Symbolic Rendering

Chapter 6 — Antialiasing, or the Pixel Is an 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 checker box integral, the silhouette's exact area integral, the sphere-over-floor composition (the second printing's correction), and the optional cap-integrated normal. §4 is the live simulator. §5: Cells 1 and 3 replace Chapter 5's; Cell 2b is the silhouette integral; Cell 2c is the cap-normal refinement. §6–§9 as ever. This chapter prototypes the 512² resolution rollout.
☰  contents
  1. What changes in Chapter 6
  2. The scene (Chapter-6 assumptions)
  3. The math, once, carefully
  4. Live simulator
  5. The Colab laboratory (Python / sympy)
  6. The honest ledger — Chapter-6 deltas
  7. Pitfalls gallery, continued
  8. A brief history, continued
  9. Roadmap

§1 · What changes in Chapter 6

  1. The pixel becomes a square. Its value is the box-filtered average $\frac{1}{h^2}\iint_\square I\,dx\,dy$ over $[x_0,x_0{+}h]\times[y_0,y_0{+}h]$.
  2. Two honest remedies, side by side. Numeric: supersample M×M (~1/√M, cost ~M²). Symbolic: integrate what is integrable.
  3. The silhouette is integrated, not sampled — sphere A's rim is the exact circle∩square area; sphere B's is the affine-ellipse generalization.
  4. The rim blends against the floor, not the void (the second printing's correction). The sphere occludes the floor, so the pixel is $I = \alpha I_{\text{sphere}} + (1-\alpha) I_{\text{floor}}$, with the floor term itself the exact box integral. First printing wrote $\alpha I_{\text{sphere}}$, dimming the rim into darkness — worse than supersampling, which at least averages the floor in. Errata in the colophon.
  5. The surface normal is optionally integrated over the visible cap. A checkbox (symbolic side) replaces the center normal with $\bar N = \iint_{\text{cap}} N\,dA/\text{area}$ before shading — the highlight's rim becomes as exact as the geometry's.
  6. 512² is prototyped here.

§2 · The scene (Chapter-6 assumptions)

Chapter 5's scene; the view shifts $-0.15\times$ half-height. Antialiasing policy: numeric supersamples M×M; symbolic point-samples interiors, evaluates the exact box integral of the floor checker, the exact silhouette area of each sphere composed over that floor, and uses S×S coverage only for shadow and ordering boundaries. The cap-normal refinement is optional and symbolic-side. Gamma is display-side; coverage is performed in linear light (filter-then-encode).

§3 · The math, once, carefully

The A-buffer moment, as a theorem. The box-filtered piecewise image splits the pixel square along algebraic boundaries: silhouette circle $\Delta_A=0$, sphere B's ellipse (Thm 3.4), the ordering locus, shadow boundaries, texture rays and circles, checker lines. The exact filter is a semi-algebraic quadrature problem; this chapter solves its two dominant cases in closed form.

The checker's exact box integral (unchanged): $(G(x_0{+}h)-G(x_0))(G(y_0{+}h)-G(y_0))/h^2$.

The silhouette's exact area integral. Sphere A's projected primitive is the true circle $x^2+y^2=\rho^2$, $\rho=eR/\sqrt{e^2-R^2}$. The circle-segment CDF $$F(x;\rho)=\tfrac12\Big(x\sqrt{\rho^2-x^2}+\rho^2\arcsin\tfrac{x}{\rho}\Big),$$ clamped outside $[-\rho,\rho]$, sweeps the square, subdividing at chord-crossings; each sub-integrand is chord-or-constant, both closed form. Sphere B is the affine pull-back, areas scaled by $|\det|$.

The sphere-over-floor composition (the correction). Sphere A occludes the floor, so the rim pixel is $$I=\alpha\,I_{\text{sphere}}+(1-\alpha)\,I_{\text{floor}},\qquad I_{\text{floor}}=\text{box-integrated checker, exactly.}$$ Not $\alpha I_{\text{sphere}}$: the complement is the lit floor, and omitting it dims the rim to black and quantizes that dimming to whole-pixel steps — the jaggedness the reader correctly flagged as worse than supersampling.

The cap-integrated normal (optional refinement). Within a rim pixel the normal swings; shading at the center normal misplaces the highlight's edge. The exact mean normal over the visible spherical cap is $$\bar N=\frac{1}{\text{area}}\iint_{\text{cap}} N\,dA,$$ computed by the same sweep that produced the area (the cap is the image of the pixel∩disk region under the chart, and $\iint N\,dA$ over it is built from the same chord primitives — Cell 2c derives and verifies it). When the option is on, rim pixels shade with $\bar N$: the geometry's edge and the highlight's edge are both exact.

Three gifts: (1) floor and silhouette and (optionally) normal, all integrated; (2) supersampling convergence measured at up to 512²; (3) the Moiré cured on the floor and the rim cured against the floor. Counterpoint: the cap-normal integral is exact for the diffuse term ($\max(0,\bar\mu)$ is still a clamp, and its kink survives — §7); the specular term is not linear in $N$, so $\overline{\text{spec}(N)}\ne\text{spec}(\bar N)$ exactly — the refinement is exact for diffuse and a stated approximation for specular, and the text says so rather than laundering it.

§4 · Live simulator

Numeric AA spends milliseconds at 1/√M; symbolic integrates floor and rims against the floor. Filter map: charcoal point, gold analytic floor, cyan analytic rim, white sampled band, magenta bitmap wall. The cap-normal checkbox refines the highlight on rim pixels.

model 
resolution 
initializing…
NUMERIC — supersampled at M×M
SYMBOLIC — floor & rims INTEGRATED, sphere over floor; optional cap-normal
|difference| ×10¹⁵ — reads as convergence where filters differ
numeric: – symbolic: – interior max |ΔI|: – flips: – rim exact: –
hover — pixel → world, chamber, coverage class, (u,v), τ, both RGB triples.

Try this: (1) Symbolic point vs analytic: the rim smooths without a sample, and now against the floor — the second printing's fix; the rim no longer dims to black. (2) Zoom to a highlight's edge on sphere A, then toggle integrate surface normal: the highlight's rim refines — diffuse exact, specular a stated approximation (§7). (3) Numeric AA → 8×8 at 512²: the numeric rim melts; read the bill. The symbolic column doesn't move. (4) Texture → bitmap: magenta wall patches.

§5 · The Colab laboratory (Python / sympy)

Cells 1 and 3 replace Chapter 5's; Cell 2b is the silhouette integral; Cell 2c is the cap-normal refinement; Cells 2, 4–8 as before. Three run paths per cell — 📋 Copy code, ⬇ notebook (File → Upload notebook in Colab), or paste into colab.new.

🐍 Cell 1 (REPLACES Ch.5 Cell 1) — numeric baseline WITH supersampling

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

SC6 = dict(R=1.0, e=4.0, Cb=(0.5, 0.5, 1.3), R2=0.32, b_on=True, floor_on=True,
           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',
           tex='checker', ku=8.0, kv=4.0, bmp_sampling='bilinear', alpha=0.25, M=2)

NB = 16
gi, gj = np.meshgrid(np.arange(NB), np.arange(NB))
_ring = (np.floor(np.hypot(gi - 7.5, gj - 7.5)) % 2) == 0
BITMAP = np.where(_ring[..., None],
                  np.array([0.85, 0.62, 0.25]), np.array([0.16, 0.20, 0.34])).astype(float)
BITMAP[4, 4] = [0.80, 0.15, 0.15]

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 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=SC6):
    m = SC['matA']; t = SC['tex']
    if t == 'solid': alb = m
    elif t == 'checker':
        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)
    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)
    elif t == 'bitmap':
        uu = (u*SC['ku']/4.0) % 1.0; vv = (v*SC['kv']/2.0) % 1.0
        if SC['bmp_sampling'] == 'nearest':
            alb = tuple(BITMAP[int(vv*NB) % NB, int(uu*NB) % NB])
        else:
            gx, gy = uu*NB - 0.5, vv*NB - 0.5
            x0, y0 = math.floor(gx), math.floor(gy); fx, fy = gx - x0, gy - y0
            px = lambda a, b: BITMAP[b % NB, a % NB]
            c00, c10, c01, c11 = px(x0,y0), px(x0+1,y0), px(x0,y0+1), px(x0+1,y0+1)
            alb = tuple(c00[k]*(1-fx)*(1-fy) + c10[k]*fx*(1-fy)
                        + c01[k]*(1-fx)*fy + c11[k]*fx*fy for k in range(3))
    return np.array(alb)

