跳转至

计算几何基础

前置知识：高中的基础解析几何知识。

比如平面/空间向量，极坐标系，直线与圆一类的。

为了计算几何花了点时间把之前的文化课内容捡起来了，有点感叹。

因为不知道有什么教程比较好，所以结构上基本是参考 OI-wiki 然后做了改进。

内容的话，OI-wiki 也有些地方写的不好，所以大部分还是自己总结的。

本文的代码目前还未经过验证，仅供参考，未验证代码会有一个 FIX 标记

约定ψ(｀∇´)ψ

能不写浮点数就不写浮点数。

不需要拘泥于模板化，能用就行，学的时候模板化一点方便以后弄

常用函数/操作ψ(｀∇´)ψ

首先浮点数比较大小这个显然是需要的：

设置一个合适的 eps 即可，通常来说 binary64 标准浮点数已经够用了。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
const double eps = 1e-8;
int sgn(double x) {
    return (x > eps) - (x < -eps);
}
bool cmp(double x, double y) {
    return sgn(x - y) < 0;
}
bool equal(double x, double y) {
    return sgn(x - y) == 0;
}

对于 \(\pi\)，可以直接背下它的后几位小数值。

更为通用的方式是：const double pi = acos(-1.0);，这是能用库函数得到的最接近 \(\pi\) 的浮点数。

然后是 double atan2(double y, double x) 这个函数，它的作用是返回 \(\arg x+yi\)，弧度制，或者说 \(\arctan(\dfrac{y}{x})\)。

相比于 atan()，它的优势在于值域，atan2() 的值域是 \([-\pi, \pi]\)；

众所周知 \(\tan\) 的值在一三，二四象限分别是相等的，因为 atan() 不知道这个值到底是在哪个象限，所以标准会规定为这个角一定不在 \(y\) 轴的左边，所以值域为 \([-\dfrac{\pi}{2},\dfrac{\pi}{2}]\)。

而 atan2 则可以分辨出这个角到底是在哪个象限，并且对边界有处理（因为你给了两个参数）。

acos(), asin() 也是一样的，他们的值域长度也都只有 \(\pi\) ，acos() 只取了 \([0, \pi]\)，asin() 只取了 \([-\dfrac{\pi}{2},\dfrac{\pi}{2}]\)。

数据记录ψ(｀∇´)ψ

对于平面/空间中的一个点，我们直接记录其坐标，如果是极坐标系就记录极角和极径。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
struct Point { // 2D only.
    double x, y;
    Point(double ix = 0, double iy = 0) : x(ix), y(iy) {}
    Point operator + (const Point &rhs) const {
        return Point(x + rhs.x, y + rhs.y);
    } // 向量加要用
    Point operator - (const Point &rhs) const {
        return Point(x - rhs.x, y - rhs.y);
    } // 初始化向量要用
    Point operator * (const double &rate) const {
        return Point(x * rate, y * rate);
    } // 向量数乘要用
} orig;

对于平面/空间中的一个向量，使用向量的坐标表示，记录其坐标表示的坐标即可。

有时候需要旋转一个向量，这个东西直接用和差角推一下就行。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
struct Vector : public Point {
    double len() { return sqrt(x * x + y * y); }
    Vector(Point st = orig, Point ed = orig) {
        Point vec = ed - st;
        x = vec.x, y = vec.y;
    }
    Vector &operator = (const Point &vec) {
        x = vec.x, y = vec.y;
        return *this;
    } // 后面有构造函数会用到
    Vector rot(double t) {
        double nx = x * cos(t) - y * sin(t);
        double ny = x * sin(t) + y * cos(t);
        x = nx, y = ny;
        return *this;
    } // 逆时针旋转
    Vector unit() {
        double l = len();
        x /= l, y /= l;
        return *this;
    } // 返回同方向的单位向量
    Vector vert() {
        swap(x, y);
        if(sgn(x) != 0) 
            x = x * -1.0;
        else y = y * -1.0;
        return *this;
    } // 返回一个垂直于此向量的非零向量
      // 若返回值为零向量则说明该向量也为零向量
    Point toPoint() { return Point(x, y); }
} zero; 
// 个人习惯是构造函数用有向线段，等号用坐标。
// 因为会有 Vector = xx 的形式，而且构造函数也可以直接用 (orig, pt) 来构造。

