二维计算几何全家桶

参考文章：范神的神仙博客

前置芝士

• 一些高中数学

• 向量的叉积：向量的点积为 $a\cdot b=|a||b|\cos <a,b>$，向量的叉积为 $a\times b=|a||b|\sin <a,b>$，在平面直角坐标系中 $a(x_1,y_1),b(x_2,y_2)$，则 $a\times b=x_1{y}_2-x_2{y}_1$。容易发现叉积不满足交换律，不要写反。

存储及基本运算

我也不知道为什么到了这就突然不压行了。

struct qwq{//点 double x,y;
};
struct qaq{//线 qwq d,a,b;//d为a->b方向的单位向量 
};
struct cir{//圆 qwq p;double r;//圆心及半径 
};
inline qwq operator + (const qwq &a,const qwq &b){return {a.x+b.x,a.y+b.y};
}
inline qwq operator - (const qwq &a,const qwq &b){return {a.x-b.x,a.y-b.y};
}
inline qwq operator * (const qwq &a,const double &b){return {a.x*b,a.y*b};
}
inline qwq operator / (const qwq &a,const double &b){return {a.x/b,a.y/b};
}
inline double operator * (const qwq &a,const qwq &b){//点积 return a.x*b.x+a.y*b.y;
}
inline double operator ^ (const qwq &a,const qwq &b){//叉积 return a.x*b.y-b.x*a.y;
}
inline double len(qwq a){return sqrt(a.x*a.x+a.y*a.y);
}
inline qaq get(qwq a,qwq b){//已知两点构造直线qwq d=(b-a)/len(b-a);return {d,a,b};
}

点线基本问题

点到直线垂足

求完垂足显然就可以用初中知识求出点到直线距离，点关于直线对称点等，不一一展开。

inline qwq cz(qaq l,qwq p){return l.a+l.d*((p-l.a)*l.d);
}

点与直线位置关系

用叉积与 $0$ 的大小关系可以判断顺时针/逆时针/在线上。

