计算几何

计算几何

计算几何

点和向量

首先二维平面的点用 (x,y) 表示,向量是一个具有方向和大小的线段,向量 (x,y) 可以理解为从原点 O(0,0)(x,y) 连接出去的一条线段,同时具有方向。

当然,向量可以平移,平移后的向量用 (A,B) 表示,A,B 均为点,此时的向量等同于之前的 (BxAx,ByAy)(即把 A 平移到原点的位置)

两点之间距离公式:(AxBx)2+(AyBy)2,这就是一个向量的大小。

向量的运算

  • 向量 (x,y) 加点 A,这个可以理解为将 A 点顺着向量的方向移动,结果就是一个新的点 (x+Ax,y+Ay)

  • 向量 (x,y) 加/减向量 (x,y),得到 (x+x,y+y),这个没有太大的几何意义。

  • 向量 (x,y) 乘/除实数 k 得到 (xk,yk),得到等比例放大/缩小的向量。

点积

|A| 表示向量 A 的大小,那么定义向量 A 和向量 B 的点积是 AB=|A||B|cosθ,其中 θA,B 之间的夹角。

点积的几何意义就是 AB 上的投影长度乘以 B 的模长。

点积并不是很常用,计算它也不需要得到 θ,只需要 AxBx+AyBy 就好,公式推导需要用到三角函数,此处不再赘述。

应用:若 AB>0,则 θ 为锐角;若 AB=0,则 θ 为直角,否则 θ 为钝角。

叉积

|A| 表示向量 A 的大小,那么定义向量 A 和向量 B 的叉积是 A×B=|A||B|sinθ,其中 θA,B 之间的夹角。

叉积的几何意义就是 AB 两个向量形成的平行四边形的有向面积,即如果 AB 的上方,它就是负数;否则是正数。(逆时针围成的图形面积)

如果 AB 在同一直线(方向可能不同),那么 A×B=0

计算叉积也不需要得到 θ,只需要 AxByAyBx 就好,公式推导需要用到三角函数,此处不再赘述。

应用

如果 A×B>0,则 BA 的逆时针方向,如果 A×B<0,则 BA 的顺时针方向;否则 A,B 共线。

如果计算三个点构成的有向平行四边形的面积,那就让其中某个点作为原点,然后计算新的向量 A×B 即可。

如果是计算三角形,答案取绝对值然后除以 2 就可以了。

因为向量可以平移,我们还可以用叉积来判断两个向量是否平行或者重合。

旋转

向量 A 逆时针旋转 θ,则得到新的向量 A=(xcosθysinθ,xsinθ+ycosθ)

如果是旋转 π2 的倍数,那么交换 x,y 再视情况取反就可以了,这个分类讨论不难。

直线、射线和线段

直线表达式常用的有两种写法 ax+by+c=0y=kx+b,两种写法都比较好用,但是我更喜欢后面的一种写法。

线段可以用两个点直接表示,同时也可以写出包含此线段的直线。

下面设直线 l(A,B) 表示包含两个点 A,B 的直线 l,则我们可以得到:

点和直线的位置关系

直线上任取两个点 A,B,与另外的点 P 构成两个向量 (B,P),(A,B),如果这两个向量的叉积 >0,代表 P 在直线右侧;<0 P 在直线左侧;否则 P 在直线上。

同理可以定义点和线段的位置关系,这个用上面的分类讨论一下就可以了。

点到直线距离

已知点 P 和直线 l(A,B),可以求出 AB 线段的长度,和 P,A,B 用叉积求出平行四边形面积,除以 AB 长度就是高,即点 Pl 的距离。

点在直线上的投影

已知直线上两个点 p1,p2,以及直线外一个点 ppl 上的投影 p0

k=|p0p1||p2p1|,那么求得 k 就知道了 p0 的位置,因为 k 是这两条线段的比值。

根据点积的概念有 A(p1,p)B(p1,p2)=A(p1,p0)B(p1,p2),所以 |p0p1|=(pp1)(p2p1)|p2p1|,带入 k 得到 k=(pp1)(p2p1)|p2p1||p2p1|,所以 p0=p1+k(p2p1)=p1+(pp1)(p2p1)|p2p1||p2p1|(p2p1)

同理可以得到点对于一条直线的对称点。

点到线段的距离

下面三个取最小值:

  • 点到线段某个端点的距离。

  • 点到这个直线的距离。(点的投影在线段上)

直线的位置关系

