计算几何学习(待完善)

本文最后更新于：2024年9月13日 早上

今日补完了昆明站赛站，自告奋勇接下了两题结果卡在了计算几何上，以为自己有着模板就直接上了，结果一直卡在段错误，因此决定这周进行计算几何专题的专门练习。
其实这些内容在寒假也学了一些，当时主要是入了计算机图形学的坑，玩OpenGL的时候顺便补了点计几。当时也找到了很多的blog进行学习，并且学习了大佬的代码风格(因为感觉这样写好棒).这里的内容主要还是把平常的几何知识封装以下，所以除了一些思维上的新东西应该都只有板子。

准备工作

建点

1
2
3
4
5
6
7
8
9
10
#define LD double
#define Vector Point
const LD eps=1e-8;//精度差
inline int dcmp(LD a){return a<-eps?-1:(a>eps?1:0);}//处理精度
inline LD ABS(LD a){return a*dcmp(a);}//取绝对值
struct Point{//名字是点，但实际上向量也是一样的用法
    LD x,y;Point(LD X=0,LD Y=0){x=X,y=Y;}//为点赋初值
    inline void in(){scanf("%lf%lf",&x,&y);}
    inline void out(){printf("%.2lf %.2lf\n",x,y);}
};

向量

模长

1
inline LD Len(Vector a){return sqrt(Dot(a,a));}

点积&&叉积

1
2
inline LD Dot(Vector a,Vector b){return a.x*b.x+a.y*b.y;}//点积
inline LD Cro(Vector a,Vector b){return a.x*b.y-a.y*b.x;}//叉积

求法线

1
inline Vector Normal_line(Vector a) {double len=Len(a);return Vector(-a.y/len,a.x/len);}//求向量A的单位法线

求夹角

1
inline LD Angle(Vector a,Vector b) {return acos(Dot(a,b)/len(a)/len(b));}

位置变换

向量位移重载运算符直接加减&&数乘

1
2
3
inline Vector operator+(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
inline Vector operator-(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
inline Vector operator*(Vector a,LD b){return Vector(a.x*b,a.y*b);}

点&&向量旋转

将向量顺时针旋转$\theta$角度(点的话是绕原点旋转)，相当于乘上了个旋转矩阵。

1
2
3
4
5
inline Point turn_P(Point a,LD theta){//点A\向量A顺时针旋转theta(弧度)
    LD x=a.x*cos(theta)+a.y*sin(theta);
    LD y=-a.x*sin(theta)+a.y*cos(theta);
    return Point(x,y);
}

点绕点旋转

点绕点旋转，实际上可以看作是先平移坐标系把中心点平移到原点的位置，然后点绕原点旋转后再平移回去

1
2
3
4
5
inline Point turn_PP(Point a,Point b,LD theta){//将点A绕点B顺时针旋转theta(弧度)
    LD x=(a.x-b.x)*cos(theta)+(a.y-b.y)*sin(theta)+b.x;
    LD y=-(a.x-b.x)*sin(theta)+(a.y-b.y)*cos(theta)+b.y;
    return Point(x,y);
}

点与线

判断点是否在线段上

点若在线段上，呢么其与线段上两端点的夹角为180度，其叉积为0，点积不大于0

1
2
3
4
5
inline bool pan_PL(Point p,Point a,Point b){//判断点P是否在线段AB上
    if(a==b)return false;//AB重合，线段不存在
    return !dcmp(Cro(p-a,b-a))&&dcmp(Dot(p-a,p-b))<=0;//注意，根据icpc2021昆明i题发现，由于求交点时本身就会产生精度误差，使用叉积来判断向量平行不再放宽精度则会放大此误差，因此假如已经确定点在线段所在的直线上后可以根据需要去除第一个条件 
    //return dcmp(Dot(p-a,p-b))<=0; 
}

判断点是否在直线上

点若在直线上，呢么就只剩下与直线上任意两其他点的叉积为0，无点积约束

1
2
3
inline bool pan_PL_(Point p,Point a,Point b){//判断点P是否在线段AB上
    return !dcmp(Cro(p-a,b-a));
}

点到线段AB距离

假如点与线段的垂线不在线段上，返回与最近的端点之间的距离。若是在线段上，呢么用与两端点围成三角形的面积除以底边长来算，运用了两向量叉积的数值等于两向量组成的平行四边形的面积的性质。

1
2
3
4
5
6
7
8
inline bool operator==(Point a,Point b){return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);}//判断点的相等
inline LD dis_PL(Point p,Point a,Point b){//点P到线段AB距离
    if(a==b)return Len(p-a);//AB重合
    Vector x=p-a,y=p-b,z=b-a;
    if(dcmp(Dot(x,z))<0)return Len(x);//P距离A更近
    if(dcmp(Dot(y,z))>0)return Len(y);//P距离B更近
    return Abs(Cro(x,z)/Len(z));//面积除以底边长
}

点到直线AB距离

到直线的距离就不用判断端点了

1
2
3
4
5
6
inline bool operator==(Point a,Point b){return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);}//判断点的相等
inline LD dis_PL_(Point p,Point a,Point b){//点P到线段AB距离
    if(a==b)return Len(p-a);//AB重合
    Vector x=p-a,y=p-b,z=b-a;
    return Abs(Cro(x,z)/Len(z));//面积除以底边长
}

