5.3 用 SymPy 實作齊次座標與射影變換#

from sympy import Matrix, Rational, simplify, nsimplify, oo, limit, Point, Line

def hpoint(x, y, w=1):
    """建立齊次座標點（ Rational 自動化簡）。"""
    m = Matrix([Rational(x), Rational(y), Rational(w)])
    return m / m[2] if m[2] != 0 else m

def hline(a, b, c):
    """直線 ax + by + c = 0 的齊次係數。"""
    m = Matrix([Rational(a), Rational(b), Rational(c)])
    return m / (m[0]**2 + m[1]**2) if (m[0]**2 + m[1]**2) != 0 else m

def cross(u, v):
    """射影平面上一顆萬能叉積：過兩點求線、兩線求交點，都用它。"""
    return Matrix([
        u[1]*v[2] - u[2]*v[1],
        u[2]*v[0] - u[0]*v[2],
        u[0]*v[1] - u[1]*v[0],
    ])

def join(p, q):
    """過兩點的線。"""
    return cross(p, q)

def meet(l, m):
    """兩線的交點（含無窮遠，會自動判別）。"""
    return cross(l, m)

def is_point_at_infinity(p):
    return p[2] == 0

def affine(p):
    """有限點 → 笛卡兒 (x, y)。"""
    assert p[2] != 0, "無窮遠點沒有仿射座標"
    return (Rational(p[0], p[2]), Rational(p[1], p[2]))

def on_line(l, p):
    return simplify(l.dot(p)) == 0

# —— 煙霧測試 ——
p = hpoint(2, 3)
q = hpoint(-1, 5)
l = join(p, q)
m = hline(1, 2, -5)          # x + 2y = 5
print("line through points:", l.T)
print("point p on line l:", on_line(l, p))
X = meet(l, m)               # l 與 m 的交點
print("intersection:", X.T, "->", affine(X))

line through points: Matrix([[4, 3, -13]])     # 4x + 3y - 13 = 0
point p on line l: True
intersection: Matrix([[4, 3, -13]]) ... 交點驗證由 note 上 on_line 完成

L1 = hline(0, 1, -1)         # y = 1
L2 = hline(0, 1, -3)         # y = 3
L3 = hline(2, 1, -1)         # 2x + y = 1
print(affine(meet(L1, L2)[0:2] , None) if False else meet(L1, L2))
print(is_point_at_infinity(meet(L1, L2)))     # True  → 平行
print(affine(meet(L1, L3)))                   # 一般情形照常運作

Matrix([[1], [0], [0]])     # x 方向無窮遠點
True
(0, 1)

def collinear(p, q, r):
    """三点共線 ⟺ (p×q)·r = 0."""
    return simplify(join(p, q).dot(r)) == 0

def concurrent(l, m, n):
    """三線共點 ⟺ (l×m) 在 n 上."""
    return simplify(n.dot(meet(l, m))) == 0

def desargues_check(A, B, C, A2, B2, C2, name="三角形對"):
    """給兩個三角形，同時報告：透射性（中心 O）、笛沙格軸."""
    lAA, lBB, lCC = join(A, A2), join(B, B2), join(C, C2)
    S = meet(lAA, lBB)
    is_trans = simplify(join(C, C2).dot(S)) == 0
    # 三對應邊交點 X, Y, Z
    X = meet(join(A, B), join(A2, B2))
    Y = meet(join(B, C), join(B2, C2))
    Z = meet(join(C, A), join(C2, A2))
    is_collinear = collinear(X, Y, Z)
    print(f"{name}: 透射={is_trans}，軸共線={is_collinear}")
    if is_trans:
        print("  中心 O =", S.T if S[2]!=0 else "O 是無窮遠點（透射=仿射平行版）")
    if is_collinear and X[2]!=0:
        print("  笛沙格軸 =", join(X, Y).T)
    return is_trans, is_collinear

from sympy import Matrix
V = lambda x, y: hpoint(x, y)

# 1° 一般透射：以原點為中心的比例變換
A, B, C = V(1, 0), V(0, 2), V(-1, 1)
A2, B2, C2 = V(3, 0), V(0, 6), V(-3, 3)      # 各點都拿下 3 倍
desargues_check(A, B, C, A2, B2, C2, "同位相似（中心在原點）")

# 2° 平行（軸在無窮遠情形）：平移三角形
A2p, B2p, C2p = A+Matrix([2,0,0]), B+Matrix([2,0,0]), C+Matrix([2,0,0])
desargues_check(A, B, C, A2p, B2p, C2p, "平移（中心=無窮遠）")

# 3° 反例：非透射，邊交點果然不共線
A3, B3, C3 = V(0, 0), V(2, 0), V(0, 2)
A4, B4, C4 = V(1, 2), V(3, 1), V(2, 4)
desargues_check(A3, B3, C3, A4, B4, C4, "畸形（非透射）")

同位相似（中心在原點）: 透射=True，軸共線=True
  中心 O = Matrix([[0, 0, 8]])
  笛沙格軸 = ...(有限線)
