计算几何，我来辣！

前置的辅助工具

进入计算几何之前先来了解一下一些前置的定义 / 函数：

typedef double db;
typedef Vector Point;
const double eps = 1e-8;
const double PI = acos(-1.0);
int dcmp(db x)
{
    if(fabs(x) < eps) return 0;
    return x > 0 ? 1 : -1;    
}

这些是啥呢？

1. double 这个数据类型的简写。这样可以减少代码量。
2. 因为向量和点都是坐标存储，所以我们其实可以把它们看做同一个东西。
   • 我们先定义 Vector 这个类型后才能使用 typedef。
3. 精度控制。一般根据题目要求进行调整。
4. \(\pi\) 的预定义，方便之后可能的使用。
5. 在精度限制下进行比较的函数。只要精度够了就认为相等。

我们进入正题吧。

向量

在高中数学，我们就很熟悉向量了。不论是平面向量，空间向量，还是解析几何，向量都一直是一个解决问题的利器。当然，解析几何里，它也是一个基本的元素。

首先是存储和一些基本运算。这里使用重载运算符的方式实现。

struct Vector {
    db x, y;
    Vector(db x = 0, db y = 0) : x(x), y(y) {}
    Vector operator + (const Vector &A) const {
        return Vector(x + A.x, y + A.y);
    }
    Vector operator - (const Vector &A) const {
        return Vector(x - A.x, y - A.y);
    }
    Vector operator * (const db &A) const {
        return Vector(x * A, y * A);
    }
    Vector operator / (const db &A) const {
        return Vector(x / A, y / A);
    }
    bool operator == (const Vector &A) const {
        return dcmp(x - A.x) == 0 && dcmp(y - A.y) == 0;
    }
};

之后我们还要对向量的另一些基本信息和基本运算进行支持！

极角

一种方法是直接调用已有的函数。不过听说好像有毒瘤出题人对着卡精度？ 不过精度要求在 \(10^{-9}\) 以下的话这个方法应该是没有问题的，出题人硬要卡到 \(10^{-18}\) 就难说了。

db PolarAngle(Vector A) {
    return atan2(A.y, A.x);
}

旋转

众所周知，向量 \((x,y)\) 逆时针旋转 \(\alpha\) 度的话坐标会变成 \((x \cos \alpha - y \sin \alpha, y \cos \alpha + x \sin \alpha)\)。

直接搞就可以了。

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

点积和叉积

也是我们熟悉的内容。

db Dot(Vector A, Vector B) {
    return A.x * B.x + A.y * B.y;
}

db Cross(Vector A, Vector B) {
    return A.x * B.y - A.y * B.x;
}

不过我们更关心它们的一些性质：

1. 向量点积可以判断向量夹角。
   • \(\vec a \cdot \vec b = 0 \Leftrightarrow \vec a \perp \vec b\)
   • \(\vec a \cdot \vec b > 0 \Leftrightarrow \text{夹角为锐角}\)
   • \(\vec a \cdot \vec b < 0 \Leftrightarrow \text{夹角为钝角}\)
2. 向量叉积是有向面积。
   • \(\vec a \times \vec b\) 是 \(\vec a\) 和 \(\vec b\) 构成平行四边形的面积。
   • 有向性。当 \(\vec b\) 在 \(\vec a\) 左边的时候为正，在右边为负。
   • 当 \(\vec a \ \ /\kern -0.8em / \ \ \vec b\) 时为 \(0\)。

所以有一张很著名的图：

对于 \(\vec a \ \text{op} \ \vec b\) 来说：

这张图提供了感性理解：点积可以判断向量的“前后”关系，叉积可以判断向量的“左右”关系。

以此我们就可以进行一些实用操作：

计算三角形面积

利用叉积定义计算即可。

db Area(Point A, Point B, Point C) {
    return fabs(Cross(B - A, C - A) / 2);
}

计算向量模长

利用点积定义计算即可。

db Length(Vector A) {
    return sqrt(Dot(A, A));
}

计算向量夹角

利用点积定义计算即可。

db Angle(Vector A, Vector B) {
    return acos(Dot(A, B) / Length(A) / Length(B));
}

点在直线哪端

我们已知直线上两个点 \(A,B\)（或者，直线上一点 \(P\) 和方向向量 \(\vec v\)，本质一样），以及直线外一点 \(Q\)，可以利用叉积判断点 \(Q\) 在直线哪一端。

