二维计算几何基础

什么是“二维计算几何”？当我们将需要解决的几何问题的范围限制在二维平面内，这样就用到了二维计算几何。要用电脑解平面几何题？这就是陷入误区了，我们并不是用计算机算数学卷子上的几何题去了，而是解决一些更加复杂的几何相关问题。

为了解决复杂且抽象的问题，我们一定要选择合适的研究方法。对于计算机来说，给它看几何图形……

我们可以把要研究的图形放在平面直角坐标系或极坐标系下，这样解决问题就会方便很多。

计算几何中坐标一般是实数，编程时使用double，不要使用精度较低的float。

在进行浮点数运算时会产生精度误差，可以设置一个偏差值eps(epsilon)来控制精度。

判断浮点数是否等于零或两个浮点数是否相等要用eps辅助判断。

#include<bits/stdc++.h>
using namespace std;
#define db double
const db pi = acos(-1.0);   //高精度圆周率
const db eps = 1e-8;        //可根据具体需要更改精度
inline int sgn(db x){		//判断是否等于零
	if(fabs(x)<eps) return 0;
	else return x<0?-1:1;
}
inline int dcmp(db x,db y){	//比较浮点数
	if(fabs(x-y)<eps)return 0;
	else return x<y?-1:1;
}

点和向量

1. 点

二位平面中的点用(x,y)表示。

struct Point{
	db x,y;
	Point(){}
	Point(db _x,db _y):x(_x),y(_y){}
};

2. 两点之间的距离

