Symbolic Rendering

Chapter 8 — Glossy Reflection, or the Lobe Integral

An explainer comparing numerical rendering with symbolic rendering,
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
Python-ready for Jupyter Notebooks in Google Colab.
Requires an internet connection for mathematical typesetting (MathJax).
❦
How to use this file. §1 what changes; §2 the scene; §3 the derivation — the lobe integral, the mirror anchor, the exact Lambert weighting, and the roughness ruler σθ = 1/√(2n). §4 is the live simulator: numeric jitters the lobe, symbolic integrates it; the gloss view paints the reflected energy; n → 4096 recovers Chapter 4's mirror. §5: Cells 1 and 3 replace Chapter 7's; Cells 2, 4–8 are new. §6–§9 as ever.
☰  contents
  1. What changes in Chapter 8
  2. The scene (Chapter-8 assumptions)
  3. The math, once, carefully
  4. Live simulator
  5. The Colab laboratory (Python / sympy)
  6. The honest ledger — Chapter-8 deltas
  7. Pitfalls gallery, continued
  8. A brief history, continued
  9. Roadmap

§1 · What changes in Chapter 8

Chapter 7's scene stands (disk key light with its conic penumbra, point fill, textures, analytic floor AA). Chapter 4's mirror returns — glossy:

  1. The floor reflects with a lobe. Its mirror weight $k_{r,f}$ now carries a roughness exponent $n_g$: the reflected radiance is the integral of Chapter 4's composed chamber expression over a cosine-power lobe around the mirror direction.
  2. The anchor is the mirror itself. At $n_g \to \infty$ the lobe collapses onto the mirror direction and this chapter must return Chapter 4's image — Cell 6 runs that regression; the simulator exposes it at the top of the shininess slider.
  3. The estimators keep their currencies. Numeric jitters within strata (variance: speckle that thickens as the lobe widens); symbolic stratifies deterministically (bias), with an exact Lambert weighting option: stratum solid-angle weights $w = 2\pi(\cos^q\theta_a - \cos^q\theta_b)/q$ in closed form — the book's running "the cognoscenti will recognize" gag, finally recognized.

§2 · The scene (Chapter-8 assumptions)

Chapter 7's scene, plus: the floor is glossy (default $k_{r,f} = 0.3$, $n_g = 16$); the spheres' micro-reflection is fixed at $k_r = 0.05$ — present, sampled with the same lobe machinery, but the stage belongs to the floor, where the blur is wide and the estimator shapes are visible. One bounce (the lobe ray's hit shades directly, per Chapter 4's discipline); reflected shadow feelers reuse the Shadow Lemma verbatim. Units follow the series convention — reflected contribution is the lobe-weighted average — with the normalized-BRDF bookkeeping derived and confessed in Cell 2.

§3 · The math, once, carefully

The lobe integral. With mirror direction $\hat m = D - 2(D\cdot N)N$ at a floor point $P$, lobe half-angle $\theta_{\max}$ (default $\pi/3$), and gloss exponent $n_g$:

$$I_{\text{refl}}(P) = k_{r,f}\cdot \frac{\displaystyle\int_{\Omega} I_{\text{scene}}(P, \omega)\, (\cos\theta)^{n_g}\,d\omega} {\displaystyle\int_{\Omega}(\cos\theta)^{n_g}\,d\omega} \;=\; k_{r,f}\,\frac{\sum_k w_k I_k}{\sum_k w_k},$$

