<模板><计算几何>半平面求交学习小记

主要是依靠这篇博文学习的 http://blog.csdn.net/accry/article/details/6070621

我的半平面交代码模板也参考自这里，个人进行了简单优化。

刚发现个很水的变换时针方法。如果你需要逆时针，而题目给的顺时针，那么读入的时候改成

for(int i=n-1; i>=0; i--)scanf("%lf%lf",&p[i].x,&p[i].y);

二分double型改变大小的时候low=mid比low=mid+eps要好一些(个人感觉)

半平面求交 及 多边形求核 模板

#include <cstdio>#include <iostream>#include <cmath>#include <algorithm>using namespace std;const double EPS=1e-11;const int NUM=1510;int DB (double x){    if (x>EPS)        return 1;    if (x<-EPS)        return -1;    return 0;}struct Point{    double x,y;     Point(){}    Point(double _x,double _y)    {        x=_x;        y=_y;    }     void get ()    {scanf("%lf%lf",&x,&y);}    void Out ()    {        printf("(%.5lf %.5lf)",x,y);    }     Point operator+(Point a)    {        return Point(x+a.x,y+a.y);    }    Point operator-(Point a)    {        return Point(x-a.x,y-a.y);    }     double operator*(Point a)     //叉积    {        return x*a.y-y*a.x;    }double operator^(Point a)     //点积    {        return x*a.y+y*a.x;    }    Point operator*(double t)    {        return Point(x*t,y*t);    }     Point operator/(double t)    {        return Point(x/t,y/t);    }     bool operator==(Point a)    {        return DB(x-a.x)==0&&DB(y-a.y)==0;    }     bool operator!=(Point a)    {        return DB(x-a.x)||DB(y-a.y);    }}points[NUM],p[NUM],q[NUM];int n;double r;int cCnt,curCnt;void getline (Point x,Point y,double &a,double &b,double &c){a = y.y - x.y;b = x.x - y.x;c = y.x * x.y - x.x * y.y;}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) );}/*半平面相交（直线切割多边形）（点标号从1开始）*/void cut (double a,double b ,double c){int i;curCnt = 0;for (i=1;i<=cCnt;i++){        if (DB(a*p[i].x + b*p[i].y + c) >= 0)q[++curCnt] = p[i];        else{            if (DB(a*p[i-1].x + b*p[i-1].y + c) > 0)                q[++curCnt] = intersect(p[i],p[i-1],a,b,c);            if (DB(a*p[i+1].x + b*p[i+1].y + c) > 0)                q[++curCnt] = intersect(p[i],p[i+1],a,b,c);        }    }    for (i=1;i<=curCnt;i++)p[i] = q[i];    p[curCnt+1] = q[1];p[0] = p[curCnt];    cCnt = curCnt;}void GuiZhengHua (){//规整化方向,逆时针变顺时针,顺时针变逆时针for (int i=1;i<(n+1)/2;i++)swap(points[i], points[n-i]);}double getArea (Point p[],int n){double ans=0;for (int i=0;i<n;i++)ans+=p[i]*p[i+1];return ans/2;}void initial (){    for (int i=1;i<=n;i++)p[i] = points[i];  p[n+1] = p[1];p[0] = p[n];cCnt = n;double s=getArea(p,n);     //自动规整化为顺时针,如果方向确定,可不执行if (DB(s)>=0)      //面积为正,说明为逆时针GuiZhengHua ();}bool Deal (){int i;//注意:默认点是顺时针,如果题目不是顺时针,规整化方向  initial();for (i=1;i<=n;i++){double a,b,c;getline (points[i],points[i+1],a,b,c);cut(a,b,c);}       //p[0]~p[cCnt-1]及p[1]~p[cCnt]均存放多边形的核的cCnt-1个顶点    /*     如果要向内推进r，用该部分代替上个函数     for (i=1;i<=n;i++){         Point ta, tb, tt;         tt.x = points[i+1].y - points[i].y;         tt.y = points[i].x - points[i+1].x;         double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);         tt.x = tt.x * k;         tt.y = tt.y * k;         ta.x = points[i].x + tt.x;         ta.y = points[i].y + tt.y;         tb.x = points[i+1].x + tt.x;         tb.y = points[i+1].y + tt.y;         double a,b,c;         getline(ta,tb,a,b,c);         cut(a,b,c);     }     */// 计算面积printf("%.2lf\n",fabs(getArea (p,cCnt)));    //顺时针处理后面积为负,但写成-1.0*会WA……if (cCnt == 0)return false;return true;} int main (){int T;scanf("%d",&T);while (T--){scanf("%d",&n);for (int i=1;i<=n;i++)points[i].get();points[0]=points[n];points[n+1] = points[1];Deal ();}return 0;}

nlogn算法,Poj2451

个人感觉比较难理解，据说该算法特殊情况有bug……

代码参考自：//http://hi.baidu.com/godwitness/item/8107cb1f6a21d7f987ad4e84

