计算几何深度实战：从凸包、最近点对与扫描线到 Voronoi 图、Delaunay 三角化与几何查询的完整工程链路

计算几何（Computational Geometry）是算法工程里最容易被低估、却在底层无处不在的一门学问。碰撞检测、GIS 空间索引、计算机图形学光栅化、机器人运动规划、地图瓦片切分、芯片布局布线，背后都是几何原语的组合。本文沿"几何基础 → 经典算法 → 空间剖分 → 工程映射"的主线，用可运行的 Python 把凸包、最近点对、扫描线、点在多边形内、旋转卡壳、Delaunay/Voronoi 逐一打通，并落到工程选型与 AI/系统场景。

一、几何基础：叉积与方向判定

一切平面几何算法的基石，是把"点在线的哪一侧""两条线段是否左转"这类直觉，翻译成有符号面积（叉积）。对向量 OA=(ax-ox, ay-oy) 与 OB=(bx-ox, by-oy)，二维叉积取其 z 分量：

cross(O,A,B) = (Ax-Ox)*(By-Oy) - (Ay-Oy)*(Bx-Ox)

• cross > 0：B 在 OA 的左侧（逆时针 / 左转）
• cross < 0：B 在 OA 的右侧（顺时针 / 右转）
• cross = 0：三点共线

这个原语既是凸包"是否该弹栈"的判据，也是 Delaunay 外接圆测试的组成部分。先把它跑通：

def cross(o, a, b):
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def ccw(a, b, c):
    return cross(a, b, c) > 0  # 严格左转

# 验证：以 O(0,0) 为支点，A(1,0) 到 B(0,1) 是左转
print("cross((0,0),(1,0),(0,1)) =", cross((0, 0), (1, 0), (0, 1)))   # 1.0 -> 左转
print("ccw left-turn:", ccw((0, 0), (1, 0), (0, 1)))                  # True
print("cross((0,0),(1,0),(1,-1)) =", cross((0, 0), (1, 0), (1, -1)))  # -1.0 -> 右转

二、凸包：Andrew 单调链

凸包（Convex Hull）是包含所有点的最小凸多边形，等价于"用一根橡皮筋套住所有点"。Andrew 单调链（Monotone Chain）先按 x（再按 y）排序，分别构造下链与上链，维护一个栈，当栈顶两段不是左转（即 cross <= 0）就弹栈——本质就是不断"向内收口"。时间复杂度 O(n log n)（瓶颈在排序）。

def cross(o, a, b):
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def convex_hull(points):
    pts = sorted(set(points))
    if len(pts) <= 1:
        return pts
    def build(hull_pts):
        res = []
        for p in hull_pts:
            while len(res) >= 2 and cross(res[-2], res[-1], p) <= 0:
                res.pop()
            res.append(p)
        return res
    lower = build(pts)
    upper = build(pts[::-1])
    return lower[:-1] + upper[:-1]  # 去掉重复端点

pts = [(0, 0), (1, 1), (2, 0), (1, 2), (0.5, 0.5), (2, 2), (0, 2), (3, 1)]
hull = convex_hull(pts)
print("hull (CCW):", hull)
print("hull size:", len(hull), "of", len(set(pts)), "unique points")

注意 cross <= 0 而非 < 0：共线点会被剔除，保证边界最简。若想保留共线边界点，改成 < 0 即可。

三、最近点对：分治

给定 n 个点求最小距离，暴力是 O(n^2)。分治法先按 x 排序，递归左右两半得到最小距离 d，再只在跨中缝、且 x 距离 < d 的"窄带"里检查——因为 y 方向上超过 d 的点不可能构成更近对，内层循环可被 break 截断，整体回到 O(n log n)。

import math, random

def dist(a, b):
    return math.hypot(a[0] - b[0], a[1] - b[1])

def closest_pair(px):
    def cp(pts):
        n = len(pts)
        if n <= 3:
            d = float('inf'); best = None
            for i in range(n):
                for j in range(i + 1, n):
                    if dist(pts[i], pts[j]) < d:
                        d = dist(pts[i], pts[j]); best = (pts[i], pts[j])
            return d, best
        m = n // 2
        mid = pts[m]
        d1, p1 = cp(pts[:m])
        d2, p2 = cp(pts[m:])
        d, best = (d1, p1) if d1 < d2 else (d2, p2)
        strip = [p for p in pts if abs(p[0] - mid[0]) < d]
        strip.sort(key=lambda q: q[1])
        for i in range(len(strip)):
            for j in range(i + 1, len(strip)):
                if strip[j][1] - strip[i][1] >= d:
                    break
                if dist(strip[i], strip[j]) < d:
                    d = dist(strip[i], strip[j]); best = (strip[i], strip[j])
        return d, best
    return cp(sorted(px, key=lambda p: (p[0], p[1])))

