//! The fixed colour science: spectra in, tristimulus out.
//!
//! Everything here is per-installation rather than per-stock — one observer,
//! a handful of illuminants, one basis — and all of it is small enough to
//! compile in. See [`crate::tables`] for the numbers themselves.

use crate::tables::{ILLUMINANT_D50, ILLUMINANT_D55, ILLUMINANT_D65, OBSERVER, SPECTRUM};

/// A spectral distribution on the crate's fixed 380–780 nm, 5 nm grid.
pub type Spectrum = [f32; SPECTRUM];

/// Linear sRGB primaries from CIE XYZ, for a D65 white.
///
/// The adaptation to whatever white the picture is actually being viewed under
/// happens before this is applied — see [`Viewing::to_srgb`].
const XYZ_TO_SRGB: [[f32; 3]; 3] = [
    [3.2404542, -1.5371385, -0.4985314],
    [-0.969266, 1.8760108, 0.041556],
    [0.0556434, -0.2040259, 1.0572252],
];

/// sRGB's own white, which is what [`XYZ_TO_SRGB`] expects to be handed.
const D65_WHITE: [f32; 3] = [0.9504559, 1.0, 1.0890578];

/// Bradford cone response, and its inverse.
///
/// A von Kries adaptation done in XYZ — dividing each channel by the white —
/// is the obvious thing and it is wrong: XYZ axes are not cone responses, so
/// the result drifts in hue. Bradford is the transform that makes "the same
/// colour, seen under a different light" mean what a viewer means by it, and
/// it matters here because a print is specified under D50 while sRGB is D65.
const BRADFORD: [[f32; 3]; 3] = [
    [0.8951, 0.2664, -0.1614],
    [-0.7502, 1.7135, 0.0367],
    [0.0389, -0.0685, 1.0296],
];

const BRADFORD_INV: [[f32; 3]; 3] = [
    [0.9869929, -0.1470543, 0.1599627],
    [0.4323053, 0.5183603, 0.0492912],
    [-0.0085287, 0.0400428, 0.9684867],
];

/// Named illuminants a profile may ask for.
///
/// The daylight ones are tabulated because the CIE defines them that way. The
/// tungsten ones are computed from Planck's law, which costs no data at all,
/// and the approximation is unusually safe here: an enlarger's lamp colour is
/// cancelled by the filtration solved in [`crate::bake::print_balance`], the
/// same way a darkroom worker dials it out on the colour head.
pub fn illuminant(name: &str) -> Spectrum {
    match name {
        "D50" => ILLUMINANT_D50,
        "D55" => ILLUMINANT_D55,
        "D65" => ILLUMINANT_D65,
        // The tungsten-halogen enlarger source, with and without the heat
        // filter that a real head carries. Both land on the same blackbody:
        // the filter's effect is a colour shift, and a colour shift ahead of
        // the balance step is by construction invisible.
        "TH-KG3" | "TH-KG3-L" | "T" => blackbody(3400.0),
        // A cinema projector's xenon short-arc lamp, which is what a release
        // print is *looked at* under — Kodak 2383 and 2393 name it.
        //
        // Approximated by D55, and named here rather than left to fall through
        // to the unknown branch, which logged a warning on every bake of a
        // Vision3 stock and implied something was broken. A xenon arc sits near
        // 6000 K, close to daylight and nothing like the tungsten above, so
        // daylight is the right family to stand in for it. It is not a
        // blackbody — the arc has line structure a Planckian curve has no way
        // to express — and the residue of that is small here because the
        // viewing step adapts the white point out either way.
        "K75P" => ILLUMINANT_D55,
        other => {
            if let Some(kelvin) = other.strip_prefix("BB").and_then(|k| k.parse::<f32>().ok()) {
                return blackbody(kelvin);
            }
            log::warn!("unknown illuminant {other:?}, falling back to D55");
            ILLUMINANT_D55
        }
    }
}

/// Planck's law on the grid, normalised to unit mean.
///
/// Normalised because only the shape matters: absolute level is set by the
/// exposure controls, and leaving it in would make the choice of units a
/// visible parameter.
pub fn blackbody(kelvin: f32) -> Spectrum {
    const H: f64 = 6.626_070_15e-34;
    const C: f64 = 2.997_924_58e8;
    const KB: f64 = 1.380_649e-23;

    let mut out = [0.0f32; SPECTRUM];
    let mut total = 0.0f64;
    for (i, slot) in out.iter_mut().enumerate() {
        let lambda =
            f64::from(crate::tables::LAMBDA_MIN + crate::tables::LAMBDA_STEP * i as f32) * 1e-9;
        let radiance = (2.0 * H * C * C)
            / (lambda.powi(5) * ((H * C / (lambda * KB * f64::from(kelvin))).exp() - 1.0));
        *slot = radiance as f32;
        total += radiance;
    }
    let mean = (total / SPECTRUM as f64) as f32;
    for slot in &mut out {
        *slot /= mean;
    }
    out
}

/// How a developed image is looked at: an illuminant, and the adaptation it
/// implies.
///
/// Built once per bake rather than per sample, because the white point and the
/// adaptation matrix depend only on the illuminant and computing them inside
/// the LUT loop would be the same work 32 768 times.
pub struct Viewing {
    illuminant: Spectrum,
    /// The illuminant's own Y, which normalises the integral so that a clear
    /// frame comes out at exactly 1.0 rather than at whatever the tabulated
    /// units happen to give.
    normalisation: f32,
    /// XYZ under this illuminant to linear sRGB, adaptation folded in.
    matrix: [[f32; 3]; 3],
}