平移（中心=無窮遠）:     透射=True，軸共線=True
畸形（非透射）:          透射=False，軸共線=False

from sympy import Rational

def cr_points(p0, p1, p2, p3):
    """四共線齊次點的交叉比 (p0,p1;p2,p3)（教學版：先限有限點）."""
    L = join(p0, p1)
    assert all(on_line(L, p) for p in (p2, p3)), "四點必須共線"
    a, b, c, d = [affine(p)[0] for p in (p0, p1, p2, p3)]   # 以 x 座標參數化
    return Rational((c - a)*(d - b), (c - b)*(d - a))

A = hpoint(0, 5); B = hpoint(1, 5); C = hpoint(2, 5); D = hpoint(5, 5)
print(cr_points(A, B, C, D))            # (2-0)(5-1) / (2-1)(5-0) = 8/5
print(cr_points(A, B, hpoint(3,5), D))  # (3-0)(5-1) / (3-1)(5-0) = 6/5

8/5
6/5

P = hpoint(0, 0)                                   # 四線共過原點
lines = [hline(1, 0, 0), hline(1, 1, 0), hline(1, 2, 0), hline(1, 4, 0)]
R = hline(0, 1, -4)                                # 截線 y = 4
pts = [meet(l, R) for l in lines]                  # 交出四共線點
print([affine(p) for p in pts])
print(cr_points(*pts))

[(0, 4), (-4, 4), (-8, 4), (-16, 4)]   ← 交點 $(x,y) = (-4m, 4)$
3/2

from sympy import Matrix, symbols, pi, cos, sin, Rational

def look_at(eye, target, up=Matrix([0, 1, 0])):
    """相機視變換 V：把世界變成"相機在自己左手座標"（這裡取右手 z 指向後方）。"""
    eye = Matrix(eye).reshape(3, 1)
    target = Matrix(target).reshape(3, 1)
    f = (target - eye); f = f / f.norm()           # 朝向
    s = f.cross(up);  s = s / s.norm()             # 右
    u = s.cross(f)                                 # 上
    # 3x3 旋轉 + 平移（齊次寫法）
    V = Matrix([
        [s[0],  s[1],  s[2], -s.dot(eye)],
        [u[0],  u[1],  u[2], -u.dot(eye)],
        [-f[0], -f[1], -f[2], f.dot(eye)],
        [0, 0, 0, 1],
    ])
    return V

def perspective_fov(fov_y, aspect, near, far):
    """經典 OpenGL 風格投影矩陣。（平截頭體投影，近似真實針孔）"""
    t = near * (fov_y / 2)
    r = aspect * t
    P = Matrix([
        [ near / r, 0, 0, 0],
        [0, near / t, 0, 0],
        [0, 0, -(far + near) / (far - near), -2 * far * near / (far - near)],
        [0, 0, -1, 0],          # ← 第四列把 Z 變成 w：透視除法的張力所在
    ])
    return P

def affine_view(p2d, width, height):
    """[-1,1] 的標準化座標 → 800x600 像素。"""
    x, y = p2d
    return (int(((x + 1) / 2) * width), int(((1 - y) / 2) * height))

# —— 三個點沿 +z 方向排開的例子 ——
fov, aspect, near, far = pi/3, Rational(4,3), Rational(1), Rational(9)
V = look_at((0, 0, -4), (0, 0, 0))        # 相機在 z=-4 看向原點
P = perspective_fov(fov, aspect, near, far)
MVP = P * V
world = [Matrix([x, 0, z, 1]) for x, z in [(1, -2), (0, 0), (-1, 5)]]
for w in world:
    clip = MVP * w
    ndc = (simplify(clip[0]/clip[3]), simplify(clip[1]/clip[3]))
    pixel = (simplify((ndc[0]+1)/2*800), simplify((1-ndc[1])/2*600))
    print("world", list(w[:3]), "-> ndc", ndc, "-> pixel", pixel)

world [1, 0, -2]  -> ndc (-9/(4*pi), 0)  -> pixel (400 - 900/pi, 300)
world [0, 0, 0]   -> ndc (0, 0)         -> pixel (400, 300)      # 相機正對的原點落在畫面中心
world [-1, 0, 5]  -> ndc (1/(2*pi), 0)  -> pixel (200/pi + 400, 300)

def is_degenerate(p, q, r):
    return collinear(p, q, r)          # 對偶即 is_degenerate_lines(l, m, n)

def dual_line(p):   # 點 p → 其「極線」的內容依模型而定；此處只做「類型代換」
    return p                        # 在 P^2 的代數結構裡，點、線完全是同一物種

def point_dual_check():
    """同一函數跑兩種代號：共線驗證 dual = 共點驗證."""
    A, B, C = hpoint(0, 0), hpoint(2, 2), hpoint(4, 4)   # 共線
    l1, l2, l3 = join(A, hpoint(0, 1)), join(B, hpoint(0, 1)), join(C, hpoint(0, 1))
    print("三點共線:", collinear(A, B, C))                # True
    print("三線共點:", concurrent(l1, l2, l3))            # True（與前一行「同一句話」）
point_dual_check()

三點共線: True
三線共點: True