计算几何

• 二维计算几何
  • 基本函数
    • eps
    • sign
    • myabs
    • mysqr
    • Point/Vector
    • 向量内积
    • 向量叉积
    • 向量旋转
    • Line
    • 判断线段相交
    • 求直线交点
    • 点到线的最短距离
    • 凸多边形面积
  • 基本算法
    • 排序算法
      • 极角排序
      • 水平序排序
    • 二维凸包
      • 定义
      • 实现
    • 旋转卡壳
      • 定义
      • 实现
    • 闵可夫斯基和
      • 定义
      • 实现
    • 半平面交
      • 定义
      • 实现
• 三维计算几何
  • 基本概念
    • 直线
    • 平面
    • 夹角
      • 两直线夹角.
      • 直线与平面的夹角
    • 交点
  • 基本定理
    • 三余弦定理(最小角定理)
    • 三正弦定理(最大角定理)

二维计算几何

参考:https://www.luogu.com.cn/blog/command-block/ji-suan-ji-he-suan-fa-hui-zong.

这一部分比较常用,知识点也比较简单,会主要聚焦于写法.

基本函数

eps

const double eps=1e-10

浮点数计算存在误差,因此大部分时候不能直接进行大小是否相等的判断,我们需要通过规定精度来减少这种误差.

sign

inline int sign(double x){
    if(x<eps&&-x<eps)return 0;
    if(x>eps)return 1;
    return -1;
}

用于判断一个浮点数是正数还是负数.

myabs

inline double myabs(double x){
    if(sign(x)==-1)return -x;
    return x;
}

用来求绝对值.

mysqr

inline double mysqr(double x){
    return x*x;
}

用来求平方.

Point/Vector

struct Point{
    double x,y;
    Point(double X,double Y){
        x=X;y=Y;
    }
    inline double len(){//判断到原点的距离
        return std::sqrt(x*x+y*y);
    }
    inline int que(){//判断点的象限
        if(sign(x)== 1&&sign(y)== 1)return 1;
        if(sign(x)==-1&&sign(y)== 1)return 2;
        if(sign(x)==-1&&sign(y)==-1)return 3;
        return 4;
    }
};
inline Point operator +(const Point &A,const Point &B){
    return Point(A.x+B.x,A.y+B.y);
}
inline Point operator -(const Point &A,const Point &B){
    return Point(A.x-B.x,A.y-B.y);
}
inline Point operator *(const double &c,const Point &A){
    return Point(c*A.x,c*A.y);
}
inline Point operator *(const Point &A,const double &c){
    return Point(c*A.x,c*A.y);
}
inline Point operator /(const Point &A,const double &c){
    return Point(A.x/c,A.y/c);
}
inline double dist(const Point &A,const Point &B){
    return std::sqrt(mysqr(A.x-B.x)+mysqr(A.y-B.y));
}

存储基本的点的信息.也可以看作一个二维向量.

向量内积

inline double operator *(const Point &A,const Point &B){
	return A.x*B.x+A.y*B.y;
}

也就是\(\vec a\sdot \vec b=|\vec a||\vec b|\cos \theta=x_ax_b+y_ay_b\).也就等于\(\vec a\)在\(b\)上的投影与\(\vec b\)的模长的乘积.

内积可以用来判断夹角:

1. 如果\(\vec a\sdot \vec b=0\),则说明\(\vec a\bot \vec b\).
2. 如果\(\vec a\sdot \vec b>0\),则说明\(\vec a\)和\(\vec b\)正方向的夹角小于\(90\degree\).
3. 如果\(\vec a\sdot \vec b<0\),则说明\(\vec a\)和\(\vec b\)正方向的夹角大于\(90\degree\).

两个向量同时旋转相同角度,其内积结果不变.

向量叉积

double operator ^(const Point &A,const Point &B){
    return A.x*B.y-B.x*A.y;
}

也就是\(\vec a\times\vec b=x_ay_b-y_ax_b\).也就等于\(\vec a,\vec b\)两个向量张成的平行四边形(有向)的面积.

叉积可以用来判断方向:

1. 如果\(\vec a\times \vec b=0\),说明二者共线.
2. 如果\(\vec a\times \vec b<0\),说明从\(\vec a\)到\(\vec b\)的方向是顺时针.
3. 如果\(\vec a\times \vec b>0\),说明从\(\vec a\)到\(\vec b\)的方向是逆时针.

向量旋转

//插入Point内部
inline Point rotate(double theta){
    double st=std::sin(theta),double ct=std::cos(theta);
    return Point(x*ct-A.y*st,A.x*st+A.y*ct);
}

也就是将这个竖向量乘左乘旋转矩阵\(\begin{bmatrix}\cos \theta&-\sin \theta\\\sin \theta &\cos \theta\end{bmatrix}\).

Line

