#!/usr/bin/env python3 # ═════════════════════════════════════════════════════════════════════════════ # PLAYER COMPUTER — Mammoth Thaw (22/32) # by Gene Kogan · 2026 · https://genekogan.com/player_computer/mammoth_thaw # # A mammoth thaws out of the permafrost, seen only in X-ray plates. # # This single file IS the piece: it draws every frame, synthesizes every sound, # and muxes them into the final video with ffmpeg. No other project files are # needed. You (or your agent) are invited to make a VARIATION of it: # # Generate a variation of this music video using only code. # Start from https://genekogan.com/player_computer/code/mammoth_thaw.py.txt # # The original render (for reference, yours should differ): # video: https://genekogan.com/player_computer/media/mammoth_thaw.mp4 # cover: https://genekogan.com/player_computer/media/mammoth_thaw.jpg # # Requirements: python3, numpy, pillow, and ffmpeg on PATH. # pip install numpy "pillow<13" # Speech/vocals (in pieces that have them) use the macOS `say` command; on # other platforms swap in espeak-ng / any TTS at the say_wav()/speak() calls, # or mute those lines. git provenance stamps degrade gracefully outside a repo. # Run: python3 mammoth_thaw.py (writes frames/, audio/, and the final mp4 # next to the script; writes ~8 GB of frames, # takes 5-15 min on a modern machine) # ═════════════════════════════════════════════════════════════════════════════ """ player_computer_final — "MAMMOTH THAW" (final delivery cut) Nordic folk fiddle (hardanger drones, open fifths) that rots into SLUDGE DOOM METAL. 62 bpm, D. 16 bars. The transition IS the arc. A woolly mammoth comes out of permafrost block 7. You never see the mammoth. You only ever see the PLATES — a radiology lab keeps shooting X-rays of the block as it thaws, and the film is those plates, one after another, on a light box. Ice is nearly transparent; ivory is off the scale; the skeleton resolves out of the milk. Somewhere in the sludge a soft-tissue shadow between the ribs starts to have an opinion. Then the heart restarts on the plate, and the animal walks off the right edge of the film and out of the beam. VISUAL SUBSTRATE (new to this repo): a real radiographic density accumulator. Every bone is a 3D primitive (tapered capsule / ellipsoid / fbm slab); for each pixel we integrate the analytic PATH LENGTH of a parallel ray through it times the material's density (ivory 1.9, bone 1.0, cartilage .34, muscle .11, ice .075, air 0). That density field goes through Beer-Lambert exposure, focal-spot unsharpness, quantum mottle that rises with attenuation, exposure blooming, detector fixed-pattern noise, a collimator, and a light box. Bone is white because bone stops photons, not because it was painted white. Composition: engine : audio-first x shot-parallel x world-camera content: audio-groove (hardanger fiddle -> downtuned sludge) x radiograph x tts-voices (Daniel, radiologist dictation) x effects-post FINAL CUT (player_computer_final). The film is unchanged; the delivery is: * 1920x1080 native. RS = H/720 = 1.5. This substrate is a physical simulation, so the rule is that WORLD scale and LENGTHS scale and COUNTS do not: the camera's pixels-per-metre goes up 1.5x (so the ray-sum is genuinely integrated at 1080p rather than blown up), and with it the focal-spot sigma, the Compton scatter radius, the exposure bloom, the detector fixed-pattern noise scale, its dead columns, the dust and the hair on the light-box glass, the plate jitter and the exposure scan bar. Quantum mottle and film grain are drawn at 1280x720 and NEAREST-blown-up so a speck of noise is the same size on screen as it always was. The burned-in chrome is authored in 1280x720 units and passes through a ScaledDraw proxy. * No renderer-debug overlays. Everything burnt into this picture is a radiograph's own annotation — plate number, kVp/mAs, grid and SID, the specimen block, the lead L/R marker, the elapsed clock, ice depth and core temperature, the OD step wedge, the calipers, the ECG and the radiologist's dictation. That is the fiction, and all of it stays. * Title moment: the "MAMMOTH · THAW" card keeps its series sub-line and gains PLAYER COMPUTER under it, in the plate's own annotation ink. Run from repo root: python3 renders/player_computer_final/mammoth_thaw/render.py --sheet python3 renders/player_computer_final/mammoth_thaw/render.py --jobs 3 """ import argparse, datetime, hashlib, math, os, subprocess, wave from pathlib import Path import numpy as np from PIL import Image, ImageDraw, ImageFont, ImageFilter NAME = "mammoth_thaw" TITLE = "MAMMOTH THAW" SUBT = "PLAYER COMPUTER" SETDIR = "player_computer_final" SETNUM = "03" AW, AH = 1280, 720 # authoring frame — the chrome lives here W, H, FPS = 1920, 1080, 24 # ── delivery scale ────────────────────────────────────────────────────────── # Named RS to keep `S` free. Lengths scale; counts do not. RS = H / AH def PX(v): return int(round(v*RS)) def PF(v): return v*RS def _scale_xy(v, s): if isinstance(v, (list, tuple)): return [_scale_xy(u, s) for u in v] return v*s class ScaledDraw: """ImageDraw proxy that multiplies geometry by RS at rasterisation time. The chrome layer is authored in 1280x720 annotation coordinates; this is what puts it on a 1920x1080 sheet without hand-scaling every literal. Fonts are scaled inside font(), so `textlength` comes back in REAL pixels — use TL() when a measured width has to go back into authoring space. """ __slots__ = ("_d", "_s") _GEOM = frozenset(("line", "rectangle", "rounded_rectangle", "ellipse", "polygon", "arc", "chord", "pieslice", "point", "text")) def __init__(self, d, s): self._d, self._s = d, s def __getattr__(self, name): f = getattr(self._d, name) if name not in self._GEOM: return f s = self._s def wrapped(xy, *a, **kw): w = kw.get("width") if w is not None: kw["width"] = max(2, int(round(w*s))) return f(_scale_xy(xy, s), *a, **kw) return wrapped def mkdraw(im): d = ImageDraw.Draw(im) return d if RS == 1.0 else ScaledDraw(d, RS) def TL(d, txt, f): """A measured text width, in AUTHORING units.""" return d.textlength(txt, font=f)/RS BPM = 62.0 BEAT = 60.0 / BPM BAR = 4 * BEAT SR = 44100 OUT = Path(__file__).resolve().parent FRAMES = OUT / "frames"; FRAMES.mkdir(exist_ok=True) AUD = OUT / "audio"; AUD.mkdir(exist_ok=True) ROOT = Path(__file__).resolve().parent # standalone: was repo root (used for git provenance) FONTS = ROOT / "fonts" SECTIONS = [ ("plate", 0, 3), # solo fiddle, ice, nothing but a ghost ("thaw", 3, 6), # drones + drips, skeleton resolving ("fracture", 6, 8), # first downtuned chord; ice cracks ("sludge", 8, 11), # full doom ("restart", 11, 13), # the heart ("walk", 13, 16), # it leaves ] N_BARS = SECTIONS[-1][2] DUR = N_BARS * BAR + 2.4 N_FRAMES = int(DUR * FPS) MUSIC_DESC = f"nordic hardanger fiddle -> sludge doom, {BPM:.0f}bpm, D, {N_BARS} bars" ENGINE_DESC = "radiographic density accumulator (ray-sum) / light box / plate stack" # ════════════════════════════════════════════════════════════════════════════ # AUDIO PRIMITIVES # ════════════════════════════════════════════════════════════════════════════ def mtof(m): return 440.0 * 2.0 ** ((m - 69) / 12.0) _PC = {"C":0,"C#":1,"Db":1,"D":2,"D#":3,"Eb":3,"E":4,"F":5,"F#":6,"Gb":6, "G":7,"G#":8,"Ab":8,"A":9,"A#":10,"Bb":10,"B":11} def nf(name): i = 2 if (len(name) > 2 and name[1] in "#b") else 1 return mtof(12 * (int(name[i:]) + 1) + _PC[name[:i]]) def adsr(n, a, d, s, r): e = np.zeros(n) ai, di, ri = max(1, int(a*SR)), max(1, int(d*SR)), max(1, int(r*SR)) ai = min(ai, n); e[:ai] = np.linspace(0, 1, ai) if ai < n: dd = min(di, n - ai) e[ai:ai+dd] = np.linspace(1, s, dd); e[ai+dd:] = s if ri < n: e[-ri:] *= np.linspace(1, 0, ri) return e def bandshape(x, lo=0.0, hi=0.0, order=4): """Exact FFT band shaping. Every noise source goes through this so nothing is a raw full-band blast (AESTHETIC 13a).""" n = len(x) if n < 8: return x X = np.fft.rfft(x); fq = np.maximum(np.fft.rfftfreq(n, 1/SR), 1e-6) g = np.ones_like(fq) if lo: g *= 1.0/np.sqrt(1.0 + (lo/fq)**order) if hi: g *= 1.0/np.sqrt(1.0 + (fq/hi)**order) return np.fft.irfft(X*g, n) def notch(x, f0, q=6.0, depth=0.85): n = len(x) if n < 8: return x X = np.fft.rfft(x); fq = np.maximum(np.fft.rfftfreq(n, 1/SR), 1e-6) g = 1.0 - depth*np.exp(-((np.log(fq/f0))**2)*(q*q)) return np.fft.irfft(X*g, n) def peaks(x, fs, gains, q=9.0): """Body resonances — a violin corpus, or a guitar cabinet.""" n = len(x) if n < 8: return x X = np.fft.rfft(x); fq = np.maximum(np.fft.rfftfreq(n, 1/SR), 1e-6) g = np.ones_like(fq) for f0, gg in zip(fs, gains): g += gg*np.exp(-((np.log(fq/f0))**2)*(q*q)) return np.fft.irfft(X*g, n) def comb_fb(x, freq, g=0.86, pad=0.6): """Feedback comb (sympathetic string) done in the frequency domain, because a per-sample IIR over 60 s of audio in Python is not a thing we do. y[n] = x[n] + g*y[n-L] -> Y = X / (1 - g*z^-L).""" L = max(2, int(round(SR/freq))) n = len(x); N = n + int(pad*SR) N = 1 << int(np.ceil(np.log2(N))) X = np.fft.rfft(x, N) k = np.arange(len(X)) Hd = 1.0 - g*np.exp(-2j*np.pi*k*L/N) return np.fft.irfft(X/Hd, N)[:n] # ---- the fiddle ----------------------------------------------------------- _VBODY = (280.0, 460.0, 730.0, 1180.0, 2600.0) _VGAIN = (1.5, 0.9, 1.3, 0.7, 0.5) def bow(freq, dur, g_noise=0.10, nh=34, a=.06, d=.25, s=.82, r=.14, scoop=0.02, grit=0.0, seed=0, tilt=5200.0): """Hardanger-ish bowed string. Additive odd+even stack with a 1/k^1.15 tilt (bowed strings are brighter than a plucked 1/k), a pitch scoop into the note from bow pressure, bow-scrape noise gated by the attack, and the violin corpus formants. Deliberately vibrato-free — Nordic fiddling doesn't wobble, it leans.""" n = int(dur*SR) if n <= 0: return np.zeros(0) t = np.arange(n)/SR f = freq*(1.0 - scoop*np.exp(-t*26.0)) ph = 2*np.pi*np.cumsum(f)/SR rng = np.random.RandomState(seed) out = np.zeros(n) for k in range(1, nh+1): fk = freq*k if fk > SR*0.45: break gk = (1.0/k**1.15) / np.sqrt(1.0 + (fk/tilt)**3) out += gk*np.sin(ph*k + rng.uniform(0, 2*np.pi)) if grit: # bow crunch: slight AM at ~48 Hz out *= 1.0 + grit*np.sin(2*np.pi*47.0*t + rng.uniform(0, 6)) env = adsr(n, a, d, s, r) scrape = bandshape(rng.randn(n), lo=1600, hi=7000)*np.exp(-t*9.0)*g_noise y = peaks(out*env + scrape*env, _VBODY, _VGAIN, q=7.0) return y/(np.max(np.abs(y))+1e-9)*0.9 def fiddle_drone(freq, dur, sway=0.55, seed=0): """An open string held under the melody, breathing with the bow arm.""" x = bow(freq, dur, a=.5, d=.6, s=.9, r=.6, g_noise=.05, seed=seed) t = np.arange(len(x))/SR return x*(0.72 + 0.28*np.sin(2*np.pi*(sway/BEAT/4)*t + seed)) # ---- the sludge ----------------------------------------------------------- def guitar(root, dur, chug=False, gain=9.0, fifth=True, octv=True, sub=True, seed=0, squeal=0.0): """Downtuned sludge guitar. Detuned saw stack (root + fifth + octave + sub), pre-EQ, asymmetric hard clip, 4x12 cabinet, second-stage soft clip. `chug` = palm mute: fast decay + a darker cab.""" n = int(dur*SR) if n <= 0: return np.zeros(0) t = np.arange(n)/SR rng = np.random.RandomState(seed) voices = [(1.0, 1.0), (1.0, 1.0)] if fifth: voices.append((1.4983, 0.85)) if octv: voices.append((2.0, 0.55)) if sub: voices.append((0.5, 0.75)) x = np.zeros(n) for vi, (mul, amp) in enumerate(voices): det = 1.0 + (vi - len(voices)/2)*0.0016 f0 = root*mul*det ph = 2*np.pi*f0*t + rng.uniform(0, 6) for k in range(1, 15): if f0*k > SR*0.45: break x += amp*np.sin(ph*k)/k x = bandshape(x, lo=52, hi=3000) x = np.tanh(x*gain + 0.22*np.tanh(x*gain*2.2)) # asymmetric grind x = notch(x, 380.0, q=3.0, depth=0.55) # scooped mids x = peaks(x, (105.0, 720.0, 2050.0), (1.1, 0.5, 0.35), q=4.0) x = bandshape(x, lo=68, hi=(2100 if chug else 3600)) # 4x12 cabinet x = np.tanh(x*1.9) if squeal: # feedback harmonic sq = np.sin(2*np.pi*root*6.02*t)*np.clip((t/dur)**2.2, 0, 1) x = x*(1-0.25*squeal) + sq*squeal*0.5 if chug: env = adsr(n, .004, .10, .16, .06) else: env = adsr(n, .010, .45, .70, .35) y = x*env return y/(np.max(np.abs(y))+1e-9)*0.95 def doom_kick(dur=.62, seed=1): n = int(dur*SR); t = np.arange(n)/SR f = 42 + 130*np.exp(-t*26) body = np.sin(2*np.pi*np.cumsum(f)/SR)*np.exp(-t*5.2) beat_ = np.sin(2*np.pi*np.cumsum(f*0.5)/SR)*np.exp(-t*3.0)*0.5 ck = bandshape(np.random.RandomState(seed).randn(n), lo=1200, hi=5000) ck *= np.exp(-t*180)*0.55 return np.tanh((body + beat_ + ck)*1.7)*0.98 def doom_snare(dur=.70, seed=2, crack=1.0): n = int(dur*SR); t = np.arange(n)/SR rng = np.random.RandomState(seed) nz = bandshape(rng.randn(n), lo=210, hi=7200) body = (np.sin(2*np.pi*168*t) + .7*np.sin(2*np.pi*242*t) + .4*np.sin(2*np.pi*331*t)) y = nz*np.exp(-t*7.5)*0.9*crack + body*np.exp(-t*13)*0.65 return np.tanh(y*1.5)*0.95 def doom_tom(f0=96, dur=.75, seed=3): n = int(dur*SR); t = np.arange(n)/SR f = f0*(1 + 0.7*np.exp(-t*14)) body = np.sin(2*np.pi*np.cumsum(f)/SR)*np.exp(-t*4.2) sk = bandshape(np.random.RandomState(seed).randn(n), lo=400, hi=2600) return np.tanh((body + sk*np.exp(-t*45)*.4)*1.4)*0.9 def ride(dur=1.1, seed=9, g=1.0): n = int(dur*SR); t = np.arange(n)/SR bell = sum(np.sin(2*np.pi*f*t) for f in (410, 663, 941, 1327, 1903)) nz = bandshape(np.random.RandomState(seed).randn(n), lo=3200, hi=9500) return (bell*np.exp(-t*9)*.09 + nz*np.exp(-t*4.5)*.20)*g def crash(dur=3.4, seed=13, g=1.0): n = int(dur*SR); t = np.arange(n)/SR nz = bandshape(np.random.RandomState(seed).randn(n), lo=900, hi=8800) return nz*(np.exp(-t*1.5) + .35*np.exp(-t*.36))*.55*g def gong(dur=5.0, seed=17): n = int(dur*SR); t = np.arange(n)/SR y = np.zeros(n) rng = np.random.RandomState(seed) for f in (54, 79, 118, 163, 227, 311, 409, 552): y += np.sin(2*np.pi*f*t*(1+0.0012*rng.randn()))/np.sqrt(f) nz = bandshape(rng.randn(n), lo=200, hi=4000)*np.exp(-t*2.2)*.25 return (y*np.exp(-t*0.7) + nz)*np.clip(t*60, 0, 1)*0.5 # ---- textures ------------------------------------------------------------- def stomp(seed=19, hard=1.0): """A boot on a floorboard. Nordic fiddling grooves on the foot, and this is also the doom kick's ancestor — the same thump, 3 bars early.""" n = int(.34*SR); t = np.arange(n)/SR rng = np.random.RandomState(seed) f = 58 + 90*np.exp(-t*40) body = np.sin(2*np.pi*np.cumsum(f)/SR)*np.exp(-t*13) board = bandshape(rng.randn(n), lo=180, hi=1400)*np.exp(-t*55)*0.45 plank = bandshape(rng.randn(n), lo=60, hi=260)*np.exp(-t*20)*0.35 return np.tanh((body + board + plank)*1.5)*0.8*hard def drip(seed=21, f0=1350): """A drop of meltwater falling into a tray. Pitch rises as the cavity shrinks — the classic drip formant sweep.""" n = int(.16*SR); t = np.arange(n)/SR f = f0*(1 + 1.5*(t/0.16)**1.6) y = np.sin(2*np.pi*np.cumsum(f)/SR)*np.exp(-t*22) tick = bandshape(np.random.RandomState(seed).randn(n), lo=2600, hi=9000) return (y + tick*np.exp(-t*260)*.5)*0.5 def ice_crack(seed=23, dur=.9): n = int(dur*SR); t = np.arange(n)/SR rng = np.random.RandomState(seed) nz = bandshape(rng.randn(n), lo=260, hi=3400)*np.exp(-t*18) ping = np.sin(2*np.pi*(720*np.exp(-t*3)+180)*t)*np.exp(-t*6)*.5 grind = bandshape(rng.randn(n), lo=90, hi=520)*np.exp(-t*4)*.55 return np.tanh((nz + ping + grind)*1.6)*0.7 def heart_beat(dur=1.5, f=44.0, seed=27, g=1.0): """lub-dub. Two thumps, the second softer and lower.""" n = int(dur*SR); y = np.zeros(n) for j, (off, amp, ff, dec) in enumerate(((0.0, 1.0, f, 9.0), (0.30, 0.62, f*0.82, 12.0))): m = int((dur-off)*SR); t = np.arange(m)/SR thump = np.sin(2*np.pi*np.cumsum(ff*(1+2.1*np.exp(-t*40)))/SR) thump *= np.exp(-t*dec) body = bandshape(np.random.RandomState(seed+j).randn(m), lo=40, hi=210)*np.exp(-t*17)*0.4 y[int(off*SR):int(off*SR)+m] += (thump + body)*amp return np.tanh(y*1.8)*0.9*g def wind(dur, seed=31): """Tundra: pitched, drifting, filtered noise — never a static hiss.""" n = int(dur*SR); rng = np.random.RandomState(seed) out = np.zeros(n); blk = 8192 for i in range(0, n, blk): m = min(blk, n-i); u = i/max(1, n) fc = 240 + 340*math.sin(u*7.1) + 300*math.sin(u*2.3) seg = rng.randn(m+512) out[i:i+m] = bandshape(seg, lo=max(60, fc*0.5), hi=fc*2.1)[:m] env = 0.5 + 0.5*np.sin(2*np.pi*np.arange(n)/SR*0.07) return out*env*0.18 def reverb(x, rt=2.4, mix=.34, seed=29, pre=0.02): n = int(rt*SR); t = np.arange(n)/SR rng = np.random.RandomState(seed) ir = rng.randn(n)*np.exp(-t*(5.0/rt)) ir[:int(pre*SR)] = 0 ir /= np.abs(ir).sum()/40.0 + 1e-9 L = 1 << int(np.ceil(np.log2(len(x)+n))) wet = np.fft.irfft(np.fft.rfft(x, L)*np.fft.rfft(ir, L))[:len(x)] wet /= np.max(np.abs(wet))+1e-9 return x*(1-mix) + wet*mix*(np.max(np.abs(x))+1e-9) def delay(x, time=.36, fb=.42, mix=.22, taps=8): d = int(time*SR); out = x.copy() for i in range(1, taps+1): g = mix*(fb**i); s = d*i if s >= len(x): break out[s:] += x[:len(x)-s]*g return out # ════════════════════════════════════════════════════════════════════════════ # VOICE — the radiologist dictating over the plates # ════════════════════════════════════════════════════════════════════════════ def _h(*p): return hashlib.md5("|".join(str(q) for q in p).encode()).hexdigest()[:16] def read_wav(p): with wave.open(str(p)) as w: ch = w.getnchannels() x = np.frombuffer(w.readframes(w.getnframes()), " macOS `say` (the canonical voices) -> espeak-ng / # espeak (Linux; language mapped from the say voice name, rate is wpm in both) # -> Windows SAPI (default voice, rate mapped from wpm) -> timed silence as # the last resort (duration from a chars/wpm heuristic, loud warning, never # cached so a later run with an engine present re-voices). A re-voiced film is # a different performance of the same score; that is by design. # Force a tier with POOP_TTS=say|espeak|sapi|none. def _tts_lang(voice): v = str(voice) if "Spanish" in v: return "es-mx" if "Mexico" in v else "es" if "Portuguese" in v: return "pt-br" if "Brazil" in v else "pt" if "English (UK)" in v: return "en-gb" return "en-us" def _tts_engine(): import shutil, platform want = os.environ.get("POOP_TTS", "").strip().lower() if want: return want if shutil.which("say"): return "say" if shutil.which("espeak-ng") or shutil.which("espeak"): return "espeak" if platform.system() == "Windows": return "sapi" return "none" def _tts_render(text, voice, rate, path): """Synthesize text -> mono 44.1k wav at `path` with the best available engine. Returns False if no engine (caller falls back to timed silence).""" import shutil, sys, base64 eng = _tts_engine() tmp = path.with_suffix(".tts.wav") try: if eng == "say": aiff = path.with_suffix(".aiff") subprocess.run(["say", "-v", voice, "-r", str(rate), "-o", str(aiff), text], check=True) subprocess.run(["ffmpeg", "-y", "-i", str(aiff), "-ar", str(SR), "-ac", "1", str(path)], check=True, capture_output=True) aiff.unlink(missing_ok=True) return True if eng == "espeak": exe = shutil.which("espeak-ng") or shutil.which("espeak") or "espeak-ng" subprocess.run([exe, "-v", _tts_lang(voice), "-s", str(int(rate)), "-w", str(tmp), str(text)], check=True) elif eng == "sapi": r = max(-10, min(10, round((int(rate) - 175) / 25))) esc = str(text).replace("'", "''") ps = ("Add-Type -AssemblyName System.Speech;" "$s=New-Object System.Speech.Synthesis.SpeechSynthesizer;" f"$s.Rate={r};$s.SetOutputToWaveFile('{tmp}');" f"$s.Speak('{esc}');$s.Dispose()") enc = base64.b64encode(ps.encode("utf-16-le")).decode() subprocess.run(["powershell", "-NoProfile", "-EncodedCommand", enc], check=True) else: return False subprocess.run(["ffmpeg", "-y", "-i", str(tmp), "-ar", str(SR), "-ac", "1", str(path)], check=True, capture_output=True) return True except Exception as e: print(f"[tts] {eng} failed ({e}) — falling back to timed silence", file=sys.stderr) return False finally: tmp.unlink(missing_ok=True) def _tts_silence(text, rate): import sys dur = max(0.6, len(str(text)) / (max(60, int(rate)) * 5.0 / 60.0)) print(f'[tts] no speech engine — timed silence ({dur:.2f}s): ' f'"{str(text)[:48]}"', file=sys.stderr) return np.zeros(int(dur * SR)) def say_wav(text, voice, rate, path): """text -> mono 44.1k voice wav (cached on disk; deterministic per engine).""" path = Path(path) if not path.exists(): if not _tts_render(text, voice, rate, path): return _tts_silence(text, rate) return read_wav(path) def dictate(text, voice="Daniel", rate=152, seed=0): """Dictaphone: band-limited, a hair of tape wobble, small room.""" x = say_wav(text, voice, rate, AUD/("say_"+_h(text, voice, rate)+".wav")) x = bandshape(x, lo=190, hi=4600) x = peaks(x, (1100.0,), (0.6,), q=3.0) x = np.tanh(x*1.5) x = reverb(x, rt=0.9, mix=.20, seed=800+seed) return x/(np.max(np.abs(x))+1e-9) # The dictation. (bar, text, on-screen annotation, gain) VO = [ (0.55, "Plate one. Permafrost block seven. Subject intact.", "PLATE 001 / BLOCK 7 / SUBJ INTACT", 1.00), (3.20, "Ivory density exceeds the scale.", "IVORY: DENSITY OFF SCALE", 1.00), (5.15, "Ice is going. Keep plating.", "ICE RETREAT 4 cm/hr — CONTINUE", 1.00), (7.10, "There is something in the thoracic cavity.", "THORACIC CAVITY: MASS?", 1.05), (10.30, "That is not artifact.", "NOT ARTIFACT", 1.10), (11.55, "That is contraction.", "CONTRACTION. 4 bpm.", 1.15), (13.15, "It is walking off the plate.", "SUBJECT AMBULATORY", 1.10), (14.90, "Keep the light box on.", "LIGHT BOX: ON", 1.05), ] # ════════════════════════════════════════════════════════════════════════════ # THE SONG # ════════════════════════════════════════════════════════════════════════════ class Song: def __init__(self, dur): self.n = int(dur*SR) self.tr = {} self.kick_t = [] def t(self, bar, step=0, swing=0.0): sw = swing*(BEAT/4) if (step % 2) else 0.0 return bar*BAR + step*(BEAT/4) + sw def put(self, track, sig, at, g=1.0, pan=0.0): b = self.tr.setdefault(track, np.zeros((self.n, 2))) i = int(at*SR); j = min(self.n, i+len(sig)) if i >= self.n or j <= i: return if i < 0: sig = sig[-i:]; i = 0; j = min(self.n, i+len(sig)) if j <= i: return th = (pan*.5+.5)*(np.pi/2) b[i:j] += np.stack([sig[:j-i]*np.cos(th), sig[:j-i]*np.sin(th)], 1)*g def bus(self, track, fn): if track in self.tr: b = self.tr[track] self.tr[track] = np.stack([fn(b[:, 0]), fn(b[:, 1])], 1) def sec_env(self, levels, glide=0.7): env = np.ones(self.n) for nm, b0, b1 in SECTIONS: i0, i1 = int(b0*BAR*SR), min(self.n, int(b1*BAR*SR)) if i1 > i0: env[i0:i1] = levels.get(nm, 1.0) env[int(SECTIONS[-1][2]*BAR*SR):] = levels.get(SECTIONS[-1][0], 1.0) k = max(1, int(glide*SR)) return np.convolve(env, np.ones(k)/k, "same") def mixdown(self, gains, pump_depth=.34, pump_rel=.30, levels=None, duck=None): mix = np.zeros((self.n, 2)) for k, b in self.tr.items(): if k == "vox": continue mix += b*gains.get(k, 1.0) if levels: mix *= self.sec_env(levels)[:, None] if self.kick_t: env = np.ones(self.n); rl = int(pump_rel*SR) shape = 1 - pump_depth*np.exp(-np.arange(rl)/(pump_rel*SR/4)) for at in self.kick_t: i = int(at*SR); j = min(self.n, i+rl) if i < self.n: env[i:j] = np.minimum(env[i:j], shape[:j-i]) env = np.convolve(env, np.ones(512)/512, "same") mix *= env[:, None] if duck is not None: mix *= (1.0 - 0.50*duck)[:, None] if "vox" in self.tr: mix += self.tr["vox"]*gains.get("vox", 1.0) a = math.exp(-2*math.pi*26.0/SR) # DC / rumble trim for c in range(2): col = mix[:, c]; lp = np.empty(self.n); z = 0.0 for i in range(0, self.n, 8192): blk = col[i:i+8192] for j in range(len(blk)): z = (1-a)*blk[j] + a*z; lp[i+j] = z mix[:, c] = col - lp mix = np.tanh(mix*1.22)/np.tanh(1.22) return mix/(np.max(np.abs(mix))+1e-9)*0.94 def write(self, path, mix): with wave.open(str(path), "w") as w: w.setnchannels(2); w.setsampwidth(2); w.setframerate(SR) w.writeframes((np.clip(mix, -1, 1)*32767).astype("= 2: s.put("fiddle", bow(f*2**(2/12.0), (BEAT/4)*0.55, g_noise=.06, a=.008, d=.05, s=.4, r=.04, seed=bar*29+j+700), s.t(bar, st, SWF)-(BEAT/4)*0.42, g=mg*acc*0.42, pan=-.10) s.put("fiddle", bow(f, ln, g_noise=.09+.10*heavy, grit=.22 if heavy else .0, a=.030 if st % 4 == 0 else .045, d=.18, s=.80, r=.10, seed=bar*29+j), s.t(bar, st, SWF), g=mg*acc, pan=-.18+.36*((j % 3)/2.0)) # the double-stop fifth below — the hardanger signature if st % 4 == 0: s.put("fiddle", bow(f/1.4983, ln*1.15, g_noise=.07, seed=bar*29+j+400), s.t(bar, st, SWF), g=mg*acc*0.58, pan=-.34) # -------- the boot --------------------------------------------------- if sec in ("plate", "thaw"): hardness = .55 if sec == "plate" else .78 for st in (0, 8): at = s.t(bar, st) s.put("drums", stomp(seed=bar*3+st, hard=1.0), at, g=hardness, pan=-.10) s.kick_t.append(at) if sec == "thaw": s.put("drums", stomp(seed=bar*3+90, hard=.55), s.t(bar, 6, SWF), g=.34, pan=-.16) if bar % 2 == 1: s.put("drums", stomp(seed=bar*3+70, hard=.45), s.t(bar, 14, SWF), g=.28, pan=.14) # -------- meltwater drips ------------------------------------------ if sec in ("plate", "thaw", "walk"): nd = {"plate": 2, "thaw": 5, "walk": 3}[sec] for q in range(nd): st = int(R.randint(0, 16)) s.put("fx", drip(seed=bar*13+q, f0=1000+700*R.rand()), s.t(bar, st) + R.rand()*0.05, g=.26, pan=-.6+1.2*R.rand()) # -------- the sludge ------------------------------------------------ if heavy or walk_: amt = {"fracture": .55, "sludge": 1.0, "restart": .85, "walk": .40}[sec] if sec == "fracture": # the guitar arrives as one enormous held chord per half-bar for st in (0, 8): s.put("gtr", guitar(TUNING*2**((0 if st == 0 else 3)/12.0), BEAT*2.1, gain=7.0, seed=bar*7+st, squeal=.25 if st else 0.0), s.t(bar, st), g=.46*amt, pan=-.1) elif sec in ("sludge", "restart"): # two guitars, detuned seven cents apart and panned hard, plus # an octave-down sub layer under the open chords — the weight # comes from the beating between the two, not from more gain. rf = (RIFF, RIFF2)[bar % 2] for j, (st, deg, ch, ln16) in enumerate(rf): if sec == "restart" and bar == 11 and st > 3: continue d = (BEAT/4)*ln16*(0.95 if ch else 1.25) f0 = TUNING*2**(deg/12.0) for side, (cents, pan, gn, dl) in enumerate(( (-0.0040, -.62, 9.6, 0.000), (+0.0040, .62, 8.9, 0.009))): s.put("gtr", guitar(f0*(1+cents), d, chug=ch, gain=gn, seed=bar*11+j+side*900, squeal=0.0 if ch else .16), s.t(bar, st)+dl, g=.44*amt, pan=pan) if not ch: # the sub under the chord s.put("gtr", guitar(f0*0.5, d*1.15, gain=6.0, fifth=False, octv=False, seed=bar*11+j+1700), s.t(bar, st), g=.26*amt, pan=0.0) else: # walk: it recedes for st in (0,): s.put("gtr", guitar(TUNING, BAR*0.8, gain=5.0, seed=bar*17, squeal=.12), s.t(bar, st), g=.30*amt, pan=0.0) # -------- the kit ---------------------------------------------------- if sec == "fracture": for st in (0, 8): at = s.t(bar, st) s.put("drums", doom_kick(), at, g=.62); s.kick_t.append(at) s.put("drums", doom_snare(crack=.8), s.t(bar, 8), g=.44, pan=.05) for st in range(0, 16, 4): s.put("drums", ride(seed=90+bar*4+st, g=.7), s.t(bar, st), g=.20, pan=.28) elif sec in ("sludge", "restart"): heavier = sec == "sludge" for st in (0, 6, 8, 14) if heavier else (0, 8): at = s.t(bar, st) s.put("drums", doom_kick(), at, g=.72 if st in (0, 8) else .5) if st in (0, 8): s.kick_t.append(at) for st in (4, 12): s.put("drums", doom_snare(crack=1.0), s.t(bar, st), g=.58 if heavier else .40, pan=.04) for st in range(0, 16, 2): s.put("drums", ride(seed=200+bar*8+st, g=.9 if st % 4 == 0 else .55), s.t(bar, st), g=.22 if heavier else .13, pan=.30) if bar % 4 == 3: # the tom fill into the next for q, st in enumerate((10, 12, 13, 14, 15)): s.put("drums", doom_tom(f0=150-q*17, seed=bar*9+q), s.t(bar, st), g=.46, pan=-.4+.2*q) elif walk_: for st in (0,): at = s.t(bar, st) s.put("drums", doom_kick(dur=.8), at, g=.36); s.kick_t.append(at) for st in range(0, 16, 4): s.put("drums", ride(seed=400+bar*4+st, g=.5), s.t(bar, st), g=.13, pan=.25) # ---- the transitions: a run into every door ----------------------------- # A fiddle pickup lifts into the fracture and into the sludge, and a tom # crescendo pushes the last bar of thaw over the edge. The old version cut # between sections; this one is pushed across them. for tgt, degs in ((6, (0, 3, 5, 7, 10, 12)), (8, (12, 10, 12, 15, 17, 19))): for q, dg2 in enumerate(degs): at = tgt*BAR - (len(degs)-q)*(BEAT/4) s.put("fiddle", bow(fid_hi*2**(dg2/12.0), (BEAT/4)*1.1, g_noise=.11, grit=.10, a=.012, d=.08, s=.6, r=.05, seed=5000+tgt*17+q), at, g=.24+0.05*q, pan=-.30+.10*q) for q in range(6): # tom crescendo into the fracture s.put("drums", doom_tom(f0=170-q*15, seed=6100+q), 6*BAR - (6-q)*(BEAT/2), g=.16+0.08*q, pan=-.42+.16*q) for q in range(4): # and a snare pickup into the doom s.put("drums", doom_snare(crack=.55+0.15*q, seed=6200+q), 8*BAR - (4-q)*(BEAT/4), g=.20+0.10*q, pan=.06) # ---- events ------------------------------------------------------------- s.put("fx", ice_crack(seed=23, dur=1.2), 3*BAR - .18, g=.55, pan=-.35) s.put("fx", ice_crack(seed=41, dur=1.6), 6*BAR - .30, g=.85, pan=.30) s.put("fx", ice_crack(seed=43, dur=1.1), 6*BAR + BEAT*2, g=.55, pan=-.2) s.put("fx", ice_crack(seed=47, dur=1.9), 8*BAR - .32, g=.95, pan=.10) s.put("drums", crash(g=1.0, seed=13), 6*BAR - .05, g=.55, pan=-.12) s.put("drums", crash(g=1.2, seed=19), 8*BAR - .05, g=.72, pan=.12) s.put("drums", crash(g=1.0, seed=21), 11*BAR - .05, g=.48, pan=-.10) s.put("drums", crash(g=1.4, seed=25), 12*BAR - .05, g=.62, pan=.08) s.put("fx", gong(seed=17), 13*BAR - .3, g=.55, pan=0.0) # the crushing chord that lands exactly with the heart steadying (bar 12) for side, (cents, pan) in enumerate(((-0.005, -.66), (0.005, .66))): s.put("gtr", guitar(TUNING*(1+cents), BAR*1.9, gain=9.8, seed=7000+side, squeal=.22), 12*BAR, g=.58, pan=pan) s.put("gtr", guitar(TUNING*0.5, BAR*2.1, gain=6.4, fifth=False, octv=False, seed=7002), 12*BAR, g=.34, pan=0.0) s.put("fx", gong(dur=6.0, seed=33), 15*BAR + BEAT*2, g=.42, pan=0.0) # ---- the heart: the whole point of the piece --------------------------- # First flutters at the end of sludge; steadies through restart; carries # the walk out under everything. hb = [(10.55, .30, 46.), (10.90, .38, 45.), (11.00, .55, 44.), (11.35, .62, 44.), (11.70, .70, 43.), (12.00, .85, 42.), (12.35, .90, 42.), (12.70, 1.0, 41.), (13.00, 1.0, 41.), (13.50, .95, 41.), (14.00, .90, 40.), (14.50, .85, 40.), (15.00, .80, 40.), (15.50, .70, 40.), (16.00, .62, 40.)] for j, (bar, g, f) in enumerate(hb): s.put("heart", heart_beat(f=f, seed=27+j), bar*BAR, g=.55*g, pan=0.0) # ---- dictation --------------------------------------------------------- duck = np.zeros(s.n) for j, (bar, txt, _ann, g) in enumerate(VO): x = dictate(txt, seed=j) at = bar*BAR s.put("vox", x, at, g=.62*g, pan=0.04) i0 = int(at*SR); i1 = min(s.n, i0+len(x)+int(.25*SR)) if i1 > i0: duck[i0:i1] = np.maximum(duck[i0:i1], 1.0) duck = np.convolve(duck, np.ones(int(.22*SR))/int(.22*SR), "same") s.bus("fiddle", lambda x: reverb(delay(x, BEAT*.75, .28, .16), rt=2.8, mix=.30, seed=901)) s.bus("drone", lambda x: reverb(x, rt=5.0, mix=.46, seed=903)) s.bus("gtr", lambda x: reverb(x, rt=1.9, mix=.20, seed=907)) s.bus("drums", lambda x: reverb(x, rt=2.6, mix=.22, seed=909)) s.bus("fx", lambda x: reverb(x, rt=3.4, mix=.42, seed=911)) s.bus("heart", lambda x: reverb(x, rt=1.4, mix=.16, seed=913)) s.bus("air", lambda x: reverb(x, rt=6.0, mix=.5, seed=915)) mix = s.mixdown(dict(fiddle=1.0, drone=1.0, gtr=1.0, drums=1.0, fx=1.0, heart=1.0, air=1.0, vox=1.0), pump_depth=.22, pump_rel=.30, duck=duck, levels=dict(plate=.52, thaw=.68, fracture=.86, sludge=1.0, restart=.92, walk=.66)) wav = AUD/"final.wav" s.write(wav, mix) return wav, mix def analyze(mix): x = mix.mean(1) hop = SR/FPS; win = int(hop*1.8) E = {k: np.zeros(N_FRAMES) for k in ("rms", "sub", "low", "mid", "high")} for f in range(N_FRAMES): i = int(f*hop); seg = x[i:i+win] if len(seg) < 16: continue E["rms"][f] = np.sqrt((seg**2).mean()) sp = np.abs(np.fft.rfft(seg*np.hanning(len(seg)))) fr = np.fft.rfftfreq(len(seg), 1/SR) E["sub"][f] = sp[fr < 70].sum() E["low"][f] = sp[(fr >= 70) & (fr < 220)].sum() E["mid"][f] = sp[(fr >= 220) & (fr < 2200)].sum() E["high"][f] = sp[fr >= 2200].sum() for k in E: p = np.percentile(E[k], 96)+1e-9 E[k] = np.clip(E[k]/p, 0, 1.3) lo = E["low"] flux = np.maximum(0, lo - np.concatenate([[0], lo[:-1]])) E["hit"] = np.clip(np.convolve(flux, [.25, .5, .25], "same") / (np.percentile(flux, 97)+1e-9), 0, 1) np.savez(AUD/"env.npz", **E) return E _ENV = {} def env(): if not _ENV: z = np.load(AUD/"env.npz") for k in z.files: _ENV[k] = z[k] return _ENV # ════════════════════════════════════════════════════════════════════════════ # NOISE # ════════════════════════════════════════════════════════════════════════════ def value_noise(h, w, scale, seed): rng = np.random.RandomState(seed) gh, gw = int(h/scale)+2, int(w/scale)+2 g = rng.rand(gh, gw) ys = np.linspace(0, gh-1-1e-3, h); xs = np.linspace(0, gw-1-1e-3, w) y0 = ys.astype(int); x0 = xs.astype(int) fy = (ys-y0)[:, None]; fx = (xs-x0)[None, :] sy = fy*fy*(3-2*fy); sx = fx*fx*(3-2*fx) g00 = g[np.ix_(y0, x0)]; g01 = g[np.ix_(y0, x0+1)] g10 = g[np.ix_(y0+1, x0)]; g11 = g[np.ix_(y0+1, x0+1)] return ((g00*(1-sx)+g01*sx)*(1-sy) + (g10*(1-sx)+g11*sx)*sy).astype(np.float32) def blur_f(a, sigma): """Separable gaussian on a float32 plane. PIL will not filter mode 'F', and the density field must stay float — quantising to uint8 before the focal-spot blur throws away exactly the low-density gradients (ice, fat, fur) that this whole substrate is about.""" if sigma <= 0.05: return a.astype(np.float32) r = max(1, int(round(sigma*2.6))) k = np.exp(-0.5*(np.arange(-r, r+1)/sigma)**2) k = (k/k.sum()).astype(np.float32) p = np.pad(a, ((0, 0), (r, r)), mode="edge") o = np.zeros_like(a, np.float32) for j, wgt in enumerate(k): o += wgt*p[:, j:j+a.shape[1]] p = np.pad(o, ((r, r), (0, 0)), mode="edge") o2 = np.zeros_like(a, np.float32) for j, wgt in enumerate(k): o2 += wgt*p[j:j+a.shape[0], :] return o2 def blur_big(a, sigma, ds=8): """Wide glow. Mean-pool, blur small, blow back up — a 15px sigma done honestly is 80 full-frame adds per frame and we have 1544 frames.""" h, w = a.shape hh, ww = h//ds, w//ds s = a[:hh*ds, :ww*ds].reshape(hh, ds, ww, ds).mean((1, 3)) s = blur_f(s.astype(np.float32), max(0.7, sigma/ds)) up = np.repeat(np.repeat(s, ds, 0), ds, 1) o = np.empty_like(a, np.float32) o[:hh*ds, :ww*ds] = up if hh*ds < h: o[hh*ds:, :ww*ds] = up[-1:, :] if ww*ds < w: o[:, ww*ds:] = o[:, ww*ds-1:ww*ds] return o def fbm(h, w, scale, seed, octaves=4): out = np.zeros((h, w), np.float32); amp = 1.0; nrm = 0.0 for o in range(octaves): out += amp*value_noise(h, w, max(2.0, scale/(2**o)), seed+o) nrm += amp; amp *= .5 return out/nrm # ════════════════════════════════════════════════════════════════════════════ # THE RADIOGRAPH — a density accumulator # # Rays travel along +z into the plate. For every primitive we add # rho * (analytic path length of the ray through it) # to a density buffer D(x,y). Nothing here is a painted shape: the ribcage is # bright where ribs overlap because two path lengths added, and the ice is # milky because a 0.075-density slab 1.4 m thick is genuinely a little opaque. # # Materials (roughly relative linear attenuation at diagnostic energies): # ════════════════════════════════════════════════════════════════════════════ RHO_IVORY = 1.95 RHO_BONE = 1.00 RHO_CART = 0.34 RHO_MUSC = 0.115 RHO_FAT = 0.075 RHO_ICE = 0.072 RHO_HIDE = 0.055 RHO_FUR = 0.018 class Plate: """Density accumulator with a world->pixel camera. cam = (world_cx, world_cy, pixels_per_world_unit). World y is UP.""" __slots__ = ("D", "w", "h", "cx", "cy", "s") def __init__(self, w, h, cam): self.w, self.h = w, h self.cx, self.cy, self.s = cam # The ONE place the delivery scale enters the physics: metres-to-pixels # goes up with RS, so the ray-sum is integrated at the delivered # resolution instead of being enlarged after the fact. self.s *= RS self.D = np.zeros((h, w), np.float32) def px(self, x, y): return ((x-self.cx)*self.s + self.w*0.5, (self.cy-y)*self.s + self.h*0.5) def _bbox(self, x0, x1, y0, y1, pad): ix0 = int(max(0, math.floor(min(x0, x1)-pad))) ix1 = int(min(self.w, math.ceil(max(x0, x1)+pad))) iy0 = int(max(0, math.floor(min(y0, y1)-pad))) iy1 = int(min(self.h, math.ceil(max(y0, y1)+pad))) if ix1 <= ix0 or iy1 <= iy0: return None return ix0, ix1, iy0, iy1 def ell(self, cx, cy, ax, ay, az, rho): """Path length through an axis-aligned ellipsoid: 2*az*sqrt(1-u^2-v^2).""" sx, sy = self.px(cx, cy) AX, AY = ax*self.s, ay*self.s bb = self._bbox(sx-AX, sx+AX, sy-AY, sy+AY, 1) if bb is None or AX < .3 or AY < .3: return ix0, ix1, iy0, iy1 = bb Y, X = np.mgrid[iy0:iy1, ix0:ix1] u = (X-sx)/AX; v = (Y-sy)/AY q = 1.0 - u*u - v*v np.maximum(q, 0.0, out=q) self.D[iy0:iy1, ix0:ix1] += (2.0*az*rho)*np.sqrt(q, dtype=np.float32) def cap(self, p0, p1, r0, r1, rho, flat=1.0): """Tapered capsule (swept sphere). `flat` squashes the z-extent, for blades like the scapula that are wide but thin.""" sx0, sy0 = self.px(*p0); sx1, sy1 = self.px(*p1) R0, R1 = r0*self.s, r1*self.s rm = max(R0, R1)+2 bb = self._bbox(min(sx0, sx1)-rm, max(sx0, sx1)+rm, min(sy0, sy1)-rm, max(sy0, sy1)+rm, 1) if bb is None or rm < 1.2: return ix0, ix1, iy0, iy1 = bb Y, X = np.mgrid[iy0:iy1, ix0:ix1] dx, dy = sx1-sx0, sy1-sy0 L2 = dx*dx + dy*dy + 1e-6 t = np.clip(((X-sx0)*dx + (Y-sy0)*dy)/L2, 0.0, 1.0) ex = X - (sx0 + t*dx); ey = Y - (sy0 + t*dy) R = R0 + (R1-R0)*t q = R*R - (ex*ex + ey*ey) np.maximum(q, 0.0, out=q) self.D[iy0:iy1, ix0:ix1] += (2.0*rho*flat/self.s)*np.sqrt(q, dtype=np.float32) def arc(self, pts, radii, rho, flat=1.0): for j in range(len(pts)-1): self.cap(pts[j], pts[j+1], radii[j], radii[j+1], rho, flat) def slab(self, x0, x1, y0, y1, thick, rho, seed, scale=0.9, soft=0.35, melt=0.0, octaves=4): """The permafrost: a block whose thickness is an fbm field, eaten away from the outside as `melt` goes 0->1.""" sx0, sy0 = self.px(x0, y1); sx1, sy1 = self.px(x1, y0) bb = self._bbox(sx0, sx1, sy0, sy1, 2) if bb is None: return ix0, ix1, iy0, iy1 = bb hh, ww = iy1-iy0, ix1-ix0 if hh < 2 or ww < 2: return n = fbm(hh, ww, max(3.0, scale*self.s), seed, octaves=octaves) Y, X = np.mgrid[iy0:iy1, ix0:ix1] u = np.clip((X-sx0)/max(1.0, sx1-sx0), 0, 1) v = np.clip((Y-sy0)/max(1.0, sy1-sy0), 0, 1) edge = (np.minimum(np.minimum(u, 1-u), np.minimum(v, 1-v))/soft) np.clip(edge, 0, 1, out=edge) edge = edge*edge*(3-2*edge) th = thick*(0.45 + 1.0*n)*edge # subtract MORE than the maximum thickness (thick*1.45), unevenly, # so the retreat is patchy but complete th -= melt*thick*(1.78 - 0.32*n) np.maximum(th, 0.0, out=th) self.D[iy0:iy1, ix0:ix1] += th.astype(np.float32)*rho def crack(self, pts, wid, rho, seed=0): """A fracture plane: AIR inside the ice, so NEGATIVE density.""" rng = np.random.RandomState(seed) for j in range(len(pts)-1): a, b = pts[j], pts[j+1] self.cap(a, b, wid*(0.6+0.8*rng.rand()), wid*(0.6+0.8*rng.rand()), -rho) def hd_curve(v, speed=1.0, contrast=1.58, k=0.445, dmax=1.30, fog=0.014): """The film's characteristic (H&D) curve: toe, straight portion, shoulder. Without this the tusks simply clip — ivory at 1.95 rho over a 33 cm chord is four times a rib, so a linear map either crushes the ribcage or fuses the tusks into one white sausage. A shoulder is how real film survives a subject with that much dynamic range, so we grow one.""" x = np.maximum(v, 0.0)*speed xg = x**contrast return fog + dmax*xg/(xg + k**contrast) def develop(D, kvp=1.0, mas=1.0, focal=1.15, seed=0, mottle=1.0, scatter=0.0): """D (density * path length) -> a plate image, bone bright. Beer-Lambert -> H&D curve -> focal-spot unsharpness -> scatter veil -> quantum mottle.""" h, w = D.shape if focal > 0.02: D = blur_f(D, focal) mu = 1.05/max(0.35, kvp) # higher kVp = more penetrating v = hd_curve(1.0 - np.exp(-mu*D), speed=mas) if scatter > 0.01: # Compton fog off a metre and a half of ice: a low-frequency veil that # eats contrast. It is why the early plates are useless, and its # retreat is the film's other clock. v = (v*(1.0 - 0.66*scatter) + blur_big(v, PF(30.0))*1.05*scatter + 0.055*scatter) if mottle: rng = np.random.RandomState(seed) sig = (0.008 + 0.040*np.clip(v, 0, 1.4))*mottle if RS == 1.0: n = rng.normal(0, 1, (h, w)).astype(np.float32) else: # mottle is a grain size, so it is drawn at 720p and blown up gn = rng.normal(0, 1, (AH, AW))*32.0 + 128.0 gi = Image.fromarray(np.clip(gn, 0, 255).astype(np.uint8)) n = (np.asarray(gi.resize((w, h), Image.NEAREST), np.float32) - 128.0)/32.0 v = v + n*sig return np.clip(v, 0.0, 1.35) # --- the light-box grade ---------------------------------------------------- def _ramp(stops, n=256): stops = np.array(stops, np.float32) xs = np.linspace(0, 1, len(stops)); g = np.linspace(0, 1, n) return np.stack([np.interp(g, xs, stops[:, c]) for c in range(3)], 1) PLATE_LUT = _ramp([(3, 7, 15), (10, 20, 36), (26, 46, 70), (66, 92, 118), (136, 158, 178), (206, 220, 232), (243, 249, 255), (255, 255, 255)]) def colorize(v): i = np.clip(v*(len(PLATE_LUT)-1)/1.0, 0, len(PLATE_LUT)-1).astype(np.int32) return PLATE_LUT[i] _FIX = {} def fixed_pattern(h, w): """Detector fixed-pattern gain + a few dead columns. Same every frame — that is exactly what makes it read as one physical detector.""" key = (h, w) if key not in _FIX: g = 0.965 + 0.070*fbm(h, w, PF(26.0), 771, octaves=3) rng = np.random.RandomState(7717) cw = max(1, int(round(RS))) # one detector column, in pixels nc = max(1, int(round(w/RS))) # ...and how many there are for c in rng.choice(nc, 3, replace=False): g[:, int(c*RS):int(c*RS)+cw] *= 0.55 for c in rng.choice(nc, 6, replace=False): g[:, int(c*RS):int(c*RS)+cw] *= 0.93 _FIX[key] = g.astype(np.float32) return _FIX[key] _DUST = {} def dust(h, w): """Specks and a hair on the light box glass. Bone-dark, plate-side.""" key = (h, w) if key not in _DUST: m = np.zeros((h, w), np.float32) im = Image.fromarray(m); d = ImageDraw.Draw(im) rng = np.random.RandomState(3391) for _ in range(46): # the same 46 specks, larger x, y = rng.rand()*w, rng.rand()*h r = PF(0.6 + 1.9*rng.rand()) d.ellipse([x-r, y-r, x+r, y+r], fill=0.55+0.4*rng.rand()) pts = [(w*0.71, h*0.10)] for _ in range(9): pts.append((pts[-1][0]+PF(rng.uniform(-26, 34)), pts[-1][1]+PF(rng.uniform(4, 26)))) d.line(pts, fill=0.5, width=max(1, int(round(RS))), joint="curve") _DUST[key] = np.asarray(im, np.float32) return _DUST[key] def to_lightbox(v, plate_rect, collim, glow=1.0, box=1.0, seed=0): """Compose the exposed field onto a sheet of film on a back-lit box. Outside the collimator the film is CLEAR, so the box shines straight through it — that blazing frame around a dark image is the single most recognisable thing about looking at an X-ray, and it is free light.""" h, w = v.shape px0, py0, px1, py1 = plate_rect cx0, cy0, cx1, cy1 = collim yy, xx = np.mgrid[0:h, 0:w] inside_p = ((xx >= px0) & (xx < px1) & (yy >= py0) & (yy < py1)) inside_c = ((xx >= cx0) & (xx < cx1) & (yy >= cy0) & (yy < cy1)) field = v*inside_c clear = (inside_p & ~inside_c).astype(np.float32) outside = (~inside_p).astype(np.float32) out = field + clear*1.02*box + outside*0.34*box if glow: # ONLY the box light haloes. Letting the exposed field bloom into # itself here washes every plate to milk — the image is a shadow, it # does not emit. out = out + blur_big(clear + outside*0.8, PF(13.0))*0.14*glow return out def bloom(rgb, thr=0.62, amt=0.55, rad=13): """Exposure blooming: the brightest bone spills into its neighbourhood.""" lum = rgb.mean(2)/255.0 m = np.clip((lum-thr)/max(1e-3, 1-thr), 0, 1) src = (rgb*m[..., None]).astype(np.float32) im = Image.fromarray(np.clip(src, 0, 255).astype(np.uint8)) b = im.resize((rgb.shape[1]//3, rgb.shape[0]//3), Image.BILINEAR) b = b.filter(ImageFilter.GaussianBlur(rad/3.0)) b = b.resize((rgb.shape[1], rgb.shape[0]), Image.BILINEAR) return np.clip(rgb + np.asarray(b, np.float32)*amt, 0, 255) # ════════════════════════════════════════════════════════════════════════════ # THE ANIMAL — built as a skeleton, posed, then integrated # # World units are metres. y is UP, ground is y=0. The animal faces +x. # ════════════════════════════════════════════════════════════════════════════ def bez(p0, p1, p2, n=7): out = [] for j in range(n): u = j/(n-1.0); m = 1-u out.append((m*m*p0[0] + 2*m*u*p1[0] + u*u*p2[0], m*m*p0[1] + 2*m*u*p1[1] + u*u*p2[1])) return out def rot(p, o, a): c, s = math.cos(a), math.sin(a) dx, dy = p[0]-o[0], p[1]-o[1] return (o[0] + dx*c - dy*s, o[1] + dx*s + dy*c) class Pose: """Everything the animal is doing at time t.""" __slots__ = ("t", "dx", "walk", "ph", "breathe", "heart", "trunk", "head", "settle") def __init__(self, t, dx=0.0, walk=0.0, breathe=0.0, heart=0.0, trunk=0.0, head=0.0, settle=0.0): self.t = t; self.dx = dx; self.walk = walk self.ph = t/(BEAT*2.0) # one stride per two beats self.breathe = breathe; self.heart = heart self.trunk = trunk; self.head = head; self.settle = settle def leg(P, hip, l_up, l_lo, a_up, a_lo, r_up, r_lo, rho, foot=True): """Two-segment FK leg + foot. Angles are from straight-down, +ve forward.""" kx = hip[0] + math.sin(a_up)*l_up ky = hip[1] - math.cos(a_up)*l_up fx = kx + math.sin(a_up+a_lo)*l_lo fy = ky - math.cos(a_up+a_lo)*l_lo P.cap(hip, (kx, ky), r_up, r_up*0.72, rho) P.ell(kx, ky, r_up*0.95, r_up*0.95, r_up*0.9, rho*0.9) # joint P.cap((kx, ky), (fx, fy), r_lo, r_lo*0.86, rho) if foot: P.ell(fx, fy-0.06, r_lo*1.5, r_lo*0.95, r_lo*1.2, rho*0.95) for j in range(4): # toe bones ox = (j-1.5)*r_lo*0.85 P.cap((fx+ox*0.5, fy-0.10), (fx+ox, max(0.045, fy-0.30)), r_lo*0.36, r_lo*0.26, rho) return (fx, fy) def mammoth(P, p, parts="all", ice=0.0, hidef=1.0): """Integrate the whole animal into the plate. `parts` selects a subset so close-ups don't pay for the ribcage they can't see.""" D = p.dx all_ = parts == "all" want = (lambda k: all_ or k in parts) bob = math.sin(p.ph*math.tau)*0.055*p.walk sh_y = 3.30 + bob - p.settle*0.35 hd = p.head # ---- soft tissue first (it sits under everything) --------------------- if want("soft"): br = 1.0 + 0.045*math.sin(p.t*0.9) * (0.4+p.breathe) P.ell(D-0.35, 2.20, 1.62, 1.02*br, 1.05, RHO_MUSC*0.85) P.ell(D+1.05, 2.55, 0.85, 0.80, 0.80, RHO_MUSC*0.7) P.ell(D-1.45, 2.30, 0.75, 0.72, 0.70, RHO_MUSC*0.7) # the shaggy coat — barely there, but it is the silhouette if hidef > 0.4: P.ell(D-0.35, 2.15, 2.05, 1.42, 1.25, RHO_FUR) P.ell(D-0.35, 1.35, 1.95, 0.95, 1.15, RHO_FUR*1.4) # ---- spine + hump ----------------------------------------------------- spine = bez((D-1.95, 2.58), (D-0.55, 3.02+bob), (D+1.28, sh_y-0.10), n=13) if want("spine"): # thin cord, fat vertebrae: a uniform sausage reads as a bone, a # beaded chain reads as a SPINE P.arc(spine, [0.070]*len(spine), RHO_BONE*0.80) for j, (sx, sy) in enumerate(spine): u = j/(len(spine)-1.0) hgt = 0.12 + 0.78*math.exp(-((u-0.80)**2)/0.050) # the hump P.cap((sx, sy), (sx-0.04, sy+hgt), 0.080, 0.026, RHO_BONE*1.05) P.ell(sx, sy, 0.118, 0.105, 0.115, RHO_BONE*0.62) # tail tl = bez((D-1.95, 2.58), (D-2.42, 2.20), (D-2.30, 1.35), n=6) P.arc(tl, [0.075, 0.062, 0.050, 0.040, 0.032, 0.024], RHO_BONE*0.9) if want("pelvis"): P.ell(D-1.62, 2.42, 0.44, 0.36, 0.30, RHO_BONE*0.85) P.cap((D-1.72, 2.66), (D-1.30, 2.30), 0.20, 0.12, RHO_BONE*0.7, flat=0.35) # ---- ribcage ---------------------------------------------------------- if want("ribs"): for i in range(11): u = i/10.0 ax = D + 1.05 - u*1.90 ay = 2.90 + 0.24*(1-u) + bob*(1-u*0.4) spread = 1.30 - 0.42*abs(u-0.35) drop = 1.02 - 0.10*u br = 1.0 + 0.075*math.sin(p.t*0.9 + u*0.6)*(0.3+p.breathe) for sgn, dmul in ((1, 1.0), (-1, 0.80)): # near + far side pts = bez((ax, ay), (ax - 0.30 + sgn*0.10, ay - spread*0.55*br), (ax + 0.16 + sgn*0.06, ay - drop*1.02*br), n=8) rr = [0.075, 0.070, 0.064, 0.058, 0.052, 0.046, 0.040, 0.033] P.arc(pts, rr, RHO_BONE*dmul) # sternum P.cap((D+0.85, 1.78), (D-0.55, 1.86), 0.10, 0.075, RHO_BONE*0.75) # ---- the heart -------------------------------------------------------- if want("heart"): hh = p.heart pulse = 1.0 + 0.16*hh*math.sin(p.t*math.tau*0.72) P.ell(D+0.42, 2.05, 0.44*pulse, 0.50*pulse, 0.42, RHO_MUSC + (RHO_CART*2.4)*hh) if hh > 0.05: # great vessels lighting up for (vx, vy, vr) in ((0.42, 2.58, 0.11), (0.10, 2.30, 0.09), (0.72, 2.34, 0.09)): P.cap((D+0.42, 2.20), (D+vx*1.15, vy), vr, vr*0.6, RHO_CART*1.6*hh) # ---- legs ------------------------------------------------------------- if want("legs"): s1 = math.sin(p.ph*math.tau); s2 = math.sin(p.ph*math.tau + math.pi) # column legs, and they must actually REACH THE GROUND (y=0): front # hip sits at ~3.0, so the two segments have to sum to ~2.85 for (hx, hy, up, lo, ru, rl, phs, far) in ( (D+1.02, sh_y-0.28, 1.54, 1.30, 0.200, 0.150, s1, 1.0), (D+0.84, sh_y-0.30, 1.52, 1.28, 0.190, 0.142, s2, 0.70), (D-1.42, 2.36, 1.28, 1.02, 0.195, 0.148, s2, 1.0), (D-1.24, 2.34, 1.26, 1.00, 0.185, 0.140, s1, 0.70)): a_up = phs*0.36*p.walk a_lo = -abs(phs)*0.30*p.walk - 0.06 leg(P, (hx, hy), up, lo, a_up, a_lo, ru, rl, RHO_BONE*far, foot=True) # scapula blades P.cap((D+1.20, sh_y+0.30), (D+0.92, sh_y-0.32), 0.30, 0.14, RHO_BONE*0.75, flat=0.30) P.cap((D+1.02, sh_y+0.28), (D+0.80, sh_y-0.30), 0.26, 0.12, RHO_BONE*0.55, flat=0.30) # ---- skull ------------------------------------------------------------ if want("skull"): kx, ky = D+1.92, sh_y+0.02 + hd*0.10 P.ell(kx, ky, 0.44, 0.46, 0.42, RHO_BONE*0.80) # braincase P.ell(kx-0.03, ky+0.40, 0.38, 0.34, 0.36, RHO_BONE*0.60) # the dome P.ell(kx+0.10, ky-0.44, 0.32, 0.34, 0.30, RHO_BONE*0.95) # maxilla P.cap((kx-0.42, ky-0.02), (D+1.34, sh_y-0.06), 0.22, 0.24, RHO_BONE*0.7) # occiput P.ell(kx+0.16, ky+0.06, 0.13, 0.15, 0.16, -RHO_BONE*0.72) # orbit P.ell(kx-0.03, ky+0.30, 0.20, 0.16, 0.22, -RHO_BONE*0.34) # air sinus P.ell(kx+0.24, ky-0.20, 0.09, 0.10, 0.10, -RHO_BONE*0.5) # nares P.cap((kx+0.02, ky-0.70), (kx-0.30, ky-0.86), 0.16, 0.11, RHO_BONE*0.8) # mandible # tusks — ivory, off the scale, and they spiral for sgn, dm in ((1, 1.0), (-1, 0.86)): base = (kx+0.22+sgn*0.05, ky-0.52) pts = bez(base, (kx+1.30, ky-1.55-sgn*0.05), (kx+2.05, ky-0.70+sgn*0.10), n=9) pts += bez((kx+2.05, ky-0.70+sgn*0.10), (kx+2.32, ky-0.05), (kx+1.70, ky+0.22), n=6)[1:] rr = np.linspace(0.165, 0.038, len(pts)).tolist() P.arc(pts, rr, RHO_IVORY*dm) P.ell(base[0], base[1], 0.19, 0.16, 0.17, RHO_IVORY*0.55*dm) # trunk: soft tissue, with the cartilage rings faintly stacked if want("trunk") or all_: tp = bez((kx+0.22, ky-0.62), (kx+0.95, ky-1.60), (kx+0.80, ky-2.35 + p.trunk*0.75), n=9) tp += bez(tp[-1], (kx+1.35, ky-2.95 + p.trunk*1.1), (kx+1.72, ky-2.35 + p.trunk*1.4), n=6)[1:] rr = np.linspace(0.32, 0.12, len(tp)).tolist() P.arc(tp, rr, RHO_MUSC*2.2) # the ring cartilage: a stack of faint discs is the ONLY thing # that makes a trunk legible on a plate, since it has no bone for j in range(len(tp)): rj = rr[j] P.ell(tp[j][0], tp[j][1], rj*1.0, rj*0.30, rj*0.95, RHO_CART*1.35) # ear (mammoth ears are small — cold is cold) P.ell(kx-0.30, ky+0.06, 0.22, 0.28, 0.05, RHO_HIDE*1.6) def ice_block(P, melt, seed=101, cracks=0.0, cx=0.0): """Permafrost block 7.""" P.slab(cx-3.2, cx+4.4, -0.15, 4.6, 1.85, RHO_ICE, seed, scale=1.05, soft=0.30, melt=melt) P.slab(cx-3.0, cx+4.2, -0.10, 4.4, 0.80, RHO_ICE*0.8, seed+9, scale=0.42, soft=0.26, melt=min(1.0, melt*1.25), octaves=3) # the specimen table rail — steel, and the only straight line in the film P.cap((cx-3.6, -0.16), (cx+4.8, -0.16), 0.055, 0.055, RHO_BONE*1.15) P.cap((cx-3.6, -0.34), (cx+4.8, -0.34), 0.030, 0.030, RHO_BONE*0.75) # meltwater in the tray — it never fully leaves, and it is the last thing # on the last plate P.slab(cx-3.1, cx+4.3, -0.13, 0.10 + 0.34*melt, 0.30 + 0.55*melt, RHO_MUSC*1.25, seed+21, scale=0.55, soft=0.18, octaves=3) if cracks > 0.01: rng = np.random.RandomState(seed+55) nc = int(3 + 12*cracks) for j in range(nc): x = cx - 3.0 + 7.2*rng.rand() pts = [(x, 4.4)] for _ in range(7): pts.append((pts[-1][0] + rng.uniform(-0.42, 0.42), pts[-1][1] - rng.uniform(0.35, 0.80))) P.crack(pts, 0.024*(0.5+cracks), RHO_ICE*3.6, seed=seed+j) # ════════════════════════════════════════════════════════════════════════════ # CHROME — plate numbers, lead markers, dictation, calipers. # Drawn crisply AFTER the grade so nothing fringes the text (AESTHETIC 13b). # ════════════════════════════════════════════════════════════════════════════ # ── portable font resolution (cross-platform; replaces the repo-only lookup) ── import warnings as _warnings _FONT_ALIASES = { "Menlo.ttc": ["Menlo.ttc", "DejaVuSansMono.ttf", "consola.ttf", "LiberationMono-Regular.ttf"], "Georgia.ttf": ["Georgia.ttf", "georgia.ttf", "DejaVuSerif.ttf", "LiberationSerif-Regular.ttf"], "Georgia Bold.ttf": ["Georgia Bold.ttf", "georgiab.ttf", "DejaVuSerif-Bold.ttf", "LiberationSerif-Bold.ttf"], "Georgia Italic.ttf": ["Georgia Italic.ttf", "georgiai.ttf", "DejaVuSerif-Italic.ttf", "LiberationSerif-Italic.ttf"], "Impact.ttf": ["Impact.ttf", "impact.ttf", "Anton-Regular.ttf", "DejaVuSans-Bold.ttf"], "Helvetica.ttc": ["Helvetica.ttc", "arial.ttf", "Arial.ttf", "DejaVuSans.ttf", "LiberationSans-Regular.ttf"], } def _font_dirs(): here = Path(__file__).resolve() dirs = [here.parent / "fonts"] + [p / "fonts" for p in list(here.parents)[1:4]] try: home = Path.home() except Exception: home = None dirs += [Path("/System/Library/Fonts"), Path("/System/Library/Fonts/Supplemental"), Path("/Library/Fonts"), Path("C:/Windows/Fonts"), Path("/usr/share/fonts"), Path("/usr/local/share/fonts")] if home: dirs += [home / "Library/Fonts", home / ".fonts", home / ".local/share/fonts"] return dirs _FONT_DIRS = _font_dirs() _FF = {} def _find_font(name): """Path of a usable font file for `name`, or None. Cached per name.""" if name in _FF: return _FF[name] found = None for cand in _FONT_ALIASES.get(name, [name]): for d in _FONT_DIRS: if not d.is_dir(): continue p = d / cand if p.is_file(): found = p; break try: found = next(iter(d.rglob(cand)), None) except OSError: found = None if found: break if found: break if found is None: _warnings.warn(f"font {name} not found in fonts/ or system font dirs; " f"using Pillow default (layout will differ)") _FF[name] = found return found def _load_font(p, size): """ImageFont for path `p` (from _find_font) at `size`; Pillow default if p is None.""" if p is None: try: return ImageFont.load_default(size=int(size)) except TypeError: return ImageFont.load_default() return ImageFont.truetype(str(p), size) _FC = {} def font(size, name="Menlo.ttc"): key = (size, name) if key not in _FC: p = _find_font(name) # Scaled ONCE, here. Call sites always pass authoring sizes. _FC[key] = _load_font(p, max(1, PX(size))) return _FC[key] INK = (196, 214, 232) INK_D = (120, 146, 172) HOT = (250, 253, 255) WARN = (255, 208, 150) def tc(sec): h = int(sec//3600); m = int(sec//60) % 60; s = sec % 60 return f"T+{h:02d}:{m:02d}:{s:05.2f}" # ════════════════════════════════════════════════════════════════════════════ # SHOTS # ════════════════════════════════════════════════════════════════════════════ class Shot: __slots__ = ("idx", "i0", "i1", "n", "engine", "section", "seed", "plate", "kvp", "mas", "wipe", "card") def __init__(self, idx, i0, i1, engine, section, plate, kvp, mas, wipe, card=None): self.idx, self.i0, self.i1 = idx, i0, i1 self.n = i1-i0 self.engine, self.section = engine, section self.seed = 62000 + idx*7919 self.plate, self.kvp, self.mas, self.wipe = plate, kvp, mas, wipe self.card = card def technique_at(t): """The lab's exposure technique, which starts timid. Plate 1 is thin because nobody has shot a mammoth before.""" return float(0.58 + 0.42*np.clip(t/(2.4*BAR), 0, 1)**0.8) def melt_at(t): """Ice retreat over the whole piece — the one monotone clock.""" u = np.clip((t - 1.2*BAR)/(9.4*BAR), 0, 1) return float(u**0.78) def heart_at(t): b = t/BAR if b < 10.4: return 0.0 if b < 12.0: return float((b-10.4)/1.6)*0.55 return float(min(1.0, 0.55 + (b-12.0)/1.4*0.45)) def walk_at(t): b = t/BAR return float(np.clip((b-13.25)/0.70, 0, 1)) def leave_dx(t): """How far the animal has moved out of the beam.""" b = t/BAR u = np.clip((b-13.60)/2.30, 0, 1) return float(10.8*(u**1.45)) BASE_CAM = (0.35, 2.15, 128.0) class Engine: """A shot. Sets a camera, integrates density, develops a plate.""" CAM = BASE_CAM # VARY = (world dx, world dy, +/- fraction of scale) drawn per shot. An # engine used five times must not be the same picture five times. VARY = (0.10, 0.05, 0.06) PARTS = "all" ICE = True def __init__(self, shot, rng): self.shot = shot; self.rng = rng self.i0 = shot.i0 self.j = float(rng.random()) self.drift = (float(rng.random())-0.5, float(rng.random())-0.5) def cam(self, t, u, e): cx, cy, s = self.CAM vx, vy, vs = self.VARY cx += (self.j-0.5)*2.0*vx cy += self.drift[1]*2.0*vy s *= 1.0 + self.drift[0]*2.0*vs # slow radiographic table drift + a push driven by the low end k = 1.0 + 0.035*u + 0.05*e["sub"] return (cx + 0.05*math.sin(t*0.35+self.j*6), cy + 0.03*math.cos(t*0.29+self.j*4), s*k) def build(self, P, t, u, e): # +x is the way it faces, so leaving is +leave_dx for EVERY engine. # (The block stays where it is; only the animal goes.) p = Pose(t, dx=leave_dx(t), walk=walk_at(t), breathe=float(e["mid"]), heart=heart_at(t), trunk=0.35*math.sin(t*0.5), head=0.1*math.sin(t*0.4)) if self.ICE: ice_block(P, melt_at(t), cracks=np.clip((t/BAR-4.6)/2.2, 0, 1)) mammoth(P, p, parts=self.PARTS) def frame(self, k, u, e): t = (self.i0+k)/FPS P = Plate(W, H, self.cam(t, u, e)) self.build(P, t, u, e) return P class Slab(Engine): """Wide. The block on the table. Almost nothing legible yet.""" CAM = (0.35, 2.10, 118.0) VARY = (0.22, 0.10, 0.07) class Full(Engine): CAM = (0.20, 2.20, 132.0) VARY = (0.32, 0.14, 0.10) class Skull(Engine): CAM = (2.30, 3.05, 268.0) VARY = (0.30, 0.18, 0.17) PARTS = ("skull", "spine", "soft", "trunk") class Tusk(Engine): """Ivory. Three of these in the film, so the framing has to move.""" CAM = (3.50, 2.55, 300.0) VARY = (0.62, 0.34, 0.30) PARTS = ("skull", "soft") class Ribs(Engine): CAM = (0.30, 2.35, 250.0) VARY = (0.52, 0.24, 0.22) PARTS = ("ribs", "spine", "soft", "heart") class Foot(Engine): """Both near feet and the steel rail. A single leg is a stick.""" CAM = (-0.10, 0.80, 205.0) VARY = (0.90, 0.26, 0.26) PARTS = ("legs", "soft", "ribs") class Trunk(Engine): """The one part of a mammoth with no bone in it at all.""" CAM = (3.48, 1.72, 232.0) VARY = (0.30, 0.22, 0.16) PARTS = ("skull", "soft", "trunk") class Heart(Engine): """The money shot: soft tissue between the ribs, taking on contrast.""" CAM = (0.55, 2.10, 400.0) VARY = (0.22, 0.14, 0.14) PARTS = ("ribs", "heart", "soft") class Walk(Engine): """It leaves. The camera does not follow.""" CAM = (1.60, 2.05, 108.0) VARY = (0.28, 0.12, 0.08) ICE = True def build(self, P, t, u, e): p = Pose(t, dx=-1.2+leave_dx(t), walk=max(0.55, walk_at(t)), breathe=float(e["mid"]), heart=heart_at(t), trunk=0.5*math.sin(t*0.8), head=0.18*math.sin(t*0.9)) ice_block(P, min(1.0, melt_at(t)*1.35), cracks=1.0) mammoth(P, p, parts="all") class Stack(Engine): """The lab's archive: six plates on one box, each a different exposure, each a different region. This is the shot that says 'they kept plating'.""" ICE = True GRID = [(0.35, 2.05, 52.0, "all"), (2.30, 3.05, 110.0, ("skull", "spine")), (3.55, 2.60, 132.0, ("skull",)), (0.30, 2.35, 100.0, ("ribs", "spine")), (1.05, 0.75, 150.0, ("legs",)), (0.55, 2.10, 160.0, ("ribs", "heart"))] def frame(self, k, u, e): t = (self.i0+k)/FPS tiles = [] cols, rows = 3, 2 tw, th = W//cols, H//rows for gi, (cx, cy, s, parts) in enumerate(self.GRID): P = Plate(tw-PX(16), th-PX(16), (cx, cy + 0.02*math.sin(t+gi), s)) p = Pose(t, dx=-leave_dx(t), walk=walk_at(t), breathe=float(e["mid"]), heart=heart_at(t), trunk=0.3*math.sin(t*0.5)) ice_block(P, melt_at(t)*(0.6+0.12*gi), cracks=np.clip((t/BAR-4.6)/2.2, 0, 1)) mammoth(P, p, parts=parts, hidef=0.3) tiles.append((gi, P)) out = Plate(W, H, BASE_CAM) for gi, P in tiles: cxx, cyy = (gi % cols)*tw+PX(8), (gi//cols)*th+PX(8) out.D[cyy:cyy+P.h, cxx:cxx+P.w] = P.D return out ENGINES = {"slab": Slab, "full": Full, "skull": Skull, "tusk": Tusk, "ribs": Ribs, "foot": Foot, "trunk": Trunk, "heart": Heart, "walk": Walk, "stack": Stack} PLAN = { "plate": (["slab", "slab", "skull", "stack", "slab", "skull"], [8, 6, 6, 5]), "thaw": (["skull", "tusk", "ribs", "stack", "foot", "slab", "ribs", "tusk"], [4, 3, 6, 4]), "fracture": (["ribs", "tusk", "full", "skull", "trunk"], [3, 4, 2, 6]), "sludge": (["full", "ribs", "skull", "tusk", "stack", "trunk", "ribs", "foot", "full", "tusk"], [2, 2, 3, 4]), "restart": (["heart", "ribs", "heart", "full", "heart"], [3, 4, 2, 6]), "walk": (["walk", "full", "walk", "walk", "ribs", "full"], [4, 6, 8]), } CARDS = {"plate": "MAMMOTH · THAW", "walk": None} def build_shots(): R = np.random.RandomState(1791) shots = []; idx = 0; last = []; plate_no = 1 for nm, b0, b1 in SECTIONS: engs, menu = PLAN[nm] t = b0*BAR; first = True while t < b1*BAR - 1e-6: step = menu[R.randint(len(menu))]*(BEAT/2) t2 = min(t+step, b1*BAR) if (b1*BAR - t2) < BEAT*0.9: t2 = b1*BAR i0, i1 = int(t*FPS), int(t2*FPS) if i1 > i0: # no engine twice within three cuts — one repeat reads as # a mistake, two in a row reads as a still frame pool = ([x for x in engs if x not in last[-2:]] or [x for x in engs if x not in last[-1:]] or list(engs)) eng = "slab" if idx == 0 else pool[R.randint(len(pool))] last.append(eng) # exposure varies plate to plate; two are deliberately wrong kvp = float(np.clip(R.normal(1.0, 0.14), 0.70, 1.32)) mas = float(np.clip(R.normal(1.0, 0.14), 0.72, 1.26)) if nm == "walk": # chasing it off the plate mas = float(np.clip(mas*1.22, 1.02, 1.42)) if idx == 9: kvp, mas = 1.40, 0.64 # under-cooked plate if idx == 22: kvp, mas = 0.74, 1.30 # blown out shots.append(Shot(idx, i0, i1, eng, nm, plate_no, kvp, mas, wipe=(R.rand() < 0.72), card=CARDS.get(nm) if first else None)) idx += 1; plate_no += 1 + (R.rand() < 0.35) first = False t = t2 if shots: shots[-1].i1 = N_FRAMES; shots[-1].n = N_FRAMES - shots[-1].i0 return shots # ════════════════════════════════════════════════════════════════════════════ # POST — tint -> vignette -> grain -> (chrome) -> letterbox # ════════════════════════════════════════════════════════════════════════════ _VIG = {} def lightbox_vig(): if "v" not in _VIG: yy, xx = np.mgrid[0:H, 0:W] nx = (xx-W/2)/(W/2); ny = (yy-H/2)/(H/2) r = np.sqrt(nx*nx + ny*ny)/1.42 v = np.clip(1.0 - 0.52*r**2.0, 0, 1) # the two fluorescent tubes behind the box tube = (0.06*np.exp(-((ny-0.42)**2)/0.012) + 0.06*np.exp(-((ny+0.42)**2)/0.012)) _VIG["v"] = (v + tube).astype(np.float32)[..., None] return _VIG["v"] def geometry(shot, i, jit=0): """Collimator + film rectangles for this plate. Every plate is loaded a little crooked, which is most of why it reads as a physical object; `jit` is the sheet shifting in its clips when the low end hits.""" r = np.random.RandomState(shot.seed ^ 0x5EED) mx = 26 + int(r.rand()*16); my = 20 + int(r.rand()*14) dy = jit; dx = int(jit*0.4) px0, py0, px1, py1 = mx+dx, my+dy, AW-mx+dx, AH-my+dy ci = 22 + int(r.rand()*26) return (px0, py0, px1, py1), (px0+ci, py0+ci, px1-ci, py1-ci) def post(P, i, e, shot): """P is a Plate (density). Everything from density to delivered frame.""" t = i/FPS k = i - shot.i0 u = k/max(1, shot.n-1) # the rects come back in authoring units — the chrome draws in those, # the numpy passes below need them in real pixels plate_rect, collim = geometry(shot, i, jit=int(round(5.0*float(e["hit"]) - 2.0*float(e["sub"])))) PR = tuple(PX(q) for q in plate_rect) CO = tuple(PX(q) for q in collim) # --- exposure: the plate is being made right now ----------------------- wipe_n = int(FPS*0.34) exposing = shot.wipe and k < wipe_n mas = shot.mas*technique_at(t)*(1.0 + 0.30*float(e["hit"]) + 0.16*float(e["sub"])) v = develop(P.D, kvp=shot.kvp, mas=mas, focal=PF(1.05 + 0.18*float(e["high"])), seed=4400+i, mottle=1.0, scatter=(1.0 - melt_at(t))**1.25*0.92) v *= fixed_pattern(H, W) v -= dust(H, W)*0.5 v = to_lightbox(v, PR, CO, glow=1.0, box=1.0, seed=i) if exposing: # the scan bar: below it the plate is still clear, above it exposed wu = k/max(1, wipe_n-1) yb = int(CO[1] + (CO[3]-CO[1])*wu) v[yb:CO[3], CO[0]:CO[2]] = 1.28 yl0, yl1 = max(0, yb-PX(9)), min(H, yb+PX(4)) v[yl0:yl1, CO[0]:CO[2]] = 1.55 rgb = colorize(np.clip(v/1.30, 0, 1)) rgb = bloom(rgb, thr=0.74, amt=0.26 + 0.34*float(e["sub"]), rad=PF(15)) # ---- 1. tint: cold plate, faintly warm light box ---------------------- lum = rgb.mean(2, keepdims=True)/255.0 rgb = rgb + (1-lum)*np.array([-10, 0, 16], np.float32) \ + lum*np.array([8, 4, -6], np.float32) # ---- 2. vignette ------------------------------------------------------ rgb = rgb*lightbox_vig() # ---- 3. grain --------------------------------------------------------- rng = np.random.RandomState(9100+i) if RS == 1.0: rgb = rgb + rng.normal(0, 2.6, (H, W, 1)).astype(np.float32) else: # grain is a look, not a resolution gn = rng.normal(0, 2.6, (AH, AW))*8.0 + 128.0 gi = Image.fromarray(np.clip(gn, 0, 255).astype(np.uint8)) rgb = rgb + ((np.asarray(gi.resize((W, H), Image.NEAREST), np.float32) - 128.0)/8.0)[..., None] rgb = np.clip(rgb, 0, 255) img = Image.fromarray(rgb.astype(np.uint8)) chrome(img, i, t, e, shot, plate_rect, collim, u) # ---- 4. letterbox ----------------------------------------------------- d = ImageDraw.Draw(img) bh = int(H*0.034) d.rectangle([0, 0, W, bh], fill=(0, 0, 0)) d.rectangle([0, H-bh, W, H], fill=(0, 0, 0)) return img def chrome(img, i, t, e, shot, plate_rect, collim, u): """Burned-in plate annotation. Everything lives INSIDE the collimated field and on its own alpha layer, for two reasons: the clear film outside the collimator is the brightest thing in the frame (pale text there is unreadable and the letterbox eats it), and a real radiograph burns its annotation into the exposed area. Composited crisply after the grade — the picture may fringe, the words never do (AESTHETIC 13b).""" cx0, cy0, cx1, cy1 = collim f9, f11, f13, f16 = font(10), font(12), font(14), font(17) bar = t/BAR hb = heart_at(t) ov = Image.new("RGBA", img.size, (0, 0, 0, 0)) d = mkdraw(ov) # every coordinate below is authoring-space SCRIM = (0, 3, 9, 122) # --- the exposed field's border (the collimator edge) ------------------ d.rectangle([cx0, cy0, cx1-1, cy1-1], outline=(40, 58, 78, 255), width=1) # --- plate identity, top left ----------------------------------------- d.rectangle([cx0+1, cy0+1, cx0+248, cy0+64], fill=SCRIM) d.text((cx0+12, cy0+8), f"PLATE {shot.plate:03d}", font=f16, fill=INK+(255,)) d.text((cx0+12, cy0+31), f"{int(56+shot.kvp*30):02d} kVp {shot.mas*8.0:4.1f} mAs", font=f11, fill=INK_D+(255,)) d.text((cx0+12, cy0+46), "AP · GRID 12:1 · SID 110", font=f9, fill=INK_D+(255,)) # --- specimen, top right ---------------------------------------------- d.rectangle([cx1-266, cy0+1, cx1-1, cy0+58], fill=SCRIM) for j, (s, ff, cc) in enumerate(( ("PERMAFROST BLOCK 7", f11, INK), ("SPEC. M-001 M. PRIMIGENIUS", f9, INK_D), ("N 71°02' E 152°18'", f9, INK_D))): d.text((cx1-12-TL(d, s, ff), cy0+9+j*16), s, font=ff, fill=cc+(255,)) # --- lead marker: lead is opaque, so it burns white -------------------- mk = "L" if (shot.idx % 2 == 0) else "R" mx, my = cx1-46, cy0+70 d.rectangle([mx-6, my-4, mx+26, my+30], fill=(14, 26, 40, 235)) d.text((mx, my), mk, font=font(24), fill=HOT+(255,)) # --- clock + the two clocks that matter, bottom left ------------------ d.rectangle([cx0+1, cy1-52, cx0+318, cy1-1], fill=SCRIM) d.text((cx0+12, cy1-46), tc(3600*3 + t*470), font=f13, fill=INK+(255,)) d.text((cx0+12, cy1-28), f"ICE {max(0.0, 148*(1-melt_at(t))):5.1f} cm " f"CORE {-19.4 + 21.0*melt_at(t):+5.1f} °C", font=f11, fill=INK_D+(255,)) # --- the step wedge, bottom right ------------------------------------- sw, sh2 = 132, 9 sx, sy = cx1-16-sw, cy1-40 d.rectangle([sx-10, cy1-50, cx1-1, cy1-1], fill=SCRIM) for q in range(11): c = tuple(int(x) for x in PLATE_LUT[int(q/10.0*(len(PLATE_LUT)-1))]) d.rectangle([sx+q*sw//11, sy, sx+(q+1)*sw//11-1, sy+sh2], fill=c+(255,)) d.rectangle([sx, sy, sx+sw, sy+sh2], outline=(60, 82, 106, 255)) d.text((sx, sy+13), "OD 0.2 ——————— 3.2", font=f9, fill=INK_D+(255,)) # --- stack engine: the archive's individual plate frames -------------- if shot.engine == "stack": tw, th = AW//3, AH//2 for gi in range(6): gx, gy = (gi % 3)*tw+8, (gi//3)*th+8 d.rectangle([gx, gy, gx+tw-17, gy+th-17], outline=(58, 84, 108, 255)) d.rectangle([gx+1, gy+1, gx+64, gy+19], fill=SCRIM) d.text((gx+7, gy+4), f"{shot.plate+gi:03d}", font=f9, fill=INK_D+(255,)) # --- the caliper, on ivory -------------------------------------------- if shot.engine == "tusk": r = np.random.RandomState(shot.seed ^ 0xCA11) ax, ay = cx0+90+int(r.rand()*70), cy0+150+int(r.rand()*90) bx, by = cx1-130-int(r.rand()*80), cy0+96+int(r.rand()*70) d.line([ax, ay, bx, by], fill=(110, 146, 182, 255), width=1) for (ex, ey) in ((ax, ay), (bx, by)): d.line([ex, ey-9, ex, ey+9], fill=(132, 168, 202, 255), width=2) lbl = f"Ø {19.0+r.rand()*5:4.1f} cm L {3.30+r.rand()*0.7:4.2f} m" lx, ly = (ax+bx)//2 - 60, (ay+by)//2 - 26 d.rectangle([lx-6, ly-3, lx+156, ly+18], fill=SCRIM) d.text((lx, ly), lbl, font=f11, fill=INK+(255,)) # --- the ECG. Wherever the heart is, not wherever the camera is. ------ if hb > 0.001 or bar > 9.6: hy = cy1-108 pts = [] for q in range(0, cx1-cx0-40, 3): ph = ((t*0.62 + q/260.0) % 1.0) spike = 0.0 if hb > 0.02: spike = (math.exp(-((ph-0.30)**2)/0.00035) - 0.30*math.exp(-((ph-0.26)**2)/0.0006) + 0.34*math.exp(-((ph-0.46)**2)/0.004)) wob = math.sin(q*0.21 + t*3.1)*0.05 pts.append((cx0+20+q, hy - (spike*hb + wob*(0.4+hb))*46)) live = hb > 0.05 col = (150, 255, 196, 255) if live else (74, 100, 118, 255) d.line(pts, fill=col, width=2, joint="curve") lbl = f"ECG {int(3+hb*9):02d} bpm" if live else "ECG —— ASYSTOLE" lw = TL(d, lbl, f11) d.rectangle([cx1-28-lw, hy+8, cx1-8, hy+28], fill=SCRIM) d.text((cx1-20-lw, hy+10), lbl, font=f11, fill=col) # --- radiologist dictation, transcribed onto the plate ---------------- for (vb, _txt, a, _g) in VO: vlen = 1.35 + 0.055*len(_txt) if not (vb <= bar < vb + vlen/BAR + 0.42): continue al = max(0.0, min(1.0, min((bar-vb)*7.0, (vb + vlen/BAR + 0.42 - bar)*6.0))) if al < 0.02: continue wdt = TL(d, "> " + a, f13) d.rectangle([cx0+12, cy1-92, cx0+34+wdt, cy1-64], fill=(2, 7, 16, int(150*al))) d.line([cx0+16, cy1-90, cx0+16, cy1-66], fill=WARN+(int(255*al),), width=2) d.text((cx0+26, cy1-86), "> " + a, font=f13, fill=WARN+(int(255*al),)) # --- title card -------------------------------------------------------- if shot.card: kk = i - shot.i0 life = FPS*3.4 if kk < life: al = max(0.0, min(1.0, min(kk/8.0, (life-kk)/16.0))) # THE TITLE MOMENT — struck in the plate's own annotation ink, # with the show name below the series line. f = font(46, "Georgia Bold.ttf") lw = TL(d, shot.card, f) d.text((AW/2-lw/2, AH*0.455), shot.card, font=f, fill=(232, 242, 252, int(255*al))) s2 = "RADIOGRAPHIC SERIES · 16 h · 241 PLATES" lw2 = TL(d, s2, f13) d.rectangle([AW/2-lw2/2-14, AH*0.455+74, AW/2+lw2/2+14, AH*0.455+98], fill=(2, 7, 16, int(140*al))) d.text((AW/2-lw2/2, AH*0.455+78), s2, font=f13, fill=INK_D+(int(255*al),)) fs = font(13, "Menlo.ttc") tr = 6.5 sw3 = sum(TL(d, ch, fs) + tr for ch in SUBT) - tr xs3 = AW/2 - sw3/2; ys3 = AH*0.455 + 112 d.line([AW/2-sw3/2-16, ys3-9, AW/2+sw3/2+16, ys3-9], fill=INK_D+(int(150*al),), width=1) for ch in SUBT: d.text((xs3, ys3), ch, font=fs, fill=HOT+(int(205*al),)) xs3 += TL(d, ch, fs) + tr # --- the last words ---------------------------------------------------- if bar > 15.72: al = max(0.0, min(1.0, min((bar-15.72)*4.4, (16.66-bar)*3.2))) if al > 0.01: f = font(21, "Georgia Bold.ttf") s2 = "PLATE 241 · FIELD EMPTY" lw = TL(d, s2, f) d.text((AW/2-lw/2, AH*0.50), s2, font=f, fill=(224, 238, 250, int(255*al))) img.paste(ov, (0, 0), ov) def render_shot(job): shot, force = job E = env() rng = np.random.default_rng(shot.seed) eng = ENGINES[shot.engine](shot, rng) made = 0 for k in range(shot.n): i = shot.i0 + k p = FRAMES/f"f{i:05d}.png" if p.exists() and not force: continue e = {kk: float(E[kk][min(i, N_FRAMES-1)]) for kk in E} u = k/max(1, shot.n-1) P = eng.frame(k, u, e) post(P, i, e, shot).save(p, compress_level=1) made += 1 return f"shot {shot.idx:02d} {shot.engine:6s} {shot.section:8s} {made}/{shot.n}" def contact_sheet(shots, frac=0.5): import time cols = 6; rows = (len(shots)+cols-1)//cols tw, th = PX(320), PX(180) sheet = Image.new("RGB", (cols*tw, rows*(th+PX(26))), (8, 10, 14)) sd = ImageDraw.Draw(sheet) E = env(); t0 = time.time() for n, sh in enumerate(shots): st = time.time() rng = np.random.default_rng(sh.seed) eng = ENGINES[sh.engine](sh, rng) kk = max(0, int(sh.n*frac)) i = sh.i0+kk e = {c: float(E[c][min(i, N_FRAMES-1)]) for c in E} P = eng.frame(kk, kk/max(1, sh.n-1), e) im = post(P, i, e, sh).resize((tw, th), Image.LANCZOS) ms = (time.time()-st)*1000 cx, cy = (n % cols)*tw, (n//cols)*(th+PX(26)) sheet.paste(im, (cx, cy)) note = (f"{sh.idx:02d} {sh.engine[:5]} {sh.section[:4]} {sh.i0/FPS:.0f}s " f"{sh.n}f kv{sh.kvp:.2f} ma{sh.mas:.2f} {ms:.0f}ms") sd.text((cx+PX(5), cy+th+PX(5)), note, font=font(12), fill=(190, 200, 214)) print(note, flush=True) p = OUT/"contact_sheet.png"; sheet.save(p) print(f"contact sheet -> {p} ({len(shots)} shots, {time.time()-t0:.1f}s)") def main(): ap = argparse.ArgumentParser() ap.add_argument("--sheet", action="store_true") ap.add_argument("--sheet-at", type=float, default=0.5) ap.add_argument("--shots", default="") ap.add_argument("--force", action="store_true") ap.add_argument("--mux-only", action="store_true") ap.add_argument("--audio-only", action="store_true") ap.add_argument("--jobs", type=int, default=min(12, os.cpu_count() or 4)) a = ap.parse_args() wav = AUD/"final.wav" if not wav.exists() or not (AUD/"env.npz").exists() or a.force: print(f"[1/3] song… {N_BARS} bars @ {BPM:.0f}bpm = {DUR:.1f}s") wav, mix = build_song(); analyze(mix) if a.audio_only: print(f"audio -> {wav}"); return shots = build_shots() if a.sheet: contact_sheet(shots, a.sheet_at); return if not a.mux_only: sel = set(int(x) for x in a.shots.split(",") if x.strip() != "") jobs = [(s, a.force) for s in shots if not sel or s.idx in sel] print(f"[2/3] frames… {len(jobs)} shots / {N_FRAMES} frames " f"on {a.jobs} workers") import multiprocessing as mp with mp.get_context("fork").Pool(a.jobs) as pool: for r in pool.imap_unordered(render_shot, jobs): print(" ", r, flush=True) print("[3/3] mux…") out = OUT/f"{NAME}.mp4" subprocess.run( ["ffmpeg", "-y", "-framerate", str(FPS), "-i", str(FRAMES/"f%05d.png"), "-i", str(wav), "-c:v", "libx264", "-preset", "medium", "-crf", "19", "-pix_fmt", "yuv420p", "-c:a", "aac", "-b:a", "256k", "-shortest", "-movflags", "+faststart", # mp4 silently drops non-standard tag keys, so provenance rides in # `comment` / `description` / `artist`, which survive the muxer "-metadata", f"title={SUBT} — {TITLE}", "-metadata", f"artist=renders/{SETDIR}/{NAME}/render.py", "-metadata", f"comment=generator=renders/{SETDIR}/{NAME}/render.py; " f"{MUSIC_DESC}; {ENGINE_DESC}", "-metadata", f"description={MUSIC_DESC}; {ENGINE_DESC}", "-metadata", f"date={datetime.date.today().isoformat()}", str(out)], check=True, capture_output=True) try: sha = subprocess.check_output(["git", "rev-parse", "--short", "HEAD"], cwd=ROOT).decode().strip() br = subprocess.check_output(["git", "rev-parse", "--abbrev-ref", "HEAD"], cwd=ROOT).decode().strip() except Exception: sha = br = "unknown" (OUT/"PROVENANCE.txt").write_text( f"generator: renders/{SETDIR}/{NAME}/render.py\n" f"git: {sha} branch: {br}\n" f"timestamp: {datetime.datetime.now().astimezone().isoformat()}\n" f"duration: {DUR:.2f}s fps: {FPS} size: {W}x{H} (16:9)\n" f"scale: RS={RS} — native re-rasterisation from {AW}x{AH} authoring units\n" f"music: {MUSIC_DESC}\n" f"sections: {' '.join(n for n, _, _ in SECTIONS)}\n" f"engines: {ENGINE_DESC} (shot-parallel, stateless per frame)\n" f"voice: macOS say — Daniel (radiologist dictation), verified via say -v '?'\n") print(f"DONE {out} ({DUR:.1f}s)") if __name__ == "__main__": main()