Chapter 3 — Shadows Are Quantifiers, Tamed
in the spirit of Utah graphics meets Macsyma (today: Python + sympy).
Kimi K3 — Assistance with Mathematical Development
☰ contents
§1 · What changes in Chapter 3
Everything from Chapter 2 stands: RGB, two lights, Blinn/Phong, the silhouette theorem for sphere A. The novelties, in escalating order:
- A second object. Sphere B (offset center, own radius, own material) — and with it visibility ordering: the nearest of two hits wins. The decision tree grows its first ordering branch.
- Cast shadows. A point is lit by light $\ell$ iff the open segment from point to light misses the other sphere. "Misses" is a statement about a quadratic having no root in $(0,1)$ — and Sturm (1829) already knew how to decide root-placement in an interval without computing any roots.
- The image becomes discontinuous in the scene parameters. Until now $\partial I/\partial\theta$ existed almost everywhere with kinks; shadows add jump discontinuities. The gradient acquires Dirac sheets living on curves (Cell 6 shows the finite-difference fringes).
- The first genuinely different algorithms. Numeric: find roots, test the interval, at runtime, per pixel, per light. Symbolic: a polynomial predicate — no square root anywhere in the test — with the occlusion order provably constant over each silhouette overlap (Theorem 3.3). The working simulator evaluates the comparison per pixel because one comparison is cheaper than the machinery to cache it; the theorem's content is that this is a choice, not a necessity.
- A new species of disagreement. With branches come branch boundaries. Numeric and symbolic now agree to ~10⁻¹⁵ inside every chamber, and may choose different chambers at isolated boundary pixels (last-ulp ordering). The honest ledger becomes two numbers: interior max |ΔI|, and a flip census. The flips trace the silhouette and shadow curves — the algebra's zero sets, made visible.
§2 · The scene (Chapter-3 assumptions)
Chapter 2's scene, plus: sphere B with center $C_b = (b_x, b_y, b_z)$, radius $R_2$, material $m_B$. Standing assumptions (stated, not smuggled): the spheres are disjoint; the camera and both lights lie outside both spheres; the blocker sits fully in front of the camera, $|e - b_z| > R_2$; no containment (neither sphere contains the other, nor any light, nor any shaded point of the other). Under these, every theorem below is exact; where a slider violates them, the simulator's diagnostics say so. Shadows are other-sphere only: a convex sphere's self-shadowing is exactly the terminator clamp $\max(0,\mu)$ — Theorem 3.2 proves the equivalence, and the boundary case $f(0)=0$ of the shadow quadratic is excluded by construction, and we say so. Hard shadows from point lights; penumbra appears only as the XOR of two independent hard masks. Soft shadows are integrals over area lights — a later chapter's problem, labeled as such.
§3 · The math, once, carefully
Two quadratics, one ordering. Sphere A as ever: $A_q = x^2+y^2+e^2$, $\Delta_A = e^2R^2 - (e^2-R^2)(x^2+y^2)$, $t_A = (e^2-\sqrt{\Delta_A})/A_q$. Sphere B, with $D = (x,y,-e)$ and $E = (0,0,e)$:
$$b_B = (E-C_b)\cdot D, \qquad c_B = |E-C_b|^2 - R_2^2, \qquad \Delta_B = b_B^2 - A_q\,c_B, \qquad t_B = \frac{-b_B - \sqrt{\Delta_B}}{A_q},$$ $$\text{on}_A = \big(\Delta_A \ge 0\big) \land \big(\Delta_B \lt 0 \;\lor\; t_A \le t_B\big), \qquad \text{on}_B = \big(\Delta_B \ge 0\big) \land \lnot\,\text{on}_A .$$The gifts generalize verbatim: for any sphere, the near root gives $-b - tA = \sqrt{\Delta}$, so $\nu = N\cdot\hat V = \sqrt{\Delta}/(R\sqrt A)$ holds for the blocker too — the closed form did not know the sphere was centered.
The shadow quadratic. Point $P$ (on either sphere), light $L$, blocker sphere $(C_b, R_b)$. The segment $P + s(L-P)$, $s \in (0,1)$, meets the blocker iff $f(s) = a_2 s^2 + 2 b_2 s + c_2 = 0$ has a root in the open interval, where
$$a_2 = (L-P)\cdot(L-P), \qquad b_2 = (P - C_b)\cdot(L-P), \qquad c_2 = |P - C_b|^2 - R_b^2, \qquad \Delta_s = b_2^2 - a_2 c_2 .$$Theorem 3.1 (The Shadow Predicate — Sturm for a Segment). Under non-containment ($c_2 \ge 0$ and $f(1) \ge 0$: $P$ and $L$ outside or on the blocker, and $a_2 \gt 0$), $$\text{shadowed} \iff \exists\, s \in (0,1):\; f(s) = 0 \iff \Delta_s \ge 0 \;\land\; b_2 \lt 0 \;\land\; a_2 + b_2 \gt 0 .$$ Proof. $a_2 \gt 0$, so $f$ opens upward with vertex at $s_v = -b_2/a_2$, where $f(s_v) = -\Delta_s/a_2$. If $\Delta_s \ge 0$ and $0 \lt s_v \lt 1$ (i.e. $b_2 \lt 0$, $a_2 + b_2 \gt 0$), then $f(s_v) \le 0$ while $f(0) \ge 0$, so a root lies in $(0, s_v] \subset (0,1)$. If $s_v \le 0$, $f$ is increasing on $(0,1)$, so $f \ge f(0) \ge 0$; if $s_v \ge 1$, $f$ is decreasing, $f \ge f(1) \ge 0$. If $\Delta_s \lt 0$, no real roots at all. ∎ The test is polynomial — three signs, no square roots. Shading needs √Δ; shadows are algebraic sign logic. This is Sturm's 1829 theorem in its smallest nontrivial instance, and the shadow region is a semi-algebraic set whose boundary is an algebraic curve. Cell 2 verifies the lemma against explicit root computation on 200,000 random configurations; the simulator's self-test re-verifies 200 more at every page load.
Theorem 3.2 (The Self-Shadow Was Already There). Run the same test against one's own sphere: $c_2 = |P|^2 - R^2 = 0$, so the roots are $s = 0$ and $s = -2b_2/a_2$ with $b_2 = P\cdot(L-P)$. Since $\operatorname{sign}(b_2) = \operatorname{sign}(\mu)$ and $|L| \gt R$ gives $-2b_2 \lt a_2$: $$\text{self-shadowed} \iff \mu \lt 0 \iff \max(0,\mu) = 0.$$ The diffuse clamp is the occlusion test against one's own convex body — the terminator was a shadow boundary all along. Cell 5 verifies the equivalence to machine precision. (This is why the simulator excludes self-blockers by construction: the work was already done in Chapter 1's clamp.)
Theorem 3.3 (The Constant Order Theorem). For two disjoint spheres under one pinhole camera, the nearer-surface relation is constant over the entire overlap of their silhouettes. Proof. Each silhouette region is a disk (convex); their overlap is convex, hence connected. On it, $t_A - t_B$ is continuous and never zero (zero would put one point on both spheres, contradicting disjointness). A continuous nonvanishing real function on a connected set has constant sign. ∎ The numeric renderer compares $t$'s per pixel; the symbolic renderer may decide the sign once per frame. The working simulator evaluates the comparison per pixel — one comparison is cheaper than the cache — and the theorem's real content is that the ordering locus $t_A = t_B$ is a single algebraic curve, not a per-pixel accident. Cell 5(a) measures the agreement either way: zero flips across the overlap.
Theorem 3.4 (The Off-Axis Silhouette Is a Conic). $\Delta_B(x,y)$ is quadratic in $(x,y)$, so the blocker's silhouette is a conic — and its conic discriminant factors: $$Q_{xx}Q_{yy} - Q_{xy}^2 \;=\; \big(|E-C_b|^2 - R_2^2\big)\,\big((e-b_z)^2 - R_2^2\big) \gt 0,$$ an ellipse under the standing assumptions (the camera is outside the blocker, which is fully in front of it). Cell 5(d) has the CAS perform the factorization; the simulator's silhouette overlay draws the conic live.
The chamber system. With lit-masks $\lambda_\ell = 1 - \text{shadowed}_\ell$ (per surface, per light, blocker = the other sphere), the image is
$$I_k = \text{on}_A \cdot I^A_k + \text{on}_B \cdot I^B_k, \qquad I^\bullet_k = k_a a_k + \sum_{\ell=1}^{2} \lambda^\bullet_\ell\, c_{\ell k}\Big(k_d\,m^\bullet_k \max(0,\mu^\bullet_\ell) + k_s\,\text{spec}^\bullet_\ell\Big),$$a disjoint union of chambers: 2 surfaces × 2² light-masks + background = 9 regions, whose boundaries are semi-algebraic curves in $(x,y)$. Two lights' masks are independent, so penumbra is the XOR chamber — visible for free, before any area-light theory.
Three gifts of Chapter 3:
- The shadow test is sqrt-free. Numeric pipelines compute both roots and interval-test them; the predicate replaces that with three polynomial signs. The simulator runs both formulations — explicit roots on the numeric side, sign logic on the symbolic side — and their agreement is the lemma, demonstrated.
- Branch flips are the theory's fingerprint, not its failure. Inside each open chamber the two pipelines agree to float-ordering noise; along boundary curves they may choose different chambers in the last ulp. The difference canvas therefore draws the silhouettes and shadow boundaries: the zero sets of $\Delta_A$, $\Delta_B$, $t_A - t_B$, and each $\Delta_s$, without any edge detector.
- No epsilon. Classic shadow rays need the ugly $\epsilon$-offset hack or they acne. The exact sign tests use strict inequalities against computed quantities; there is no fudge factor to tune, and no acne to suppress.
Counterpoint. Totality has a price — the symbolic form evaluates both spheres and both shadow quadratics at every pixel (no early exit), a constant-factor cost, still vectorized-class (Cell 4). And the predicate is quantifier elimination for one quadratic: Tarski (1951) guarantees QE for all polynomial scenes, Collins' CAD (1975) makes it mechanical, and both are hyper-exponential in general. Quadrics are the sweet spot where the elimination fits on one line. Beyond them, the frontier from Chapter 1 stands. A further note on compilation: the monolithic expression is no longer the right target — compile the chamber decomposition, not the monolith: lambdify the masks and the ungated per-light terms, and assemble in numpy. This is CSE-by-hand at the level of logic, and it is what every real engine does without saying so.
§4 · Live simulator
Two spheres, two lights, shadows both ways. Numeric shoots explicit-root shadow rays; symbolic evaluates the Shadow Predicate's sign logic. The ledger reports interior max |ΔI| and the flip census; the chambers view paints the nine regions, the difference panel draws the boundary curves in white, and the silhouette overlay draws Theorem 3.4's conic.
Try this: (1) Uncheck cast shadows — the shadow lifts; agreement is untouched: the symbolic path just substituted mask → 1. (2) Uncheck light 2 — the gold fill dies and the shadowed patch plunges to ambient: penumbra → umbra by substitution c₂ → 0. (3) Drag the blocker's b_x: the shadow slides across the sphere while the ellipse overlay migrates — the conic of Theorem 3.4, tracked live. (4) Switch the view to visibility chambers: nine regions, flips in white — the algebraic boundaries with no edge detector. (5) Heaviside multiply (raw): the NaNs come from the other chamber — a pixel on sphere A still evaluates I_B, whose √Δ_B is NaN where B was missed, and 0·NaN = NaN poisons across chambers. The regularized mode restores the identity.
§5 · The Colab laboratory (Python / sympy)
Paste the cells in order, run in order; three run paths per cell —
📋 Copy code, ⬇ notebook (File → Upload notebook in Colab), or paste into
colab.new.
🐍 Cell 1 — the numerical baseline: two spheres, explicit-root shadow rays
import numpy as np
import matplotlib.pyplot as plt
import time
# Scene defaults match the live simulator exactly.
SC3 = dict(R=1.0, e=4.0,
Cb=(0.5, 0.5, 1.3), R2=0.32,
lights=(((5., 5., 10.), (1., 1., 1.)),
((-6., 2., 8.), (0.9, 0.65, 0.25))),
ambient=(0.15, 0.15, 0.16),
matA=(0.63, 0.19, 0.16), matB=(0.33, 0.44, 0.62),
ka=0.1, kd=0.7, ks=0.5, shin=32.0, model='blinn')
def render_numeric_2s(width, height, SC=SC3):
"""Chapter 3 numerical renderer: nearest of two spheres; per light, an
EXPLICIT-ROOT shadow ray against the other sphere. Branches at runtime."""
R, e = SC['R'], SC['e']
Cb = np.array(SC['Cb']); R2 = SC['R2']
img = np.zeros((height, width, 3))
code = np.zeros((height, width), dtype=np.uint8) # chamber codes
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])
Aq = float(D @ D)
# ---- sphere A ----
discA = e*e*R*R - (e*e - R*R)*(x*x + y*y)
hitA = discA >= 0.0
tA = (e*e - np.sqrt(max(discA, 0.0))) / Aq if hitA else np.inf
# ---- sphere B ----
EC = E - Cb
bB = float(EC @ D); cB = float(EC @ EC) - R2*R2
discB = bB*bB - Aq*cB
hitB = discB >= 0.0
tB = (-bB - np.sqrt(max(discB, 0.0))) / Aq if hitB else np.inf
# ---- nearest-surface branch ----
if hitA and tA <= tB:
Cen, Rr, mat, onA = np.zeros(3), R, np.array(SC['matA']), True
t = tA
elif hitB:
Cen, Rr, mat, onA = Cb, R2, np.array(SC['matB']), False
t = tB
else:
continue # background
P = E + t * D
N = (P - Cen) / Rr
BCen, BRr = (Cb, R2) if onA else (np.zeros(3), R) # blocker = other
v_hat = (E - P) / np.linalg.norm(E - P)
col = SC['ka'] * np.array(SC['ambient'])
cd = 1 if onA else 2
for li, (pos, color) in enumerate(SC['lights']):
if color is None:
continue
Lp = np.array(pos); seg = Lp - P
# ---- shadow ray: EXPLICIT roots, interval test (classic) ----
a2 = float(seg @ seg)
b2 = float((P - BCen) @ seg)
c2 = float((P - BCen) @ (P - BCen)) - BRr*BRr
dsc = b2*b2 - a2*c2
shadowed = False
if dsc >= 0.0:
sq = np.sqrt(dsc)
s1 = (-b2 - sq) / a2; s2 = (-b2 + sq) / a2
if (0.0 < s1 < 1.0) or (0.0 < s2 < 1.0):
shadowed = True
if shadowed:
continue # this light contributes nothing
w_hat = seg / np.sqrt(a2)
mu = float(N @ w_hat)
if SC['model'] == 'phong':
r = 2.0 * mu * N - w_hat
spec = max(0.0, float(r @ v_hat)) ** SC['shin']
else:
h = w_hat + v_hat; h = h / np.linalg.norm(h)
spec = max(0.0, float(N @ h)) ** SC['shin']
col += np.array(color) * (SC['kd'] * mat * max(0.0, mu)
+ SC['ks'] * spec)
cd |= 4 << li # lit bit for this light
img[j, i] = col
code[j, i] = cd
return np.clip(img, 0.0, 1.0), code
t0 = time.perf_counter()
img_num, code_num = render_numeric_2s(256, 256)
t1 = time.perf_counter()
print(f"numeric render (two spheres, shadow rays): {t1 - t0:.3f} s")
plt.figure(figsize=(4, 4)); plt.imshow(img_num)
plt.title("Numeric: two spheres, two lights, shadows"); plt.axis('off'); plt.show()🐍 Cell 2 — the Shadow Predicate: quantifier elimination for a quadratic, verified 200,000×
import sympy as sp
# --- The symbolic statement ----------------------------------------------------
s, a2, b2, c2 = sp.symbols('s a_2 b_2 c_2', real=True)
f = a2*s**2 + 2*b2*s + c2
print("shadow quadratic: f(s) =", f, " on the OPEN segment s in (0, 1)")
print("roots: s =", *sp.solve(sp.Eq(f, 0), s))
print("""
Shadow Predicate (Thm 3.1; non-containment: c2 >= 0 and f(1) >= 0):
exists s in (0,1) with f(s) = 0
<=> disc_s >= 0 AND b2 < 0 AND a2 + b2 > 0
The existential over a continuum collapses to three polynomial signs.
No square root appears: shading needs sqrt(disc); shadows do not.
""")
# --- Numerical verification against EXPLICIT root computation ------------------
rng = np.random.default_rng(7)
M = 200_000
P = rng.normal(size=(M, 3)) * 2.0 # surface points
L = rng.normal(size=(M, 3)) * 4.0 + [0, 0, 6.0] # lights
Cb = rng.normal(size=(M, 3)) * 1.5 # blocker centers
Rb = rng.uniform(0.2, 1.5, size=M) # blocker radii
seg = L - P
A2 = np.einsum('ij,ij->i', seg, seg)
B2 = np.einsum('ij,ij->i', P - Cb, seg)
C2 = np.einsum('ij,ij->i', P - Cb, P - Cb) - Rb**2
DSC = B2**2 - A2*C2
keep = (C2 >= 0) & (A2 + 2*B2 + C2 >= 0) # non-containment of P and L
keep &= np.abs(DSC) > 1e-9 # skip near-tangent degenerates
# formulation 1 (classic numeric): explicit roots, interval membership
sq = np.sqrt(np.maximum(DSC, 0))
s1 = (-B2 - sq) / A2; s2 = (-B2 + sq) / A2
shadow_roots = (DSC >= 0) & (((s1 > 0) & (s1 < 1)) | ((s2 > 0) & (s2 < 1)))
# near-boundary samples where the explicit test itself is ulp-sensitive:
marginal = (np.abs(s1) < 1e-7) | (np.abs(s1 - 1) < 1e-7) | \
(np.abs(s2) < 1e-7) | (np.abs(s2 - 1) < 1e-7)
keep &= ~marginal
# formulation 2 (symbolic): the Shadow Predicate's sign logic
shadow_signs = (DSC >= 0) & (B2 < 0) & (A2 + B2 > 0)
agree = (shadow_roots == shadow_signs)[keep]
print(f"verified on {keep.sum()} random configurations "
f"(skipped {M - keep.sum()} degenerate/marginal)")
print(f"agreements: {agree.sum()} / {keep.sum()} "
f"-> predicate holds: {agree.all()} ✓")🐍 Cell 3 — compile the chambered image; interior agreement + the flip census
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 render_symbolic_2s(X, Y, SC=SC3, regularized=True):
"""Chapter 3 symbolic renderer: chamber masks x total shading expressions.
Shadows via the Shadow Predicate's SIGN LOGIC (no square roots)."""
R, e = SC['R'], SC['e']
bx, by, bz = SC['Cb']; R2 = SC['R2']
sq = (lambda v: np.sqrt(np.maximum(v, 0))) if regularized else np.sqrt
Aq = X**2 + Y**2 + e**2
# sphere A
discA = e*e*R*R - (e*e - R*R)*(X**2 + Y**2)
tA = (e*e - sq(discA)) / Aq
# sphere B
bB = -bx*X - by*Y - (e - bz)*e
cB = bx**2 + by**2 + (e - bz)**2 - R2**2
discB = bB**2 - Aq*cB
tB = (-bB - sq(discB)) / Aq
# chambers
onA = (discA >= 0) & ((discB < 0) | (tA <= tB))
onB = (discB >= 0) & ~onA
code = onA.astype(np.uint8) + 2*onB.astype(np.uint8)
v = np.sqrt(Aq); vx, vy, vz = -X/v, -Y/v, e/v # pixel-only view unit vector
img = np.zeros(X.shape + (3,))
amb = SC['ka'] * np.array(SC['ambient'])
for surf, on, t, Rr, Cen, mat in (('A', onA, tA, R, (0., 0., 0.), SC['matA']),
('B', onB, tB, R2, SC['Cb'], SC['matB'])):
Px, Py, Pz = t*X, t*Y, e*(1 - t)
Nx, Ny, Nz = (Px - Cen[0])/Rr, (Py - Cen[1])/Rr, (Pz - Cen[2])/Rr
nu = sq(discA if surf == 'A' else discB) / (Rr*np.sqrt(Aq))
BCen, BRr = (SC['Cb'], R2) if surf == 'A' else ((0., 0., 0.), R)
col = np.zeros(X.shape + (3,)) + amb
for li, (pos, color) in enumerate(SC['lights']):
wx, wy, wz = pos[0]-Px, pos[1]-Py, pos[2]-Pz
a2 = wx*wx + wy*wy + wz*wz
b2 = (Px-BCen[0])*wx + (Py-BCen[1])*wy + (Pz-BCen[2])*wz
c2 = (Px-BCen[0])**2 + (Py-BCen[1])**2 + (Pz-BCen[2])**2 - BRr**2
dsc = b2*b2 - a2*c2
lit = ~((dsc >= 0) & (b2 < 0) & (a2 + b2 > 0)) # Shadow Predicate
code = code | ((on & lit).astype(np.uint8) << (2 + li))
wn = np.sqrt(a2)
mu = (Nx*wx + Ny*wy + Nz*wz) / wn
if SC['model'] == 'phong':
wv = (wx*vx + wy*vy + wz*vz) / wn
spec = np.maximum(0, 2*mu*nu - wv) ** SC['shin']
else:
hx, hy, hz = wx/wn + vx, wy/wn + vy, wz/wn + vz
hn = np.sqrt(hx*hx + hy*hy + hz*hz)
spec = np.maximum(0, (Nx*hx + Ny*hy + Nz*hz)/hn) ** SC['shin']
add = SC['kd']*np.array(mat)*np.maximum(0, mu)[..., None] \
+ SC['ks']*spec[..., None]
col += np.array(color) * add * lit[..., None]
img += col * on[..., None] # chamber select by multiply
return np.clip(np.nan_to_num(img), 0, 1), code
X, Y = grid(256, 256)
t0 = time.perf_counter()
img_sym, code_sym = render_symbolic_2s(X, Y)
t1 = time.perf_counter()
print(f"symbolic render: {t1 - t0:.4f} s")
same = code_num == code_sym
flips = ~same
d = np.abs(img_num - img_sym)
print(f"interior max |num - sym| = {d[same].max():.3e} "
f"(over {same.sum()} matching-chamber pixels)")
print(f"branch flips (chamber disagreement): {flips.sum()} pixels "
f"-- the boundary curves themselves")
fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(img_num); axs[0].set_title("Numeric")
axs[1].imshow(img_sym); axs[1].set_title("Symbolic (chambers)")
ov = img_sym.copy(); ov[flips] = [1, 1, 1]
axs[2].imshow(ov); axs[2].set_title("flips in white: the algebraic boundaries")
for ax in axs: ax.axis('off')
plt.show()🐍 Cell 4 — honest timings: totality has a price, and it is a constant factor
def render_numeric_vec_2s(width, height, SC=SC3):
"""Chapter-3 math, numpy-vectorized: masks instead of ifs. Still numerical."""
R, e = SC['R'], SC['e']
bx, by, bz = SC['Cb']; R2 = SC['R2']
half = 1.25*e*R/np.sqrt(e*e - R*R)
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)
Aq = X**2 + Y**2 + e**2
discA = e*e*R*R - (e*e - R*R)*(X**2 + Y**2)
tA = (e*e - np.sqrt(np.maximum(discA, 0))) / Aq
bB = -bx*X - by*Y - (e - bz)*e
cB = bx**2 + by**2 + (e - bz)**2 - R2**2
discB = bB**2 - Aq*cB
tB = (-bB - np.sqrt(np.maximum(discB, 0))) / Aq
onA = (discA >= 0) & ((discB < 0) | (tA <= tB))
onB = (discB >= 0) & ~onA
v = np.sqrt(Aq); vx, vy, vz = -X/v, -Y/v, e/v
img = np.zeros(X.shape + (3,))
amb = SC['ka'] * np.array(SC['ambient'])
for on, t, Rr, Cen, mat in ((onA, tA, R, (0., 0., 0.), SC['matA']),
(onB, tB, R2, SC['Cb'], SC['matB'])):
Px, Py, Pz = t*X, t*Y, e*(1 - t)
Nx, Ny, Nz = (Px - Cen[0])/Rr, (Py - Cen[1])/Rr, (Pz - Cen[2])/Rr
BCen, BRr = (SC['Cb'], R2) if Rr == R else ((0., 0., 0.), R)
col = np.zeros(X.shape + (3,)) + amb
for pos, color in SC['lights']:
wx, wy, wz = pos[0]-Px, pos[1]-Py, pos[2]-Pz
a2 = wx*wx + wy*wy + wz*wz
b2 = (Px-BCen[0])*wx + (Py-BCen[1])*wy + (Pz-BCen[2])*wz
c2 = (Px-BCen[0])**2 + (Py-BCen[1])**2 + (Pz-BCen[2])**2 - BRr**2
dsc = b2*b2 - a2*c2
lit = ~((dsc >= 0) & (b2 < 0) & (a2 + b2 > 0))
wn = np.sqrt(a2)
mu = (Nx*wx + Ny*wy + Nz*wz) / wn
hx, hy, hz = wx/wn + vx, wy/wn + vy, wz/wn + vz
hn = np.sqrt(hx*hx + hy*hy + hz*hz)
spec = np.maximum(0, (Nx*hx + Ny*hy + Nz*hz)/hn) ** SC['shin']
col += np.array(color) * lit[..., None] * (
SC['kd']*np.array(mat)*np.maximum(0, mu)[..., None]
+ SC['ks']*spec[..., None])
img += col * on[..., None]
return np.clip(img, 0, 1)
print(f"{'grid':>9} {'numeric loop':>14} {'numeric vec':>13} {'symbolic chambers':>19}")
for W in (64, 128, 256):
Xg, Yg = grid(W, W)
t0 = time.perf_counter(); render_numeric_2s(W, W); t1 = time.perf_counter()
render_numeric_vec_2s(W, W); t2 = time.perf_counter()
render_symbolic_2s(Xg, Yg); t3 = time.perf_counter()
print(f"{W:>4}x{W:<4} {t1-t0:>13.3f}s {t2-t1:>12.3f}s {t3-t2:>18.3f}s")
print("\nThe symbolic form evaluates BOTH spheres and BOTH shadow quadratics at")
print("every pixel: no early exit is possible inside a total expression. That is")
print("roughly 2x the one-sphere work -- a constant factor, not a complexity class.")
print("Early exit is the numeric renderer's privilege; totality is the symbolic")
print("object's price of admission, and it buys inspectability.")🐍 Cell 5 — superpowers: the chamber map, penumbra as XOR, boundaries as curves, the self-shadow equivalence, the conic factorization
# (a) THE CHAMBER MAP: 2 surfaces x 2^2 light-masks + background = 9 regions.
PAL = {0:(0,0,0), # background
1+4+8:(217,77,66), 1+4:(158,51,46), 1+8:(128,41,36), 1:(77,26,26), # A
2+4+8:(115,153,217), 2+4:(84,115,166), 2+8:(69,92,133), 2:(38,51,77)} # B
cmap_img = np.zeros(code_sym.shape + (3,), dtype=np.uint8)
for k, rgb in PAL.items():
cmap_img[code_sym == k] = rgb
fig, axs = plt.subplots(1, 2, figsize=(10, 4.4))
axs[0].imshow(cmap_img); axs[0].set_title("chamber map: 9 semi-algebraic regions")
axs[0].axis('off')
# (b) PENUMBRA IS XOR: exactly one light shadowed -> mixed lighting.
m1 = (code_sym >> 2) & 1 # lit by light 1
m2 = (code_sym >> 3) & 1 # lit by light 2
surf = (code_sym & 3) > 0
penumbra = surf & (m1 ^ m2 == 1)
print(f"pixels: surface {surf.sum()}, penumbra (XOR) {penumbra.sum()}, "
f"umbra (both shadowed) {(surf & (m1 + m2 == 0)).sum()}")
pen_img = np.zeros(code_sym.shape + (3,), dtype=np.uint8)
pen_img[surf & (m1 & m2 == 1)] = [217, 77, 66] # fully lit
pen_img[penumbra] = [184, 134, 11] # exactly one light: gold
pen_img[surf & (m1 + m2 == 0)] = [31, 58, 95] # umbra: navy
axs[1].imshow(pen_img); axs[1].set_title("lit / penumbra (XOR) / umbra")
axs[1].axis('off'); plt.show()
# (c) THE SHADOW BOUNDARY IS A CURVE: trace the zero set of the lit-mask
# transition of the compiled symbolic expression itself.
half = 1.25*SC3['e']*SC3['R']/np.sqrt(SC3['e']**2 - SC3['R']**2)
fig, ax = plt.subplots(figsize=(4.6, 4.4))
ax.imshow(img_sym, extent=[-half, half, -half, half])
ax.contour(X, Y, (m1 ^ 1).astype(float), levels=[0.5],
colors='white', linewidths=0.8)
ax.set_title("light-1 shadow boundary: zero set of the mask")
ax.axis('off'); plt.show()
# (d) THE SELF-SHADOW EQUIVALENCE (Thm 3.2) and the CONIC FACTORIZATION (Thm 3.4):
# the terminator clamp IS the self-occlusion test; the CAS factors the
# blocker's silhouette discriminant.
import sympy as sp
print("\nThm 3.2 check: self-shadow <=> mu < 0, on 200k sphere points:")
rng = np.random.default_rng(3)
Pn = rng.normal(size=(200_000, 3)); Pn /= np.linalg.norm(Pn, axis=1, keepdims=True)
Lv = np.array(SC3['lights'][0][0])
segv = Lv - Pn
a2v = np.einsum('ij,ij->i', segv, segv)
b2v = np.einsum('ij,ij->i', Pn, segv) # c2 = 0 exactly on the sphere
# exact-root self test: roots 0 and -2*b2/a2; blocked iff the second lies in (0,1)
s_other = -2*b2v/a2v
blocked_self = (s_other > 0) & (s_other < 1)
mu_v = np.einsum('ij,ij->i', Pn, segv) / np.sqrt(a2v)
equiv = (blocked_self == (mu_v < 0))
print(f" equivalence holds: {equiv.sum()}/{len(equiv)} -> {equiv.all()} ✓")
print("\nThm 3.4: the blocker's silhouette conic, factored by the CAS:")
xx, yy = sp.symbols('x y', real=True)
ee, RR = SC3['e'], SC3['R']
bx_, by_, bz_, R2_ = SC3['Cb'][0], SC3['Cb'][1], SC3['Cb'][2], SC3['R2']
Aq_s = xx**2 + yy**2 + ee**2
bB_s = -bx_*xx - by_*yy - (ee - bz_)*ee
cB_s = bx_**2 + by_**2 + (ee - bz_)**2 - R2_**2
discB_s = sp.expand(bB_s**2 - Aq_s*cB_s)
Q = sp.Poly(discB_s, xx, yy)
Qxx, Qyy = Q.coeff_monomial(xx**2), Q.coeff_monomial(yy**2)
Qxy = Q.coeff_monomial(xx*yy)/2
detQ = sp.simplify(Qxx*Qyy - Qxy**2)
print(f" Qxx*Qyy - Qxy^2 = {detQ}")
factored = (bx_**2 + by_**2 + (ee - bz_)**2 - R2_**2)*((ee - bz_)**2 - R2_**2)
print(f" factored form = {sp.nsimplify(detQ)} vs product {factored}")
print(f" equal: {abs(float(detQ) - factored) < 1e-9}, positive (ellipse): {float(detQ) > 0}")🐍 Cell 6 — derivatives in the presence of visibility: smooth part + boundary delta
# How does the image respond to moving the blocker (sphere B) in x?
# INSIDE a chamber, the shadow mask is constant: dI/d(b_x) = (smooth shading
# change on sphere B only; sphere A's shading does not contain b_x at all).
# AT the shadow boundary, the mask jumps: the derivative is a Dirac delta
# living on the boundary curve -- the visibility gradient of Chapter 1's ledger.
# Symbolic rendering makes the split EXPLICIT: differentiate the smooth field,
# account for the boundary separately.
# (a) Smooth part: sphere B's shading as a function of its center -- exact.
# Here, evaluated numerically by the compiled form at two nearby centers.
dBx = 1e-4
SC_lo = {**SC3, 'Cb': (SC3['Cb'][0] - dBx, SC3['Cb'][1], SC3['Cb'][2])}
SC_hi = {**SC3, 'Cb': (SC3['Cb'][0] + dBx, SC3['Cb'][1], SC3['Cb'][2])}
img_lo, code_lo = render_symbolic_2s(X, Y, SC_lo)
img_hi, code_hi = render_symbolic_2s(X, Y, SC_hi)
dI = (img_hi - img_lo) / (2*dBx)
Gvis = 0.5 + 0.5 * dI / (np.abs(dI).max() + 1e-30)
plt.figure(figsize=(4.5, 4)); plt.imshow(Gvis)
plt.title("∂I/∂b_x: smooth on B's face, spiky on the shadow boundary")
plt.axis('off'); plt.colorbar(); plt.show()
# (b) The boundary delta: where did the CHAMBER assignment move?
moved = code_lo != code_hi
print(f"chamber-boundary pixels swept by a ±{dBx} blocker move: {moved.sum()}")
print("Those pixels carry the Dirac part. Smooth field: differentiable by sympy.")
print("Boundary: a moving algebraic curve. Both are exact; they are different")
print("kinds of object, and the symbolic form keeps them apart.")
ov = img_sym.copy(); ov[moved] = [1, 1, 1]
plt.figure(figsize=(4.5, 4)); plt.imshow(ov)
plt.title("the swept shadow/silhouette boundary (white)")
plt.axis('off'); plt.show()🐍 Cell 7 — pitfalls, live: cross-chamber NaNs, and the two-number reporting rule
# (a) RAW Heaviside assembly: I = onA*I_A + onB*I_B with RAW sqrt.
# Where sphere B is missed, sqrt(disc_B) is NaN; I_B is NaN; and the
# selector 0 * NaN = NaN poisons pixels that belong to SPHERE A.
img_raw, code_raw = render_symbolic_2s(X, Y, SC3, regularized=False)
nan_px = np.isnan(img_raw).any(axis=-1)
on_any_surface = (code_sym & 3) > 0
print(f"NaN pixels, raw assembly: {nan_px.sum()} of {nan_px.size}")
print(f"NaN pixels that lie ON a sphere (cross-chamber poisoning): "
f"{(nan_px & on_any_surface).sum()}")
print("Chapter 1's background NaNs have learned to walk: they now cross chambers.")
fig, axs = plt.subplots(1, 2, figsize=(9, 4))
axs[0].imshow(np.nan_to_num(img_raw)); axs[0].set_title("raw Heaviside: NaN census")
axs[0].axis('off')
axs[1].imshow(nan_px, cmap='autumn'); axs[1].set_title("NaN support")
axs[1].axis('off'); plt.show()
# (b) Regularization annihilates them, and the two-number rule reports honestly:
img_reg, code_reg = render_symbolic_2s(X, Y, SC3, regularized=True)
print(f"NaNs after regularization: {np.isnan(img_reg).sum()}")
same = code_num == code_reg
d = np.abs(img_num - img_reg)
print(f"interior max |num - sym|: {d[same].max():.3e} "
f"flip pixels: {(~same).sum()}")
print("Rule: with branches, NEVER report a single max. The max lives on a")
print("boundary pixel by construction; the interior max is the agreement;")
print("the flip census is the boundary. Two numbers, or you are lying.")🐍 Cell 8 — the frontier, updated: k blockers are a product; reflection is a composition (teaser for Ch. 4)
# (a) THE FOR-ALL IS A PRODUCT. "Lit" = "no blocker occludes": with k blocker
# spheres, the lit mask is the PRODUCT of k per-blocker lit masks. The
# universal quantifier over a finite set is multiplication. Swell, measured:
import sympy as sp
xx, yy = sp.symbols('x y', real=True)
R_, e_ = sp.symbols('R e', positive=True)
Aq_ = xx**2 + yy**2 + e_**2
discA_ = e_**2*R_**2 - (e_**2 - R_**2)*(xx**2 + yy**2)
tA_ = (e_**2 - sp.sqrt(discA_)) / Aq_
P_ = sp.Matrix([tA_*xx, tA_*yy, e_*(1 - tA_)])
L_ = sp.Matrix(sp.symbols('l_x l_y l_z', real=True))
def shadow_ops(cx, cy, cz, rb):
Cb_ = sp.Matrix([cx, cy, cz])
seg_ = L_ - P_
a2_ = seg_.dot(seg_)
b2_ = (P_ - Cb_).dot(seg_)
c2_ = (P_ - Cb_).dot(P_ - Cb_) - rb**2
return sp.expand(b2_**2 - a2_*c2_)
d1 = shadow_ops(*sp.symbols('b_1x b_1y b_1z', real=True), sp.Symbol('r_1'))
print(f"one blocker's shadow discriminant: {sp.count_ops(d1)} ops")
print("k blockers: k such masks ANDed -- linear growth in k, exact throughout.")
print("The quantifier over finitely many quadric blockers is tamed.")
# (b) REFLECTION IS COMPOSITION -- and composition is swell.
# Reflect the view direction about the sphere-A normal, then intersect
# sphere B: substitute P into B's quadratic. Measure, don't trust.
N_ = P_ / R_
D_ = sp.Matrix([xx, yy, -e_])
rdir = D_ - 2*D_.dot(N_)*N_ # reflected direction (unnormalized)
print(f"\nreflected ray direction: {sp.count_ops(rdir)} ops per component")
Cb_ = sp.Matrix(sp.symbols('b_x b_y b_z', real=True))
R2_ = sp.Symbol('R_2', positive=True)
Er = P_ - Cb_
ar = rdir.dot(rdir)
br = Er.dot(rdir)
cr = Er.dot(Er) - R2_**2
disc_refl = sp.expand(br**2 - ar*cr)
print(f"reflected-ray discriminant against sphere B: {sp.count_ops(disc_refl)} ops")
print("One bounce, one substitution: the discriminant is already a monster.")
print("Chapter 4 will live here: reflection = expression composition, and the")
print("question is not 'can we' but 'how much swell can CSE claw back'.")§6 · The honest ledger — Chapter-3 deltas
| Dimension | Chapter 2 | Chapter 3 |
|---|---|---|
| Objects | one sphere | two spheres; nearest-hit ordering is a chamber condition on $t_A - t_B$, constant over each overlap (Thm 3.3) |
| Shadows | — | per light, per surface: the Shadow Predicate — three sign tests, no square root; numeric keeps explicit roots, and the two formulations cross-validate |
| Self-shadow | the terminator clamp | …is the self-occlusion test, proven (Thm 3.2); excluded by construction, confessed |
| Silhouette | the exact circle ρ | + the blocker's ellipse, with factoring discriminant (Thm 3.4), drawn by the overlay |
| Quantifiers | — | ∃-over-(0,1) eliminated by hand; ∀-over-blockers is a product of masks (Cell 8) |
| Agreement | max |ΔI| ~1e-15 | two numbers: interior max |ΔI| ~1e-15 and a flip census on boundary curves (Cells 3, 7) |
| Boundaries | one silhouette | silhouettes, the ordering curve $t_A = t_B$, and each light's shadow curve — all algebraic, all drawn by the flip census (Cell 5) |
| Epsilon hacks | none needed yet | none needed still — exact strict sign tests replace the shadow-ray ε-offset; no acne by construction |
| Derivatives | smooth field + silhouette delta | smooth part per chamber; the moving shadow boundary carries the Dirac part (Cell 6) |
| Totality tax | — | both spheres, both shadow quadratics, every pixel — a constant factor, honestly reported (Cell 4) |
| Failure modes | per-channel NaNs | cross-chamber NaN poisoning in raw Heaviside assembly (Cell 7a) |
§7 · Pitfalls gallery, continued
- Branch flips are not errors — but hiding them is. On boundary curves the two pipelines may legitimately choose different chambers in the last ulp. Report interior max and flip census as two numbers; a single "max |ΔI|" now lies by construction.
- Cross-chamber NaNs. In raw Heaviside assembly, $0 \cdot \mathrm{NaN}$ from the other surface's expression poisons pixels this surface owns. Regularize every √, or keep select semantics.
- The self-shadow knife edge. Testing a sphere against itself starts the segment on the blocker's own surface: $f(0)=0$ sits exactly on the root boundary, where strict and non-strict tests disagree in the last ulp. Convex self-shadowing is the terminator clamp $\max(0,\mu)$ — Theorem 3.2 makes this a proof; exclude self-blockers by construction, and say so.
- Penumbra ≠ soft shadow. Two point lights give a genuine penumbra as XOR of independent masks, but true soft shadows are integrals over area sources — Chapter-later territory, not a bug.
- Ordering ties. $t_A = t_B$ on the intersection curve of two surfaces; the $\le$ choice is a convention. Numeric and symbolic must share the convention or the flip census lights up the entire overlap curve.
- The constant-order theorem has hypotheses. Disjoint, convex silhouettes, one camera. Interpenetrating spheres break the nonvanishing condition and the once-per-frame reading dies with them — which is why the working simulator pays the per-pixel comparison and treats the theorem as what it is: a proof about the scene, not an optimization mandate.
§8 · A brief history, continued
Shadow feelers are as old as rendering: Arthur Appel's 1968 "Some Techniques for Shading Machine Renderings of Solids" already shot rays toward the light — before ray tracing had a name. Whitted's 1980 recursive ray tracer made shadows automatic and bought the shadow-acne that the ε-offset hack has treated ever since; Williams' shadow mapping (1978) took the approximate depth-buffer route that games still live on. This chapter is the other tradition: the shadow test as an exact sign predicate, with no ε anywhere. The quantifier-elimination lineage is older than graphics: Sturm (1829) counted real roots in intervals by sign variations; Tarski (1951) proved every polynomial scene admits QE; Collins' cylindrical algebraic decomposition (1975) made it mechanical — and expensive. The Shadow Predicate is the quadratic case of Sturm, small enough to hold in one hand: the whole frontier of Chapter 1's Cell 8, collapsed to three signs, for exactly as long as the world stays quadric. The self-shadow equivalence (Thm 3.2) is older still in spirit: Lambert's cosine law was a visibility statement before visibility had a name.
§9 · Roadmap
- Ch. 4 — Planar and spherical reflection: the reflected ray is a substituted expression (Cell 8's monster, tamed by CSE); recursion becomes composition; one bounce, honestly measured.
- Ch. 5 — Texture maps as symbolic functions u(x,y); the sRGB gamma question faced at last (we still render linear and show it raw).
- Ch. 6 — Symbolic antialiasing, now over chambers: integrate the piecewise image over the pixel square; the chamber boundaries are algebraic curves, and exact integration meets them head-on.
- Ch. 7 — Area lights and soft shadows: penumbra as an integral over the source; the XOR chamber becomes a quadrature.
Section ids are stable (s1…s9, sim, lab, cell1…cell8) — cite the id when requesting revisions.
Merged from index3a.html (theorems 3.1–3.4, standing assumptions) and index3b.html (working simulator, cells) — the audit is on record.