double dot(Vector a, Vector b) {
    return (a.x * b.x) + (a.y * b.y);
}
double cross(Vector a, Vector b) {
    return (a.x * b.y) - (b.x * a.y);
}

对于平面/空间中的一条线段，直接记录两个端点就行。

1
2
3
4
5
struct Segment {
    Point st, ed;
    Segment(Point a = orig, Point b = orig) : st(a), ed(b) {}
    Line toLine() { return Line(st, Vector(st, ed)); }
};

对于平面/空间中的一条直线，我们没有必要记录解析式，有向量就够了，我们只需要记录直线上一点以及它的一个方向向量。

1
2
3
4
5
6
struct Line {
    Point pt;
    Vector vec;
    Line(Point p = orig, Vector v = orig) : pt(p), vec(v) {}
    void adjust() { if(vec.y < 0) vec = zero - vec; } // 调整方向
};

对于平面/空间中的一条射线，我们类似直线的方式来记录，不过记录的是起始点和与射线方向相同的方向向量。

对于平面中的一个多边形，我们直接记录多边形的所有顶点即可，一般按照逆时针顺序记录。

1
using Polygon = std::vector<Point>;

对于平面中的一类曲线，直接记录解析式，圆可以直接记录圆心和半径。

1
2
3
4
5
struct Circle {
    Point o;
    double r;
    Circle(Point io = orig, double ir = 0) : o(io), r(ir) {}
};

一些公式ψ(｀∇´)ψ

正余弦定理ψ(｀∇´)ψ

正弦定理ψ(｀∇´)ψ

对于 \(\Delta ABC\)，存在 \(\dfrac{a}{\sin A} = \dfrac{b}{\sin B} = \dfrac{c}{\sin C} = 2R\)。

其中 \(R\) 为 \(\Delta ABC\) 的外接圆的半径。

余弦定理ψ(｀∇´)ψ

对于 \(\Delta ABC\)，存在 \(a^2 = b^2 + c^2 + 2bc\cos A\)。

其它的两个公式通过轮换可以得到。

其它公式ψ(｀∇´)ψ

对于到角公式和夹角公式，因为我们一般使用向量去记录直线，所以没有必要记了，直接用向量求就行。

直线的到角定义是 \(l1 \to l2\) 逆时针转过的最小角，夹角定义为两个到角中的最小角。

而向量的夹角定义不太一致，它的取值是 \([0, \pi]\)，因为向量是具有方向的（这个夹角需要是两个向量形成的射线所构成的角）。

  • 向量夹角：直接用数量积就可以求，刚好 acos() 值域也对应的上。
  • 直线到角：考虑把两个方向向量调整到 \(y > 0\)，然后求一下叉积看一下方向，如果叉积为正，说明两向量夹角就是到角，否则如果为负，则到角为两向量夹角的补角。
  • 直线夹角：直接利用定义求出两个到角返回小的一个即可。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
double angle(Vector a, Vector b) {
    return acos(dot(a, b) / (a.len() * b.len()));
}
double rad(Line a, Line b) {
    double cs = cross(a.vec, b.vec);
    double theta = angle(a.vec, b.vec);
    if(sgn(cs) == 0) {
        if(sgn(a.vec.y) == sgn(b.vec.y))
            return 0.0;
        else return pi;
    }
    if(sgn(cs) > 0) return theta;
    if(sgn(cs) < 0) return pi - theta;
}
double angle(Line a, Line b) {
    double ab = rad(a, b);
    double ba = rad(b, a);
    if(cmp(ab, ba)) 
        return ab;
    else return ba;
}
// FIX

对于点到直线的距离公式，也没有必要记，我们可以通过向量来构造，这个等下会说。

三角恒等变换ψ(｀∇´)ψ

其实就是和差角，和差化积，万能公式什么的，高中应该都会涉及，我不想写了。

基本操作ψ(｀∇´)ψ

点点关系ψ(｀∇´)ψ

欧几里得距离ψ(｀∇´)ψ

1
2
3
double edis(Point a, Point b) {
    return Vector(b, a).len();
}

曼哈顿距离ψ(｀∇´)ψ

1
2
3
double mdis(Point a, Point b) {
    return fabs(a.x - b.x) + fabs(a.y - b.y);
}

切比雪夫距离ψ(｀∇´)ψ

1
2
3
double cdis(Point a, Point b) {
    return max(fabs(a.x - b.x), fabs(a.y - b.y));
}

转化ψ(｀∇´)ψ

