ph-color-bake 0.1.0

//! RGB↔XYZ and RGB↔RGB matrix derivation on the host.

use crate::BakeError;
use crate::Primaries;

/// Q4.28 encoding of 1.0. This host-side scale is checked against
/// `ph_color::Q4_28::ONE` by the bake tests.
pub const COEF_ONE: i32 = 1 << 28;

/// 3×3 row-major `f64` matrix.
pub type Mat3 = [[f64; 3]; 3];

/// Convert a real coefficient to saturating Q4.28 with round-half-away-from-zero.
#[must_use]
pub fn to_q428(value: f64) -> i32 {
    let scaled = value * f64::from(COEF_ONE);
    let rounded = if scaled >= 0.0 {
        (scaled + 0.5).floor()
    } else {
        (scaled - 0.5).ceil()
    };
    rounded.clamp(f64::from(i32::MIN), f64::from(i32::MAX)) as i32
}

/// Identity Q4.28 matrix.
#[must_use]
pub fn identity_q428() -> [[i32; 3]; 3] {
    [[COEF_ONE, 0, 0], [0, COEF_ONE, 0], [0, 0, COEF_ONE]]
}

/// Derive linear RGB→RGB from source and destination primaries.
///
/// Same primaries and white yield an exact identity Q4.28 matrix. Otherwise
/// the path is RGB_src → XYZ → RGB_dst using `f64` inversion on the host.
pub fn rgb_to_rgb(src: Primaries, dst: Primaries) -> Result<[[i32; 3]; 3], BakeError> {
    let src_to_xyz = rgb_to_xyz(src)?;
    if src == dst {
        return Ok(identity_q428());
    }
    let dst_to_xyz = rgb_to_xyz(dst)?;
    let xyz_to_dst = invert(dst_to_xyz)?;
    let composed = mul(xyz_to_dst, src_to_xyz);
    Ok([
        [
            to_q428(composed[0][0]),
            to_q428(composed[0][1]),
            to_q428(composed[0][2]),
        ],
        [
            to_q428(composed[1][0]),
            to_q428(composed[1][1]),
            to_q428(composed[1][2]),
        ],
        [
            to_q428(composed[2][0]),
            to_q428(composed[2][1]),
            to_q428(composed[2][2]),
        ],
    ])
}

fn xy_to_xyz(xy: crate::Xy) -> Result<[f64; 3], BakeError> {
    if !xy.y.is_finite() || xy.y.abs() < 1e-12 {
        return Err(BakeError::SingularMatrix);
    }
    let y = 1.0;
    let x = xy.x / xy.y;
    let z = (1.0 - xy.x - xy.y) / xy.y;
    if ![x, y, z].iter().all(|c| c.is_finite()) {
        return Err(BakeError::SingularMatrix);
    }
    Ok([x * y, y, z * y])
}

fn rgb_to_xyz(p: Primaries) -> Result<Mat3, BakeError> {
    let r = xy_to_xyz(p.r)?;
    let g = xy_to_xyz(p.g)?;
    let b = xy_to_xyz(p.b)?;
    let w = xy_to_xyz(p.white)?;
    let prim = [[r[0], g[0], b[0]], [r[1], g[1], b[1]], [r[2], g[2], b[2]]];
    let inv = invert(prim)?;
    let s = mulv(inv, w);
    let scaled = [
        [prim[0][0] * s[0], prim[0][1] * s[1], prim[0][2] * s[2]],
        [prim[1][0] * s[0], prim[1][1] * s[1], prim[1][2] * s[2]],
        [prim[2][0] * s[0], prim[2][1] * s[1], prim[2][2] * s[2]],
    ];
    let _ = invert(scaled)?;
    Ok(scaled)
}

fn mul(a: Mat3, b: Mat3) -> Mat3 {
    let mut out = [[0.0; 3]; 3];
    for i in 0..3 {
        for j in 0..3 {
            out[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
        }
    }
    out
}

fn mulv(a: Mat3, v: [f64; 3]) -> [f64; 3] {
    [
        a[0][0] * v[0] + a[0][1] * v[1] + a[0][2] * v[2],
        a[1][0] * v[0] + a[1][1] * v[1] + a[1][2] * v[2],
        a[2][0] * v[0] + a[2][1] * v[1] + a[2][2] * v[2],
    ]
}

fn invert(m: Mat3) -> Result<Mat3, BakeError> {
    let det = m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
        - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
        + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
    if !det.is_finite() || det.abs() < 1e-12 {
        return Err(BakeError::SingularMatrix);
    }
    let inv_det = 1.0 / det;
    Ok([
        [
            (m[1][1] * m[2][2] - m[1][2] * m[2][1]) * inv_det,
            (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * inv_det,
            (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * inv_det,
        ],
        [
            (m[1][2] * m[2][0] - m[1][0] * m[2][2]) * inv_det,
            (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * inv_det,
            (m[0][2] * m[1][0] - m[0][0] * m[1][2]) * inv_det,
        ],
        [
            (m[1][0] * m[2][1] - m[1][1] * m[2][0]) * inv_det,
            (m[0][1] * m[2][0] - m[0][0] * m[2][1]) * inv_det,
            (m[0][0] * m[1][1] - m[0][1] * m[1][0]) * inv_det,
        ],
    ])
}

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

    #[test]
    fn same_space_is_identity() {
        let m = rgb_to_rgb(Primaries::srgb(), Primaries::srgb()).expect("sRGB identity");
        assert_eq!(m, identity_q428());
    }

    #[test]
    fn white_on_a_primary_is_singular() {
        let mut bad = Primaries::srgb();
        bad.white = bad.r;
        assert_eq!(
            rgb_to_rgb(bad, Primaries::srgb()),
            Err(crate::BakeError::SingularMatrix)
        );
    }

    #[test]
    fn duplicate_primaries_are_singular() {
        let mut bad = Primaries::srgb();
        bad.g = bad.r;
        assert_eq!(
            rgb_to_rgb(Primaries::srgb(), bad),
            Err(crate::BakeError::SingularMatrix)
        );
    }

    #[test]
    fn identical_invalid_primaries_are_singular() {
        let mut bad = Primaries::srgb();
        bad.g = bad.r;
        assert_eq!(rgb_to_rgb(bad, bad), Err(crate::BakeError::SingularMatrix));
    }
}