返回点到直线AB的垂足

求出垂足到a，a到b的长度的比例，然后加到a点上即可。

1
2
3
4
5
inline Point Foot_PL(Point p,Point a,Point b){//点P到直线AB的垂足
    Vector x=p-a,y=p-b,z=b-a;
    LD len1=Dot(x,z)/Len(z),len2=-1.0*Dot(y,z)/Len(z);//分别计算AP,BP在AB,BA上的投影
    return a+z*(len1/(len1+len2));//点A加上向量AF
}

返回点P关于AB的对称点

将点p到垂足的距离延长一倍即可

1
2
3
inline Point Symmetry_PL(Point p,Point a,Point b){//点P关于直线AB的对称点
    return p+(Foot_PL(p,a,b)-p)*2;
}

线与线

求两直线交点

1
2
3
4
inline Point cross_LL(Point a,Point b,Point c,Point d){//两直线AB,CD的交点
    Vector x=b-a,y=d-c,z=a-c;
    return a+x*(Cro(y,z)/Cro(x,y));//点A加上向量AF
}

判断直线AB与线段CD是否相交

实际上就是看直线与线段所在的直线交点是否在线段上

1
2
3
inline bool pan_cross_L_L(Point a,Point b,Point c,Point d){//判断直线AB与线段CD是否相交
    return pan_PL(cross_LL(a,b,c,d),c,d);
}

判断两线段是否相交

若是端点在另一条直线上直接返回true，不在的话根据叉积的性质来判定，一个线段两端点端点和另一个线段的两端点组成的向量叉积假如符号相同的话就代表在同一侧.

1
2
3
4
5
6
inline bool pan_cross_LL(Point a,Point b,Point c,Point d){//判断两线段AB,CD是否相交
    if(pan_PL(a,c,d)||pan_PL(b,c,d)||pan_PL(c,a,b)||pan_PL(d,a,b))return true;
    LD c1=Cro(b-a,c-a),c2=Cro(b-a,d-a);
    LD d1=Cro(d-c,a-c),d2=Cro(d-c,b-c);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(d1)*dcmp(d2)<0;
}

点与多边形

判断点是否在任意多边形内

使用射线法，时间复杂度为O(n)。
射线法（Ray casting algorithm）是一种判断点是否在多边形内部的一种简单方法。即从该点做一条射线，计算它跟多边形边界的交点个数，如果交点个数为奇数，那么点在多边形内部，否则点在多边形外部。
附上简单易懂的证明以及细节

1
2
3
4
5
6
7
8
9
10
inline int PIP(Point *P,int n,Point a){//射线法判断点A是否在任意多边形P[]以内
    int cnt=0;LD tmp;
    for(int i=1;i<=n;i++){
        int j=i<n?i+1:1;
        if(pan_PL(a,P[i],P[j]))return 2;//点在多边形上
        if(a.y>=min(P[i].y,P[j].y)&&a.y<max(P[i].y,P[j].y))//纵坐标在该线段两端点之间
            tmp=P[i].x+(a.y-P[i].y)/(P[j].y-P[i].y)*(P[j].x-P[i].x),cnt+=dcmp(tmp-a.x)>0;//交点在A右方
    }
    return cnt&1;//返回2表明在多边形边上，返回1代表在多边形内，返回0代表在多边形外
}

