"""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")