Merge remote-tracking branch 'origin/master' into worktree-spot-removal

# Conflicts:
#	docs/traceability.md

core/dr-film/src/boolean_grain.rs

@@ -0,0 +1,294 @@
+//! TRACES: FR-DEV-3f
+//! Grain as a Boolean model — grains that overlap, rather than noise that does
+//! not.
+//!
+//! # Why the counting model was not enough
+//!
+//! [`crate::grain`] gets the *variance* of a developed density right: a count
+//! of independent yes/no events, `D(Dmax − uD)/N`, calibrated from published
+//! granularity. What it cannot get right is the **structure**, because it
+//! treats every pixel as an independent draw.
+//!
+//! Measured at 35 mm, a pixel of a 5472-wide frame covers about 6.6 µm and
+//! holds some 300 crystals of ~370 nm. Three hundred independent events per
+//! pixel average almost flat, and the little that survives has no spatial
+//! extent — which is exactly why it reads as sensor noise rather than as film.
+//!
+//! Real grain is visible because it *clumps*. A crystal is far smaller than a
+//! pixel, but crystals overlap into structures that are not, and those survive
+//! the filtering that averages independent noise away.
+//!
+//! # The model
+//!
+//! Newson, Delon & Galerne, *A Stochastic Film Grain Model for
+//! Resolution-Independent Rendering* (Computer Graphics Forum, 2017).
+//!
+//! Grain centres are a Poisson process of local intensity `λ(y)`; each centre
+//! carries a disc. The developed film is the **union** of those discs, and a
+//! point is opaque exactly when some disc covers it. Because a disc covers a
+//! whole neighbourhood, nearby points are *correlated* — and that correlation
+//! is the clumping, which arrives for free rather than being added.
+//!
+//! # Why it couples to density and not to a grey level
+//!
+//! The paper drives `λ` from an image's grey level, because it renders grain
+//! onto a finished picture. We are not doing that: we have a *density* per
+//! layer, from a measured characteristic curve, and a dye that absorbs through
+//! it.
+//!
+//! The two meet exactly. In a Boolean model the chance a point is left
+//! uncovered is
+//!
+//! ```text
+//! P(uncovered) = exp(−λ · E[A])
+//! ```
+//!
+//! which is Beer–Lambert. So the model's coverage *is* optical density, and
+//!
+//! ```text
+//! λ = D · ln(10) / E[A]
+//! ```
+//!
+//! puts our measured densities straight into it — with `E[A]`, the mean grain
+//! area, already computed from published RMS granularity in
+//! [`crate::grain`]. Nothing here is tuned by eye.
+//!
+//! # What it costs
+//!
+//! Monte Carlo per pixel, against the counting model's two hashes and a square
+//! root. The shader form stores no grains: space is cut into cells, a
+//! generator is seeded from each cell's index, and only the cells a sample
+//! could reach are visited. That keeps it inside the fused pass — no
+//! neighbouring *pixel* is read — but it is emphatically not free, and
+//! [`BooleanGrain::samples_for`] is where that trade is made explicit.
+
+use crate::grain::Grain;
+
+/// Mean grain area, in µm², recovered from the counting model's calibration.
+///
+/// The two models are the same emulsion seen two ways, so they must not
+/// disagree about how big a crystal is: `Grain` already inverts published RMS
+/// granularity for exactly this number, and taking it from there is what stops
+/// a Boolean render and a counting render describing different films.
+pub fn mean_grain_area_um2(grain: &Grain, pixel_size_um: f32) -> [f32; 3] {
+    let pixel_area = (pixel_size_um * pixel_size_um).max(1e-6);
+    grain.particles.map(|n| pixel_area / n.max(1e-6))
+}
+
+/// TRACES: FR-DEV-3f
+/// What the shader needs to render the Boolean model.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct BooleanGrain {
+    /// Grain radius per layer, in *pixels* at the current sampling scale.
+    ///
+    /// In pixels rather than micrometres because that is the unit the shader
+    /// works in, and converting once here keeps the conversion out of the
+    /// inner loop.
+    pub radius_px: [f32; 3],
+    /// `ln(10) / E[A]`, per layer: the factor taking a density to a Poisson
+    /// intensity. Precomputed because it is constant per bake and the shader
+    /// would otherwise recompute a logarithm per pixel per layer.
+    pub lambda_per_density: [f32; 3],
+    /// The density each layer saturates at, as in the counting model.
+    pub density_max: [f32; 3],
+    /// Monte Carlo samples per pixel.
+    pub samples: u32,
+    /// Standard deviation of the sampling kernel, in pixels.
+    ///
+    /// The pixel's own footprint: what a scanner or an eye integrates over.
+    /// Too small and the render is binary salt and pepper; too large and the
+    /// grain is blurred out of existence.
+    pub sigma_px: f32,
+}
+
+impl BooleanGrain {
+    /// Derive the parameters for a stock at a given sampling scale.
+    pub fn new(grain: &Grain, pixel_size_um: f32, samples: u32) -> Self {
+        let area = mean_grain_area_um2(grain, pixel_size_um);
+        let pixel_area = (pixel_size_um * pixel_size_um).max(1e-6);
+
+        let mut radius_px = [0.0f32; 3];
+        let mut lambda_per_density = [0.0f32; 3];
+        for l in 0..3 {
+            // A disc of this area, expressed as a fraction of a pixel.
+            let area_px = area[l] / pixel_area;
+            radius_px[l] = (area_px / std::f32::consts::PI).sqrt();
+            // Beer-Lambert, read backwards: coverage exp(-lambda*E[A]) is
+            // transmittance 10^-D, so lambda = D * ln(10) / E[A].
+            lambda_per_density[l] = std::f32::consts::LN_10 / area_px.max(1e-9);
+        }
+
+        Self {
+            radius_px,
+            lambda_per_density,
+            density_max: grain.density_max,
+            samples: samples.max(1),
+            // Half a pixel: the footprint of one sample of a sensor whose
+            // pixels abut. Wider would be a soft scanner, narrower a sharper
+            // one than exists.
+            sigma_px: 0.5,
+        }
+    }
+
+    /// How many Monte Carlo samples a given quality asks for.
+    ///
+    /// The estimator's own noise falls as `1/sqrt(N)`, so this trades one kind
+    /// of grain against another: too few samples and the *sampling* shows as a
+    /// second, wrong texture on top of the film's.
+    pub fn samples_for(quality: Quality) -> u32 {
+        match quality {
+            Quality::Preview => 16,
+            Quality::Export => 64,
+        }
+    }
+
+    /// Expected coverage at a density — what the render must average to.
+    ///
+    /// The Boolean model's mean is analytic even though its texture is not,
+    /// which is what makes it testable without rendering anything: whatever
+    /// the grain does locally, across a flat patch it has to come back to the
+    /// density the characteristic curve asked for.
+    pub fn expected_coverage(&self, layer: usize, density: f32) -> f32 {
+        let d = density.clamp(0.0, self.density_max[layer]);
+        1.0 - 10f32.powf(-d)
+    }
+}
+
+/// How hard to work at the Monte Carlo.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Quality {
+    /// Interactive. Some sampling noise, which at preview scale is hidden
+    /// under the grain it is sampling.
+    Preview,
+    /// Final render, where the sampling noise must be well below the grain.
+    Export,
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::profile::Profile;
+
+    fn portra() -> Profile {
+        Profile::parse(include_str!("../profiles/kodak_portra_400.yaml")).unwrap()
+    }
+
+    fn model(pixel_size_um: f32) -> BooleanGrain {
+        let g = Grain::for_pixel_size(&portra(), pixel_size_um);
+        BooleanGrain::new(&g, pixel_size_um, 16)
+    }
+
+    #[test]
+    fn coverage_is_beer_lambert() {
+        // The identity the whole coupling rests on: a Boolean model's uncovered
+        // fraction is exp(-lambda E[A]), and transmittance is 10^-D, so the
+        // model's coverage *is* the film's opacity. If this drifts, the grain
+        // is no longer rendering the density the curve asked for.
+        let m = model(6.6);
+        // Inside the layer's own range. Past its Dmax the coverage clamps —
+        // correctly, since a film cannot develop denser than its maximum — and
+        // an earlier version of this test probed 2.0 against a layer that
+        // reaches 1.798, then blamed the model for the clamp.
+        for d in [0.0f32, 0.3, 1.0, 1.7] {
+            let coverage = m.expected_coverage(1, d);
+            let transmittance = 1.0 - coverage;
+            assert!(
+                (transmittance - 10f32.powf(-d)).abs() < 1e-5,
+                "at density {d}: transmittance {transmittance}, expected {}",
+                10f32.powf(-d)
+            );
+        }
+    }
+
+    #[test]
+    fn clear_film_has_no_grains_and_fully_developed_film_is_nearly_solid() {
+        let m = model(6.6);
+        assert!(m.expected_coverage(1, 0.0) < 1e-6);
+        // At its own maximum, not at some density it never reaches: Portra's
+        // green layer tops out near 1.8, which transmits about 1.6% — dense,
+        // and not opaque. A film that went fully black would be one whose
+        // shadows carried no detail at all.
+        let dmax = m.density_max[1];
+        assert!(
+            m.expected_coverage(1, dmax) > 0.98,
+            "{}",
+            m.expected_coverage(1, dmax)
+        );
+    }
+
+    #[test]
+    fn coverage_clamps_at_the_layers_own_maximum() {
+        // The property the two tests above tripped over, asserted directly:
+        // asking for more density than the emulsion has gives the emulsion's
+        // own ceiling rather than extrapolating one.
+        let m = model(6.6);
+        let dmax = m.density_max[1];
+        assert_eq!(
+            m.expected_coverage(1, dmax),
+            m.expected_coverage(1, dmax + 5.0)
+        );
+    }
+
+    #[test]
+    fn the_two_models_describe_the_same_crystal() {
+        // The counting model and this one are one emulsion seen two ways. If
+        // they disagreed about grain size they would render as different
+        // films, and the difference would look like a modelling choice rather
+        // than the bug it is.
+        let px = 6.6;
+        let g = Grain::for_pixel_size(&portra(), px);
+        let area = mean_grain_area_um2(&g, px);
+        // Portra's green layer: ~0.14 um^2, about 370 nm across.
+        assert!(
+            (0.10..0.20).contains(&area[1]),
+            "grain area {} um^2 is not what the counting model calibrated",
+            area[1]
+        );
+        let m = BooleanGrain::new(&g, px, 16);
+        // And the radius in pixels must match that area at this scale.
+        let area_px = area[1] / (px * px);
+        let expect_r = (area_px / std::f32::consts::PI).sqrt();
+        assert!((m.radius_px[1] - expect_r).abs() < 1e-6);
+    }
+
+    #[test]
+    fn zooming_in_makes_the_grains_bigger_in_pixels() {
+        // Resolution independence, which is the paper's headline claim and the
+        // thing the counting model can only approximate: a grain is a fixed
+        // size *on the film*, so looking closer must resolve it, not merely
+        // reduce the variance.
+        let close = model(2.0);
+        let far = model(12.0);
+        assert!(
+            close.radius_px[1] > far.radius_px[1] * 3.0,
+            "close {} far {}",
+            close.radius_px[1],
+            far.radius_px[1]
+        );
+    }
+
+    #[test]
+    fn a_finer_stock_has_smaller_grains() {
+        let mut fine = portra();
+        let mut coarse = portra();
+        fine.rms_granularity = [4.0; 3];
+        coarse.rms_granularity = [16.0; 3];
+
+        let f = BooleanGrain::new(&Grain::for_pixel_size(&fine, 6.6), 6.6, 16);
+        let c = BooleanGrain::new(&Grain::for_pixel_size(&coarse, 6.6), 6.6, 16);
+        assert!(
+            c.radius_px[1] > f.radius_px[1],
+            "coarse {} is not larger than fine {}",
+            c.radius_px[1],
+            f.radius_px[1]
+        );
+    }
+
+    #[test]
+    fn export_samples_more_than_preview() {
+        assert!(
+            BooleanGrain::samples_for(Quality::Export)
+                > BooleanGrain::samples_for(Quality::Preview)
+        );
+    }
+}

