计算几何
1.向量
向量的定义
在二维及以上维度里,既有大小又有方向的量叫做向量
向量可以用有向线段表示
如图,从A到B的有向线段表示了一个向量,记作
注意: 向量可以用有向线段表示,但向量并不等同于有向线段,所有方向相同且长度相等的有向线段都可以表示同一个向量。
正因如此,我们可以将有向线段的起点平移至平面直角坐标系的原点,并用其终点的坐标来表示向量。
如图,
因为坐标系中的点坐标可以表示向量,所以我们可以将坐标系中的点
向量的模
向量的长度叫做向量的模,记作
如上图,
零向量
长度为
零向量的起点和终点重合,故零向量的方向不确定,也就是说零向量的方向是任意的,但零向量的大小是确定的。
相反向量
若两个向量大小相等,方向相反,则称这两个向量互为相反向量,记作
特别的,我们规定,零向量的相反向量为它本身。
表示成坐标形式就是
共线向量
方向相同或相反的非零向量叫做共线向量(或平行向量),记作
由于零向量的方向不确定,我们规定,零向量与任意向量平行,即
2.向量的运算
以下默认
向量的加法
向量既有大小,又有方向,其加减运算不能简单的将大小相加。
向量的加法满足平行四边形定则,即以参与运算的两个向量作为邻边作平行四边形,这两个向量的和即为该平行四边形的对角线,如下图:
因为平行四边形对边平行且相等,所以我们可以将向量的加法理解为将
表示成坐标形式就是
向量的加法满足交换律和结合律,即:
加法代码:
Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}
向量的减法
与实数运算中减去一个数等于加上这个数的相反数类似,减去一个向量等于加上这个向量的相反向量,即
由此我们可以将减法转化成加法。
从图中我们可以看出,
所以
根据这个几何意义,
减法代码:
Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}
向量的数乘
向量的数乘即向量与一个实数相乘,其实质是向量的缩放。
实数
当
当
当
表示成坐标形式就是
数乘满足以下运算律:
设
结合律:
向量对于数的分配律(第一分配律):
数对于向量的分配律(第二分配律):
数乘代码:
Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}
向量的数量积
又称点积、内积,记作
定义:
用坐标表示为:
坐标表示的推导:
如图,
根据余弦定理,
即
点积满足以下运算律:
设
交换律:
关于数乘的结合律:
分配律:
点积可用于判定两向量是否垂直:
点积代码:
double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
向量的向量积
又称叉积、外积,记作
叉积的模长:
叉积的方向垂直于参与运算的两向量所在的平面,且方向与参与运算的两向量的相对位置有关。具体来说:
如图,对于
若若
若
向量的叉积同样可以用坐标表示:
用坐标求出的值若为正数,则
若为负数,则
若为
利用这个性质我们可以判断两个向量的位置关系。
向量的叉积满足以下运算律:
设
关于数乘的结合律:
左分配律:
右分配律:
叉积代码:
double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}
3.求直线交点
如图,求直线
代码:
Vec getcrp(Line A,Line B) //crp=cross point 交点
{
double k=cross(B.a-A.a,B.b-B.a)/cross(A.b-A.a,B.b-B.a); //A.a和A.b表示直线A的两个端点,直线B同理
return A.a+(A.b-A.a)*k;
}
4.求多边形面积
如图,对于平面上的一个多边形,我们可以从某个点向其他所有顶点连线,将其分割成若干个三角形,三角形面积的和即为多边形面积。
那么三角形面积怎么求呢?用向量。
如图,从某一点向其他所有点连向量,因为向量的叉积可以求两向量作为邻边围成的平行四边形的面积,而三角形的面积是平行四边形的面积的一半,这样我们就可以用叉积求出多边形的面积。
那么可能有人会问了:如果多边形不是凸多边形怎么办呢?
确实,如果多边形非凸,我们会算上一些多余部分。
但是我们可以按顺序依次算叉积,即计算
因为用坐标计算的叉积,若第二个向量在第一个向量的顺时针方向,算出来的结果是负数,那么加上这个负数就相当于减去的多余部分。
如图,当我们加上
因此,只要我们依次计算叉积,就可以保证求出来的面积就是多边形的面积,也就可以保证这样求面积的方法是普遍适用的。(当然,很少有题目会让你求非凸多边形的面积)
5.判断点与线、线与线的位置关系
直线的左半平面与右半平面
在计算几何中,直线由有顺序的两个点确定,如下图为直线
那么我们定义,直线
相似的,直线
特别的,直线
点与线的位置关系(向量与向量的位置关系)
判断点与线的位置关系,其实就是判断点在直线的左半平面,直线上还是右半平面。
如图,判断点
前面在讲叉积时说过,用坐标算出来的叉积结果可以判断参与运算的两个向量之间的位置关系,那么我们可以用同样的办法判断点与直线的位置关系。
计算
若为正数,则
若为负数,则
若为
线与线的位置关系
将直线上的两个点表示成向量,用向量的运算判断。
平行
相交
两条直线不是平行就是相交
垂直
6.二维凸包
定义:能包含平面上给定的所有点的最小凸多边形。
性质:凸包上的边按逆时针顺序排列,第
求法:Graham扫描法
首先选取
由于凸包第
怎么判断新加入的点是否满足性质呢?
我们可以把加入栈的点之间依次连接向量,将栈顶的点与新点连接向量,然后判断两个向量之间的位置关系。
如图,
那么我们判断
若
否则边
最后在将栈顶和栈底的点连边,就可以构成一个完整的凸包。
代码(Luogu P2742):
#include <bits/stdc++.h>
#define maxn 1000005
using namespace std;
struct Vec{double x,y;}v[maxn];
Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}
Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}
Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}
double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}
double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
double angle(Vec a){return atan2(a.y,a.x);}
double dist(Vec a){return sqrt(a.x*a.x+a.y*a.y);}
double sqrdist(Vec a){return a.x*a.x+a.y*a.y;}
Vec O;//极点
bool cmp(Vec a,Vec b)
{
double k1=angle(a-O);
double k2=angle(b-O);
if (k1!=k2) return k1<k2;
if (a.x!=b.x) return a.x<b.x;
return a.y<b.y;
}
int n;
int main()
{
scanf("%d",&n);
for (int i=1;i<=n;i++)
{
scanf("%lf%lf",&v[i].x,&v[i].y);
if (i==1) O=v[i];
if (v[i].y<O.y||(v[i].y==O.y&&v[i].x<O.x)) O=v[i];
}
sort(v+1,v+1+n,cmp);
Vec s[maxn];
int top=0;
for (int i=1;i<=n;i++)
{
while(top>1&&cross(s[top]-s[top-1],v[i]-s[top])<=0) top--;
//top>1是因为若栈中只剩一个点是无法判断的,此时新点必须入栈才能进一步判断之后的点
s[++top]=v[i];
}
double ans=0.0;
for (int i=1;i<=top;i++)
{
int tar=i%top+1;//逆时针方向的下一个点
ans+=dist(s[tar]-s[i]);
}
printf("%.2lf",ans);
return 0;
}
7.判断点是否在凸包内
- 做法一:
O(n) 枚举
如图,判断点
将凸包上的边按逆时针顺序连接向量,并将凸包所有节点与
若
否则,
- 做法二:二分
如图,判断点
首先将点
然后,令
接下来,我们判断点
最后,
此时点
Code:
//d为需判断的点,E为凸包顶点坐标
//1为在凸包内,0反之
if (cross(E[N-1]-E[0],d-E[0])>0){puts("0");continue;}
if (cross(E[1]-E[0],d-E[0])<0){puts("0");continue;}
int l=1,r=N-1;
while(l+1<r)
{
int mid=(l+r)>>1;
if (cross(E[mid]-E[0],d-E[0])>0) l=mid;
else r=mid;
}
if (cross(E[r]-E[l],d-E[l])>=0) puts("1");
else puts("0");
8.闵可夫斯基和
定义:两个多边形
其中
此处我们只讨论凸包的闵可夫斯基和。
其实就是将其中一个多边形的某个顶点放在另一个多边形上绕一圈,如图:
红色凸包与蓝色凸包的闵可夫斯基和即为图中整个多边形。
由图中可看出,闵可夫斯基和对应的边就是将两个凸包的边按极角排序后一次收尾相接,根据这个性质我们可以求出闵可夫斯基和的形状。
那么闵可夫斯基和的位置怎么求呢?
因为闵可夫斯基和的定义是向量相加,所以闵可夫斯基和
最后记得将相邻的且极角相等的边合并,不然会有奇怪的错误
例题:P4557 [JSOI2018] 战争
给出两个凸包
我们考虑将这个式子移一下项,可得:
即:
我们发现,原条件等价于
回看一下闵可夫斯基和的定义,我们发现,就是将
Code:
#include <bits/stdc++.h>
#define maxn 100005
using namespace std;
const double pi=acos(-1);//圆周率
struct Vec{double x,y;};//向量结构体
bool operator == (Vec a,Vec b){return a.x==b.x&&a.y==b.y;}
bool operator != (Vec a,Vec b){return !(a==b);}
Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}
Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}
Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}
double angle(Vec a){return atan2(a.y,a.x);}//极角
double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}//叉积
double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
double dist(Vec a){return sqrt(a.x*a.x+a.y*a.y);}
double sqrdist(Vec a){return a.x*a.x+a.y*a.y;}
Vec O;
bool cmp(Vec a,Vec b)//求凸包时的排序
{
double k1=angle(a-O);
double k2=angle(b-O);
if (k1!=k2) return k1<k2;
if (a.x!=b.x) return a.x<b.x;
return a.y<b.y;
}
void hull(Vec *v,int n,Vec *e,int &num)//求凸包
{
sort(v+1,v+1+n,cmp);
static Vec st[maxn];
int top=0;
for (int i=1;i<=n;i++)
{
while(top>1&&cross(st[top]-st[top-1],v[i]-st[top])<=0) top--;
st[++top]=v[i];
}
for (int i=1;i<=top;i++)
{
int tar=(i%top)+1;
e[++num]=st[tar]-st[i];//记录凸包的边
}
}
int n,m,q;
Vec v1[maxn],v2[maxn];
Vec e[maxn<<1];//存储两个凸包的边,记得开两倍
Vec E[maxn<<1];
int num,N;
bool cmp2(Vec a,Vec b)
{
double k1=angle(a);
double k2=angle(b);
//因为atan2函数求出的极角是-π到π之间的,而排序时要按照极角在0到2π之间的值排序,所以将小于0的极角加上2π
if (k1<0) k1+=pi*2;
if (k2<0) k2+=pi*2;
return k1<k2;
}
int main()
{
scanf("%d%d%d",&n,&m,&q);
for (int i=1;i<=n;i++)
{
scanf("%lf%lf",&v1[i].x,&v1[i].y);
if (i==1) O=v1[i];
if (v1[i].y<O.y||(v1[i].y==O.y&&v1[i].x<O.x)) O=v1[i];
}
hull(v1,n,e,num);
for (int i=1;i<=m;i++)
{
scanf("%lf%lf",&v2[i].x,&v2[i].y);
v2[i].x*=-1;v2[i].y*=-1;
if (i==1) O=v2[i];
if (v2[i].y<O.y||(v2[i].y==O.y&&v2[i].x<O.x)) O=v2[i];
}
hull(v2,m,e,num);
sort(e+1,e+1+num,cmp2);//按极角排序
for (int i=1;i<=num;i++)
{
if (i==1||angle(e[i])!=angle(e[i-1])) E[++N]=e[i];
else E[N]=E[N]+e[i];//合并极角相等的边
}
E[0]=v1[1]+v2[1];//闵可夫斯基和的y坐标最小的点
for (int i=1;i<=N;i++) E[i]=E[i-1]+E[i];//求出各顶点的坐标
for (int i=1;i<=q;i++)
{
Vec d;
scanf("%lf%lf",&d.x,&d.y);
if (cross(E[N-1]-E[0],d-E[0])>0){puts("0");continue;}
if (cross(E[1]-E[0],d-E[0])<0){puts("0");continue;}
int l=1,r=N-1;
while(l+1<r)
{
int mid=(l+r)>>1;
if (cross(E[mid]-E[0],d-E[0])>0) l=mid;
else r=mid;
}
if (cross(E[r]-E[l],d-E[l])>=0) puts("1");
else puts("0");
}
return 0;
}
练习:CF87E Mogohu-Rea Idol
题解link
9.旋转卡壳
一种重要算法,可用于求凸包直径,也就是求凸包上任意两点间距离的最大值。
做法:先求出凸包,然后逆时针枚举凸包上的边,对于每条边找到凸包上距离该边最远的点,设为
由于逆时针旋转后,新的
如何判断
如图,
要比较
因为
面积怎么求?用叉积算出平行四边形面积在除以二即可(其实不除也行,因为比较面积的两倍跟比较面积是等价的)
代码(Luogu P1452):
#include <bits/stdc++.h>
#define maxn 1000005
using namespace std;
struct Vec{double x,y;}v[maxn];
Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}
Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}
Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}
double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}
double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
double angle(Vec a){return atan2(a.y,a.x);}
double dist(Vec a){return sqrt(a.x*a.x+a.y*a.y);}
double sqrdist(Vec a){return a.x*a.x+a.y*a.y;}
Vec O;
bool cmp(Vec a,Vec b)
{
double k1=angle(a-O);
double k2=angle(b-O);
if (k1!=k2) return k1<k2;
if (a.x!=b.x) return a.x<b.x;
return a.y<b.y;
}
int n;
int main()
{
scanf("%d",&n);
for (int i=1;i<=n;i++)
{
scanf("%lf%lf",&v[i].x,&v[i].y);
if (i==1) O=v[i];
if (v[i].y<O.y||(v[i].y==O.y&&v[i].x<O.x)) O=v[i];
}
sort(v+1,v+1+n,cmp);
Vec s[maxn];
int top=0;
for (int i=1;i<=n;i++)
{
while(top>1&&cross(s[top]-s[top-1],v[i]-s[top])<=0) top--;
s[++top]=v[i];
}
double ans=0.0;
int k=2;
for (int i=1;i<=top;i++)
{
int tar=i%top+1;
while(cross(s[tar]-s[i],s[k%top+1]-s[i])>cross(s[tar]-s[i],s[k]-s[i])) k=k%top+1;
ans=max(ans,sqrdist(s[k]-s[i]));
ans=max(ans,sqrdist(s[k]-s[tar]));
}
printf("%.0lf",ans);
return 0;
}
10.半平面交
求若干个半平面的交集,这里的有效半平面定义为直线的左半平面。
求法:
首先将所有的直线按极角排序,极角相同的只保留最左侧的那条直线。
显然最后的交是一个凸多边形,我们使用deque(双端队列)来维护。
考虑新加入一个半平面(一条直线)
- 情况一:
(注:图中直线
如图,若deque中最后两条直线的交点在新加入的直线的右侧,显然deque尾端的直线对于半平面的交已经没有了影响,可以弹出。
- 情况二:
如图,若deque中最后两条直线的交点在新加入的直线的左侧,那么此时deque尾端的直线显然都还对半平面的交有影响,不必弹出,此时直接将新的直线入队。
这样看起来,我们用单调栈貌似也能实现目的,为什么我们要用deque呢?
- 情况三
(注:图中直线
如图,deque中的前两条直线的交点也有可能在新加入的直线的右侧,此时显然deque头端的直线对半平面的交没有了影响,所以我们要弹出队头直线。
- 情况四
(注:此处图示中直线的含义与情况三图示相同)
如图,若deque中的前两条直线的交点在新加入的直线的左侧,同情况二,直接将新的直线入队。
- 最后
由于deque中前两条直线的交点可能在队尾直线的右侧,此时队头直线是多余的。
并且最后两条直线的交点也有可能在队头直线的右侧,此时队尾直线是多余的。
所以,我们还要将队头和队尾多余的直线弹出。
例题:Luogu P4196
这题说是让我们求凸多边形的交,但实际上,每个凸多边形可以看成若干个半平面的交,多边形的每条边所在的直线代表着一个半平面,这样我们就将原问题转化为了半平面交。
代码(注意,若要使用本代码解决其他半平面交问题,请确保构造的直线的有效半平面为左半平面):
#include <bits/stdc++.h>
#define maxn 1005
using namespace std;
struct Vec{double x,y;}v[maxn];
struct Line{Vec a,b;double k;}L[maxn],l[maxn];
Vec operator + (Vec a,Vec b){return (Vec){a.x+b.x,a.y+b.y};}
Vec operator - (Vec a,Vec b){return (Vec){a.x-b.x,a.y-b.y};}
Vec operator * (Vec a,double b){return (Vec){a.x*b,a.y*b};}
double cross(Vec a,Vec b){return a.x*b.y-a.y*b.x;}
double dot(Vec a,Vec b){return a.x*b.x+a.y*b.y;}
double angle(Vec a){return atan2(a.y,a.x);}
double dist(Vec a){return sqrt(a.x*a.x+a.y*a.y);}
double sqrdist(Vec a){return a.x*a.x+a.y*a.y;}
Vec getcrp(Line A,Line B)
{
double k=cross(B.a-A.a,B.b-B.a)/cross(A.b-A.a,B.b-B.a);
return A.a+(A.b-A.a)*k;
}
bool cmp(Line A,Line B)
{
if (A.k!=B.k) return A.k<B.k;
return cross(B.a-A.a,B.b-A.a)>0;//若斜率相等,就将在最左侧的直线排到最前面
}
bool check(Line l,Vec p){return cross(l.b-l.a,p-l.a)<0;}//判断点是否在直线的右侧
int n,N,num;
int main()
{
scanf("%d",&n);
for (int i=1;i<=n;i++)
{
int m;
scanf("%d",&m);
for (int j=1;j<=m;j++) scanf("%lf%lf",&v[j].x,&v[j].y);
for (int j=1;j<=m;j++)
{
int tar=j%m+1;
N++;
L[N].a=v[j];L[N].b=v[tar];L[N].k=angle(v[tar]-v[j]);
}
}
sort(L+1,L+1+N,cmp);
for (int i=1;i<=N;i++)
if (i==1||L[i].k!=L[i-1].k) l[++num]=L[i];
Line q[maxn*2];
int h=1,t=0;
Vec crp[maxn];
for (int i=1;i<=num;i++)
{
while(h<t&&check(l[i],crp[t])) t--;//情况一
while(h<t&&check(l[i],crp[h+1])) h++;//情况三
q[++t]=l[i];
if (h<t) crp[t]=getcrp(q[t-1],q[t]);
}
while(h<t&&check(q[t],crp[h+1])) h++;//弹出多余直线
while(h<t&&check(q[h],crp[t])) t--;//弹出多余直线
crp[h]=getcrp(q[h],q[t]);
double ans=0.0;
for (int i=h+1;i<t;i++) ans+=cross(crp[i]-crp[h],crp[i+1]-crp[h]);//求面积
printf("%.3lf",ans/2);
return 0;
}