where $\theta$ is measured from $\hat m$, the denominator (series units) is $2\pi\big(1 - \cos^{n_g+1}\theta_{\max}\big)/(n_g+1)$ — a closed form the symbolic side uses as its normalizer — and $I_{\text{scene}}$ is Chapter 4's composed expression: sphere re-hits by the strict $t' \gt 0$ test, floor self-hits excluded by construction ($D'_z \lt 0$ from a floor point hits the floor again only behind the ray; horizon and miss both contribute background, which is zero — the boundary-of-domain honesty is in §7), shadows by the lemma, signs only.

Theorem (Mirror anchor). As $n_g \to \infty$ with $\theta_{\max}$ fixed, the lobe mass concentrates at $\theta = 0$: $I_{\text{refl}} \to k_{r,f}\,I_{\text{scene}}(P, \hat m)$ — Chapter 4's mirror term exactly. Proof: the weight ratio becomes a point mass, $\cos^{n_g}$ being a nascent delta on $[0, \theta_{\max}]$; Cell 6 runs the scene through it. ∎

The exact Lambert weighting (the cognoscenti clause, honored). On stratum $[\theta_a, \theta_b] \times \Delta\varphi$, the lobe weight integrates exactly: $$w_k = \Delta\varphi \int_{\theta_a}^{\theta_b} \cos^{n_g}\theta \sin\theta\,d\theta = \frac{\Delta\varphi}{n_g+1}\Big(\cos^{n_g+1}\theta_a - \cos^{n_g+1}\theta_b\Big).$$ The numeric estimator can't use this without abandoning its jittered creed; the symbolic estimator adopts it freely — stratum geometry exact, integrand evaluated at centroids — turning one factor of the quadrature error into zero. Cell 2 verifies the weights against fine quadrature and measures the improvement on the scene.

The roughness ruler. A cosine-power lobe has angular spread $\sigma_\theta \approx 1/\sqrt{2n_g}$: halving the blur costs 4× the exponent. The lobe is Chapter 5's texture-frequency story retold in angles — and the ruler is printed in the simulator's readout, so blur has units.

Three gifts of Chapter 8: (1) a composed integrand, exact and differentiable, doing real work inside an integral — the role the symbolic object has been rehearsing since Cell 6 of Chapter 1; (2) the mirror regression as a one-slider theorem; (3) the estimator contrast at its most visible — numeric speckle thickens with the blur, symbolic bias sharpens into banding at low K, and the difference canvas shows both at once. Counterpoint: the lobe crosses chambers — the integrand is piecewise along algebraic curves inside the integration domain, so quadrature loses its textbook order at the lobe's interior boundaries, and neither estimator may claim otherwise. And the energy-normalized units (the $(n_g+2)/2\pi$ divisor of the BRDF textbooks) are stated, derived, and then set aside per house convention — Cell 2 shows the renormalization is one line and the images agree after it.

§4 · Live simulator

The floor blurs its mirror; the spheres keep a faint glossy kiss. Numeric jitters the lobe (hashed per-pixel — reproducible, but visibly speckled); symbolic stratifies it (deterministic; exact Lambert weights optional). The ledger keeps RMS and max; the gloss view paints the reflected term alone.

model 
resolution 
initializing…
NUMERIC — jittered lobe samples; speckle thickens with blur
SYMBOLIC — stratified lobe quadrature; exact weights, deterministic
|difference| ×10² — speckle vs banding, the two currencies of estimator error
numeric: – symbolic: – RMS |ΔI|: – max |ΔI|: – σθ = –
hover — pixel → world, surface, gloss fraction, and both RGB triples.

Try this: (1) gloss n to 4096: the lobe collapses — both panels return Chapter 4's mirror, and the difference snaps back toward machine black (the anchor regression, live). (2) n = 4, K = 2: numeric speckle is obvious; symbolic shows smooth bands. Two currencies, one integrand. (3) Flip symbolic weights to centroid: a faint systematic shift appears in the diff — the price of midpoint weights, made visible. (4) gloss view: the reflected term alone — watch it widen and fade as n drops; that is σθ = 1/√(2n), read directly off the readout.

§5 · The Colab laboratory (Python / sympy)

Cumulative: Cells 1 and 3 replace Chapter 7's; Cells 2, 4–8 are new. Paste 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 (REPLACES Ch.7 Cell 1) — numeric baseline WITH a glossy floor: jittered lobe sampling

· · paste into colab.new → Run
import numpy as np
import matplotlib.pyplot as plt
import math, time

SC8 = dict(R=1.0, e=4.0, Cb=(0.5, 0.5, 1.3), R2=0.32, b_on=True, floor_on=True,
           L1=(5., 5., 10.), rho=0.45, c1=(1., 1., 1.),
           L2=(-6., 2., 8.), c2=(0.9, 0.65, 0.25), l2_on=True,
           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',
           tex='checker', ku=8.0, kv=4.0, alpha=0.25,
           krf=0.30, gloss_n=16.0, K=4, seed=7, thmax=math.pi/3)

def srgb_encode(c):
    c = np.maximum(np.asarray(c, dtype=float), 0.0)
    return np.where(c <= 0.0031308, 12.92*c, 1.055*np.power(c, 1/2.4) - 0.055)

def disk_basis(L, target=(0., 0., 0.)):
    w = np.array(target, dtype=float) - np.array(L, dtype=float)
    w /= np.linalg.norm(w)
    up = np.array([0., 1., 0.]) if abs(w[1]) < 0.9 else np.array([1., 0., 0.])
    e1 = np.cross(w, up); e1 /= np.linalg.norm(e1)
    e2 = np.cross(w, e1)
    return e1, e2

def chartA1(P, R):
    Nx, Ny, Nz = P[0]/R, P[1]/R, P[2]/R
    return ((math.atan2(Ny, Nx) + math.pi)/(2*math.pi),
            math.acos(max(-1.0, min(1.0, Nz)))/math.pi)

def texture1(u, v, SC=SC8):
    m = SC['matA']; t = SC['tex']
    if t == 'solid': alb = m
    elif t == 'stripes':
        b = 1.0 if math.sin(math.pi*u*SC['ku']) >= 0 else 0.2
        alb = (m[0]*b, m[1]*b, m[2]*b)
    elif t == 'marble':
        mm = 0.5 + 0.5*math.sin(2*math.pi*(SC['ku']*u + 0.35*math.sin(3*math.pi*SC['kv']*v)))
        b = 0.55 + 0.45*mm; alb = (m[0]*b, m[1]*b, m[2]*b)
    else:
        s = math.sin(math.pi*u*SC['ku'])*math.sin(math.pi*v*SC['kv'])
        b = 1.0 if s >= 0 else 0.2; alb = (m[0]*b, m[1]*b, m[2]*b)
    return np.array(alb)

def floor_tex1(x, y, SC=SC8):
    s = 1.0 if math.sin(math.pi*x)*math.sin(math.pi*y) >= 0 else -1.0
    base = 0.85 if s > 0 else 0.15
    return np.array((base*(1-SC['alpha']), base*(1-SC['alpha']), base))

def shadowed1(P, Q, C, Rr):
    seg = Q - P
    a2 = float(seg @ seg); b2 = float((P - C) @ seg)
    c2 = float((P - C) @ (P - C)) - Rr*Rr
    if a2 < 1e-18: return False
    dsc = b2*b2 - a2*c2
    if dsc >= 0:
        sq = math.sqrt(dsc); s1, s2 = (-b2 - sq)/a2, (-b2 + sq)/a2
        if (0 < s1 < 1) or (0 < s2 < 1): return True
    return False

def light_point(P, N, Q, v_hat, alb, SC, blockers):
    if any(shadowed1(P, Q, C, Rr) for C, Rr in blockers):
        return np.zeros(3), 0.0
    w = Q - P; wn = np.linalg.norm(w); w_hat = w/wn
    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']
    return SC['kd']*alb*max(0.0, mu) + SC['ks']*spec, 1.0

def hit_sphere1(Q0, Dv, C, Rr):
    """Nearest t > 0 (strict); None if missed. The bounce's exclusion test."""
    a = float(Dv @ Dv); b = float((Q0 - C) @ Dv)
    c = float((Q0 - C) @ (Q0 - C)) - Rr*Rr
    dsc = b*b - a*c
    if dsc < 0: return None
    sq = math.sqrt(dsc)
    for t in sorted(((-b - sq)/a, (-b + sq)/a)):
        if t > 0: return t
    return None

def scene_at(Q0, Dv, SC=SC8):
    """Chapter 4's composed scene, one level: what direction Dv from Q0 sees."""
    R, e = SC['R'], SC['e']; Cb = np.array(SC['Cb']); R2 = SC['R2']
    tA = hit_sphere1(Q0, Dv, np.zeros(3), R)
    tB = hit_sphere1(Q0, Dv, Cb, R2) if SC['b_on'] else None
    tF = -Q0[2]/Dv[2] if (SC['floor_on'] and Dv[2] < 0 and abs(Q0[2]) > 1e-12) else None
    cand = [(t, k) for t, k in ((tA, 0), (tB, 1), (tF, 2)) if t is not None]
    if not cand: return np.zeros(3)
    t, kind = min(cand)
    Q = Q0 + t*Dv
    if kind == 0:
        N = Q/R; alb = texture1(*chartA1(Q, R), SC)
        blockers = [(Cb, R2)] if SC['b_on'] else []
    elif kind == 1:
        N = (Q - Cb)/R2; alb = np.array(SC['matB']); blockers = [(np.zeros(3), R)]
    else:
        N = np.array([0., 0., 1.]); alb = floor_tex1(Q[0], Q[1], SC)
        blockers = [(np.zeros(3), R)] + ([(Cb, R2)] if SC['b_on'] else [])
    v_hat = (np.array([0., 0., e]) - Q); v_hat /= np.linalg.norm(v_hat)
    col = SC['ka']*np.array(SC['ambient'])
    # lights: disk (strata of K, cheap: reuse K=2) + point
    e1, e2 = disk_basis(SC['L1'])
    Kd = 2
    acc = np.zeros(3)
    for sj in range(Kd):
        for si in range(Kd):
            u, v = (si + .5)/Kd, (sj + .5)/Kd
            rr = SC['rho']*math.sqrt(u); th = 2*math.pi*v
            Qp = np.array(SC['L1']) + rr*(math.cos(th)*e1 + math.sin(th)*e2)
            c, _ = light_point(Q, N, Qp, v_hat, alb, SC, blockers)
            acc += np.array(SC['c1'])*c
    col += acc/(Kd*Kd)
    if SC['l2_on']:
        c, _ = light_point(Q, N, np.array(SC['L2']), v_hat, alb, SC, blockers)
        col += np.array(SC['c2'])*c
    return col

def render_numeric_gloss(width, height, SC=SC8):
    """Numeric: primary shading as Ch.7 + glossy floor by JITTERED lobe samples."""
    rng = np.random.default_rng(SC['seed'])
    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), np.uint8)
    gloss = np.zeros((height, width, 3))
    half = 1.25*e*R/np.sqrt(e**2 - R**2)
    E = np.array([0.0, 0.0, e])
    e1, e2 = disk_basis(SC['L1'])
    K = SC['K']
    thm = SC['thmax']
    for j in range(height):
        y = half*(1 - 2*(j + 0.5)/height) - 0.15*half
        for i in range(width):
            x = half*(2*(i + 0.5)/width - 1)
            Aq = x*x + y*y + e*e
            discA = e*e*R*R - (e*e - R*R)*(x*x + y*y)
            tA = (e*e - math.sqrt(discA))/Aq if discA >= 0 else math.inf
            tB = math.inf
            if SC['b_on']:
                bB = -Cb[0]*x - Cb[1]*y - (e - Cb[2])*e
                cB = float(Cb @ Cb) - R2*R2
                discB = bB*bB - Aq*cB
                if discB >= 0: tB = (-bB - math.sqrt(discB))/Aq
            tf = 1.0 if SC['floor_on'] else math.inf
            tm = min(tA, tB, tf)
            if tm == math.inf: continue
            P = np.array([tm*x, tm*y, e*(1 - tm)])
            if tm == tA:
                N = P/R; alb = texture1(*chartA1(P, R), SC)
                blockers = [(Cb, R2)] if SC['b_on'] else []; kind = 1
            elif tm == tB:
                N = (P - Cb)/R2; alb = np.array(SC['matB'])
                blockers = [(np.zeros(3), R)]; kind = 2
            else:
                N = np.array([0., 0., 1.]); alb = floor_tex1(P[0], P[1], SC)
                blockers = [(np.zeros(3), R)] + ([(Cb, R2)] if SC['b_on'] else []); kind = 3
            v_hat = (E - P)/np.linalg.norm(E - P)
            col = SC['ka']*np.array(SC['ambient'])
            acc = np.zeros(3)
            for sj in range(2):
                for si in range(2):
                    u, v = (si + rng.random())/2, (sj + rng.random())/2
                    rr = SC['rho']*math.sqrt(u); th = 2*math.pi*v
                    Qp = np.array(SC['L1']) + rr*(math.cos(th)*e1 + math.sin(th)*e2)
                    c, _ = light_point(P, N, Qp, v_hat, alb, SC, blockers)
                    acc += np.array(SC['c1'])*c
            col += acc/4
            if SC['l2_on']:
                c, _ = light_point(P, N, np.array(SC['L2']), v_hat, alb, SC, blockers)
                col += np.array(SC['c2'])*c
            # ---- GLOSSY FLOOR: jittered lobe samples ----
            g = np.zeros(3)
            if kind == 3 and SC['krf'] > 0:
                Dv = np.array([x, y, -e]); Dv /= np.linalg.norm(Dv)
                m = Dv - 2.0*float(Dv @ N)*N              # mirror direction
                # frame about m
                up = np.array([0., 1., 0.]) if abs(m[1]) < 0.9 else np.array([1., 0., 0.])
                u1 = np.cross(m, up); u1 /= np.linalg.norm(u1)
                u2 = np.cross(m, u1)
                wsum = 0.0
                for sj in range(K):
                    for si in range(K):
                        uu = (si + rng.random())/K; vv = (sj + rng.random())/K
                        ct = 1 - uu*(1 - math.cos(thm))
                        st = math.sqrt(max(0.0, 1 - ct*ct)); ph = 2*math.pi*vv
                        Dp = ct*m + st*(math.cos(ph)*u1 + math.sin(ph)*u2)
                        w = ct**SC['gloss_n']
                        g += w*scene_at(P, Dp, SC); wsum += w
                if wsum > 0: g /= wsum
                col = col + SC['krf']*g
            img[j, i] = col; code[j, i] = kind; gloss[j, i] = g
    return np.clip(img, 0, 1), code, gloss