#include <cstdio>#include <iostream>#include <cmath>#include <algorithm>using namespace std;const double STD=1e-10,big=100000.0;const int NUM = 20010;struct Point{double x,y;void get (){scanf("%lf%lf",&x,&y);}}node[NUM];struct polygon   //存放最后半平面交中相邻边的交点,就是一个多边形的所有点{int n;           //多边形顶点数Point p[NUM];    //从p[0]开始保存}pg;struct Line      //半平面,这里是线段{Point a,b;}data[NUM];        //半平面集合double at2[NUM];   //每条线段对应坐标系的角度int ord[NUM];int dq[NUM+1];     //双端栈int n,lnum;int DB (double x)  {      if (x>STD)          return 1;      if (x<-STD)          return -1;      return 0;  } //叉积>0代表在左边,<0代表在右边,=0代表共线//e是否在o->s的左边onleft(DB(cross))>=0double cross (Point o, Point s, Point e) {return (s.x-o.x)*(e.y-o.y)-(e.x-o.x)*(s.y-o.y);}//直线求交点Point isIntersected (Point s1, Point e1, Point s2, Point e2){double dot1,dot2;Point pp;dot1 = cross(s2,e1,s1); dot2 = cross(e1,e2,s1); pp.x = (s2.x * dot2 + e2.x * dot1) / (dot2 + dot1); pp.y = (s2.y * dot2 + e2.y * dot1) / (dot2 + dot1); return pp;} //象限排序bool cmp (int u,int v){if (DB(at2[u]-at2[v]) == 0)return DB(cross(data[v].a,data[v].b,data[u].b))>=0;return at2[u]<at2[v];} //判断半平面的交点在当前半平面外bool Judgein (int x,int y,int z){Point pnt = isIntersected(data[x].a, data[x].b, data[y].a, data[y].b); //求交点return DB(cross(data[z].a,data[z].b,pnt)) < 0; //判断交点位置,如果在右面,返回true,如果要排除三点共线,改成<=}//半平面交void HalfPlaneIntersection (polygon &pg){ //预处理int n=lnum,tmpn,i; /* 对于atan2(y,x)结果为正表示从 X 轴逆时针旋转的角度，结果为负表示从 X 轴顺时针旋转的角度.atan2(a, b) 与 atan(a/b)稍有不同,atan2(a,b)的取值范围介于 -pi 到 pi 之间（不包括 -pi）,而atan(a/b)的取值范围介于-pi/2到pi/2之间*/for (i=0;i<n;i++) { //atan2(y,x)求出每条线段对应坐标系的角度at2[i] = atan2( data[i].b.y - data[i].a.y, data[i].b.x - data[i].a.x);ord[i] = i;}sort (ord,ord+n,cmp);for (i=1,tmpn=1;i<n;i++) //处理重线的情况if( DB(at2[ord[i-1]] - at2[ord[i]]) != 0 )ord[tmpn++] = ord[i];n = tmpn; //圈地int bot = 1,top = bot + 1; //双端栈,bot为栈底,top为栈顶dq[bot] = ord[0]; dq[top] = ord[1]; //先压两根线进栈for (i=2;i<n;i++)    { //bot < top 表示要保证栈里至少有2条线段,如果剩下1条,就不继续退栈//Judgein,判断如果栈中两条线的交点如果在当前插入先的右边,就退栈while ( bot < top && Judgein(dq[top-1] , dq[top] , ord[i]) )top--;//对栈顶要同样的操作while ( bot < top && Judgein(dq[bot+1] , dq[bot] , ord[i]) )bot++; dq[++top] = ord[i]; } //最后还要处理一下栈里面存在的栈顶的线在栈底交点末尾位置,或者栈顶在栈尾两条线的右边while ( bot < top && Judgein(dq[top-1] , dq[top] , dq[bot]) )top--; while ( bot < top && Judgein(dq[bot+1] , dq[bot] , dq[top]) )bot++; //最后一条线是重合的dq[--bot] = dq[top]; //求多边形pg.n = 0; for (i=bot+1;i<=top;i++) //求相邻两条线的交点pg.p[pg.n++] = isIntersected(data[dq[i-1]].a, data[dq[i-1]].b, data[dq[i]].a,data[dq[i]].b);}void add (double a,double b,double c,double d)  //添加线段{data[lnum].a.x = a; data[lnum].a.y = b; data[lnum].b.x = c; data[lnum].b.y = d;lnum++;}void init ()           //计算时要求多边形的核一定位于半平面左侧{ int i;scanf("%d",&n);/*double a,b,c,d;for (i=0;i<n;i++){//输入代表一条向量(x = (c - a),y = (d - b)),核在其左侧;scanf ("%lf%lf%lf%lf",&a,&b,&c,&d);add (a,b,c,d);}*/for (i=0;i<n;i++)       //以多边形顶点输入,要求逆时针node[i].get();data[n]=data[i];for (i=0;i<n;i++)add(node[i].x,node[i].y,node[i+1].x,node[i+1].y);//下面是构造一个大矩形边界add(0,0,big,0);       //downadd(big,0,big,big);   //rightadd(big,big,0,big);   //upadd(0,big,0,0);       //left}double Deal (){init ();HalfPlaneIntersection(pg); //求半平面交//最后多边形的各个点保存在pg里面//计算半平面交的面积double area = 0;n = pg.n;for (int i=0;i<n;i++) area += pg.p[i].x * pg.p[(i+1)%n].y - pg.p[(i+1)%n].x * pg.p[i].y;    //x1 * y2 - x2 * y1用叉积求多边形面积 area=fabs(area)/2.0;           //所有面积应该是三角形面积之和，而叉积求出来的是四边形的面积和,所以要除2return area;}