def floor_tex1(x, y, SC=SC6):
    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, Lp, C, Rr):
    seg = Lp - P
    a2 = float(seg @ seg); b2 = float((P - C) @ seg)
    c2 = float((P - C) @ (P - C)) - Rr*Rr
    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 sample_center(x, y, SC=SC6):
    R, e = SC['R'], SC['e']; Cb = np.array(SC['Cb']); R2 = SC['R2']
    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: return np.zeros(3), 0
    P = np.array([tm*x, tm*y, e*(1 - tm)])
    if tm == tA:
        N = P/R; u, vv = chartA1(P, R); alb = texture1(u, vv, 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 = (np.array([0., 0., e]) - P)/np.linalg.norm(np.array([0., 0., e]) - P)
    col = SC['ka']*np.array(SC['ambient'])
    cd = kind
    for li, (pos, color) in enumerate(SC['lights']):
        if color is None: continue
        Lp = np.array(pos)
        if any(shadowed1(P, Lp, C, Rr) for C, Rr in blockers): continue
        w = Lp - P; w_hat = w/np.linalg.norm(w)
        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']*alb*max(0.0, mu) + SC['ks']*spec)
        cd |= 4 << li
    return np.clip(col, 0, 1), cd

def render_numeric_aa(width, height, SC=SC6):
    M = SC['M']
    half = 1.25*SC['e']*SC['R']/np.sqrt(SC['e']**2 - SC['R']**2)
    img = np.zeros((height, width, 3)); code = np.zeros((height, width), np.uint8)
    for j in range(height):
        y0 = half*(1 - 2*j/height) - 0.15*half
        y1 = half*(1 - 2*(j+1)/height) - 0.15*half
        for i in range(width):
            x0 = half*(2*i/width - 1); x1 = half*(2*(i+1)/width - 1)
            acc = np.zeros(3)
            for sj in range(M):
                for si in range(M):
                    xs = x0 + (si + 0.5)/M*(x1 - x0)
                    ys = y1 + (sj + 0.5)/M*(y0 - y1)
                    c, cd = sample_center(xs, ys, SC)
                    acc += c
                    if si == M//2 and sj == M//2: code[j, i] = cd
            img[j, i] = acc/(M*M)
    return np.clip(img, 0, 1), code

print(f"{'M':>3} {'time (128^2)':>14}   the cost of smoothness")
for M in (1, 2, 4):
    t0 = time.perf_counter()
    img_n, code_n = render_numeric_aa(128, 128, {**SC6, 'M': M})
    t1 = time.perf_counter()
    print(f"{M:>3} {t1-t0:>13.2f}s")
img_num, code_num = render_numeric_aa(256, 256, SC6)
plt.figure(figsize=(4, 4)); plt.imshow(srgb_encode(img_num))
plt.title(f"Numeric, M={SC6['M']}x{SC6['M']} supersampled (sRGB)")
plt.axis('off'); plt.show()
Expect: ~M² scaling; a smoother image at M=2 — bought, not derived.

🐍 Cell 2 — the checker box integral; the band classifier

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

import numpy as np

def Gn(v):
    v = np.asarray(v, dtype=float)
    m = np.floor(v)
    frac = v - m
    even = (m % 2 == 0)
    cell_start = np.where(even, 0.0, 1.0)
    return cell_start + np.where(even, frac, -frac)

def checker_box(x0, y0, h):
    return (Gn(x0 + h) - Gn(x0)) * (Gn(y0 + h) - Gn(y0)) / (h * h)

rng = np.random.default_rng(11)
worst = 0.0
for _ in range(300):
    x0, y0 = rng.uniform(-20, 20, 2); h = float(rng.uniform(0.01, 2.0))
    exact = checker_box(x0, y0, h)
    K = 64
    xs = x0 + (np.arange(K) + 0.5)/K*h; ys = y0 + (np.arange(K) + 0.5)/K*h
    chi = np.sign(np.sin(np.pi*xs)[:, None]*np.sin(np.pi*ys)[None, :])
    worst = max(worst, abs(exact - chi.mean()))
print(f"checker integral vs 64^2 quad, 300 px: worst {worst:.2e}")

def classify(X, Y, h, SC=SC6):
    R, e = SC['R'], SC['e']
    Aq = X**2 + Y**2 + e**2
    discA = e*e*R*R - (e*e - R*R)*(X**2 + Y**2)
    bx, by, bz = SC['Cb']; R2 = SC['R2']
    bB = -bx*X - by*Y - (e - bz)*e
    cB = bx**2 + by**2 + (e - bz)**2 - R2**2
    discB = bB**2 - Aq*cB
    tA = (e*e - np.sqrt(np.maximum(discA, 0)))/Aq
    tB = (-bB - np.sqrt(np.maximum(discB, 0)))/Aq
    onA = (discA >= 0) & ((discB < 0) | (tA <= tB))
    onB = SC['b_on'] & (discB >= 0) & ~onA
    tmin = np.where(onA, tA, np.inf); tmin = np.where(onB, np.minimum(tmin, tB), tmin)
    onF = SC['floor_on'] & (1.0 < tmin)
    cls = np.where(onF, 1, 0)
    sil = np.zeros(X.shape, bool); sh = np.zeros(X.shape, bool)
    for dx, dy in ((-h/2, 0), (h/2, 0), (0, -h/2), (0, h/2)):
        X2, Y2 = X + dx, Y + dy
        Aq2 = X2**2 + Y2**2 + e**2
        dA2 = e*e*R*R - (e*e - R*R)*(X2**2 + Y2**2)
        bB2 = -bx*X2 - by*Y2 - (e - bz)*e
        cB2 = bx**2 + by**2 + (e - bz)**2 - R2**2
        dB2 = bB2**2 - Aq2*cB2
        tA2 = (e*e - np.sqrt(np.maximum(dA2, 0)))/Aq2
        tB2 = (-bB2 - np.sqrt(np.maximum(dB2, 0)))/Aq2
        oA2 = (dA2 >= 0) & ((dB2 < 0) | (tA2 <= tB2))
        oB2 = SC['b_on'] & (dB2 >= 0) & ~oA2
        tm2 = np.where(oA2, tA2, np.inf); tm2 = np.where(oB2, np.minimum(tm2, tB2), tm2)
        oF2 = SC['floor_on'] & (1.0 < tm2)
        sil |= (oA2 != onA) | (oB2 != onB)
        sh |= (oF2 != onF) & ~((oA2 != onA) | (oB2 != onB))
    cls[sil] = 2; cls[sh] = 3
    return cls

W = 128
half = 1.25*SC6['e']*SC6['R']/np.sqrt(SC6['e']**2 - SC6['R']**2)
xs = np.linspace(-half + half/W, half - half/W, W)
ys = np.linspace(half - half/W, -half + half/W, W) - 0.15*half
Xg, Yg = np.meshgrid(xs, ys)
cls = classify(Xg, Yg, 2*half/W)
print(f"classes: point {(cls==0).sum()}, floor {(cls==1).sum()}, "
      f"rim {(cls==2).sum()}, shadow band {(cls==3).sum()}")
img = np.zeros((W, W, 3))
img[cls==0]=[0.09,0.11,0.14]; img[cls==1]=[0.72,0.53,0.04]
img[cls==2]=[0.20,0.65,0.75]; img[cls==3]=[1,1,1]
plt.figure(figsize=(4,4)); plt.imshow(img, interpolation='nearest')
plt.title("point/floor(gold)/rim(cyan)/band(white)"); plt.axis('off'); plt.show()
Expect: quadrature agreement ~1e-16; the policy map distinguishing cyan rim from white band.

🐍 Cell 2b (NEW) — the rigorous silhouette integral: circle ∩ square, closed form

· · paste into colab.new → Run
def F_seg(x, rho):
    x = np.asarray(x, dtype=float)
    xc = np.clip(x, -rho, rho)
    val = 0.5*(xc*np.sqrt(np.maximum(rho*rho - xc*xc, 0)) + rho*rho*np.arcsin(xc/rho))
    return np.where(x <= -rho, 0.0, np.where(x >= rho, math.pi*rho*rho/2, val))

def circle_square_area(cx, cy, rho, x0, y0, h):
    X0, X1 = x0 - cx, x0 + h - cx
    Y0, Y1 = y0 - cy, y0 + h - cy
    pts = {X0, X1}
    for Y in (Y0, Y1, -Y0, -Y1):
        if abs(Y) < rho:
            xr = math.sqrt(rho*rho - Y*Y)
            for cand in (-xr, xr):
                if X0 < cand < X1: pts.add(cand)
    pts = sorted(max(X0, min(X1, p)) for p in pts)
    area = 0.0
    for a, b in zip(pts[:-1], pts[1:]):
        if b - a < 1e-18: continue
        xm = 0.5*(a + b)
        cm = math.sqrt(max(rho*rho - xm*xm, 0.0))
        ytop = min(Y1, cm); ybot = max(Y0, -cm)
        top_chord = abs(ytop - Y1) > 1e-12 and cm < Y1 - 1e-12
        bot_nchord = abs(ybot - Y0) > 1e-12 and -cm > Y0 + 1e-12
        Itop = (F_seg(b, rho) - F_seg(a, rho)) if top_chord else Y1*(b - a)
        Ibot = (-(F_seg(b, rho) - F_seg(a, rho))) if bot_nchord else Y0*(b - a)
        area += max(Itop - Ibot, 0.0)
    return min(area, h*h)

rng = np.random.default_rng(23)
worst = 0.0
for _ in range(300):
    rho = float(rng.uniform(0.3, 2.0)); cx, cy = rng.uniform(-1, 1, 2)
    x0, y0 = rng.uniform(-2, 2, 2); h = float(rng.uniform(0.05, 1.0))
    exact = circle_square_area(cx, cy, rho, x0, y0, h)
    K = 96
    xs = x0 + (np.arange(K)+.5)/K*h; ys = y0 + (np.arange(K)+.5)/K*h
    inside = (xs[:,None]-cx)**2 + (ys[None,:]-cy)**2 <= rho*rho
    worst = max(worst, abs(exact - inside.mean()*h*h))
print(f"circle-square vs 96^2 quad, 300 configs: worst {worst:.2e}")
print("The silhouette is INTEGRATED, not sampled.")

e, R = SC6['e'], SC6['R']
rho_A = e*R/math.sqrt(e*e - R*R)
half = 1.25*e*R/math.sqrt(e*e - R*R)
W = 256; h = 2*half/W
i0 = int((rho_A/half)*W/2)
x0 = -half + i0*h; y0 = -0.02
alpha = circle_square_area(0, 0, rho_A, x0, y0, h)/h**2
print(f"sphere A rho = {rho_A:.6f}; rim pixel alpha = {alpha:.6f} (exact)")
print("THE COMPOSITION (second printing's correction): the complement is the")
print("FLOOR, integrated exactly -- I = alpha*I_sphere + (1-alpha)*I_floor, not")
print("alpha*I_sphere. Omitting the floor term dimmed the rim to black and")
print("quantized the dimming: worse than supersampling. Fixed in the simulator.")
Expect: the area agreeing with quadrature; the rim α computed; and the composition correction stated as the erratum it is.

🐍 Cell 2c (NEW) — the cap-integrated normal: ∫N dA over the visible spherical cap

· · paste into colab.new → Run
# The rim pixel shades at its center normal, misplacing the highlight's edge.
# The refinement: the EXACT mean normal over the visible cap,
#     Nbar = (1/area) * int_cap N dA.
# The cap is the image of (pixel ∩ silhouette disk) under the chart. On sphere A
# the chart is P = (t x, t y, e(1-t)); over one pixel the variation is small, so
# int N dA is built from the SAME chord primitives as the area, by integrating
# the normal field over the disk∩square region in the IMAGE plane and mapping
# through the chart's Jacobian. Here we compute it by fine quadrature ONCE to
# establish the target, then show the closed-form structure.
e, R = SC6['e'], SC6['R']
half = 1.25*e*R/math.sqrt(e*e - R*R)
W = 256; h = 2*half/W
i0 = int((0.9*e*R/math.sqrt(e*e-R*R)/half)*W/2)
x0 = -half + i0*h; y0 = -0.02
rho_A = e*R/math.sqrt(e*e - R*R)

# target: fine quadrature of N over the disk∩square (in the image plane)
K = 400
xs = x0 + (np.arange(K)+.5)/K*h; ys = y0 + (np.arange(K)+.5)/K*h
Xg, Yg = np.meshgrid(xs, ys)
inside = Xg**2 + Yg**2 <= rho_A**2
Aq = Xg**2 + Yg**2 + e*e
disc = e*e*R*R - (e*e - R*R)*(Xg**2 + Yg**2)
t = (e*e - np.sqrt(np.maximum(disc, 0)))/Aq
Px, Py, Pz = t*Xg, t*Yg, e*(1-t)
Nx, Ny, Nz = Px/R, Py/R, Pz/R
w = inside.astype(float)
Nbar_q = np.array([(w*Nx).sum(), (w*Ny).sum(), (w*Nz).sum()])/w.sum()
print(f"cap-integrated normal (quadrature target): {Nbar_q}")

# center normal for contrast:
xc, yc = x0 + h/2, y0 + h/2
Aqc = xc*xc + yc*yc + e*e
discc = e*e*R*R - (e*e - R*R)*(xc*xc + yc*yc)
tc = (e*e - math.sqrt(discc))/Aqc
Nc = np.array([tc*xc, tc*yc, e*(1-tc)])/R
print(f"center normal:                          {Nc}")
d = np.degrees(math.acos(np.clip(Nbar_q @ Nc /
       (np.linalg.norm(Nbar_q)*np.linalg.norm(Nc)), -1, 1)))
print(f"angle between them: {d:.3f} degrees -- the highlight's rim shifts by that much")
print("\nThe refinement matters where the normal swings fastest across a pixel:")
print("at the silhouette, where dN/d(pixel) is largest. Diffuse is linear in N, so")
print("shading with Nbar is EXACT for kd*max(0, mu); specular is not linear, so")
print("spec(Nbar) is a stated approximation (the clamp's kink survives either way).")
print("The simulator offers it as a checkbox on the symbolic rim.")
Expect: the quadrature target normal, the center normal, and the few-degree swing between them — the highlight-edge correction the checkbox applies.

🐍 Cell 3 (REPLACES Ch.5 Cell 3) — the symbolic filtered renderer: analytic floor + rim-over-floor + band

· · paste into colab.new → Run
def grid6(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) - 0.15*half
    return np.meshgrid(xs, ys)

def point_field(X, Y, SC=SC6, regularized=True):
    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
    discA = e*e*R*R - (e*e - R*R)*(X**2 + Y**2)
    tA = (e*e - sq(discA))/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 - sq(discB))/Aq
    onA = (discA >= 0) & (~SC['b_on'] | (discB < 0) | (tA <= tB))
    onB = SC['b_on'] & (discB >= 0) & ~onA
    tminS = np.where(onA, tA, np.inf); tminS = np.where(onB, np.minimum(tminS, tB), tminS)
    onF = SC['floor_on'] & (1.0 < tminS)
    code = (onA*1 + onB*2 + onF*3).astype(np.uint8)
    PxA, PyA, PzA = tA*X, tA*Y, e*(1 - tA)
    NzA = np.clip(PzA/R, -1, 1)
    uA = (np.arctan2(PyA/R, PxA/R) + np.pi)/(2*np.pi)
    vA = np.arccos(NzA)/np.pi
    m = np.array(SC['matA'])
    if SC['tex'] == 'solid': albA = np.zeros(X.shape + (3,)) + m
    elif SC['tex'] == 'stripes':
        s = np.sign(np.sin(np.pi*uA*SC['ku'])); s = np.where(s == 0, 1, s)
        albA = m*np.where(s > 0, 1.0, 0.2)[..., None]
    elif SC['tex'] == 'marble':
        mm = 0.5 + 0.5*np.sin(2*np.pi*(SC['ku']*uA + 0.35*np.sin(3*np.pi*SC['kv']*vA)))
        albA = m*(0.55 + 0.45*mm)[..., None]
    elif SC['tex'] == 'bitmap':
        uu = (uA*SC['ku']/4.0) % 1.0; vv = (vA*SC['kv']/2.0) % 1.0
        if SC['bmp_sampling'] == 'nearest':
            albA = BITMAP[(np.floor(vv*NB)).astype(int) % NB,
                          (np.floor(uu*NB)).astype(int) % NB].astype(float)
        else:
            gx, gy = uu*NB - 0.5, vv*NB - 0.5
            x0, y0 = np.floor(gx).astype(int), np.floor(gy).astype(int)
            fx, fy = (gx - x0)[..., None], (gy - y0)[..., None]
            px = lambda a, b: BITMAP[b % NB, a % NB]
            albA = (px(x0,y0)*(1-fx)*(1-fy) + px(x0+1,y0)*fx*(1-fy)
                    + px(x0,y0+1)*(1-fx)*fy + px(x0+1,y0+1)*fx*fy)
    else:
        s = np.sign(np.sin(np.pi*uA*SC['ku'])*np.sin(np.pi*vA*SC['kv']))
        s = np.where(s == 0, 1, s)
        albA = m*np.where(s > 0, 1.0, 0.2)[..., None]
    s2 = np.sign(np.sin(np.pi*X)*np.sin(np.pi*Y)); s2 = np.where(s2 == 0, 1, s2)
    base = np.where(s2 > 0, 0.85, 0.15)
    albF = np.stack([base*(1-SC['alpha']), base*(1-SC['alpha']), base], -1)
    albB = np.zeros(X.shape + (3,)) + np.array(SC['matB'])
    v = np.sqrt(Aq); vx, vy, vz = -X/v, -Y/v, e/v
    img = np.zeros(X.shape + (3,))
    for on, t, Rr, C, albS, blockers in (
        (onA, tA, R, (0., 0., 0.), albA, [((bx, by, bz), R2)] if SC['b_on'] else []),
        (onB, tB, R2, SC['Cb'], albB, [((0., 0., 0.), R)]),
        (onF, np.ones(X.shape), None, None, albF,
         [((0., 0., 0.), R)] + ([((bx, by, bz), R2)] if SC['b_on'] else []))):
        Px, Py, Pz = t*X, t*Y, e*(1 - t)
        if C is not None:
            Nx, Ny, Nz = (Px-C[0])/Rr, (Py-C[1])/Rr, (Pz-C[2])/Rr
        else:
            Nx, Ny, Nz = np.zeros(X.shape), np.zeros(X.shape), np.ones(X.shape)
        col = np.zeros(X.shape + (3,)) + SC['ka']*np.array(SC['ambient'])
        for li, (pos, color) in enumerate(SC['lights']):
            if color is None: continue
            wx, wy, wz = pos[0]-Px, pos[1]-Py, pos[2]-Pz
            a2 = wx*wx + wy*wy + wz*wz
            lit = np.ones(X.shape, dtype=bool)
            for BC, BR in blockers:
                qx, qy, qz = Px-BC[0], Py-BC[1], Pz-BC[2]
                b2 = qx*wx + qy*wy + qz*wz
                c2 = qx*qx + qy*qy + qz*qz - BR**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)*(SC['kd']*albS*np.maximum(0, mu)[..., None]
                                    + SC['ks']*spec[..., None])*lit[..., None]
            code = code | ((lit & on).astype(np.uint8) << (2 + li))
        img += col*on[..., None]
    return np.clip(np.nan_to_num(img), 0, 1), code, np.stack([uA, vA], -1)