t0 = time.perf_counter()
img_num, code_num, gloss_num = render_numeric_gloss(256, 256)
t1 = time.perf_counter()
print(f"numeric render (glossy floor, K={SC8['K']}x{SC8['K']} jittered): {t1-t0:.2f} s")
fig, axs = plt.subplots(1, 2, figsize=(9, 4))
axs[0].imshow(srgb_encode(img_num)); axs[0].set_title("Numeric: blurred mirror")
axs[1].imshow(srgb_encode(np.clip(gloss_num*2, 0, 1)))
axs[1].set_title("gloss term alone (x2)")
for ax in axs: ax.axis('off')
plt.show()
Expect: the mirror image of the spheres smeared into a soft smudge; the gloss view shows the blurred reflections with visible speckle at K=4.

🐍 Cell 2 (NEW) — the normalizer in closed form, the exact Lambert weights, and the units confession

· · paste into colab.new → Run
import sympy as sp

n, th, thm = sp.symbols('n theta theta_max', positive=True)
# --- the series-units normalizer, in closed form ---------------------------------
W = 2*sp.pi*sp.integrate(sp.cos(th)**n * sp.sin(th), (th, 0, thm))
print("lobe normalizer: 2*pi*int cos^n sin =", sp.simplify(W))
# -> 2*pi*(1 - cos^(n+1)(th_max))/(n+1)

