//! TRACES: FR-DEV-3j | FR-DEV-2
//! The view transform — the one stage that maps scene-linear colour to a
//! display range (D19, ARCH §6.14).
//!
//! # What it is
//!
//! A log-logistic sigmoid, per channel:
//!
//! ```text
//! f(x) = w · r / (1 + r),   r = (x / k)^n
//! ```
//!
//! `n` is the contrast — the slope in log-log terms, before the shoulder
//! bends it. `k` and `w` are solved from two conditions rather than set:
//! scene middle grey lands on display middle grey, and the scene white the
//! photographer chose lands on display white. So the curve has a toe, a
//! midtone slope and a shoulder that approaches `w` — a hair above 1.0 —
//! without ever reaching it. Everything the shoulder has not reached by the
//! white point is clipped by the output transform, which is the last moment
//! and the only place a clip belongs.
//!
//! # Why per channel, and why the middle channel is put back
//!
//! Per channel is what makes a bright saturated colour desaturate as it
//! approaches white — a blown sky rolls toward white rather than toward a
//! saturated corner of the gamut, which is what film and every camera JPEG
//! do. It also bends hue: the three channels sit at different places on the
//! curve, so their ratios change, and an orange flame drifts toward yellow.
//! So after the curve the middle channel is moved back to where it sat
//! *between the other two* before it — the same fraction of the way from the
//! smallest to the largest. The smallest and largest keep what the curve gave
//! them, which keeps the desaturation; the hue, which is decided by that
//! fraction, survives. It is the "preserve hue" step of darktable's sigmoid,
//! at full strength.
//!
//! # Why these defaults
//!
//! [`SCENE_GREY`] is where the retired default base curve put middle grey
//! (FR-DEV-3e): linear sensor data from a correctly exposed frame has it
//! near 13% of saturation, and a camera JPEG shows it at 18%. The contrast
//! and white defaults were chosen against that same retired curve: at 1.4 and
//! 4 stops the midtones stay within a quarter of a stop of it between scene
//! 0.03 and 1.0, while a highlight a stop past sensor saturation still rolls
//! into white rather than stopping dead at it. The upper midtones come out a
//! little darker than the curve had them, which is the price of that
//! headroom and what the white slider is for.

/// Scene-linear middle grey: where the retired default curve placed it.
pub const SCENE_GREY: f32 = 0.13;

/// Display-linear middle grey — what a camera JPEG shows a grey card as.
pub const DISPLAY_GREY: f32 = 0.18;

/// The default contrast: the sigmoid's log-log slope `n`, and for the DNG
/// reference curve a power of `DEFAULT_CONTRAST / REFERENCE_CONTRAST` about grey.
///
/// Fitted, not chosen: on Lightroom 6 exports whose look settings were neutral,
/// the DNG reference curve matched the exports best with the input bent by about
/// 1.08 (held-out MSE 224 at 1.4, ~150 at 1.5). The sigmoid, which is no longer
/// the default, fitted best near 1.7 and is better at 1.5 than at 1.4.
pub const DEFAULT_CONTRAST: f32 = 1.5;

/// The contrast at which each curve is its own reference: the sigmoid's match
/// to the retired base curve (see the module note), and the reference table
/// untouched.
pub const REFERENCE_CONTRAST: f32 = 1.4;

/// The default white point, in stops above [`SCENE_GREY`].
pub const DEFAULT_WHITE: f32 = 4.0;

/// The contrast range a photographer is offered.
pub const CONTRAST_RANGE: (f32, f32) = (1.0, 3.0);

/// The white point range, in stops above middle grey.
///
/// The floor is not taste. The two conditions `k` and `w` are solved from
/// have a solution only while `2^(white · n)` exceeds `1 / DISPLAY_GREY`,
/// and at the lowest contrast that needs `white` above about 2.47 stops.
pub const WHITE_RANGE: (f32, f32) = (2.5, 10.0);

/// The curve's three numbers, solved from the photographer's two.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Sigmoid {
    /// Contrast: the exponent.
    pub n: f32,
    /// `1 / k`, so the shader multiplies rather than divides.
    pub inv_k: f32,
    /// The asymptote the shoulder approaches, a little above 1.0.
    pub w: f32,
}

