DarkRoom/tools/film-profiles/convert.py

#!/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"

# Every stock spektrafilm publishes. A print paper is a stock like any other,
# distinguished only by `support: paper`; a cine print film likewise, by
# `stage: printing`. The renderer treats all three the same and the *picker*
# decides what a photographer is offered, which is why they are all converted
# here rather than filtered at the source.
STOCKS = [
    "fujifilm_c200",
    "fujifilm_crystal_archive_typeii",
    "fujifilm_pro_400h",
    "fujifilm_provia_100f",
    "fujifilm_velvia_100",
    "fujifilm_xtra_400",
    "kodak_2383",
    "kodak_2393",
    "kodak_ektachrome_100",
    "kodak_ektacolor_edge",
    "kodak_ektar_100",
    "kodak_endura_premier",
    "kodak_gold_200",
    "kodak_kodachrome_64",
    "kodak_portra_160",
    "kodak_portra_400",
    "kodak_portra_800",
    "kodak_portra_800_push1",
    "kodak_portra_800_push2",
    "kodak_portra_endura",
    "kodak_supra_endura",
    "kodak_ultra_endura",
    "kodak_ultramax_400",
    "kodak_verita_200d",
    "kodak_vision3_50d",
    "kodak_vision3_200t",
    "kodak_vision3_250d",
    "kodak_vision3_500t",
    # Black and white. Upstream ships these on `dev` only, and they are the
    # only open, measured monochrome profiles in existence -- everything else
    # is a datasheet graph nobody has digitised.
    "kodak_trix",
    "kodak_doublex",
    "kodak_2302",
]

# Stocks that live on upstream's `dev` branch rather than `main`.
#
# The black-and-white ones, which is the whole reason for the exception: there
# is no monochrome stock on `main` at all, and these are the only open,
# measured B&W profiles that exist anywhere. Pinned per stock rather than
# moving everything to `dev`, so the colour stocks stay on the released branch.
DEV_BRANCH = {"kodak_trix", "kodak_doublex", "kodak_2302"}

# Which development time to ship, for a stock whose curves are measured at
# several. Index into the profile's own `development_time` array.
#
# The middle of the published range, which is the manufacturer's standard
# process: 6.5 minutes for Double-X, 5 for 2302. The others are push and pull,
# and are the obvious next control -- the data for them ships in the upstream
# file and is dropped here only because nothing can yet ask for it.
DEVELOPMENT_INDEX = 2

# ---------------------------------------------------------------------------
# Constructed black-and-white stocks.
# ---------------------------------------------------------------------------
#
# **These are not measurements, and the profiles say so.** Ilford publish
# spectral sensitivity and characteristic curves as *graphs*, and no
# granularity figure at all -- the word does not appear in their technical
# information. Nobody has digitised them, so a profile has to be constructed.
#
# What each part actually rests on:
#
#   speed        Published and exact. ISO 400/27 for HP5 Plus is a fact.
#   contrast     ISO 6:1993 defines normal development as an average gradient
#                of 0.62 over 1.30 log-E. A standard, not a guess.
#   spectral     *Borrowed* from Kodak Double-X, which is a measured
#                panchromatic negative -- peak near 430nm, the dip at 500, the
#                cutoff at 660 -- shifted by the speed difference. Conventional
#                panchromatic sensitisation is much alike across black-and-white
#                stocks, and this is far better founded than reading pixels off
#                a printed curve.
#   silver       Neutral, which is not an approximation: developed silver
#                absorbs flat across the visible band, and Double-X's measured
#                density is flat 1.0.
#   granularity  **Estimated.** Ordered by each film's known relative grain and
#                nothing more.
#
# So these render as an Ilford-*shaped* stock -- right speed, right contrast,
# right relative grain -- and not as a measurement of one. Anyone who digitises
# the real curves should replace them outright.
PARAMETRIC_BW = [
    # stock, display name, ISO, Dmax, granularity estimate
    ("ilford_pan_f_plus", "Ilford Pan F Plus", 50, 2.5, 5.0),
    ("ilford_delta_100", "Ilford Delta 100", 100, 2.4, 6.0),
    ("ilford_fp4_plus", "Ilford FP4 Plus", 125, 2.5, 7.0),
    ("ilford_delta_400", "Ilford Delta 400", 400, 2.4, 9.0),
    ("ilford_hp5_plus", "Ilford HP5 Plus", 400, 2.6, 11.0),
]

# The stock whose measured panchromatic response the constructed ones borrow.
BW_REFERENCE = "kodak_doublex"
BW_REFERENCE_ISO = 250

# ISO 6:1993's normal-development contrast.
BW_GAMMA = 0.62

