Chapter 5 — Texture, or Albedo Becomes a Function

import numpy as np
import matplotlib.pyplot as plt
import math, time

SC5 = 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',
           legacy=False, alpha=0.25)

# --- the shared procedural bitmap: IDENTICAL formula in the JS simulator ------
NB = 16
ii, jj = np.meshgrid(np.arange(NB), np.arange(NB))
_ring = (np.floor(np.hypot(ii - 7.5, jj - 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]      # orientation marker dot

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)

# --- numeric realizations: scalar, lazily evaluated (runtime branches) --------
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=SC5):
    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))
    if SC['legacy']:
        alb = tuple(srgb_encode(alb))          # wrong on purpose: encode BEFORE lighting
    return np.array(alb)

def floor_tex1(x, y, SC=SC5):
    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
    alb = (base*(1-SC['alpha']), base*(1-SC['alpha']), base)
    return np.array(srgb_encode(alb)) if SC['legacy'] else np.array(alb)

def shadowed1(P, Lp, C, Rr):
    """Shadow feeler, EXPLICIT roots (numeric formulation)."""
    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 render_numeric_tex(width, height, SC=SC5):
    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)
    uvm = np.full((height, width, 2), np.nan)
    alb_out = np.zeros((height, width, 3))
    half = 1.25*e*R/np.sqrt(e**2 - R**2)
    E = np.array([0.0, 0.0, e])
    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; u, vv = chartA1(P, R)
                alb = texture1(u, vv, SC)
                uvm[j, i] = (u, vv)
                blockers = [(Cb, R2)] if SC['b_on'] else []
                kind = 1
            elif tm == tB:
                N = (P - Cb)/R2; alb = np.array(SC['matB'])
                if SC['legacy']: alb = srgb_encode(alb)
                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)
            amb = np.array(SC['ambient'])
            if SC['legacy']: amb = srgb_encode(amb)
            col = SC['ka']*amb
            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
            img[j, i] = col; code[j, i] = cd; alb_out[j, i] = alb
    return np.clip(img, 0, 1), code, uvm, alb_out

t0 = time.perf_counter()
img_num, code_num, uvm_num, alb_num = render_numeric_tex(256, 256)
t1 = time.perf_counter()
print(f"numeric textured render: {t1 - t0:.3f} s")
fig, axs = plt.subplots(1, 2, figsize=(9, 4))
axs[0].imshow(srgb_encode(img_num)); axs[0].set_title("Numeric, sRGB display")
axs[1].imshow(img_num); axs[1].set_title("Numeric, linear raw")
for ax in axs: ax.axis('off')
plt.show()

import sympy as sp

x, y = sp.symbols('x y', real=True)
R, e = sp.symbols('R e', positive=True)
Aq = x**2 + y**2 + e**2
disc = e**2*R**2 - (e**2 - R**2)*(x**2 + y**2)
tA = (e**2 - sp.sqrt(disc))/Aq
Px, Py, Pz = tA*x, tA*y, e*(1 - tA)
Nx, Ny, Nz = Px/R, Py/R, Pz/R

# --- THEOREM: atan2(Ny, Nx) = atan2(y, x) exactly (positive common scaling) ----
lhs = sp.atan2(Ny, Nx)
rhs = sp.atan2(y, x)
num_check = []
rng = np.random.default_rng(3)
for _ in range(20000):
    xv, yv = rng.normal(size=2)*0.8
    dA = 16 - 15*(xv*xv + yv*yv)
    if dA < 0: continue
    tv = (16 - np.sqrt(dA))/(xv*xv + yv*yv + 16)
    num_check.append(abs(np.arctan2(tv*yv, tv*xv) - np.arctan2(yv, xv)) < 1e-12)
print(f"longitude == pixel polar angle: {sum(num_check)}/{len(num_check)} ✓ (exact)")

# --- COROLLARY: v depends only on r^2 ------------------------------------------
r2 = sp.symbols('r2', positive=True)
Nz_r = sp.simplify(Nz.xreplace({x**2: r2 - y**2}))
print("Nz as a function of r^2 only:", Nz_r)
print("=> constant-v curves are CIRCLES, constant-u curves are RAYS.")
print("=> seam (Ny=0, Nx<0) renders as the slit y=0, x<0.")
print("=> pole (Nz=1) renders at (x,y)=(0,0): the pinwheel at window center.")