impl Sigmoid {
    /// Solve the curve for a contrast and a white point in stops.
    ///
    /// Out-of-range inputs are clamped to [`CONTRAST_RANGE`] and
    /// [`WHITE_RANGE`] rather than trusted: they arrive from a sidecar, which
    /// may have been written by a build with other limits, and outside them
    /// the solution below divides by something that is no longer positive.
    ///
    /// With `r_g` the value of `r` at scene grey and `q = 2^(white · n)`, the
    /// two conditions `f(grey) = display grey` and `f(grey · 2^white) = 1`
    /// are `w·r_g/(1+r_g) = g` and `w·q·r_g/(1+q·r_g) = 1`. Dividing one by
    /// the other eliminates `w` and leaves `r_g = (g·q − 1) / (q·(1 − g))`.
    pub fn new(contrast: f32, white: f32) -> Self {
        let n = contrast.clamp(CONTRAST_RANGE.0, CONTRAST_RANGE.1) as f64;
        let white = white.clamp(WHITE_RANGE.0, WHITE_RANGE.1) as f64;
        let g = f64::from(DISPLAY_GREY);
        let q = (white * n).exp2();
        let r_grey = (g * q - 1.0) / (q * (1.0 - g));
        let w = g * (1.0 + r_grey) / r_grey;
        // r = (x / k)^n, and r at scene grey is r_grey, so
        // k = grey / r_grey^(1/n).
        let k = f64::from(SCENE_GREY) / r_grey.powf(1.0 / n);
        Self {
            n: n as f32,
            inv_k: (1.0 / k) as f32,
            w: w as f32,
        }
    }

    /// The default curve.
    pub fn default_curve() -> Self {
        Self::new(DEFAULT_CONTRAST, DEFAULT_WHITE)
    }

    /// One channel through the curve. The CPU reference the shader is
    /// tested against.
    pub fn channel(&self, x: f32) -> f32 {
        let r = (x.max(0.0) * self.inv_k).powf(self.n);
        self.w * r / (1.0 + r)
    }

    /// A colour through the curve, with the middle channel put back between
    /// the other two. See the module documentation.
    pub fn apply(&self, c: [f32; 3]) -> [f32; 3] {
        let x = c.map(|v| v.max(0.0));
        let y = x.map(|v| self.channel(v));
        let lo = x[0].min(x[1]).min(x[2]);
        let hi = x[0].max(x[1]).max(x[2]);
        if hi - lo <= 1e-9 {
            return y;
        }
        let (y_lo, y_hi) = (self.channel(lo), self.channel(hi));
        x.map(|v| y_lo + (y_hi - y_lo) * (v - lo) / (hi - lo))
    }
}