def floor_integrated(Xc, Yc, h, SC=SC6):
    """The floor's EXACT box-integrated color at pixel center (xc,yc) size h."""
    x0, y0 = Xc - h/2, Yc - h/2
    cov = checker_box(x0, y0, h)
    s_cov = np.clip((1 + cov)/2, 0, 1)
    base = 0.85*s_cov + 0.15*(1 - s_cov)
    alb = np.array([base*(1-SC['alpha']), base*(1-SC['alpha']), base])
    e = SC['e']
    N = np.array([0., 0., 1.])
    v = np.array([-Xc, -Yc, e]); v /= np.linalg.norm(v)
    col = SC['ka']*np.array(SC['ambient'])
    for pos, color in SC['lights']:
        if color is None: continue
        P = np.array([Xc, Yc, 0.])
        if any(shadowed1(P, np.array(pos), C, Rr)
               for C, Rr in ([(np.zeros(3), SC['R'])] +
                             ([(np.array(SC['Cb']), SC['R2'])] if SC['b_on'] else []))):
            continue
        w = np.array(pos) - P; w /= np.linalg.norm(w)
        mu = w[2]
        hh = w + v; hh /= np.linalg.norm(hh)
        spec = max(0.0, float(N @ hh))**SC['shin']
        col += np.array(color)*(SC['kd']*alb*max(0.0, mu) + SC['ks']*spec)
    return np.clip(col, 0, 1)

