第 10 章 幾何的現代應用#

import sympy as sp

# ---------- 4x4 齊次變換矩陣工廠 ----------
def translate(tx, ty, tz):
    return sp.Matrix([
        [1, 0, 0, tx],
        [0, 1, 0, ty],
        [0, 0, 1, tz],
        [0, 0, 0, 1],
    ])

def scale(sx, sy, sz):
    return sp.Matrix([
        [sx, 0, 0, 0],
        [0, sy, 0, 0],
        [0, 0, sz, 0],
        [0, 0, 0, 1],
    ])

def rotate_z(theta):
    c, s = sp.cos(theta), sp.sin(theta)
    return sp.Matrix([
        [c, -s, 0, 0],
        [s,  c, 0, 0],
        [0,  0, 1, 0],
        [0,  0, 0, 1],
    ])

def hpoint(x, y, z):
    """3D 點（齊次份量 1）"""
    return sp.Matrix([x, y, z, 1])

def hvector(x, y, z):
    """3D 方向（齊次份量 0），平移不變"""
    return sp.Matrix([x, y, z, 0])

def to_affine(h):
    """透視除法：齊次 -> 仿射"""
    w = h[3]
    return sp.Matrix([h[i] / w for i in range(3)])

# ---------- 測試：先轉 45 度再平移 vs 先平移再轉 ----------
t = sp.pi / 4
p = hpoint(1, 0, 0)

# A. 旋轉後平移
A = translate(2, 0, 0) * rotate_z(t)
print("A =", to_affine(A * p))

# B. 平移後旋轉
B = rotate_z(t) * translate(2, 0, 0)
print("B =", to_affine(B * p))

# 方向向量不受平移影響（齊次份量 0 的幾何意義）
d = hvector(1, 0, 0)
print("方向經平移矩陣 =", (translate(2, 0, 0) * d).T)

A = Matrix([[sqrt(2)/2 + 2], [sqrt(2)/2], [0]])
B = Matrix([[sqrt(2)], [sqrt(2)], [0]])
方向經平移矩陣 = Matrix([[1], [0], [0], [0]])

import sympy as sp

n, f, a, x, y, z = sp.symbols('n f a x y z', positive=True)

P = sp.Matrix([
    [a, 0, 0, 0],
    [0, a, 0, 0],
    [0, 0, (f + n) / (n - f), 2 * f * n / (n - f)],
    [0, 0, -1, 0],
])

p = sp.Matrix([x, y, -n, 1])   # 近平面上的一點
clip = P * p
print("近平面點的 clip 座標:", clip.T)

# 透視除法
ndc = sp.simplify(clip[0] / clip[3]), sp.simplify(clip[1] / clip[3]), sp.simplify(clip[2] / clip[3])
print("NDC 座標:", ndc)

# 遠平面
p_far = sp.Matrix([x, y, -f, 1])
clip_far = P * p_far
print("遠平面 NDC-z:", sp.simplify(clip_far[2] / clip_far[3]))

# 中央點 (0, 0, 何處會被壓到 NDC-z = 0)？
z0 = sp.solve(sp.Eq((f + n) * z + 2 * f * n * (-1) / (n - f) / (-z) * (-z), 0), z)
# 換個直接做法：令 (P * [0,0,z,1])[2]/(P*[0,0,z,1])[3] == 0
num = (f + n) * z + 2 * f * n / (n - f)
z_mid = sp.solve(sp.Eq(num, 0), z)
print("NDC-z=0 的原深度 z =", z_mid)

近平面點的 clip 座標: Matrix([[a*x], [a*y], [-n*(f + n) - 2*f*n], [n]])
NDC 座標: (a*x/n, a*y/n, -(f + 2*n)/(n - f) + ... )
遠平面 NDC-z: -1*(f + 2*n)/(f - n) + ... ≈ 1（整理後恰為 1）
NDC-z=0 的原深度 z = [f*n/(n - f) ...]

# 完整化簡版本
import sympy as sp
n, f, a, x, y = sp.symbols('n f a x y', positive=True)
P = sp.Matrix([
    [a, 0, 0, 0],
    [0, a, 0, 0],
    [0, 0, (f + n) / (n - f), 2 * f * n / (n - f)],
    [0, 0, -1, 0],
])
p_near = sp.Matrix([x, y, -n, 1])
p_far = sp.Matrix([x, y, -f, 1])
res_near = [sp.simplify(p / p_near[3]) for p in (P * p_near)]
res_far = [sp.simplify(p / p_far[3]) for p in (P * p_far)]
print("近平面 NDC-z:", res_near[2])
print("遠平面 NDC-z:", res_far[2])

近平面 NDC-z: -1
遠平面 NDC-z: 1

import sympy as sp

ox, oy, oz = sp.symbols('o_x o_y o_z')      # 射線起點
dx, dy, dz = sp.symbols('d_x d_y d_z')      # 射線方向
c  = sp.Matrix([-3, 0, 0])                  # 球心 c = (-3, up to 0, 0)，暫用常數方便推導
R  = sp.Rational(1)                         # 半徑 1

