//! 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)
        );
    }
}