# Where mid-grey sits above the speed point, in log exposure. At gamma 0.62
# this puts an 18% card near 0.68 above fog, which is what a correctly exposed
# negative reads.
BW_SPEED_POINT = -1.1

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
        branch = "dev" if stock in DEV_BRANCH else "main"
        print(f"fetching {stock} ({branch})", file=sys.stderr)
        urllib.request.urlretrieve(
            f"https://raw.githubusercontent.com/andreavolpato/spektrafilm/{branch}"
            f"/src/spektrafilm/data/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 columns(rows, want=3):
    """Normalise a table to `want` columns per row.

    Two shapes arrive. A colour stock gives three columns and passes through.
    A black-and-white stock gives **one** -- it has one emulsion -- and is
    spread across all three, which is exact rather than approximate: three
    layers with identical sensitivity and identical curves respond identically,
    which is what one layer does.

    The dye is the exception and is handled by `split_dye`.
    """
    out = []
    for row in rows:
        if len(row) >= want:
            out.append(row[:want])
        else:
            out.append([row[0]] * want)
    return out


def split_dye(rows):
    """A single emulsion's dye, divided across the three layers.

    The renderer sums the three layers' contributions -- Beer-Lambert, so
    densities add -- and a monochrome stock has one dye, not three. Replicating
    it unchanged would treat one emulsion as three stacked copies of itself and
    render everything three times too dense. A third each reconstructs the
    single layer exactly, and stays correct anywhere the three densities differ,
    which is where the baked lookup interpolates.
    """
    out = []
    for row in rows:
        if len(row) >= 3:
            out.append(row[:3])
        else:
            v = row[0]
            out.append([None if v is None else v / 3.0] * 3)
    return out


def pick_development(rows, index):
    """One development time's column, for a stock measured at several."""
    return [[row[min(index, len(row) - 1)]] if isinstance(row, list) else [row] for row in rows]


def convert(stock):
    d = json.loads((CACHE / f"{stock}.json").read_text())
    info, data = d["info"], d["data"]
    monochrome = info.get("channel_model") == "bw"

    # A monochrome stock may be measured at several development times. Take one
    # column before anything else, so everything below sees the usual shape.
    if monochrome:
        curves = data["density_curves"]
        if curves and isinstance(curves[0], list) and len(curves[0]) > 1:
            data["density_curves"] = pick_development(curves, DEVELOPMENT_INDEX)
        base = data.get("base_density")
        if base and isinstance(base[0], list):
            data["base_density"] = [
                row[min(DEVELOPMENT_INDEX, len(row) - 1)] for row in base
            ]
    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"stage: {info['stage']}        # filming | printing",
        f"monochrome: {'true' if monochrome else 'false'}",
        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"], columns(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"], split_dye(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 columns(data["density_curves"]):
        out.append(f"  - {row(triple, 5)}")
    return "\n".join(out) + "\n"


def build_parametric_bw(stock, name, iso, density_max, granularity):
    """A constructed black-and-white profile. See PARAMETRIC_BW for its basis."""
    import math

    ref = json.loads((CACHE / f"{BW_REFERENCE}.json").read_text())["data"]
    wavelengths = ref["wavelengths"]
    n = len(wavelengths)

    # The measured panchromatic response, moved by the speed difference. A
    # faster film is more sensitive at every wavelength, which is a shift of
    # the log curve rather than a change in its shape.
    speed_shift = math.log10(iso / BW_REFERENCE_ISO)
    sensitivity = []
    for row_in in ref["log_sensitivity"]:
        v = row_in[0]
        sensitivity.append(None if v is None or v != v else v + speed_shift)

    # Developed silver, which absorbs neutrally. A third each, because the
    # renderer sums the three layers -- the same reason `split_dye` exists.
    dye = [[1.0 / 3.0] * 3 for _ in range(n)]
    base = [0.1] * n

    # A straight line of slope `gamma`, softened into a toe below and a
    # shoulder above.
    log_exposure = [-3.0 + 7.0 * i / 255.0 for i in range(256)]
    toe, shoulder = 0.35, 0.45
    curves = []
    for x in log_exposure:
        d = (
            BW_GAMMA * toe * math.log10(1 + 10 ** ((x - BW_SPEED_POINT) / toe))
            - BW_GAMMA
            * shoulder
            * math.log10(
                1 + 10 ** ((x - BW_SPEED_POINT - density_max / BW_GAMMA) / shoulder)
            )
        )
        curves.append([d] * 3)

    g = num(granularity)
    out = [
        "# Generated by tools/film-profiles/convert.py.",
        "#",
        "# *** CONSTRUCTED, NOT MEASURED. ***",
        "#",
        "# Ilford publish this film's spectral sensitivity and characteristic",
        "# curve as graphs, and no granularity figure at all, so this profile is",
        "# built rather than extracted. Its speed is the published ISO rating,",
        "# and its contrast is ISO 6:1993's normal development (average gradient",
        "# 0.62). Its spectral response is borrowed from Kodak Double-X, a",
        "# *measured* panchromatic negative, shifted by the speed difference.",
        "# Its granularity is an estimate, ordered against the other films here.",
        "#",
        "# So it renders as a film of this speed and contrast, not as a",
        "# measurement of this emulsion. Replace it outright if anyone digitises",
        "# the real curves.",
        "",
        "version: 'constructed-1'",
        "stock: " + stock,
        "name: '" + name + "'",
        "kind: negative          # negative | positive",
        "support: film    # film | paper",
        "stage: filming        # filming | printing",
        "monochrome: true",
        "reference_illuminant: D55",
        "viewing_illuminant: D50",
        "",
        "rms_granularity: [" + g + ", " + g + ", " + g + "]",
        "",
        "# log10 spectral sensitivity per layer, 380-780nm at 5nm.",
        "log_sensitivity:",
    ]
    for wl, v in zip(wavelengths, sensitivity):
        cell = "-9" if v is None else num(v)
        out.append("  - [" + cell + ", " + cell + ", " + cell + "]   # %.0fnm" % wl)

    out += ["", "# Developed silver: neutral across the band.", "dye_density:"]
    for wl, triple in zip(wavelengths, dye):
        out.append("  - " + row(triple) + "   # %.0fnm" % wl)

    out += [
        "",
        "# Base plus fog. No orange mask on a black-and-white film.",
        "base_density: " + row(base),
        "",
        "# The characteristic curve, parametric.",
        "log_exposure_min: " + num(log_exposure[0]),
        "log_exposure_max: " + num(log_exposure[-1]),
        "density_curves:",
    ]
    for triple in curves:
        out.append("  - " + 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):
        """A float literal an `f32` can actually hold.

        Seven significant figures, not eight: `f32` carries about 7.2 decimal
        digits, so an eighth is noise that rounds away at compile time -- and
        clippy's `excessive_precision` says so, which under CI's `-D warnings`
        is a failed build rather than a note.

        Always a literal, too: `%g` renders an exact zero as `0`, which is an
        integer in Rust and will not compile in an `[f32; _]`.
        """
        s = f"{v:.7g}"
        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")
    ap.add_argument(
        "--tables",
        action="store_true",
        help="also regenerate src/tables.rs (needs colour-science; the tables are "
        "fixed colour science and do not change when a stock is added)",
    )
    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")

    for stock, name, iso, dmax, gran in PARAMETRIC_BW:
        dest = PROFILE_DIR / (stock + ".yaml")
        dest.write_text(build_parametric_bw(stock, name, iso, dmax, gran))
        print(f"{dest.relative_to(ROOT)}  {dest.stat().st_size / 1024:.1f} kB  (constructed)")

    lic = CACHE / "SPEKTRAFILM_LICENSE.txt"
    if lic.exists():
        (PROFILE_DIR / "LICENSE-PROFILES.txt").write_text(lic.read_text())

    if args.tables:
        tables = ROOT / "core/dr-film/src/tables.rs"
        emit_tables(tables)
        print(f"{tables.relative_to(ROOT)}  {tables.stat().st_size / 1024:.1f} kB")

    registry = ROOT / "core/dr-film/src/built_in.rs"
    emit_registry(registry)
    print(f"{registry.relative_to(ROOT)}  {registry.stat().st_size / 1024:.1f} kB")


def emit_registry(dest):
    """The compiled-in stock list.

    Generated rather than hand-maintained, because it has to agree exactly with
    what is in `profiles/` -- a stock converted but never listed is a file that
    ships and cannot be chosen, which looks like a bug in the picker.
    """
    lines = [
        "//! Generated by tools/film-profiles/convert.py. Do not edit.",
        "//!",
        "//! The stocks compiled in as a floor. A floor rather than the whole",
        "//! story: a stock is a file, and the point of the format is that",
        "//! anyone can add one without a release.",
        "",
        "/// Each stock's id and its YAML, in the order the converter ran.",
        "pub static BUILT_IN: &[(&str, &str)] = &[",
    ]
    for stock in STOCKS + [p[0] for p in PARAMETRIC_BW]:
        lines.append(f'    ("{stock}", include_str!("../profiles/{stock}.yaml")),')
    lines.append("];")
    dest.write_text("\n".join(lines) + "\n")


if __name__ == "__main__":
    main()