计算几何学习笔记
_ANIG_
·
·
算法·理论
说的比较简略，不保证全对。
基本运算
二维向量
加减
(x_1,y_1)+(x_2,y_2)=(x_1+y_1,x_2+y_2)
(x_1,y_1)-(x_2,y_2)=(x_1-y_1,x_2-y_2)
对于平面上的两个点 A(x_1,y_1),B(x_2,y_2)，则向量 \overrightarrow{AB}=B-A。
数乘
k(x,y)=(kx,ky)
模长
|(x,y)|=\sqrt{x^2+y^2}
点乘
结果是数值：
(x_1,y_1)\cdot(x_2,y_2)=x_1x_2+y_1y_2
设向量 \vec a,\vec b 的夹角为 \theta，则：
若 \theta=0，则 \vec a\cdot\vec b=|a||b|。
若 0<\theta<90^\circ，则 \vec a\cdot \vec b>0。
若 \theta=90^\circ，则 \vec a\cdot\vec b=0。
若 90^\circ<\theta<180^\circ，则 \vec a\cdot\vec b<0。
若 \theta=180^\circ，则 \vec a\cdot\vec b=-|a||b|。
叉乘
结果是向量，其模长为：
|(x_1,y_1)\times(x_2,y_2)|=x_1y_2-x_2y_1
这个结果的绝对值等于 \vec a 和 \vec b 围成的平行四边形的面积。
若 \vec a\times\vec b<0，则 \vec a 在 \vec b 的逆时针方向。
若 \vec a\times\vec b=0，则 \vec a 与 \vec b 共线。
若 \vec a\times \vec b>0，则 \vec a 在 \vec b 的顺时针方向。
下文叉乘代表结果的模长。
旋转
### 两点距离
$AB=|\overrightarrow{AB}|
线
用两个点 A,B 表示有向线段 AB。
点在直线哪一侧
若 \overrightarrow{AB}\times \overrightarrow{AC}>0，则 C 在 AB 左侧。
若 \overrightarrow{AB}\times \overrightarrow{AC}=0，则 C 在 AB 上。
若 \overrightarrow{AB}\times \overrightarrow{AC}<0，则 C 在 AB 右侧。
点是否在线段上
### 两线段是否相交
先判断是否有一个线段端点在另一个线段上。
$AB$ 和 $CD$ 相交的条件为 $A,B$ 在直线 $CD$ 的不同侧且 $CD$ 在 $AB$ 的不同侧。
### 两直线交点
要求直线 $AB$ 和直线 $CD$ 的交点。
直线 $AB$ 上的点可表示为 $A+k\overrightarrow{AB}$。
则 $A+k_1\overrightarrow{AB}=C+k_2\overrightarrow{CD}$。
两边同时叉乘 $\overrightarrow{CD}$ 得：
$$A\times \overrightarrow{CD}+k_1\overrightarrow{AB}\times \overrightarrow{CD}=C\times \overrightarrow{CD}$$
解得 $k_1=\frac{\overrightarrow{AC}\times\overrightarrow{CD}}{\overrightarrow{AB}\times \overrightarrow{CD}}$。
交点即为 $A+k_1\overrightarrow{AB}$。
## 图形
### 三角形面积
面积等于对应梯形面积的一半，而梯形面积可以叉乘求出。
$$S_{\triangle ABC}=|\frac{\overrightarrow{AB}\times \overrightarrow{AC}}2|$$
### 凸多边形面积
把多边形剖成三角形，面积相加。
### 判断点是否在凸多边形内部
将这个点与凸多边形的边构成的三角形面积相加，判断是否等于凸多边形面积，单次 $O(n)$。
一个更快的做法是将所有顶点按极角排序，二分出来要查询的点位于哪个三角形，复杂度 $O(\log n)$。
## 三维向量
### 加减
$$(x_1,y_1,z_1)+(x_2,y_2,z_2)=(x_1+x_2,y_1+y_2,z_1+z_2)$$
$$(x_1,y_1,z_1)-(x_2,y_2,z_2)=(x_1-x_2,y_1-y_2,z_1-z_2)$$
### 点乘
$$(x_1,y_1,z_1)\cdot(x_2,y_2,z_2)=x_1x_2+y_1y_2+z_1z_2$$
几何意义同二维向量。
### 叉乘
叉乘结果是一个向量。
$$(x_1,y_1,z_1)\times (x_2,y_2,z_2)=(y_1z_2-y_2z_1,z_1x_2-z_2x_1,x_1y_2-x_2y_1)$$
模长为对应平行四边形面积，方向用右手定则判断。
### 平面
可以用三个不共线的点确定一个平面。
#### 法向量
设 $A,B,C$ 为平面上逆时针排列的三个点，则 $A,B,C$ 所在平面的法向量为：$\overrightarrow{AB}\times \overrightarrow{AC}$，方向垂直于这个平面向上。
#### 判断一个点在平面的哪一侧
要判断点 $D$ 在 $A,B,C$ 所在平面的哪一侧。
先求出平面的法向量 $\overrightarrow{AP}$，用点乘算出 $AP$ 和 $AD$ 的夹角是锐角还是钝角，若锐角，则在平面上方。
### 判断一个点是否在三个点确定的圆内
要判断点 $D$ 是否在经过 $A,B,C$ 的圆内，先把所有点 $(x,y)$ 变换为三维向量 $(x,y,x^2+y^2)$。
若 $D$ 在 $A,B,C$ 所在平面的下方，在 $D$ 在 $A,B,C$ 确定的圆内部。
# 极角排序
$(x,y)$ 的极角可以使用函数 ```atan2(y,x)``` 求出。
按极角排序相当于个所有点按逆时针顺序排序。
[[CQOI2017] 小Q的草稿](https://www.luogu.com.cn/problem/P3699)
以每个点为中心，统计这个点的贡献。
极角排序后，类似扫描线，每次加入一条线段，删除一条线段或查询一个点，set 维护即可。
[Bob and stages](https://www.luogu.com.cn/problem/CF852H)
枚举起点，设 $f_{k,x}$ 表示选到第 $k$ 个点，最后一个点是 $x$ 的最大面积。
但是要求凸性，所以需要记倒数第二个点是什么。
可以提前把所有两两之间的连线按极角排序，从小到大转移，这样，凸包的限制自动满足。
然后就是判断一个三角形内是否有其他点。
可以把所有点按极角排序，枚举三角形的一条边，另一条逆时针扫一圈。
# 坐标旋转
有些题，为了更方便做，可以把所有点绕原点旋转一个角度，或平移一段距离。
[[SDOI2018] 物理实验](https://www.luogu.com.cn/problem/P4605)
把直线导轨转到与 $y$ 轴重合的位置，把所有端点投到 $y$ 轴上离散化，扫描线维护一下，这样每一段区间能照到的线段是唯一的，双指针维护即可。
# 二维凸包
先找到最靠下的点，若有多个取最靠左的。
这个点一定在凸包上，以这个点为中心极角排序，若极角相同按照到这个点的距离排序。
然后用栈维护当前凸包。加入一个点，从栈顶弹出不合法的点，再把当前的点加进去。
最后，栈里面的点就是凸包顶点。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e5+5;
const double eps=1e-12,pi=acos(-1);
struct node{
double x,y;
friend node operator-(node a,node b){
return {a.x-b.x,a.y-b.y};
}
}g[N];
double jj(node x){
return atan2(x.y,x.x);
}
double dis(node a,node b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double xl(node x,node y){
auto tmp=jj(y-x);
if(tmp<0)tmp+=2*pi;
return tmp;
}
int n,q[N],tl;
signed main(){
cin>>n;
int wz=0;
for(int i=1;i<=n;i++){
cin>>g[i].x>>g[i].y;
if(!wz||g[i].y<g[wz].y||(g[i].y==g[wz].y&&g[i].x<g[wz].x))wz=i;
}
swap(g[wz],g[1]);
sort(g+2,g+n+1,[&](node a,node b){
return jj(a-g[1])<jj(b-g[1]);
});
q[++tl]=1;
for(int i=2;i<=n;i++){
while(tl>1&&xl(g[q[tl]],g[i])<xl(g[q[tl-1]],g[q[tl]])-eps)tl--;
q[++tl]=i;
}
double res=0;
for(int i=2;i<=tl;i++)res+=dis(g[q[i]],g[q[i-1]]);
res+=dis(g[q[1]],g[q[tl]]);
printf("%.2lf",res+eps);
}
```
另一种做法是以 $x$ 为第一关键字，$y$ 为第二关键字排序。
然后从左往右跑一边上面类似的做法。
这样会跑出来一个下凸壳，倒着再跑一遍即可得到完整凸包。
此做法精度误差更小。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=2e5+5;
const double eps=1e-12,pi=acos(-1);
struct node{
double x,y;
friend node operator-(node a,node b){
return {a.x-b.x,a.y-b.y};
}
}g[N];
double cro(node a,node b){
return a.x*b.y-a.y*b.x;
}
double dis(node a,node b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int n,q[N],tl;
signed main(){
cin>>n;
for(int i=1;i<=n;i++){
cin>>g[i].x>>g[i].y;
}
sort(g+1,g+n+1,[&](node a,node b){
if(a.x==b.x)return a.y<b.y;
return a.x<b.x;
});
for(int i=1;i<=n;i++){
while(tl>1&&cro(g[q[tl-1]]-g[q[tl]],g[i]-g[q[tl]])>eps)tl--;
q[++tl]=i;
}
for(int i=n;i>=1;i--){
while(tl>1&&cro(g[q[tl-1]]-g[q[tl]],g[i]-g[q[tl]])>eps)tl--;
q[++tl]=i;
}
double res=0;
for(int i=2;i<=tl;i++)res+=dis(g[q[i]],g[q[i-1]]);
printf("%.2lf",res+eps);
}
```
[[ICPC 2021 Nanjing R] Paimon Polygon](https://www.luogu.com.cn/problem/P9845)
可以发现，划分出来的两个凸包要么包含，要么相离。
若包含，则大凸包就是整个点集的凸包，直接求出来整个点集的凸包，判断剩余的点集是否合法即可。
若相离，则极角排序后，两个凸包包含的点对于两个区间。
对于每个点，判断这个点能否作为左端点，右端点和中间点，均摊下来可以 $O(n)$ 统计答案。
# 闵可夫斯基和
对于两个凸包 $\{(x,y)|(x,y)\in A\}$ 和 $\{(x,y)|(x,y)\in B\}$，闵可夫斯基和的结果为 $\{(x_1+x_2,y_1+y_2)|(x_1,y_1)\in A,(x_2,y_2)\in B\}$ 的凸包。
具体求法是把凸包的边变成向量，按极角排序再合并起来。
归并排序可以做到 $O(n)$。
[[JSOI2018] 战争](https://www.luogu.com.cn/problem/P4557)
问题等价于求出点集 $A$ 和 $-B$ 的闵可夫斯基和，然后多次询问一个点是否在结果凸包内。
用上文提到的方法二分查询，总复杂度 $O(n\log n)$。
# 半平面交
给定若干条有向直线，每条直线表示保留这条直线左侧的半平面。求所有半平面的交。
把直线按极角排序，增量构造。
维护双端队列，表示当前有用的直线。
若队尾的两条直线的交点位于当前处理的直线的右侧，则最后一条直线就没用了。
由于最后会成为一个凸包，所以要给队头做同样的判断。
最后还要用队头更新队尾。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e4+5;
const double eps=1e-12;
struct node{
double x,y;
friend node operator-(node a,node b){
return {a.x-b.x,a.y-b.y};
}
friend node operator*(double a,node b){
return {a*b.x,a*b.y};
}
friend node operator+(node a,node b){
return {a.x+b.x,a.y+b.y};
}
};
double jj(node x){
return atan2(x.y,x.x);
}
bool chk(double x){
return abs(x)<=eps;
}
double cro(node a,node b){
return a.x*b.y-a.y*b.x;
}
struct vec{
node a,b;
friend bool operator<(vec a,vec b){
if(chk(jj(a.b-a.a)-jj(b.b-b.a)))return cro(a.a-b.a,b.b-b.a)<-eps;
return jj(a.b-a.a)<jj(b.b-b.a);
}
};
node jd(vec a,vec b){
double k=cro(a.a-b.a,b.b-b.a)/cro(b.b-b.a,a.b-a.a);
return a.a+k*(a.b-a.a);
}
bool chk(vec a,vec b,vec c){
auto x=jd(a,b);
return cro(x-c.a,c.b-c.a)>-eps;
}
double S(node a,node b,node c){
return cro(b-a,c-a)/2;
}
vector<vec>q;
vec qs[N];
node g[N];
int n,hd,tl;
signed main(){
cin>>n;
for(int i=1;i<=n;i++){
vector<node>jl;
int m;
cin>>m;
while(m--){
double x,y;
cin>>x>>y;
jl.push_back({x,y});
}
for(int i=1;i<jl.size();i++)q.push_back({jl[i-1],jl[i]});
q.push_back({jl.back(),jl[0]});
}
sort(q.begin(),q.end());
vector<vec>tmp;
double lst=1e9;
for(auto c:q){
if(chk(jj(c.b-c.a)-lst))continue;
lst=jj(c.b-c.a);
tmp.push_back(c);
}
q=tmp;
hd=1;tl=0;
for(auto c:q){
while(hd<tl&&chk(qs[tl-1],qs[tl],c))tl--;
while(hd<tl&&chk(qs[hd+1],qs[hd],c))hd++;
qs[++tl]=c;
}
while(hd<tl&&chk(qs[tl-1],qs[tl],qs[hd]))tl--;
int idx=0;
for(int i=hd+1;i<=tl;i++)g[++idx]=jd(qs[i],qs[i-1]);
g[++idx]=jd(qs[hd],qs[tl]);
double res=0;
for(int i=3;i<=idx;i++)res+=S(g[1],g[i-1],g[i]);
printf("%.3lf",res+eps);
}
```
[[SCOI2015] 小凸想跑步](https://www.luogu.com.cn/problem/P4250)
[[ZJOI2019] 浙江省选](https://www.luogu.com.cn/problem/P5328)
# 反演
> 反演变换适用于题目中存在多个圆/直线之间的相切关系的情况。利用反演变换的性质，在反演空间求解问题，可以大幅简化计算。
对于一个反演半径 $r$ 和反演中心 $O$，$P$ 的反演点 $P'$ 满足 $O,P,P'$ 共线，且 $OP\times OP'=r^2$。
有两个比较有用的性质：
- 过点 $O$ 的圆反演图形为一条直线。
- 在原图相切的图形反演后也相切。
[Drawing Circles is Fun](https://www.luogu.com.cn/problem/CF372E)
以 $O$ 为反演中心，任意长度为反演半径进行反演变换。
原图中两个经过 $O$ 的相切的圆反演后就是两条平行线。
于是原图的限制变成了对边平行，也就是任意两条线段作为对角线都能得到平行四边形。于是这个问题就变成了容易的问题。
# 半平面数据结构
给定 $n$ 条直线，$q$ 次询问，每次询问一个点下方有几条直线。
把直线每 $T$ 个分成一块，分开处理。
对于一组直线，枚举直线两两之间的交点。
扫描线的过程中维护直线的相对顺序，遇到交点就调整，二分查询。
总复杂度 $O((T^2+q\log T)\frac nT)$。
取 $T=\sqrt q$，总复杂度 $O(n\sqrt q)$。
[【UNR #4】己酸集合](https://uoj.ac/problem/553)
# 剖分
## 梯形剖分
多边形不是凸的了。
从前往后扫描线，处理出来每条边的方向。
剖出来的竖直的先必须满足向左的边在上，向右的边在下，这样就能避免剖到多边形外面。
[Spider](https://www.luogu.com.cn/problem/CF223D)
## Delaunay 三角剖分
把平面上若干个点剖成三角形。
先把所有点随机扰动，使得不存在四点共圆或三点共线。
把所有点按 $x$ 排序，分治构造。
假设现在已经构造出了左半边和右半边。
先找到最靠下的连接左右两边的边。具体实现可以求出来这个点集的下凸壳，找到连接左右两半的那条。
假设找到的点为 $A,B$。
找到右侧与 $B$ 相连的点 $C$，使得：
- $C$ 在 $\overrightarrow{BA}$ 的顺时针方向。
- 不存在与 $B$ 相连的其他点 $D$，使得 $D$ 在 $A,B,C$ 构成的圆内部。
左半边的点同理。
通过打擂台的方式可以找到唯一的一个点。
假设这个点位于右半边，为 $C$。
连接 $AC$，删掉与 $B$ 相连的与 $AC$ 有交的边。
以 $AC$ 为新的最底部的边重复这个过程。
总复杂度 $O(n\log n)$，我实现的似乎并不是这个复杂度，但在 $10^5$ 的数据范围下能跑。
[平面欧几里得最小生成树](https://www.luogu.com.cn/problem/P6362)
先跑出来三角剖分，只保留三角剖分的边跑最小生成树即可。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e5+5;
const double eps=1e-9;
struct vec{
double x,y,z;
friend vec operator-(vec a,vec b){
return {a.x-b.x,a.y-b.y,a.z-b.z};
}
};
vec cro(vec a,vec b){
return {a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};
}
double dot(vec a,vec b){
return a.x*b.x+a.y*b.y+a.z*b.z;
}
bool chk(vec a,vec b,vec c,vec x){
return dot(cro(b-a,c-a),x-a)<=0;
}
struct node{
double x,y;
void rnd(){
x+=((int)rand()-RAND_MAX/2)*1.0/RAND_MAX*eps;
y+=((int)rand()-RAND_MAX/2)*1.0/RAND_MAX*eps;
}
friend bool operator<(node a,node b){
return a.x<b.x;
}
friend node operator-(node a,node b){
return {a.x-b.x,a.y-b.y};
}
}p[N];
vec to(node x){
return {x.x,x.y,x.x*x.x+x.y*x.y};
}
bool in(node a,node b,node c,node x){
return chk(to(a),to(b),to(c),to(x));
}
double mc(node x){
return x.x*x.x+x.y*x.y;
}
double Dis(node a,node b){
return mc(b-a);
}
double cro(node a,node b){
return a.x*b.y-a.y*b.x;
}
bool in(node a,node b,node x){
if(cro(a-x,b-x))return 0;
if(Dis(a,x)>Dis(a,b))return 0;
if(Dis(b,x)>Dis(a,b))return 0;
return 1;
}
bool xj(node a,node b,node c,node d){
if(in(c,d,a)||in(c,d,b)||in(a,b,c)||in(a,b,d))return 0;
return ((cro(a-c,d-c)>0)^(cro(b-c,d-c)>0))&((cro(c-a,b-a)>0)^(cro(d-a,b-a)>0));
}
unordered_set<int>e[N];
int n,fa[N],qs[N],tl;
void solve(int l,int r){
if(r-l+1<=3){
for(int i=l;i<=r;i++){
for(int j=l;j<=r;j++){
if(i==j)continue;
e[i].insert(j);
}
}
return;
}
int mid=l+r>>1;
solve(l,mid);solve(mid+1,r);
int a=0,b=0,s1=0,s2=0;
tl=0;
for(int i=l;i<=r;i++){
while(tl>1&&cro(p[qs[tl-1]]-p[qs[tl]],p[i]-p[qs[tl]])>=0)tl--;
qs[++tl]=i;
}
for(int i=1;i<tl;i++){
if(qs[i]<=mid&&qs[i+1]>mid){
a=qs[i];b=qs[i+1];
break;
}
}
auto fd1=[&](){
int x=0;
for(auto c:e[a]){
if(cro(p[a]-p[b],p[c]-p[b])>=0)continue;
if(!x){
x=c;
continue;
}
// cout<<c<<" "<<x<<" "<<in(p[a],p[b],p[x],p[c])<<" "<<in(p[a],p[b],p[c],p[x])<<endl;
if(in(p[a],p[b],p[x],p[c]))x=c;
}
return x;
};
auto fd2=[&](){
int x=0;
for(auto c:e[b]){
if(cro(p[b]-p[a],p[c]-p[a])<=0)continue;
if(!x){
x=c;
continue;
}
if(in(p[a],p[b],p[x],p[c]))x=c;
}
return x;
};
while(1){
e[a].insert(b);e[b].insert(a);
s1=fd1();s2=fd2();
if(!s1&&!s2)break;
if(!s2||(s1&&in(p[a],p[b],p[s2],p[s1]))){
vector<int>jl;
for(auto c:e[a]){
if(xj(p[a],p[c],p[b],p[s1]))jl.push_back(c);
}
for(auto c:jl){
e[a].erase(e[a].find(c));
e[c].erase(e[c].find(a));
}
a=s1;
s1=0;
}else{
vector<int>jl;
for(auto c:e[b]){
if(xj(p[b],p[c],p[a],p[s2]))jl.push_back(c);
}
for(auto c:jl){
e[b].erase(e[b].find(c));
e[c].erase(e[c].find(b));
}
b=s2;
s2=0;
}
}
}
struct edge{
int x,y;double z;
friend bool operator<(edge a,edge b){
return a.z<b.z;
}
};
vector<edge>q;
int fd(int x){
if(x==fa[x])return x;
return fa[x]=fd(fa[x]);
}
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
cin>>n;
for(int i=1;i<=n;i++)cin>>p[i].x>>p[i].y,p[i].rnd();
sort(p+1,p+n+1);
solve(1,n);
for(int i=1;i<=n;i++){
for(auto c:e[i])q.push_back({i,c,sqrt(Dis(p[i],p[c]))});
}
for(int i=1;i<=n;i++)fa[i]=i;
double res=0;
sort(q.begin(),q.end());
for(auto[x,y,z]:q){
x=fd(x);y=fd(y);
if(x==y)continue;
fa[x]=y;
res+=z;
}
printf("%.6lf",res);
}
```
# 三维凸包
增量构造。
先随机扰动，这样凸包上的每个面都是三角形。
直接存下来所有三角形。
每次加入一个点，把这个点能看到的面删掉，再把边界处的边和这个点连起来，构成新的凸包。
可以直接枚举所有面，只有边界上的边会有恰好一个邻面被删了。
总复杂度 $O(n^2)$，代码很好写。
[【模板】三维凸包](https://www.luogu.com.cn/problem/P4724)
```cpp
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=2005;
const double eps=1e-8;
struct node{
double x,y,z;
void rnd(){
x+=(rand()-RAND_MAX/2)*1.0/RAND_MAX*eps;
y+=(rand()-RAND_MAX/2)*1.0/RAND_MAX*eps;
z+=(rand()-RAND_MAX/2)*1.0/RAND_MAX*eps;
}
friend node operator-(node a,node b){
return {a.x-b.x,a.y-b.y,a.z-b.z};
}
}p[N];
double mc(node x){
return sqrt(x.x*x.x+x.y*x.y+x.z*x.z);
}
double dot(node a,node b){
return a.x*b.x+a.y*b.y+a.z*b.z;
}
node cro(node a,node b){
return {a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};
}
struct face{
int a,b,c;
};
vector<face>q;
bool chk(node a,node b,node c,node d){
return dot(cro(b-a,c-a),d-a)>0;
}
int n;
set<pair<int,int> >jl;
void add(int x,int y){
if(jl.count({y,x}))jl.erase(jl.find({y,x}));
else jl.insert({x,y});
}
signed main(){
cin>>n;
for(int i=1;i<=n;i++)cin>>p[i].x>>p[i].y>>p[i].z,p[i].rnd();
q.push_back({1,2,3});
q.push_back({3,2,1});
for(int i=4;i<=n;i++){
vector<face>tmp;
jl.clear();
for(auto [a,b,c]:q){
if(chk(p[a],p[b],p[c],p[i])){
add(a,b);add(b,c);add(c,a);
}else tmp.push_back({a,b,c});
}
for(auto [a,b]:jl)tmp.push_back({a,b,i});
q=tmp;
}
double res=0;
for(auto [a,b,c]:q)res+=mc(cro(p[b]-p[a],p[c]-p[a]))/2;
printf("%.3lf",res+eps);
}
```