两条直线的位置关系直接使用叉积判断即可。

两条线段的位置关系同理可以应用叉积的判断:

  • 如果每条线段的两个端点都在另外一条线段的直线的两侧,那么这两条线段相交。

于是用最简单的叉积相乘 <0 判断即可。

直线交点

可以通过两条直线的 y=kx+b 联立得到交点,但是这样的话精度和特判有点麻烦,我们可以考虑更加简单的做法。

如果第一条直线上有两个点 A,B,第二条直线上有两个点 C,DPA,B,C,D 构成的四边形中间,也就是对角线交点,那么就有相似 |DP||CP|=SABDSABC=(A,D)×(A,B))(A,B)×(A,C)

注意:AB 的下方,CD 的左方。(有的时候不需要调整顺序)

那么就有 |DP||CP|=DxPxPxCx=DyPyPyCy,联立两个方程,得到:

Px=SABDCx+SABCDxSABD+SABC  Py=SABDCy+SABCDySABD+SABC

于是用叉积求三角形面积即可。

线段的交点也就比较简单了,判断是否有交即可,然后求直线交点。

多边形

点和多边形的关系

一种方法就是暴力枚举每条边,然后对于每条边的两个端点按照顺序与这个点 P 构成的向量作叉乘,这样的话我们会得到 n 个结果,这 n 个结果要么都是负数,要么都是正数,如果出现其它情况那就说明这个点在多边形外面。

上面这种方法只能用于凸多边形,更加常用的方法是引出来一条射线,然后判断这条射线与多少条边相交,如果是奇数,那么就在多边形内部。

但是需要注意的是我们可能需要进行多次随机取值,因为如果划到了端点就不好判断了,时间复杂度依然是 O(n)

多边形的面积

直接对于每个边的两个端点与原点构成的向量按照顺序作叉乘就可以了。

最后需要除以 2 并且取绝对值,因为每次求出的相当于是一个三角形的面积,这些三角形因为叉乘有正负,多加的会被容斥掉,所以这么做对凸多边形和凹多边形都可以。

凸包

例题:P2742 [USACO5.1] 圈奶牛Fencing the Cows /【模板】二维凸包

大意:给定平面上若干个点,求出包含这些点的最小凸多边形的长度,这个凸多边形被称为凸包。

解法如下:

首先我们将所有点按照 x 坐标排序,x 坐标相同的按照 y 坐标排序,这两个维度都是从小到大排序,分上下凸包考虑,此处以下凸包为例,设 p1,p2,,pk 是当前凸包的点集,容易发现排序之后最左边的点和最右边的点都在凸包上。

接下来考虑加入一个点 x,那么如果 (pk1,pk)×(pk1,x)<0,说明 (pk1,x)(pk1,pk) 下方,应该弹出 pk

所以每个元素最多被弹出一次,最后处理完上下凸包之后拼起来就可以了。(上凸包需要倒序枚举进栈,并且最后一个点不需要重复进栈)

