DarkRoom/core/dr-pipeline/src/ops/film_sim.rs

//! TRACES: FR-DEV-3f
//! Film simulation — the stock renders the picture.
//!
//! # Why this one replaces the base curve
//!
//! [`crate::ops`]' other nodes adjust a picture. This one *makes* it. The base
//! curve exists because sensor data is scene-referred and nothing anybody looks
//! at is (FR-DEV-3e); a film stock's characteristic curve does the same job,
//! from measurements, with a toe and a shoulder that were coated onto acetate
//! rather than drawn. Running both renders the image twice — the camera's
//! JPEG-ish rendering, and then a film's rendering of that — which is not what
//! either is for and looks like neither.
//!
//! So this node declares [`Operation::renders`], and the composer answers by
//! emitting neither the base curve nor the camera matrix. Both jobs move here:
//! the fragment takes camera RGB, converts it to linear sRGB itself with the
//! matrix already in the uniform block, and returns linear sRGB. That is a
//! contract worth stating plainly, because a node that got half of it wrong
//! would produce a picture that renders perfectly and is wrong everywhere.
//!
//! # Why the tables are not parameters
//!
//! For the same reason [`crate::ops::vignetting`]'s coefficients are not: they
//! are measurements of a physical thing, not something a slider moves. The
//! sliders here are exposure and print exposure, which are what a photographer
//! and a printer actually control. `dr-film` turns a stock plus those two
//! numbers into [`FilmTables`]; this node knows only the layout.
//!
//! Declared as a plain struct here rather than imported, so that dr-pipeline
//! keeps its no-dependency property (ARCH §6.5a) exactly as `vignetting` does
//! with `Pa`.
use std::sync::{Arc, LazyLock};

use crate::descriptor::{Attribute, LocalizedKey, OpDescriptor, OpId, ParamDescriptor, ParamId};
use crate::operation::{Operation, Uniform};

pub const ID: OpId = OpId("film_sim");
pub const EXPOSURE: ParamId = ParamId("exposure");
pub const PRINT_EXPOSURE: ParamId = ParamId("print_exposure");
pub const PUSH: ParamId = ParamId("push");
pub const FORMAT: ParamId = ParamId("format");

/// TRACES: FR-DEV-3f
/// The frames a photograph can be simulated on, smallest first.
///
/// A genuinely fixed list, unlike the stocks: nobody invents a film format, so
/// this is a declared `enum` parameter and gets its control, its place in the
/// sidecar and its undo step for free. The *sizes* live in `dr_film::Format`;
/// this crate carries only the names, in the same order.
static FORMATS: [LocalizedKey; 6] = [
    LocalizedKey("param.film_sim.format.35mm"),
    LocalizedKey("param.film_sim.format.645"),
    LocalizedKey("param.film_sim.format.6x6"),
    LocalizedKey("param.film_sim.format.6x7"),
    LocalizedKey("param.film_sim.format.4x5"),
    LocalizedKey("param.film_sim.format.8x10"),
];

/// How many samples a characteristic curve carries.
///
/// Must agree with `dr_film::profile::CURVE_SAMPLES`. Restated rather than
/// imported because importing it is exactly the dependency this crate does not
/// take; [`FilmTables::is_well_formed`] is what stops the two drifting.
pub const CURVE_SAMPLES: usize = 256;

/// The uniform field names the fragment reads the exposure matrix from.
///
/// A table rather than a formatted string, because a `Uniform`'s name is
/// `&'static str`: building one per composition would mean leaking a string
/// every time a slider moved.
static MATRIX_FIELDS: [[&str; 3]; 3] = [
    ["m00", "m01", "m02"],
    ["m10", "m11", "m12"],
    ["m20", "m21", "m22"],
];

