DarkRoom/tools/film-profiles/proto.py

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