8.4 SymPy 與 NumPy 的計算幾何實作#

學習動機：讀完 8.1–8.3，你已經知道凸包、Voronoi、Delaunay 的原理。但真要動手，你會撞上兩個現實問題。第一，SymPy 有沒有現成的凸包函式？如果有，它用了哪個演算法？跟你自己寫的 Graham scan 差在哪？第二，SymPy 的精確算術很慢——當你手上是 NumPy 的十萬個浮點座標時，正確做法是「用 NumPy 打架、用 SymPy 驗屍」。本節就是把兩者接起來：SymPy 給你正確性，NumPy 給你速度，matplotlib 給你眼睛。

8.4.1 SymPy 有 `ConvexHull` 嗎？#

有，但名字不叫 ConvexHull。它在 sympy.geometry.util，叫 convex_hull：

from sympy.geometry import Point
from sympy.geometry.util import convex_hull

基本用法#

# 08.4_sympy_convex_hull.py
from sympy.geometry import Point, Segment
from sympy.geometry.util import convex_hull

pts = [Point(0,0), Point(2,0), Point(2,2), Point(0,2), Point(1,1), Point(1,-1)]
h = convex_hull(*pts)
print(type(h))
print(h)
# 預期輸出：
# <class 'sympy.geometry.polygon.Polygon'>
# Polygon(Point2D(0, 0), Point2D(1, -1), Point2D(2, 0), Point2D(2, 2), Point2D(0, 2))

幾個重點：

參數是一串點（*args），不是一個列表；
回傳的是 Polygon——已經是逆時針（counterclockwise）排列好的頂點；
頂點沒有重複、也不包含尾端點；
逆時針是 SymPy 內部的慣例，本書所有幾何慣例一致。

退化情況怎麼處理？#

點數太少或有共線點時，SymPy 會「降維」回傳更基本的幾何物件：

# 08.4_convex_hull_degenerate.py
from sympy.geometry import Point, Segment
from sympy.geometry.util import convex_hull

# 三點共線：凸包退化為線段
print(convex_hull(Point(0,0), Point(1,0), Point(2,0), Point(1,1)))
# 預期輸出：
# Polygon(Point2D(0, 0), Point2D(2, 0), Point2D(1, 1))

# 完全共線的 4 點：convex_hull 會回傳 Segment
print(convex_hull(Point(0,0), Point(1,0), Point(2,0), Point(3,0)))
# 預期輸出：
# Segment2D(Point2D(0, 0), Point2D(3, 0))

# 只有一個點
print(convex_hull(Point(5,5)))
# 預期輸出：
# Point2D(5, 5)

“降維” 的行為值得記住：2D 的凸包只會回傳 Polygon、Segment 或 Point，型別由退化程度決定。

它是不是「Graham scan」？#

不是。SymPy 的 convex_hull 內部使用「gift wrapping（Jarvis march）」的變體（原始碼位置：sympy.geometry.util._find_hull_and_maybe_normals_...），複雜度為 \(O(nh)\) （ \(h\) 為凸包頂點數）。這個選擇對 SymPy 是合理的——它每屏一次都做精確運算，代價在正確性；但 \(n\) 型資料量大時會慢。

因此本節的答案是：SymPy 有 convex_hull，但它不夠快。所以以下兩條路線：

自己寫 Graham scan / Monotone chain（8.2 已寫過），用純 Python 或 NumPy 向量化；
用 convex_hull 做小案情或驗證（確認給定輸出無誤）。

8.4.2 自己寫 Graham scan，與 SymPy 交叉驗證#

# 08.4_cross_check.py
import math
import random
from sympy.geometry import Point
from sympy.geometry.util import convex_hull

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

def monotone_chain(points):
    """Andrew 1979 monotone chain，O(n log n)"""
    pts = sorted(set(points))
    if len(pts) <= 2:
        return pts
    lower, upper = [], []
    for p in pts:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        lower.append(p)
    for p in reversed(pts):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)
    return lower[:-1] + upper[:-1]

