二维计算几何复习

二维计算几何

 

声明:

由于本人较弱,并不能保证以下内容的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;
}

  

 

 

 

 

 

 

 

 

 

 

 

 
posted @ 2018-04-26 17:27  ws_zzy  阅读(189)  评论(0编辑  收藏  举报