判断点是否在凸多边形内

使用二分法，时间复杂度为O(logn)
首先选择多边形其中一个点为起点，连接其它点作射线，将整个凸多边形分为多个小三角形。然后在射线围成的区域里判断点到底是在左半边还是右半边，最后确定的区域肯定会在一个三角形内。之后判断这个点是否在这个三角形内即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
inline int judge(Point a,Point L,Point R){//判断AL是否在AR右边
    return dcmp(Cro(L-a,R-a))>0;//必须严格以内
}
inline int PIP_(Point *P,Re n,Point a){//二分法判断点A是否在凸多边形P以内
    //点按逆时针给出
    if(judge(P[1],a,P[2])||judge(P[1],P[n],a))return 0;//在P[1_2]或P[1_n]外
    if(pan_PL(a,P[1],P[2])||pan_PL(a,P[1],P[n]))return 2;//在P[1_2]或P[1_n]上
    Re l=2,r=n-1;
    while(l<r){//二分找到一个位置pos使得P[1]_A在P[1_pos],P[1_(pos+1)]之间
        Re mid=l+r+1>>1;
        if(judge(P[1],P[mid],a))l=mid;
        else r=mid-1;
    }
    if(judge(P[l],a,P[l+1]))return 0;//在P[pos_(pos+1)]外
    if(pan_PL(a,P[l],P[l+1]))return 2;//在P[pos_(pos+1)]上
    return 1;
}

线与多边形

判断线段是否在任意多边形内

用射线法判断两端点是否在多边形内，之后再判断线段和多边形的其他线段是否有交点。

判断线段是否在任意多边形内

只需要判断两端点都在多边形内即可，不过再返回都为2时(及线段两端点都在多边形上)，需要进行特判看是否线段在多边形上。

多边形相关

判断任意两个多边形是否相离

属于不同多边形的任意两边都不相交 且一个多边形上的任意点都不被另一个多边形所包含。

1
2
3
4
5
6
7
8
9
10
11
inline int judge_PP(Point *A,int n,Point *B,int m){//【判断多边形A与多边形B是否相离】 
    for(int i1=1;i1<=n;++i1){
        int j1=i1<n?i1+1:1;
        for(int i2=1;i2<=m;++i2){
            int j2=i2<m?i2+1:1;
            if(pan_cross_LL(A[i1],A[j1],B[i2],B[j2]))return 0;//两线段相交
            if(PIP(B,m,A[i1])||PIP(A,n,B[i2]))return 0;//点包含在内
        }
    }
    return 1;
}

求任意多边形面积

使用鞋带定理(Shoelace formula) ，用来求任意多边形面积,这里是证明

1
2
3
4
5
inline LD P_Area(Point *P,Re n){//鞋带定理求任意多边形P的面积
    LD S=0;
    for(int i=1;i<=n;++i)S+=Cro(P[i],P[i<n?i+1:1]);
    return S/2.0;
}

求圆的面积并

首先得学习一下自适应辛普森法，这是一种用二次函数来逼近曲线的一种方法(有点像泰勒，怪)，不过肯定是有精度误差的，因此需要将积分区域减少以减少误差。

自适应辛普森法
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
inline LD f(LD x)
{
    /*函数表达式*/
}
inline LD simpson(LD l,LD r)
{
    LD mid=(l+r)/2.0;
    return (f(l)+4.0*f(mid)+f(r))*(r-l)/6.0;
} 
inline LD solve(LD l,LD r,LD eps)
{
    LD mid=(l+r)/2.0;
    LD s=simpson(l,r),sl=simpson(l,mid),sr=simpson(mid,r);
    if(dcmp(sl+sr-s)) return (sl+sr+(sl+sr-s));//注意，这里对精度要求较高，可能需要更改dcmp的精度范围
    return solve(l,mid,eps/2.0)+solve(mid,r,eps/2.0);
}

