hdu 4327 简单立体几何 + 半平面交

http://acm.hdu.edu.cn/showproblem.php?pid=4327

题意：给你一些点( x, y ), 0 < x < 1 , 0 < y < 1.再给你一个概率分布函数,表示每个坐标点的概率分布，显然整个矩形的分布概率之和为1

另外有一个规则：每个点都只管辖离自己最近的区域，即和其他所有点的中垂线所切割成的区域。

求每块区域分布的概率是多少

解法：显然，每个点所管辖的区域可以用这个点与其他点得中垂线去切割平面，这个点所在的一侧为有效半平面，所以要先判断点与直线的方位再去切割平面，这里我写了两个切割的函数，

容易看出，题目所给的概率分布函数是空间中的一个平面，所以可以想象以半平面切割后的区域为底，最顶上是一个斜面的立体（这个空间立体要仔细想想）

那么这个立体的体积就是我们要求的答案。

具体在求的时候可以先求出底下棱柱的体积，最上面的那个角单独计算。

咋一看，最上面的部分什么形状都没有，没办法算啊，仔细一看就知道了，可以将它切割成若干个四棱锥，锥顶始终是最矮的那个点，这个需要大家发挥空间想象能力。

这个题搞了我一个下午+晚上的时间啊，最后发现错在半平面交 模板，我的模板是n^2的，以前写的时候都没有注意一个很小的细节，在这道题上暴露出来了，发现问题，及时解决总是好的，go on AC！

ps：半平面交学习的话可以去这里，写的很好

#include <cmath>
#include <cstdio>
#include<algorithm>
using namespace std;
const  double INF = 1000000000.0;
const  int maxn = 1010;
const double eps = 1e-12;
inline double sgn(double x) {return fabs(x)<eps?0:(x>0?1:-1);}
struct Point{
    double x,y;
    Point(double tx=0,double ty=0){x=tx;y=ty;}
    bool operator == (const Point& t) const {
        return sgn(x-t.x)==0 && sgn(y-t.y)==0;
    }
     Point operator - (const Point& t) const {
        Point tmp;
        tmp.x = x - t.x;
        tmp.y = y - t.y;
        return tmp;
    }
}p[maxn],tmp[maxn],pp[maxn],GP;
struct Seg{Point s,e;};
struct Line {
    double a, b, c;
};//中垂线
double cross(Point a,Point b,Point c){return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);}
inline Point intersect(Point x,Point y,double a,double b,double c){  
    double u = fabs(a * x.x + b * x.y + c);  
    double v = fabs(a * y.x + b * y.y + c);  
    return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );  
}  
int n,ban_tot;
void CUT1(double a,double b,double c){//points must be clockwise    a*x+b*y+c>=0表示点在半平面内   
    int i,tot=0;
    for(int i = 1; i <= ban_tot; ++i){  
        if(a*p[i].x + b*p[i].y + c >= eps) pp[++tot] = p[i];  
        else {  
            if(a*p[i-1].x + b*p[i-1].y + c > eps){  
                pp[++tot] = intersect(p[i],p[i-1],a,b,c);  
            }  
            if(a*p[i+1].x + b*p[i+1].y + c > eps){  
                pp[++tot] = intersect(p[i],p[i+1],a,b,c);  
            }  
        }  
    }  ban_tot=tot;
    pp[tot+1]=pp[1];pp[0]=pp[tot];
    memcpy(p,pp,sizeof(pp));
}
void CUT2(double a,double b,double c){//a*x+b*y+c<=0表示点在半平面内        
    int i,tot=0;
    for(int i = 1; i <= ban_tot; ++i){  
        if(!(a*p[i].x + b*p[i].y + c > eps) ) pp[++tot] = p[i]; //点在半平面内或者边界上 
        else {  
            if(a*p[i-1].x + b*p[i-1].y + c < -eps){  //上个点完全在半平面内
                pp[++tot] = intersect(p[i],p[i-1],a,b,c);  
            }  
            if(a*p[i+1].x + b*p[i+1].y + c < -eps){  //下个点完全在半平面内
                pp[++tot] = intersect(p[i],p[i+1],a,b,c);  
            }  
        }  
    }  ban_tot=tot;
    pp[tot+1]=pp[1];pp[0]=pp[tot];//两端都扩展一个点
    memcpy(p,pp,sizeof(pp));
}
Line Turn(Point s, Point e) {                                        // 线段转直线表达式
    Line ln;
    ln.a = s.y - e.y;
    ln.b = e.x - s.x;
    ln.c = s.x*e.y - e.x*s.y;
    return ln;
}
Line make(Point a,Point b)//求a,b中垂线的直线方程
{
    double x0=(a.x+b.x)/2;
    double y0=(a.y+b.y)/2;
    Line tmp=Turn(a,b);
    Line ans;
    ans.a=tmp.b;
    ans.b=-tmp.a;
    ans.c=tmp.a*y0-tmp.b*x0;
    return ans;
}
Line ln[maxn];
inline double PPdis(Point a, Point b) {                                // 点点距离
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline double PLdis(Point p,Point l1,Point l2){                        // 点线距离
    return fabs(cross(p,l1,l2))/PPdis(l1,l2);
}
double calc(Point *p,int n)
{
      if(n<3) return 0;
      double area=0,V=0;
      for(int i=0;i<n;i++) area+=p[(i+1)%n].y*(p[i].x-p[(i+2)%n].x);
      area/=2;
      area=fabs(area);
      double H=10;
      int pos=0;
      for(int i=0;i<n;i++)
      {
          if(2-p[i].x-p[i].y<H)
          {
                  H=2-p[i].x-p[i].y;
                  pos=i;
          }
      }
      V+=area*H;
      for(int i=pos+1;;i++)
      {  
          int id1=i%n;
          int id2=(i+1)%n;
          if(id2==pos) break;
          double h=PLdis(p[pos],p[id1],p[id2]);
          double s=((2-p[id1].x-p[id1].y-H) + (2-p[id2].x-p[id2].y-H)) * PPdis(p[id1],p[id2])/2;
          V+=s*h/3;
      }
      return V;
}
double ans[110];
int main(){
    int t,ca=1,n;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d",&n);
        for(int i=1;i<=n;i++) scanf("%lf%lf",&tmp[i].x,&tmp[i].y);
        for(int i=1;i<=n;i++)
        {
            p[1]=Point(0,0);
            p[2]=Point(0,1);
            p[3]=Point(1,1);
            p[4]=Point(1,0);
            p[0]=p[4];
            p[5]=p[1];//首尾都要加上一个点，保证边界的合法点不会漏掉
			ban_tot=4;
            double x0=tmp[i].x,y0=tmp[i].y;
            for(int j=1;j<=n;j++)
            {
                if(i==j) continue;
                Line mid_line = make(tmp[i],tmp[j]);
                double a=mid_line.a,b=mid_line.b,c=mid_line.c;
                if(a*x0+b*y0+c >= eps)
                {
                    CUT1(a,b,c);
                }
                else 
                {
                    CUT2(a,b,c);
                }
            }
            double tmpv=calc(p,ban_tot);
            ans[i]=tmpv;
        }
        printf("Case #%d:\n",ca++);
        for(int i=1;i<=n;i++)
        {
            printf("%.6lf\n",ans[i]);
        }
    }
    return 0;
}