二维计算几何
声明:
由于本人较弱,并不能保证以下内容的100%正确
欢迎大佬来挑错
基本定义
1 struct Point{ 2 double x,y; 3 Point(){ 4 x=y=0.0; 5 } 6 Point(double xx,double yy){ 7 x=xx;y=yy; 8 } 9 }; 10 typedef Point Vec; 11 12 Vec operator + (Vec v1,Vec v2){ 13 return Vec(v1.x+v2.x,v1.y+v2.y); 14 } 15 Vec operator - (Point p1,Point p2){ 16 return Vec(p1.x-p2.x,p1.y-p2.y); 17 } 18 Vec operator * (Vec v,double k){ 19 return Vec(v.x*k,v.y*k); 20 } 21 Vec operator / (Vec v,double k){ 22 return Vec(v.x/k,v.y/k); 23 } 24 25 bool operator < (const Point &a,const Point &b){ 26 if(a.x==b.x)return a.y<b.y; 27 else return a.x<b.x; 28 } 29 30 double Myabs(double x){ 31 if(x<0)return -x; 32 else return x; 33 } 34 35 const double eps=1e-10; 36 int dcmp(double x){ 37 if(Myabs(x)<eps)return 0; 38 if(x>0)return 1; 39 else return -1; 40 } 41 42 bool operator == (const Point &a,const Point &b){ 43 return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0; 44 }
点积
1 double Dt(Vec v1,Vec v2){ 2 return v1.x*v2.x+v1.y*v2.y; 3 } 4 double Getlen(Vec v){ 5 return sqrt(Dt(v,v)); 6 } 7 double Getang(Vec v1,Vec v2){ 8 return acos(Dt(v1,v2)/Getlen(v1)/Getlen(v2)); 9 }
叉积
1 double Crs(Vec v1,Vec v2){ 2 return v1.x*v2.y-v2.x*v1.y; 3 } 4 double GetS2(Point A,Point B,Point C){ 5 return Crs(B-A,C-A); 6 }
基本运算
1 Vec Rotate(Vec v,double rad){ 2 return Vec(v.x*cos(rad)-v.y*sin(rad),v.x*sin(rad)+v.y*cos(rad)); 3 } 4 Vec Getn(Vec v){ 5 double L=Getlen(v); 6 return Vec(-v.y/L,v.x/L); 7 }
1 Point GetLineIntersection(Point P,Vec v,Point Q,Vec w){ 2 Vec u=P-Q; 3 double t=Crs(w,u)/Crs(v,w); 4 return P+v*t; 5 } 6 double DistanceToLine(Point P,Point A,Point B){ 7 Vec v1=B-A,v2=P-A; 8 return Myabs(Crs(v1,v2))/Getlen(v1); 9 } 10 Point GetLineProjection(Point P,Point A,Point B){ 11 Vec v=B-A; 12 double t=Dt(v,P-A)/Dt(v,v); 13 return A+v*t; 14 } 15 16 17 bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){ 18 double c1=Crs(a2-a1,b1-a1),c2=Crs(a2-a1,b2-a1); 19 double c3=Crs(b2-b1,a1-b1),c4=Crs(b2-b1,a2-b1); 20 return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0; 21 } 22 bool OnSegment(Point p,Point a1,Point a2){ 23 return dcmp(Crs(a1-p,a2-p))==0&&dcmp(Dt(a1-p,a2-p))<0; 24 } 25 26 double PolygonArea(Point* p,int n){ 27 double area=0; 28 for(int i=1;i<n-1;++i){ 29 area+=Crs(p[i]-p[0],p[i+1]-p[0]); 30 } 31 return Myabs(area)/2; 32 }
二维计算几何常用算法
点在多边形内的判定
1 int isPointInPolygon(Point p,vector<Point>poly){ 2 int wn=0; 3 int n=poly.size(); 4 for(int i=0;i<n;++i){ 5 if(OnSegment(p,poly[i],poly[(i+1)%n]))return -1; 6 int k=dcmp(Crs(poly[(i+1)%n]-poly[i],p-poly[i])); 7 int d1=dcmp(poly[i].y-p.y); 8 int d2=dcmp(poly[(i+1)%n].y-p.y); 9 if(k>0 && d1<=0 &&d2>0)wn++; 10 if(k<0 && d2<=0 &&d1>0)wn--; 11 } 12 if(wn!=0)return 1; 13 else return 0; 14 }
二维凸包
注意
输入不能有重复点
精度高时使用dcmp比较
1 int ConvexHull(Point *p,int n,Point *ch){ 2 sort(p,p+n); 3 int m=0; 4 for(int i=0;i<n;++i){ 5 while(m>1&&Crs(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)--m; 6 ch[m++]=p[i]; 7 } 8 int k=m; 9 for(int i=n-2;i>=0;--i){ 10 while(m>k&&Crs(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)--m; 11 ch[m++]=p[i]; 12 } 13 if(n>1)m--; 14 return m; 15 }
旋转卡壳求直径
1 double RotatingCaliper(Point* ch,int m){ 2 int q=1; 3 ch[m]=ch[0]; 4 double ans=0; 5 for(int p=0;p<m;++p){ 6 while(Myabs(Crs(ch[q+1]-ch[p+1],ch[p]-ch[p+1]))>Myabs(Crs(ch[q]-ch[p+1],ch[p]-ch[p+1])))q=(q+1)%m; 7 ans=max(ans,max(Getlen(ch[p]-ch[q]),Getlen(ch[p+1]-ch[q+1]))); 8 ans=max(ans,max(Getlen(ch[p]-ch[q+1]),Getlen(ch[p+1]-ch[q]))); 9 } 10 return ans; 11 }
半平面交
O(n2)
半平面交注意判断无解和无界
1 typedef vector<Point> Polygon; 2 Polygon CutPolygon(Polygon poly,Point A,Point B){ 3 Polygon newpoly; 4 int n=poly.size(); 5 for(int i=0;i<n;++i){ 6 Point C=poly[i]; 7 Point D=poly[(i+1)%n]; 8 if(dcmp(Crs(B-A,C-A))>=0)newpoly.push_back(C); 9 if(dcmp(Crs(B-A,C-D))!=0){ 10 Point ip=GetLineIntersection(A,B-A,C,D-C); 11 if(OnSegment(ip,C,D))newpoly.push_back(ip); 12 } 13 } 14 return newpoly; 15 }
O(nlogn)
1 bool OnLeft(Line L,Point p){
2 return Crs(L.v,p-L.P)>0;
3 }
4 Point GetLineIntersection(Line a,Line b){
5 Vec u=a.P-b.P;
6 double t=Crs(b.v,u)/Crs(a.v,b.v);
7 return a.P+a.v*t;
8 }
9
10 Point p[maxn];
11 Line q[maxn];
12 int HalfPlaneIntersection(Line *L,int n,Point *poly){
13 sort(L,L+n);
14
15 int h,t;
16 q[h=t=1]=L[0];
17 for(int i=1;i<n;++i){
18 while((h<t)&&(!OnLeft(L[i],p[t-1])))--t;
19 while((h<t)&&(!OnLeft(L[i],p[h])))++h;
20 q[++t]=L[i];
21 if(dcmp(Crs(q[t].v,q[t-1].v))==0){
22 --t;
23 if(OnLeft(q[t],L[i].P))q[t]=L[i];
24 }
25 if(h<t)p[t-1]=GetLineIntersection(q[t-1],q[t]);
26 }
27 while((h<t)&&(!OnLeft(q[h],p[t-1])))--t;
28 if(t-h<=1)return 0;
29 p[t]=GetLineIntersection(q[t],q[h]);
30
31 int m=0;
32 for(int i=h;i<=t;++i)poly[m++]=p[i];
33 return m;
34 }
bool OnLeft(Line L,Point p){ return Crs(L.v,p-L.P)>0; } Point GetLineIntersection(Line a,Line b){ Vec u=a.P-b.P; double t=Crs(b.v,u)/Crs(a.v,b.v); return a.P+a.v*t; } Point p[maxn]; Line q[maxn]; int HalfPlaneIntersection(Line *L,int n,Point *poly){ sort(L,L+n); int h,t; q[h=t=1]=L[0]; for(int i=1;i<n;++i){ while((h<t)&&(!OnLeft(L[i],p[t-1])))--t; while((h<t)&&(!OnLeft(L[i],p[h])))++h; q[++t]=L[i]; if(dcmp(Crs(q[t].v,q[t-1].v))==0){ --t; if(OnLeft(q[t],L[i].P))q[t]=L[i]; } if(h<t)p[t-1]=GetLineIntersection(q[t-1],q[t]); } while((h<t)&&(!OnLeft(q[h],p[t-1])))--t; if(t-h<=1)return 0; p[t]=GetLineIntersection(q[t],q[h]); int m=0; for(int i=h;i<=t;++i)poly[m++]=p[i]; return m; }