random.seed(7)
data = [(random.random() * 10, random.random() * 10) for _ in range(40)]
d, pair = closest_pair(data)
print(f"min dist = {d:.4f}, pair = {pair}")

窄带内按 y 排序后，每个点只需与后面 y 距离 < d 的点比较，这是算法压到 O(n log n) 的关键。

四、扫描线：矩形面积并

"扫描线"（Sweep Line）用一个假想线扫过平面，把二维问题降成一维事件序列。求多个轴对齐矩形的并集面积：沿 x 排序所有竖边，逐段处理，对当前 x 区间里"被覆盖的 y 区间"做合并计数，再乘以 dx 累加。

def rectangle_union_area(rects):
    # rects: list of (x1, y1, x2, y2)，约定 x1<x2, y1<y2
    xs = sorted(set([r[0] for r in rects] + [r[2] for r in rects]))
    area = 0.0
    for i in range(len(xs) - 1):
        xL, xR = xs[i], xs[i + 1]
        active = [(min(r[1], r[3]), max(r[1], r[3]))
                  for r in rects if r[0] <= xL and r[2] >= xR]
        if not active:
            continue
        active.sort()
        ycov = 0.0; cur_s, cur_e = active[0]
        for s, e in active[1:]:
            if s <= cur_e:
                cur_e = max(cur_e, e)      # 区间合并
            else:
                ycov += cur_e - cur_s
                cur_s, cur_e = s, e
        ycov += cur_e - cur_s
        area += (xR - xL) * ycov
    return area

rects = [(0, 0, 4, 4), (2, 2, 6, 6), (5, 0, 8, 3)]
print("union area:", rectangle_union_area(rects))  # 36.0

A=16、B=16、C=9，两两重叠 4 与 1，并集 = 16+16-4 + 9-1 = 36。扫描线把"二维覆盖"变成"一维区间合并"，同一思路还能解线段求交、天际线等问题。

五、点在多边形内：射线法

判断点是否在多边形内部，最稳健的通用法是射线法（Ray Casting）：从点向右引一条水平射线，统计它与多边形边相交的次数——奇数次在内部，偶数次在外部。处理水平边与顶点穿越时，用"yi>y != yj>y"的开闭约定避免重复计数。

def point_in_polygon(p, poly):
    x, y = p
    n = len(poly)
    inside = False
    j = n - 1
    for i in range(n):
        xi, yi = poly[i]
        xj, yj = poly[j]
        if ((yi > y) != (yj > y)) and (x < (xj - xi) * (y - yi) / (yj - yi) + xi):
            inside = not inside
        j = i
    return inside

poly = [(0, 0), (4, 0), (4, 4), (0, 4), (0, 0)]  # 正方形
print("center (2,2):", point_in_polygon((2, 2), poly))   # True
print("outside (5,5):", point_in_polygon((5, 5), poly))  # False
print("near-right-edge (3.9,2):", point_in_polygon((3.9, 2), poly)) # True

若多边形是凸的，有更快的"所有边同侧"判定；射线法胜在对任意简单多边形都正确，是点选、命中测试、地理围栏的标准实现。

六、旋转卡壳：凸包直径

得到凸包后，最远点对（凸包直径）可由旋转卡壳（Rotating Calipers）在 O(n) 内求出：把一对"对踵点"沿凸包逆时针同步推进，总步数 O(n)。它还能顺手解出最小外接矩形、最远点对、凸多边形宽度等。

def dist(a, b):
    return math.hypot(a[0] - b[0], a[1] - b[1])

def rotating_calipers_diameter(hull):
    n = len(hull)
    if n == 1:
        return 0.0
    if n == 2:
        return dist(hull[0], hull[1])
    best = 0.0
    j = 1
    for i in range(n):
        while True:
            nj = (j + 1) % n
            if dist(hull[i], hull[nj]) > dist(hull[i], hull[j]):
                j = nj
            else:
                break
        best = max(best, dist(hull[i], hull[j]))
    return best

pts = [(0, 0), (1, 1), (2, 0), (1, 2), (2, 2), (0, 2), (3, 1)]
hull = convex_hull(pts)
print("hull:", hull)
print("diameter:", rotating_calipers_diameter(hull))

七、Delaunay 三角化与 Voronoi 图

Delaunay 三角化把点集连成三角形，约束是"任意三角形的外接圆内不含其他点"（空圆性质）。它最大化最小角、避免狭长三角形，是网格生成、插值、地形建模的默认选择。其对偶就是 Voronoi 图：每个站点（site）的 Voronoi 胞腔，是所有"离它最近的点"构成的区域，是最近邻查询、细胞分割、设施选址的几何基础。