static DESCRIPTOR: LazyLock<Arc<OpDescriptor>> = LazyLock::new(|| {
    Arc::new(OpDescriptor {
        // Tone and colour both, and not `Effect`: a stock is not something applied
        // on top of a photograph, it is what the photograph was made on.
        attributes: vec![Attribute::Tone, Attribute::Colour],
        id: ID,
        label: LocalizedKey("op.film_sim"),
        params: vec![
            ParamDescriptor::stops("exposure", "param.film_sim.exposure", -3.0, 3.0),
            ParamDescriptor::stops("print_exposure", "param.film_sim.print_exposure", -3.0, 3.0),
            // TRACES: FR-DEV-3f
            // Development, in stops of push. Bounded by what the manufacturers
            // actually published: Double-X's measured axis spans about -1 to +2,
            // and beyond a range like that a curve would have to be invented.
            ParamDescriptor::stops("push", "param.film_sim.push", -1.0, 3.0),
            // TRACES: FR-DEV-3f
            // Which frame this was taken on — the half of the enlargement a
            // photograph cannot supply. A crystal is a fixed size in micrometres,
            // so how grainy a picture looks is film size against output size, and
            // the same emulsion on 4x5 renders about three times smoother than on
            // 35mm at the same print.
            ParamDescriptor::choice("format", "param.film_sim.format", FORMATS.to_vec()),
        ],
    })
});

/// A stock reduced to what a shader runs, as `dr-film` bakes it.
///
/// Layout is the contract between the two crates, so it is written down here
/// and checked rather than assumed:
///
/// - `exposure_matrix[l][c]` — layer `l`'s response to linear sRGB channel `c`.
/// - `curves` — `CURVE_SAMPLES` density triples, uniform over
///   `[curve_log_min, curve_log_max]`.
/// - `lut` — `lut_size³` linear sRGB triples, uniform over `[0, density_max]`
///   on each axis, with the **red axis varying fastest**: index
///   `(b * size + g) * size + r`. That is the order a 3D texture upload
///   expects, so the consumer hands the slice straight to the driver. Filling
///   it the other way round transposes red and blue in the finished picture —
///   which is a plausible photograph of the wrong colour, and which the unit
///   tests on both sides of this seam happily pass, because each side is
///   internally consistent. `dr-film` pins it; `dr-gpu`'s `film_sim` test
///   catches it end to end.
#[derive(Debug, Clone, PartialEq)]
pub struct FilmTables {
    pub exposure_matrix: [[f32; 3]; 3],
    pub curves: Vec<[f32; 3]>,
    pub curve_log_min: f32,
    pub curve_log_max: f32,
    pub lut: Vec<[f32; 3]>,
    pub density_max: f32,
    pub lut_size: usize,
    /// TRACES: FR-DEV-3f
    /// Grains in one pixel's patch of film, per layer, with the density
    /// ceiling and uniformity the variance is taken against. Zero particles
    /// means no grain, which is how the control is turned off.
    pub grain_particles: [f32; 3],
    pub grain_density_max: [f32; 3],
    pub grain_uniformity: f32,
}

impl FilmTables {
    /// Whether these tables are the shape the shader will index them at.
    ///
    /// Checked on the way in, because the failure otherwise is a shader
    /// sampling past the end of a texture: undefined, silent, and different on
    /// every driver.
    pub fn is_well_formed(&self) -> bool {
        self.curves.len() == CURVE_SAMPLES
            && self.lut_size >= 2
            && self.lut.len() == self.lut_size.pow(3)
            && self.density_max > 0.0
            && self.curve_log_max > self.curve_log_min
    }
}

/// TRACES: FR-DEV-3f
#[derive(Debug, Default, Clone)]
pub struct FilmSim {
    exposure: f32,
    print_exposure: f32,
    push: f32,
    /// Index into `FORMATS`. Zero is 35 mm, which is the neutral choice.
    format: f32,
    tables: Option<FilmTables>,
}

impl FilmSim {
    pub fn new() -> Self {
        Self::default()
    }