如果直线是向量，那可以进一步用点积判断点在向量上/延长线上/反向延长线上。

	if(fabs(l.d^(p-l.a))<eps){if(((p-l.a)*(p-l.b))<eps) puts("ON_SEGMENT");//注意点与端点重合的情况，实质上是<=0else if((l.d*(p-l.a))<-eps) puts("ONLINE_BACK");else puts("ONLINE_FRONT");}else{if((l.d^(p-l.a))>eps) puts("COUNTER_CLOCKWISE");else puts("CLOCKWISE");}

直线与直线位置关系

共线用叉积判断，垂直用点积判断。

	if(fabs(l1.d^l2.d)<eps) puts("2");//共线/平行else if(fabs(l1.d*l2.d)<eps) puts("1");//垂直else puts("0");//其它

判断两线段是否相交

• 快速排斥实验：定义一条线段占的区域为各边平行坐标轴且以该线段为对角线的矩形。当两个矩形不相交是两线段一定不相交。
• 跨立实验：若线段 $a$ 两端点在线段 $b$ 所在直线两侧且线段 $b$ 两端点在线段 $a$ 所在直线两侧，则通过了跨立实验。

发现通过了这两个实验之一都无法保证两线段相交，而将它们结合起来用就是充要条件了。

不要忘了跨立实验有端点在线段上的情况。（似乎写的有点丑）

inline bool check(qwq p1,qwq p2,qwq p3,qwq p4){//快速排斥实验 if(max(p1.x,p2.x)<min(p3.x,p4.x)) return 0;if(max(p3.x,p4.x)<min(p1.x,p2.x)) return 0;if(max(p1.y,p2.y)<min(p3.y,p4.y)) return 0;if(max(p3.y,p4.y)<min(p1.y,p2.y)) return 0;return 1;
}
inline int side(qwq p,qaq l){//跨立实验 double tmp=l.d^(p-l.a);if(fabs(tmp)<eps) return 0;if(tmp<-eps) return 1;return 2;
}
inline void solve(){l1=get(p1,p2),l2=get(p3,p4);if(!check(p1,p2,p3,p4)) puts("0");else{int a1=side(p1,l2),a2=side(p2,l2),a3=side(p3,l1),a4=side(p4,l1);if((a1!=a2||!a1||!a2)&&(a3!=a4||!a3||!a4)) puts("1");else puts("0"); }
}

求两直线交点

设交点为 $p$，则代入第一个直线的方程式得 $p=a_1+k\cdot d_1$，且同时 $p$ 要满足在第二条直线上，即 $(p-a_2)\times d_2=0$，联立一下：$(a_1-a_2+k\times d_1)\times d_2=0$，叉积拆开即为 $k=\frac{(a_2-a_1)\times d_2}{d_1\times d2}$。

inline qwq jd(qaq l1,qaq l2){double k=((l2.a-l1.a)^l2.d)/(l1.d^l2.d);return l1.a+(l1.d*k);
} 

求两线段距离

特殊情况是相交，答案为 $0$；其余情况就是端点到另一线段距离了，我是用垂足在不在线段上来判的，似乎又写丑了。

	l1=get(p1,p2),l2=get(p3,p4);if(xj()){puts("0.0000000000");continue;}qwq p=cz(p1,l2);double ans=min(dis(p1,p3),min(dis(p1,p4),min(dis(p2,p3),dis(p2,p4)))); if(on(p,l2)) ans=min(ans,dis(p,p1));p=cz(p2,l2);if(on(p,l2)) ans=min(ans,dis(p,p2));p=cz(p3,l1);if(on(p,l1)) ans=min(ans,dis(p,p3));p=cz(p4,l1);if(on(p,l1)) ans=min(ans,dis(p,p4));printf("%.10lf\n",ans);

公共投影问题

给定 $n$ 条线段，求是否存在一条直线使得这些线段在这条直线上的投影有公共点。

问题转化：假设存在，那么过公共点作直线垂线一定与每条线段都相交。而这条垂线一定可以通过平移使得它经过所有端点中的两个。于是枚举两个端点判断即可。

直线的旋转

向量 $(x,y)$ 逆时针旋转 $\theta$ 度得到的向量为 $(x\cos \theta -y\sin \theta ,y\cos \theta +x\sin \theta )$，详见数学必修三。

inline qwq rotate(qwq a,double t){double si=sin(t),co=cos(t);return {a.x*co-a.y*si,a.y*co+a.x*si};
}

多边形问题

下面默认所有顶点都按逆时针顺序存储。

求多边形面积

任选一个点，把多边形拆成若干三角形即可根据叉积的几何意义计算：$S=\frac{1}{2}|\sum _{i=1}^n{p}_i\times p_{i\%n+1}|$，显然任意点直接取原点是最方便的。

判断凹/凸多边形

忽略相邻三点共线情况，凸多边形即相邻两边所成向量之间的叉积正负性始终相同。

	p[0]=p[n],p[n+1]=p[1];int pre=0;rep(i,1,n){double tmp=(p[i+1]-p[i])^(p[i]-p[i-1]);if(fabs(tmp)<eps) continue;int now=(tmp>eps)?1:-1;if(!pre) pre=now;else if(pre!=now) return puts("0"),0;}

判断一个点是否在多边形内

这里是回转数算法。

把待判断的点 $a$ 与多边形每一个顶点依次相连，相邻顶点之间会形成若干有方向的夹角。画图得知当且仅当这些夹角和为 $0$ 时点在多边形外部。

注意 $eps$ 不要过小，容易炸。

inline bool on(qwq p,qaq l){if(p==l.a||p==l.b) return 1;double tmp=(p-l.a)^l.d;if(fabs(tmp)>eps) return 0;return p.x>=min(l.a.x,l.b.x)&&p.x<=max(l.a.x,l.b.x)&&p.y>=min(l.a.y,l.b.y)&&p.y<=max(l.a.y,l.b.y);
}
inline double jiao(qwq a,qwq b){double ans=acos(min(1.0,(a^b)/(len(a)*len(b))));//取min是为了避免一些精度误差带来的玄学问题 return ans;
}
inline bool pd(qaq l,qwq p){//判断正角还是负角return (l.d^(p-l.a))>eps;
}
inline void solve(){double now=0;rep(i,1,n){if(on(a,get(p[i],p[i+1]))) return puts("1"),void();//在多边形上 if(pd(get(a,p[i]),p[i+1])) now+=jiao(a-p[i],a-p[i+1]);else now-=jiao(a-p[i],a-p[i+1]);}if(fabs(now)<eps) puts("0");//在多边形外 else puts("2");
}

二维凸包

定义：包含给定的所有点的周长最小的凸多边形。

求凸包可以用 Andrew 算法在 $O(n\log n)$ 的时间内求解。

先把所有点以 $x$ 为第一关键字，$y$ 为第二关键字升序排序。那么排完序的第一个点和最后一个点一定在凸包上，且以它们为分界线可以将凸包分为上下两部分，而两部分转弯的顺逆时针情况都是相同的。所以可以用单调栈维护，先求出下凸壳，再在没使用过的点中求出上凸壳。

实现的时候注意因为第一个点要用来更新上凸壳，所以开始不能记 $vis[1]=1$；最后栈内一定重复加了一个 $1$，所以要 $top--$。

inline void solve(){sort(a+1,a+n+1);s[++top]=1;rep(i,2,n){while(top>1&&((a[s[top]]-a[s[top-1]])^(a[i]-a[s[top]]))<0) vis[s[top]]=0,--top;vis[i]=1,s[++top]=i;}int pos=top;for(int i=n;i;--i) if(!vis[i]){while(top>pos&&((a[s[top]]-a[s[top-1]])^(a[i]-a[s[top]]))<0)vis[s[top]]=0,--top;vis[i]=1,s[++top]=i;}--top;
}

旋转卡壳

常用于求平面上最远点对。

首先显然最远点对的两个点一定都在凸包上，而逆时针遍历凸包的每一条边的过程中，离当前边所在直线距离最远的点也一定沿着逆时针方向旋转，于是维护这个最远点的位置，用这个点到当前边两端点的距离更新答案即可。

inline double dis(qwq a,qwq b){double x=a.x-b.x,y=a.y-b.y;return x*x+y*y;
}
inline double disl(qwq p,qaq l){return fabs(l.d^(p-l.a));
}
inline double kk(){p[tot+1]=p[1],p[0]=p[tot],nxt[tot]=1;rep(i,1,tot-1) nxt[i]=i+1;int now=2;double ans=0;rep(i,1,tot){qaq l=get(p[i],p[i+1]);while(disl(p[nxt[now]],l)>disl(p[now],l)) now=nxt[now];ans=max(ans,max(dis(p[i],p[now]),dis(p[i+1],p[now])));}return sqrt(ans);
}

半平面交

顾名思义，求解多个凸多边形的面积交问题。

先将所有边进行极角排序，即用 $atan2(y,x)$ 求出 $\tan \alpha =\frac{y}{x},\alpha \in (-\pi ,\pi ]$ 的角 $\alpha$ ，不同的角按角的大小排序，相同的将位置更靠里的放在前面。

然后用一个 $deque$ 维护对答案可能有贡献的边及相邻边的交点。具体地，每加入一条边，先不断从队尾弹出已经不合法的交点，再不断从队首弹出已经不合法的交点，最后插入当前边已经当前边和之前的队尾的交点。

插入进行到最后时，可能队尾的一些点被队首卡掉了，而队首的点一直被后来加入的队尾们约束着，不可能被卡掉，所以要不断弹出队尾直至合法。

最后别忘了加入一下队尾和队首的交点。

构成一个多边形至少要三条边，若最终队列长度小于 $3$，说明面积交为 $0$。

inline bool in(qwq p,qaq l){return (l.d^(p-l.a))>eps;}
inline bool out(qwq p,qaq l){return (l.d^(p-l.a))<-eps;}
inline bool cmp(qaq l1,qaq l2){double p1=atan2(l1.d.y,l1.d.x),p2=atan2(l2.d.y,l2.d.x);if(fabs(p1-p2)>eps) return p1<p2;return in(l1.a,l2);
}
inline qwq jiao(qaq l1,qaq l2){double k=((l2.a-l1.a)^l2.d)/(l1.d^l2.d);return l1.a+l1.d*k;
}
inline void half(){sort(a+1,a+tot+1,cmp);rep(i,1,tot){if(l<=r&&fabs(q[r].d^a[i].d)<eps) continue;//判掉无用的平行边while(l<r&&out(jd[r],a[i])) --r;//当队列长为2且唯一交点不合法时先弹出队首会出错，必须先弹队尾while(l<r&&out(jd[l+1],a[i])) ++l;//交点是从队首的下一个开始产生的q[++r]=a[i];if(l<r) jd[r]=jiao(a[i],q[r-1]);}while(l<r&&out(jd[r],q[l])) --r;if(l<r) jd[r+1]=jiao(q[l],q[r]),++r;n=r-l;rep(i,1,n) p[i]=jd[i+l];p[n+1]=p[1];
}

圆相关问题

三角形内切圆

求出两条角平分线，交点即为圆心。半径即为点到直线距离。

inline qaq pf(qaq l1,qaq l2){//求角平分线qwq d=(l1.d+l2.d)/len(l1.d+l2.d);return {d,l1.a,l1.a+d};
}
inline cir nqy(qwq p1,qwq p2,qwq p3){cir ans;qaq l1=pf(get(p1,p2),get(p1,p3));qaq l2=pf(get(p2,p1),get(p2,p3));ans.p=jd(l1,l2),ans.r=dis(ans.p,get(p1,p2));return ans; 
}

三角形外接圆

求两条中垂线交点即可。

inline qaq zcx(qaq l){//中垂线qwq mid={(l.a.x+l.b.x)/2,(l.a.y+l.b.y)/2};return {{l.d.y,-l.d.x},mid};
}
inline cir wjy(qwq p1,qwq p2,qwq p3){qaq l1=zcx(get(p1,p2)),l2=zcx(get(p2,p3));qwq p=jd(l1,l2);return {p,dis(p,p1)};
}

圆和直线交点

在保证有交点的前提下，求出圆心到直线垂足后勾股定理即可。

inline void query(qaq l){vector<qwq> ans;qwq cz=l.d*(l.d*(c.p-l.a))+l.a;double d1=len(cz-c.p),d2=sqrt(c.r*c.r-d1*d1);ans.push_back(cz+l.d*d2),ans.push_back(cz-l.d*d2);sort(ans.begin(),ans.end());for(qwq tmp:ans) printf("%.8lf %.8lf ",tmp.x,tmp.y);putchar('\n');
}

圆与圆交点

发现以两圆心和其中一个交点为顶点的三角形三边均已知，就可以用余弦定理求出其中一个角，旋转经过两圆心的直线即可。

inline void query(cir c1,cir c2){vector<qwq> ans;double dis=len(c1.p-c2.p);qaq l=get(c1.p,c2.p);double t=acos((c1.r*c1.r-c2.r*c2.r+dis*dis)/(2*c1.r*dis));qwq d1=rotate(l.d,t),d2=rotate(l.d,2*pi-t);ans.push_back(c1.p+d1*c1.r),ans.push_back(c1.p+d2*c1.r);sort(ans.begin(),ans.end());for(qwq tmp:ans) printf("%.8lf %.8lf ",tmp.x,tmp.y);putchar('\n');
}

过圆外一点作圆的切线

圆外点 $p$，圆心和任意切点组成的直角三角形两边已知，可以求解其中一锐角并将直线旋转。

inline void query(qwq p,cir c){vector<qwq> ans;qaq l=get(p,c.p);double dis=len(p-c.p),t=asin(c.r/dis),d=sqrt(dis*dis-c.r*c.r);qwq d1=rotate(l.d,t),d2=rotate(l.d,pi*2-t);ans.push_back(p+d1*d),ans.push_back(p+d2*d);sort(ans.begin(),ans.end());for(qwq tmp:ans) printf("%.8lf %.8lf\n",tmp.x,tmp.y);
}

求两圆公切线

分内含/内切/相交/外切/外离五种情况考虑。也是利用一些简单的几何性质求解直角三角形，然后旋转得到交点。

代码求的是在圆 $c_1$ 上的交点。

inline void solve(cir c1,cir c2){int flag=0;if(c1.r<c2.r) swap(c1,c2),flag=1;double dis=len(c1.p-c2.p);if(dis+c2.r<c1.r) return;//内含qaq l=get(c1.p,c2.p);if(fabs(c1.r-c2.r-dis)<eps){//内切qwq ans=c1.p+l.d*c1.r;return printf("%.10lf %.10lf\n",ans.x,ans.y),void();}vector<qwq> ans;double t=acos((c1.r-c2.r)/dis);qwq d1=rotate(l.d,t),d2=rotate(l.d,2*pi-t);if(flag) ans.pb(c2.p+d1*c2.r),ans.pb(c2.p+d2*c2.r);else ans.pb(c1.p+d1*c1.r),ans.pb(c1.p+d2*c1.r);if(fabs(dis-c1.r-c2.r)<eps) ans.pb(c1.p+l.d*c1.r);//外切else if(dis>c1.r+c2.r){//外离if(flag) swap(c1,c2),l=get(c1.p,c2.p);double xie=dis*c1.r/(c1.r+c2.r),t=acos(c1.r/xie);d1=rotate(l.d,t),d2=rotate(l.d,2*pi-t);ans.pb(c1.p+d1*c1.r),ans.pb(c1.p+d2*c1.r);}sort(ans.begin(),ans.end());for(qwq tmp:ans) printf("%.10lf %.10lf\n",tmp.x,tmp.y);
}

扇形面积与弓形面积

根据初中公式求解即可。注意这里算的都是圆心角 $\le \pi$ 的部分的面积。

inline double angle(qwq p1,qwq p2){return acos((p1*p2)/len(p1)/len(p2));
}
inline double sshan(cir c,qwq a,qwq b){double t=angle(a-c.p,b-c.p);return t/2*c.r*c.r;
}
inline double sgong(cir c,qwq a,qwq b){return sshan(c,a,b)-fabs((b-c.p)^(a-c.p)/2);
}

圆与圆面积交

特判内含/外离之后考虑处理相交。发现当圆心角在 $\pi$ 以内的时候对答案的贡献为扇形面积减三角形面积，否则答案为扇形面积加三角形面积，但此时 $\sin \theta$ 为负，依然正确。

inline double query(cir c1,cir c2){if(c1.r<c2.r) swap(c1,c2);double dis=len(c1.p-c2.p);if(dis>=c1.r+c2.r) return 0;if(c2.r+dis<=c1.r) return pi*c2.r*c2.r;double r1=c1.r*c1.r,r2=c2.r*c2.r;double t1=2*acos((dis*dis+r1-r2)/(2*dis*c1.r));double t2=2*acos((dis*dis+r2-r1)/(2*dis*c2.r));return t1/2*r1+t2/2*r2-r1*sin(t1)/2-r2*sin(t2)/2;
}

圆与多边形面积交

把多边形拆成多个由多边形上相邻两顶点和圆心组成的三角形，根据圆心与多边形的位置关系将每部分面积交标上正负符号，然后分类：

• （当前正在讨论的多边形上的）两点全在圆内，答案为三角形面积；

• 两点全在圆外：

  • 两点所成线段与圆没有交点，答案为扇形面积；
  • 两点所成线段与圆有交点，答案为扇形面积减多余的弓形面积；

• 一点在圆内一点在圆外：答案为一个三角形面积 + 一个扇形面积。

inline qwq cirline(cir c,qaq l){qwq cz=l.d*((c.p-l.a)*l.d)+l.a;double dis=len(c.p-cz);return cz+l.d*sqrt(c.r*c.r-dis*dis);
}
inline bool check(qwq a,qwq b,cir c){return max(a.x,b.x)>=c.p.x-eps&&min(a.x,b.x)<=c.p.x+eps||max(a.y,b.y)>=c.p.y-eps&&min(a.y,b.y)<=c.p.y+eps;
}
inline double calc(cir c,qwq a,qwq b){double tmp=(a-c.p)^(b-c.p);if(fabs(tmp)<eps) return 0;int f=tmp>0?1:-1,f1=len(a-c.p)<=c.r,f2=len(b-c.p)<=c.r;if(f1&&f2) return tmp/2;qaq l=get(a,b);qwq cz=l.d*((c.p-l.a)*l.d)+l.a;double jl=len(c.p-cz);if(jl<c.r) tmp=sqrt(c.r*c.r-jl*jl);if(!f1&&!f2){qaq l1=get(c.p,a),l2=get(c.p,b);double s=sshan(c,l1.a+l1.d*c.r,l2.a+l2.d*c.r);if(jl<c.r&&check(a,b,c)) s-=sgong(c,cz-l.d*tmp,cz+l.d*tmp);return s*f;}if(!f1) swap(a,b),l=get(a,b);qwq p1=cirline(c,l),p2=cirline(c,get(c.p,b));double s=fabs((c.p-a)^(p1-a)/2)+sshan(c,p1,p2);return s*f;
}
inline double query(){a[n+1]=a[1];double ans=0;rep(i,1,n) ans+=calc(c,a[i],a[i+1]);return fabs(ans);
}

杂项

曼哈顿距离与切比雪夫距离

平面上两点 $(x_1,y_1)$ 和 $(x_2,y_2)$，它们的曼哈顿距离为 $|x_1-x_2|+|y_1-y_2|$，切比雪夫距离为 $\max (|x_1-x_2|,|y_1-y_2|)$。

• 曼哈顿距离 $\Rightarrow$ 切比雪夫距离：$(x,y)\Rightarrow (x+y,x-y)$；
• 切比雪夫距离 $\Rightarrow$ 曼哈顿距离：$(x,y)\Rightarrow (\frac{x+y}{2},\frac{x-y}{2})$。