半平面交模板
妹的,一直没有想清楚无解的情况到底是如何判断的。
偷来一个模板。
半平面交的结果:1.凸多边形(后面会讲解到)2.无界,因为有可能若干半平面没有形成封闭3.直线,线段,点,空(属于特殊情况吧)
算法:1:根据上图可以知道,运用给出的多边形每相邻两点形成一条直线来切割原有多边形,如果多边形上的点i在有向直线的左边或者在直线上即保存起来,否则判断此点的前一个点i-1和后一个点i+1是否在此直线的左边或线上,在的话分别用点i和点i-1构成的直线与此时正在切割的直线相交求出交点,这个交点显然也要算在切割后剩下的多边形里,同理点i和点i+1。原多边形有n条边,每条边都要进行切割,所以时间复杂度为O(n^2)。
2:第二种就是训练指南上面详细讲解的运用双端队列的半平面交算法,时间复杂度为O(nlogn)。仔细阅读代码应该能理解。
代码实现:以 poj3130 为例,裸的模板。
#include<iostream> #include<cstdio> #include<cstring> #include<cstdlib> #include<cmath> #include<algorithm> #define inf 0x7fffffff #define exp 1e-10 #define PI 3.141592654 using namespace std; const int maxn=111; struct Point { double x,y; Point (double x=0,double y=0):x(x),y(y){} }an[maxn],bn[maxn],cn[maxn]; ///an:记录最开始的多边形;bn:临时保存新切割的多边形;cn:保存新切割出的多边形 typedef Point Vector; Vector operator + (Vector A,Vector B) {return Vector(A.x+B.x , A.y+B.y); } Vector operator - (Vector A,Vector B) {return Vector(A.x-B.x , A.y-B.y); } Vector operator * (Vector A,double p) {return Vector(A.x*p , A.y*p); } Vector operator / (Vector A,double p) {return Vector(A.x/p , A.y/p); } int dcmp(double x) {if (fabs(x)<exp) return 0;return x>0 ? 1 : -1; } double cross(Vector A,Vector B) { return A.x*B.y-B.x*A.y; } double A,B,C; int n,m; void getline(Point a,Point b)///获取直线 Ax + By + C = 0 { A=b.y-a.y; B=a.x-b.x; C=b.x*a.y-a.x*b.y; } ///getline()函数得到的直线和点a和点b所连直线的交点 Point intersect(Point a,Point b) { double u=fabs(A*a.x+B*a.y+C); double v=fabs(A*b.x+B*b.y+C); Point ans; ans.x=(a.x*v+b.x*u)/(u+v); ans.y=(a.y*v+b.y*u)/(u+v); return ans; } void cut()///切割,原多边形的点为顺时针存储 { int cnt=0; for (int i=1 ;i<=m ;i++) { if (A*cn[i].x + B*cn[i].y + C>=0) bn[++cnt]=cn[i]; else { if (A*cn[i-1].x + B*cn[i-1].y + C > 0) bn[++cnt]=intersect(cn[i-1],cn[i]); if (A*cn[i+1].x + B*cn[i+1].y + C > 0) bn[++cnt]=intersect(cn[i+1],cn[i]); } } for (int i=1 ;i<=cnt ;i++) cn[i]=bn[i]; cn[0]=bn[cnt]; cn[cnt+1]=bn[1]; m=cnt;///新切割出的多边形的点数 } void solve() { for (int i=1 ;i<=n ;i++) cn[i]=an[i]; an[n+1]=an[1]; cn[n+1]=an[1]; cn[0]=an[n]; m=n; for (int i=1 ;i<=n ;i++) { getline(an[i],an[i+1]); cut(); } } int main() { while (scanf("%d",&n)!=EOF && n) { for (int i=1 ;i<=n ;i++) scanf("%lf%lf",&an[i].x,&an[i].y); reverse(an+1,an+n+1); solve(); if (m) puts("1"); else puts("0"); } return 0; }
View Code
poj3335 ,给出两种方法
1.O(n^2)
#include<iostream> #include<cstdio> #include<cstring> #include<cstdlib> #include<cmath> #include<algorithm> #define inf 0x7fffffff #define exp 1e-10 #define PI 3.141592654 using namespace std; const int maxn=111; struct Point { double x,y; Point (double x=0,double y=0):x(x),y(y){} }an[maxn],bn[maxn],cn[maxn]; typedef Point Vector ; Vector operator + (Vector A,Vector B) {return Vector(A.x+B.x , A.y+B.y); } Vector operator - (Vector A,Vector B) {return Vector(A.x-B.x , A.y-B.y); } Vector operator * (Vector A,double p) {return Vector(A.x*p , A.y*p); } Vector operator / (Vector A,double p) {return Vector(A.x/p , A.y/p); } int dcmp(double x) {if (fabs(x)<exp) return 0;return x>0 ? 1 : -1; } double cross(Vector A,Vector B) { return A.x*B.y-B.x*A.y; } int n,m; double A,B,C; void getline(Point a,Point b) { A=b.y-a.y; B=a.x-b.x; C=b.x*a.y-a.x*b.y; } Point intersect(Point a,Point b) { double u=fabs(A*a.x+B*a.y+C); double v=fabs(A*b.x+B*b.y+C); Point ans; ans.x=(a.x*v+b.x*u)/(u+v); ans.y=(a.y*v+b.y*u)/(u+v); return ans; } void cut() { int cnt=0; for (int i=1 ;i<=m ;i++) { if (A*cn[i].x+B*cn[i].y+C>=0) bn[++cnt]=cn[i]; else { if (A*cn[i-1].x+B*cn[i-1].y+C>0) bn[++cnt]=intersect(cn[i-1],cn[i]); if (A*cn[i+1].x+B*cn[i+1].y+C>0) bn[++cnt]=intersect(cn[i+1],cn[i]); } } for (int i=1 ;i<=cnt ;i++) cn[i]=bn[i]; cn[0]=bn[cnt]; cn[cnt+1]=bn[1]; m=cnt; } void solve() { for (int i=1 ;i<=n ;i++) cn[i]=an[i]; an[n+1]=an[1]; cn[n+1]=cn[1]; cn[0]=cn[n]; m=n; for (int i=1 ;i<=n ;i++) { getline(an[i],an[i+1]); cut(); } } int main() { int t; scanf("%d",&t); while (t--) { scanf("%d",&n); for (int i=1 ;i<=n ;i++) scanf("%lf%lf",&an[i].x,&an[i].y); solve(); if (!m) printf("NO\n"); else printf("YES\n"); } return 0; }
View Code
2.O(nlogn)
#include<iostream> #include<cstdio> #include<cstring> #include<cstdlib> #include<cmath> #include<algorithm> #define inf 0x7fffffff #define exp 1e-10 #define PI 3.141592654 using namespace std; const int maxn=111; struct Point { double x,y; Point (double x=0,double y=0):x(x),y(y){} }an[maxn]; typedef Point Vector; struct Line { Point p; Vector v; double ang; Line (){} Line (Point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x); } //Line (Point p,Vector v):p(p),v(v) {ang=atan2(v.y,v.x); } friend bool operator < (Line a,Line b) { return a.ang<b.ang; } }bn[maxn]; Vector operator + (Vector A,Vector B) {return Vector(A.x+B.x , A.y+B.y); } Vector operator - (Vector A,Vector B) {return Vector(A.x-B.x , A.y-B.y); } Vector operator * (Vector A,double p) {return Vector(A.x*p , A.y*p); } Vector operator / (Vector A,double p) {return Vector(A.x/p , A.y/p); } int dcmp(double x) {if (fabs(x)<exp) return 0;return x>0 ? 1 : -1; } double cross(Vector A,Vector B) { return A.x*B.y-B.x*A.y; } bool OnLeft(Line L,Point p) { return cross(L.v,p-L.p)>=0;///点P在有向直线L的左边(>=0说明在线上也算) } Point GetIntersection(Line a,Line b) { Vector u=a.p-b.p; double t=cross(b.v,u)/cross(a.v,b.v); return a.p+a.v*t; } //Point GetIntersection(Line l1, Line l2) { // Point p; // double dot1,dot2; // //dot1 = multi(l2.a, l1.b, l1.a); // dot1=cross(l1.b-l2.a , l1.a-l2.a); // //dot2 = multi(l1.b, l2.b, l1.a); // dot2=cross(l2.b-l1.b , l1.a-l1.b); // p.x = (l2.a.x * dot2 + l2.b.x * dot1) / (dot2 + dot1); // p.y = (l2.a.y * dot2 + l2.b.y * dot1) / (dot2 + dot1); // return p; //} int HalfplaneIntersection(Line *L,int n,Point *poly) { sort(L,L+n); int first,last; Point *p=new Point[n]; Line *q=new Line[n]; q[first=last=0]=L[0]; for (int i=1 ;i<n ;i++) { while (first<last && !OnLeft(L[i],p[last-1])) last--; while (first<last && !OnLeft(L[i],p[first])) first++; q[++last]=L[i]; if (fabs(cross(q[last].v , q[last-1].v))<exp) { last--; if (OnLeft(q[last] , L[i].p)) q[last]=L[i]; } if (first<last) p[last-1]=GetIntersection(q[last-1],q[last]); } while (first<last && !OnLeft(q[first],p[last-1])) last--; if (last-first<=1) return 0; p[last]=GetIntersection(q[last],q[first]); int m=0; for (int i=first ;i<=last ;i++) poly[m++]=p[i]; return m; } void calPolygon(Point *p, int n, double &area, bool &shun) { p[n] = p[0]; area = 0; double tmp; for (int i = 0; i < n; i++) area += p[i].x * p[i + 1].y - p[i].y * p[i + 1].x; area /= 2.0; if (shun = area < 0) area = -area; } bool calCore(Point *ps, int n) { Line l[maxn]; ps[n] = ps[0]; bool shun; double area; calPolygon(ps, n, area, shun); if (shun) for (int i = 0; i < n; i++) bn[i] = Line(ps[i], ps[i] - ps[i + 1]); else for (int i = 0; i < n; i++) bn[i] = Line(ps[i], ps[i + 1] - ps[i]); Point pp[maxn]; return HalfplaneIntersection(bn, n, pp); } int main() { int t; int n; scanf("%d",&t); while (t--) { scanf("%d",&n); Point cn[maxn]; for (int i=0 ;i<n ;i++) { scanf("%lf%lf",&cn[i].x,&cn[i].y); } // reverse(cn,cn+n); // for (int i=0 ;i<n ;i++) // { // bn[i].p=cn[i]; // bn[i].v=cn[(i+1)%n]-cn[i]; // bn[i].ang=atan2(bn[i].v.y , bn[i].v.x); // } // int m=HalfplaneIntersection(bn,n,an); if (!calCore(cn,n)) puts("NO"); else puts("YES"); } return 0; }
练习:
后续:欢迎提出宝贵的意见。