    /// Load a baked stock, or clear it.
    ///
    /// Malformed tables are refused rather than stored: an operation that is
    /// active but cannot be indexed is worse than one that is off, because the
    /// first renders garbage and the second renders the photograph.
    pub fn set_tables(&mut self, tables: Option<FilmTables>) {
        match tables {
            Some(t) if !t.is_well_formed() => {
                log::error!(
                    "film_sim: refusing malformed tables ({} curve samples, {} lut entries at size {})",
                    t.curves.len(),
                    t.lut.len(),
                    t.lut_size
                );
                self.tables = None;
            }
            other => self.tables = other,
        }
    }

    /// The loaded stock's tables, for whoever has to upload them.
    pub fn tables(&self) -> Option<&FilmTables> {
        self.tables.as_ref()
    }
}

impl Operation for FilmSim {
    fn descriptor(&self) -> Arc<OpDescriptor> {
        DESCRIPTOR.clone()
    }

    fn set_param(&mut self, id: ParamId, value: f32) {
        match id {
            EXPOSURE => self.exposure = value,
            PRINT_EXPOSURE => self.print_exposure = value,
            PUSH => self.push = value,
            FORMAT => self.format = value,
            _ => log::warn!("film_sim: unknown parameter {id}"),
        }
    }

    fn param(&self, id: ParamId) -> f32 {
        match id {
            EXPOSURE => self.exposure,
            PRINT_EXPOSURE => self.print_exposure,
            PUSH => self.push,
            FORMAT => self.format,
            _ => 0.0,
        }
    }

    /// Active exactly when a stock is loaded.
    ///
    /// Not "when a slider has moved", which is the rule everywhere else and
    /// would be wrong here: a stock at zero exposure compensation is the whole
    /// point of choosing it, and a node that went quiet at its defaults would
    /// mean picking a film did nothing until you also nudged something.
    fn is_active(&self) -> bool {
        self.tables.is_some()
    }

    /// This node renders; the camera's own rendering must not also run.
    fn renders(&self) -> bool {
        true
    }

    fn set_film_tables(&mut self, tables: Option<&FilmTables>) {
        self.set_tables(tables.cloned());
    }

    fn uniforms(&self) -> Vec<Uniform> {
        let Some(t) = &self.tables else {
            return Vec::new();
        };
        let m = t.exposure_matrix;
        // Exposure rides in the matrix on the CPU when the stock is baked, so
        // what is left here is the *shader's* copy of the same nine numbers.
        // Spelled out one at a time because a uniform is a named `f32` in this
        // pipeline and a matrix would be a second kind of thing for one caller.
        let mut out = Vec::with_capacity(MATRIX_FIELDS.len() + 5);
        for (l, row) in m.iter().enumerate() {
            for (c, v) in row.iter().enumerate() {
                out.push(Uniform {
                    name: MATRIX_FIELDS[l][c],
                    value: *v,
                });
            }
        }
        for (l, name) in ["gn0", "gn1", "gn2"].into_iter().enumerate() {
            out.push(Uniform {
                name,
                value: t.grain_particles[l],
            });
        }
        for (l, name) in ["gd0", "gd1", "gd2"].into_iter().enumerate() {
            out.push(Uniform {
                name,
                value: t.grain_density_max[l],
            });
        }
        out.push(Uniform {
            name: "grain_u",
            value: t.grain_uniformity,
        });
        out.push(Uniform {
            name: "log_min",
            value: t.curve_log_min,
        });
        out.push(Uniform {
            name: "log_max",
            value: t.curve_log_max,
        });
        out.push(Uniform {
            name: "density_max",
            value: t.density_max,
        });
        out.push(Uniform {
            name: "lut_size",
            value: t.lut_size as f32,
        });
        out.push(Uniform {
            name: "print_exposure",
            value: self.print_exposure,
        });
        out
    }