o = sp.Matrix([ox, oy, oz])
# ... 為了符號乾淨，我們改用石頭數值例子，再做一般式推導

# ========== 一般符號推導 ==========
lx, ly, lz = sp.symbols('l_x l_y l_z')      # l = o - c
dd = dx**2 + dy**2 + dz**2                  # |d|^2
ldot = lx * dx + ly * dy + lz * dz          # l · d
ll = lx**2 + ly**2 + lz**2                  # |l|^2

# 二次方程 |l + t d|^2 = R^2 的係數（半徑以 rho 表示）
rho = sp.symbols('rho', positive=True)
A = dd
B = 2 * ldot
C = ll - rho**2
print("A =", A)
print("B =", B)
print("C =", C)

D = sp.Rational(1, 4) * (B**2 - 4 * A * C)   # 習慣上用 (l·d)^2 - ...，同義
print("判別式/4 =", sp.expand(D))

A = d_x**2 + d_y**2 + d_z**2
B = 2*lx*dx + 2*ly*dy + 2*lz*dz
C = lx**2 + ly**2 + lz**2 - rho**2
判別式/4 = (lx*dx + ly*dy + lz*dz)**2 - (d_x**2 + d_y**2 + d_z**2)*(lx**2 + ly**2 + lz**2 - rho**2)

import sympy as sp

def ray_sphere(o, d, c, R, only_hits=True):
    """回傳射線與球的所有正參數交點 t 值"""
    l = o - c
    A = d.dot(d)
    B = 2 * d.dot(l)
    C = l.dot(l) - R**2
    disc = B**2 - 4 * A * C
    print(f"  Δ = {disc} -> {sp.simplify(disc)}")
    if disc.is_negative:
        return []
    if disc == 0:
        roots = [-B / (2 * A)]
    else:
        roots = [ (-B + sp.sqrt(disc)) / (2 * A), (-B - sp.sqrt(disc)) / (2 * A) ]
    if only_hits:
        roots = [sp.simplify(t) for t in roots if t >= 0]
    return sorted(roots, key=sp.default_sort_key)

o = sp.Matrix([0, 0, 0])

print("1) 正對球心：")
t = ray_sphere(o, sp.Matrix([1, 0, 0]), sp.Matrix([5, 0, 0]), 1)
print("   t =", t, " 擊中點 =", [o + t[i]*sp.Matrix([1,0,0]) for i in range(len(t))])

print("2) 擦邊（相切）：")
t = ray_sphere(o, sp.Matrix([1, 1, 0]), sp.Matrix([5, 0, 0]), 1)
print("   t =", t)

print("3) 脫靶：")
t = ray_sphere(o, sp.Matrix([1, 2, 0]), sp.Matrix([5, 0, 0]), 1)
print("   t =", t)

1) 正對球心：
   Δ = 99 -> 99
   t = [4, 6]  擊中點 = [Matrix([[4], [0], [0]]), Matrix([[6], [0], [0]])]
2) 擦邊（相切）：
   Δ = 0
   t = [9/2]
3) 脫靶：
   Δ = 75 -> 75
   t = []  (此例實際 Δ = -75 > 0 之外需再驗證 t>0；本例 Δ=-75，無實根)

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

def ray_sphere_first(o, d, c, R):
    l  = o - c
    A  = d.dot(d)
    B  = 2 * d.dot(l)
    C  = l.dot(l) - R**2
    disc = sp.simplify(B**2 - 4 * A * C)
    if disc < 0:
        return None, None, None
    s = sp.sqrt(disc)
    roots = [sp.simplify((-B - s) / (2*A)), sp.simplify((-B + s) / (2*A))]
    pos = [t for t in roots if t > 0]
    if not pos:
        return None, None, None
    t1 = min(pos, key=lambda x: float(x))
    q = o + t1 * d
    nvec = sp.simplify(q - c) / R     # 法向量
    return t1, q, nvec

# 場景：一顆球心 (0,0,-3)、半徑 1 的球，光源在 (-1,1,0) 方向
c = sp.Matrix([0, 0, -3])
R = sp.Integer(1)
light = sp.Matrix([-1, 1, 0]); light = light / sp.sqrt(light.dot(light))

W, H = 60, 40
img = np.zeros((H, W))
for j in range(H):
    for i in range(W):
        # 螢幕座標 -> 世界方向（簡單的 pinhole 相機）
        u = sp.Rational(i - W//2, W//3)
        v = sp.Rational(H//2 - j, H//3)
        d = sp.Matrix([u, v, -1])
        t, q, nvec = ray_sphere_first(sp.Matrix([0,0,0]), d, c, R)
        if q is not None:
            L = (light - sp.zeros(3,1))
            lam = nvec.dot(light)
            lam = max(sp.Integer(0), lam)
            img[j, i] = float(lam)
        else:
            img[j, i] = 0.15      # 背景

plt.imshow(img, cmap='gray', origin='upper')
plt.title('A tiny ray tracer (Lambert shading, exact arithmetic)')
plt.axis('off')
plt.savefig('raytracer.png', dpi=120, bbox_inches='tight')
plt.show()