DarkRoom/core/dr-film/src/grain.rs

//! TRACES: FR-DEV-3f
//! Grain — the silver that did or did not develop.
//!
//! # Why this is not noise
//!
//! An emulsion is a suspension of silver-halide crystals. Light sensitises
//! some of them; development turns a sensitised crystal into an opaque grain,
//! all-or-nothing. So the density a patch of film reaches is a *count* of
//! developed grains, and a count of independent yes/no events has a variance
//! whether or not anybody wanted texture.
//!
//! Writing it down gives the look for free, rather than as a slider:
//!
//! ```text
//! p        = D / Dmax                    the chance one grain develops
//! N                                      grains in this pixel's patch of film
//! mean     = D
//! variance = D · (Dmax − u · D) / N
//! ```
//!
//! The variance peaks near the middle of the density range and vanishes at
//! both ends — clear film has nothing to develop, fully black film has nothing
//! left undeveloped. That is why grain lives in the midtones, and it is a
//! consequence here rather than a "midtone bias" control.
//!
//! # Why the zoom problem is not a problem
//!
//! `N` is grains *per pixel*, so it scales with how much film a pixel covers.
//! Zoomed out, each pixel averages more grains and the variance falls, which
//! is exactly what happens when you look at a print from further away. Nothing
//! has to be super-sampled and nothing has to be filtered: the model is already
//! a function of scale, and [`Grain::for_pixel_size`] is where the scale
//! enters.
//!
//! # Where crystal habit lives
//!
//! In [`crate::profile::Profile::rms_granularity`]. It measures exactly the
//! thing that differs between a traditional cubic emulsion and a tabular one:
//! for the same speed, tabular crystals present more area per unit of silver,
//! so fewer, flatter grains cover the frame and the film reads finer. A
//! stock's grain character therefore needs no new model and no new code — it
//! is one number in the profile.
//!
//! # Getting that number is the hard part
//!
//! It is **not** simply "what the datasheet says", and assuming so is a trap
//! worth naming:
//!
//! - **Kodak colour negatives publish Print Grain Index, not RMS.** PGI is a
//!   perceptual scale from viewer surveys — 25 is roughly the threshold of
//!   visibility and four units is one just-noticeable difference — and Kodak
//!   state it *cannot* be compared to RMS granularity. There is no published
//!   conversion, so a PGI figure cannot be dropped into this field.
//! - **Ilford publish no granularity figure at all.** Their technical sheets
//!   describe grain only in words ("fine grain", "finest grain"); the word
//!   granularity does not appear in them. Spectral sensitivity and the
//!   characteristic curves are there, but as graphs.
//! - RMS granularity proper is published mostly for black-and-white, reversal
//!   and motion-picture stocks.
//!
//! So the shipped defaults are plausible rather than measured, and every stock
//! currently carries the same ones — which means **grain does not yet tell one
//! stock from another**. Fixing that is per-stock data, not code, and for the
//! films whose figures are unpublished it needs a defensible estimate rather
//! than a number copied off a page.

use crate::profile::Profile;

/// The aperture RMS granularity is defined against: a 48 µm circle.
///
/// Fixed by the measurement, not by us. Every published granularity figure is
/// a standard deviation of density read through this aperture, so inverting it
/// for a particle area needs the same one.
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;

/// TRACES: FR-DEV-3f
/// The frame a photograph is being simulated on.
///
/// **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
/// What the shader needs to add grain to a density.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Grain {
    /// Grains in one pixel's patch of film, per layer. The whole scale
    /// dependence is in here.
    pub particles: [f32; 3],
    /// Each layer's maximum density — the ceiling `p = D/Dmax` is taken
    /// against.
    pub density_max: [f32; 3],
    /// How uniformly grains develop. Just below 1: a real emulsion's grains
    /// are not identical, and at exactly 1 the variance expression collapses
    /// at the top of the curve.
    pub uniformity: f32,
}

impl Grain {
    /// Derive the parameters for a stock at a given sampling scale.
    ///
    /// `pixel_size_um` is the film distance one rendered pixel covers. It is
    /// the only argument that changes with zoom, and it is what makes a
    /// preview show less grain than a 100% view without anything being
    /// filtered.
    pub fn for_pixel_size(profile: &Profile, pixel_size_um: f32) -> Self {
        let density_max = profile.max_density_per_layer();
        let density_min = profile.min_density_per_layer();
        let uniformity = profile.grain_uniformity;

        let pixel_area = (pixel_size_um * pixel_size_um).max(1e-6);
        let mut particles = [0.0f32; 3];
        for layer in 0..3 {
            // Invert the granularity definition for the area of one grain:
            //
            //     sigma48^2 = D_ref (Dmax - u D_ref) a / A48
            //
            // then count how many of those fit in a pixel. A stock with a
            // *lower* published granularity has smaller grains, so more of
            // them per pixel, so less variance -- which is the whole content
            // of the number.
            let sigma48 = profile.rms_granularity[layer] / 1000.0;
            let d_ref = RMS_REFERENCE_NET_DENSITY + density_min[layer];
            let denom = (d_ref * (density_max[layer] - uniformity * d_ref)).max(1e-6);
            let area = (sigma48 * sigma48 * RMS_APERTURE_AREA_UM2 / denom).max(1e-9);

            // At least one grain per pixel. Below that the model is describing
            // a pixel smaller than a single crystal, where "how many developed"
            // stops being a useful question and the variance would run away.
            particles[layer] = (pixel_area / area).max(1.0);
        }

        Self {
            particles,
            density_max,
            uniformity,
        }
    }