# --- exact stratum weights: w = dphi*(cos^(n+1) a - cos^(n+1) b)/(n+1) ------------
a, b = sp.symbols('a b', positive=True)
w_exact = sp.integrate(sp.cos(th)**n * sp.sin(th), (th, a, b))
print("stratum weight:", sp.simplify(w_exact))

# verify against quadrature for a sample lobe
def wq(nv, av, bv, K=20000):
    tt = np.linspace(av, bv, K)
    return np.trapezoid(np.cos(tt)**nv*np.sin(tt), tt)
for nv, av, bv in ((16, 0.0, 0.5), (16, 0.5, 1.0), (64, 0.1, 0.2)):
    we = float(w_exact.subs({n: nv, a: av, b: bv}))
    print(f"  n={nv} [{av},{bv}]: exact {we:.10e}  quad {wq(nv, av, bv):.10e}")

# --- the midpoint-weight tax, measured -------------------------------------------
def weights(K, nv, thm_v, exact=True):
    ws = []
    for k in range(K):
        ta = k/K*thm_v; tb = (k + 1)/K*thm_v
        if exact:
            ws.append((math.cos(ta)**(nv + 1) - math.cos(tb)**(nv + 1))/(nv + 1))
        else:
            tm = (ta + tb)/2
            ws.append(math.cos(tm)**nv*math.sin(tm)*(tb - ta))
    return np.array(ws)
for K in (2, 4, 8):
    tot_ex = weights(K, 16, math.pi/3, True).sum()
    tot_mp = weights(K, 16, math.pi/3, False).sum()
    tot_true = (1 - math.cos(math.pi/3)**17)/17
    print(f"K={K}: exact-weight sum error {abs(tot_ex - tot_true):.2e}, "
          f"midpoint sum error {abs(tot_mp - tot_true):.2e}")

# --- the units confession ---------------------------------------------------------
print("""
Series units: the reflected term is the lobe-weighted AVERAGE of the scene.
BRDF units (normalized Phong): multiply by (n+2)/(2*pi) and weight by the
surface cosine. The renorm is one line; the images agree after it:
""")
for nv in (4, 16, 64):
    print(f"  n={nv}: (n+2)/(2*pi) = {(nv + 2)/(2*math.pi):.4f}   "
          f"(energy normalization, derived and set aside per house convention)")
Expect: the closed-form normalizer and stratum weights, quadrature agreement to ~1e-10, and the midpoint tax printed per K.

🐍 Cell 3 (REPLACES Ch.7 Cell 3) — the symbolic renderer: stratified lobe quadrature with exact weights; the estimator census