def np_to_sympy(hull):
    return convex_hull(*[Point(*p) for p in hull])

def same_hull(h1, h2):
    """兩個凸包是否含同一組頂點（忽略順序，無法完全相等時為 degenerate）"""
    a = {tuple(p) if hasattr(p, "__iter__") else (float(p.x), float(p.y)) for p in h1}
    b = {(float(p.x), float(p.y)) if hasattr(p, "x") else tuple(p) for p in h2}
    return a == b

random.seed(42)
N = 40
pts = [(random.randint(-20, 20), random.randint(-20, 20)) for _ in range(N)]
my_hull  = monotone_chain(pts)
sympy_hull = np_to_sympy(my_hull)
still = convex_hull(*[Point(*p) for p in pts])
print("我的 Monotone chain:", my_hull)
print("SymPy 的 convex_hull（由同一結果建出的多邊形）:", sympy_hull)
print("SymPy 原始輸入:", still)
print("頂點集合一致嗎?" , same_hull(my_hull, still.vertices if hasattr(still, "vertices") else still))
# 預期輸出（最多不同是「點的順序」，集合應當一致）：
# 頂點集合一致嗎? True

8.4.3 NumPy 向量化：一行的 cross 矩陣#

做幾何到大數據階段，重點是讓 Python 的迴圈消失。觀察：對所有點做「判斷 \(p\) 相對於直線 \((a, b)\) 的方向」這件事，本質上就是矩陣乘法。

考慮直線 \(a \to b\) 與所有點 \(p_i\) ：

\[ t_i = (b_x - a_x)(p_{i,y} - a_y) - (b_y - a_y)(p_{i,x} - a_x) \]

這可以寫成 NumPy 的 dot product：

# 08.4_numpy_cross.py
import numpy as np

def cross_all(a, b, pts):
    """對每個 pts 回傳 (b-a) × (p-a)。1 維矩陣版。"""
    a = np.asarray(a, dtype=float); b = np.asarray(b, dtype=float)
    pts = np.asarray(pts, dtype=float)
    d = b - a                       # (2,)
    return d[0] * (pts[:, 1] - a[1]) - d[1] * (pts[:, 0] - a[0])   # (n,)

pts = np.array([[0,0],[4,0],[0,4],[4,4],[2,2]])
print(cross_all([0,0], [4,0], pts))
# 預期輸出：
# [  0。   0。  16。  16。   8。]

上面的輸出是作者機器的順序（第一位是 0.0；留意空白）。它告訴你：「相對於一行」的狀態可以一個運算完成，這是 8.2 節單點迴圈的 \(n\) 倍速度等等。

這在處理半平面判定、在線段同側的判斷、還有Voronoi cell 頂點檢查都有用：你可以一次把 \(n\) 個候選點的 side 批次得出。

向量化的凸包（等幅 Monotone chain 的 NumPy 版本）#

# 08.4_numpy_monotone.py
import numpy as np

def monotone_chain_np(pts):
    """Andrew monotone chain 的 NumPy 輔助版（主體仍是 Python 迴圈，但排序與 side 檢查為 NumPy）"""
    pts = np.asarray(pts, dtype=float)
    pts = pts[np.lexsort((pts[:,1], pts[:,0]))]
    # 去重：保留第一個不同 x,y
    if len(pts) <= 2:
        return pts
    cross = lambda o, a, b: (a[0]-o[0])*(b[1]-o[1]) - (a[1]-o[1])*(b[0]-o[0])
    lower, upper = [], []
    for p in pts:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        lower.append(p)
    for p in pts[::-1]:
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)
    return np.array(lower[:-1] + upper[:-1])

random_pts = np.random.default_rng(1).integers(-50, 50, size=(100, 2))
print(monotone_chain_np(random_pts)[:5])
# 預期輸出（前 5 列，順序為逆時針）：
# ...數值隨亂數種子而定

當你想真的快，直接上 SciPy：scipy.spatial.ConvexHull 背後是以 QuickHull 為核心的 qhull；那是第 10 章工業級工具的主軸。