具体地，计算 \(\vec {PQ} \times \vec v\)，之后：

1. \(\vec {PQ} \times \vec v > 0\)，说明 \(Q\) 在直线下方。
2. \(\vec {PQ} \times \vec v = 0\)，说明 \(Q\) 在直线上。
3. \(\vec {PQ} \times \vec v < 0\)，说明 \(Q\) 在直线上方。

这个由叉积易知。

点到直线的距离

已知直线外一点 \(P\)，直线上两点 \(A,B\)。

做法大概是拿叉积得到平行四边形面积，再把底边长度除掉，就得到高，也就是点到直线距离。

db Dis_Point_to_Line(Point P, Point A, Point B)
{
    Vector u = B - A, v = P - A;
    return fabs(Cross(u, v)) / Length(u);
}

点到线段的距离

这个就需要讨论一下了。

考虑在点到直线距离的基础上，如果平行四边形高在外部，就不合理了，这个时候距离是到最近的线段端点的距离。

所以怎么做？

利用点积的性质判断一下夹角大小（前后）即可。

db Dis_Point_to_Seg(Point P, Point A, Point B)
{
    if(A == B) return Length(P - A); 
    // Case 1 : A 与 B 重合
    Vector v1 = B - A, v2 = P - A, v3 = P - B;
    if(dcmp(Dot(v1, v2)) < 0) return Length(v2);
    // Case 2 : 距离 A 近
    else if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
    // Case 3 : 距离 B 近
    return fabs(Cross(v1, v2)) / Length(v1);
    // Case 4 : 转化为点到直线距离
}

上图给出了前三种特殊情况的示意。

判断线段相交

首先是快速排斥实验：

如果两条线段“离得太远了”，根本没有相交的机会怎么办？我们把线段看成一个矩形，先看看这两个矩形是否相交？要是矩形都不相交，线段肯定没有相交的机会的！

这个就是“快速排斥实验”。但是，未通过快速排斥实验 是 两条线段不相交 的 充分不必要条件。也就是说，通过快速排斥实验 的线段 也可能不相交。

既然它只是一个充分不必要条件，所以我们一般应用的少。

更普遍的一个方法是 跨立实验。

之前我们知道了如何判断点在直线的哪一端；又由于两条线段相交，那么某线段的两个端点一定分别在另一线段端点两端；所以我们直接分别对于两条线段，各自判断一次它的两个端点和另一条线段所在直线的位置关系即可。

bool Is_Intersect(Point A, Point B, Point C, Point D)
{
    db c = Cross(B - A, C - A), d = Cross(B - A, D - A);
    db a = Cross(D - C, A - C), b = Cross(D - C, B - C);
    return dcmp(c) * dcmp(d) < 0 && dcmp(a) * dcmp(b) < 0;
} // 判断线段 AB 和 CD 是否相交

注意：

这里没有对于特殊情况进行讨论。如果题目不保证 / 有可能出现特殊情况，需要在上面的代码加入特判。

1. 如果认为 「两线段只有一个公共端点」 也算相交，那么需要特判。
2. 如果两线段会出现重合 / 部分重合 / 交点在其中一条线段上的情况，特判三点共线情况即可。

求两条直线的交点

设直线 \(l_1,l_2\)，线上各有一点 \(x_1,x_2\)，方向向量 \(\vec {v_1}, \vec {v_2}\)。

设交点为 \(p\)，那么 \(p\) 一定在 \(l_1\) 上，又在 \(l_2\) 上。这句话真的不是废话！因为我们接下来的公式推导全都是基于这个基本事实的。

首先，由 \(p\) 在 \(l_1\) 上，可以知道，\(p\) 可以表示成 \(x_1 + k \cdot \vec {v_1}\) 的形式，\(k\) 是一个实数。之后，由于 \(p\) 在 \(l_2\) 上，可以由叉积的性质知道：\((p - x_1) \times \vec {v_2} = 0\)。

把 \(p\) 代进去，由于叉积具有分配率，稍微化简一下：

\(\begin{aligned} (p - x_2) \times v_2 & = 0 \\ (x_1 + k \cdot v_1 - x_2) \times v_2 & = 0 \\ k \cdot v_1 \times v_2 & = (x_2 - x_1) \times v_2 \\ k & = \dfrac{(x_2 - x_1) \times v_2}{v_1 \times v_2} \end{aligned}\)