(哇，博客园一直在审核，好多东西都搜不到，待更。)

扫描线法

扫描线好像求⚪的面积并很难搞，只不过和自适应辛普森差不多都是求面积的所以先放到一块，交给你了！后来的林泽音。(　ﾟ 3ﾟ)

三角形面积并

跟⚪的面积并大概差不多？只不过都是直线应该用扫描线更好一点？(・ω・)

圆相关

给点求圆

拥有初中学历的人都知道，三点才能确定一个圆。
好像可以暴力带入公式啥的，不过太麻烦了就先旷过。
拥有初中学历的人还知道，圆的圆心在任意两条不平行弦的垂直平分线的交点上。

1
2
3
4
5
6
struct Circle{Point O;LD r;Circle(Point P,LD R=0){O=P,r=R;}};
inline Circle getcircle(Point A,Point B,Point C){//三点确定一圆向量垂心法
    Point P1=(A+B)*0.5,P2=(A+C)*0.5;
    Point O=cross_LL(P1,P1+Normal(B-A),P2,P2+Normal(C-A));
    return Circle(O,Len(A-O));
}

最小覆盖圆算法

使用随机增量法，把点全部打乱然后用同样的方法添加
用随机化降低不完全枚举造成有影响后果的概率以达到降低复杂度的目的. . .印象中上一个这个干的算法好像是Pollard_Rho因数分解。
具体说明在这里

1
2
3
4
5
6
7
8
9
10
11
12
13
14
struct Circle{Point O;LD r;Circle(Point P,LD R=0){O=P,r=R;}};
inline int PIC(Circle C,Point a){return dcmp(Len(a-C.O)-C.r)<=0;}//判断点A是否在圆C内
inline void Random(Point *P,Re n){for(Re i=1;i<=n;++i)swap(P[i],P[rand()%n+1]);}//随机一个排列
inline Circle MinCover(Point *P,Re n){//求点集P的最小覆盖圆
    Random(P,n);Circle C=Circle(P[1],0);
    for(int i=2;i<=n;++i)if(!PIC(C,P[i])){
        C=Circle(P[i],0);
        for(int j=1;j<i;++j)if(!PIC(C,P[j])){
            C.O=(P[i]+P[j])*0.5,C.r=Len(P[j]-C.O);
            for(int k=1;k<j;++k)if(!PIC(C,P[k]))C=getcircle(P[i],P[j],P[k]);
        }
    }
    return C;
}

凸包相关

给点集求凸包

我学习的博客上说是Graham扫描法，但是他实际上给的是Andrew算法的代码…
感觉Andrew要比Graham要简单，就只介绍下Andrew吧。
首先先按按照x优先的顺序排序（坐标从小到大），之后易知第一个点和最后一个点肯定是在凸包上。求凸包时需要扫描一遍所有点，遇到比栈顶向量右边的就替换栈顶向量，左边的直接入栈。 扫描到头会发现. . .这种算法只会算出一半的凸包，呢就反过来再来一遍，也就是第一遍算下凸包，第二遍算上凸包。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
inline bool cmp1(Vector a,Vector b){return a.x==b.x?a.y<b.y:a.x<b.x;};//按坐标排序
inline int ConvexHull(Point *P,int n,Point *cp){//cp存的是凸包元素
    sort(P+1,P+n+1,cmp1);
    int t=0;
    for(int i=1;i<=n;++i){
        while(t>1&&dcmp(Cro(cp[t]-cp[t-1],P[i]-cp[t-1]))<=0)--t;//假如P[i]此时在栈顶两点组成直线的的右边，呢么此时栈顶两点就必不可能是凸包上的点。
        cp[++t]=P[i];
    }
    int St=t;
    for(int i=n-1;i>=1;--i){
        while(t>St&&dcmp(Cro(cp[t]-cp[t-1],P[i]-cp[t-1]))<=0)--t;
        cp[++t]=P[i];
    }
    return --t;//要减一
}

旋转卡壳

