2017-05-24

【算法笔记】漫话计算几何

本文为博主原创，未经许可不得转载

二维几何基础

向量和点

关于向量，在数学课上早已详细讲过，这里主要说一下向量和点的区别：

点-点=向量
向量+向量=向量
点+向量=点
点+点是没有意义的

在代码实现上，两者都是同样的定义，但是在题目分析时，一定要区分开来

下面是它们的常用定义：

注：极角指从ｘ轴正半轴旋转到该向量方向所需要的角度

struct Point{
	double x,y;
	Point(double a=0,double b=0):x(a),y(b){}
	inline void read(){scanf("%lf%lf",&x,&y);}//读入
	inline void print(){printf("%.6lf %.6lf ",x,y);}//输出
	inline Point operator+(const Point &b){return Point(x+b.x,y+b.y);}//向量加法
	inline Point operator-(const Point &b){return Point(x-b.x,y-b.y);}//向量减法
	inline Point operator*(const double p){return Point(x*p,y*p);}//数乘向量
	inline Point operator/(const double p){return Point(x/p,y/p);}//数除向量
	inline bool operator==(const Point &b)const{return fabs(x-b.x)<eps&&fabs(y-b.y)<eps;}//向量相等
	inline bool operator< (const Point &b)const{return x<b.x||(x==b.x&&y<b.y);}//排序
	inline double Angle(){return atan2(y,x);}//向量极角
};

基本运算

点积（定义式）：$v·w=|v||w|cos\theta$

点积（坐标式）：$v·w=v_x w_x + v_y w_y$

利用点积可以计算向量长度和两向量夹角：

1
2
3

inline double Dot(Point A,Point B){return A.x*B.x+A.y*B.y;}//向量点积
inline double Lenth(Point A){return sqrt(Dot(A,A));}//向量长度
inline double Angle(Point A,Point B){return acos(Dot(A,B)/Lenth(A)/Lenth(B));}//向量夹角

叉积（定义式）：$v\times w=|v||w|sin\theta$

叉积（坐标式）：$v\times w=v_x w_y-v_y w_x$

叉积的几何意义：两个向量$v$和$w$的叉积表示$v$和$w$组成的平行四边形的有向面积

叉积的性质：$v\times w=-w\times v$

1	inline double Cross(Point A,Point B){return A.xB.y-B.xA.y;}//向量叉积

向量旋转：（$\alpha$为逆时针旋转的角度）

$x’=x\cos\alpha -y\sin\alpha$

$y’=x\sin\alpha +y\cos\alpha$

证明：
　　

　　
由上图可知：
　　
$\frac{x}{\cos\beta}=\frac{y}{\sin\beta}=\frac{x’}{\cos(\alpha+\beta)}=\frac{y’}{\sin(\alpha+\beta)}=R$
　　
$x’=R\cos(\alpha+\beta)=R\cos\alpha\cos\beta-R\sin\alpha\sin\beta=x\cos\alpha -y\sin\alpha$
　　
$y’=R\sin(\alpha+\beta)=R\sin\alpha\cos\beta+R\cos\alpha\sin\beta=x\sin\alpha +y\cos\alpha$

1	inline Point rotate(double rad){return Point(xcos(rad)-ysin(rad),xsin(rad)+ycos(rad));}//向量旋转

单位法线：左旋90度，然后把长度化为1即可

1	inline Point Normal(Point A){return Point(-A.y/Lenth(A),A.x/Lenth(A));}

点和直线

直线的参数方程

除了我们熟知的一般式$Ax+By+C=0$外，还有一种在高中数学中比较冷门的表示方法：参数式

但是在计算几何中，由于参数式与向量密切相关，所以应用广泛

直线可以用直线上一点$P_0$和方向向量$v$表示，记做：$P_0+vt$ ($t$为参数)

如果已知直线上两点$A$和$B$，则参数方程就是$A+(B-A)t$

参数方程最方便的地方在于直线、射线、线段的方程形式是一样的，区别仅仅在于参数$t$的值域

直线：没有限制

射线：$t>0$

线段：$0<t<1$

直线的交点

设两直线分别为$P+tv$和$Q+tw$，交点在第一条直线的参数为$t_1$，第二条为$t_2$，则：