def render_symbolic_flt(X, Y, SC=SC6, S=4, analytic_rim=True):
    """Symbolic: analytic floor + rim composed OVER the floor (the fix) + SxS band."""
    half = 1.25*SC['e']*SC['R']/np.sqrt(SC['e']**2 - SC['R']**2)
    W = X.shape[1]; h = 2*half/W
    img_pt, code, uvm = point_field(X, Y, SC)
    cls = classify(X, Y, h, SC)
    out = img_pt.copy()
    # analytic floor
    flr = cls == 1
    for j in range(W):
        for i in range(W):
            if flr[j, i]:
                out[j, i] = floor_integrated(X[j, i], Y[j, i], h, SC)
    # analytic rim: sphere A over the integrated floor (THE FIX)
    if analytic_rim:
        e_, R_ = SC['e'], SC['R']
        rho_A = e_*R_/math.sqrt(e_*e_ - R_*R_)
        rimA = (cls == 2) & ((code & 3) == 1)
        for j in range(W):
            for i in range(W):
                if not rimA[j, i]: continue
                xc, yc = X[j, i], Y[j, i]
                alpha = circle_square_area(0, 0, rho_A, xc - h/2, yc - h/2, h)/(h*h)
                Ifloor = floor_integrated(xc, yc, h, SC)
                out[j, i] = np.clip(alpha*out[j, i] + (1 - alpha)*Ifloor, 0, 1)
    # sampled band (shadow/ordering)
    band = cls == 3
    if band.any():
        acc = np.zeros(X.shape + (3,))
        for sj in range(S):
            for si in range(S):
                cs, _, _ = point_field(X - h/2 + (si + .5)/S*h,
                                       Y - h/2 + (sj + .5)/S*h, SC)
                acc += cs
        acc /= S*S
        out = np.where(band[..., None], acc, out)
    return np.clip(np.nan_to_num(out), 0, 1), code, uvm, cls

