平面计算几何全家桶

平面计算几何全家桶

准备工作

const double eps=1e-8;
const double inf=1e18;
inline int dcmp(double a){return a<-eps?-1:(a>eps?1:0);}//-eps!!! 

计算几何题一般用的是浮点数,所以要格外注意精度问题,封装了dcmp函数后，把所有含比较的地方的变量加上dcmp,如x==0换成dcmp(x)==0

点与向量

向量的线性运算

struct vec{
	double x,y;
	vec(){}
	vec(double _x,double _y){x=_x;y=_y;}
	friend vec operator + (vec p,vec q){return vec(p.x+q.x,p.y+q.y);}
	friend vec operator - (vec p,vec q){return vec(p.x-q.x,p.y-q.y);} 
	friend vec operator * (vec p,double k){return vec(p.x*k,p.y*k);}
	friend vec operator / (vec p,double k){return p*(1.0/k);}
};

本质上点和向量用的是同一个类型，但还是定义两个名字清晰一些

typedef vec point;

点+向量=点

向量旋转

点积

$$\bm{a} \cdot \bm{b}= |\bm{a}||\bm{b}|\cos <\bm{a},\bm{b}>=x_1x_2+y_1y_2$$

double dot(vec p,vec q){return p.x*q.x+p.y*q.y;}
double len(vec p){return sqrt(dot(p,p));}
double dist(point P,point Q){return sqrt(dot(P-Q,P-Q));}
double dist2(point P,point Q){return dot(P-Q,P-Q);}//距离的平方

$|\bm{a}|^2=\bm{a} \cdot \bm{a}$,因此用点积表示向量长度与两点距离

用点积可求向量投影

叉积

$$\bm{a} \times \bm{b}= |\bm{a}||\bm{b}|\sin <\bm{a},\bm{b}>=x_1y_2-x_2y_1$$

注意叉积不满足交换律$\bm{a} \times \bm{b} = - (\bm{b} \times \bm{a})$

二维向量叉积的结果是一个实数,其几何意义是$\bm{a} ,\bm{b}$围成的平行四边形的面积(有符号)。符号的判定可用右手定则:若 P × Q > 0, 则P在Q的顺时针方向(拇指向上)。 若 P × Q < 0, 则P在Q的逆时针方向。

double cross(vec p,vec q){return p.x*q.y-p.y*q.x;}

直线

直线方程

在高中我们学习了五种直线方程的表示方法，计算几何中常用两点表示直线(射线/线段),还可以用点向式:

$Q=P+ t \bm{v}$

其中$P$为直线上一点,$\bm{v}$为直线的方向向量,$Q$为直线上任意一点,$t$为参数。 给出$P,\bm{v}$,则直线唯一确定

其实直线的参数方程$\begin{cases} x=x_0+t\cos \theta \\ y=y_0+t\sin \theta \end{cases}$就是这种表示法的特殊情况

struct line{
	point P;vec v;//P+tv
	line(){}
        line(point A,point B){P=A,v=B-A;}
	line(point _P,vec _v){P=_P;v=_v;}
};

直线的交点

如何用点向式方程求直线的交点?

设两直线分别为 $(P,\bm{v}),(Q,\bm{m})$ ,交点$A = P + t \bm{v}$.

又因为$A$ 在另一条直线上, 故$\vec{QA}\times \bm{m}=0$

$(P + t \bm{v} - Q) \times \bm{m}=0$

叉积关于向量加法运算具有分配律，展开后得

$t (\bm{v} \times \bm{m}) = (Q-P)\times \bm{m}$

解的$t = \frac{(Q-P)\times \bm{m}}{\bm{v} \times \bm{m}}$

故 $$A=P+\frac{(Q-P)\times \bm{m}}{\bm{v} \times \bm{m}}\bm{v}$$

point interSec(line l,line m){
	return l.P+l.v*(cross(m.P-l.P,m.v)/cross(l.v,m.v)); //((P-Q)+tv)m=0，解出t,交点为P+tv 
}

注意两直线平行时无交点，该函数会返回(nan,inf),因此要特判

点、线段、直线的位置关系

位置关系的判断主要依赖于点积和叉积的符号性质。如利用叉积可以判断"上下"关系(指沿法向),点积可判断"左右"(指沿投影的方向)。叉积为0表示平行，点积为0表示垂直

请读者自行画图观察(有空再补上，先放个板子。

点与直线

点P在直线AB上:用叉积判断，若共线则叉积为0

bool isOnLine(point P,point A,point B){return dcmp(cross(P-A,B-A))==0;}

点P在直线AB上的投影(垂足)。

设垂足为H,$\vec{AH}=t \vec{AB},\vec{PH}=\vec{AP}-\vec{AH},\vec{PH}\cdot \vec{AB}=0$,解出$t$即可
(用长度算要根号，误差太大)

point footPnt(point P,point A,point B){
	vec AP=P-A,AB=B-A;
	return A+AB*(dot(AP,AB)/dot(AB,AB));
}

点P在AB上的对称点:关于垂足对称即可.

point symPnt(point P,point A,point B){
	return P+(footPnt(P,A,B)-P)*2; 
}

线段

点P在线段AB上。(叉积为0表示共线,点积为负数表示位置在中间)

bool isOnSeg(point P,point A,point B){return dcmp(cross(P-A,B-A))==0&&dcmp(dot(P-A,P-B))<=0;}

直线AB与线段CD相交

bool isCrossLS(point A,point B,point C,point D){
return dcmp(cross(B-A,C-D))!=0&&isOnSeg(interSec(line(A,B),line(C,D)),C,D);}//注意排除平行的 

线段AB与线段CD相交

多边形

求凸包

用Graham算法(本质是单调栈)

void convexHull(point *a,int n){
	static point s[MAXN+5];
    for(int i=1;i<=n;i++) if(a[i].x<a[1].x||(a[i].x==a[1].x&&a[i].y<a[1].y)) swap(a[i],a[1]);
    sort(a+2,a+1+n,[=](point p,point q)->bool{//极角排序 
        double ang=cross(p-a[1],q-a[1]);
        if(fabs(ang)<eps) return dist2(p,a[1])<dist2(q,a[1]);
        else return ang>eps;
    });
	int top=0;
    for(int i=1;i<=n;i++){//单调栈 
   		while(top>1&&cross(s[top]-s[top-1],a[i]-s[top-1])<=eps) top--;
   		s[++top]=a[i]; 
    }
}

判断点是否在多边形内

任意多边形:射线法

从某点出发任意引出一条射线，若射线与多边形的边相交奇数次则该点在多边形内。否则在多边形外。

int isInPoly(poly &v,point P){
	int cnt=0;
	line l=line(P,P+point(1,0));
	for(int i=0;i<(int)v.size();i++){
		point A=v[i],B=(i+1<(int)v.size())?v[i+1]:v[0];
		if(A.y>B.y) swap(A,B);
		if(isOnSeg(P,A,B)) return 2;//在边上
		if(P.y>=A.y&&P.y<B.y) cnt+=dcmp(interSec(l,line(A,B)).x-P.x)>0; //注意穿过顶点的处理，看两条边是否在射线同侧
	}
	return cnt%2==1;
}

凸多边形:二分法