· · paste into colab.new → Run
def render_symbolic_gloss(width, height, SC=SC8, exact_w=True):
    """Symbolic: SAME integrand, deterministic strata; exact Lambert stratum
       weights (exact_w) or centroid weights (the midpoint tax, Cell 2)."""
    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), np.uint8)
    gloss = np.zeros((height, width, 3))
    half = 1.25*e*R/np.sqrt(e**2 - R**2)
    E = np.array([0.0, 0.0, e])
    e1, e2 = disk_basis(SC['L1'])
    K = SC['K']; nv = SC['gloss_n']; thm = SC['thmax']
    # strata: rings x sectors, centroids + exact weights
    strata = []
    for sj in range(K):
        for si in range(K):
            ta = si/K*thm; tb = (si + 1)/K*thm
            pa, pb = sj/K*2*math.pi, (sj + 1)/K*2*math.pi
            tm, pm = (ta + tb)/2, (pa + pb)/2
            if exact_w:
                w = ((math.cos(ta)**(nv + 1) - math.cos(tb)**(nv + 1))/(nv + 1)) \
                    * (pb - pa)
            else:
                w = math.cos(tm)**nv*math.sin(tm)*(tb - ta)*(pb - pa)
            strata.append((tm, pm, w))
    Wtot = sum(w for _, _, w in strata)
    for j in range(height):
        y = half*(1 - 2*(j + 0.5)/height) - 0.15*half
        for i in range(width):
            x = half*(2*(i + 0.5)/width - 1)
            Aq = x*x + y*y + e*e
            discA = e*e*R*R - (e*e - R*R)*(x*x + y*y)
            tA = (e*e - math.sqrt(discA))/Aq if discA >= 0 else math.inf
            tB = math.inf
            if SC['b_on']:
                bB = -Cb[0]*x - Cb[1]*y - (e - Cb[2])*e
                cB = float(Cb @ Cb) - R2*R2
                discB = bB*bB - Aq*cB
                if discB >= 0: tB = (-bB - math.sqrt(discB))/Aq
            tf = 1.0 if SC['floor_on'] else math.inf
            tm2 = min(tA, tB, tf)
            if tm2 == math.inf: continue
            P = np.array([tm2*x, tm2*y, e*(1 - tm2)])
            if tm2 == tA:
                N = P/R; alb = texture1(*chartA1(P, R), SC)
                blockers = [(Cb, R2)] if SC['b_on'] else []; kind = 1
            elif tm2 == tB:
                N = (P - Cb)/R2; alb = np.array(SC['matB'])
                blockers = [(np.zeros(3), R)]; kind = 2
            else:
                N = np.array([0., 0., 1.]); alb = floor_tex1(P[0], P[1], SC)
                blockers = [(np.zeros(3), R)] + ([(Cb, R2)] if SC['b_on'] else []); kind = 3
            v_hat = (E - P)/np.linalg.norm(E - P)
            col = SC['ka']*np.array(SC['ambient'])
            acc = np.zeros(3)
            for sj in range(2):
                for si in range(2):
                    u, v = (si + .5)/2, (sj + .5)/2
                    rr = SC['rho']*math.sqrt(u); th = 2*math.pi*v
                    Qp = np.array(SC['L1']) + rr*(math.cos(th)*e1 + math.sin(th)*e2)
                    c, _ = light_point(P, N, Qp, v_hat, alb, SC, blockers)
                    acc += np.array(SC['c1'])*c
            col += acc/4
            if SC['l2_on']:
                c, _ = light_point(P, N, np.array(SC['L2']), v_hat, alb, SC, blockers)
                col += np.array(SC['c2'])*c
            # ---- GLOSSY FLOOR: stratified lobe quadrature, exact weights ----
            g = np.zeros(3)
            if kind == 3 and SC['krf'] > 0:
                Dv = np.array([x, y, -e]); Dv /= np.linalg.norm(Dv)
                m = Dv - 2.0*float(Dv @ N)*N
                up = np.array([0., 1., 0.]) if abs(m[1]) < 0.9 else np.array([1., 0., 0.])
                u1 = np.cross(m, up); u1 /= np.linalg.norm(u1)
                u2 = np.cross(m, u1)
                for tm_, pm, w in strata:
                    ct, st = math.cos(tm_), math.sin(tm_)
                    Dp = ct*m + st*(math.cos(pm)*u1 + math.sin(pm)*u2)
                    g += w*scene_at(P, Dp, SC)
                g /= Wtot
                col = col + SC['krf']*g
            img[j, i] = col; code[j, i] = kind; gloss[j, i] = g
    return np.clip(img, 0, 1), code, gloss

t0 = time.perf_counter()
img_sym, code_sym, gloss_sym = render_symbolic_gloss(256, 256)
t1 = time.perf_counter()
print(f"symbolic render (stratified lobe, exact weights): {t1-t0:.2f} s")

same = code_num == code_sym
d = np.abs(img_num - img_sym)[same]
print(f"interior max  |num - sym| = {d.max():.3e}")
print(f"interior RMS  |num - sym| = {np.sqrt((d**2).mean()):.3e}   <- estimator error")
print(f"flips: {(~same).sum()} px")
print("Reading: disagreement concentrates in the glossy floor where the numeric")
print("speckle lives; the symbolic side is smooth (deterministic bias only).")

fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(srgb_encode(img_num)); axs[0].set_title("Numeric (jittered lobe)")
axs[1].imshow(srgb_encode(img_sym)); axs[1].set_title("Symbolic (strata, exact w)")
axs[2].imshow(np.clip(np.abs(img_num - img_sym)*30, 0, 1))
axs[2].set_title("|delta| x30: speckle vs smooth bias")
for ax in axs: ax.axis('off')
plt.show()
Expect: a cleaner blur on the symbolic side; the ×30 diff showing numeric speckle concentrated in the floor's reflection; RMS dominated by variance.

🐍 Cell 4 (NEW) — honest timings: the lobe multiplies Chapter 7's source loop

· · paste into colab.new → Run
print(f"{'K (lobe)':>9} {'numeric 64^2':>14} {'symbolic 64^2':>15} {'sym centroid-w':>15}")
for K in (2, 4, 8):
    t0 = time.perf_counter(); render_numeric_gloss(64, 64, {**SC8, 'K': K})
    t1 = time.perf_counter(); render_symbolic_gloss(64, 64, {**SC8, 'K': K}, exact_w=True)
    t2 = time.perf_counter(); render_symbolic_gloss(64, 64, {**SC8, 'K': K}, exact_w=False)
    t3 = time.perf_counter()
    print(f"{K:>9} {t1-t0:>13.2f}s {t2-t1:>14.2f}s {t3-t2:>14.2f}s")

print("""
Reading the table:
 * Lobe samples multiply the scene_at cost: K^2 x (scene re-trace incl. a small
   disk-light loop). Both pipelines pay; exact vs centroid weights differ by a
   rounding error in TIME -- the weight tax is in ACCURACY (Cell 2), not speed.
 * The honest accounting of the whole series so far:
     Chapter 6: pixel integral   (M^2 on numeric; band-confined S^2 on symbolic)
     Chapter 7: source integral  (K^2 on both; currencies differ)
     Chapter 8: lobe integral    (K^2 on both, MULTIPLIED by the source loop)
   Integrals are stacking. The book's closing chapters live in the product.
""")
Expect: both columns growing ~K²; centroid vs exact rows indistinguishable in time — the weight choice is an accuracy decision, not a speed one.

🐍 Cell 5 (NEW) — superpower study: gloss convergence in both currencies, against a K=16 exact-weight truth

· · paste into colab.new → Run
W = 96
print("computing K=16 exact-weight reference (patience)...")
img_true, code_true, _ = render_symbolic_gloss(W, W, {**SC8, 'K': 16}, exact_w=True)

print(f"\n{'K':>4} {'numeric RMS (s7)':>17} {'numeric RMS (s99)':>18} "
      f"{'sym RMS (exact-w)':>18} {'sym RMS (centroid)':>19}")
rows = []
for K in (2, 4, 8):
    row = [K]
    for seed in (7, 99):
        im, cn, _ = render_numeric_gloss(W, W, {**SC8, 'K': K, 'seed': seed})
        same = cn == code_true
        d = np.abs(im - img_true)[same]
        row.append(np.sqrt((d**2).mean()))
    im, cn, _ = render_symbolic_gloss(W, W, {**SC8, 'K': K}, exact_w=True)
    same = cn == code_true
    d = np.abs(im - img_true)[same]
    row.append(np.sqrt((d**2).mean()))
    im, cn, _ = render_symbolic_gloss(W, W, {**SC8, 'K': K}, exact_w=False)
    d = np.abs(im - img_true)[same]
    row.append(np.sqrt((d**2).mean()))
    rows.append(row)
    print(f"{K:>4} {row[1]:>17.4e} {row[2]:>18.4e} {row[3]:>18.4e} {row[4]:>19.4e}")

Ks = [r[0] for r in rows]
plt.figure(figsize=(6.5, 3.8))
plt.plot(Ks, [r[1] for r in rows], 'o-', color='#7a1f1f', label='numeric (jittered)')
plt.plot(Ks, [r[3] for r in rows], 's-', color='#1f3a5f', label='symbolic (exact w)')
plt.plot(Ks, [r[4] for r in rows], '^--', color='#b8860b', label='symbolic (centroid)')
plt.yscale('log'); plt.xlabel('K'); plt.ylabel('interior RMS vs K=16 truth')
plt.title('the lobe integral, priced in both currencies'); plt.legend(); plt.show()

print("\nThe exact weights buy a visible margin at low K (where the strata are fat")
print("and the midpoint assumption is weakest) -- Cell 2's ±2pi-improvement, spent.")
print("Numeric keeps its unbiasedness certificate; symbolic keeps its determinism.")
print("Same integrand. Same lobe. Two honest estimators, each with a receipt.")
Expect: three decaying curves; exact-weight symbolic below centroid at low K; numeric noisy but unbiased; the receipt columns match the story.

🐍 Cell 6 (NEW) — the mirror anchor regression: n → ∞ returns Chapter 4, in both pipelines

· · paste into colab.new → Run
# As n_g -> infinity the lobe becomes the mirror direction. This chapter must
# return Chapter 4's mirror image. (Mirror-mode: thmax tiny, n huge, K=1
# centroid = the mirror ray.)
SC_mir = {**SC8, 'gloss_n': 1e6, 'K': 1, 'thmax': 1e-6}
img_ms, code_ms, gloss_ms = render_symbolic_gloss(128, 128, SC_mir)

# reference: scene_at along the EXACT mirror direction, no lobe at all
def mirror_floor(P, x, y, SC=SC8):
    e = SC['e']
    Dv = np.array([x, y, -e]); Dv /= np.linalg.norm(Dv)
    N = np.array([0., 0., 1.])
    m = Dv - 2.0*float(Dv @ N)*N
    return scene_at(P, m, SC)