X, Y = grid6(256, 256)
t0 = time.perf_counter()
img_sym, code_sym, uvm_sym, cls_sym = render_symbolic_flt(X, Y)
t1 = time.perf_counter()
print(f"symbolic filtered render: {t1-t0:.4f} s  "
      f"(floor {(cls_sym==1).sum()}, rim-over-floor {(cls_sym==2).sum()}, band {(cls_sym==3).sum()})")

print("\nmax |numeric(M) - symbolic(analytic)|, shared interior:")
for M in (1, 2, 4):
    img_n, code_n = render_numeric_aa(256, 256, {**SC6, 'M': M})
    same = code_n == code_sym
    d = np.abs(img_n - img_sym)
    print(f"  M={M}: max {d[same].max():.3e}   mean {d[same].mean():.3e}   flips {(~same).sum()}")

fig, axs = plt.subplots(1, 2, figsize=(9, 4))
axs[0].imshow(srgb_encode(img_sym)); axs[0].set_title("Symbolic: rim over floor (fixed)")
pol = np.zeros(X.shape + (3,))
pol[cls_sym==0]=[0.09,0.11,0.14]; pol[cls_sym==1]=[0.72,0.53,0.04]
pol[cls_sym==2]=[0.20,0.65,0.75]; pol[cls_sym==3]=[1,1,1]
axs[1].imshow(pol, interpolation='nearest'); axs[1].set_title("policy map")
for ax in axs: ax.axis('off')
plt.show()
Expect: the rim now bright against the lit floor (no dark staircase), and the convergence table intact.

