//! Turning a stock into something a shader can run.
//!
//! # The decomposition
//!
//! A physically-honest film simulation looks like it needs a spectral
//! integration per pixel, and vkdt's does exactly that. It does not have to,
//! and the reason is worth writing down because it is what makes this cheap
//! enough to run on a phone:
//!
//! 1. **Exposure is a 3×3 matrix.** A layer's exposure is
//!    `∫ S(λ)·L(λ) dλ`, and the scene spectrum `L` reconstructed from an sRGB
//!    triple is *linear* in that triple — that is what a spectral basis is. So
//!    the whole integral collapses into nine numbers, computed once, exactly.
//!    No approximation is involved.
//!
//! 2. **The characteristic curve is three 1D functions.** Sampled exactly, at
//!    [`crate::profile::CURVE_SAMPLES`].
//!
//! 3. **Everything after that is a function of three densities.** The dye
//!    transmittance, the print exposure through the negative, the paper's own
//!    curves and dyes, the viewing illuminant, the adaptation — all of it takes
//!    three numbers in and gives three numbers out. So it bakes into one small
//!    3D lookup, and the per-pixel cost is a matrix multiply, three curve taps
//!    and one texture fetch. A print is two: the film's lookup ends at the
//!    paper's log exposure, where the enlarger's exposure is an addition, and
//!    the paper's curve and lookup take it from there — see [`Paper`].
//!
//! Splitting 2 from 3 rather than baking a single LUT over exposure is
//! deliberate and measured: the curve carries all of the sharp shape and the
//! dye mixing is smooth, so putting the curve in the 3D LUT would force it
//! three times larger for the same error.

use crate::profile::{Profile, CURVE_SAMPLES};
use crate::spectrum::{illuminant, Spectrum, Viewing};
use crate::tables::{SPECTRUM, SRGB_BASIS};

/// The mid-grey a photographic exposure is reckoned from.
///
/// 18.4% rather than 18%: it is the value the upstream profiles are calibrated
/// against, and a profile calibrated at one grey and rendered at another is
/// off by a fraction of a stop everywhere.
pub const MID_GREY: f32 = 0.184;

/// The edge length of the baked density lookup.
///
/// 32 holds the worst-case interpolation error to about 0.003 in linear sRGB,
/// which is below one 8-bit code value, in 384 kB. Doubling it buys a factor
/// of four in error for eight times the memory, and there is nothing to spend
/// that on: the error is already under what the output can represent.
pub const LUT_SIZE: usize = 32;

/// TRACES: FR-DEV-3f
/// The most development times a stock may measure: one curve row, and one
/// push station, each. Every stock shipped measures five; the ceiling is what
/// the shader's fixed uniform block can hold.
pub const MAX_CURVE_ROWS: usize = 8;

/// What to develop: the materials, and where the enlarger is balanced.
///
/// **Not how far, and not how bright.** Push, print exposure and camera
/// exposure are [`Settings`], evaluated per pixel against these tables, so
/// that a mask layer can hold its own and a pixel under it can take the
/// weighted average of everyone's (FR-DEV-3f). What is left here is what a
/// photograph has one of.
pub struct Recipe<'a> {
    /// The stock the picture was taken on.
    pub film: &'a Profile,
    /// The paper it is printed on. `None` views the film directly, which is
    /// what a reversal stock wants and what makes a negative come out orange
    /// and inverted — that being what a negative actually looks like.
    pub print: Option<&'a Profile>,
    /// The camera exposure the enlarger is balanced at, in stops. Ignored
    /// without a `print`.
    ///
    /// The *photograph's* exposure, never a region's. An enlarger has one
    /// filtration for the whole print: a negative exposed a stop brighter in
    /// one corner prints a stop darker there, and that difference is the
    /// picture — balancing it away per pixel would erase every local exposure
    /// change a layer made.
    pub exposure_ev: f32,
}

impl<'a> Recipe<'a> {
    /// The straightforward reading of a stock: reversal viewed directly,
    /// negative printed on the paper its datasheet names.
    pub fn new(film: &'a Profile, print: Option<&'a Profile>) -> Self {
        Self {
            film,
            print,
            exposure_ev: 0.0,
        }
    }
}

