Count the silver instead of adding noise

An emulsion is a suspension of crystals. Light sensitises some; development
turns a sensitised one opaque, all or nothing. So a patch of film's density
is a *count* of developed grains, and a count of independent yes/no events
has a variance whether or not anyone wanted texture:

    mean     = D
    variance = D * (Dmax - u * D) / N

That expression is the whole feature. It peaks in the middle of the density
range and vanishes at both ends -- clear film has nothing developed to vary,
black film has nothing left to develop -- so grain lives in the midtones as a
consequence rather than as a "midtone bias" slider.

I was wrong earlier that this needs the detail stage. Nothing in it reads a
neighbouring pixel; the only reason to move it was that grain must be fixed in
film space rather than screen space, and that solves itself: N is grains *per
pixel*, so it scales with the film a pixel covers. Zoom out, each pixel
averages more grains, less variance -- correct, with nothing super-sampled and
nothing filtered. It stays in the fused pass.

Grain goes on the density and *before* the dye, which is the physical order
and not cosmetic. Perturbing the finished colour -- what an effect does --
tints highlights wrong, because that noise never passes through the dye.

Crystal habit lives in `rms_granularity`, the number every datasheet
publishes, now a profile field. It measures exactly what differs between a
cubic emulsion and a tabular one: at equal speed, tabular crystals present
more area per unit silver, so the film reads finer. Delta 100 is quoted near 9
where HP5 is near 12, and that gap *is* the habit. Adding a stock whose grain
is its whole reputation is therefore editing one line, not writing a model.

Three things this cost, all of them worth writing down:

  - The default granularity is a colour negative's, blue coarsest. Applied to
    Tri-X it put *colour* speckle on a black and white photograph. Monochrome
    stocks collapse it at parse, where every other per-layer table is already
    replicated from the one measured channel.
  - Helpers cannot read uniforms. The composer prefixes a uniform with its
    operation's id and rewrites references inside a fragment body only;
    helpers are shared and deduplicated, so a bare `gn0` names nothing.
    `film_lut` already took its size as an argument for this reason, and now
    says so.
  - The end-to-end test compares the shader against the CPU model, and grain
    is stochastic, so that comparison now runs with grain off. Which means a
    grain that never left the CPU would look exactly like a passing suite --
    hence a second test that grain off is bit-identical, one grain per pixel
    moves it, and ten thousand move it less.

Not here, deliberately: no grain slider. The parameters are physical and
`rms_granularity` is the honest place to scale one from, but its range wants
choosing rather than guessing. Nor a film format -- 35 mm is assumed, and
medium format at the same stock is far less grainy per unit of picture.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

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

@@ -93,6 +93,13 @@ pub struct FilmTables {
     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 {
@@ -207,6 +214,22 @@ impl Operation for FilmSim {
                 });
             }
         }
+        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,
@@ -265,10 +288,23 @@ let log_exposure = log10(max(exposure, vec3<f32>(0.0)) + 1e-10);
 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(density / density_max, vec3<f32>(0.0), vec3<f32>(1.0)), lut_size);"
+c = film_lut(clamp(grained / density_max, vec3<f32>(0.0), vec3<f32>(1.0)), lut_size);"
             .into()
     }
 
@@ -277,7 +313,85 @@ c = film_lut(clamp(density / density_max, vec3<f32>(0.0), vec3<f32>(1.0)), lut_s
     }
 }
 
-static HELPERS: [crate::operation::Helper; 3] = [
+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: "\
@@ -349,6 +463,9 @@ mod tests {
             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,
         }
     }