db Dist(Point &a,Point &b){
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

3. 向量

有大小，有方向的量称为向量(矢量) ，只有大小没有方向的量称为标量。

在平面上两个点可确定一个向量，由于向量不是一个有向线段，及向量没有起点与终点的概念，所有向量平移后不变。大多时候为了方便可直接将向量平移到原点处。向量的表示在形式上与点一样，可以用点的数据结构表示向量。

typedef Point Vector;

4. 向量的运算

在struct Point中，对向量运算重载运算符。

+：点加点没有意义，点加向量得到另一个点，向量加向量得到另一个向量。

Point operator+(Point &a){return Point(x+a.x,y+a.y);}

-：两个点的差得到一个向量，向量A减B得到有B指向A的向量。

Point operator-(Point &a){return Point(x-a.x,y-a.y);}

*：向量乘实数等比例放大。

Point operator*(db a){return Point(x*a,y*a);}

/:向量与实数相除等比例缩小。

Point operator/(db a){return Point(x/a,y/a);}

==:相等

bool operator==(Point &a){return sgn(x-a.x)==0&&sgn(y-a.y)==0;}

点积和叉积

向量的基本运算是点积和叉积，计算几何的各种操作几乎都基于这两种操作。

1. 点积(Dot product)

其中θ为A，B之间的夹角。点积的几何意义是A在B上投影的长度乘以B的模长。

在编程中并不需要知道θ。如果已知A=(A.x,A.y),B=(B.x,B.y)，那么有：

A⋅B=A.x∗B.x+A.y∗B.y

A.x∗B.x+A.y∗B.y

=(|A|cosθ1∗|B|cosθ2)+(|A|sinθ1∗|B|sinθ2)

|A||B|(cosθ1cosθ2+sinθ1sinθ2)

|A||B|cos(θ1−θ2)

|A||B|cosθ

db Dot(Vector &a,Vector &b){return a.x*b.x+a.y*b.y;}

2. 点积的应用

（1）判断A,B之间的角是钝角锐角还是直角。

（2）求A的长度(模)

db len(Vector &a){return sqrt(Dot(a,a));}

（3）求A,B夹角的大小

db Angle(Vector &a,Vector &b){return acos(Dot(a,b)/len(a)/len(b));}

3. 叉积(Cross product)

叉积是比点积等常用的集合概念。

θ表示A旋转到B所经过的夹角。

两个向量的叉积的结果是一个带符号的数值，其几何意义为A,B形成平行四边形的“有向”面积，其正负可通过“右手定则”判断。

db Cross(Vector &a,Vector &b){return a.x*b.y-a.y*b.x;}

A.x∗B.y−A.y∗B.x

=|A||B|(cosθ1sinθ2−sinθ1cosθ2)

=|A||B|sin(θ2−θ1)

当θ2>θ1时，结果为正，且B在A的逆时针方向。

故可通过正负来判断A，B的相对位置。

4. 叉积的基本应用

（1）判断A，B的方向关系。

（2）计算两向量构成的平行四边形的“有向面积”

db Area2(Point a,Point b,Point c){return Cross(b-a,c-a);}

（3）计算3个点构成的三角形面积

（4）向量旋转

是向量(x,y)绕起点逆时针旋转θ，旋转后的向量为(x′,y′)

Vector Rotate(Vector a,db rad){
	return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}

（5）用叉积检查向量是否平行或重合

bool Parallel(Vector a,Vector b){return sgn(Cross(a,b))==0;}

点和线

1. 直线的表示

直线的表示方法很多，编程时可灵活选择：

（1）用直线上的两个点表示。

（2）ax+by+c=0，普通式。

（3）y=kx+b，斜街式，注意斜率不存在的情况。

（4）P=P0+vt，点向式。用P0和v来表示直线P，t是任意值，P0是直线上一点，v是方向向量，v=A−B。

t无限制时，P是直线。

t∈[0,1]，P是A，B之间的线段。

t>0,P是射线。

struct line{
	Point p1,p2;
	line(){}
	line(Point _p1,Point _p2):p1(_p1),p2(_p2){}
	//一个点和角(angle)确定直线，0<=angle<pi
	line(Point p,db angle){
		p1=p;
		if(sgn(angle-pi/2)==0){p1=p+Point(0,1);}
		else p1=p+Point(1,tan(angle));
	}
	//ax+by+c=0;
	line(db a,db b,db c){
		if(sgn(a)==0){
			p1=Point(0,-c/b);
			p2=Point(1,-c/b);
		}
		else if(sgn(b)==0){
			p1=Point(-c/a,0);
			p2=Point(-c/a,1);
		}
		else{
			p1=Point(0,-c/b);
			p2=Point(-c/a,0);
		}
	}
};

2. 线段的表示

可以用两个点表示线段，直接用直线的数据结构表示线段。

typedef line Segment;

3. 在二维平面上，点和直线有三种位置关系，在直线上，左侧和右侧。用叉积的正负可判断方向。

int Point_line_relation(Point p,line v){
	int c=sgn(Cross(p-v.p1,v.p2-v.p1));
	if(c<0)return 1;		//顺时针
	else if(c>0)return 2;	        //逆时针
	return 0;			//在直线上
}

4. 点和线段的关系

先用叉积判断是否共线，再用点积判断点与线段两端的角是否为180°。

int Point_on_seg(Point p,line v){
	return sgn(Cross(v.p1-p,v.p2-p))==0&&sgn(Dot(v.p1-p,v.p2-p))<=0;
}

5. 点到直线的距离

已知点p到直线v(p1,p2)，求p到v的距离。首先用叉积求p,p1,p2构成的平行四边形的面积，然后除以线段(p1,p2)的长度即可。

db Dis_point_line(Point p,line v){
    return fabs(Cross(p-v.p1,v.p2-v.p1)/Dist(v.p1,v.p2));//注意叉积有正负
}

6. 点在直线上的投影

已知直线上的两点p1，p2以及直线外一点p,求投影点p0，我们先来看下上图左边的那种情形。

令，即k是线段p0p1和p2p1长度的比值，那么有。

根据点积的概念有，

及

带入即：

那么：

容易证明上图右边那种情况也满足此式。

Point Point_line_proj(Point p,line v){
	double k=Dot(p-v.p1,v.p2-p)/len(v.p2-v.p1)/len(v.p2-v.p1);
	return v.p1+(v.p2-v.p1)*k;
}

7. 点关于直线的对称点

求点p关于直线的对称点，先求投影点，在求对称点。

Point Point_line_symmetry(Point p,line v){
	Point p1=Point_line_proj(p,v);
	return Point(p1.x*2-p.x,p1.y*2-p.y);
}

8. 点到线段的距离

点p到线段AB的距离，在p到A的距离，p到B的距离，p到直线AB的距离(如果垂足在线段AB上)。

db Dis_point_seg(Point p,Segment v){
	if(Dot(p-v.p1,v.p2-v.p1)<0||Dot(p-v.p2,v.p1-v.p2)<0)return min(Dist(p,v.p1),Dist(p,v.p2));
	return Dis_point_line(p,v);
}

9. 两条直线的位置关系

int line_relation(line v1,line v2){
	if(sgn(Cross(v1.p1-v1.p2,v2.p1-v2.p2))==0){
		if(Point_line_relation(v1.p1,v2)==0)return 1;//重合
		else return 0;//平行
	}
	return 2;//相交
}

10. 求两个直线的交点

可以联立方程求解，也可以用简单的叉积来实现。

如图，AB与CD的交点为P。

容易得到以下两个等式：

联立上面两个方程，即可得到点P的坐标：

Point Cross_point(line a,line b){
	db s1=Cross(b.p1-a.p1,a.p2-a.p1);
	db s2=Cross(a.p2-a.p1,b.p2-a.p1);//注意叉积的正负
	return Point(s1*b.p2.x+s2*b.p1.x,s1*b.p2.y+s2*b.p1.y)/(s1*s2);
}

注意函数要对(s1+s2)做除法，在调用前要保证两直线不平行且不重合。

11. 判断两线段是否有交点

bool Cross_segment(Segment a,Segment b){
	db c1=Cross(a.p1-a.p2,b.p1-a.p2),c2=Cross(a.p1-a.p2,b.p2-a.p2);
	db d1=Cross(b.p1-b.p2,a.p1-b.p2),d2=Cross(b.p1-b.p2,a.p2-b.p2);
	return c1*c2<=0&&d1*d2<=0;
}

12. 求两个线段的交点

先判断有没有交点，有交点转化为直线求就行了。

多边形

1. 判断点在多边形内部

给定一个点P和一个多边形，判断点P在多边形的内部还是外部。

从点P引出一条直线(与多边形交于两点)，检查P与多边形每条边的相交情况，检查的方向不影响结果，但不能中途改变方向。

检查以下3个参数：

显然u,v是用来判断是否相交，c用来判断p与边的位置关系。

最后当num>0时p在多边形内部。

int Point_in_polygon(Point pt,Point *p,int n){
	//注意P[]中的点为严格逆时针或顺时针
	for(int i=0;i<n;++i){
		if(p[i]==pt)return 3;//点p在多边形的顶点上
	}
	for(int i=0;i<n;++i){
		line v=line(p[i],p[(i+1)%n]);
		if(Point_on_seg(pt,v))return 2;//点在多边形边上
	}
	int num=0;
	for(int i=0;i<n;++i){
		int j=(i+1)%n;
		int c=sgn(Cross(pt-p[j],p[i]-p[j]));
		int u=sgn(p[i].y-pt.y);
		int v=sgn(p[j].y-pt.y);
		if(c>0&&u<0&&v>=0)num++;
		if(c<0&&u>=0&&v<0)num--;//注意这里不能用n*v<0判断
	}
	return num>0;//1在内部，0在外部 
}

3. 求多边形面积

对于一个多边形，可以在其内部找一点p，将多边形划分为n个三角形，计算n个三角形的面积之和即可。

实际上p的位置可任意选取，叉积的正负可将多余的面积抵消掉，所以一般将p选在原点。

db Polygon_area(Point*p,int n){
	db ans=0;
	for(int i=0;i<n;++i){
		ans+=Cross(p[i],p[(i+1)/n]);
	}
	return ans/2;//注意面积有正负，可根据需要判断是否要去绝对值
}