以上的三个距离是可以相互转化的，不过暂时先咕咕掉。

在 Gym100739E-Life as a Monster 这个题里会用到，有比较详细的推导。

点线关系ψ(｀∇´)ψ

点到直线的距离ψ(｀∇´)ψ

一定要注意的一个事情是我们是在使用坐标来表示向量，相当于是，把所有向量都变成 \((0, 0)\) 起始再做运算。

假设点为 \(Q\), 直线为 \(l:\{P, \textbf{n}\}\)。

我们考虑构造出这个距离，尝试将其作为某个图形的高，想到平面向量的叉积的几何意义就是一个平行四边形的面积。

于是构造出一个底边在 \(l\) 上的平行四边形，这个底边的邻边为 \(PQ\)，于是可以知道这个平行四边形的面积为 \(\vec{PQ}\times \textbf{n}\)。

然后我们要求的就是这个高，对应底边长就是 \(|\textbf{n}|\)，于是就求出来了。

1
2
3
4
double dis(Point p, Line l) {
    double s = fabs(cross(Vector(l.pt, p), l.vec));
    return s / l.vec.len();
}

点到线段的距离ψ(｀∇´)ψ

\(P \to AB\)。

需要考虑有没有垂足存在，可以考察构成的三角形中，两个底角是否存在钝角。

然后转化成点到点的距离和点到线段的距离就行。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
double dis(Point p, Segment l) {
    double ret = 0;
    Point A = l.st, B = l.ed;
    if(cmp(pi / 2.0, angle(Vector(A, p), Vector(A, B))))
        return edis(p, A);
    if(cmp(pi / 2.0, angle(Vector(B, p), Vector(B, A))))
        return edis(p, B);
    return dis(p, l.toLine());
}
// FIX

点在直线的哪一侧ψ(｀∇´)ψ

还是可以用叉积的性质，构造出 \(\vec{PQ}\)，然后计算 \(s = \vec{PQ}\times \textbf{n}\)。

如果 \(s > 0\) 就是逆时针转，如果 \(s < 0\) 就是顺时针转，否则就是在线上。

不过需要注意方向向量的方向，这个会对结果有影响，一般来说我们希望方向向量 \(y \ge 0\)，在实现时请注意，这对应着上面实现中的 void adjust()。

1
2
3
int side(Point p, Line l) {
    return sgn(cross(Vector(l.pt, p), l.vec));
}

线线关系ψ(｀∇´)ψ

求夹角ψ(｀∇´)ψ

在上面已经说过了。

判断平行/垂直ψ(｀∇´)ψ

（直线）

这个直接使用向量就行了！

1
2
3
4
5
6
7
8
bool vertical(Vector a, Vector b) {
    return sgn(dot(a, b)) == 0;
}
bool parallel(Vector a, Vector b) {
    return equal(a.x / b.x, a.y / b.y);
}

// vertical/parallel(a.vec, b.vec);

求到角也是直接用向量就行了，如果你要夹角（\(\in (0, \pi)\)）的话判断一下就好。

因为是平面内的直线所以相交没有什么好判的，只要不是平行或者垂直就是相交。

（线段）

平行和垂直直接把线段变成向量，然后比较就行。

然后至于相交的话，因为是有端点存在的，发现相交实际上就是两个端点被另一条线段分开在两边。

这个可以通过上面点和直线的关系进行判断，然后你发现对于一条线段满足还不是充分条件。

可以发现应该是对于两条线段都满足这个限制，但是需要特判线段是一个点的情况。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
bool intersect(Segment a, Segment b) {
    Vector va = Vector(a.st, a.ed);
    Vector vb = Vector(b.st, b.ed);
    Line la = a.toLine(), lb = b.toLine();

    if(va == zero) return side(va, lb) == 0;
    if(vb == zero) return side(vb, la) == 0;

    if(parallel(va, vb) || vertical(va, vb)) 
        return false;
    return (side(a.st, lb) != side(a.ed, lb)) 
        && (side(b.st, la) != side(b.ed, la));
}
// FIX

求交点ψ(｀∇´)ψ

（直线）

网上看到的大部分的博客都直接假定成线段来做。

但实际上只知道方向向量和直线上一点并不一定能平移出两个相交的线段。

考虑求出某一个点到交点的距离。

1

假定线上的两个点为 \(A,B\)，我们考虑找出 \(AC\) 这个距离。

