Symbolic Rendering
Chapter 2 — Half-Vectors, Three Channels, Two Lights
An explainer comparing numerical rendering with
symbolic rendering,
in the spirit of Utah graphics meets Macsyma (today: Python + sympy).
in the spirit of Utah graphics meets Macsyma (today: Python + sympy).
Dr. L. Van Warren — Original Idea & Mathematical Development
Kimi K3 — Assistance with Mathematical Development
Kimi K3 — Assistance with Mathematical Development
❦
How to use this file.
§1 states exactly what changes relative to Chapter 1 — and what provably does not.
§2–§3 extend the derivation: the half-vector, three channels, two lights.
§4 is the live simulator: numeric and symbolic RGB renderers side by side, difference
underneath, with per-channel accounting.
§5 is the laboratory: eight Python cells (sympy + numpy), each with its three run paths.
§6–§9: the ledger deltas, new pitfalls, the history, and the road to shadows.
☰ contents
§1 · What changes in Chapter 2
The sphere, the camera, the quadratic, the hit point, the normal, the silhouette — unchanged. Chapter 1's Δ, its closed-form $t$, and the exact circle $\rho = eR/\sqrt{e^2-R^2}$ carry over verbatim, because geometry is upstream of shading. Three extensions, in order of increasing symbolism:
- Specular model: Phong's reflection term $r\cdot\hat V$ is joined by Blinn's half-vector term $N\cdot\hat H$ — a one-line change in the derivation, a different highlight in the image.
- One channel → three: the image is now three expressions $I_r, I_g, I_b$
sharing all of their geometry — and the sharing itself becomes a measurable object
(
sympy.cse). - One light → a symbolic Σ: two point lights with positions and colors as free parameters. Disabling a light is not a branch — it is the substitution $c_2 \to 0$.
§2 · The scene (Chapter-2 assumptions)
Chapter 1's scene, plus: material diffuse color $m \in [0,1]^3$; ambient tint $a \in [0,1]^3$; specular color white (bakeable later); two point lights $L_1, L_2$ with colors $c_1, c_2 \in [0,1]^3$; a model switch $\beta \in \{\text{Phong}, \text{Blinn}\}$; linear-space rendering with the display clamp $\min(1, \cdot)$ applied at the end (the gamma question is deferred, honestly, to Chapter 5). Background remains 0; no shadows — that is Chapter 3's quantifier.
§3 · The math, once, carefully
Intersection exactly as before: $A = x^2+y^2+e^2$, $\Delta = e^2R^2 - (e^2-R^2)(x^2+y^2)$, $t = (e^2-\sqrt{\Delta})/A$, $P = (tx, ty, e(1-t))$, $N = P/R$, $\hat V = (-x,-y,e)/\sqrt A$, $\nu = \sqrt{\Delta}/(R\sqrt A)$. For each light $\ell \in \{1,2\}$: $\hat w_\ell = (L_\ell - P)/\|L_\ell - P\|$ and $\mu_\ell = N\cdot\hat w_\ell$. The two specular geometries:
$$\text{Phong: } \quad \hat r\cdot\hat V = 2\mu_\ell\,\nu - \hat w_\ell\cdot\hat V, \qquad \text{spec}_\ell = \max(0,\; \hat r\cdot\hat V)^n \quad \text{(no new square root)}$$ $$\text{Blinn: } \quad \hat H_\ell = \frac{\hat w_\ell + \hat V}{\|\hat w_\ell + \hat V\|}, \qquad \text{spec}_\ell = \max(0,\; N\cdot\hat H_\ell)^n$$The image, one expression per channel $k \in \{r,g,b\}$:
$$I_k(x, y; \theta) = \begin{cases} k_a\,a_k + \displaystyle\sum_{\ell=1}^{2} c_{\ell k}\Big(k_d\,m_k\max(0, \mu_\ell) + k_s\,\text{spec}_\ell\Big), & \text{if } \Delta(x,y) \ge 0\\[6pt] 0, & \text{otherwise}\end{cases}$$with $\min(1, I_k)$ applied at display — a clamp the symbolic form carries explicitly and numeric pipelines apply silently.
Chapter 2's three gifts:
- One intersection serves three channels. All geometry — Δ, t, P, N, every
$\mu_\ell$ — is shared by $I_r, I_g, I_b$.
sympy.csemakes the sharing explicit and countable: the swell meter in Cell 2 shows the naive triple collapse back to barely more than one channel's cost. Common subexpression elimination is the symbolic renderer's oldest friend — it is Macsyma'soptimize(), reborn. - Disabling a light is substitution, not branching. The numeric renderer writes
if not enabled: continue; the symbolic renderer substitutes $c_2 \to 0$ into the same expression. Cell 5 proves the two agree to machine precision — you get a checkbox with no new logic. - The image is exactly affine in each light's color. $\partial^2 I_k / \partial c_{\ell k}^2 = 0$, provable in the CAS, and $\partial I_k / \partial c_{\ell k}$ is literally that light's illumination geometry. Color grading is linear algebra performed on one derivation, forever.
Counterpoint: three channels triple the naive swell, and Blinn's $\hat H$ costs one square root per light that Phong's closed form $2\mu_\ell\nu - \hat w_\ell\cdot\hat V$ avoids. But Blinn's classic win — $\hat H$ is constant for a directional light and orthographic viewer — is, in this framework, a provable special case: $\partial \hat H/\partial x = 0$, after which CSE hoists it out of the per-pixel work entirely. The CAS doesn't just compute the model; it can certify its simplifications.
§4 · Live simulator
The same sphere, now in color, under two lights. The symbolic panel contains no light loop — the Σ was unrolled at derivation time; the light-2 checkbox performs the substitution $c_2 \to 0$. Silhouette ρ = eR/√(e²−R²) in red, as ever.
model
resolution
initializing…
NUMERIC — the algorithm per pixel (runtime ifs, per light)
SYMBOLIC — the derived Σ evaluated per pixel
|difference| per channel, amplified ×10¹⁵ — structural agreement, now in triplicate
numeric: –
symbolic: –
max |ΔI| (any channel): –
ρ = –
hover the canvases — the probe shows pixel → world (x,y), Δ, branch, and both RGB triples.
Try this: (1) Uncheck light 2 — the gold fill vanishes and the difference canvas stays black: you just performed $c_2 \to 0$, a substitution, and the agreement never noticed. (2) Flip Phong ↔ Blinn — the highlight migrates; both panels move together, because both are the same kind of object. (3) Switch symbolic mode to Heaviside multiply (raw) — the background fills with magenta, per channel: $0 \cdot \mathrm{NaN} = \mathrm{NaN}$, in triplicate.
§5 · The Colab laboratory (Python / sympy)
Same discipline as Chapter 1: paste the cells in order, run in order;
every cell carries its three run paths — 📋 Copy code, ⬇ notebook
(upload via File → Upload notebook), or paste into
colab.new.
🐍 Cell 1 — the numerical baseline: RGB, two lights, Phong or Blinn
import numpy as np
import matplotlib.pyplot as plt
import time
def render_numeric_rgb(width, height, R=1.0, e=4.0,
lights=(((5., 5., 10.), (1., 1., 1.)),
((-6., 2., 8.), (0.9, 0.65, 0.25))),
ambient=(0.15, 0.15, 0.16), material=(0.63, 0.19, 0.16),
ka=0.1, kd=0.7, ks=0.5, shin=32.0, model='blinn'):
"""Classic numerical renderer, Chapter 2: RGB, two lights, Phong or Blinn.
Branches execute at runtime: per pixel, per light, per channel."""
img = np.zeros((height, width, 3))
E = np.array([0.0, 0.0, e])
half = 1.25 * e * R / np.sqrt(e**2 - R**2)
for j in range(height):
y = half * (1.0 - 2.0 * (j + 0.5) / height)
for i in range(width):
x = half * (2.0 * (i + 0.5) / width - 1.0)
D = np.array([x, y, -e])
A = float(D @ D); b = float(E @ D); c = e * e - R * R
disc = b * b - A * c
if disc < 0.0: # ---- branch: ray misses
continue
t = (-b - np.sqrt(disc)) / A
P = E + t * D
N = P / R
v_hat = (E - P) / np.linalg.norm(E - P)
col = ka * np.array(ambient, dtype=float)
for pos, color in lights:
if color is None: # ---- branch: light disabled
continue
cL = np.array(color, dtype=float)
w = np.array(pos, dtype=float) - P
w_hat = w / np.linalg.norm(w)
mu = float(N @ w_hat)
if model == 'phong':
r = 2.0 * mu * N - w_hat
spec = max(0.0, float(r @ v_hat)) ** shin
else: # blinn half-vector
h = w_hat + v_hat
h = h / np.linalg.norm(h)
spec = max(0.0, float(N @ h)) ** shin
col += cL * (kd * np.array(material) * max(0.0, mu) + ks * spec)
img[j, i] = col
return np.clip(img, 0.0, 1.0)
t0 = time.perf_counter()
img_num = render_numeric_rgb(256, 256)
t1 = time.perf_counter()
print(f"numeric RGB render: {t1 - t0:.3f} s")
plt.figure(figsize=(4, 4))
plt.imshow(img_num); plt.title("Numeric RGB — Blinn, two lights")
plt.axis('off'); plt.show()Expect: a crimson sphere keyed white from the upper right, filled gold from the left; ~0.5–2 s at 256² in Colab.
🐍 Cell 2 — the symbolic derivation: half-vector, three channels, symbolic Σ — and the CSE that tames it
import sympy as sp
# --- Symbols -----------------------------------------------------------------
x, y = sp.symbols('x y', real=True)
R, e = sp.symbols('R e', positive=True)
ka, kd, ks, n = sp.symbols('k_a k_d k_s n', positive=True)
L1v = sp.symbols('l_1x l_1y l_1z', real=True); L1 = sp.Matrix(L1v)
L2v = sp.symbols('l_2x l_2y l_2z', real=True); L2 = sp.Matrix(L2v)
c1v = sp.symbols('c_1r c_1g c_1b', real=True); c1 = sp.Matrix(c1v)
c2v = sp.symbols('c_2r c_2g c_2b', real=True); c2 = sp.Matrix(c2v)
av = sp.symbols('a_r a_g a_b', real=True); amb = sp.Matrix(av)
mv = sp.symbols('m_r m_g m_b', real=True); mat = sp.Matrix(mv)
# --- Shared geometry: exactly Chapter 1, unchanged ----------------------------
E = sp.Matrix([0, 0, e]); D = sp.Matrix([x, y, -e])
A = D.dot(D)
disc = sp.expand(E.dot(D)**2 - A*(e**2 - R**2))
t_hit = (-E.dot(D) - sp.sqrt(disc)) / A
P = E + t_hit * D
N = P / R
v_hat = -D / sp.sqrt(A)
nu = sp.sqrt(disc) / (R * sp.sqrt(A)) # N.v_hat, closed form (Ch.1 gift)
def light_geo(L):
w = L - P
w_hat = w / sp.sqrt(w.dot(w))
mu = N.dot(w_hat) # diffuse term, this light
phi = 2*mu*nu - w_hat.dot(v_hat) # Phong r.v, closed form (no new sqrt)
h = w_hat + v_hat
psi = N.dot(h) / sp.sqrt(h.dot(h)) # Blinn N.h (one new sqrt per light)
return mu, phi, psi
mu1, phi1, psi1 = light_geo(L1)
mu2, phi2, psi2 = light_geo(L2)
def channel(k, spec1, spec2):
return sp.Piecewise(
(ka*amb[k] + c1[k]*(kd*mat[k]*sp.Max(0, mu1) + ks*sp.Max(0, spec1)**n)
+ c2[k]*(kd*mat[k]*sp.Max(0, mu2) + ks*sp.Max(0, spec2)**n),
disc >= 0),
(0, True))
I_blinn = sp.Matrix([channel(k, psi1, psi2) for k in range(3)])
I_phong = sp.Matrix([channel(k, phi1, phi2) for k in range(3)])
print("--- swell meter (ops per channel, unshared) ---")
print("blinn:", [sp.count_ops(expr) for expr in I_blinn])
print("phong:", [sp.count_ops(expr) for expr in I_phong])
# --- The swell-tamer: common subexpression elimination ------------------------
# Three channels share ALL geometry; two lights share N, t, v_hat.
subs, reduced = sp.cse(list(I_blinn))
total = sum(sp.count_ops(rhs) for _, rhs in subs) + \
sum(sp.count_ops(expr) for expr in reduced)
print("blinn after cse:", total, "ops total (all three channels)")
print("shared subexpressions hoisted:", len(subs))Expect: naive swell roughly triples Chapter 1's — then cse collapses it to little more than one channel's cost. The sharing is not a hope; it is counted.
🐍 Cell 3 — compile the RGB image, evaluate everywhere, verify against numeric
params = (x, y, R, e, *L1v, *L2v, *c1v, *c2v, *av, *mv, ka, kd, ks, n)
f_blinn = sp.lambdify(params, list(I_blinn), modules='numpy')
f_phong = sp.lambdify(params, list(I_phong), modules='numpy')
SC2 = dict(R=1.0, e=4.0,
lights=(((5., 5., 10.), (1., 1., 1.)),
((-6., 2., 8.), (0.9, 0.65, 0.25))),
ambient=(0.15, 0.15, 0.16), material=(0.63, 0.19, 0.16),
ka=0.1, kd=0.7, ks=0.5, shin=32.0)
def grid(W, H, R=1.0, e=4.0):
half = 1.25*e*R/np.sqrt(e*e - R*R)
xs = np.linspace(-half + half/W, half - half/W, W)
ys = np.linspace( half - half/H, -half + half/H, H)
return np.meshgrid(xs, ys)
def eval_sym(f, X, Y, SC, lights=None):
(p1, col1), (p2, col2) = SC['lights'] if lights is None else lights
with np.errstate(invalid='ignore', divide='ignore'):
chans = f(X, Y, SC['R'], SC['e'], *p1, *p2, *col1, *col2,
*SC['ambient'], *SC['material'],
SC['ka'], SC['kd'], SC['ks'], SC['shin'])
return np.clip(np.nan_to_num(np.stack(chans, axis=-1)), 0, 1)
X, Y = grid(256, 256)
t0 = time.perf_counter(); img_sym = eval_sym(f_blinn, X, Y, SC2); t1 = time.perf_counter()
print(f"symbolic RGB eval: {t1 - t0:.4f} s")
img_num = render_numeric_rgb(256, 256, **SC2) # model='blinn' by default
diff = np.abs(img_num - img_sym)
print(f"max |num - sym| over ALL channels = {diff.max():.3e} (structural agreement)")
fig, axs = plt.subplots(1, 3, figsize=(12, 4))
for ax, im, ttl in zip(axs, [img_num, img_sym, diff / diff.max()],
["Numeric", "Symbolic (compiled)", "|difference| (auto-scaled)"]):
ax.imshow(im); ax.set_title(ttl); ax.axis('off')
plt.show()Expect: identical images; max difference ~1e-15 across all three channels — the Phong-identity and half-vector closed forms validated simultaneously.
🐍 Cell 4 — honest timings, RGB edition: loop vs vectorized vs compiled Blinn vs compiled Phong
def render_numeric_vec_rgb(width, height, R=1.0, e=4.0,
lights=(((5., 5., 10.), (1., 1., 1.)),
((-6., 2., 8.), (0.9, 0.65, 0.25))),
ambient=(0.15, 0.15, 0.16), material=(0.63, 0.19, 0.16),
ka=.1, kd=.7, ks=.5, shin=32., model='blinn'):
"""Chapter-2 math, numpy-vectorized: masks instead of ifs. Still numerical."""
half = 1.25*e*R/np.sqrt(e**2 - R**2)
xs = np.linspace(-half + half/width, half - half/width, width)
ys = np.linspace( half - half/height, -half + half/height, height)
X, Y = np.meshgrid(xs, ys)
A = X**2 + Y**2 + e**2
disc = e**2*R**2 - (e**2 - R**2)*(X**2 + Y**2)
hit = disc >= 0
sq = np.sqrt(np.maximum(disc, 0))
t = (e**2 - sq) / A
Px, Py, Pz = t*X, t*Y, e*(1 - t)
Nx, Ny, Nz = Px/R, Py/R, Pz/R
va = np.sqrt(A); vx, vy, vz = -X/va, -Y/va, e/va
m = np.array(material)
img = np.zeros((height, width, 3)) + ka*np.array(ambient)
for pos, color in lights:
if color is None:
continue
wx, wy, wz = pos[0]-Px, pos[1]-Py, pos[2]-Pz
wn = np.sqrt(wx**2 + wy**2 + wz**2)
mu = (Nx*wx + Ny*wy + Nz*wz) / wn
if model == 'phong':
rv = 2*mu*(Nx*vx + Ny*vy + Nz*vz) - (wx*vx + wy*vy + wz*vz)/wn
spec = np.maximum(0, rv)**shin
else:
hx, hy, hz = wx/wn + vx, wy/wn + vy, wz/wn + vz
hn = np.sqrt(hx**2 + hy**2 + hz**2)
spec = np.maximum(0, (Nx*hx + Ny*hy + Nz*hz)/hn)**shin
img += np.array(color)*(kd*np.maximum(0, mu)[..., None]*m + ks*spec[..., None])
return np.clip(np.where(hit[..., None], img, 0.0), 0, 1)
print(f"{'grid':>9} {'loop':>10} {'vectorized':>12} {'c-blinn':>10} {'c-phong':>10}")
for W in (64, 128, 256):
Xg, Yg = grid(W, W)
t0 = time.perf_counter(); render_numeric_rgb(W, W, **SC2); t1 = time.perf_counter()
render_numeric_vec_rgb(W, W, **SC2); t2 = time.perf_counter()
eval_sym(f_blinn, Xg, Yg, SC2); t3 = time.perf_counter()
eval_sym(f_phong, Xg, Yg, SC2); t4 = time.perf_counter()
print(f"{W:>4}x{W:<4} {t1-t0:>9.3f}s {t2-t1:>11.3f}s {t3-t2:>9.3f}s {t4-t3:>9.3f}s")
print("\nBlinn vs Phong compiled: nearly equal here -- both cost one norm per light.")
print("Blinn's classic win (constant h for distant light + distant viewer) is a")
print("SPECIAL CASE the CAS can prove: dh/dx = 0 under directional light, then CSE")
print("hoists it out of the per-pixel work. Certification, not folklore.")Expect: loop ≫ vectorized ≈ compiled-blinn ≈ compiled-phong. (512² is left to the patient — the loop is quadratic and unapologetic.)
🐍 Cell 5 — superpowers: substitution is a theorem, color is affine, the silhouette endures
# (a) Disabling light 2 IS the substitution c2 -> 0. Verify against omission:
img_sub = eval_sym(f_blinn, X, Y, SC2,
lights=(SC2['lights'][0], ((-6., 2., 8.), (0., 0., 0.))))
img_omit = render_numeric_rgb(256, 256, **{**SC2, 'lights': (SC2['lights'][0],)})
print("disable-by-substitution == omit-the-light:",
np.allclose(img_sub, img_omit, atol=1e-12))
# (b) The image is EXACTLY AFFINE in each light color -- a theorem, not a hope:
I_unc = [ka*amb[k] + c1[k]*(kd*mat[k]*mu1 + ks*psi1**n)
+ c2[k]*(kd*mat[k]*mu2 + ks*psi2**n) for k in range(3)]
print("d2I/dc1k^2 =", [sp.diff(I_unc[k], c1v[k], 2) for k in range(3)], " (all zero)")
print("dI_r/dc_1r =", sp.diff(I_unc[0], c1v[0]))
print(" -- sensitivity to a light's redness IS that light's illumination geometry.")
# (c) The silhouette theorem is untouched (geometry identical to Chapter 1):
half = 1.25*SC2['e']*SC2['R']/np.sqrt(SC2['e']**2 - SC2['R']**2)
rho = SC2['e']*SC2['R']/np.sqrt(SC2['e']**2 - SC2['R']**2)
fig, ax = plt.subplots(figsize=(4.5, 4.5))
ax.imshow(img_sym, extent=[-half, half, -half, half])
ax.add_patch(plt.Circle((0, 0), rho, color='red', fill=False, lw=1))
ax.set_title("RGB render, exact analytic silhouette overlaid")
plt.show()Expect: True; [0, 0, 0]; a geometric expression for ∂I/∂c; and the red circle threading the pixel boundary as if Chapter 1 never ended.
🐍 Cell 6 — exact derivatives and partial evaluation, now per channel
# (a) Exact per-channel sensitivity to the gold fill light's x position.
# Differentiate the UNCLAMPED field (the Max kinks are handled by reasoning,
# as in Chapter 1 §7).
dI_dl2x = [sp.diff(I_unc[k], L2v[0]) for k in range(3)]
print("ops per channel of dI/dl_2x:", [sp.count_ops(expr) for expr in dI_dl2x])
g = sp.lambdify(params, dI_dl2x, modules='numpy')
(p1, col1), (p2, col2) = SC2['lights']
with np.errstate(invalid='ignore', divide='ignore'):
G = g(X, Y, SC2['R'], SC2['e'], *p1, *p2, *col1, *col2,
*SC2['ambient'], *SC2['material'],
SC2['ka'], SC2['kd'], SC2['ks'], SC2['shin'])
discN = SC2['e']**2*SC2['R']**2 - (SC2['e']**2 - SC2['R']**2)*(X**2 + Y**2)
G = np.where(discN[..., None] >= 0, np.nan_to_num(np.stack(G, axis=-1)), 0.0)
Gvis = 0.5 + 0.5 * G / np.abs(G).max()
plt.figure(figsize=(4.5, 4)); plt.imshow(Gvis)
plt.title("∂I/∂l₂ₓ per channel (darker/lighter = dims/brightens)")
plt.axis('off'); plt.show()
# (b) Partial evaluation: bake the entire scene, keep only (x, y) free.
subs_scene = {R: SC2['R'], e: SC2['e'],
**{L1v[i]: p1[i] for i in range(3)},
**{L2v[i]: p2[i] for i in range(3)},
**{c1v[i]: col1[i] for i in range(3)},
**{c2v[i]: col2[i] for i in range(3)},
**{av[i]: SC2['ambient'][i] for i in range(3)},
**{mv[i]: SC2['material'][i] for i in range(3)},
ka: SC2['ka'], kd: SC2['kd'], ks: SC2['ks'], n: SC2['shin']}
I_baked = [expr.subs(subs_scene) for expr in I_blinn]
print("ops: free-params:", [sp.count_ops(expr) for expr in I_blinn],
" baked:", [sp.count_ops(expr) for expr in I_baked])
f_baked = sp.lambdify((x, y), I_baked, modules='numpy')
with np.errstate(invalid='ignore'):
img_baked = np.clip(np.nan_to_num(np.stack(f_baked(X, Y), axis=-1)), 0, 1)
print("baked matches compiled render:", np.allclose(img_baked, img_sym, atol=1e-12))Expect: a per-channel sensitivity field (strongest where the gold fill grazes), and a baked triple of expressions that matches to 1e-12. Inverse relighting = gradient descent on G.
🐍 Cell 7 — pitfalls, live: saturation is a region; means hide error; NaNs come per channel
# (a) SATURATION IS A REGION. The display clamp min(1, I) is silent in numeric
# pipelines; symbolically it is another kink, and its domain is computable.
with np.errstate(invalid='ignore', divide='ignore'):
chans = f_blinn(X, Y, SC2['R'], SC2['e'], *p1, *p2, *col1, *col2,
*SC2['ambient'], *SC2['material'],
SC2['ka'], SC2['kd'], SC2['ks'], SC2['shin'])
I_raw = np.nan_to_num(np.stack(chans, axis=-1))
sat = (I_raw > 1.0) & (discN[..., None] >= 0)
print("saturated pixels per channel:", sat.sum(axis=(0, 1)),
" (raise ks in the simulator and watch this grow)")
fig, axs = plt.subplots(1, 2, figsize=(8, 4))
axs[0].imshow(np.clip(I_raw, 0, 1)); axs[0].set_title("clipped image (Min kink applied)")
axs[1].imshow(sat.any(axis=-1), cmap='autumn')
axs[1].set_title("saturation region { I_k > 1 }")
for ax in axs: ax.axis('off')
plt.show()
# (b) MEAN-OVER-CHANNELS HIDES ERROR. Always report the per-channel max.
d = np.abs(img_num - img_sym)
print(f"max over channels: {d.max():.3e} "
f"mean-over-channels then max: {d.mean(axis=-1).max():.3e} "
f"(up to 3x smaller: error hiding)")
# (c) The Heaviside temptation, in triplicate: H(disc)*I NaNs PER CHANNEL.
I_heavi = [sp.Heaviside(disc, 0) *
(ka*amb[k] + c1[k]*(kd*mat[k]*sp.Max(0, mu1) + ks*sp.Max(0, psi1)**n)
+ c2[k]*(kd*mat[k]*sp.Max(0, mu2) + ks*sp.Max(0, psi2)**n))
for k in range(3)]
heavi = {'Heaviside': lambda u, h0: (u >= 0).astype(float)}
fH = sp.lambdify(params, I_heavi, modules=['numpy', heavi])
with np.errstate(invalid='ignore'):
bad = fH(X, Y, SC2['R'], SC2['e'], *p1, *p2, *col1, *col2,
*SC2['ambient'], *SC2['material'], SC2['ka'], SC2['kd'], SC2['ks'], SC2['shin'])
print("NaNs per channel (raw Heaviside):", [int(np.isnan(bk).sum()) for bk in bad])
I_regH = [expr.subs(sp.sqrt(disc), sp.sqrt(sp.Max(0, disc))) for expr in I_heavi]
fR = sp.lambdify(params, I_regH, modules=['numpy', heavi])
with np.errstate(invalid='ignore'):
good = fR(X, Y, SC2['R'], SC2['e'], *p1, *p2, *col1, *col2,
*SC2['ambient'], *SC2['material'], SC2['ka'], SC2['kd'], SC2['ks'], SC2['shin'])
print("NaNs after regularization:", [int(np.isnan(gk).sum()) for gk in good])
print("matches Piecewise render:",
np.allclose(np.nan_to_num(np.stack(good, axis=-1)), img_sym, atol=1e-12))Expect: a small saturation cap near the key highlight; the mean ~⅓ of the max; then tens of thousands of NaNs per channel, annihilated by regularization.
🐍 Cell 8 — the frontier escalates: shadows with two lights are independent quantifiers (teaser for Chapter 3)
# A blocker sphere B (center Cb, radius R2). A surface point P is lit by light L
# iff the open segment P -> L misses B: the quadratic |P + s(L-P) - Cb|^2 = R2^2
# has NO root s in (0,1). "No root in (0,1)" is the quantifier.
R2 = sp.symbols('R_2', positive=True)
Cb = sp.Matrix(sp.symbols('c_bx c_by c_bz', real=True))
def shadow_quadratic(L):
seg = L - P # segment direction: surface point -> light
a2 = seg.dot(seg)
b2 = (P - Cb).dot(seg)
c2 = (P - Cb).dot(P - Cb) - R2**2
return sp.expand(b2**2 - a2*c2) # <0: the line misses B entirely
sdisc1 = shadow_quadratic(L1)
sdisc2 = shadow_quadratic(L2)
print("blocker discriminant, light 1:", sp.count_ops(sdisc1), "ops")
print("blocker discriminant, light 2:", sp.count_ops(sdisc2), "ops")
# sdisc < 0 certainly means lit; sdisc >= 0 requires the root-in-(0,1) test --
# decidable per pixel (Sturm sequences / CAD), but NOT a closed form:
# the shadow boundary is a semi-algebraic set.
# With TWO lights the masks are INDEPENDENT: penumbra = mask1 XOR mask2, for free.
I_shadowy = sp.Piecewise(
(sp.Symbol('I_both'), (disc >= 0) & (sdisc1 < 0) & (sdisc2 < 0)),
(sp.Symbol('I_L1'), (disc >= 0) & (sdisc1 < 0) & (sdisc2 >= 0)),
(sp.Symbol('I_L2'), (disc >= 0) & (sdisc1 >= 0) & (sdisc2 < 0)),
(sp.Symbol('I_amb'), (disc >= 0)),
(0, True))
print(I_shadowy)
print("\nFour visibility chambers, counted by the algebra itself (naive version:")
print("the sdisc >= 0 chambers still need the root-in-(0,1) refinement).")
print("Chapter 3: make that test exact with Sturm/CAD, and watch the shadow")
print("boundaries emerge as algebraic curves in (x, y). The frontier, measured.")Expect: a four-chamber Piecewise. Discussion: independent per-light quantifiers, penumbra for free, and where the closed forms run out.
§6 · The honest ledger — Chapter-2 deltas
| Dimension | Chapter 1 | Chapter 2 |
|---|---|---|
| Channels | one expression | three expressions sharing all geometry; the sharing is explicit and countable via cse (Cell 2) |
| Lights | one point light | symbolic Σ; disable = substitution $c_2 \to 0$, verified against omission to 1e-12 (Cell 5a) |
| Color | — | image exactly affine in light colors; $\partial I/\partial c$ is the illumination geometry (Cell 5b) |
| Specular | Phong $r\cdot\hat V$, closed form $2\mu\nu - \hat w\cdot\hat V$ | Blinn $N\cdot\hat H$: one new sqrt per light; the constant-$\hat H$ win is a provable special case (Cell 4) |
| Clamp | clip at display, unremarked | $\min(1,\cdot)$ is a kink with a computable domain: the saturation region $\{I_k \ge 1\}$ (Cell 7a) |
| Error metric | max |ΔI| | max over channels; the mean hides up to 3× (Cell 7b) |
| Failure modes | 0·NaN, Max kinks, swell | the same, per channel; raw Heaviside now poisons three channels at once (Cell 7c) |
§7 · Pitfalls gallery, continued
- Saturation is a region, not an event. The clamp $\min(1, I_k)$ kinks on the boundary of $\{I_k \ge 1\}$ — a semi-algebraic set you can compute and plot. Numeric pipelines clip silently at display and lose the information.
- Mean-over-channels hides error. Report the per-channel max, always.
- Blinn moves the kinks. The $\max(0, N\cdot\hat H)$ terminator is a different curve from Phong's $\max(0, r\cdot\hat V)$ terminator; derivatives of the two models disagree along different loci. (And ungated Blinn can highlight surface facing away from the light — the classic artifact; gate with $\mu_\ell > 0$ if it bites.)
- NaNs are per-channel. One poisoned channel poisons the pixel's color; the raw Heaviside failure in RGB is magenta with structure.
- Swell triples unless you CSE. Measure with
count_opsbefore and after; the sharing is the whole game.
§8 · A brief history, continued
Jim Blinn, SIGGRAPH 1977 — "Models of Light Reflection for Computer Synthesized
Pictures" — imported the Torrance–Sparrow microfacet theory (1967) into graphics and,
almost in passing, replaced the reflection vector with the half-vector. The fixed-function
hardware of the 1990s enshrined the choice: OpenGL's lighting was Blinn's, so for a
decade the model you computed with was decided by silicon. The color side runs deeper:
CIE 1931's standard observer made RGB itself a linear algebra, and sRGB (HP/Microsoft,
1996) made it a household one — with a gamma nonlinearity this chapter deliberately
defers (we render in linear light and show it raw; Chapter 5 will face it).
Meanwhile, in the other world: common subexpression elimination is compiler folklore
formalized by Cocke (1970), and Macsyma's optimize() performed it for
FORTRAN code generation — the direct ancestor of the sympy.cse that tamed
Cell 2. The same decade gave us both halves of this chapter; the file merely reunites
them again.
§9 · Roadmap
- Ch. 3 — Shadows: per-light quantified blockers; the root-in-(0,1) test made exact via Sturm sequences / CAD; two lights give penumbra as independent masks; shadow boundaries as algebraic curves.
- Ch. 4 — Planar reflection: the reflected ray is a substituted expression; recursion becomes composition; watch the swell with three channels.
- Ch. 5 — Texture maps as symbolic functions u(x,y); the sRGB gamma question faced honestly.
- Ch. 6 — Symbolic antialiasing: integrate I(x,y) over the pixel square — the exact box filter, now over clamped, saturated, three-channel integrands.
❦
Set by the Press · single HTML5 file · the expression is the ground truth.
Section ids are stable (s1…s9, sim, lab, cell1…cell8) — cite the id when requesting revisions.
Section ids are stable (s1…s9, sim, lab, cell1…cell8) — cite the id when requesting revisions.