struct Line{
    Point a,b;
    double len(){
        return dist(a,b);//计算线段长度,直线和射线基本没用
    }
};
bool onl0(const Point &A,const Line &L){
    return sign((A-L.a)^(A-L.b))==0;
}//判断点是否在直线上
bool onl2(const Point &A,const Line &L){
    return (sign((A-L.a)^(A-L.b))==0)&&
           (sign(A.x-std::min(L.a.x,L.b.x))>=0)&&
           (sign(A.x-std::max(L.a.x,L.b.x))>=0)&&
           (sign(A.y-std::min(L.a.y,L.b.y))>=0)&&
           (sign(A.y-std::max(L.a.y,L.b.y))<=0);
        
}//判断点是否在线段上

用来维护直线/射线/线段的基本信息.

判断线段相交

bool isinter(const Line &L1,const Line &L2){
    if(std::max(L1.a.x,L1.b.x)<std::min(L2.a.x,L2.b.x))return 0;
    if(std::max(L2.a.x,L2.b.x)<std::min(L1.a.x,L1.b.x))return 0;
    if(std::max(L1.a.y,L1.b.y)<std::min(L2.a.y,L2.b.y))return 0;
    if(std::max(L2.a.y,L2.b.y)<std::min(L1.a.y,L1.b.y))return 0;
    if(((L1.a-L2.b)^(L2.a-L2.b))*((L1.b-L2.b)^(L2.a-L2.b))<0)return 0;
    if(((L2.a-L1.b)^(L1.a-L1.b))*((L2.b-L1.b)^(L1.a-L1.b))<0)return 0;
    return 1;
}

前四行被称为快速排斥实验,后两行被称为跨立实验,也就是相交的两线段,对于其中一个线段而言,它的两个端点必然在另一个线段的两侧.但我们注意到共线线段也会通过跨立实验,因此拿快速排斥实验来判掉这种情况.

注意,线段交点不同于直线交点,两条共线的交点是有可能在端点上的.这个需要特判一下.

当然一般判断线段是否相交可以直接求直线交点,然后判断在不在线段上.

求直线交点

Point inter(const Line &L1,const Line &L2){
    double ls=(L1.b-L1.a)^(L2.a-L1.a);
    double rs=(L1.b-L1.a)^(L2.b-L1.a);
    return L2.a+(L2.b-L2.a)*ls/(ls-rs);
}

原理在于,画一下\(x_{L_1},x_{L_2},y_{L_1},y_{L_2}\)围成的四边形,计算面积后用等高不等底计算.注意\(ls\)和\(rs\)所代表的面积一正一负(不一定谁是正),因此需要减一下.

点到线的最短距离

inline double disl0(const Point &A,const Line &L){
    return myabs(((L.a-A)^(L.b-A))/dis(L.a,L.b));
}
inline double disl2(const Point &A,const Line &L){
    if(sign((A-L.a)*(L.b-L.a))<0)return dist(A,L.a);
    if(sign((A-L.b)*(L.a-L.b))<0)return dist(A,L.b);
    return disl0(A,L);
}

直线好做,直接算垂线段长度,用面积除以底长.

线段用点积判一下夹角即可.

凸多边形面积

利用叉乘,任取平面上一点\(O\),则\(S=\frac1 2\sum_{i=1}^n\overrightarrow {OP_i}\times \overrightarrow{OP_{i+1}}\).证明的话考虑分\(O\)在内部和\(O\)在外部两种情况分类讨论.注意此时的\(P\)必须逆时针排列.

另外有皮克定理:在一个平面直角坐标系内,以整点为顶点的简单多边形(任两边不交叉),它内部整点数为\(a\),它的边上(包括顶点)的整点数为\(b\),则它的面积\(S=a+\frac{b}{2}-1\).

基本算法

排序算法

极角排序

定义原点\(O\)并建立坐标系,所有点按照和\(O\)所连直线与\(x\)轴正方向的夹角排序.

极角排序通常使用叉乘来实现,因为叉乘可以快速计算两个向量的方向.但是注意到需要先判断象限,再判断叉乘.

bool cmp(Point A,Point B){//判断A能否在B前面
    A=A-O,B=B-O;//O=(0,0)的时候可以省略
    if(A.qua()!=B.qua()){
        return A.qua()<B.qua();
    }
    if(sign(A^B)==0)return A.x<B.x;
    return sign(B^A)==1;
}

水平序排序

bool cmp(Point A,Point B){
    if(sign(A.x-B.x)==0)return sign(A.y-B.y)==-1;
    return sign(A.x-B.x)==-1;
}

也就是\(x\)相同比\(y\),否则比\(x\).

二维凸包

定义

包住平面上某个点集的周长最小的简单多边形,一定是凸多边形.

实现

用水平序排序,然后从左往右扫一遍得到上凸壳,从右往左扫一遍得到下凸壳.

用一下以下函数:

inline bool anticlock(const Point &A,const Point &B,const Point &C){
    return sign((A-B)^(B-C))==-1;
}
inline bool isline(const Point &A,const Point &B,const Point &C){
    return sign((A-B)^(B-C))==0;
}