$P_x+t_1 v_x=Q_x+t_2 w_x$

$P_y+t_1 v_y=Q_y+t_2 w_y$

联立可得：$(v_x w_y-w_x v_y)t_1=w_x(P_y-Q_y)-w_y(P_x-Q_x)$

设$u=P-Q$可得：$t_1=\frac{Cross(w,u)}{Cross(v,w)}$

inline Point Intersection(Point P,Point v,Point Q,Point w){//直线交点
	Point u=P-Q;
	double t1=Cross(w,u)/Cross(v,w);
	return P+v*t1;
}

点到直线的距离

直接用叉积计算，即用平行四边形的面积除以底长

1	inline double Line_Distance(Point A,Point B,Point P){return fabs(Cross(B-A,P-A)/Lenth(B-A));}//点到直线距离

点到线段的距离

设投影点为$Q$，如果$Q$在线段$AB$上，则所求距离就是该点到直线$AB$的距离

如果$Q$在射线$BA$上，则所求距离为$PA$距离

否则就是$PB$距离

判断$Q$的位置时用点积进行

inline int dcmp(double x){
	if(fabs(x)<eps) return 0;
	else return x<0?-1:1;
}
inline double Setment_Distance(Point A,Point B,Point P){//点到线段距离
	if(A==B) return Lenth(P-A);
	else if(dcmp(Dot(B-A,P-A))<0) return Lenth(P-A);
	else if(dcmp(Dot(A-B,P-B))<0) return Lenth(P-B);
	else return Line_Distance(A,B,P);
}

求投影点

求点$P$在直线$AB$上的投影点$Q$

设直线$AB$的参数方程为$A+vt$，点$Q$所对应的参数为$t_0$，那么$Q=A+vt_0$

因为$PQ$垂直$AB$，所以$Dot(v,P-(A+vt_0))=0$

根据分配率展开得到：$t_0=\frac{Dot(v,P-A)}{Dot(v,v)}$

inline Point Projection(Point A,Point B,Point P){//点在直线上的投影点
	Point v=B-A;
	return P+v*(Dot(v,P-A)/Dot(v,v));
}

线段相交判定

线段相交的充要条件：每条线段的两个端点都在另一条线段的两侧

用叉积判断即可

1
2
3

bool Segment_Jiao(Point A,Point B,Point C,Point D){
	return dcmp(Cross(B-A,C-A)*Cross(B-A,D-A))<0&&dcmp(Cross(D-C,A-C)*Cross(D-C,B-C))<0;
}

判断点是否在线段上

点在直线上的充要条件：叉积等于$0$，点积小于$0$

1
2
3

bool On_Segment(Point P,Point A,Point B){
	return dcmp(Cross(A-P,B-P))==0&&dcmp(Dot(A-P,B-P))<0;
}

二维几何常用算法

求多边形面积

我们可以从起点出发把多边形分成$n-2$个三角形，然后把面积加起来

三角形的面积用叉积计算即可

值得一提的是这个方法对非凸多边形同样适用

因为我们计算的是三角形的有向面积，在外面的部分可以抵消掉

inline double Area(){
	double ret=0;
	for(int i=2;i<n;++i) ret+=Cross(A[i]-A[1],A[i+1]-A[1]); 
	return ret/2;
}

点在多边形内判定

我们作出一条从该点出发水平向右的射线，统计多边形穿过这条射线正反多少次

我们把这个次数记做绕数$wn$，逆时针穿过时，$wn++$，顺时针穿过时，$wn–$

若$wn!=0$就说明点在多边形内

具体实现时，用叉积判断穿过方向，用坐标判断是否穿过

inline int check(){
	int wn=0;
	for(int i=1;i<=n;++i){
		if(On_Segment(P,A[i],A[i+1])) return -1;//在边界上
		int k=dcmp(Cross(A[i+1]-A[i],P-A[i]));
		int d1=dcmp(A[i].y-P.y);
		int d2=dcmp(A[i+1].y-P.y);
		if(k>0&&d1<=0&&d2>0) ++wn;
		if(k<0&&d2<=0&&d1>0) --wn;
	}
	return wn?1:0;
}

凸包

凸包是把给定点包围在内部的，面积最小的凸多边形