    /// The standard deviation of the developed density at density `d`.
    ///
    /// The shader evaluates this per pixel; it is here so the property can be
    /// asserted without a device.
    pub fn sigma(&self, layer: usize, d: f32) -> f32 {
        let dmax = self.density_max[layer];
        let d = d.clamp(0.0, dmax);
        let var = d * (dmax - self.uniformity * d) / self.particles[layer];
        var.max(0.0).sqrt()
    }
}

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

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

    fn trix() -> Profile {
        Profile::parse(include_str!("../profiles/kodak_trix.yaml")).unwrap()
    }

    #[test]
    fn grain_vanishes_at_both_ends_and_peaks_between() {
        // The property that makes this a model rather than a texture: clear
        // film has nothing developed to vary, and fully developed film has
        // nothing left to develop. Everything interesting is in the middle.
        let g = Grain::for_pixel_size(&portra(), 6.6);
        let dmax = g.density_max[1];
        let clear = g.sigma(1, 0.0);
        let mid = g.sigma(1, dmax * 0.5);
        let full = g.sigma(1, dmax);
        assert!(clear < mid && full < mid, "{clear} {mid} {full}");
        assert!(clear < 1e-6, "clear film is not grainless: {clear}");
    }

    #[test]
    fn a_smaller_pixel_sees_more_grain() {
        // The scale dependence, and the reason nothing needs super-sampling: a
        // pixel covering less film averages fewer grains, so it is noisier. A
        // model that got this backwards would show *more* grain in a zoomed-out
        // preview than in the export.
        let p = portra();
        let close = Grain::for_pixel_size(&p, 3.0);
        let far = Grain::for_pixel_size(&p, 12.0);
        let d = close.density_max[1] * 0.5;
        assert!(
            close.sigma(1, d) > far.sigma(1, d),
            "close {} far {}",
            close.sigma(1, d),
            far.sigma(1, d)
        );
    }

    #[test]
    fn halving_the_pixel_halves_the_grain() {
        // Not merely monotone: variance goes as 1/N and N goes as area, so
        // sigma goes as 1/pixel_size. Asserted because the exponent is what
        // decides whether a preview at half scale looks like the export.
        let p = portra();
        let a = Grain::for_pixel_size(&p, 4.0);
        let b = Grain::for_pixel_size(&p, 8.0);
        let d = a.density_max[1] * 0.5;
        let ratio = a.sigma(1, d) / b.sigma(1, d);
        assert!(
            (ratio - 2.0).abs() < 0.05,
            "sigma ratio {ratio}, expected 2"
        );
    }

    #[test]
    fn a_coarser_stock_is_grainier_at_the_same_scale() {
        // What the granularity number is for, and what will carry crystal
        // habit when a tabular stock arrives: the same pixel, the same
        // density, a coarser emulsion, more noise.
        let mut fine = portra();
        let mut coarse = portra();
        fine.rms_granularity = [4.0; 3];
        coarse.rms_granularity = [16.0; 3];

        let d = fine.max_density_per_layer()[1] * 0.5;
        let f = Grain::for_pixel_size(&fine, 6.6).sigma(1, d);
        let c = Grain::for_pixel_size(&coarse, 6.6).sigma(1, d);
        assert!(c > f * 2.0, "fine {f} coarse {c}");
    }

    #[test]
    fn a_monochrome_stock_grains_every_layer_alike() {
        // Its three layers are one emulsion spread three ways, so a difference
        // between them would be grain the film does not have — and it would
        // show as colour speckle on a black and white photograph, which is the
        // most obvious way this could look wrong.
        let g = Grain::for_pixel_size(&trix(), 6.6);
        let d = g.density_max[0] * 0.5;
        let s: Vec<f32> = (0..3).map(|l| g.sigma(l, d)).collect();
        assert!(
            (s[0] - s[1]).abs() < 1e-5 && (s[1] - s[2]).abs() < 1e-5,
            "monochrome grain differs by layer: {s:?}"
        );
    }

    #[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
        // the variance runs away and the picture fills with salt and pepper.
        let g = Grain::for_pixel_size(&portra(), 0.001);
        assert!(g.particles.iter().all(|n| *n >= 1.0), "{:?}", g.particles);
    }
}