core/dr-film/src/grain.rs

@@ -76,13 +76,65 @@ const RMS_APERTURE_AREA_UM2: f32 = std::f32::consts::PI * 24.0 * 24.0;
 /// The net density the granularity figure is quoted at.
 const RMS_REFERENCE_NET_DENSITY: f32 = 1.0;
 
-/// A 35 mm frame's width, in micrometres.
+/// TRACES: FR-DEV-3f
+/// The frame a photograph is being simulated on.
 ///
-/// What turns a pixel count into a grain size. A photograph has no inherent
-/// film format, so simulating one means choosing what the frame *would have
-/// been*; 35 mm is the choice that makes the numbers mean what a photographer
-/// expects, since published granularity and every intuition about how grainy a
-/// stock looks come from 35 mm.
+/// **Grain is a function of enlargement, and this is the half of it the
+/// photograph cannot supply.** A crystal is a fixed size in micrometres, so
+/// how grainy a picture looks depends entirely on how much the frame was
+/// magnified to make it — and that is film size against output size.
+///
+/// The same emulsion on 4x5 packs about 3,800 crystals into the pixel that
+/// holds 300 on 35 mm, so it renders roughly 3.5 times smoother at the same
+/// output size. Treating everything as 35 mm, as this did, made every
+/// photograph as grainy as the smallest common format.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Format {
+    Mm35,
+    Format645,
+    Format6x6,
+    Format6x7,
+    Sheet4x5,
+    Sheet8x10,
+}
+
+impl Format {
+    /// The formats, in the order the picker offers them.
+    ///
+    /// Smallest first, so index zero is 35 mm — the commonest frame, and the
+    /// one whose grain every published figure and every photographer's
+    /// intuition is calibrated against.
+    pub const ALL: [Format; 6] = [
+        Format::Mm35,
+        Format::Format645,
+        Format::Format6x6,
+        Format::Format6x7,
+        Format::Sheet4x5,
+        Format::Sheet8x10,
+    ];
+
+    /// The frame's width in micrometres — the *image* area, not the sheet.
+    pub fn width_um(self) -> f32 {
+        match self {
+            Format::Mm35 => 36_000.0,
+            // 6x4.5 and 6x6 share a 56 mm gate; only the other axis differs,
+            // and grain scales with the linear magnification of the axis being
+            // enlarged.
+            Format::Format645 | Format::Format6x6 => 56_000.0,
+            Format::Format6x7 => 70_000.0,
+            // The image area of a sheet, which is smaller than the nominal
+            // inches: a "4x5" exposes about 121 x 97 mm.
+            Format::Sheet4x5 => 121_000.0,
+            Format::Sheet8x10 => 248_000.0,
+        }
+    }
+
+    pub fn from_index(i: usize) -> Self {
+        Self::ALL.get(i).copied().unwrap_or(Format::Mm35)
+    }
+}
+
+/// A 35 mm frame's width, in micrometres. The default format.
 pub const FRAME_WIDTH_UM: f32 = 36_000.0;
 
 /// TRACES: FR-DEV-3f