旋转卡壳，本质是用凸包上的直线做平行线去找到对踵点，然后由于凸包的壕多性质都可以用对踵点来获得，所以可以用这个算法求很多东西：凸包直径(两点距离最长)，凸包的宽(两点距离最短)，凸包间的最大/小距离，凸多边形最小面积外接矩形，凸多边形最小周长外接矩形
比较细的讲解可以学习自这篇博文

求凸包直径/宽
1
2
3
4
5
6
double maxAns=Len(cp[2]-cp[1]),minAns=Len(cp[2]-cp[1]);
for(int i=1,j=3;i<=cnt;++i){
    while(dcmp(Cro(cp[i+1]-cp[i],cp[j]-cp[i])-Cro(cp[i+1]-cp[i],cp[j+1]-cp[i]))<0)j=j<cnt?j+1:1;//找到对踵点。注意是<0，如果写<=0的话可能会被两个点的数据卡掉
    maxAns=max(maxAns,max(Len(cp[j]-cp[i]),Len(cp[j]-cp[i+1])));//求凸包直径
    minAns=min(minAns,min(Len(cp[j]-cp[i]),Len(cp[j]-cp[i+1])));//求凸包的宽
}
求面积最小矩形覆盖

主要思路和基本的旋转卡壳差不多，但是要同时旋转三个点，以及定了矩形框架后用交线来
打好的代码本来是交给luogu的. . .但是一直40pts还不给看数据，怀疑是莫得spj的问题但是有spj的BZOJ也早就凉了，所以到现在还莫得过，所以先把锅推给莫得spj，自己先用着代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
double minArea=10000000.0;
	for(int i=1,j=2,z=2,k=2;i<=n;++i){
		while(dcmp(Cro(tu[i]-tu[i+1],tu[j]-tu[i+1])-Cro(tu[i]-tu[i+1],tu[j+1]-tu[i+1]))>0)j=j<n?j+1:1;
		while(dcmp(Dot(tu[i]-tu[i+1],tu[z]-tu[i+1])-Dot(tu[i]-tu[i+1],tu[z+1]-tu[i+1]))>0)z=z<n?z+1:1;//找到宽
		k=(i==1)?j:k;
		while(dcmp(Dot(tu[i+1]-tu[i],tu[k]-tu[i])-Dot(tu[i+1]-tu[i],tu[k+1]-tu[i]))>0)k=k<n?k+1:1;//找到宽
        LD dis=Len(tu[i]-tu[i+1]),kuan,gao;
        LD l=Abs(Dot(tu[z]-tu[i],tu[i+1]-tu[i])/dis);
		LD r=Abs(Dot(tu[k]-tu[i+1],tu[i]-tu[i+1])/dis);
        LD h=Abs(Cro(tu[i]-tu[j],tu[i+1]-tu[j])/dis);
		gao=Abs(h);kuan=Abs(r+l-dis);
		if(gao*kuan<minArea)
        {
			//cout<<gao<<' '<<kuan<<' '<<minArea<<endl;
            minArea=gao*kuan;
            Square[1]=cross_LL(tu[j],tu[j]+tu[i+1]-tu[i],tu[z],tu[z]+Normal_line(tu[i+1]-tu[i]));
            Square[2]=cross_LL(tu[j],tu[j]+tu[i+1]-tu[i],tu[k],tu[k]+Normal_line(tu[i+1]-tu[i]));
            Square[3]=cross_LL(tu[i],tu[i+1],tu[k],tu[k]+Normal_line(tu[i+1]-tu[i]));
            Square[4]=cross_LL(tu[i],tu[i+1],tu[z],tu[z]+Normal_line(tu[i+1]-tu[i]));
            //保证逆时针的存点
        }
	}
	printf("%.5lf\n",minArea);

求周长最小矩形覆盖

面积最小都算出来了，周长最小也是很简单的，就是每次判断一下l+r的大小就行。
代码几乎还是上面的就不搞了

求半平面交