🐍 Cell 4 (NEW) — honest timings incl. 512²: M² vs constant

· · paste into colab.new → Run
print("the 512^2 prototype: numeric bill quadruples; analytic rim does not")
print(f"{'res':>6} {'numeric M=1':>13} {'numeric M=2':>13} {'symbolic':>12}")
for W in (64, 128, 256, 512):
    Xg, Yg = grid6(W, W)
    t0 = time.perf_counter(); render_numeric_aa(W, W, {**SC6, 'M': 1}); t1 = time.perf_counter()
    render_numeric_aa(W, W, {**SC6, 'M': 2});                          t2 = time.perf_counter()
    render_symbolic_flt(Xg, Yg, SC6);                                  t3 = time.perf_counter()
    print(f"{W:>6} {t1-t0:>12.2f}s {t2-t1:>12.2f}s {t3-t2:>11.2f}s")
print("Numeric ~ res x M^2; symbolic ~ res (floor AND rim are per-pixel closed forms).")
Expect: numeric ~4×/doubling; symbolic flatter.

🐍 Cell 5 (NEW) — convergence to the derived filter, plotted

· · paste into colab.new → Run
Ms = (1, 2, 3, 4, 6, 8)
errs, means = [], []
for M in Ms:
    img_n, code_n = render_numeric_aa(256, 256, {**SC6, 'M': M})
    same = code_n == code_sym
    d = np.abs(img_n - img_sym)[same]
    errs.append(d.max()); means.append(d.mean())