时间复杂度 O(nlogn),代码如下:

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
#include<bits/stdc++.h>
#define eps 1e-8
#define ll long long
#define db double
#define N 100005
using namespace std;
struct point{db x,y;}p[N],sta[N];
inline point operator-(point a,point b){return (point){a.x-b.x,a.y-b.y};}
inline ll sgn(double x){
if(fabs(x)<eps) return -1;
if(x>eps) return 1;
return 0;
}
inline bool cmp(point a,point b){
if(sgn(a.x-b.x)==-1) return sgn(b.y-a.y);
else return sgn(b.x-a.x);
}
inline ll cross(point a,point b){return a.x*b.y-a.y*b.x;}
inline db lenth(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
ll n,i,top,pos;
db ans;
int main(){
ios::sync_with_stdio(false);
cin>>n,cout<<fixed<<setprecision(2);
for(i=1;i<=n;i++) cin>>p[i].x>>p[i].y;
sort(p+1,p+n+1,cmp);
for(i=1;i<=n;i++){
while(top>=2&&cross(sta[top]-sta[top-1],p[i]-sta[top-1])<=0) top--;
sta[++top]=p[i];
}
pos=top;
for(i=n-1;i>=1;i--){
while(top>=pos+1&&cross(sta[top]-sta[top-1],p[i]-sta[top-1])<=0) top--;
sta[++top]=p[i];
}
for(i=1;i<top;i++) ans+=lenth(sta[i],sta[i+1]);
cout<<ans<<endl;
return 0;
}

这个算法被称为 Andrew 算法。

点是否在凸包内部

凸包内最左边的那个点作为原点,剩下的点按照逆时针通过叉积排序就可以了。

查找的时候二分出来它应该插入到哪个位置,得到前后两个点,然后根据叉积判断位置关系就可以了。

还要特判边界情况,即在凸包最后一个点的外侧或者第一个点的外侧或者和原点重合了。

时间复杂度 O(logn),甚至不需要扫描每一条边。

闵可夫斯基和

给定两个凸包 A,B,构建一个凸包 CA,B 闵可夫斯基和,即满足 aA,bB|a+bC,其中 a,b 均为向量,aA 为在凸包 A 中的向量,a+b 为向量加法。

首先有一个性质就是两个凸包的闵可夫斯基和还是凸包,并且 C 的边的条数恰好是 A,B 两个凸包的边的数量之和,更明确一点,C 的边集构成的向量集合就是 A,B 两个凸包构成的向量集合的并集。

因此我们可以直接对于所有 A,B 中的边按照逆时针排序,这个可以用叉积解决,然后固定一下 C 中的一个点,这里一般选择 A 的第一个点和 B 的第一个点加起来的结果。(可以看成是一个向量)

特别注意:C 中的第一个点连接的边也得是 A 中的第一条边和 B 中的第一条边作为向量靠下面(逆时针排序后的第一个)的那一个,因此我们不能直接 sort 排序,而是使用归并排序来合并 A,B 两个向量集合。

最后我们就在 O(n) 的时间复杂度内求出了这个凸包 C

最后还要再对闵可夫斯基和求一遍凸包才能通过,因为有可能存在三点共线的情况。

例题:P4557 [JSOI2018] 战争

这道题让你求出每次 A 沿着向量 w 平移之后与 B 有没有公共点,AB 是我们根据给定的点集算出来的凸包。

那么若有公共点则存在 aA,bB 满足 b+w=a,移项得 w=ab,那么就是求 w 在不在 AB 的闵可夫斯基和内部,直接用 O(logn) 求点是否在闵可夫斯基和内部即可。

时间复杂度 O(nlogn)

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
#include<bits/stdc++.h>
#define eps 1e-8
#define ll long long
#define db double
#define N 200005
using namespace std;
struct point{ll x,y;}p1[N],p2[N],sta[N],g[N],p3[N],p4[N],st;
inline point operator+(point a,point b){return (point){a.x+b.x,a.y+b.y};}
inline point operator-(point a,point b){return (point){a.x-b.x,a.y-b.y};}
inline ll sgn(double x){
if(fabs(x)<eps) return -1;
if(x>eps) return 1;
return 0;
}
inline bool cmp(point a,point b){
if(sgn(a.x-b.x)==-1) return sgn(b.y-a.y);
else return sgn(b.x-a.x);
}
inline ll cross(point a,point b){return a.x*b.y-a.y*b.x;}
inline ll lenth(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
inline bool cmps(point a,point b){return cross(a,b)>0;}
ll n,m,q,i,j,top,pos,x,y,tot;
inline void andrew(){
sort(p1+1,p1+n+1,cmp),sort(p2+1,p2+m+1,cmp);
top=0;
for(i=1;i<=n;i++){
while(top>=2&&cross(sta[top]-sta[top-1],p1[i]-sta[top-1])<=0) top--;
sta[++top]=p1[i];
}
pos=top;
for(i=n-1;i>=1;i--){
while(top>=pos+1&&cross(sta[top]-sta[top-1],p1[i]-sta[top-1])<=0) top--;
sta[++top]=p1[i];
}
n=top-1;
for(i=1;i<top;i++) p1[i]=sta[i];
top=0;
for(i=1;i<=m;i++){
while(top>=2&&cross(sta[top]-sta[top-1],p2[i]-sta[top-1])<=0) top--;
sta[++top]=p2[i];
}
pos=top;
for(i=m-1;i>=1;i--){
while(top>=pos+1&&cross(sta[top]-sta[top-1],p2[i]-sta[top-1])<=0) top--;
sta[++top]=p2[i];
}
m=top-1;
for(i=1;i<top;i++) p2[i]=sta[i];
}
inline void Minkovski(){
g[++tot] = (point){p1[1].x+p2[1].x,p1[1].y+p2[1].y};
for(i=1;i<=n;i++) p3[i]=p1[i%n+1]-p1[i];
for(i=1;i<=m;i++) p4[i]=p2[i%m+1]-p2[i];
for(i=1,j=1;i<=n&&j<=m;) tot++,g[tot]=g[tot-1]+(cross(p3[i],p4[j])>0?p3[i++]:p4[j++]);
for(;i<=n;i++) tot++,g[tot]=g[tot-1]+p3[i];
for(;j<=m;j++) tot++,g[tot]=g[tot-1]+p4[j];
tot--,sort(g+1,g+tot+1,cmp);
top=0;
for(i=1;i<=tot;i++){
while(top>=2&&cross(sta[top]-sta[top-1],g[i]-sta[top-1])<=0) top--;
sta[++top]=g[i];
}
pos=top;
for(i=tot-1;i>=1;i--){
while(top>=pos+1&&cross(sta[top]-sta[top-1],g[i]-sta[top-1])<=0) top--;
sta[++top]=g[i];
}
tot=top-1;
for(i=1;i<=tot;i++) g[i]=sta[i];
st=g[1];
for(i=1;i<=tot;i++) g[i]=g[i]-st;
sort(g+2,g+tot+1,cmps);
}
inline ll query(point x){
if(tot==1){
if(x.x==g[1].x&&x.y==g[1].y) return 1;
return 0;
}
if(cross(x,g[2])>0||cross(g[tot],x)>0) return 0;
ll res = lower_bound(g+2,g+tot+1,x,cmps)-g-1;
return cross(g[res+1]-g[res],x-g[res])>=0||(cross(g[res+1],x)==0&&((g[res+1].x>0&&x.x>0&&x.x<=g[res+1].x)||(g[res+1].x<0&&x.x<0&&x.x>=g[res+1].x)));
}
int main(){
ios::sync_with_stdio(false);
cin>>n>>m>>q;
for(i=1;i<=n;i++) cin>>p1[i].x>>p1[i].y;
for(i=1;i<=m;i++) cin>>p2[i].x>>p2[i].y,p2[i].x=-p2[i].x,p2[i].y=-p2[i].y;
andrew();
Minkovski();
while(q--){
cin>>x>>y;
cout<<query((point){x,y}-st)<<endl;
}
return 0;
}

注意上面那份代码最后还要求一次凸包,原因在于可能用重复计算的点。

例题 2:here

题目大意:给定长度为 n 的序列 A,定义长度为 k 的子序列的权值为奇数位上的和减去偶数位上的和,求出所有长度为 1kn 的子序列的最大权值。

其中:1n1050Ai2×104

首先如果 k 是偶数情况的话,我们可以令 BiAiAi+1,然后选出连续的 k 段区间使得区间和最大,这个可以用费用流解决,所以设 f(i) 表示 k=i 的答案,那么对于 imod2=0 的情况 f 明显是凸函数,因为费用流每次增广的流量是递减的,奇数的情况同理。

于是就可以设 (x,y) 坐标表示当 i=x 的时候存在 f(i)=y,于是我们发现所有 imod2=1 对于所有区间都有凸性,所以我们合并的时候用闵可夫斯基和就好。

因为凸包和凸包的闵可夫斯基和也是凸包,并且根据它的定义我们知道这么做一定是正确的。

具体来说,线段树每个节点维护一个平面直角坐标系上的凸包,凸包每个点的定义就是上面说的那样,然后在 pushup 的过程中根据题目的要求(一正一负)记录一下正负的情况合并就可以了。

记得特判掉边界情况。

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#include<bits/stdc++.h>
#pragma GCC optimize("Ofast")
#define ll long long
#define N 100005
using namespace std;
struct point{ll x,y;}g[N],p3[N],p4[N];
vector<point> tr[N<<2][2][2],st;
inline point operator+(point a,point b){return (point){a.x+b.x,a.y+b.y};}
inline point operator-(point a,point b){return (point){a.x-b.x,a.y-b.y};}
inline ll lenth(point a){return a.x*a.x+a.y*a.y;}
inline bool cmp(point a,point b){
if(a.x==b.x) return a.y<b.y;
else return a.x<b.x;
}
inline ll cross(point a,point b){return a.x*b.y-a.y*b.x;}
ll n,i,j,a[N],tot,ans[N],temp[N];
inline vector<point> Minkovski(vector<point> as,vector<point> bs){
if(as.size()==0) return bs;
if(bs.size()==0) return as;
ll i,j;
tot=0,st.clear();
g[++tot] = (point){as[0].x+bs[0].x,as[0].y+bs[0].y};
for(i=1;i<as.size();i++) p3[i]=as[i]-as[i-1];
for(i=1;i<bs.size();i++) p4[i]=bs[i]-bs[i-1];
for(i=1,j=1;i<as.size()&&j<bs.size();){
tot++;
if(cross(p3[i],p4[j])<0) g[tot]=g[tot-1]+p3[i++];
else g[tot]=g[tot-1]+p4[j++];
}
for(;i<as.size();i++) tot++,g[tot]=g[tot-1]+p3[i];
for(;j<bs.size();j++) tot++,g[tot]=g[tot-1]+p4[j];
for(i=1;i<=tot;i++) st.push_back(g[i]);
return st;
}
inline vector<point> max(vector<point> a,vector<point> b,vector<point> c,vector<point> d,ll len){
ll beg = (len%2==0?2:1);
st.clear();
for(ll i=beg;i<=len;i+=2) temp[i]=-0x3f3f3f3f3f3f3f3f;
temp[0]=0;
for(ll i=0;i<a.size();i++) temp[a[i].x]=max(temp[a[i].x],a[i].y);
for(ll i=0;i<b.size();i++) temp[b[i].x]=max(temp[b[i].x],b[i].y);
for(ll i=0;i<c.size();i++) temp[c[i].x]=max(temp[c[i].x],c[i].y);
for(ll i=0;i<d.size();i++) temp[d[i].x]=max(temp[d[i].x],d[i].y);
for(ll i=beg;i<=len;i+=2) st.push_back((point){i,temp[i]});
return st;
}
inline void build(ll s,ll t,ll p){
if(s==t){
tr[p][0][0].push_back((point){1,a[s]});
tr[p][1][1].push_back((point){1,-a[s]});
return ;
}
build(s,(s+t)/2,2*p),build((s+t)/2+1,t,2*p+1);
ll len1,len2;
if((t-s+1)%2==0) len1=t-s,len2=t-s+1;
else len1=t-s+1,len2=t-s;
tr[p][0][0] = max(Minkovski(tr[2*p][0][1],tr[2*p+1][0][0]),Minkovski(tr[2*p][0][0],tr[2*p+1][1][0]),tr[2*p][0][0],tr[2*p+1][0][0],len1);
tr[p][0][1] = max(Minkovski(tr[2*p][0][1],tr[2*p+1][0][1]),Minkovski(tr[2*p][0][0],tr[2*p+1][1][1]),tr[2*p][0][1],tr[2*p+1][0][1],len2);
tr[p][1][0] = max(Minkovski(tr[2*p][1][1],tr[2*p+1][0][0]),Minkovski(tr[2*p][1][0],tr[2*p+1][1][0]),tr[2*p][1][0],tr[2*p+1][1][0],len2);
tr[p][1][1] = max(Minkovski(tr[2*p][1][1],tr[2*p+1][0][1]),Minkovski(tr[2*p][1][0],tr[2*p+1][1][1]),tr[2*p][1][1],tr[2*p+1][1][1],len1);
}
int main(){
ios::sync_with_stdio(false);
cin>>n;
for(i=1;i<=n;i++) cin>>a[i];
build(1,n,1);
memset(ans,-0x3f,sizeof(ans));
for(i=0;i<tr[1][0][0].size();i++) ans[tr[1][0][0][i].x]=max(ans[tr[1][0][0][i].x],tr[1][0][0][i].y);
for(i=0;i<tr[1][0][1].size();i++) ans[tr[1][0][1][i].x]=max(ans[tr[1][0][1][i].x],tr[1][0][1][i].y);
for(i=1;i<=n;i++) cout<<ans[i]<<" ";
return 0;
}

极坐标系

一个点的位置在极坐标系上可以用下面的方式表示:

  • 原点 (0,φ)

A 为平面上一点。

极点 OA 之间的距离 |OA| 称为 极径,记为 ρ
以极轴为始边,OA 为终边的角 xOA 称为 极角,记为 φ

A=(ρ,φ)

我们可以很简单的将所有点按照极角排序,我们可以直接使用 atan 还原极角,但是这样的话分类讨论的情况较多,我们可以直接使用系统自带的函数 atan2 解决。

即设 A=(x,y)(平面直角坐标系),A 的极角就是 atan2(y,x)A 的极径就是 x2+y2

有时候极角排序比用叉积排序更加方便。