ph-color-bake 0.1.0

//! Ottosson Oklab matrices and cube-root LUT for host baking.

use crate::lut::bake_lut;
use crate::matrix::to_q428;

/// Linear sRGB → LMS (Ottosson).
#[must_use]
pub fn m1() -> [[f64; 3]; 3] {
    [
        [0.4122214708, 0.5363325363, 0.0514459929],
        [0.2119034982, 0.6806995451, 0.1073969566],
        [0.0883024619, 0.2817188376, 0.6299787005],
    ]
}

/// LMS' → Lab (Ottosson).
#[must_use]
pub fn m2() -> [[f64; 3]; 3] {
    [
        [0.2104542553, 0.7936177850, -0.0040720468],
        [1.9779984951, -2.4285922050, 0.4505937099],
        [0.0259040371, 0.7827717662, -0.8086757660],
    ]
}

/// LMS → linear sRGB (Ottosson).
#[must_use]
pub fn m1_inv() -> [[f64; 3]; 3] {
    [
        [4.0767416621, -3.3077115913, 0.2309699292],
        [-1.2684380046, 2.6097574011, -0.3413193965],
        [-0.0041960863, -0.7034186147, 1.7076147010],
    ]
}

/// Lab → LMS' (Ottosson).
#[must_use]
pub fn m2_inv() -> [[f64; 3]; 3] {
    [
        [1.0, 0.3963377774, 0.2158037573],
        [1.0, -0.1055613458, -0.0638541728],
        [1.0, -0.0894841775, -1.2914855480],
    ]
}

/// Map a real 3×3 to Q4.28.
#[must_use]
pub fn to_q428_mat(m: [[f64; 3]; 3]) -> [[i32; 3]; 3] {
    [
        [to_q428(m[0][0]), to_q428(m[0][1]), to_q428(m[0][2])],
        [to_q428(m[1][0]), to_q428(m[1][1]), to_q428(m[1][2])],
        [to_q428(m[2][0]), to_q428(m[2][1]), to_q428(m[2][2])],
    ]
}

/// Affine addend so `M2_inv` can consume offset `a`/`b` stored as `value + 0.5`.
#[must_use]
pub fn m2_inv_offset_q428() -> [i32; 3] {
    let m = m2_inv();
    [
        to_q428(m[0][1] * -0.5 + m[0][2] * -0.5),
        to_q428(m[1][1] * -0.5 + m[1][2] * -0.5),
        to_q428(m[2][1] * -0.5 + m[2][2] * -0.5),
    ]
}

/// Cube root on the unit interval.
#[must_use]
pub fn cbrt(x: f64) -> f64 {
    x.clamp(0.0, 1.0).cbrt()
}