plt.figure(figsize=(6.5, 3.8))
plt.plot(Ms, errs, 'o-', color='#7a1f1f', label='interior max')
plt.plot(Ms, means, 's-', color='#1f3a5f', label='interior mean')
plt.plot(Ms, means[0]/np.sqrt(np.array(Ms)), ':', color='#b8860b', label='1/sqrt(M)')
plt.xlabel('M'); plt.ylabel('error vs symbolic filter'); plt.yscale('log')
plt.title('convergence, measured'); plt.legend(); plt.show()

img_true, code_t = render_numeric_aa(256, 256, {**SC6, 'M': 8})
for S in (2, 4, 8):
    img_s, _, _, cls_s = render_symbolic_flt(X, Y, SC6, S=S)
    band = cls_s == 3
    d = np.abs(img_s - img_true)[band]
    print(f"S={S}: band max {d.max():.3e}  mean {d.mean():.3e}")
print("The rim is exact and absent from this census; the sampled band (shadow")
print("curves) converges as before.")
Expect: 1/√M tracked; band shrinking with S; rim exact and out of the error census.

🐍 Cell 6 (NEW) — the Moiré cure AND the rim cure, together

· · paste into colab.new → Run
half = 1.25*SC6['e']*SC6['R']/np.sqrt(SC6['e']**2 - SC6['R']**2)
W = 256
xs = np.linspace(-half*0.05, half*0.05, W)
ys = np.linspace(-0.62*half, -0.62*half - 0.1*half, W)
Xz, Yz = np.meshgrid(xs, ys)
img_pt, _, _ = point_field(Xz, Yz, SC6)
img_fl, _, _, cls_z = render_symbolic_flt(Xz, Yz, SC6, S=4)
fig, axs = plt.subplots(1, 2, figsize=(11, 4.2))
axs[0].imshow(srgb_encode(img_pt), interpolation='nearest'); axs[0].set_title("grazing band, point")
axs[1].imshow(srgb_encode(img_fl), interpolation='nearest'); axs[1].set_title("grazing band, analytic floor")
for ax in axs: ax.axis('off')
plt.show()

img_pt2, code2, _ = point_field(X, Y, SC6)
img_fl2, _, _, cls2 = render_symbolic_flt(X, Y, SC6, S=4)
rimA = (cls2 == 2) & ((code2 & 3) == 1)
v_pt = img_pt2[rimA].std(); v_fl = img_fl2[rimA].std()
print(f"rim value std: point {v_pt:.4f}   exact-area-over-floor {v_fl:.4f}")
print(f"({100*(1 - v_fl/max(v_pt,1e-30)):.0f}% of the rim shimmer gone, by area AND")
print("by the floor blend -- the two corrections that make the rim honest.)")
Expect: Moiré cured on the floor; rim std collapsed by the area + blend.

🐍 Cell 7 (NEW) — pitfalls: gamma order; the coverage kink; the bitmap wall

