平面计算几何全家桶
平面计算几何全家桶
准备工作
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;
点+向量=点
向量旋转
点积
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{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;
}