考虑以 \(\textbf{b}\) 为底构造两个三角形，可以求出高构造出一个相似。

如图所示，\(s1,s2\) 都是可以直接叉积求出来的，然后 \(|TU|,|AS|\) 的比值就是 \(|\textbf{a}|\) 和 \(|AC|\) 的比值。

于是我们平移就可以了，因为叉积自带方向所以不需要额外关心了。

1
2
3
4
5
Point intersecion(Line a, Line b) {
    Vector t = Vector(a.pt, b.pt);
    double rate = (cross(b.vec, t)) / (cross(b.vec, a.vec));
    return a.pt + a.vec * rate;
}

（线段）

这个比直线好理解一点，但是代码可能略麻烦，建议直接转成直线求。

我们有向量，有点，可以直接算很多东西，比如面积一类的。

2

仍旧是考虑通过某个已知的点平移求出交点。

考虑构造定比分点的形式，这样能搞出某个点到交点的向量。

显然存在 \(k\in \mathbb{R}\text{ s.t. } k\vec{CA} + (1 - k)\vec{CB} = \vec{CE}\)。

考虑怎么求出这个 \(k\)，实际上就是求 \(T = \dfrac{|AE|}{|EB|}\)。

可以有 \(T = \dfrac{S_{\Delta AED}}{S_{\Delta DEB}} = \dfrac{S_{\Delta AEC}}{S_{\Delta{CEB}}}\)。

两个比例显然可以直接合并，于是 \(T = \dfrac{S_{\Delta ACD}}{S_{\Delta BCD}}\)。

显然这两个用向量叉积可计算（平行四边形面积一半，记得加绝对值），然后比值就可以计算了。

\(T = \dfrac{|\vec{AC}\times\vec{AD}|}{|\vec{BC}\times\vec{BD}|}\)。

然后 \(T = \dfrac{k}{1-k}\Rightarrow k = \dfrac{T}{T+1}\)。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
Point intersection(Segment a, Segment b) {
    return intersection(a.toLine(), b.toLine());
}
Point intersection(Segment a, Segment b) {
    assert(intersect(a, b) == true);
    Point A = a.st, B = a.ed;
    Point C = b.st, D = b.ed;
    double T = fabs(cross(Vector(A, C), Vector(A, D)))
            / fabs(cross(Vector(B, C), Vector(B, D)));
    double k = T / (T + 1.0);
    Vector vec = Vector(C, A) * k + Vector(C, B) * (1.0 - k);
    return C + vec.toPoint();
}
// FIX

圆相关ψ(｀∇´)ψ

线与圆的位置关系ψ(｀∇´)ψ

考虑 \(r\) 和 \(dis(o, l)\) 即可。

1
2
3
4
5
6
7
int side(Line l, Circle C) {
    double d = dis(C.o, l);
    if(cmp(r, d)) return 0; // 相离
    if(equal(r, d)) return 1; // 相切
    if(cmp(d, r)) return 2; // 相交
}
// FIX

线与圆的交点ψ(｀∇´)ψ

先判断位置关系。

如果相交，有圆心和线的距离，考虑勾股定理求出弦长。

然后可以构造出圆心和垂足所在的直线，平移即可得到两个交点（代码中切点会被返回两次）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
double sqr(double val) { return val * val; }
std::pair<Point, Point> intersections(Line L, Circle C) {
    if(side(L, C) == 0) return make_pair(orig, orig);
    double l = (sqr(C.r) - sqr(dis(C.r, L))) / 2.0;
    Point mid = intersection(L, Line(C.o, L.vec.vert()));
    Point A = mid + (L.vec.unit() * l).toPoint();
    Point B = mid * 2.0 - A;
    return make_pair(A, B);
}
// FIX;

圆与圆的位置关系ψ(｀∇´)ψ

考虑两个 \(r\) 以及 \(o\) 之间的距离的关系即可。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
int side(Circle C1, Circle C2) {
    int odis = edis(C1.o, C2.o);
    if(cmp(C1.r + C2.r, odis)) return 0; // 相离
    if(equal(C1.r + C2.r, odis)) return 1; // 外切
    if(cmp(odis, C1.r + C2.r)) return 2; // 相交
    if(equal(C1.r, C2.r + odis) || equal(C2.r, C1.r + odis))
        return 3; // 内切
    return 4 // 内含
}
// FIX

圆与圆的交点ψ(｀∇´)ψ