# --- magnification: the chart Jacobian at pole vs mid-latitude ------------------
u_expr = (sp.atan2(Ny, Nx) + sp.pi)/(2*sp.pi)
v_expr = sp.acos(Nz)/sp.pi
print(f"\nops: u = {sp.count_ops(u_expr)}, v = {sp.count_ops(v_expr)}")
print("texture composition adds:", sp.count_ops(u_expr) + sp.count_ops(v_expr),
      "ops before tau itself")
print("At the pole the u-gradient is unbounded (angular coordinate); at the rim")
print("the v-gradient blows up (the sphere turns edge-on). Both are visible as")
print("texture squeezing -- and both are computed, not asserted, in Cell 7.")

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) - 0.15*half
    return np.meshgrid(xs, ys)

# --- symbolic realization: TOTAL, vectorized; evaluated everywhere, selected ---
def textureV(u, v, SC=SC5):
    m = np.array(SC['matA']); t = SC['tex']
    if t == 'solid':
        alb = np.zeros(u.shape + (3,)) + m
    elif t == 'checker':
        s = np.sign(np.sin(np.pi*u*SC['ku']) * np.sin(np.pi*v*SC['kv']))
        s = np.where(s == 0, 1, s)
        alb = m*np.where(s > 0, 1.0, 0.2)[..., None]
    elif t == 'stripes':
        s = np.where(np.sign(np.sin(np.pi*u*SC['ku'])) == 0, 1,
                     np.sign(np.sin(np.pi*u*SC['ku'])))
        alb = m*np.where(s > 0, 1.0, 0.2)[..., None]
    elif t == 'marble':
        mm = 0.5 + 0.5*np.sin(2*np.pi*(SC['ku']*u + 0.35*np.sin(3*np.pi*SC['kv']*v)))
        alb = m*(0.55 + 0.45*mm)[..., None]
    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 = 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]
            alb = (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)
    if SC['legacy']:
        alb = srgb_encode(alb)
    return alb

def render_symbolic_tex(X, Y, SC=SC5, 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)
    # chart evaluated EVERYWHERE (totality); acos clamped total
    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
    albA = textureV(uA, vA, SC)
    s = np.sign(np.sin(np.pi*X)*np.sin(np.pi*Y)); s = np.where(s == 0, 1, s)
    base = np.where(s > 0, 0.85, 0.15)
    albF = np.stack([base*(1-SC['alpha']), base*(1-SC['alpha']), base], -1)
    if SC['legacy']: albF = srgb_encode(albF)
    albB = np.zeros(X.shape + (3,)) + np.array(SC['matB'])
    if SC['legacy']: albB = srgb_encode(albB)
    alb = albA*onA[..., None] + albB*onB[..., None] + albF*onF[..., None]
    amb = np.array(SC['ambient'])
    if SC['legacy']: amb = srgb_encode(amb)
    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 []))):
        t3 = t if C is not None else np.ones(X.shape)
        Px, Py, Pz = t3*X, t3*Y, e*(1 - t3)
        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']*amb
        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))      # Shadow Lemma
            wn = np.sqrt(a2)
            mu = (Nx*wx + Ny*wy + Nz*wz)/wn
            if SC['model'] == 'phong':
                wv = (wx*vx + wy*vy + wz*vz)/wn
                nu = np.clip((Nx*vx + Ny*vy + Nz*vz), -1, 1)
                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']
            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]
    uvm = np.stack([uA, vA], -1)
    return np.clip(np.nan_to_num(img), 0, 1), code, uvm, alb

X, Y = grid(256, 256)
t0 = time.perf_counter()
img_sym, code_sym, uvm_sym, alb_sym = render_symbolic_tex(X, Y)
t1 = time.perf_counter()
print(f"symbolic textured render: {t1 - t0:.4f} s")

same = code_num == code_sym
d = np.abs(img_num - img_sym)
onA_both = (code_num & 3 == 1) & same
duv = np.abs(uvm_num[onA_both] - uvm_sym[onA_both])
print(f"interior max |num - sym|  = {d[same].max():.3e}   over {same.sum()} px")
print(f"flips: {(~same).sum()} px")
print(f"max |du,dv| on sphere A   = {duv.max():.3e}   (the chart agrees exactly)")
dalb = np.abs(alb_num[onA_both] - alb_sym[onA_both])
print(f"max |d_tau| on sphere A   = {dalb.max():.3e}   (bitmap included)")

fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(srgb_encode(img_num)); axs[0].set_title("Numeric (sRGB)")
axs[1].imshow(srgb_encode(img_sym)); axs[1].set_title("Symbolic (sRGB)")
ov = srgb_encode(img_sym).copy(); ov[~same] = [1, 1, 1]
axs[2].imshow(ov); axs[2].set_title("flips in white")
for ax in axs: ax.axis('off')
plt.show()

print(f"{'texture':>14} {'numeric loop 128':>18} {'symbolic vec 256':>20}")
for tex in ('solid', 'checker', 'stripes', 'marble', 'bitmap'):
    for smp in (['nearest', 'bilinear'] if tex == 'bitmap' else ['bilinear']):
        SCt = {**SC5, 'tex': tex, 'bmp_sampling': smp}
        t0 = time.perf_counter(); render_numeric_tex(128, 128, SCt);  t1 = time.perf_counter()
        render_symbolic_tex(X, Y, SCt);                               t2 = time.perf_counter()
        tag = tex + ('/' + smp if tex == 'bitmap' else '')
        print(f"{tag:>14} {t1-t0:>17.3f}s {t2-t1:>19.3f}s")

print("""
Reading the table:
 * Closed-form textures cost a few transcendentals per pixel; the bitmap costs
   indexing (nearest) or four texel reads and three lerps (bilinear).
 * The symbolic form pays the chart and tau for EVERY pixel (totality), so its
   texture cost is constant across the frame; the numeric loop pays only on hits.
 * Bilinear vs nearest is not a quality knob here -- it is a choice between two
   EXACT functions: piecewise-constant vs piecewise-polynomial (Cell 6).
""")

cc = sp.symbols('c', positive=True)
E_expr = sp.Piecewise((12.92*cc, cc <= sp.Rational(7803, 10**6)),
                      (1.055*cc**sp.Rational(5, 12) - 0.055, True))
a = sp.Rational(1, 2)
diff = sp.N(E_expr.subs(cc, 2*a) - 2*E_expr.subs(cc, a))
print(f"E(0.5 + 0.5) - [E(0.5) + E(0.5)] = {diff}")
print("nonzero => lighting in gamma space is WRONG, by a computable amount.")
print("Correct order: shade in linear, encode at display. (Theorem, not taste.)")

# --- measure the legacy error over the whole frame -----------------------------
img_correct, _, _, _ = render_symbolic_tex(X, Y, {**SC5, 'legacy': False})
img_correct = srgb_encode(img_correct)
img_legacy, _, _, _  = render_symbolic_tex(X, Y, {**SC5, 'legacy': True})
img_naive            = np.clip(img_correct / 1.0, 0, 1)   # linear raw (no encode)
img_naive            = np.clip(np.nan_to_num(
    render_symbolic_tex(X, Y, {**SC5, 'legacy': False})[0]), 0, 1)

surf = (code_sym & 3) > 0
err = np.abs(img_legacy - img_correct)
print(f"\nlegacy vs correct: mean |delta| on surfaces = {err[surf].mean():.4f}, "
      f"max = {err[surf].max():.4f}")
lin = np.nan_to_num(render_symbolic_tex(X, Y, SC5)[0])
mid = surf & (lin.mean(axis=-1) > 0.1) & (lin.mean(axis=-1) < 0.6)
print(f"midtone mean |delta| = {err[mid].mean():.4f}   "
      f"(the washout lives in the midtones)")

fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(img_correct); axs[0].set_title("correct: shade linear, encode last")
axs[1].imshow(img_legacy);  axs[1].set_title("legacy: lighting in gamma space")
axs[2].imshow(np.clip(err*4, 0, 1)); axs[2].set_title("|delta| x4: the measured error")
for ax in axs: ax.axis('off')
plt.show()

# --- the sRGB piecewise itself, with its kink -----------------------------------
cg = np.linspace(0, 1, 500)
plt.figure(figsize=(6, 3.4))
plt.plot(cg, srgb_encode(cg), color='#1f3a5f', lw=2, label='sRGB encode E(c)')
plt.plot(cg, cg, color='grey', ls=':', label='identity (Chapters 1-4 display)')
plt.axvline(0.0031308, color='#b8860b', lw=1)
plt.text(0.006, 0.09, 'kink: 12.92c | 1.055c^(1/2.4)-0.055', fontsize=9)
plt.title("the display encode is a Piecewise"); plt.legend(); plt.show()