@@ -245,6 +297,41 @@ mod tests {
         );
     }
 
+    #[test]
+    fn a_larger_format_is_less_grainy_at_the_same_output_size() {
+        // TRACES: FR-DEV-3f
+        // The point of the whole control, and a fact about photography rather
+        // than about this code: enlarge 35 mm and 4x5 to the same print and the
+        // sheet is visibly smoother, because each of its pixels averages far
+        // more crystals. Treating every frame as 35 mm made a large-format
+        // photograph as grainy as a small one.
+        let p = portra();
+        let out_px = 5472.0;
+        let small = Grain::for_pixel_size(&p, Format::Mm35.width_um() / out_px);
+        let large = Grain::for_pixel_size(&p, Format::Sheet4x5.width_um() / out_px);
+
+        let d = small.density_max[1] * 0.5;
+        let ratio = small.sigma(1, d) / large.sigma(1, d);
+        // Linear magnification is 121/36, so the crystal count per pixel goes
+        // as its square and sigma as its reciprocal: about 3.4x.
+        assert!(
+            (2.5..4.5).contains(&ratio),
+            "35mm is {ratio:.2}x grainier than 4x5, which is not the enlargement"
+        );
+    }
+
+    #[test]
+    fn the_formats_are_ordered_smallest_first() {
+        // Index zero has to be the neutral choice — 35 mm, which is what every
+        // published granularity figure is calibrated against.
+        let widths: Vec<f32> = Format::ALL.iter().map(|f| f.width_um()).collect();
+        assert_eq!(widths[0], FRAME_WIDTH_UM);
+        assert!(
+            widths.windows(2).all(|w| w[0] <= w[1]),
+            "formats are not ordered by size: {widths:?}"
+        );
+    }
+
     #[test]
     fn a_pixel_never_holds_less_than_one_grain() {
         // Past this the model describes a pixel smaller than a crystal, where

core/dr-film/src/lib.rs

@@ -38,6 +38,7 @@
 //! auditable. The sRGB reflectance basis is Mallett & Yuksel (2019).
 
 pub mod bake;
+pub mod boolean_grain;
 mod built_in;
 pub mod grain;
 pub mod profile;
@@ -45,7 +46,8 @@ pub mod spectrum;
 pub mod tables;
 
 pub use bake::{bake, Baked, Recipe};
-pub use grain::Grain;
+pub use boolean_grain::BooleanGrain;
+pub use grain::{Format, Grain};
 pub use profile::{Kind, Profile, Stage, Support};
 
 use built_in::BUILT_IN;

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

@@ -37,6 +37,23 @@ 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.
 ///
@@ -70,6 +87,13 @@ static DESCRIPTOR: OpDescriptor = OpDescriptor {
         // 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),
     ],
 };
 
@@ -129,6 +153,8 @@ 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>,
 }
 
@@ -173,6 +199,7 @@ impl Operation for FilmSim {
             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}"),
         }
     }
@@ -182,6 +209,7 @@ impl Operation for FilmSim {
             EXPOSURE => self.exposure,
             PRINT_EXPOSURE => self.print_exposure,
             PUSH => self.push,
+            FORMAT => self.format,
             _ => 0.0,
         }
     }