/// TRACES: FR-DEV-3f
/// What a pixel is developed with, against a [`Baked`] stock.
///
/// The shader's uniforms, as the CPU sees them: every field is linear in what
/// the tables are indexed by, which is what lets the composer blend several
/// layers' settings into one before the fragment runs.
#[derive(Debug, Clone, Copy, Default, PartialEq)]
pub struct Settings {
    /// Camera exposure, in stops: a gain on the scene.
    pub exposure_ev: f32,
    /// Development, in stops of push. Positive develops longer. Nothing for a
    /// stock measured at one process, of which there are many.
    pub push_stops: f32,
    /// Enlarger exposure, in stops. Nothing without a print.
    pub print_exposure_ev: f32,
}

/// A recipe reduced to three tables.
///
/// Plain `f32` with a documented layout, and no notion of a texture: what to
/// bind this to is dr-gpu's decision, and keeping it out of here is what lets
/// the whole model be tested on the CPU.
#[derive(Debug, Clone)]
pub struct Baked {
    /// Linear sRGB to the three layers' exposure, before the log — row `l`,
    /// column `c` is layer `l`'s response to sRGB channel `c`. At unit gain:
    /// [`Settings::exposure_ev`] is applied per pixel.
    pub exposure_matrix: [[f32; 3]; 3],
    /// The characteristic curves: `curve_rows` rows of `CURVE_SAMPLES`
    /// samples, row after row, each uniform over
    /// `[curve_log_min, curve_log_max]`.
    ///
    /// Row `r` is the stock as measured at its `r`th development time, which
    /// is push [`Self::push_stations`]`[r]`. The rows are the measurements
    /// themselves rather than a resampling: between two, density is linear in
    /// push (development is interpolated in log time, and push is log time),
    /// so interpolating the rows by push reproduces
    /// [`Profile::curves_at_push`] exactly. A stock measured at one process
    /// has one row.
    pub curves: Vec<[f32; 3]>,
    pub curve_rows: usize,
    /// The push each row was developed to, ascending, one per row.
    pub push_stations: Vec<f32>,
    pub curve_log_min: f32,
    pub curve_log_max: f32,
    /// Film density to what comes next, `LUT_SIZE³` entries uniform over
    /// `[0, density_max]` on each axis: linear sRGB when the film is viewed
    /// directly, and the paper's log₁₀ exposure through it, per layer, when it
    /// is printed.
    ///
    /// **The red axis varies fastest**, then green, then blue — that is,
    /// `lut[(b * size + g) * size + r]`. Stated because it is not the order
    /// this loop reads most naturally, and it is not arbitrary: it is the
    /// order a 3D texture upload expects, so the consumer can hand the slice
    /// straight to the driver. Filling it the other way round renders a
    /// picture with red and blue transposed, which looks like a plausible
    /// photograph of the wrong colour.
    pub lut: Vec<[f32; 3]>,
    /// The paper, when there is one. See [`Paper`].
    pub paper: Option<Paper>,
    /// The deepest density any row develops to, so one lookup covers every
    /// push.
    pub density_max: f32,
    pub lut_size: usize,
}

/// TRACES: FR-DEV-3f
/// The print half of a baked stock: enlarger to paper to viewing.
///
/// Split from the film's lookup at the paper's log exposure, for the reason
/// the film is split from its own curve. The enlarger's exposure is a shift
/// *in that log exposure*, the same stops on all three layers, so a print
/// exposure is an addition between the two lookups — exact at any value and
/// free per pixel. Baking it into one lookup instead needs a slice per
/// setting, and interpolating between slices misses by several code values,
/// because the paper's curve is the sharpest thing in the print.
#[derive(Debug, Clone)]
pub struct Paper {
    /// The enlarger's filtration, per layer, in log₁₀ exposure: what makes a
    /// mid-grey scene print neutral at the photograph's exposure. See
    /// [`Recipe::exposure_ev`].
    pub balance: [f32; 3],
    /// The paper's characteristic curves, `CURVE_SAMPLES` samples uniform
    /// over `[log_min, log_max]`.
    pub curves: Vec<[f32; 3]>,
    pub log_min: f32,
    pub log_max: f32,
    /// Paper density to linear sRGB, laid out as [`Baked::lut`] is, uniform
    /// over `[0, density_max]`.
    pub lut: Vec<[f32; 3]>,
    pub density_max: f32,
}

/// Where `push` falls among the rows: the lower row and the fraction toward
/// the next. Clamped at both ends, as `curves_at_push` clamps to the first and
/// last measured process.
fn push_row(stations: &[f32], push: f32) -> (usize, f32) {
    if stations.len() < 2 {
        return (0, 0.0);
    }
    let last = stations.len() - 1;
    let hi = stations
        .iter()
        .position(|p| *p >= push)
        .unwrap_or(last)
        .max(1);
    let lo = hi - 1;
    let f = (push - stations[lo]) / (stations[hi] - stations[lo]).max(1e-6);
    (lo, f.clamp(0.0, 1.0))
}