# The bitmap SAMPLES an ideal ring function at texel centers. Compare its two
# exact reconstructions against the ideal, then sweep frequency to find aliasing.
uu, vv = np.meshgrid(np.linspace(0, 1, 512, endpoint=False),
                     np.linspace(0, 1, 512, endpoint=False))
ideal = np.where((np.floor(np.hypot((uu*NB - 8) % NB - 7.5 + 0*uu,
                                    (vv*NB - 8) % NB - 7.5)) % 2) == 0,
                 0.85, 0.15)
# (ideal is defined texel-wise: the function the bitmap sampled at centers)
gi, gj = np.meshgrid(np.arange(NB), np.arange(NB))
ideal_centers = np.where((np.floor(np.hypot(gi - 7.5, gj - 7.5)) % 2) == 0, 0.85, 0.15)

def recon(uu, vv, mode):
    if mode == 'nearest':
        return ideal_centers[(np.floor(vv*NB)).astype(int) % NB,
                             (np.floor(uu*NB)).astype(int) % NB]
    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, gy - y0
    px = lambda a, b: ideal_centers[b % NB, a % NB]
    return (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)

fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(ideal_centers, cmap='gray', interpolation='nearest')
axs[0].set_title("the bitmap (what was sampled)")
axs[1].imshow(recon(uu, vv, 'nearest'), cmap='gray', interpolation='nearest')
axs[1].set_title("nearest: piecewise-constant function")
axs[2].imshow(recon(uu, vv, 'bilinear'), cmap='gray', interpolation='nearest')
axs[2].set_title("bilinear: piecewise-polynomial function")
for ax in axs: ax.axis('off')
plt.show()

# --- Nyquist, measured: checker frequency vs reconstruction error ---------------
freqs = np.arange(1, 17)
errs_near, errs_bil = [], []
for k in freqs:
    f = np.sign(np.sin(np.pi*k*uu)*np.sin(np.pi*k*vv)) > 0
    # sample f into a k-scaled bitmap grid of NB texels, reconstruct, compare
    fb = (np.sign(np.sin(np.pi*k*(np.arange(NB)+0.5)/NB)
                  [None, :]*np.ones((NB, 1))) > 0)
    fb2 = fb * (np.sign(np.sin(np.pi*k*(np.arange(NB)+0.5)/NB))[:, None] > 0)
    rn = fb2[(np.floor(vv*NB)).astype(int), (np.floor(uu*NB)).astype(int)]
    errs_near.append(np.abs(rn - f).mean())
    # bilinear on the boolean pattern as floats
    fbf = fb2.astype(float)
    gx, gy = uu*NB - 0.5, vv*NB - 0.5
    x0, y0 = np.floor(gx).astype(int), np.floor(gy).astype(int)
    fx2, fy2 = gx - x0, gy - y0
    px = lambda a, b: fbf[b % NB, a % NB]
    rb = (px(x0,y0)*(1-fx2)*(1-fy2) + px(x0+1,y0)*fx2*(1-fy2)
          + px(x0,y0+1)*(1-fx2)*fy2 + px(x0+1,y0+1)*fx2*fy2)
    errs_bil.append(np.abs(np.clip(rb, 0, 1) - f).mean())

plt.figure(figsize=(6.5, 3.6))
plt.plot(freqs, errs_near, 'o-', color='#7a1f1f', label='nearest')
plt.plot(freqs, errs_bil, 's-', color='#1f3a5f', label='bilinear')
plt.axvline(NB/2, color='#b8860b', lw=1.5)
plt.text(NB/2 + 0.2, 0.05, 'NB/2 = 8: the bitmap Nyquist', fontsize=9)
plt.xlabel('checker frequency k'); plt.ylabel('mean reconstruction error')
plt.title('sampling a symbolic function has a Nyquist limit'); plt.legend(); plt.show()
print("Past k = NB/2 the error stops growing gracefully and starts oscillating:")
print("aliasing is a theorem about frequency content, symbolic or not.")