    fn wgsl_body(&self) -> String {
        // Filtered by hand rather than through a sampler, which is what the
        // framing prologue already does for the source: this pipeline has no
        // sampler binding, and adding one to interpolate two lookups would
        // cost a binding in every shader whether or not a film is loaded.
        "\
// Camera RGB to linear sRGB. The film's exposure matrix is defined against
// sRGB primaries, and this node has taken over the conversion the composer
// would otherwise have emitted at the end — see `Operation::renders`.
let scene = vec3<f32>(
    dot(u.cam_to_srgb_0.rgb, c),
    dot(u.cam_to_srgb_1.rgb, c),
    dot(u.cam_to_srgb_2.rgb, c),
);

// What each emulsion layer was exposed to. A matrix, exactly: the scene
// spectrum reconstructed from an sRGB triple is linear in that triple, so the
// integral over wavelength collapsed into these nine numbers when the stock
// was baked.
let exposure = vec3<f32>(
    dot(vec3<f32>(m00, m01, m02), scene),
    dot(vec3<f32>(m10, m11, m12), scene),
    dot(vec3<f32>(m20, m21, m22), scene),
);
// 1e-10 rather than a clamp to zero: a black pixel has to land somewhere on
// the curve, and the toe is where it belongs.
let log_exposure = log10(max(exposure, vec3<f32>(0.0)) + 1e-10);

// The characteristic curve: what density each layer develops to. Clamped, not
// extrapolated — past the shoulder a real emulsion stops responding, and
// extrapolating would turn a blown highlight into a colour cast that grows the
// more it is overexposed.
let density = film_curve(clamp((log_exposure - log_min) / (log_max - log_min),
                               vec3<f32>(0.0), vec3<f32>(1.0)));

// TRACES: FR-DEV-3f
// Grain, on the density and before the dye.
//
// That order is the physical one and it is not cosmetic: grain is silver that
// did or did not develop, so it perturbs *density*, and the dye absorbs
// through whatever density resulted. Adding noise to the finished colour --
// which is what an effect does -- tints the highlights wrong, because that
// noise never passes through the dye at all.
let grained = film_grain(density, source_px,
                         vec3<f32>(gn0, gn1, gn2),
                         vec3<f32>(gd0, gd1, gd2),
                         grain_u);

// Dye absorption, the print through the negative, the paper, the viewing
// illuminant and the chromatic adaptation — all of which take exactly three
// numbers in, which is why they fit in one lookup.
c = film_lut(clamp(grained / density_max, vec3<f32>(0.0), vec3<f32>(1.0)), lut_size);"
            .into()
    }

    fn helpers(&self) -> &'static [crate::operation::Helper] {
        &HELPERS
    }
}