impl Baked {
    /// Look a colour up the way the shader will, at the stock's own settings.
    pub fn apply(&self, rgb: [f32; 3]) -> [f32; 3] {
        self.apply_at(rgb, &Settings::default())
    }

    /// Look a colour up the way the shader will, for tests and for previews.
    pub fn apply_at(&self, rgb: [f32; 3], settings: &Settings) -> [f32; 3] {
        let gain = 2f32.powf(settings.exposure_ev);
        let mut log_exposure = [0.0f32; 3];
        for (l, slot) in log_exposure.iter_mut().enumerate() {
            let m = self.exposure_matrix[l];
            let e = gain * (m[0] * rgb[0] + m[1] * rgb[1] + m[2] * rgb[2]);
            *slot = (e.max(0.0) + 1e-10).log10();
        }
        let density = self.sample_curves(log_exposure, settings.push_stops);
        let through = sample_cube(&self.lut, self.lut_size, density, self.density_max);
        let Some(paper) = &self.paper else {
            return through;
        };
        let shift = settings.print_exposure_ev * 2f32.log10();
        let paper_log = [0, 1, 2].map(|l| through[l] + paper.balance[l] + shift);
        let paper_density = sample_curve(&paper.curves, paper.log_min, paper.log_max, paper_log);
        sample_cube(&paper.lut, self.lut_size, paper_density, paper.density_max)
    }

    fn sample_curves(&self, log_exposure: [f32; 3], push_stops: f32) -> [f32; 3] {
        let (row, g) = push_row(&self.push_stations, push_stops);
        let lo = self.sample_curve_row(log_exposure, row);
        if self.curve_rows < 2 {
            return lo;
        }
        let hi = self.sample_curve_row(log_exposure, row + 1);
        [0, 1, 2].map(|c| lo[c] * (1.0 - g) + hi[c] * g)
    }

    fn sample_curve_row(&self, log_exposure: [f32; 3], row: usize) -> [f32; 3] {
        let samples = self.curves.len() / self.curve_rows;
        let curve = &self.curves[row * samples..(row + 1) * samples];
        sample_curve(curve, self.curve_log_min, self.curve_log_max, log_exposure)
    }
}

/// Three curves sampled uniformly over `[log_min, log_max]`, read at a log
/// exposure per layer. Clamped at both ends, as the shader's is.
fn sample_curve(curve: &[[f32; 3]], log_min: f32, log_max: f32, at: [f32; 3]) -> [f32; 3] {
    let last = curve.len() - 1;
    let span = log_max - log_min;
    let mut out = [0.0f32; 3];
    for (c, slot) in out.iter_mut().enumerate() {
        let t = ((at[c] - log_min) / span).clamp(0.0, 1.0) * last as f32;
        let i = (t.floor() as usize).min(last - 1);
        let f = t - i as f32;
        *slot = curve[i][c] * (1.0 - f) + curve[i + 1][c] * f;
    }
    out
}

/// A cube of `n³` triples over `[0, max]` per axis, red fastest, read
/// trilinearly.
fn sample_cube(lut: &[[f32; 3]], n: usize, density: [f32; 3], max: f32) -> [f32; 3] {
    let mut base = [0usize; 3];
    let mut frac = [0f32; 3];
    for c in 0..3 {
        let t = (density[c] / max).clamp(0.0, 1.0) * (n - 1) as f32;
        base[c] = (t.floor() as usize).min(n - 2);
        frac[c] = t - base[c] as f32;
    }
    let mut out = [0.0f32; 3];
    for dx in 0..2 {
        for dy in 0..2 {
            for dz in 0..2 {
                let w = if dx == 0 { 1.0 - frac[0] } else { frac[0] }
                    * if dy == 0 { 1.0 - frac[1] } else { frac[1] }
                    * if dz == 0 { 1.0 - frac[2] } else { frac[2] };
                let e = lut[((base[2] + dz) * n + base[1] + dy) * n + base[0] + dx];
                for c in 0..3 {
                    out[c] += w * e[c];
                }
            }
        }
    }
    out
}