# (a) THE SEAM TEAR: use a chart with the WRONG wrap (u in [-0.5, 0.5) via bare
#     arctan instead of atan2) and the texture discontinuity becomes visible as
#     a one-pixel tear along y = 0, x < 0 -- exactly the slit of section 3.
u_bad = np.arctan(np.where(np.abs(X) > 1e-12, Y/np.maximum(np.abs(X), 1e-30), 0))
u_bad = (u_bad/np.pi) % 1.0          # wrong chart: loses quadrant information
alb_bad = textureV(u_bad, uvm_sym[..., 1], {**SC5, 'tex': 'checker'})
alb_ok  = textureV(uvm_sym[..., 0], uvm_sym[..., 1], {**SC5, 'tex': 'checker'})
onA = (code_sym & 3) == 1
tear = onA & (np.abs(alb_bad - alb_ok).sum(axis=-1) > 0.5)
print(f"wrong-chart tear pixels on sphere A: {tear.sum()} "
      f"(they lie on the seam slit y = 0, x < 0)")
fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(np.where(onA[..., None], alb_bad, 0.0))
axs[0].set_title("wrong chart (bare arctan): the tear")
axs[1].imshow(np.where(onA[..., None], alb_ok, 0.0))
axs[1].set_title("correct chart (atan2)")

# (b) THE POLE PINCH, MEASURED: |du| between horizontally adjacent pixels.
du = np.abs(np.diff(uvm_sym[..., 0], axis=1))
du = np.minimum(du, 1 - du)          # circular distance across the wrap
cx, cy = du.shape[0]//2, du.shape[1]//2
ring_near = du[cx-8:cx+8, cy-8:cy+8].mean()
ring_far  = du[cx-8:cx+8, du.shape[1]//4*3-8:du.shape[1]//4*3+8].mean()
print(f"\nmean |du| per pixel step near the pole (center): {ring_near:.5f}")
print(f"mean |du| per pixel step mid-sphere:               {ring_far:.5f}")
print(f"magnification ratio ~ {ring_near/max(ring_far,1e-30):.1f}x -- "
      f"the pinch, measured")
axs[2].imshow(onA & (du > np.percentile(du[onA], 99)), cmap='autumn')
axs[2].set_title("top 1% |du|: the pole pinch, located")
for ax in axs: ax.axis('off')
plt.show()

# (c) RAW HEAVISIDE + TEXTURE: the census, again
img_raw, _, _, _ = render_symbolic_tex(X, Y, SC5, regularized=False)
print(f"\nraw assembly NaN pixels: {np.isnan(img_raw).any(axis=-1).sum()}")
print("tau itself is total (atan2, acos-clamped, floors): the NaNs still come")
print("from the chamber sqrts. Regularize; the texture rides along unharmed.")

# (a) WHICH TEXTURE BOUNDARIES ARE ALGEBRAIC?
# Checker edges: sin(pi*ku*u) = 0  <=>  u = m/ku  <=>  atan2(y, x) = const
#   => RAYS from the window center (algebraic!).
#               sin(pi*kv*v) = 0  <=>  v = m/kv  <=>  r^2 = const
#   => CIRCLES (algebraic!).
# So for THIS chart + closed-form tau, Chapter 3's CAD machinery still applies.
print("checker boundaries = rays u=const and circles v=const: semi-algebraic.")
print("The transcendental chart washed out -- atan2 cancelled, Nz was r^2.")
print("But for a BITMAP tau, boundaries are images of arbitrary curves through")
print("the chart: no Sturm, no CAD. Exact integration over textured pixels")
print("(Chapter 6) splits along curves it can only sample. Honest frontier.")
r2s = sp.symbols('r2', positive=True)
Nz_r = sp.simplify(Nz.xreplace({x**2: r2s - y**2}))
print("\nNz(r^2) =", Nz_r)

# (b) BUMP MAPPING IS THE NEXT COMPOSITION: N' = normalize(N + beta * grad tau).
#     It perturbs the normal by the TEXTURE'S DERIVATIVE -- which symbolic
#     rendering computes exactly (the image was already differentiable in Ch.1).
su, sv = sp.symbols('u v', real=True)
tau_b = sp.sin(sp.pi*8*su)*sp.sin(sp.pi*4*sv)      # checker, smoothed
grad = (sp.diff(tau_b, su), sp.diff(tau_b, sv))
print(f"\ntau(u,v) = {tau_b}")
print(f"dtau/du = {grad[0]}   ({sp.count_ops(grad[0])} ops)")
print(f"dtau/dv = {grad[1]}   ({sp.count_ops(grad[1])} ops)")
print("N' = (N + beta*(tau_u * Bu + tau_v * Bv)) / |...| : pure composition,")
print("the derivative exact. Chapter 9's headline, priced in advance.")