学习笔记——计算几何
前言
为了测试你板子的正确性，可以交到 Aizu 上去。我的代码都是 Shared 的。
记得开 double！记得开 double！记得开 double！记得开 double！记得开 double！记得开 double！记得开 double！记得开 double！记得开 double！记得开 double！记得开 double！
在此之前
为了一些精度的方便判断，我们定义一些前置的东西：
const double eps=1e-8;
const double PI=acos(-1);
int sgn(double x){
if(fabs(x)<eps) return 0;
if(x<0) return -1;else return 1;
}//判断 x 的正负
向量及其基本运算
先上代码：
struct Node{
double x,y;Node(){}
Node(double _x,double _y){x=_x;y=_y;}
void input(){cin>>x>>y;}
Node friend operator+(Node a,Node b){return Node(a.x+b.x,a.y+b.y);}
Node friend operator-(Node a,Node b){return Node(a.x-b.x,a.y-b.y);}
Node friend operator*(Node a,double b){return Node(a.x*b,a.y*b);}
double friend operator*(Node a,Node b){return a.x*b.x+a.y*b.y;}
double friend operator^(Node a,Node b){return a.x*b.y-b.x*a.y;}
void rot(double a){
double x1=x*cos(a)-y*sin(a);
double y1=y*cos(a)+x*sin(a);
x=x1;y=y1;
}
};
定义一个向量
根据向量的坐标表示，可以用 \((x,y)\) 来表示一个向量。
struct Node{
double x,y;Node(){}
Node(double _x,double _y){x=_x;y=_y;}
void input(){scanf("%lf%lf",&x,&y);}
};
向量加减
由于我们用坐标表示一个向量，那么两个向量相加减实际上就是对应坐标的加减。
Node friend operator+(Node a,Node b){
return Node(a.x+b.x,a.y+b.y);
}
Node friend operator-(Node a,Node b){
return Node(a.x-b.x,a.y-b.y);
}
向量数乘
大家都会。
Node friend operator*(Node a,double b){
return Node(a.x*b,a.y*b);
}
向量点积
我们都知道 \(\vec{a}\cdot\vec{b}=|\vec{a}||\vec{b}|\cos\theta\)。其中 \(|\vec{a}|=\sqrt{x_a^2+y_a^2}\)，\(|\vec{b}|=\sqrt{x_b^2+y_b^2}\)，然后根据余弦定理：\(d_{AB}^2=|\vec{a}|^2+|\vec{b}|^2-2|\vec{a}||\vec{b}|\cos\theta=(x_a-x_b)^2+(y_a-y_b)^2\)，可以得到：\(-2|\vec{a}||\vec{b}|\cos\theta=-2x_ax_b-2y_ay_b\) 即 \(|\vec{a}||\vec{b}|\cos\theta=x_ax_b+y_ay_b\)。
所以 \(\vec{a}\cdot\vec{b}=x_ax_b+y_ay_b\)。
如果你是高中生，那应该知道向量的证明方法，\(\vec{a}=x_a\vec{i}+y_a\vec{j},\vec{b}=x_b\vec{i}+y_b\vec{j}\)。然后点乘起来就是 \((x_a\vec{i}+y_a\vec{j})\cdot(x_b\vec{i}+y_b\vec{j})\)，而 \(\vec{i}\cdot\vec{j}=0\)，所以点乘的结果就是 \(x_a\vec{i}\cdot x_b\vec{i}+y_a\vec{j}\cdot y_b\vec{j}=x_ax_b+y_ay_b\)。
double friend operator*(Node a,Node b){
return a.x*b.x+a.y*b.y;
}
向量叉积
我们都知道 \(|\vec{a}\times\vec{b}|=|\vec{a}||\vec{b}|\sin\theta\)。
我们都知道叉积有几何意义就是两条向量围成的三角形面积的两倍。
看张图，可以发现所谓三角形面积的两倍就是图中黄色的面积。而黄色的面积就是 \(x_by_a-x_ay_b\)，由于 \(\vec{a}\times\vec{b}\) 的结果根据右手螺旋定则是负的，所以叉积的表达式就是 \(|\vec{a}\times\vec{b}|=x_ay_b-x_by_a\)。
叉积的用处还是非常多的。
double friend operator^(Node a,Node b){
return a.x*b.y-b.x*a.y;
}
向量旋转
嫖一张图：
展开得：
void rot(double a){
double x1=x*cos(a)-y*sin(a);
double y1=y*cos(a)+x*sin(a);
x=x1;y=y1;
}
直线及定义
先上代码：
struct Line{
Node s,e;Line(){}
Line(Node _s,Node _e){s=_s;e=_e;}
void input(){s.input();e.input();}
pair<int,Node> friend operator&(Line a,Line b){
Node res=a.s;
if(sgn((a.e-a.s)^(b.e-b.s))==0){
if(sgn((b.e-a.s)^(b.e-b.s))==0) return mkp(0,res);
else return mkp(1,res);
}
double tmp=((a.s-b.s)^(b.e-b.s))/((a.e-a.s)^(b.s-b.e));
res.x+=(a.e.x-a.s.x)*tmp;res.y+=(a.e.y-a.s.y)*tmp;return mkp(2,res);
}
};
就用一条直线的两个点来表示这条直线。同时一般用来表示线段，就是线段的两个端点。
求直线交点
首先判断没有交点的情况。那当且仅当两直线平行的时候，会出现一些奇怪的情况。平行直接用叉积判一下就行了。然后判断一下是不是重合的（当然不判也行）。
来张图：
那我们要求 \(E\) 点的坐标，考虑用 \(a_s\) 的向量去加上 \(\vec{a_sE}\)，然后考虑用一个优美的方法来求 \(\vec{a_sE}\)，\(\vec{a_sE}=\dfrac{|\vec{a_sE}|}{|\vec{a_sa_e}|}\vec{a_sa_e}\)，所以直接求 \(\vec{a_sE}\) 和 \(\vec{a_sa_e}\) 的模长的比值就行了。这个比值我们可以用 \(\dfrac{S_{\vartriangle a_sb_sb_e}}{S_{a_sb_ea_eb_s}}\) 来表示，那分子面积的两倍可以用 \(\vec{b_sa_s}\times\vec{b_sb_e}\) 以及分母的面积两倍可以用 \(\vec{a_sa_e}\times\vec{b_eb_s}\) 来表示，这样上下一比两倍就消掉了。
pair<int,Node> friend operator&(Line a,Line b){
Node res=a.s;
if(sgn((a.e-a.s)^(b.e-b.s))==0){
if(sgn((b.e-a.s)^(b.e-b.s))==0) return mkp(0,res);//重合
else return mkp(1,res);//单纯的平行
}//平行
double tmp=((a.s-b.s)^(b.e-b.s))/((a.e-a.s)^(b.s-b.e));
res.x+=(a.e.x-a.s.x)*tmp;res.y+=(a.e.y-a.s.y)*tmp;return mkp(res,2);
}
一些基础操作
欧几里得距离
不用多说。
double dist(Node a,Node b){return sqrt((a-b)*(a-b));}
点到直线距离
为了让这个东西更有普适性，我们可以直接求出点到直线的垂足。然后求这两个点的距离就行了。
点到直线的垂足
我们假设垂足是 \(R\)，那么同样我们只要知道 \(\vec{L_sR}\) 就行了。 然后我们注意到我们对 \(R\) 点的唯一一点依据就是 \(\vec{L_sP}\cdot\vec{L_sL_e}\)。为了知道 \(\vec{L_sR}\)，考虑求出 \(\dfrac{|\vec{L_sR}|}{|\vec{L_sL_e}|}\)，其实就是 \(\dfrac{\vec{L_sP}\cdot\vec{L_sL_e}}{\vec{L_sL_e}^2}\)。然后把 \(\vec{L_sL_e}\) 等比例缩放一下就是 \(\vec{L_sR}\) 了。
Node PTL(Node P,Line L){//Point To Line
Node res=L.s;double tmp=((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
res.x+=(L.e.x-L.s.x)*tmp;res.y+=(L.e.y-L.s.y)*tmp;return res;
}
点到线段距离
多判断一下垂足不在线段上的情况就行啦。直接看 \(tmp\) 的值是不是在 \(0\) 到 \(1\) 之间就可以了。
如果不是的话那一定是和线段首尾中的一个是最近的，两个比较一下就行了。
判断点是否在线段上
判断叉积为零且点到线段两端形成的向量点积为负就行了。
bool IPOS(Node P,Line L){return sgn((L.s-P)^(L.e-P))==0&&sgn((L.s-P)*(L.e-P))<=0;}
//Is Point On Segment
注意到这段代码是包含点在线段端点上的情况的，所以如果要去掉这种情况，把小于等于号变成小于就行了。
判断线段相交
本来觉得这个可以用求个直线交点然后判断交点是否在线段上来代替。最后在 Aizu 教了一法 G 了，大概是掉精的锅，所以还是要写一下快速排斥和跨立实验。
这东西其实是用了一个很耍赖的东西，就是判两次来得到是否相交。首先两条线段各自所形成的矩形一定相交，所以 \(4\) 个顶点各自判断一下。
然后进行跨立。也就是一条线段的两个点被另一条形成的直线分在左右两边，这样可以用叉积一正一负来判断。
所以总共是 \(6\) 个 if 就完了。
bool ISI(Line l1,Line l2){//Is Segment Intersection
if(max(l1.s.x,l1.e.x)<min(l2.s.x,l2.e.x)) return 0;
if(max(l2.s.x,l2.e.x)<min(l1.s.x,l1.e.x)) return 0;
if(max(l1.s.y,l1.e.y)<min(l2.s.y,l2.e.y)) return 0;
if(max(l2.s.y,l2.e.y)<min(l1.s.y,l1.e.y)) return 0;
//上面是快速排斥实验
if(sgn((l2.s-l1.s)^(l1.e-l1.s))*sgn((l2.e-l1.s)^(l1.e-l1.s))>0) return 0;
if(sgn((l1.s-l2.s)^(l2.e-l2.s))*sgn((l1.e-l2.s)^(l2.e-l2.s))>0) return 0;
return 1;//上面是跨立实验
}
关于多边形的计算
我们可以用一个 vector<Node> poly 来保存一个多边形。
计算多边形面积
这里主要利用叉积，我们假设点是逆时针给出（这样叉积是正的）。
那么就随便选一个点（比如说源点）然后直接相邻两个叉积一下就可以了。凹多边形是没关系的，因为凹进去的地方叉积是负的，会抵消掉，看图就明白了。
红色的是负的，绿色的是正的。（根据叉积的正负来推）
double CPA(vector<Node> poly){//Calculate Polygon Area
double ret=0;int n=poly.size();
rep(i,0,n-1) ret+=(poly[i]^poly[(i+1)%n])/2;
return ret;
}
判断点在多边形内
由于回转数太慢还丢精，我们用光线投射法。
为了处理一些特殊情况，直接算当前边两点的纵坐标与给出的点的差就行了。
int IPIP(Node P,vector<Node> poly){//Is Point In Polygon
int cnt=0,n=poly.size();
rep(i,0,n-1){
Node p1=poly[i],p2=poly[(i+1)%n];
if(IPOS(P,Line(p1,p2))) return -1;
int k=sgn((p2-p1)^(P-p1)),d1=sgn(p1.y-P.y),d2=sgn(p2.y-P.y);
if(k>0&&d1<=0&&d2>0) cnt++;
if(k<0&&d2<=0&&d1>0) cnt--;
}return cnt!=0?1:0;
}
判断多边形是否为凸
\(\texttt{A}\color{red}{\texttt{xDea}}\) 说你的板子可以牛逼一点。于是我们让多边形的点顺逆时针给出都行。考虑相邻两边的叉积是否都相同，那么存一下然后看是不是正负都有。
bool IPC(vector<Node> poly){
int n=poly.size();bool vis[3];memset(vis,0,sizeof(vis));
rep(i,1,n-1){
vis[sgn((poly[(i+1)%n]-poly[i])^(poly[(i+2)%n]-poly[(i+1)%n]))+1]=1;
if(vis[0]&&vis[2]) return 0;
}return 1;
}
极角排序
按照一个点的极角排序，第一个点为基准。显然如果叉积是正的，那就是逆时针的，那 \(p_1\) 就应该在 \(p_2\) 的前面
vector<Node> PST(vector<Node> poly){//Polar SorTing
Node p0=poly[0];
sort(poly.begin()+1,poly.end(),[&](Node p1,Node p2){
double crs=(p1-p0)^(p2-p0);
if(sgn(crs)>0) return 1;
else if(sgn(crs)==0&&sgn(dist(p1,p0)-dist(p2,p0))<=0) return 1;
else return 0;
});return poly;
}
\(\color{red}{\text{P.s.}}\) 由于上面这个极角排序是以第一个点为基准点的，在后面介绍闵可夫斯基和的时候所需要的极角排序是必须严格按照极角从小到大排的，所以这里给出一种严格的极角排序方法。实际上也很简单，首先判断是不是在 \(y\) 轴的同一方向，然后再判断叉积就可以了。
sort(ans.begin(),ans.end(),[&](Node a,Node b){
int A=(a.y>0||(a.y==0&&a.x>0));
int B=(b.y>0||(b.y==0&&b.x>0));
if(A!=B) return A>B;return (a^b)>0;
});
凸包相关
多个点，我们找出最外围的点，叫做凸包。
求凸包
有两种求法。Andrew 和 Graham。
Graham
这里先介绍 Graham。
这个算法的流程其实就是一个单调栈一样的东西。首先找到纵坐标最小的点，这个点一定在凸包上，然后以它为基准点做极角排序。然后把接下来一个点入栈，然后每次进来一个点，看图：
如果说栈中是 \(ABC\) 三个点，此时跑进来一个 \(D\)，发现没有构成凸包，判断方法就是 \(\vec{v}\times\vec{a}\le0\)，那么我们就把 \(C\) 弹出，以此类推。
\(\color{red}{\texttt{Important:}}\) 在用 Graham 的时候必须用以第一个点为基准点的极角排序，我们必须保证排的时候前几个点在凸包上！
vector<Node> Graham(vector<Node> poly){
poly=PST(poly);
Node p0=poly[0];int n=poly.size(),id=0;
rep(i,1,n-1)
if(p0.y>poly[i].y||(p0.y==poly[i].y&&p0.x>poly[i].x))
p0=poly[i],id=i;
swap(poly[id],poly[0]);
poly=PST(poly);
vector<Node> ret;
if(n<=2) return poly;
ret.pb(poly[0]);ret.pb(poly[1]);
rep(i,2,n-1){
while(ret.size()>1&&sgn((ret[ret.size()-1]-ret[ret.size()-2])^
(poly[i]-ret[ret.size()-2]))<=0) ret.pop_back();
ret.pb(poly[i]);
}return ret;
}
\(\color{red}{\text{P.s.}}\) 如果有的诡异题需要保留在边上的点，那么就把单调栈那一段的 while 中等号去掉，然后最后倒着扫一遍，把在最后一条边上的点加入到答案里去。如下：
ExCode
vector<Node> Graham(vector<Node> poly){
poly=PST(poly);
Node p0=poly[0];int n=poly.size(),id=0;
rep(i,1,n-1)
if(p0.y>poly[i].y||(p0.y==poly[i].y&&p0.x>poly[i].x))
p0=poly[i],id=i;
swap(poly[id],poly[0]);
poly=PST(poly);
vector<Node> ret;
if(n<=2) return poly;
ret.pb(poly[0]);ret.pb(poly[1]);
rep(i,2,n-1){
while(ret.size()>1&&sgn((ret[ret.size()-1]-ret[ret.size()-2])^
(poly[i]-ret[ret.size()-2]))<0) ret.pop_back();
ret.pb(poly[i]);
}per(i,n-1,0) if(IPOS(poly[i],Line(ret[ret.size()-1],ret[0])))
ret.pb(poly[i]);//注意这里的 IPOS 不包含端点
return ret;
}
Andrew
自己在做闵可夫斯基和的时候发现了一个问题，就是说闵可夫斯基和需要一个直接用极角来排序的极角排序，与用 Graham 时用的以第一个点为基准点的极角排序不一样。况且 Graham 也是必须用这种以凸包上的点为基准点的极角排序。所以为了避免这样繁琐的选择极角排序的过程，我们来看一种不需要极角排序的凸包求法——Andrew。
咕咕咕
旋转卡壳
然后就是这个不会读的算法了。(\(xu\acute{a}n\,zhu\check{a}n\,qi\check{a}\,qi\grave{a}o\))
来点喜闻乐见图：
接下来介绍一下简单的对踵点：
定义：如果过两个点能作一对平行线把整个多边形夹在里面，那么我们说这是一对对踵点。
显然，平面最远点对一定在对踵点中产生。
接下来我们来看怎么判断一对点是不是对踵点。
假如说我们已经求出了凸包（就是说点是逆时针给出的），那我如果想求 \(G\) 的对踵点怎么搞？你看这个 \(B\)，显然不是对吧，那怎么判断呢？欸，\(\vec{GF}\times \vec{BC}<0\)。你看这个 \(C\)，显然是对吧，欸，\(\vec{GF}\times \vec{CD}>=0\)。所以直接用叉积判断一下就可以了。
那在旋转卡壳的时候怎么搞呢？可以用双指针维护一个点的对踵点。然后每次取一对对踵点的最大值。
double Rot_Calip(vector<Node> poly){
double ans=0;Node p;int n=poly.size(),cur=1;
rep(i,0,n-1){ p=poly[i]-poly[(i+1)%n];
while(sgn(p^(poly[(cur+1)%n]-poly[cur]))<0) cur=(cur+1)%n;
ans=max(ans,dis(poly[i],poly[cur]));
}return ans;
}
这样我们就求出了平面最远点对。
那如果是平面最近点对呢？
分治。每次把到 p[mid] 的距离小于上一层分治的答案的点给抠出来做一次 \(O(k^2)\) 的暴力。最终复杂度还是比较优秀的。
闵可夫斯基和
求两个凸包的和，听起来比较奇怪。通俗地讲，就是两个凸包中所有点构成的向量相加所形成的点集。比如把凸包 \(A\) 和 \(B\) 相加，我们只要从原点向 \(A\) 中每个点引一条向量，然后把 \(B\) 顺着这些向量移动所覆盖的点集就是这两个凸包的闵可夫斯基和。
闵可夫斯基和有一些牛逼的性质。（你会证吗？我不会。）
两个凸包的闵可夫斯基和还是一个凸的东西。但是可能会包含一些比如说三个点一条直线的情况，如果有必要还是重新求一次凸包。
这东西怎么求呢？把两个凸包的边极角排序然后顺次相连。具体来说，把相邻两个点构成的向量丢到某个容器里然后极角排序，最后顺次相连就行了。然后我们注意到如果原来就是凸包，点是顺时针给出的，那么我们直接归并就可以了。这样就能提升一点效率了。
对了，上面这个题可以直接把所有凸包求出来然后暴力求闵可夫斯基和，效率比分治加归并的闵可夫斯基和要高。（虽然但是我为了测试归并的闵可夫斯基和还是用后者写了）
\(\color{red}{\texttt{Important:}}\) 闵可夫斯基和一定要用严格的极角排序！
int get_left_down(vector<Node> poly){
Node p0=poly[0];int n=poly.size(),id=0;
rep(i,1,n-1)
if(p0.y>poly[i].y||(p0.y==poly[i].y&&p0.x>poly[i].x))
p0=poly[i],id=i;
return id;
}
vector<Node> Minkowski(vector<Node> A,vector<Node> B){
vector<Node> e,ret;int n=A.size(),m=B.size();
if(!n) return B;if(!m) return A;
rep(i,0,n-1) e.pb(A[(i+1)%n]-A[i]);
rep(i,0,m-1) e.pb(B[(i+1)%m]-B[i]);
int ida=get_left_down(A),idb=get_left_down(B);
ret.pb(A[ida]+B[idb]);
sort(e.begin(),e.end(),[&](Node a,Node b){
int A=(a.y>0||(a.y==0&&a.x>0));
int B=(b.y>0||(b.y==0&&b.x>0));
if(A!=B) return A>B;return (a^b)>0;
});
int tot=e.size();
rep(i,0,tot-1) ret.pb(ret.back()+e[i]);
return ret;
}
这是我写的应该是对的代码。
（坑：但是不知道为什么，我用分治写上面的例题的时候却挂掉了，求蝳蝳～）