/// Linear sRGB to the three layers' exposure, mid-grey normalised.
///
/// Normalised on the *green* layer alone, one shared scalar for all three.
/// Doing it per layer is tempting and wrong: it would silently flatten the
/// film's own channel balance, which is a large part of what distinguishes one
/// stock from another. Where the balance genuinely has to come out — printing a
/// negative — it is [`print_balance`]'s job, which is also where it belongs
/// physically.
pub fn exposure_matrix(film: &Profile) -> [[f32; 3]; 3] {
    let reference = illuminant(&film.reference_illuminant);
    let sensitivity = film.sensitivity();

    let mut m = [[0.0f32; 3]; 3];
    let mut mid_grey = [0.0f32; 3];
    for i in 0..SPECTRUM {
        for layer in 0..3 {
            let s = sensitivity[i][layer] * reference[i];
            mid_grey[layer] += s * MID_GREY;
            for channel in 0..3 {
                m[layer][channel] += s * SRGB_BASIS[i][channel];
            }
        }
    }

    let scale = 1.0 / mid_grey[1];
    for row in &mut m {
        for v in row.iter_mut() {
            *v *= scale;
        }
    }
    m
}

/// The enlarger head's filtration, solved rather than dialled.
///
/// Returns the per-layer log exposure offsets that make a mid-grey scene print
/// as a neutral mid-grey. This is also where a colour negative's orange mask
/// goes: the mask is a fixed density, so balancing mid-grey to neutral cancels
/// it — which is why a printed negative looks like a photograph while a scanned
/// one looks orange.
pub fn print_balance(film: &Profile, paper: &Profile, exposure_ev: f32) -> [f32; 3] {
    let matrix = exposure_matrix(film);
    let scene = MID_GREY * 2f32.powf(exposure_ev);
    let mut log_exposure = [0.0f32; 3];
    for (l, slot) in log_exposure.iter_mut().enumerate() {
        let m = matrix[l];
        *slot = ((m[0] + m[1] + m[2]) * scene + 1e-10).log10();
    }
    let mid_raw = paper_exposure(film, paper, film.density_at(log_exposure));

    // Where on the paper's curve mid-grey belongs: the density that reflects
    // 18%, read off the average of the three curves. Averaged because the
    // point of the balance is that the three end up at the same place.
    let target_density = -MID_GREY.log10() - mean(&paper.base_density);
    let target = invert_mean_curve(paper, target_density);

    let mut offsets = [0.0f32; 3];
    for (l, slot) in offsets.iter_mut().enumerate() {
        *slot = target - (mid_raw[l] + 1e-10).log10();
    }
    offsets
}

/// The paper's three layer exposures, printing through a negative at these
/// densities.
///
/// The one genuinely spectral step left in the chain — the negative's
/// transmittance is `10^-D`, so this is not a matrix and cannot be made into
/// one. It takes exactly three numbers in, which is what lets the whole
/// negative-and-print chain still bake into a 3D lookup.
fn paper_exposure(film: &Profile, paper: &Profile, density: [f32; 3]) -> [f32; 3] {
    let enlarger = illuminant(&paper.reference_illuminant);
    let sensitivity = paper.sensitivity();
    let transmittance = film.transmittance(density);

    let mut raw = [0.0f32; 3];
    for i in 0..SPECTRUM {
        let light = transmittance[i] * enlarger[i];
        for layer in 0..3 {
            raw[layer] += light * sensitivity[i][layer];
        }
    }
    raw
}