W = 128
half = 1.25*SC8['e']*SC8['R']/np.sqrt(SC8['e']**2 - SC8['R']**2)
ref = np.zeros((W, W, 3))
for j in range(W):
    y = half*(1 - 2*(j + 0.5)/W) - 0.15*half
    for i in range(W):
        x = half*(2*(i + 0.5)/W - 1)
        if code_ms[j, i] == 3:
            ref[j, i] = mirror_floor(np.array([x, y, 0.0]), x, y, SC8)

floor_px = code_ms == 3
d = np.abs(gloss_ms - ref)[floor_px]
print(f"lobe(n=1e6, K=1) vs exact mirror direction, floor pixels:")
print(f"  max |delta| = {d.max():.3e}   mean {d.mean():.3e}")
print("The anchor holds: gloss closes on the mirror, as the theorem requires.")

# and the smooth sweep: watch the lobe converge as n grows
print("\ngloss at a fixed floor point vs n (approaching the mirror value):")
P0 = np.array([0.35, -0.62*half + 0.0, 0.0])
x0, y0 = 0.35, -0.62*half
Dv = np.array([x0, y0, -SC8['e']]); Dv /= np.linalg.norm(Dv)
N = np.array([0., 0., 1.])
m = Dv - 2.0*float(Dv @ N)*N
mir = scene_at(P0, m, SC8)
for nv in (2, 8, 32, 128, 2048):
    K = 16
    thm_v = SC8['thmax']
    acc = np.zeros(3); Wt = 0.0
    up = np.array([0., 1., 0.]) if abs(m[1]) < 0.9 else np.array([1., 0., 0.])
    u1 = np.cross(m, up); u1 /= np.linalg.norm(u1); u2 = np.cross(m, u1)
    for sj in range(K):
        for si in range(K):
            ta, tb = si/K*thm_v, (si + 1)/K*thm_v
            pm = (sj + .5)/K*2*math.pi
            w = ((math.cos(ta)**(nv + 1) - math.cos(tb)**(nv + 1))/(nv + 1))*(2*math.pi/K)
            tm = (ta + tb)/2
            Dp = math.cos(tm)*m + math.sin(tm)*(math.cos(pm)*u1 + math.sin(pm)*u2)
            acc += w*scene_at(P0, Dp, SC8); Wt += w
    est = acc/Wt
    print(f"  n={nv:>6}: |lobe - mirror| = {np.abs(est - mir).max():.4e}")
Expect: machine-level agreement at the anchor; and the sweep showing the lobe estimate gliding onto the mirror value as n grows — the nascent delta, watched.

🐍 Cell 7 (NEW) — pitfalls, live: the horizon boundary, lobe-chamber crossings, and weight-noise when n is large