8.4.4 用 matplotlib 畫 Voronoi 圖#

本節的主菜：自己用 NumPy + matplotlib 畫 Voronoi 圖，而不是依賴現成套件。

步驟#

對每對站 \((p_i, p_j)\) ，算垂直平分線（SymPy 精確做一次即可）；
對站 \(p_i\) 的 cell，做半平面裁剪：把當前的多邊形（初始為超大的矩形），裁掉「不屬於站 \(i\) 」的那一側；
顯示 cell 與站。

演算法：多邊形的半平面裁剪（Sutherland–Hodgman 本質）#

虛擬碼：

clip(poly, a, b):       # 保留位於 (a,b) 左側的部分
    out = []
    for p in poly 依序：
        t = side(a, b, p)
        if t >= 0:
            # p 在內側
            if prev_t < 0:
                out.append(線段 (a,b) 與 (prev, p) 的交點)
            out.append(p)
        else:
            if prev_t >= 0:
                out.append(線段 (a,b) 與 (prev, p) 的交點)
        prev_t, prev = t, p
    return out

重複對每個「其他站」的 bisector 呼叫 clip，得到 cell。

完整可執行程式#

# 08.4_plot_voronoi.py
import math
import numpy as np
import matplotlib
matplotlib.use("Agg")        # 無視窗環境可存檔
import matplotlib.pyplot as plt

def side(a, b, p):
    """回傳 (b - a) × (p - a)，>0 在左、=0 在線上、<0 在右"""
    return (b[0]-a[0]) * (p[1]-a[1]) - (b[1]-a[1]) * (p[0]-a[0])

def seg_line_intersect(a, b, p1, p2):
    """直線 (a,b) 與線段 (p1, p2) 的交點（假設跨過）"""
    x1, y1 = p1; x2, y2 = p2
    dx1, dy1 = b[0]-a[0], b[1]-a[1]
    dx2, dy2 = x2-x1, y2-y1
    den = dx1*dy2 - dy1*dx2
    t = ((x1-a[0])*dy2 - (y1-a[1])*dx2) / den
    return (a[0] + t*dx1, a[1] + t*dy1)

def clip(poly, a, b):
    """保留 poly 中位於直線 (a,b) 左側的部分（半平面裁剪）"""
    out = []
    n = len(poly)
    for i in range(n):
        cur, prev = poly[i], poly[i-1]
        sc, sp = side(a, b, cur), side(a, b, prev)
        if sc >= 0:
            if sp < 0:
                out.append(seg_line_intersect(a, b, cur, prev))
            out.append(cur)
        elif sp >= 0:
            out.append(seg_line_intersect(a, b, cur, prev))
    return out

def voronoi_cell_np(sites, i, box):
    """站 i 的 Voronoi cell：從大矩形 box 開始，被所有其他站裁剪"""
    p = sites[i]
    poly = list(box)
    for j, q in enumerate(sites):
        if j == i: continue
        mx, my = (p[0]+q[0])/2, (p[1]+q[1])/2     # 中點
        dx, dy = q[0]-p[0], q[1]-p[1]             # p→q 方向
        # 平分線上的兩點 a, b（a = 中點, b = 中點 + 旋轉 90°），
        # 「左側」= 靠 p 的那一半平面
        a = (mx, my); b = (mx - dy, my + dx)
        if side(a, b, p) < 0:                     # 保證 p 在左側
            a, b = b, a
        poly = clip(poly, a, b)
        if not poly:
            break
    return poly

# 站集
sites = [(0.0, 0.0), (4.0, 0.0), (0.0, 4.0), (4.0, 4.0), (2.0, 2.0),
         (7.0, 2.0), (2.0, 7.0)]
box = [(-2, -2), (9, -2), (9, 9), (-2, 9)]        # 視野用的巨矩形