impl Viewing {
    pub fn new(illuminant_name: &str) -> Self {
        let illuminant = illuminant(illuminant_name);

        let mut white = [0.0f32; 3];
        let mut normalisation = 0.0f32;
        for i in 0..SPECTRUM {
            normalisation += illuminant[i] * OBSERVER[i][1];
            for c in 0..3 {
                white[c] += illuminant[i] * OBSERVER[i][c];
            }
        }
        for c in &mut white {
            *c /= normalisation;
        }

        Self {
            illuminant,
            normalisation,
            matrix: mat3_mul(XYZ_TO_SRGB, bradford_adaptation(white, D65_WHITE)),
        }
    }

    /// Integrate a transmittance (or reflectance) against the illuminant and
    /// the observer, and convert to linear sRGB.
    pub fn to_srgb(&self, transmittance: &Spectrum) -> [f32; 3] {
        let mut xyz = [0.0f32; 3];
        for i in 0..SPECTRUM {
            let light = transmittance[i] * self.illuminant[i];
            for c in 0..3 {
                xyz[c] += light * OBSERVER[i][c];
            }
        }
        for c in &mut xyz {
            *c /= self.normalisation;
        }
        mat3_apply(self.matrix, xyz)
    }
}

/// The Bradford transform taking `from` white to `to` white.
fn bradford_adaptation(from: [f32; 3], to: [f32; 3]) -> [[f32; 3]; 3] {
    let src = mat3_apply(BRADFORD, from);
    let dst = mat3_apply(BRADFORD, to);
    let scale = [
        [dst[0] / src[0], 0.0, 0.0],
        [0.0, dst[1] / src[1], 0.0],
        [0.0, 0.0, dst[2] / src[2]],
    ];
    mat3_mul(BRADFORD_INV, mat3_mul(scale, BRADFORD))
}

pub(crate) fn mat3_apply(m: [[f32; 3]; 3], v: [f32; 3]) -> [f32; 3] {
    [
        m[0][0] * v[0] + m[0][1] * v[1] + m[0][2] * v[2],
        m[1][0] * v[0] + m[1][1] * v[1] + m[1][2] * v[2],
        m[2][0] * v[0] + m[2][1] * v[1] + m[2][2] * v[2],
    ]
}

pub(crate) fn mat3_mul(a: [[f32; 3]; 3], b: [[f32; 3]; 3]) -> [[f32; 3]; 3] {
    let mut out = [[0.0f32; 3]; 3];
    for (r, row) in out.iter_mut().enumerate() {
        for (c, slot) in row.iter_mut().enumerate() {
            *slot = (0..3).map(|k| a[r][k] * b[k][c]).sum();
        }
    }
    out
}

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

    #[test]
    fn a_clear_frame_is_white_under_every_viewing_illuminant() {
        // The property the whole viewing step exists to hold: an unexposed
        // transparency on a light table is white, whatever the light table is.
        // Normalising on luminance alone would leave the illuminant's own
        // colour in the result, and a D50 print would render warm.
        for name in ["D50", "D55", "D65", "TH-KG3"] {
            let rgb = Viewing::new(name).to_srgb(&[1.0; SPECTRUM]);
            for c in rgb {
                assert!(
                    (c - 1.0).abs() < 0.01,
                    "{name}: clear frame gave {rgb:?}, which is not neutral white"
                );
            }
        }
    }

    #[test]
    fn adapting_a_white_to_itself_changes_nothing() {
        let m = bradford_adaptation(D65_WHITE, D65_WHITE);
        for r in 0..3 {
            for c in 0..3 {
                let expected = if r == c { 1.0 } else { 0.0 };
                assert!(
                    (m[r][c] - expected).abs() < 1e-5,
                    "{m:?} is not the identity"
                );
            }
        }
    }

    #[test]
    fn every_illuminant_a_shipped_profile_names_is_known() {
        // The warning that found this said "unknown illuminant K75P, falling
        // back to D55" on every bake of a Vision3 stock. The fallback was
        // reasonable and the silence was not: a profile naming a light nobody
        // implemented should be a failing test rather than a line in a log.
        for name in ["D50", "D55", "D65", "TH-KG3", "TH-KG3-L", "T", "K75P"] {
            let rgb = Viewing::new(name).to_srgb(&[1.0; SPECTRUM]);
            assert!(
                rgb.iter().all(|c| (c - 1.0).abs() < 0.01),
                "{name} is not a working illuminant: clear frame gave {rgb:?}"
            );
        }
    }

    #[test]
    fn a_hotter_blackbody_is_bluer() {
        // Cheap, but it is the one property that catches Planck's law written
        // with a sign or a reciprocal wrong, which otherwise produces a
        // plausible-looking curve pointing the wrong way.
        let cool = blackbody(2800.0);
        let hot = blackbody(9000.0);
        let blue = 10; // 430nm
        let red = 70; // 730nm
        assert!(hot[blue] / hot[red] > cool[blue] / cool[red]);
    }

    #[test]
    fn the_observer_integrates_to_a_plausible_white() {
        // Guards the generated table against a transposed or misaligned paste:
        // D65 through the 1931 observer must land on sRGB's own white.
        let v = Viewing::new("D65");
        let rgb = v.to_srgb(&[1.0; SPECTRUM]);
        assert!(rgb.iter().all(|c| (c - 1.0).abs() < 0.01), "{rgb:?}");
    }
}