/// Bake the 257-knot cube-root table.
#[must_use]
pub fn bake_cbrt257() -> crate::lut::BakedLut<257> {
    bake_lut::<257>(cbrt)
}

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

    #[test]
    fn shipped_tables_match_bake() {
        use ::ph_color::Q4_28;

        let cbrt257 = bake_cbrt257();
        assert_eq!(cbrt257.max_err_lsb, ::ph_color::CBRT_MAX_ERR_LSB);
        assert_eq!(::ph_color::CBRT.max_err_lsb(), cbrt257.max_err_lsb);
        assert_eq!(
            Q4_28::mat3_from_raw(to_q428_mat(m1())),
            ::ph_color::OKLAB_M1
        );
        assert_eq!(
            Q4_28::mat3_from_raw(to_q428_mat(m2())),
            ::ph_color::OKLAB_M2
        );
        assert_eq!(
            Q4_28::mat3_from_raw(to_q428_mat(m1_inv())),
            ::ph_color::OKLAB_M1_INV
        );
        assert_eq!(
            Q4_28::mat3_from_raw(to_q428_mat(m2_inv())),
            ::ph_color::OKLAB_M2_INV
        );
        assert_eq!(
            Q4_28::row_from_raw(m2_inv_offset_q428()),
            ::ph_color::OKLAB_M2_INV_OFFSET
        );
        for x in 0..=u16::MAX {
            assert_eq!(
                ::ph_color::CBRT
                    .lookup(::ph_color::Q0_16::from_raw(x))
                    .to_raw(),
                crate::lut::lookup(&cbrt257.knots, x)
            );
        }
    }

    fn unit_to_u16(x: f64) -> u16 {
        let x = x.clamp(0.0, 1.0);
        let scaled = x * 65535.0;
        let trunc = scaled.trunc();
        let frac = scaled - trunc;
        let rounded = if frac >= 0.5 { trunc + 1.0 } else { trunc };
        if rounded >= 65535.0 {
            65535
        } else {
            rounded as u16
        }
    }

    fn mul(mat: [[f64; 3]; 3], v: [f64; 3]) -> [f64; 3] {
        let [r0, r1, r2] = mat;
        let [x, y, z] = v;
        [
            r0[0].mul_add(x, r0[1].mul_add(y, r0[2] * z)),
            r1[0].mul_add(x, r1[1].mul_add(y, r1[2] * z)),
            r2[0].mul_add(x, r2[1].mul_add(y, r2[2] * z)),
        ]
    }

    fn ottosson_srgb_to_oklab(ch: [u16; 3]) -> [u16; 3] {
        let [r, g, b] = ch;
        let v = [
            f64::from(r) / 65535.0,
            f64::from(g) / 65535.0,
            f64::from(b) / 65535.0,
        ];
        let lms = mul(m1(), v);
        let [l, m, s] = lms;
        let lp = [cbrt(l), cbrt(m), cbrt(s)];
        let [l_star, a, b_star] = mul(m2(), lp);
        [
            unit_to_u16(l_star),
            unit_to_u16(a + 0.5),
            unit_to_u16(b_star + 0.5),
        ]
    }

    fn forward_err(ch: [u16; 3]) -> u16 {
        let got = ::ph_color::srgb_to_oklab(::ph_color::Color::<
            ::ph_color::Srgb,
            ::ph_color::Linear,
        >::new(::ph_color::Q0_16::array_from_raw(ch)))
        .ch;
        let want = ottosson_srgb_to_oklab(ch);
        let [a0, a1, a2] = got;
        let [b0, b1, b2] = want;
        a0.to_raw()
            .abs_diff(b0)
            .max(a1.to_raw().abs_diff(b1))
            .max(a2.to_raw().abs_diff(b2))
    }

    #[test]
    fn forward_stays_within_asserted_bound() {
        const GRID: [u16; 9] = [0, 8192, 16384, 24576, 32768, 40960, 49152, 57344, 65535];
        let mut max = 0u16;
        for x in GRID {
            for y in GRID {
                for z in GRID {
                    max = max.max(forward_err([x, y, z]));
                }
            }
        }
        for x in 0..=u16::MAX {
            max = max.max(forward_err([x, 0, 0]));
            max = max.max(forward_err([0, x, 0]));
            max = max.max(forward_err([0, 0, x]));
        }
        max = max.max(forward_err([620, 0, 0]));
        max = max.max(forward_err([580, 0, 340]));
        let mut i = 0u16;
        loop {
            let x = i.saturating_mul(64);
            let mut j = 0u16;
            loop {
                let y = j.saturating_mul(64);
                max = max.max(forward_err([x, y, 0]));
                max = max.max(forward_err([x, 0, y]));
                max = max.max(forward_err([0, x, y]));
                if j == 1024 {
                    break;
                }
                j = j.saturating_add(1);
            }
            if i == 1024 {
                break;
            }
            i = i.saturating_add(1);
        }
        assert!(
            max <= ::ph_color::OKLAB_FORWARD_MAX_LSB,
            "forward max {max} > bound {}",
            ::ph_color::OKLAB_FORWARD_MAX_LSB
        );
        assert!(max > 0, "bound should be tight enough to be non-vacuous");
    }
}