fig, ax = plt.subplots(figsize=(6, 6))
colors = plt.cm.tab10.colors
for i, s in enumerate(sites):
    cell = voronoi_cell_np(sites, i, box)
    if len(cell) >= 3:
        xs = [c[0] for c in cell] + [cell[0][0]]
        ys = [c[1] for c in cell] + [cell[0][1]]
        ax.fill(xs, ys, alpha=0.35, color=colors[i % 10])
        ax.plot(xs, ys, color="k", lw=1.5)
    ax.plot(*s, "k*", markersize=12, zorder=5)

ax.set_xlim(-2.5, 9.5); ax.set_ylim(-2.5, 9.5)
ax.set_title("Voronoi diagram（半平面裁剪法繪製）")
ax.set_aspect("equal")
plt.savefig("08.4_voronoi.png", dpi=120)
print("saved 08.4_voronoi.png")
# 預期輸出：
# saved 08.4_voronoi.png
# 圖中 7 個站，中央站 (2,2) 與 (7,2)、(2,7) 有界 cell，
# 其他站（貼邊、在框上）的 cell 被矩形 box 截斷。

為什麼用「巨矩形 box」？ 因為 cell 通常無界（8.3 節講過）。畫圖時的標準做法：先給一個「全世界只要我畫得下就夠大」的矩形，再對這個矩形裁剪。這就是 bounded clip。

為什麼值增加比對 8.3 版本#

8.3 節	8.4 節
用砨組密封公用的檢驗法	用裁剪（一般用法，能畫圖、能上電腦圖學）
只得出 cell 的「部分頂點」	得出視野範圍內完整的 cell
SymPy 精確、慢	NumPy 浮點、快
未使用 matplotlib 畫 face	完整可視化

兩種版本都是正確的，則適用在「用數學思考」與「讓程式跑起來」的不同場景。

8.4.5 用 SciPy 把 Voronoi / Delaunay 一次到位（工程首選）#

真實項目裡，若已進入 SciPy 生態，不必自己寫：

# 08.4_scipy_quick.py
import numpy as np
from scipy.spatial import Voronoi, Delaunay

rng = np.random.default_rng(3)
pts = rng.random((30, 2))
vd = Voronoi(pts)
dt = Delaunay(pts)
print("Voronoi vertices:", vd.vertices.shape)
print("Voronoi regions 數:", len(vd.regions), "（某些是空 list，表示無界 cell）")
print("Delaunay simplices:", dt.simplices.shape)
# 預期輸出（30 個隨機站、數值會浮動）：
# Voronoi vertices: (100+, 2)（包括 box 外的潛在頂點）
# Voronoi regions 數: 32（最後一項是∞區）...
# Delaunay simplices: (60左右, 3)
# 注意：當輸入沒有共圓時，地面簡單的 simplices 介於 2n-2-h 與其他細節關聯

但 SciPy 是浮點、C 程式；SymPy 是精確、Python。工程上對「真地要做」的東西加上 SciPy；對「要證明」的東西加上 SymPy。

8.4.6 用「定義」驗證「程式」：500 次隨機測試#

程式寫完最怕的是「看起來對」。驗證的黃金法則是：用數學定義（最近的那站）去檢驗程式輸出（cell 的歸屬）。做法：隨機撒 500 個測試點，每個點分別用兩種判斷——

定義判斷：用 SymPy 精確算它到 7 個站的距離，取最小者；
程式判斷：用 point_in_poly（同側法，適用凸多邊形）測它落在哪個 cell。

兩者必須一致。先補一個 point_in_poly：

# 08.4_verify_voronoi.py  （共用的 point_in_poly）
def point_in_poly(poly, p, eps=1e-9):
    """凸多邊形內含檢查：True＝內部，False＝外部，字串＝邊上"""
    sgn = 0
    n = len(poly)
    for i in range(n):
        s = side(poly[i], poly[(i + 1) % n], p)     # side 來自 8.4.4
        if abs(s) < eps:
            return "on-edge"
        cur = 1 if s > 0 else -1
        if sgn == 0:
            sgn = cur
        elif sgn != cur:
            return False                              # 不同側 → 外部
    return True