· · paste into colab.new → Run
# (a) THE HORIZON: from a floor point, directions with D'_z >= 0 see nothing.
#     The lobe's lower rim crosses the horizon at grazing view angles: the
#     integrand truncates mid-lobe, and both estimators must (and do) integrate
#     background = 0 there. Watch the gloss dim toward the horizon:
W = 128
img_s, code_s, gl_s = render_symbolic_gloss(W, W, {**SC8, 'K': 8})
profile = gl_s[:, W//2].mean(axis=-1)
print("gloss strength along the center column (top to bottom):")
for frac in (0.1, 0.3, 0.5, 0.7, 0.9):
    jj = int(frac*(W - 1))
    print(f"  row {jj:>3}: gloss = {profile[jj]:.4f}")
print("The floor's reflection dies toward the far edge: the lobe crosses the")
print("horizon, and the integrand's domain truncates. Not a bug -- a boundary.")

# (b) LOBE-CHAMBER CROSSINGS: inside one lobe, rays may span a silhouette, a
#     shadow boundary, AND the horizon -- quadrature order is lost at each.
#     Count crossings per lobe at a grazing floor point:
P0 = np.array([0.9, -0.35, 0.0])
x0, y0 = 0.9, -0.35
e = SC8['e']
Dv = np.array([x0, y0, -e]); Dv /= np.linalg.norm(Dv)
m = Dv - 2.0*float(Dv @ np.array([0., 0., 1.]))*np.array([0., 0., 1.])
up = np.array([0., 1., 0.]) if abs(m[1]) < 0.9 else np.array([1., 0., 0.])
u1 = np.cross(m, up); u1 /= np.linalg.norm(u1); u2 = np.cross(m, u1)
def what_hit(Dp):
    tA = hit_sphere1(P0, Dp, np.zeros(3), SC8['R'])
    tB = hit_sphere1(P0, Dp, np.array(SC8['Cb']), SC8['R2']) if SC8['b_on'] else None
    tF = -P0[2]/Dp[2] if Dp[2] < 0 and abs(P0[2]) > 1e-12 else None
    cand = [(t, k) for t, k in ((tA, 0), (tB, 1), (tF, 2)) if t is not None]
    return min(cand)[1] if cand else -1
prev = None; crossings = 0; Ks = 200
for k in range(Ks):
    th = k/Ks*SC8['thmax']
    Dp = math.cos(th)*m + math.sin(th)*u1
    h = what_hit(Dp)
    if prev is not None and h != prev: crossings += 1
    prev = h
print(f"\nlobe cross-section at a grazing point: {crossings} chamber crossings")
print("Each crossing is a place where quadrature's textbook order dies.")

# (c) WEIGHT NOISE AT HIGH n: with n huge and K small, most strata carry ~zero
#     weight; the estimate becomes 'whichever stratum survived' -- a bias spike.
for nv, K in ((4096, 2), (4096, 8)):
    thm_v = SC8['thmax']
    acc = np.zeros(3); Wt = 0.0
    for sj in range(K):
        for si in range(K):
            ta, tb = si/K*thm_v, (si + 1)/K*thm_v
            pm = (sj + .5)/K*2*math.pi
            w = ((math.cos(ta)**(nv + 1) - math.cos(tb)**(nv + 1))/(nv + 1))*(2*math.pi/K)
            tm = (ta + tb)/2
            Dp = math.cos(tm)*m + math.sin(tm)*(math.cos(pm)*u1 + math.sin(pm)*u2)
            acc += w*scene_at(P0, Dp, SC8); Wt += w
    print(f"  n=4096, K={K}: |est - mirror| = {np.abs(acc/Wt - mir_val if False else np.abs(acc/Wt - scene_at(P0, m, SC8))).max():.3e}")
print("At mirror-grade n, only strata containing theta=0 matter; the rest are")
print("dead weight. Importance-sampling the lobe IS the lobe -- the estimator")
print("must follow the integrand or starve. (Chapter 10's lesson, early.)")
Expect: the horizon dimming printed as a profile; chamber crossings counted at grazing angle; and the high-n/low-K bias spike exhibited with numbers.

🐍 Cell 8 (NEW) — the frontier, updated: the stacked integrals, and where the symbolic role settles

· · paste into colab.new → Run
print("""
The integrals have now stacked three deep, and the pattern is the theory:

  pixel square (Ch.6)  x  light source (Ch.7)  x  reflection lobe (Ch.8)

Each layer: the SAME chambered, sign-masked, differentiable integrand; an exact
anchor (checker integral / on-axis form / mirror direction); a regression
(M=1 / K=1 & rho->0 / n->inf); and two estimators with named currencies.

  Ch. 9   bump mapping: NO new integral -- the exact gradient of tau, composed
          (Chapter 4 machinery x Chapter 5 functions). A rest stop, earned.
  Ch. 10  motion blur & depth of field: time x aperture -- the aperture is
          literally this chapter's disk with the camera on it; the shutter is
          one more interval. The product integral becomes 5-dimensional, and
          the honest statement from Chapter 6 -- exact where the boundaries
          are algebraic, sampled where they are not -- becomes the whole law.
  Appendix path space: the integrand stays exact and differentiable forever;
          'evaluate it closed-form' died at Chapter 4's bounce 2, and that is
          all right. The symbolic object ends the book as what it was always
          becoming: the ground truth you can hold, differentiate, and sample --
          not the answer, the thing the answers are estimates OF.
""")
Expect: the roadmap rewritten as a product of integrals, and the symbolic object's final job title settled in print.

§6 · The honest ledger — Chapter-8 deltas

DimensionChapter 7Chapter 8
New integralover the light sourceover the reflection lobe (solid angle) — multiplied by the source loop
Anchoron-axis μ̄ closed formthe mirror direction: n → ∞ regression, live on the slider (Cell 6)
Weightsstrata uniformexact Lambert stratum weights — closed form, verified; centroid tax measured (Cells 2, 5)
Normalizerπρ² (trivial)2π(1−cos^{n+1}θ_max)/(n+1), closed form; BRDF units confessed (Cell 2)
Error lookpenumbra speckle vs bandingblur-thickened speckle vs weight-banding; three curves priced together (Cell 5)
New boundariesconic on the sourcethe horizon (domain truncation) and lobe-chamber crossings — quadrature order dies there, counted (Cell 7)
Starvation—at mirror-grade n, low-K strata starve; the estimator must follow the integrand (Cell 7c)
Roughness—σθ = 1/√(2n): blur has units, printed in the readout

§7 · Pitfalls gallery, continued

  1. The horizon truncates the lobe. From a floor point, half the lobe can face the void at grazing angles; background-zero must be integrated, not skipped — else the far floor glows. The profile in Cell 7a is the correct dim.
  2. Lobes cross chambers. One lobe can span silhouette, shadow boundary, and horizon simultaneously; quadrature order is lost at each crossing. Both estimators inherit this; only the symbolic side can count the crossings.
  3. Large n starves dumb strata. When the lobe is a needle, most strata carry no weight; the estimate is whichever stratum survived. Importance-sampling the lobe is not optional at high gloss (Cell 7c).
  4. Weights are accuracy, not speed. Exact vs centroid weights cost identical time (Cell 4) and different answers (Cell 5). Spend accuracy deliberately.
  5. The units confession matters at interop time. Series units average; BRDF units normalize. They differ by a factor with a name — write down which side of it your pipeline is on before comparing images with anyone, including your past self.

§8 · A brief history, continued

Glossy reflection is where the two traditions of this book merged in production decades before this book: Cook–Torrance (1982) brought microfacet theory's lobes into graphics; Ward's anisotropic BRDF (1992) made lobes respectable; and Kajiya's rendering equation (1986) wrote the entire field as one integral — the umbrella under which Chapters 6, 7, and 8 of this volume all stand. The stratified-versus- jittered argument is Cook's 1986 stochastic sampling paper in miniature; the exact Lambert stratum weights are the kind of closed form the radiosity era (Goral 1984, and the hemicube's critics) collected by the shelf. Blinn's microfacet work (1977) already knew the half-vector; this chapter merely integrates him. In the other world: the nascent delta sequence cosⁿθ is nineteenth-century analysis; the solid- angle strata with exact weights are spherical geometry homework from the age of astronomical tables. The lobe was always there. We finally priced it.

§9 · Roadmap

❦
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.