我们一般使用$Andrew$算法求凸包

首先把点按照$x$从小到大排序，采用在线加点的方法，使用单调栈维护上凸包和下凸包

然后把上凸包和下凸包合并起来就是完整的凸包

由于有排序操作，所以时间复杂度为$O(n\log n)$

int ConvexHull(){
	int tail=0;
	sort(A+1,A+cnt+1);
	for(int i=1;i<=cnt;++i){
		while(tail>1&&Cross(q[tail]-q[tail-1],A[i]-q[tail-1])<=0) --tail;
		q[++tail]=A[i];
	}
	int lim=tail;
	for(int i=cnt-1;i;--i){
		while(tail>lim&&Cross(q[tail]-q[tail-1],A[i]-q[tail-1])<=0) --tail;
		q[++tail]=A[i];
	}
	return tail;
}

半平面交

半平面交问题就是给定若干个半平面，求他们的公共部分

半平面一般用有向线段表示，它的左侧就是它所代表的半平面

有向直线的定义如下：

struct Line{
	Point P,v;  double ang;  Line(){} 
	Line(Point P,Point v):P(P),v(v){ang=atan2(v.y,v.x);}
	inline bool operator< (const Line &b)const{return ang<b.ang;}
};

求半平面交的方法与凸包类似，先按照极角排序，然后采用在线添加的方式，用双端队列维护

时间复杂度也是$O(n\log n)$，具体细节见代码

bool Halfplane_Jiao(){
	sort(L+1,L+n+1);  int head=1,tail=0; q[++tail]=L[1];
	for(int i=2;i<=n;++i){
		while(head<tail&&!On_Left(P[tail-1],L[i])) --tail;
		while(head<tail&&!On_Left(P[head],L[i])) ++head;
		q[++tail]=L[i];
		if(fabs(q[tail].ang-q[tail-1].ang)<eps){
			--tail;
			if(On_Left(L[i].P,q[tail])) q[tail]=L[i];
		}
		if(head<tail) P[tail-1]=Line_Jiao(q[tail-1],q[tail]);
	}
	while(head<tail&&!On_Left(P[tail-1],q[head])) --tail;
	return tail-head>1;
}

旋转卡壳

旋转卡壳算法用于求对踵点

对踵点：若分别经过$P_i$和$P_j$存在两条平行直线，将整个凸包卡在中间，则称$P_i$和$P_j$为一对对踵点

我们逆向思考一下：

如果$P_i,P_j$是凸包上最远两点，必然可以分别过$P_i,P_j$画出一对平行线

通过旋转这对平行线，我们可以让它和凸包上的一条边重合，假设这条边为$P_j P_k$

可以注意到，$P_i$是凸包上离直线$P_j P_k$最远的点

于是我们的思路就是枚举凸包上的所有边，对每一条边找出凸包上离该边最远的顶点

计算这个顶点到该边两个端点的距离，并记录最大的值

直观上这是一个$O(n^2)$的算法，和直接枚举任意两个顶点一样了

但是注意到当我们逆时针枚举边的时候，最远点的变化也是逆时针的

这样就可以不用从头计算最远点，而可以紧接着上一次的最远点继续计算

至于判断最远点，用叉积计算对应的三角形面积即可

于是我们得到了$O(n)$的算法，而且代码十分简洁

这里吐个槽：老师看到我在搞旋转卡壳，问我这东西好像得用平衡树维护吧，2333

int pos=2;  A[n+1]=A[1];
for(int i=1;i<=n;++i){
	while(fabs(Cross(A[pos%n+1]-A[i+1],A[i]-A[i+1]))>fabs(Cross(A[pos]-A[i+1],A[i]-A[i+1]))) 
		pos%=n,++pos;
	//此时，i和pos就是一对对踵点
}

参考文献

《算法竞赛入门经典——训练指南》——刘汝佳
《计算几何初步》——cydiater 传送门
《旋转卡壳算法》——Coding 传送门
《谈谈我对旋转卡壳算法本质的理解》——金字塔传送门

本文标题:【算法笔记】漫话计算几何

文章作者:chty

发布时间:2017年05月24日 - 18时07分

最后更新:2018年07月09日 - 20时18分

许可协议:本文为博主原创，未经许可不得转载。