static HELPERS: [crate::operation::Helper; 5] = [
    crate::operation::Helper {
        name: "film_hash",
        source: "\
// A hash, not a random number generator: the same pixel of the same frame has
// to grain the same way every time it is drawn, or the picture would crawl
// while nobody was editing it. Seeded from a position, so it is reproducible
// by construction rather than by holding state between frames.
//
// Two decorrelated uniforms come out, which is what a Gaussian needs.
fn film_hash(p: vec2<f32>, layer: u32) -> vec2<f32> {
    var h = u32(i32(floor(p.x))) * 73856093u
          ^ u32(i32(floor(p.y))) * 19349663u
          ^ (layer + 1u) * 83492791u;
    h = h ^ (h >> 16u);
    h = h * 2246822519u;
    h = h ^ (h >> 13u);
    h = h * 3266489917u;
    let a = h ^ (h >> 16u);
    var g = a * 747796405u + 2891336453u;
    g = ((g >> ((g >> 28u) + 4u)) ^ g) * 277803737u;
    let b = g ^ (g >> 22u);
    // Open interval: a zero would send the logarithm below to infinity.
    return vec2<f32>(
        max(f32(a) * 2.3283064e-10, 1e-7),
        max(f32(b) * 2.3283064e-10, 1e-7)
    );
}",
    },
    crate::operation::Helper {
        name: "film_grain",
        source: "\
// Developed density, with the variance a count of silver grains actually has.
//
//     mean     = D
//     variance = D * (Dmax - u * D) / N
//
// N is grains *per pixel*, so the entire scale dependence sits in that uniform
// and none of it is here: a zoomed-out pixel covers more film, averages more
// grains, and comes out smoother with nothing filtered.
//
// A Gaussian with the exact first two moments, rather than the exact compound
// Poisson-Binomial the silver actually follows. The two agree wherever grain
// is visible; the real one is skewed only in the deep toe, where the density
// is near zero and so is its variance. Sampling it properly would cost tens of
// draws per layer per pixel to change nothing anyone can see.
// Takes its parameters rather than reading uniforms, and must: the composer
// prefixes a uniform with its operation's id and rewrites the references
// *inside a fragment body only*. Helpers are shared between operations and
// deduplicated by name, so a bare `gn0` here is an identifier that exists in
// no shader. `film_lut` below takes its size for the same reason.
fn film_grain(
    density: vec3<f32>,
    at: vec2<f32>,
    n: vec3<f32>,
    dmax: vec3<f32>,
    uniformity: f32,
) -> vec3<f32> {
    var out = density;
    for (var l = 0u; l < 3u; l = l + 1u) {
        if (n[l] <= 0.0) {
            continue;
        }
        let d = clamp(density[l], 0.0, dmax[l]);
        let variance = d * (dmax[l] - uniformity * d) / n[l];
        if (variance <= 0.0) {
            continue;
        }
        let u = film_hash(at, l);
        // Box-Muller. Half the pair is discarded rather than carried: the next
        // layer wants a seed of its own, not this one's leftover.
        let z = sqrt(-2.0 * log(u.x)) * cos(6.2831853 * u.y);
        // Clamped, not wrapped: a negative density is not a colour, and the
        // ceiling is the most silver this emulsion has to develop.
        out[l] = clamp(d + z * sqrt(variance), 0.0, dmax[l]);
    }
    return out;
}",
    },
    crate::operation::Helper {
        name: "log10",
        source: "\
// WGSL has no log10, and `log2(x) * log10(2)` is the cheap identity for it.
fn log10(v: vec3<f32>) -> vec3<f32> {
    return log2(v) * 0.30103;
}",
    },
    crate::operation::Helper {
        name: "film_curve",
        source: "\
// Three characteristic curves, sampled from a 256-wide texture and
// interpolated by hand. `t` is already normalised to the curve's domain.
fn film_curve(t: vec3<f32>) -> vec3<f32> {
    let samples = u32(textureDimensions(film_curves).x);
    let last = f32(samples - 1u);
    var out = vec3<f32>(0.0);
    for (var ch = 0u; ch < 3u; ch = ch + 1u) {
        let x = t[ch] * last;
        let i = min(u32(floor(x)), samples - 2u);
        let f = x - f32(i);
        let a = textureLoad(film_curves, vec2<i32>(i32(i), 0), 0);
        let b = textureLoad(film_curves, vec2<i32>(i32(i) + 1, 0), 0);
        out[ch] = mix(a[ch], b[ch], f);
    }
    return out;
}",
    },
    crate::operation::Helper {
        name: "film_lut",
        source: "\
// Trilinear interpolation of the density lookup, by hand for the same reason
// the curve above is: there is no sampler bound, and the eight loads are
// cache-neighbours.
fn film_lut(t: vec3<f32>, size: f32) -> vec3<f32> {
    let n = i32(size);
    let x = t * (size - 1.0);
    let base = min(vec3<i32>(floor(x)), vec3<i32>(n - 2));
    let f = x - vec3<f32>(base);

    var out = vec3<f32>(0.0);
    for (var dx = 0; dx < 2; dx = dx + 1) {
        let wx = select(1.0 - f.x, f.x, dx == 1);
        for (var dy = 0; dy < 2; dy = dy + 1) {
            let wy = select(1.0 - f.y, f.y, dy == 1);
            for (var dz = 0; dz < 2; dz = dz + 1) {
                let wz = select(1.0 - f.z, f.z, dz == 1);
                let p = base + vec3<i32>(dx, dy, dz);
                out = out + wx * wy * wz
                    * textureLoad(film_lut_texture, p, 0).rgb;
            }
        }
    }
    return out;
}",
    },
];

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

    fn tables() -> FilmTables {
        FilmTables {
            exposure_matrix: [[5.0, 0.5, 0.2], [0.1, 5.0, 0.3], [0.2, 0.5, 4.0]],
            curves: vec![[0.0, 0.0, 0.0]; CURVE_SAMPLES],
            curve_log_min: -3.0,
            curve_log_max: 4.0,
            lut: vec![[0.5, 0.5, 0.5]; 32 * 32 * 32],
            density_max: 3.0,
            lut_size: 32,
            grain_particles: [0.0; 3],
            grain_density_max: [3.0; 3],
            grain_uniformity: 0.97,
        }
    }

    #[test]
    fn it_starts_inactive() {
        assert!(!FilmSim::new().is_active());
    }

    #[test]
    fn loading_a_stock_is_what_turns_it_on() {
        // Not a moved slider, which is the rule for every other node. Choosing
        // a film has to do something on its own, or picking one would appear
        // to be broken until you also nudged the exposure.
        let mut op = FilmSim::new();
        op.set_tables(Some(tables()));
        assert!(op.is_active());
        op.set_tables(None);
        assert!(!op.is_active());
    }

    #[test]
    fn malformed_tables_are_refused_rather_than_stored() {
        // The alternative is a shader indexing past the end of a texture,
        // which is undefined, silent, and different on every driver.
        let mut op = FilmSim::new();
        let mut bad = tables();
        bad.lut.truncate(10);
        op.set_tables(Some(bad));
        assert!(!op.is_active(), "malformed tables were accepted");
    }

    #[test]
    fn a_short_curve_is_refused_too() {
        let mut op = FilmSim::new();
        let mut bad = tables();
        bad.curves.truncate(CURVE_SAMPLES - 1);
        op.set_tables(Some(bad));
        assert!(!op.is_active());
    }

    #[test]
    fn it_declares_itself_a_rendering_transform() {
        // The whole reason the composer skips the base curve and the camera
        // matrix. If this ever returned false the picture would be rendered
        // twice and converted twice, which looks like a colour management bug
        // a long way from here.
        assert!(FilmSim::new().renders());
    }

    #[test]
    fn the_matrix_reaches_the_shader_in_the_order_the_fragment_reads_it() {
        // `m01` must be layer 0's response to sRGB green. A transposed matrix
        // compiles, runs, and swaps the picture's colours.
        let mut op = FilmSim::new();
        op.set_tables(Some(tables()));
        let uniforms = op.uniforms();
        let named = |n: &str| uniforms.iter().find(|u| u.name == n).unwrap().value;
        assert_eq!(named("m01"), 0.5);
        assert_eq!(named("m10"), 0.1);
        assert_eq!(named("m22"), 4.0);
    }

    #[test]
    fn an_inactive_node_publishes_no_uniforms() {
        assert!(FilmSim::new().uniforms().is_empty());
    }

    #[test]
    fn the_fragment_converts_out_of_camera_space_itself() {
        // It has to: it has taken over the conversion the composer would
        // otherwise emit at the end.
        let mut op = FilmSim::new();
        op.set_tables(Some(tables()));
        let wgsl = op.wgsl_body();
        assert!(wgsl.contains("cam_to_srgb_0"), "{wgsl}");
    }

    #[test]
    fn every_helper_defines_the_function_it_names() {
        for h in HELPERS {
            assert!(h.source.contains(&format!("fn {}(", h.name)), "{}", h.name);
        }
    }
}