dtourolle/DarkRoom
1
0
Fork 0
Code Issues 39 Pull Requests Actions Packages Projects Releases 20 Wiki Activity

Simulate a film stock from its measurements, not from someone's grade

Browse Source 
FR-DEV-3f asks for look emulation and proposes HaldCLUT import to inherit
the free film-simulation ecosystem. This takes the other road for the
stocks where the measurements exist: run the physics.

A stock here is its manufacturer's own datasheet -- spectral sensitivity,
characteristic curves, dye densities. Light exposes three emulsion layers,
the layers develop to densities, the densities are dyes that absorb, and
what is left is what reaches the eye. A colour negative comes out orange
and upside down because that is what a colour negative is; it becomes a
photograph when a paper profile prints it, with the enlarger's filtration
solved rather than dialled.

What that buys over a LUT is that the parameters stay physical. Opening up
a stop moves the picture along the film's real characteristic curve,
shoulder and all, instead of scaling a number baked at one exposure. The
data cost runs the other way too: a stock is 17 kB of published
measurements where one HaldCLUT is 800 kB of one person's grade.

It looks like it needs a spectral integration per pixel. It does not, and
that is the whole design:

  - Exposure is a 3x3 matrix. The reconstructed scene spectrum is linear
    in the sRGB triple, so the integral collapses into nine numbers,
    exactly -- no approximation.
  - The characteristic curve is three 1D functions, sampled exactly.
  - Everything after that -- dye absorption, the print through the
    negative, the paper, the viewing illuminant, the adaptation -- takes
    exactly three numbers in, so it bakes into one 32^3 lookup.

Per pixel: a matrix multiply, three curve taps, one fetch. Splitting the
curve out of the 3D lookup rather than baking one LUT over exposure is
measured, not assumed: the curve carries the sharp shape and the dye
mixing is smooth, so folding them together would need three times the
resolution for the same error. At 32^3 the worst error is 0.003 in linear
sRGB, under one 8-bit code value, and a test says so.

No wgpu dependency, deliberately, and the same isolation argument dr-lens
makes: the model is plain f32 with a documented layout, so every property
worth asserting is asserted on the CPU. Binding it to a texture is dr-gpu's
job and is not done here yet.

The expected values in tests/ came from a Python prototype running against
a different colour-science stack. Agreement to three decimals is evidence
about the model rather than about one implementation of it -- a transposed
matrix or a mispasted observer row would pass every unit test and fail
that one.