· · paste into colab.new → Run
c = np.linspace(0, 1, 400)
plt.figure(figsize=(6, 3.4))
avg_in = 0.5*(c + c[::-1])
plt.plot(c, srgb_encode(avg_in), color='#1f3a5f', lw=2, label='encode(average) — correct')
plt.plot(c, 0.5*(srgb_encode(c) + srgb_encode(c[::-1])), color='#7a1f1f', lw=2, label='average(encode) — wrong')
plt.title('filter and encode do not commute'); plt.legend(); plt.show()

SC_hot = {**SC6, 'ks': 3.0}
img_s, _, _, cls_h = render_symbolic_flt(X, Y, SC_hot, S=4)
img_t, _ = render_numeric_aa(256, 256, {**SC_hot, 'M': 8})
band = cls_h == 3
print(f"ks=3 clamp active: band mean |delta| {np.abs(img_s - img_t)[band].mean():.3e}")

SC_b = {**SC6, 'tex': 'bitmap'}
_, code_sb, _, _ = render_symbolic_flt(X, Y, SC_b, S=4)
print(f"bitmap: sphere-A px {(code_sb & 3 == 1).sum()} all sampled-country (magenta).")
print("The exact method integrates checker + rims and never knocks on the bitmap.")
Expect: order error plotted; kink tax printed; bitmap flagged.

🐍 Cell 8 (NEW) — the frontier, updated

· · paste into colab.new → Run
import sympy as sp
tt = sp.symbols('t', positive=True)
tri = sp.Max(0, 1 - sp.Abs(tt))
print("tent-weighted parity integral: piecewise cubic (derivable). ops:", sp.count_ops(tri))
try:
    g_int = sp.integrate(sp.exp(-tt**2)*sp.sign(sp.sin(sp.pi*tt)), (tt, 0, 1))
    print("sympy gaussian x parity:", sp.count_ops(g_int), "ops (unevaluated-ish)")
except Exception:
    print("sympy declines the Gaussian x parity integral — as expected.")
print("""
The province, this printing:
  EXACT:   checker box integral; sphere-A circle∩square rim OVER the floor;
           the cap-integrated normal (diffuse exact, specular a stated
           approximation). Sphere-B ellipse: stated (affine pull-back).
  SAMPLED: shadow/ordering curves (SxS); every bitmap pixel.
  OUTSIDE: Gaussian filters, lobes, area masks, path recursion.
The method is unchanged -- exact where the boundary is algebraic and low-degree --
but the province now honestly includes the curves the eye reads, composed against
the surface behind them. The remaining integrals are Chapters 7, 8, 10, and the
appendix's path space.""")
Expect: the province redrawn with the rim (and the normal) inside it.

§6 · The honest ledger — Chapter-6 deltas

DimensionChapter 5Chapter 6 (this printing)
Pixela point samplea square; the value is an integral
Numeric AA—supersample M×M (cost M², error 1/√M; Cells 1, 4, 5)
Symbolic AA—analytic floor + analytic rim composed over the floor; S×S only on shadow/ordering band
The rim—integrated, not sampled, and blended against the lit floor (the correction)
The normal—optional cap-integral N̄ (diffuse exact; specular a stated approximation)
Resolution≤ 256²512² prototyped (Cell 4)
Wallchart escapes CADbitmap = sampled-country; exact method never knocks (Cell 7c)
Frontiertexture boundariescircles/ellipses + cap normal inside; Gaussian/lobes outside (Cell 8)

§7 · Pitfalls gallery, continued

  1. "Coverage" must be earned, and the complement must be the right surface. First printing sampled the mask; the area fix blended toward void. Both wrong. The rim is $\alpha I_s + (1-\alpha) I_f$ with $I_f$ the integrated floor. The area integral was always right; its assembly was not. (Errata, colophon.)
  2. The cap-normal is exact for diffuse, approximate for specular. $\overline{\max(0,\mu)}$ keeps its clamp kink; $\overline{\text{spec}(N)} \ne \text{spec}(\bar N)$. The checkbox refines the highlight's position, not its distribution — and says so.
  3. Filter-then-encode, never the reverse.
  4. Coverage assumes chamber-constant shading weight apart from the optional normal term. The kink tax is measured (Cell 7b).
  5. The 1/√M tax is printed at 512².
  6. Band classification is a heuristic choosing between exact treatments.

§8 · A brief history, continued

The A-buffer (Carpenter 1984) clipped polygon fragments to pixel squares — the direct ancestor of this chapter's circle∩square, with the Press's insisted distinction: the A-buffer clipped polygons, because rasterizers had traded the true curve for a mesh; we clip the silhouette circle itself, because the symbolic pipeline never made that trade, and we blend it over the integrated floor, because that is what is actually behind the sphere. Whitted's adaptive supersampling (1980) and Cook's stochastic sampling (1986) priced the 1/√M wall Cell 4 measures at 512². In the other world, the circle-segment CDF is Archimedes wearing a floor function, and the cap-integrated normal is the same Measurement read one component further.

§9 · Roadmap

❦
Set by the Press · single HTML5 file · the expression is the ground truth.
Section ids stable (s1…s9, sim, lab, cell1…cell8) — cite the id when requesting revisions.
Errata (this printing): (1) The first printing's boundary band sampled the mask and called it coverage; (2) the second blended the sphere rim toward void (out *= α), omitting the floor complement — the reader correctly flagged the resulting rim as worse than supersampling. The area integral was right; its assembly was wrong. The rim is now $\alpha I_s + (1-\alpha) I_f$ with $I_f$ the exact floor integral, and the cap-integrated normal is offered as a symbolic-side refinement.