多边形的核（半平面交算法）

//看代码前最好看下zzy的论文！！！！王飞的资料里有 #include<stdio.h> #include<math.h> #include<algorithm> #include<string.h> using namespace std; const double eps = 1e-8; const int MaxN = 50005; typedef struct{double x, y;}Point; typedef struct{double k, b;}Line; typedef struct{Point p[MaxN]; int n;}Polygon; Line l[MaxN]; double ang[MaxN]; int n, ord[MaxN], dq[MaxN]; int dblcmp(double t) { if(fabs(t) < eps) return 0; return t < 0? -1 : 1; } double cross(Point p0, Point p1, Point p2) //叉积 { return (p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x); } bool cmp(int u, int v) //极角排序 { if(!dblcmp(ang[u]-ang[v]))//当极角一样是，靠近目标的放前面 return l[u].b > l[v].b; return ang[u] < ang[v]; } void Intersection(int x, int y, Point &pp) //y=kx+b,判断两直线的交点pp { pp.x = (l[x].b-l[y].b)/(l[y].k-l[x].k); pp.y = l[x].k*pp.x + l[x].b; } /*void Intersection(Point s1, Point e1, Point s2, Point e2, Point &pp) { double dot1,dot2; 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); }两个点求两直线的交点pp*/ int Judgein(int x, int y, int z) //用叉积判断点是靠近目标的还是远离目标的 { Point tmp, t1, t2; Intersection(y, z, tmp); t1.x = 1; t1.y = l[x].k + l[x].b; t2.x = 0; t2.y = l[x].b; return cross(t2, t1, tmp) <= 0; } int HalfPlaneIntersection() { int cnt, i; for(i = 0; i < n; i++) //求出与X轴夹角，用于极角排序 { ang[i] = atan2(l[i].k, 1); ord[i] = i; } sort(ord, ord+n, cmp); for(i = cnt = 1; i < n; i++) if(dblcmp(ang[ord[i]]-ang[ord[i-1]])) ord[cnt++] = ord[i]; int bot = 1, top = 2; dq[1] = ord[0]; // 双端队列 dq[2] = ord[1]; for(i = 2; i < n; i++) { while(bot < top && Judgein(ord[i], dq[top], dq[top-1])) top--; while(bot < top && Judgein(ord[i], dq[bot], dq[bot+1])) bot++; dq[++top] = ord[i]; } while(bot < top && Judgein(dq[bot], dq[top], dq[top-1])) top--; while(bot < top && Judgein(dq[top], dq[bot], dq[bot+1])) bot++; //得到平面的交点 }