再把兩種判斷接起來（請將 8.4.4 的 side、seg_line_intersect、clip、voronoi_cell_np 與下面放同一檔案）：

# 08.4_verify_voronoi.py （主體）
import numpy as np
from sympy.geometry import Point

sites = [(0.0,0.0), (4.0,0.0), (0.0,4.0), (4.0,4.0),
         (2.0,2.0), (7.0,2.0), (2.0,7.0)]
box = [(-2,-2), (9,-2), (9,9), (-2,9)]
cells = [voronoi_cell_np(sites, i, box) for i in range(len(sites))]
sps = [Point(*s) for s in sites]

rng = np.random.default_rng(0)
ok = bad = 0
for _ in range(500):
    x, y = rng.uniform(-2, 9, size=2)
    p = Point(float(x), float(y))
    nearest = min(range(len(sps)), key=lambda i: float(sps[i].distance(p)))
    r = point_in_poly(cells[nearest], (x, y))
    if r is True or r == "on-edge":
        ok += 1
    else:
        bad += 1
print("500 次隨機測試：一致 =", ok, "；不一致 =", bad)
# 預期輸出：
# 500 次隨機測試：一致 = 500 ；不一致 = 0

為什麼「不一致」恆為 0？因為 Voronoi cell 的定義就是「以該站為最近」的點集——只要裁剪與內含判斷都正確，兩種判斷必然吻合。這種「定義 vs 程式」的一致性測試是本書一貫的技倆：程式可以錯，定義不會錯。任何改動之後重跑 500 次測試，是計算幾何程式的最低門檻。

8.4.7 你應該記下的函數表#

目的	SymPy	NumPy/SciPy
凸包	`sympy.geometry.util.convex_hull(*pts)`	`scipy.spatial.ConvexHull`
Delaunay	自建（無內建）	`scipy.spatial.Delaunay`
Voronoi	自建（無內建）	`scipy.spatial.Voronoi`
垂直平分線	`Line(p,q).perpendicular_line(m)`	手算 \((b-a) \to (-dy, dx)\)
半平面交	自建 poly clip	`scipy.spatial.HalfspaceIntersection`
精確運算，字面數學	✓	✗
大數據速度	✗	✓

8.4.8 習題#

（API） convex_hull 的回傳值為什麼可能不是 Polygon？請用至少 2 種退化輸入舉例（例如全部共線、只有一點）並展示。
（程式 A） 把 8.2 節的 monotone_chain 改成 NumPy 進一步向量化 的版本：把「while 彈出堆疊」依然用 Python，但排序與資料形態改為 np.lexsort 得出後的地圖結果；與本節 monotone_chain_np 在 \(n = 20000\) 隨機點上計時，寫出結果。
（程式 B） 在 08.4_plot_voronoi.py 基礎上，加入在每個 cell 中央顯示站名的功能：把 sites 改成 dict（站名 → 座標），並在畫完 cell 之後 ax.text 展示幾個站名。觀察站名與 cell 的關聯。
（演算法） 在 clip 函式中：輸入多邊形有 \(n\) 條邊；證明 clip 本身是 \(O(n)\) ；這對計算站 \(i\) 的 cell 至多 \(O(n^2)\) 作證——確認當站數總體很大時，這不是產品級做法。
（對偶實作） 取 Voronoi 的對偶：已知每個站對之間的 bisector 被共享 ⟺ 對應 Delaunay 邊。從你 8.4.4 的 voronoi_cell_np 輸出，訪察「這對站孰站在 cell 邊上出現」以推回 Delaunay 邊集合；與 scipy.spatial.Delaunay 的 simplices 集合進行驗證（小位置集即可）。

到這裡，第 8 章「計算幾何」結束。你已能：解釋這門學科為何在 1970 年代誕生（8.1）、實作四種凸包演算法（8.2）、理解 Voronoi–Delaunay 對偶與 regular triangulation（8.3）、並把 SymPy 精確性與 NumPy 速度組合起來（8.4）。下一章（第 9 章）將走入微積分的幾何：高斯的內蘊幾何與曲率。