/// The WGSL twin of [`Sigmoid::apply`], as a helper function.
pub const VIEW_SIGMOID_WGSL: &str = "\
// The view transform (FR-DEV-3j): a log-logistic sigmoid per channel, then the
// middle channel put back between the other two so that the hue survives the
// shoulder. See `dr_pipeline::view` for the derivation and the defaults.
fn view_sigmoid(c: vec3<f32>, n: f32, inv_k: f32, w: f32) -> vec3<f32> {
    // Negative components are colours outside the working primaries. They are
    // floored here, at the last stage, which is the one place a gamut clip
    // belongs.
    let x = max(c, vec3<f32>(0.0));
    let lo = min(x.r, min(x.g, x.b));
    let hi = max(x.r, max(x.g, x.b));
    let r_lo = pow(lo * inv_k, n);
    let r_hi = pow(hi * inv_k, n);
    let y_lo = w * r_lo / (1.0 + r_lo);
    let y_hi = w * r_hi / (1.0 + r_hi);
    // Each channel's place between the smallest and the largest. A neutral
    // has no spread, and every channel then takes the one value there is.
    let spread = hi - lo;
    let t = select((x - vec3<f32>(lo)) / max(spread, 1e-9), vec3<f32>(0.0), spread <= 1e-9);
    return vec3<f32>(y_lo) + (y_hi - y_lo) * t;
}
";

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

    /// The retired default base curve, for the acceptance comparison: five
    /// points through the unit square, sampled here by straight lines between
    /// them in log-log terms — close enough to the monotone spline that drew
    /// it for a tolerance measured in quarters of a stop. Below its second
    /// point it was a straight line from the origin.
    fn retired_default(x: f32) -> f32 {
        const P: [(f32, f32); 4] = [(0.04, 0.043), (0.13, 0.175), (0.45, 0.690), (1.0, 1.0)];
        if x < 0.04 {
            return x * (0.043 / 0.04);
        }
        let x = x.min(1.0);
        let i = P.windows(2).position(|w| x <= w[1].0).unwrap_or(2);
        let ((x0, y0), (x1, y1)) = (P[i], P[i + 1]);
        let t = (x.ln() - x0.ln()) / (x1.ln() - x0.ln());
        (y0.ln() + t * (y1.ln() - y0.ln())).exp()
    }

    #[test]
    fn middle_grey_lands_on_display_grey() {
        // TRACES: FR-DEV-3j
        for (contrast, white) in [(1.0, 3.0), (1.4, 4.0), (2.5, 8.0), (3.0, 10.0)] {
            let s = Sigmoid::new(contrast, white);
            let got = s.channel(SCENE_GREY);
            assert!(
                (got - DISPLAY_GREY).abs() < 0.01,
                "contrast {contrast}, white {white}: grey went to {got}"
            );
        }
    }

    #[test]
    fn the_white_point_reaches_display_white() {
        // TRACES: FR-DEV-3j
        // The whole meaning of the slider: the scene value it names is where
        // the picture reaches white, and not before.
        for (contrast, white) in [(1.0, 3.0), (1.4, 4.0), (2.5, 8.0)] {
            let s = Sigmoid::new(contrast, white);
            let at = SCENE_GREY * white.exp2();
            assert!((s.channel(at) - 1.0).abs() < 1e-4, "{}", s.channel(at));
            assert!(s.channel(at * 0.9) < 1.0);
            assert!(s.w > 1.0, "the shoulder must approach a value above white");
        }
    }

    #[test]
    fn the_curve_is_monotone_and_keeps_going_past_one() {
        // TRACES: FR-DEV-3j | FR-DEV-2
        // What the base curve got wrong: it was flat past 1.0, so every
        // recovered highlight left it as the same number.
        let s = Sigmoid::default_curve();
        let mut last = -1.0;
        for i in 0..=2000 {
            let x = i as f32 * 0.004;
            let y = s.channel(x);
            assert!(y > last || (x == 0.0 && y == 0.0), "not increasing at {x}");
            last = y;
        }
        assert!(s.channel(2.0) > s.channel(1.0));
    }

    #[test]
    fn the_reference_stays_close_to_the_retired_curve() {
        // TRACES: FR-DEV-3j | FR-DEV-3e
        // D19's promise, now kept by the sigmoid at its reference contrast
        // rather than by the default (D21 moved the default to the DNG reference
        // curve): the midtones do not move by more than a third of a stop
        // from the retired base curve.
        let s = Sigmoid::new(REFERENCE_CONTRAST, DEFAULT_WHITE);
        let mut x = 0.03_f32;
        while x <= 1.0 {
            let ev = (s.channel(x) / retired_default(x)).log2();
            assert!(ev.abs() < 0.3, "at scene {x} the default moved {ev:+.2} EV");
            x *= 1.1;
        }
    }

    #[test]
    fn a_neutral_stays_neutral() {
        // TRACES: FR-DEV-3j
        let s = Sigmoid::default_curve();
        for v in [0.0, 0.01, 0.13, 1.0, 7.0] {
            let [r, g, b] = s.apply([v, v, v]);
            assert_eq!(r, g);
            assert_eq!(g, b);
        }
    }

    #[test]
    fn the_hue_survives_the_shoulder() {
        // TRACES: FR-DEV-3j
        // The middle channel's place between the other two is what decides
        // the hue. Without the correction an orange at the shoulder drifts
        // toward yellow as the red channel saturates first.
        let s = Sigmoid::default_curve();
        let orange = [2.0, 0.8, 0.1];
        let out = s.apply(orange);
        let before = (orange[1] - orange[2]) / (orange[0] - orange[2]);
        let after = (out[1] - out[2]) / (out[0] - out[2]);
        assert!((before - after).abs() < 1e-5, "{before} became {after}");
        // And the extremes keep what the curve gave them — the desaturation
        // toward white is the point of working per channel.
        assert!((out[0] - s.channel(2.0)).abs() < 1e-6);
        assert!((out[2] - s.channel(0.1)).abs() < 1e-6);
    }

    #[test]
    fn out_of_range_settings_are_clamped_not_trusted() {
        // A sidecar from another build may carry anything, and outside the
        // range the solution divides by a value that is no longer positive.
        let s = Sigmoid::new(0.0, 0.0);
        assert!(s.n.is_finite() && s.inv_k.is_finite() && s.w.is_finite());
        assert_eq!(s, Sigmoid::new(CONTRAST_RANGE.0, WHITE_RANGE.0));
    }
}