下面用纯 Python 实现 Bowyer-Watson 增量算法构造 Delaunay（不依赖外部库），再给出基于 scipy 的 Voronoi 胞腔提取作为对照。

import math

def circumcircle(a, b, c):
    ax, ay = a; bx, by = b; cx, cy = c
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if abs(d) < 1e-12:
        return None
    ux = ((ax * ax + ay * ay) * (by - cy) + (bx * bx + by * by) * (cy - ay) +
          (cx * cx + cy * cy) * (ay - by)) / d
    uy = ((ax * ax + ay * ay) * (cx - bx) + (bx * bx + by * by) * (ax - cx) +
          (cx * cx + cy * cy) * (bx - ax)) / d
    r = math.hypot(ax - ux, ay - uy)
    return (ux, uy, r)

def in_circumcircle(p, tri):
    cc = circumcircle(*tri)
    if cc is None:
        return False
    return math.hypot(p[0] - cc[0], p[1] - cc[1]) < cc[2] - 1e-9

def bowyer_watson(points, margin=10.0):
    if not points:
        return []
    minx = min(p[0] for p in points); maxx = max(p[0] for p in points)
    miny = min(p[1] for p in points); maxy = max(p[1] for p in points)
    dmax = max(maxx - minx, maxy - miny) * margin or margin
    midx = (minx + maxx) / 2; midy = (miny + maxy) / 2
    st = [(midx - 2 * dmax, midy - dmax), (midx + dmax, midy - dmax), (midx, midy + 2 * dmax)]
    triangles = [tuple(st)]
    for p in points:
        bad = [t for t in triangles if in_circumcircle(p, t)]
        edges = []
        for t in bad:
            edges += [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]
        boundary = [e for i, e in enumerate(edges)
                    if not any(i != j and set(e) == set(edges[j]) for j in range(len(edges)))]
        triangles = [t for t in triangles if t not in bad]
        for ea, eb in boundary:
            triangles.append((ea, eb, p))
    stset = set(st)
    return [t for t in triangles if not (set(t) & stset)]

random.seed(3)
pts = [(random.random() * 5, random.random() * 5) for _ in range(15)]
tris = bowyer_watson(pts)
print("Delaunay triangles:", len(tris))
# 每个 Delaunay 三角形的外接圆心，正是 Voronoi 图的顶点
vverts = [circumcircle(*t)[:2] for t in tris if circumcircle(*t)]
print("Voronoi vertices (circumcenters):", len(vverts))

scipy.spatial 直接给出 Voronoi 胞腔，验证上面"对偶"关系的正确性：

import numpy as np
from scipy.spatial import Delaunay, Voronoi, voronoi_plot_2d

random.seed(3)
pts = np.array([(random.random() * 5, random.random() * 5) for _ in range(15)], float)
tri = Delaunay(pts)
v = Voronoi(pts)
print("Delaunay simplices:", tri.simplices.shape)
print("Voronoi points (vertices):", v.vertices.shape, "regions:", len(v.regions))

Delaunay 三角形数 ≈ 2n，Voronoi 顶点数 ≈ 2n-5——两者共享同一套"最近邻"拓扑结构，只是边被翻转成了对偶。

八、工程映射：计算几何在 AI 与系统里的位置

几个值得记住的"映射"：

• 向量数据库的近邻查询 ≈ 最近点对问题的工程化，只是维度从 2D 升到高维，必须用 KD-tree/球树/量化索引替代暴力。
• 碰撞检测 是凸包 + 旋转卡壳 + 包围盒（BVH/AABB）的组合，物理引擎每帧都在跑。
• GPU 光栅化 本质是扫描线思想在硬件上的并行实现；Z-buffer 就是一条逐像素的"深度扫描线"。
• 地图瓦片与地理围栏 依赖点在多边形内 + 空间剖分，是 LBS 服务的基本面。

九、工程选型

十、边界与陷阱

• 精度：叉积/外接圆用浮点会触发共线、共圆退化；生产环境应引进"鲁棒几何谓词"（如 Shewchuk 的精确算术），否则可能漏判相交。
• 退化输入：多点共线、三点共圆会让"左转/空圆"判据边界抖动，需约定 <= 与 < 的取舍。
• 维度灾难：平面上优雅的 Delaunay/Voronoi，到了高维后三角形数爆炸、数值不稳定，高维近邻应改用近似索引（HNSW/IVF）。
• 坐标尺度：super_triangle 与尺度相关；点集跨数量级时先归一化再计算。

结语

计算几何把"形状"变成了可计算、可验证、可并行的对象。掌握叉积判定、凸包、分治最近点对、扫描线与 Delaunay/Voronoi 这条主线，你就拥有了从碰撞检测到向量检索、从地图围栏到网格生成的一把通用钥匙。它不像深度学习那样吸睛，却是工程底座里最恒久、最可靠的那一层。