就能求出 \(k\) 了。之后再拿 \(p = x_1 + k \cdot \vec {v_1}\) 这个式子把 \(p\) 表示出来就得到了交点。

Vector Intersect(Point A, Vector u, Point B, Vector v)
{
    db k = Cross(B - A, v) / Cross(u, v);
    return A + u * k;
}

这里没有特判平行 / 重合之类的情况，所以可能出现 nan 或者 inf 之类的结果也不奇怪（笑）。

这种方法还有一种几何意义的解释，但是并不直观。

多边形

我们按顺时针或者逆时针存储所有顶点，再记录顶点个数即可。

Point a[N]; int n;

如果需要，可以封装为一个 Poly 的结构体。不过因为多边形的存储结构本身很简单所以一般是没有必要的。

判断点是否在多边形内部

1. “光线投射算法”
   • 先特判掉一些特殊情况。
     • 比如，类似快速排斥实验，可以判断出 「点离多边形太远了」 的情况。
     • 或者，可以直接比较点是否是多边形顶点。
     • 利用叉积也可以判断点是否在多边形某条边上。
   • 考虑这样的事实：我们考虑以该点为端点引出一条射线，数交点个数。
   • 如果这条射线与多边形有奇数个交点，则该点在多边形内部，否则该点在多边形外部，也就是 奇内偶外。

具体实现很麻烦，需要考虑一些特殊情况（并非指特判，而是引出射线和多边形交点判断部分的细节）：

int Is_inPoly(Point P, Point* a, int n)
{
    int cnt = 0; int dir, d1, d2;
    for(int i = 1; i <= n; i ++)
    {
        if(dcmp(Length(P - a[i])) == 0) return -1;
        // 点和多边形顶点重合
        dir = dcmp(Cross(a[i % n + 1] - a[i], P - a[i]));
        // 点和当前这个向量的位置关系
        if(dir == 0) return -1;
        // 点在多边形边上
        d1 = dcmp(a[i].y - P.y); 
        d2 = dcmp(a[i % n + 1].y - P.y); 
        if(dir > 0 && d1 <= 0 && d2 > 0) ++ cnt;
        if(dir < 0 && d2 <= 0 && d1 > 0) -- cnt;
    }
    if(cnt) return 1;
    else    return 0;
}

2. 回转数算法
   • 把这个点和多边形所有顶点连接起来，计算相邻两边夹角的和。
   • 若结果为 \(0\)，那么在多边形外。否则在多边形内。
   • 不是很常用。

计算多边形周长

这个很容易！直接计算即可。

db Calc_C_Poly(Point* a, int n)
{
    db C = 0;
    for(int i = 1; i <= n; i ++) 
        C += Length(a[i % n + 1] - a[i]);
    return C;
}

计算多边形面积

对于多边形进行一个三角剖分，之后拿叉积计算即可。

问题在于，我们要搞所谓“三角剖分”，还得选一个在多边形内的点？好像还得做很多麻烦事？

实际上不是。因为叉积可以算出负面积，所以平面内任选一个点进行三角剖分即可。这里选择原点，最简便。

也就是说，计算公式：

\(S = \dfrac{1}{2} \displaystyle \sum_{i = 1}^n (\overrightarrow{OA}_i \times \overrightarrow{OA}_{i \bmod n + 1})\)

db Calc_S_Poly(Point *a, int n)
{
    db S = 0;
    for(int i = 1; i <= n; i ++) 
        S += Cross(a[i], a[i % n + 1]);
    return S / 2;
}

圆

记录下圆心和半径即可。

struct Circle {
    Point O;
    db r;
    Circle(db a, db b, db c) : O(a, b), r(c) {}
}

求直线与圆的交点

这个就和咱解析几何干的事差不多（

判断圆心到直线距离。

• 大于半径，没得交点。
• 等于半径，有一个交点。找一个过圆心的与已知直线垂直的直线，大力求交点即可。
• 小于半径，有两个交点。同样，找一个过圆心的与已知直线垂直的直线，求个交点；之后我们算一下交点到圆心的距离，再拿勾股定理求个半弦长；交点沿着已知直线的方向向量的正反方向分别挪动半弦长即可。

这个实现相对容易。

圆和圆的位置关系

直接比较圆心距和半径关系即可。