/// Bake a recipe into the tables a shader runs.
pub fn bake(recipe: &Recipe) -> Baked {
    let film = recipe.film;
    // At unit gain. Camera exposure is a scalar on a linear quantity, so the
    // shader applies it for the price of one multiply — and has to, since a
    // layer may hold its own.
    let matrix = exposure_matrix(film);

    // TRACES: FR-DEV-3f
    // Every measured process, not the one the slider is at: the shader
    // interpolates between rows per pixel, so a layer can push a region.
    // Resampled to one length because the rows share a texture.
    let measured = film.development_curves.len() >= 2
        && film.development_times.len() == film.development_curves.len();
    let (curves, push_stations): (Vec<[f32; 3]>, Vec<f32>) = if measured {
        let rows = film.development_curves.len().min(MAX_CURVE_ROWS);
        (
            film.development_curves[..rows]
                .iter()
                .flat_map(|c| resample(c))
                .collect(),
            film.development_times[..rows]
                .iter()
                .map(|t| 2.0 * (t / film.development_normal).log2())
                .collect(),
        )
    } else {
        (resample(&film.density_curves), vec![0.0])
    };
    let curve_rows = push_stations.len();
    // The ceiling of the deepest row, so one lookup covers every push.
    let density_max = ceiling(&curves);

    let n = LUT_SIZE;
    // Blue outermost and red innermost, so the red axis varies fastest. See
    // `Baked::lut`: this is the layout a 3D texture upload wants, and getting
    // it backwards transposes red and blue in the finished picture.
    let cube = |max: f32, f: &dyn Fn([f32; 3]) -> [f32; 3]| {
        let mut out = Vec::with_capacity(n * n * n);
        for b in 0..n {
            for g in 0..n {
                for r in 0..n {
                    let step = max / (n - 1) as f32;
                    out.push(f([r as f32 * step, g as f32 * step, b as f32 * step]));
                }
            }
        }
        out
    };

    let (lut, paper) = match recipe.print {
        None => {
            let viewing = Viewing::new(&film.viewing_illuminant);
            (
                cube(density_max, &|d| viewing.to_srgb(&film.transmittance(d))),
                None,
            )
        }
        Some(paper) => {
            let viewing = Viewing::new(&paper.viewing_illuminant);
            let curves = resample(&paper.density_curves);
            let paper_max = ceiling(&curves);
            let lut = cube(density_max, &|d| {
                paper_exposure(film, paper, d).map(|raw| (raw + 1e-10).log10())
            });
            let paper = Paper {
                balance: print_balance(film, paper, recipe.exposure_ev),
                log_min: paper.log_exposure_min,
                log_max: paper.log_exposure_max,
                lut: cube(paper_max, &|d| viewing.to_srgb(&paper.transmittance(d))),
                density_max: paper_max,
                curves,
            };
            (lut, Some(paper))
        }
    };

    Baked {
        exposure_matrix: matrix,
        curves,
        curve_rows,
        push_stations,
        curve_log_min: film.log_exposure_min,
        curve_log_max: film.log_exposure_max,
        lut,
        paper,
        density_max,
        lut_size: n,
    }
}

/// A curve at `CURVE_SAMPLES`, uniform over the same domain it came in on.
fn resample(curve: &[[f32; 3]]) -> Vec<[f32; 3]> {
    if curve.len() == CURVE_SAMPLES {
        return curve.to_vec();
    }
    (0..CURVE_SAMPLES)
        .map(|i| {
            let at = i as f32 / (CURVE_SAMPLES - 1) as f32;
            sample_curve(curve, 0.0, 1.0, [at; 3])
        })
        .collect()
}

/// The deepest density in a set of curves, floored so a lookup over it has
/// a width.
fn ceiling(curves: &[[f32; 3]]) -> f32 {
    curves
        .iter()
        .flat_map(|row| row.iter())
        .fold(0.0f32, |a, &b| a.max(b))
        .max(1e-3)
}

fn mean(s: &Spectrum) -> f32 {
    s.iter().sum::<f32>() / SPECTRUM as f32
}