前者判断是否\(ABC\)三点是一个上凸的(注意\(ABC\)三点横坐标应该不降).后者判断三点共线.

旋转卡壳

定义

定义凸包上的对踵点对,也就是拿两条平行直线卡着凸包转,这两条直线会卡住凸包的两个点,这些点对组成的集合就是对踵点对集合.

实现

逆时针枚举边,然后看对面有什么被卡住了.由于凸包的凸性,在边逆时针转的时候,另一边的对踵点也会逆时针转.只要找到距离这条边所在直线最远的点即可.但这样的点可能有俩.可以简单特判,也可以加上微小随机扰动量使得这种情况可以忽略.

闵可夫斯基和

一般只讨论凸包的闵可夫斯基和.

定义

两个区域\(A,B\)的闵可夫斯基和定义为\(\{a+b\mid a\in A,b\in B\}\).

实现

事实上,新的区域所形成的凸包,一定是原本\(A,B\)的凸包的边按照某种顺序连接起来得到的结果.

我们考虑旋转一下\(A,B\),使得\(B\)有一条边成为最右边的直上直下的一条边,然后考虑答案区域的最右边的边,这条边一定是\(B\)这个边加上\(A\)的最右边的点.这样这条边必定还在最终的凸包上.就算\(A\)最右边的是一条边,你也会发现最终的凸包最右边也一定是由\(A\)的这条边和\(B\)的这条边拼起来的.

显然,逆时针转一遍整个凸壳,将每条边改为向量(按照逆时针转的方向)然后极角排序,最后顺次链接就是答案.

半平面交

定义

定义半平面为满足\(ax+by+c>0\)或\(ax+by+c\geq 0\)的点对\((x,y)\)组成的集合,感性理解就是一条直线的一侧.

实现

首先直线不好描述左右侧,我们把直线改成向量,这样方便描述左右侧.并不妨假设所有的向量所表示的半平面在向量的左侧.如果两条向量方向相同,则取更靠左的那一条,也就是所在直线截距更大的那个,另一个直接删了就是了.当然你不想删可以把那个废物放前面,这样根据下面的操作过程中它会被弹掉.

我们这么实现:按照上面的顺序一个一个插入,维护一个单调队列.如果前两条向量的交点不在当前这条向量的控制范围,不难发现上一条向量是废物,弹掉它.这样做到最后,队首和队尾可能都有一些废物向量,把它们判掉弹掉即可.

三维计算几何

这一部分知识点比较困难,而几乎用不到,因此只讲简单知识点,就当是高中立体几何知识补档!

基本概念

直线

使用直线的方向向量\(\vec s=(n,m,p)\)和直线上一点\(M_0=(x_0,y_0,z_0)\).那么方程显然为: \[ \frac{x-x_0}{n}=\frac{y-y_0}{m}=\frac{z-z_0}{p} \] 如果换元,我们还有参数方程: \[ \begin{cases}x=x_0+nt\\y=y_0+mt\\z=z_0+pt\end{cases} \]

平面

使用平面上的一点\(P_0(x_0,y_0,z_0)\)和该平面的法向量\(\vec n\)来表示一个平面,不妨设\(\vec n=(A,B,C)\),则该平面的方程显然为: \[ A(x-x_0)+B(y-y_0)+C(z-z_0)=0 \] 如果我们令\(D=-(Ax_0+By_0+Cz_)\),那么平面方程为: \[ Ax+By+Cz+D=0 \]

夹角

两直线夹角.

直接求方向向量的夹角,然后取正值.

对于方向向量分别是\(\vec {s_1}=(n_1,m_1,p_1),\vec{s_2}=(n_2,m_2,p_2)\),也就有\(\varphi=\arccos(\frac{|\vec s_1\sdot \vec s_2|}{|\vec s_1||\vec s_2|})\\\).

直线与平面的夹角

同样使用向量,不妨设方向向量\(\vec s=(n,m,p)\),法向量\(\vec f=(a,b,c)\),那么\(\varphi=\arcsin(\frac{|\vec s\sdot \vec f|}{|\vec s||\vec f|})\).

另外,由上面这个式子,不难得到一些特殊情况下的判定标准:

1. 若直线与平面平行,则\(am+bn+cp=0\).
2. 若直线与平面垂直,则\(\frac{a}{m}=\frac{b}{n}=\frac{c}{p}\).注意这里分母可能除以\(0\),我们实际上应该是三个形如\(a=mt\)的参数方程,这里简化了.

交点

联立方程硬解.

基本定理

参考:https://zhuanlan.zhihu.com/p/401766934

三余弦定理(最小角定理)

这个定理说明直线与平面的夹角,是所有包含直线的平面与这个平面形成的夹角中最小的那一个.并且偏移量决定了差距.

三正弦定理(最大角定理)

这个定理说明二面角是另一个平面上的直线与平面的夹角中最大的那个,并且偏移量决定了差距.