lamp-testbench/crates/tm30/src/data.rs

//! Embedded spectral tables (extracted by `tools/extract_data.py`).
//!
//! All tables are the *effective* arrays colour-science 0.4.7 works with on
//! the TM-30 path (after its Sprague re-alignment), stored as CSV with
//! shortest round-trip float formatting, parsed once on first use.

use std::sync::OnceLock;

/// First wavelength of the TM-30 working grid (nm).
pub(crate) const WL_START: usize = 380;
/// Last wavelength of the TM-30 working grid (nm).
pub(crate) const WL_END: usize = 780;
/// Number of samples of the 380..=780 nm / 1 nm working grid.
pub(crate) const N_WL: usize = WL_END - WL_START + 1;
/// First wavelength of colour's default spectral shape (Planckian table).
pub(crate) const WL_DEFAULT_START: usize = 360;
/// Number of samples of colour's default spectral shape 360..=780 nm / 1 nm.
pub(crate) const N_WL_DEFAULT: usize = WL_END - WL_DEFAULT_START + 1;
/// Number of CIE 2017 test colour samples.
pub(crate) const N_TCS: usize = 99;

const CMF_1931_360_780: &str = include_str!("../data/cmf_cie1931_2deg_360_780.csv");
const CMF_1931_380_780: &str = include_str!("../data/cmf_cie1931_2deg_380_780.csv");
const CMF_1964_380_780: &str = include_str!("../data/cmf_cie1964_10deg_380_780.csv");
const D_SERIES_380_780: &str = include_str!("../data/d_series_basis_380_780.csv");
const TCS_380_780: &str = include_str!("../data/tcs_cie2017_380_780.csv");
const CIECAM02_MATRICES: &str = include_str!("../data/ciecam02_matrices.csv");

pub(crate) type Mat3 = [[f64; 3]; 3];

pub(crate) struct Tables {
    /// CIE 1931 2° CMFs on 360..=780 nm (used for the Ohno Planckian table).
    pub cmf2_default: Vec<[f64; 3]>,
    /// CIE 1931 2° CMFs on 380..=780 nm (test CCT, reference normalisation).
    pub cmf2: Vec<[f64; 3]>,
    /// CIE 1964 10° CMFs on 380..=780 nm (all CIECAM02 colorimetry).
    pub cmf10: Vec<[f64; 3]>,
    /// CIE D-series basis functions S0, S1, S2 on 380..=780 nm.
    pub d_basis: Vec<[f64; 3]>,
    /// CIE 2017 TCS reflectances, `tcs[sample][wavelength]`.
    pub tcs: Vec<Vec<f64>>,
    /// CAT02 matrix.
    pub cat02: Mat3,
    /// `MATRIX_XYZ_TO_HPE @ inv(CAT02)` exactly as computed by numpy.
    pub hpe_cat02_inv: Mat3,
}

fn parse_rows(csv: &str, name: &str) -> Vec<Vec<f64>> {
    csv.lines()
        .skip(1)
        .filter(|l| !l.trim().is_empty())
        .map(|line| {
            line.split(',')
                .map(|field| {
                    field
                        .trim()
                        .parse::<f64>()
                        .unwrap_or_else(|e| panic!("{name}: bad number {field:?}: {e}"))
                })
                .collect()
        })
        .collect()
}

/// Parse a `wavelength,c1,c2,c3` table and check the wavelength column.
fn parse_tristimulus(csv: &str, name: &str, start: usize, n: usize) -> Vec<[f64; 3]> {
    let rows = parse_rows(csv, name);
    assert_eq!(rows.len(), n, "{name}: unexpected row count");
    rows.iter()
        .enumerate()
        .map(|(i, r)| {
            assert_eq!(r.len(), 4, "{name}: unexpected column count");
            assert_eq!(r[0], (start + i) as f64, "{name}: unexpected wavelength");
            [r[1], r[2], r[3]]
        })
        .collect()
}

fn parse_tcs() -> Vec<Vec<f64>> {
    let rows = parse_rows(TCS_380_780, "tcs");
    assert_eq!(rows.len(), N_WL, "tcs: unexpected row count");
    let mut tcs = vec![vec![0.0; N_WL]; N_TCS];
    for (w, r) in rows.iter().enumerate() {
        assert_eq!(r.len(), N_TCS + 1, "tcs: unexpected column count");
        assert_eq!(r[0], (WL_START + w) as f64, "tcs: unexpected wavelength");
        for (s, sample) in tcs.iter_mut().enumerate() {
            sample[w] = r[s + 1];
        }
    }
    tcs
}

fn parse_matrices() -> (Mat3, Mat3) {
    let mut cat02 = None;
    let mut hpe = None;
    for line in CIECAM02_MATRICES.lines().skip(1) {
        let mut fields = line.split(',');
        let Some(name) = fields.next() else { continue };
        let v: Vec<f64> = fields
            .map(|f| {
                f.trim()
                    .parse::<f64>()
                    .expect("ciecam02 matrices: bad number")
            })
            .collect();
        assert_eq!(v.len(), 9, "ciecam02 matrices: expected 9 values");
        let m = [[v[0], v[1], v[2]], [v[3], v[4], v[5]], [v[6], v[7], v[8]]];
        match name.trim() {
            "CAT02" => cat02 = Some(m),
            "HPE_CAT02_INV" => hpe = Some(m),
            other => panic!("ciecam02 matrices: unknown matrix {other}"),
        }
    }
    (
        cat02.expect("CAT02 missing"),
        hpe.expect("HPE_CAT02_INV missing"),
    )
}

pub(crate) fn tables() -> &'static Tables {
    static TABLES: OnceLock<Tables> = OnceLock::new();
    TABLES.get_or_init(|| {
        let (cat02, hpe_cat02_inv) = parse_matrices();
        Tables {
            cmf2_default: parse_tristimulus(
                CMF_1931_360_780,
                "cmf 1931 360-780",
                WL_DEFAULT_START,
                N_WL_DEFAULT,
            ),
            cmf2: parse_tristimulus(CMF_1931_380_780, "cmf 1931 380-780", WL_START, N_WL),
            cmf10: parse_tristimulus(CMF_1964_380_780, "cmf 1964 380-780", WL_START, N_WL),
            d_basis: parse_tristimulus(D_SERIES_380_780, "d-series basis", WL_START, N_WL),
            tcs: parse_tcs(),
            cat02,
            hpe_cat02_inv,
        }
    })
}