187 lines
6.7 KiB
Rust
187 lines
6.7 KiB
Rust
//! Small colour-science helpers: 3×3 matrices, CIE XYZ ↔ L*a*b*, chromatic adaptation.
|
||
//!
|
||
//! All formulas are from public standards: CIE 15 (L*a*b*), ICC.1:2010 Annex A (PCS, D50),
|
||
//! and the Bradford cone response matrix published by Lam (1985).
|
||
|
||
/// The ICC profile connection space illuminant (D50), as encoded in s15Fixed16 in every
|
||
/// ICC header: X = 0.9642, Y = 1.0, Z = 0.8249.
|
||
pub const D50: [f64; 3] = [0.9642, 1.0, 0.8249];
|
||
|
||
/// CIE D65 white point as XYZ (Y = 1), from the chromaticity (0.3127, 0.3290) used by sRGB,
|
||
/// Display P3 and Adobe RGB-compatible spaces.
|
||
pub const D65_XY: [f64; 2] = [0.3127, 0.3290];
|
||
|
||
/// Row-major 3×3 matrix.
|
||
pub type Mat3 = [[f64; 3]; 3];
|
||
|
||
pub const IDENTITY: Mat3 = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
|
||
|
||
pub fn mul(a: &Mat3, b: &Mat3) -> Mat3 {
|
||
let mut o = [[0.0; 3]; 3];
|
||
for (i, row) in o.iter_mut().enumerate() {
|
||
for (j, v) in row.iter_mut().enumerate() {
|
||
*v = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
|
||
}
|
||
}
|
||
o
|
||
}
|
||
|
||
pub fn apply(m: &Mat3, v: [f64; 3]) -> [f64; 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 fn det(m: &Mat3) -> f64 {
|
||
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])
|
||
}
|
||
|
||
/// Inverse, or `None` when (nearly) singular.
|
||
pub fn invert(m: &Mat3) -> Option<Mat3> {
|
||
let d = det(m);
|
||
if !d.is_finite() || d.abs() < 1e-12 {
|
||
return None;
|
||
}
|
||
let inv = 1.0 / d;
|
||
Some([
|
||
[(m[1][1] * m[2][2] - m[1][2] * m[2][1]) * inv, (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * inv, (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * inv],
|
||
[(m[1][2] * m[2][0] - m[1][0] * m[2][2]) * inv, (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * inv, (m[0][2] * m[1][0] - m[0][0] * m[1][2]) * inv],
|
||
[(m[1][0] * m[2][1] - m[1][1] * m[2][0]) * inv, (m[0][1] * m[2][0] - m[0][0] * m[2][1]) * inv, (m[0][0] * m[1][1] - m[0][1] * m[1][0]) * inv],
|
||
])
|
||
}
|
||
|
||
/// XYZ (Y = 1) of a chromaticity `(x, y)`.
|
||
pub fn xy_to_xyz(xy: [f64; 2]) -> [f64; 3] {
|
||
let [x, y] = xy;
|
||
[x / y, 1.0, (1.0 - x - y) / y]
|
||
}
|
||
|
||
/// Bradford cone response matrix (Lam 1985).
|
||
const BRADFORD: Mat3 = [[0.8951, 0.2664, -0.1614], [-0.7502, 1.7135, 0.0367], [0.0389, -0.0685, 1.0296]];
|
||
|
||
/// Bradford chromatic adaptation matrix from white `src` to white `dst` (both XYZ).
|
||
pub fn bradford(src: [f64; 3], dst: [f64; 3]) -> Mat3 {
|
||
let inv = invert(&BRADFORD).expect("Bradford matrix is invertible");
|
||
let s = apply(&BRADFORD, src);
|
||
let d = apply(&BRADFORD, dst);
|
||
let scale = [[d[0] / s[0], 0.0, 0.0], [0.0, d[1] / s[1], 0.0], [0.0, 0.0, d[2] / s[2]]];
|
||
mul(&inv, &mul(&scale, &BRADFORD))
|
||
}
|
||
|
||
/// RGB → XYZ matrix (relative to `white`, Y = 1) for primaries given as chromaticities.
|
||
pub fn rgb_to_xyz_matrix(r: [f64; 2], g: [f64; 2], b: [f64; 2], white: [f64; 3]) -> Mat3 {
|
||
let (xr, xg, xb) = (xy_to_xyz(r), xy_to_xyz(g), xy_to_xyz(b));
|
||
let m = [[xr[0], xg[0], xb[0]], [xr[1], xg[1], xb[1]], [xr[2], xg[2], xb[2]]];
|
||
let s = apply(&invert(&m).expect("primaries are independent"), white);
|
||
[[m[0][0] * s[0], m[0][1] * s[1], m[0][2] * s[2]], [m[1][0] * s[0], m[1][1] * s[1], m[1][2] * s[2]], [m[2][0] * s[0], m[2][1] * s[1], m[2][2] * s[2]]]
|
||
}
|
||
|
||
const EPS: f64 = 216.0 / 24389.0; // (6/29)^3
|
||
const KAPPA: f64 = 24389.0 / 27.0; // (29/3)^3
|
||
|
||
#[inline]
|
||
fn lab_f(t: f64) -> f64 {
|
||
if t > EPS { t.cbrt() } else { (KAPPA * t + 16.0) / 116.0 }
|
||
}
|
||
|
||
#[inline]
|
||
fn lab_finv(f: f64) -> f64 {
|
||
let f3 = f * f * f;
|
||
if f3 > EPS { f3 } else { (116.0 * f - 16.0) / KAPPA }
|
||
}
|
||
|
||
/// CIE XYZ → L*a*b* relative to `white` (CIE 15).
|
||
pub fn xyz_to_lab(xyz: [f64; 3], white: [f64; 3]) -> [f64; 3] {
|
||
let fx = lab_f(xyz[0] / white[0]);
|
||
let fy = lab_f(xyz[1] / white[1]);
|
||
let fz = lab_f(xyz[2] / white[2]);
|
||
[116.0 * fy - 16.0, 500.0 * (fx - fy), 200.0 * (fy - fz)]
|
||
}
|
||
|
||
/// CIE L*a*b* → XYZ relative to `white`.
|
||
pub fn lab_to_xyz(lab: [f64; 3], white: [f64; 3]) -> [f64; 3] {
|
||
let fy = (lab[0] + 16.0) / 116.0;
|
||
let fx = fy + lab[1] / 500.0;
|
||
let fz = fy - lab[2] / 200.0;
|
||
[lab_finv(fx) * white[0], lab_finv(fy) * white[1], lab_finv(fz) * white[2]]
|
||
}
|
||
|
||
/// f32 versions for the hot paths (D50 PCS only).
|
||
#[inline]
|
||
pub fn xyz_to_lab_f32(xyz: [f32; 3]) -> [f32; 3] {
|
||
const W: [f32; 3] = [0.9642, 1.0, 0.8249];
|
||
const E: f32 = 216.0 / 24389.0;
|
||
const K: f32 = 24389.0 / 27.0;
|
||
let f = |t: f32| if t > E { t.cbrt() } else { (K * t + 16.0) / 116.0 };
|
||
let fx = f(xyz[0] / W[0]);
|
||
let fy = f(xyz[1] / W[1]);
|
||
let fz = f(xyz[2] / W[2]);
|
||
[116.0 * fy - 16.0, 500.0 * (fx - fy), 200.0 * (fy - fz)]
|
||
}
|
||
|
||
#[inline]
|
||
pub fn lab_to_xyz_f32(lab: [f32; 3]) -> [f32; 3] {
|
||
const W: [f32; 3] = [0.9642, 1.0, 0.8249];
|
||
const E: f32 = 216.0 / 24389.0;
|
||
const K: f32 = 24389.0 / 27.0;
|
||
let fy = (lab[0] + 16.0) / 116.0;
|
||
let fx = fy + lab[1] / 500.0;
|
||
let fz = fy - lab[2] / 200.0;
|
||
let finv = |f: f32| {
|
||
let f3 = f * f * f;
|
||
if f3 > E { f3 } else { (116.0 * f - 16.0) / K }
|
||
};
|
||
[finv(fx) * W[0], finv(fy) * W[1], finv(fz) * W[2]]
|
||
}
|
||
|
||
/// CIE 1976 colour difference.
|
||
pub fn delta_e76(a: [f64; 3], b: [f64; 3]) -> f64 {
|
||
((a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2) + (a[2] - b[2]).powi(2)).sqrt()
|
||
}
|
||
|
||
/// Quantizes to ICC s15Fixed16Number and back (what a profile can actually store).
|
||
pub fn s15f16_round(v: f64) -> f64 {
|
||
(v * 65536.0).round() / 65536.0
|
||
}
|
||
|
||
/// The D50 illuminant exactly as an ICC header stores it (s15Fixed16).
|
||
pub fn d50_quantized() -> [f64; 3] {
|
||
[s15f16_round(D50[0]), s15f16_round(D50[1]), s15f16_round(D50[2])]
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
|
||
#[test]
|
||
fn lab_roundtrip() {
|
||
for xyz in [[0.2, 0.3, 0.1], [0.9642, 1.0, 0.8249], [0.001, 0.002, 0.0005]] {
|
||
let back = lab_to_xyz(xyz_to_lab(xyz, D50), D50);
|
||
for i in 0..3 {
|
||
assert!((back[i] - xyz[i]).abs() < 1e-12);
|
||
}
|
||
}
|
||
let w = xyz_to_lab(D50, D50);
|
||
assert!((w[0] - 100.0).abs() < 1e-12 && w[1].abs() < 1e-12 && w[2].abs() < 1e-12);
|
||
}
|
||
|
||
#[test]
|
||
fn bradford_maps_white() {
|
||
let d65 = xy_to_xyz(D65_XY);
|
||
let m = bradford(d65, D50);
|
||
let w = apply(&m, d65);
|
||
for i in 0..3 {
|
||
assert!((w[i] - D50[i]).abs() < 1e-12);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn invert_roundtrip() {
|
||
let m = [[0.4, 0.3, 0.2], [0.2, 0.7, 0.1], [0.0, 0.1, 0.8]];
|
||
let p = mul(&m, &invert(&m).unwrap());
|
||
for (i, row) in p.iter().enumerate() {
|
||
for (j, v) in row.iter().enumerate() {
|
||
assert!((v - if i == j { 1.0 } else { 0.0 }).abs() < 1e-12);
|
||
}
|
||
}
|
||
}
|
||
}
|