如果圆是内切或者外切直接转化成求线和圆的交点。

如果是相交的话，考虑连接某个交点和两个圆心。

可以通过余弦定理求出该交点与某个圆心连线和圆心连线的夹角。

于是把圆心连线的方向向量旋转到对应位置，长度变为半径，然后把圆心平移过去就行了。

可以发现其实对相切的情况同样适用。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
std::pair<Point, Point> intersections(Circle C1, Cirle C2) {
    int side_c = side(C1, C2);
    if(side_c == 0 || side_c == 4) 
        return make_pair(orig, orig);
    double d = Vector(C1.o, C2.o).len();
    double t = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d));
    Point A = C1.o + (Vector(C1.o, C2.o).unit().rot(t) * C1.r).toPoint();
    Point B = C1.o + (Vector(C1.o, C2.o).unit().rot(t * -1.0) * C1.r).toPoint() 
    return make_pair(A, B);
}
// FIX

多边形相关ψ(｀∇´)ψ

周长和面积ψ(｀∇´)ψ

周长直接对于每一条边求和就行。

面积的话用叉积的几何意义就行：\(S = \dfrac{1}{2}|\sum\limits_{i = 0}^{n - 1}\textbf{v}_i\times\textbf{v}_{(i + 1) \bmod n}|\)。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
double circum(Polygon p) {
    double ret = 0.0;
    int sz = (int)p.size();
    for(int i = 0; i < sz; ++i) {
        ret += Vector(p[i], p[(i + 1) % sz]).len();
    }
    return ret;
}
double area(Polygon p) {
    double ret = 0.0;
    int sz = (int)p.size();
    auto v = [&](int idx) -> Vector {
        return Vector(p[idx], p[(idx + 1) % sz]);   
    };
    for(int i = 0; i < sz; ++i) {
        ret += cross(v(i), v((i + 1) % sz));
    }
    return ret / 2.0;
}
// FIX

Point in Polygonψ(｀∇´)ψ

就是判一个点是不是在多边形内部。

有两种常见做法，这里用的是更简单的一种回转数算法。

回转数被定义为闭合曲线绕过某个点的次数，若回转数为 \(0\) 则说明点在闭合曲线之外。

考虑把这个点和多边形的所有点连起来，计算这些边相邻两个边的夹角之和，这个夹角是有方向的，符号和叉积符号一致。

如果和为 \(0\)，那么点在多边形外，否则在多边形内，回转数实际上就是代码中的 ret 除以 \(2\pi\)。

另外敬请注意，这里的 Polygon 并不要求以什么顺序存储，但为了规范，我们一般规定以逆时针顺序储存。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
bool inside(Point p, Polygon g) {
    double ret = 0.0;
    int sz = (int)g.size();
    auto v = [&](int idx) -> Vector {
        return Vector(p, g[idx]);
    };
    for(int i = 0; i < sz; ++i) {
        ret += angle(v(i), v((i + 1) % sz))
                * sgn(cross(v(i), v((i + 1) % sz)));
    }
    return sgn(ret) != 0;
}
// FIX

当然如果要多判在边上或者是点上也是很容易的，这里略去。

极角排序ψ(｀∇´)ψ

对于给定的一个极坐标系，对所有点按极角排序。

通常来说都是从直角坐标系转换而来。

一种比较暴力的方式就是直接利用 atan2 求出极角，然后排序。

1
2
3
4
5
6
7
8
double theta(auto v) { return atan2(v.y, v.x); } 
void polarSort(std::vector<Point> &p, Point c = origin) {
    sort(p.begin(), p.end(), [&](int x, int y) -> {
        return cmp(theta(p[x] - c), theta(p[y] - c))
            && cmp(Vector(c, p[x]).len(), Vector(c, p[y]).len());
    });
}
// FIX

但这样有时候会被恶心出题人卡精度，很烦。

另外注意到叉积可以判断方向，而极角排序实际上就是从极轴出发逆时针转。

所以我们也可以直接看叉积的符号确定它们的关系。

1
2
3
4
5
6
7
void polarSort(std::vector<Point> &p, Point c = orig) {
    sort(p.begin(), p.end(), [&](int x, int y) -> {
        return sgn(cross(Vector(c, p[x]), Vector(c, p[y]))) > 0
            && cmp(Vector(c, p[x]).len(), Vector(c, p[y]).len());
    });
}
// FIX

后者常数略大。

最后更新: November 8, 2023