Profiles are converted from spektrafilm by Andrea Volpato, CC BY-SA 4.0.
The converter is in the tree and runnable, so what was changed from
upstream is auditable rather than taken on trust; profiles/CHANGELOG.txt
records it, including the one deliberate deviation -- Mallett & Yuksel's
1 kB basis instead of Hanatos's 4 MB table, which costs accuracy at the
gamut edge and saves four megabytes.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
dtourolle
2026-08-25 14:44:54 +02:00
co-authored by Claude Opus 5
parent 940058c78a
commit b6a95e1965
19 changed files with 3931 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tools/film-profiles/convert.py
+229
View File
@@ -0,0 +1,229 @@
#!/usr/bin/env python3
"""Convert spektrafilm film profiles into DarkRoom's own compact format, and
generate the colour-science tables `dr-film` compiles in.
python3 tools/film-profiles/convert.py --fetch
Kept in the tree, and kept runnable, so that the conversion from upstream is
reproducible and auditable rather than a one-off paste. CC BY-SA 4.0 requires
that a modified copy say it was modified; this script *is* the statement of
what was done, and `core/dr-film/profiles/CHANGELOG.txt` records it in prose.
The upstream profiles are published on the same 380-780nm, 5nm grid as the
Mallett 2019 sRGB basis and the CIE 1931 observer, so nothing here resamples
anything -- the conversion is a trim and a reformat, not an interpolation.
"""
import argparse
import json
import pathlib
import sys
import urllib.request
ROOT = pathlib.Path(__file__).resolve().parents[2]
PROFILE_DIR = ROOT / "core/dr-film/profiles"
UPSTREAM = "https://raw.githubusercontent.com/andreavolpato/spektrafilm/main"
UPSTREAM_PROFILES = f"{UPSTREAM}/src/spektrafilm/data/profiles"
# The stocks we ship. Each is a colour stock whose datasheet spektrafilm has
# already digitised; a print paper is a stock like any other, distinguished
# only by `support: paper`.
STOCKS = [
"kodak_portra_400",
"kodak_kodachrome_64",
"kodak_portra_endura",
]
CACHE = pathlib.Path(__file__).parent / "upstream"
def fetch():
CACHE.mkdir(exist_ok=True)
for stock in STOCKS:
dest = CACHE / f"{stock}.json"
if dest.exists():
continue
print(f"fetching {stock}", file=sys.stderr)
urllib.request.urlretrieve(f"{UPSTREAM_PROFILES}/{stock}.json", dest)
lic = CACHE / "SPEKTRAFILM_LICENSE.txt"
if not lic.exists():
urllib.request.urlretrieve(f"{UPSTREAM}/SPEKTRAFILM_LICENSE.txt", lic)
def num(v, places=6):
"""A null becomes an explicit 0, so the YAML has no holes to interpret."""
if v is None or v != v:
return "0"
s = f"{v:.{places}g}"
return "0" if s in ("-0", "-0.0") else s
def row(values, places=6):
return "[" + ", ".join(num(v, places) for v in values) + "]"
def convert(stock):
d = json.loads((CACHE / f"{stock}.json").read_text())
info, data = d["info"], d["data"]
n = len(data["wavelengths"])
assert data["wavelengths"][0] == 380.0 and data["wavelengths"][-1] == 780.0 and n == 81, (
f"{stock}: unexpected wavelength grid; the conversion assumes 380-780nm at 5nm"
)
out = [
"# Generated by tools/film-profiles/convert.py from spektrafilm.",
"# Do not edit by hand: re-run the converter instead.",
"#",
"# spektrafilm by Andrea Volpato, https://github.com/andreavolpato/spektrafilm",
"# Licensed CC BY-SA 4.0. Modified for DarkRoom: trimmed to the fields the",
"# renderer uses and reformatted; see profiles/CHANGELOG.txt.",
"",
f"version: {d['metadata']['version']!r}",
f"stock: {info['stock']}",
f"name: {info['name']!r}",
f"kind: {info['type']} # negative | positive",
f"support: {info['support']} # film | paper",
f"reference_illuminant: {info['reference_illuminant']}",
f"viewing_illuminant: {info['viewing_illuminant']}",
]
if info.get("target_print"):
out.append(f"target_print: {info['target_print']}")
out += [
"",
"# log10 spectral sensitivity per layer, 380-780nm at 5nm, in R,G,B layer",
"# order. A null upstream means the datasheet has no reading there, which is",
"# blindness, so it is written as the sentinel the loader reads as such.",
"log_sensitivity:",
]
for wl, triple in zip(data["wavelengths"], data["log_sensitivity"]):
vals = ["-9" if (v is None or v != v) else num(v) for v in triple]
out.append(f" - [{', '.join(vals)}] # {wl:.0f}nm")
out += [
"",
"# Spectral density of each layer's dye at unit density, same grid and order.",
"dye_density:",
]
for wl, triple in zip(data["wavelengths"], data["channel_density"]):
out.append(f" - {row(triple)} # {wl:.0f}nm")
base = data.get("base_density")
out += [
"",
"# The support's own density -- film base plus, for a colour negative, the",
"# orange mask. Flat zero where the datasheet does not give it.",
f"base_density: {row(base) if base else row([0.0] * n)}",
"",
"# The characteristic curves: density against log10 exposure, sampled",
f"# uniformly over [{data['log_exposure'][0]:g}, {data['log_exposure'][-1]:g}].",
f"log_exposure_min: {num(data['log_exposure'][0])}",
f"log_exposure_max: {num(data['log_exposure'][-1])}",
"density_curves:",
]
for triple in data["density_curves"]:
out.append(f" - {row(triple, 5)}")
return "\n".join(out) + "\n"
TABLE_HEADER = '''//! Generated by tools/film-profiles/convert.py. Do not edit.
//!
//! The fixed colour science: the observer, the illuminants and the spectral
//! basis. None of it is per-stock, all of it is published data, and together it
//! is under 6 kB of source -- which is the point. A film simulation's data cost
//! is dominated by whatever it uses to turn a pixel back into a spectrum, and a
//! basis is three curves where a coefficient table is megabytes.
/// The lowest wavelength sampled, in nanometres.
pub const LAMBDA_MIN: f32 = 380.0;
/// The spacing between samples, in nanometres.
pub const LAMBDA_STEP: f32 = 5.0;
/// How many wavelengths every spectral table carries.
///
/// The profiles, the observer and the basis all arrive on this grid already, so
/// nothing in this crate resamples anything.
pub const SPECTRUM: usize = 81;
'''
def emit_tables(dest):
import numpy as np
import colour
grid = colour.SpectralShape(380, 780, 5)
cmf = colour.MSDS_CMFS["CIE 1931 2 Degree Standard Observer"].copy().align(grid).values
basis = colour.recovery.MSDS_BASIS_FUNCTIONS_sRGB_MALLETT2019.copy().align(grid).values
def lit(v):
"""Always a float literal: `%g` renders an exact zero as `0`, which is
an integer in Rust and will not compile in an `[f32; _]`."""
s = f"{v:.8g}"
return s if any(ch in s for ch in ".eE") else s + ".0"
def table(name, doc, values):
lines = [f"\n{doc}\npub static {name}: [[f32; 3]; SPECTRUM] = ["]
for wl, triple in zip(grid.wavelengths, values):
cells = ", ".join(lit(v) for v in triple)
lines.append(f" [{cells}], // {wl:.0f}nm")
lines.append("];")
return "\n".join(lines)
def flat(name, doc, values, wavelengths):
lines = [f"\n{doc}\npub static {name}: [f32; SPECTRUM] = ["]
for i in range(0, len(values), 6):
chunk = ", ".join(lit(v) for v in values[i:i + 6])
lines.append(f" {chunk},")
lines.append("];")
return "\n".join(lines)
parts = [TABLE_HEADER]
parts.append(table(
"OBSERVER",
"/// CIE 1931 2-degree standard observer, x-bar/y-bar/z-bar per wavelength.",
cmf,
))
parts.append(table(
"SRGB_BASIS",
"/// Mallett & Yuksel (2019) sRGB reflectance basis: the three smooth,\n"
"/// non-negative spectra that reconstruct any sRGB colour exactly.\n"
"///\n"
"/// This is what makes the exposure step a 3x3 matrix rather than a\n"
"/// per-pixel spectral integration -- see [`crate::bake`].",
basis,
))
for name in ("D50", "D55", "D65"):
sd = colour.SDS_ILLUMINANTS[name].copy().align(grid).values
parts.append(flat(
f"ILLUMINANT_{name}",
f"/// CIE standard illuminant {name}, normalised to unit mean.",
sd / sd.mean(),
grid.wavelengths,
))
dest.write_text("\n".join(parts) + "\n")
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--fetch", action="store_true", help="download upstream profiles first")
args = ap.parse_args()
if args.fetch:
fetch()
if not CACHE.exists():
sys.exit("no upstream cache; run with --fetch")
PROFILE_DIR.mkdir(parents=True, exist_ok=True)
for stock in STOCKS:
dest = PROFILE_DIR / f"{stock}.yaml"
dest.write_text(convert(stock))
print(f"{dest.relative_to(ROOT)} {dest.stat().st_size / 1024:.1f} kB")
lic = CACHE / "SPEKTRAFILM_LICENSE.txt"
if lic.exists():
(PROFILE_DIR / "LICENSE-PROFILES.txt").write_text(lic.read_text())
tables = ROOT / "core/dr-film/src/tables.rs"
emit_tables(tables)
print(f"{tables.relative_to(ROOT)} {tables.stat().st_size / 1024:.1f} kB")
if __name__ == "__main__":
main()
tools/film-profiles/proto.py
+244
View File
@@ -0,0 +1,244 @@
"""Prototype of the spectral film chain, to be ported to core/dr-film.
linear sRGB -> reflectance (Mallett 2019 basis) -> x reference illuminant ->
layer exposures -> densities -> dye transmittance -> XYZ -> linear sRGB.
Everything runs on the profiles' own 380-780nm @5nm, 81-sample grid, which is
also the grid the Mallett basis is published on, so nothing is resampled.
Conventions follow spektrafilm's own reference implementation rather than being
invented here, because the profile data is calibrated against them:
- exposure is normalised by the *green* layer's mid-grey response, one shared
scalar for all three layers. The residual channel imbalance is the film's
real one, and for a negative it is the print stage's job to balance it out.
- the viewing step chromatically adapts from the viewing illuminant to the
output space's white, rather than dividing XYZ channelwise.
"""
import json
import numpy as np
import colour
GRID = colour.SpectralShape(380, 780, 5)
WL = GRID.wavelengths
N = len(WL)
CMF = colour.MSDS_CMFS["CIE 1931 2 Degree Standard Observer"].copy().align(GRID).values # (81,3)
BASIS = colour.recovery.MSDS_BASIS_FUNCTIONS_sRGB_MALLETT2019.copy().align(GRID).values # (81,3)
MID_GREY = 0.184
def illuminant(name):
sd = colour.SDS_ILLUMINANTS[name].copy().align(GRID).values
return sd / sd.mean()
def load(path):
d = json.load(open(path))
data = d["data"]
p = {
"info": d["info"],
"log_sensitivity": np.array(data["log_sensitivity"], dtype=float), # (81,3)
"dye_density": np.array(data["channel_density"], dtype=float), # (81,3)
"log_exposure": np.array(data["log_exposure"], dtype=float), # (256,)
"density_curves": np.array(data["density_curves"], dtype=float), # (256,3)
}
base = data.get("base_density")
p["base_density"] = (
np.array([0.0 if v is None else v for v in base], dtype=float)
if base is not None else np.zeros(N)
)
# A sample is null where the datasheet has no data. Sensitivity there means
# "blind at this wavelength"; density there means "no absorption".
p["log_sensitivity"] = np.nan_to_num(p["log_sensitivity"], nan=-9.0)
p["dye_density"] = np.nan_to_num(p["dye_density"], nan=0.0)
return p
def sensitivity(profile):
return 10.0 ** profile["log_sensitivity"] # (81, 3 layers)
def exposure_matrix(profile):
"""3x3: linear sRGB -> the three layers' exposures, mid-grey normalised.
Exposure is an integral of sensitivity against the scene spectrum, and the
scene spectrum is linear in the sRGB coefficients, so the whole step is a
matrix. This is what makes a per-pixel spectral integration unnecessary on
the way *in* -- the only place the spectrum is genuinely needed is the dye
transmittance on the way out, which is a function of three densities and so
bakes into a small 3D LUT.
"""
ill = illuminant(profile["info"].get("reference_illuminant", "D55"))
sens = sensitivity(profile)
m = sens.T @ (BASIS * ill[:, None]) # (3 layers, 3 sRGB)
# Normalised on the green layer's mid-grey response, matching spektrafilm.
# A per-layer normalisation would silently absorb the film's own channel
# balance, which is a large part of what distinguishes one stock's look
# from another's.
mid_grey_raw = (ill * MID_GREY) @ sens # (3 layers,)
return m / mid_grey_raw[1]
def densities(profile, log_exposure):
"""Sample the tabulated characteristic curves. (...,3) -> (...,3)."""
out = np.empty_like(log_exposure)
for c in range(3):
out[..., c] = np.interp(
log_exposure[..., c], profile["log_exposure"], profile["density_curves"][:, c]
)
return out
def transmittance(profile, dens):
"""Dye densities -> spectral transmittance. (...,3) -> (...,81)."""
return 10.0 ** (-(dens @ profile["dye_density"].T + profile["base_density"]))
def view(profile, trans, view_illuminant=None, output_space="sRGB"):
"""Spectral transmittance -> linear output RGB, viewed on a light table.
Adapted from the viewing illuminant to the output space's own white, so a
clear frame comes out neutral instead of carrying the light table's colour.
"""
ill = illuminant(view_illuminant or profile["info"].get("viewing_illuminant", "D50"))
norm = (ill * CMF[:, 1]).sum()
white_xyz = (ill[:, None] * CMF).sum(axis=0) / norm
xyz = ((trans * ill)[..., None] * CMF).sum(axis=-2) / norm
return colour.XYZ_to_RGB(
xyz,
colourspace=output_space,
illuminant=colour.XYZ_to_xy(white_xyz),
apply_cctf_encoding=False,
)
def develop(profile, rgb, exposure_ev=0.0):
"""Full chain for a positive (reversal) stock, viewed directly."""
raw = np.asarray(rgb, dtype=float) @ exposure_matrix(profile).T * (2.0 ** exposure_ev)
log_raw = np.log10(np.maximum(raw, 0.0) + 1e-10)
return view(profile, transmittance(profile, densities(profile, log_raw)))
def report(name):
p = load(name)
print(f"=== {p['info']['name']} ({p['info']['type']}) ===")
print("exposure matrix (rows = layers, cols = R,G,B):")
print(np.array2string(exposure_matrix(p), precision=4, suppress_small=True))
clear = view(p, np.ones(N))
print(f"clear frame -> [{clear[0]:.4f} {clear[1]:.4f} {clear[2]:.4f}]"
f" (should be neutral)")
ramp = np.array([[v] * 3 for v in (0.02, 0.09, 0.184, 0.4, 0.8)])
out = develop(p, ramp)
print("neutral ramp in -> linear sRGB out:")
for v, o in zip(ramp[:, 0], out):
print(f" {v:6.3f} -> [{o[0]:8.4f} {o[1]:8.4f} {o[2]:8.4f}]"
f" spread {o.max() - o.min():+.4f}")
print("primaries in -> out:")
for i, lbl in enumerate("RGB"):
rgb = np.full(3, 0.05)
rgb[i] = 0.5
o = develop(p, rgb)
print(f" {lbl}: [{o[0]:7.4f} {o[1]:7.4f} {o[2]:7.4f}]")
print()
# ---------------------------------------------------------------------------
# The print stage, for negative stocks.
# ---------------------------------------------------------------------------
def blackbody(temperature_k):
"""Planck's law on the grid, normalised. Analytic, so it costs no data.
Stands in for the enlarger's tungsten-halogen lamp. The lamp's own colour
is very nearly irrelevant to the result because the filtration below is
*solved* rather than specified: whatever cast the source has, the balance
step removes it, exactly as a darkroom worker dials it out on the head.
"""
wl = WL * 1e-9
h, c, kb = 6.62607015e-34, 2.99792458e8, 1.380649e-23
radiance = (2 * h * c**2) / (wl**5 * (np.exp(h * c / (wl * kb * temperature_k)) - 1))
return radiance / radiance.mean()
ENLARGER = blackbody(3400.0)
def paper_exposure_operator(neg, paper):
"""Return f(neg_densities) -> the paper's three layer exposures.
Not a matrix: the negative's transmittance is 10^-D, so the paper's
exposure is exponential in the negative's densities. This is the one
genuinely spectral step left, and it takes exactly three numbers in -- which
is what lets the whole negative-plus-print chain bake into one 3D LUT.
"""
sens = sensitivity(paper) # (81, 3)
def f(dens):
trans = transmittance(neg, dens) # (...,81)
return (trans * ENLARGER) @ sens # (...,3)
return f
def print_balance(neg, paper, print_exposure_ev=0.0):
"""Per-layer log offsets that make a mid-grey scene print neutral mid-grey.
This is the enlarger's colour head, solved instead of dialled. It is also
where a colour negative's orange mask goes: the mask is a fixed density, so
balancing mid-grey to neutral removes it, which is why a printed negative
looks like a photograph and a scanned one looks orange.
"""
expose = paper_exposure_operator(neg, paper)
mid_dens = densities(neg, np.log10((np.full(3, MID_GREY) @ exposure_matrix(neg).T) + 1e-10))
mid_raw = expose(mid_dens) # (3,)
# Where on the paper's curve mid-grey should land: the point whose density
# reflects 18%, read off the average of the three curves.
target_density = -np.log10(MID_GREY) - np.nanmean(paper["base_density"])
curve_av = np.nanmean(paper["density_curves"], axis=1)
target_log_e = np.interp(target_density, curve_av, paper["log_exposure"])
return target_log_e - np.log10(mid_raw) + print_exposure_ev * np.log10(2.0)
def print_develop(neg, paper, rgb, exposure_ev=0.0, print_exposure_ev=0.0):
"""Full negative -> print chain, viewed as a reflection print."""
raw = np.asarray(rgb, dtype=float) @ exposure_matrix(neg).T * (2.0 ** exposure_ev)
neg_dens = densities(neg, np.log10(np.maximum(raw, 0.0) + 1e-10))
offsets = print_balance(neg, paper, print_exposure_ev)
paper_raw = paper_exposure_operator(neg, paper)(neg_dens)
paper_dens = densities(paper, np.log10(paper_raw + 1e-10) + offsets)
return view(paper, transmittance(paper, paper_dens))
def report_print(neg_name, paper_name):
neg, paper = load(neg_name), load(paper_name)
print(f"=== {neg['info']['name']} printed on {paper['info']['name']} ===")
print("print balance (per-layer logE offset):",
np.array2string(print_balance(neg, paper), precision=4))
ramp = np.array([[v] * 3 for v in (0.02, 0.09, 0.184, 0.4, 0.8)])
out = print_develop(neg, paper, ramp)
print("neutral ramp in -> linear sRGB out:")
for v, o in zip(ramp[:, 0], out):
print(f" {v:6.3f} -> [{o[0]:8.4f} {o[1]:8.4f} {o[2]:8.4f}]"
f" spread {o.max() - o.min():+.4f}")
print("primaries in -> out:")
for i, lbl in enumerate("RGB"):
rgb = np.full(3, 0.05)
rgb[i] = 0.5
o = print_develop(neg, paper, rgb)
print(f" {lbl}: [{o[0]:7.4f} {o[1]:7.4f} {o[2]:7.4f}]")
print()
if __name__ == "__main__":
for name in ("kodachrome64.json", "portra400.json"):
report(name)
report_print("portra400.json", "portra_endura.json")