二维计算几何复习

二维计算几何

 

声明：

由于本人较弱，并不能保证以下内容的100%正确

欢迎大佬来挑错

 

 

基本定义

struct Point{
    double x,y;
    Point(){
        x=y=0.0;
    }
    Point(double xx,double yy){
        x=xx;y=yy;
    }
};
typedef Point Vec;

Vec operator + (Vec v1,Vec v2){
    return Vec(v1.x+v2.x,v1.y+v2.y);
}
Vec operator - (Point p1,Point p2){
    return Vec(p1.x-p2.x,p1.y-p2.y);
}
Vec operator * (Vec v,double k){
    return Vec(v.x*k,v.y*k);
}
Vec operator / (Vec v,double k){
    return Vec(v.x/k,v.y/k);
}

bool operator < (const Point &a,const Point &b){
    if(a.x==b.x)return a.y<b.y;
    else return a.x<b.x;
}

double Myabs(double x){
    if(x<0)return -x;
    else return x;
}

const double eps=1e-10;
int dcmp(double x){
    if(Myabs(x)<eps)return 0;
    if(x>0)return 1;
    else return -1;
}

bool operator == (const Point &a,const Point &b){
    return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}

 

点积

double Dt(Vec v1,Vec v2){
    return v1.x*v2.x+v1.y*v2.y;
}
double Getlen(Vec v){
    return sqrt(Dt(v,v));
}
double Getang(Vec v1,Vec v2){
    return acos(Dt(v1,v2)/Getlen(v1)/Getlen(v2));
}

 

叉积

double Crs(Vec v1,Vec v2){
    return v1.x*v2.y-v2.x*v1.y;
}
double GetS2(Point A,Point B,Point C){
    return Crs(B-A,C-A);
}

 

基本运算

 

Vec Rotate(Vec v,double rad){
    return Vec(v.x*cos(rad)-v.y*sin(rad),v.x*sin(rad)+v.y*cos(rad));
}
Vec Getn(Vec v){
    double L=Getlen(v);
    return Vec(-v.y/L,v.x/L);
}

 

Point GetLineIntersection(Point P,Vec v,Point Q,Vec w){
    Vec u=P-Q;
    double t=Crs(w,u)/Crs(v,w);
    return P+v*t;
}
double DistanceToLine(Point P,Point A,Point B){
    Vec v1=B-A,v2=P-A;
    return Myabs(Crs(v1,v2))/Getlen(v1);
}
Point GetLineProjection(Point P,Point A,Point B){
    Vec v=B-A;
    double t=Dt(v,P-A)/Dt(v,v);
    return A+v*t;
}


bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
    double c1=Crs(a2-a1,b1-a1),c2=Crs(a2-a1,b2-a1);
    double c3=Crs(b2-b1,a1-b1),c4=Crs(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
bool OnSegment(Point p,Point a1,Point a2){
    return dcmp(Crs(a1-p,a2-p))==0&&dcmp(Dt(a1-p,a2-p))<0;
}

double PolygonArea(Point* p,int n){
    double area=0;
    for(int i=1;i<n-1;++i){
        area+=Crs(p[i]-p[0],p[i+1]-p[0]);
    }
    return Myabs(area)/2;
}

 

二维计算几何常用算法

 

点在多边形内的判定

int isPointInPolygon(Point p,vector<Point>poly){
    int wn=0;
    int n=poly.size();
    for(int i=0;i<n;++i){
        if(OnSegment(p,poly[i],poly[(i+1)%n]))return -1;
        int k=dcmp(Crs(poly[(i+1)%n]-poly[i],p-poly[i]));
        int d1=dcmp(poly[i].y-p.y);
        int d2=dcmp(poly[(i+1)%n].y-p.y);
        if(k>0 && d1<=0 &&d2>0)wn++;
        if(k<0 && d2<=0 &&d1>0)wn--;
    }
    if(wn!=0)return 1;
    else return 0;
}

 

二维凸包

注意

输入不能有重复点

精度高时使用dcmp比较

int ConvexHull(Point *p,int n,Point *ch){
    sort(p,p+n);
    int m=0;
    for(int i=0;i<n;++i){
        while(m>1&&Crs(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)--m;
        ch[m++]=p[i];
    }
    int k=m;
    for(int i=n-2;i>=0;--i){
        while(m>k&&Crs(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)--m;
        ch[m++]=p[i];
    }
    if(n>1)m--;
    return m;
}

 

旋转卡壳求直径

double RotatingCaliper(Point* ch,int m){
    int q=1;
    ch[m]=ch[0];
    double ans=0;
    for(int p=0;p<m;++p){
        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;
        ans=max(ans,max(Getlen(ch[p]-ch[q]),Getlen(ch[p+1]-ch[q+1])));
        ans=max(ans,max(Getlen(ch[p]-ch[q+1]),Getlen(ch[p+1]-ch[q])));
    }
    return ans;
}

 

半平面交

 O(n²)

半平面交注意判断无解和无界

typedef vector<Point> Polygon;
Polygon CutPolygon(Polygon poly,Point A,Point B){
    Polygon newpoly;
    int n=poly.size();
    for(int i=0;i<n;++i){
        Point C=poly[i];
        Point D=poly[(i+1)%n];
        if(dcmp(Crs(B-A,C-A))>=0)newpoly.push_back(C);
        if(dcmp(Crs(B-A,C-D))!=0){
            Point ip=GetLineIntersection(A,B-A,C,D-C);
            if(OnSegment(ip,C,D))newpoly.push_back(ip);
        }
    }
    return newpoly;
}

O(nlogn)

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;
}

 

 

 

 

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;
}