/// The log exposure at which the paper's average curve reaches `density`.
fn invert_mean_curve(paper: &Profile, density: f32) -> f32 {
    let curves = &paper.density_curves;
    let last = curves.len() - 1;
    let span = paper.log_exposure_max - paper.log_exposure_min;
    let at = |i: usize| (curves[i][0] + curves[i][1] + curves[i][2]) / 3.0;
    let log_at = |i: usize| paper.log_exposure_min + span * i as f32 / last as f32;

    // A paper is negative-working, so its mean curve rises. Walk it rather
    // than binary-search: 256 samples is nothing, and a linear scan is correct
    // even where the curve is flat, which a bisection is not.
    let ascending = at(last) >= at(0);
    for i in 0..last {
        let (lo, hi) = (at(i), at(i + 1));
        let brackets = if ascending {
            lo <= density && density <= hi
        } else {
            hi <= density && density <= lo
        };
        if brackets && (hi - lo).abs() > f32::EPSILON {
            let f = (density - lo) / (hi - lo);
            return log_at(i) + (log_at(i + 1) - log_at(i)) * f;
        }
    }
    // Off the end of the curve: the nearest end is the honest answer, and it
    // keeps a badly-scaled contributed profile from producing a NaN that would
    // propagate silently through the whole LUT.
    if (density <= at(0)) == ascending {
        paper.log_exposure_min
    } else {
        paper.log_exposure_max
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::profile::CURVE_SAMPLES;

    fn profile(yaml: &str) -> Profile {
        Profile::parse(yaml).unwrap()
    }

    fn portra() -> Profile {
        profile(include_str!("../profiles/kodak_portra_400.yaml"))
    }

    fn endura() -> Profile {
        profile(include_str!("../profiles/kodak_portra_endura.yaml"))
    }

    fn kodachrome() -> Profile {
        profile(include_str!("../profiles/kodak_kodachrome_64.yaml"))
    }

    fn spread(rgb: [f32; 3]) -> f32 {
        rgb.iter().cloned().fold(f32::MIN, f32::max) - rgb.iter().cloned().fold(f32::MAX, f32::min)
    }

    #[test]
    fn the_exposure_matrix_is_diagonally_dominant() {
        // Each layer must respond most strongly to its own primary. A matrix
        // that failed this would mean the sensitivity table had been pasted in
        // the wrong channel order, which produces a picture that renders
        // perfectly and has its colours swapped.
        for film in [portra(), kodachrome()] {
            let m = exposure_matrix(&film);
            for layer in 0..3 {
                for channel in 0..3 {
                    if channel != layer {
                        assert!(
                            m[layer][layer] > m[layer][channel] * 3.0,
                            "{}: layer {layer} responds to channel {channel} too strongly: {m:?}",
                            film.stock
                        );
                    }
                }
            }
        }
    }

    #[test]
    fn a_reversal_stock_renders_a_positive() {
        let film = kodachrome();
        let baked = bake(&Recipe::new(&film, None));
        let shadow = baked.apply([0.02; 3]);
        let mid = baked.apply([MID_GREY; 3]);
        let highlight = baked.apply([0.8; 3]);
        assert!(
            shadow[1] < mid[1] && mid[1] < highlight[1],
            "{shadow:?} {mid:?} {highlight:?}"
        );
    }

    #[test]
    fn a_scanned_negative_is_inverted_and_orange() {
        // Not a defect: it is what a negative looks like, and rendering it any
        // other way would mean the print stage was silently applied.
        let film = portra();
        let baked = bake(&Recipe::new(&film, None));
        let shadow = baked.apply([0.02; 3]);
        let highlight = baked.apply([0.8; 3]);
        assert!(
            shadow[1] > highlight[1],
            "not inverted: {shadow:?} -> {highlight:?}"
        );
        let mid = baked.apply([MID_GREY; 3]);
        assert!(mid[0] > mid[2] * 4.0, "no orange mask: {mid:?}");
    }

    #[test]
    fn printing_a_negative_restores_the_picture() {
        // The property the whole print stage exists for: through the paper,
        // the same negative is the right way up and neutral again.
        let film = portra();
        let paper = endura();
        let baked = bake(&Recipe::new(&film, Some(&paper)));

        let shadow = baked.apply([0.02; 3]);
        let mid = baked.apply([MID_GREY; 3]);
        let highlight = baked.apply([0.8; 3]);
        assert!(
            shadow[1] < mid[1] && mid[1] < highlight[1],
            "print is not a positive: {shadow:?} {mid:?} {highlight:?}"
        );
        for grey in [shadow, mid, highlight] {
            assert!(
                spread(grey) < 0.06,
                "print of a neutral is not neutral: {grey:?}"
            );
        }
    }

    #[test]
    fn exposure_moves_the_print_the_way_it_moves_a_photograph() {
        // The photograph's exposure: the enlarger balanced at it, and the
        // scene brighter by it. Mid-grey stays where the balance puts it —
        // that is what the balance is for — so what a stop more does to a
        // print is lift everything either side of it along the paper's curve.
        let film = portra();
        let paper = endura();
        let brighter = bake(&Recipe {
            exposure_ev: 1.0,
            ..Recipe::new(&film, Some(&paper))
        });
        let base = bake(&Recipe::new(&film, Some(&paper)));
        let one_stop = Settings {
            exposure_ev: 1.0,
            ..Settings::default()
        };
        for v in [0.02f32, 0.6] {
            assert!(
                brighter.apply_at([v; 3], &one_stop)[1] > base.apply([v; 3])[1],
                "{v} did not print brighter a stop up"
            );
        }
        let (a, b) = (
            brighter.apply_at([MID_GREY; 3], &one_stop)[1],
            base.apply([MID_GREY; 3])[1],
        );
        assert!(
            (a - b).abs() < 1.0 / 255.0,
            "the balance let mid-grey move: {a} vs {b}"
        );
    }

    #[test]
    fn a_region_exposed_brighter_prints_brighter_than_the_enlarger_expects() {
        // TRACES: FR-DEV-3f
        // A layer's exposure is the scene's, not the enlarger's: the balance
        // stays where the photograph put it, so the region prints lighter by
        // more than the whole photograph would, which is what dodging at the
        // camera is.
        let film = portra();
        let paper = endura();
        let base = bake(&Recipe::new(&film, Some(&paper)));
        let rebalanced = bake(&Recipe {
            exposure_ev: 1.0,
            ..Recipe::new(&film, Some(&paper))
        });
        let one_stop = Settings {
            exposure_ev: 1.0,
            ..Settings::default()
        };
        let local = base.apply_at([MID_GREY; 3], &one_stop)[1];
        let global = rebalanced.apply_at([MID_GREY; 3], &one_stop)[1];
        assert!(local > base.apply([MID_GREY; 3])[1], "not brighter at all");
        assert!(
            local > global,
            "a region was rebalanced as though it were the whole print: {local} vs {global}"
        );
    }

    #[test]
    fn more_light_through_the_enlarger_darkens_the_print() {
        // TRACES: FR-DEV-3f
        // Paper is negative-working. Opening the enlarger a stop is burning
        // in, and a slider that brightened would be the wrong way round for
        // anyone who has printed.
        let film = portra();
        let paper = endura();
        let baked = bake(&Recipe::new(&film, Some(&paper)));
        let at = |stops: f32| {
            baked.apply_at(
                [MID_GREY; 3],
                &Settings {
                    print_exposure_ev: stops,
                    ..Settings::default()
                },
            )[1]
        };
        assert!(at(1.0) < at(0.0) && at(0.0) < at(-1.0));
    }

    #[test]
    fn a_push_on_a_row_is_the_measured_curve() {
        // TRACES: FR-DEV-3f
        // The rows are the measured processes, so at a row the table must be
        // that curve exactly, and between rows — density being linear in push
        // there — it must be `curves_at_push` to rounding.
        let film = profile(include_str!("../profiles/kodak_doublex.yaml"));
        let baked = bake(&Recipe::new(&film, None));
        assert_eq!(
            baked.curve_rows, 5,
            "Double-X measures five development times"
        );

        let span = film.log_exposure_max - film.log_exposure_min;
        let mut worst = 0.0f32;
        let stations = baked.push_stations.clone();
        let mut pushes: Vec<(f32, bool)> = stations.iter().map(|p| (*p, true)).collect();
        for k in 0..=16 {
            pushes.push((-1.0 + 4.0 * k as f32 / 16.0, false));
        }
        for (push, on_row) in pushes {
            let exact = film.curves_at_push(push);
            for i in (0..exact.len()).step_by(7) {
                let log = film.log_exposure_min + span * i as f32 / (exact.len() - 1) as f32;
                let got = baked.sample_curves([log; 3], push);
                for c in 0..3 {
                    let err = (got[c] - exact[i][c]).abs();
                    if on_row {
                        assert!(err < 1e-4, "push {push} is a row but misses it by {err}");
                    }
                    worst = worst.max(err);
                }
            }
        }
        assert!(worst < 1e-3, "between rows the density is off by {worst}");
    }

    #[test]
    fn a_print_exposure_is_exact_at_any_setting() {
        // TRACES: FR-DEV-3f
        // The enlarger's exposure is added between the two lookups rather than
        // baked into either, so no setting is nearer the tables than another.
        // Compared against the chain evaluated spectrally, end to end, at
        // settings chosen off every half and whole stop.
        let film = portra();
        let paper = endura();
        let baked = bake(&Recipe::new(&film, Some(&paper)));
        let offsets = print_balance(&film, &paper, 0.0);
        let viewing = Viewing::new(&paper.viewing_illuminant);

        let mut worst = 0.0f32;
        for stops in [-2.3f32, -0.6, 0.0, 0.35, 1.7] {
            for i in 0..14 {
                let v = 0.004 * 2f32.powf(i as f32 * 0.6);
                let rgb = [v, v * 0.8, v * 1.1];
                let mut log_exposure = [0.0f32; 3];
                for (l, slot) in log_exposure.iter_mut().enumerate() {
                    let m = baked.exposure_matrix[l];
                    *slot =
                        ((m[0] * rgb[0] + m[1] * rgb[1] + m[2] * rgb[2]).max(0.0) + 1e-10).log10();
                }
                let raw = paper_exposure(&film, &paper, film.density_at(log_exposure));
                let paper_log =
                    [0, 1, 2].map(|l| (raw[l] + 1e-10).log10() + offsets[l] + stops * 2f32.log10());
                let exact = viewing.to_srgb(&paper.transmittance(paper.density_at(paper_log)));
                let approx = baked.apply_at(
                    rgb,
                    &Settings {
                        print_exposure_ev: stops,
                        ..Settings::default()
                    },
                );
                for c in 0..3 {
                    worst = worst.max((exact[c] - approx[c]).abs());
                }
            }
        }
        assert!(
            worst < 1.0 / 255.0,
            "the print misses the spectral chain by {worst}"
        );
    }

    #[test]
    fn the_lut_holds_its_error_under_a_code_value() {
        // The claim LUT_SIZE is chosen on. Compared against the same chain
        // evaluated exactly, so it measures interpolation error and nothing
        // else.
        let film = kodachrome();
        let baked = bake(&Recipe::new(&film, None));
        let viewing = Viewing::new(&film.viewing_illuminant);

        let mut worst = 0.0f32;
        for i in 0..40 {
            for j in 0..40 {
                let rgb = [
                    i as f32 / 39.0,
                    j as f32 / 39.0,
                    ((i + j) % 40) as f32 / 39.0,
                ];
                let mut log_exposure = [0.0f32; 3];
                for (l, slot) in log_exposure.iter_mut().enumerate() {
                    let m = baked.exposure_matrix[l];
                    *slot =
                        ((m[0] * rgb[0] + m[1] * rgb[1] + m[2] * rgb[2]).max(0.0) + 1e-10).log10();
                }
                let exact = viewing.to_srgb(&film.transmittance(film.density_at(log_exposure)));
                let approx = baked.apply(rgb);
                for c in 0..3 {
                    worst = worst.max((exact[c] - approx[c]).abs());
                }
            }
        }
        assert!(
            worst < 1.0 / 255.0,
            "worst LUT error {worst} exceeds one code value"
        );
    }

    #[test]
    fn the_lut_stores_red_along_its_fastest_axis() {
        // The layout a 3D texture upload expects, and the one bug this whole
        // decomposition is most exposed to: fill it the other way round and
        // the picture comes back with red and blue transposed -- entirely
        // plausible-looking, and wrong. Asserted here rather than only in the
        // end-to-end GPU test, because that one needs a device and this one
        // does not.
        let film = kodachrome();
        let baked = bake(&Recipe::new(&film, None));
        let n = baked.lut_size;

        // Step one along each axis from the origin, and check that the entry
        // found is the one the *density* moved along that axis should give.
        let viewing = Viewing::new(&film.viewing_illuminant);
        let step = baked.density_max / (n - 1) as f32;
        for (axis, offset) in [(0usize, 1usize), (1, n), (2, n * n)] {
            let mut density = [0.0f32; 3];
            density[axis] = step;
            let expected = viewing.to_srgb(&film.transmittance(density));
            let stored = baked.lut[offset];
            for c in 0..3 {
                assert!(
                    (stored[c] - expected[c]).abs() < 1e-4,
                    "axis {axis} is not at stride {offset}: stored {stored:?}, \
                     the density one step along that axis gives {expected:?}"
                );
            }
        }
    }

    #[test]
    fn the_lut_is_the_size_it_says_it_is() {
        let film = kodachrome();
        let baked = bake(&Recipe::new(&film, None));
        assert_eq!(baked.lut.len(), LUT_SIZE * LUT_SIZE * LUT_SIZE);
        assert_eq!(baked.curves.len(), CURVE_SAMPLES);
        assert_eq!(baked.curve_rows, 1);
        assert!(baked.paper.is_none());

        // A print has a second lookup and a curve of its own; a development
        // series a row per push. Neither is inferred from the other.
        let negative = portra();
        let paper = endura();
        let printed = bake(&Recipe::new(&negative, Some(&paper)));
        let print = printed
            .paper
            .as_ref()
            .expect("a printed negative has a paper");
        assert_eq!(print.lut.len(), LUT_SIZE.pow(3));
        assert_eq!(print.curves.len(), CURVE_SAMPLES);

        let pushable = profile(include_str!("../profiles/kodak_doublex.yaml"));
        let rows = bake(&Recipe::new(&pushable, None));
        assert_eq!(rows.curve_rows, pushable.development_times.len());
        assert_eq!(rows.push_stations.len(), rows.curve_rows);
        assert_eq!(rows.curves.len(), rows.curve_rows * CURVE_SAMPLES);
    }
}