半平面是指一条直线会把一个平面分成两个部分，其中任意一个部分都是一个半平面，通一般情况下我们需要的是直线左边的半平面。
半平面交就是若干半平面的交集，多边形各边的半平面交就是多边形的核（一个可以看到多边形的任何一个角落的区域）
具体写法在这里

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
struct Line{
    Point a,b;LD k;Line(Point A=Point(0,0),Point B=Point(0,0)){a=A,b=B,k=atan2(b.y-a.y,b.x-a.x);}
    inline bool operator<(const Line &O)const{return dcmp(k-O.k)?dcmp(k-O.k)<0:judge(O.a,O.b,a);}//
}L[N],Q[N];
inline Point cross(Line L1,Line L2){return cross_LL(L1.a,L1.b,L2.a,L2.b);}//获取直线L1,L2的交点 
inline int judge(Line L,Point a){return dcmp(Cro(a-L.a,L.b-L.a))>0;}//判断点a是否在直线L的右边
inline int halfcut(Line *L,Re n,Point *P){//【半平面交】 
    sort(L+1,L+n+1);Re m=n;n=0;
    for(Re i=1;i<=m;++i)if(i==1||dcmp(L[i].k-L[i-1].k))L[++n]=L[i];
    Re h=1,t=0;
    for(Re i=1;i<=n;++i){
        while(h<t&&judge(L[i],cross(Q[t],Q[t-1])))--t;//当队尾两个直线交点不是在直线L[i]上或者左边时就出队
        while(h<t&&judge(L[i],cross(Q[h],Q[h+1])))++h;//当队头两个直线交点不是在直线L[i]上或者左边时就出队
        Q[++t]=L[i];
    }
    while(h<t&&judge(Q[h],cross(Q[t],Q[t-1])))--t;
    while(h<t&&judge(Q[t],cross(Q[h],Q[h+1])))++h;
    n=0;
    for(Re i=h;i<=t;++i)P[++n]=cross(Q[i],Q[i<t?i+1:h]);
    return n;
}

球

球缺

球被平面截下的一部分叫做球缺。截面叫做球缺的底面，垂直于截面的直径被截后，剩下的线段长叫做球缺的高(H)。
球缺曲面部分的面积（球冠面积）：
球缺体积公式:

球的体积并

就是两个球的体积和减去两球被截面所截的球缺体积

好玩的东西

求平面最近点对

平面上有n个二维坐标点，求最近的点对。
朴素想法肯定是n方比较，但是神奇的分治保证了可以递归处理问题，使得复杂度达到O(nlognlogn)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <bits/stdc++.h>
using namespace std;
const double inf = 1e20;
const int maxn = 100005;
struct Point
{
    double x, y;
}point[maxn];
int n, mpt[maxn];
//以x为基准排序
bool cmpxy(Point a,Point b)
{
    if(a.x!=b.x)return a.x < b.x;
    return a.y<b.y;
}
bool cmpy(int a,int b){return point[a].y < point[b].y;}
double min(double a, double b){return a < b ? a : b;}
double dis(int i, int j){return sqrt((point[i].x-point[j].x)*(point[i].x-point[j].x)+(point[i].y-point[j].y)*(point[i].y-point[j].y));}
double Closest_Pair(int left, int right)
{
    double d=inf;
    if(left == right)return d;
    if(left+1==right)return dis(left, right);
    int mid=(left+right)>>1;
    double d1=Closest_Pair(left,mid);
    double d2=Closest_Pair(mid+1,right);
    d = min(d1, d2);
    int k = 0;
    //分离出宽度为d的区间
    for(int i=left;i<=right;i++)
      if(fabs(point[mid].x-point[i].x)<=d)
        mpt[k++] = i;
    sort(mpt, mpt + k, cmpy);
    //线性扫描
    for(int i=0;i<k;i++)
        for(int j=i+1;j<k&&point[mpt[j]].y-point[mpt[i]].y<d;j++)
          d=min(d,dis(mpt[i],mpt[j]));
    return d;
}

int main()
{
  while(~scanf("%d", &n) && n)
  {
    for(int i=0;i<n;i++)scanf("%lf %lf",&point[i].x,&point[i].y);
    sort(point, point + n, cmpxy);
    printf("%.2lf\n", Closest_Pair(0,n-1)/2);
  }
  return 0;
}