zenpixels_convert/
gamut.rs

1#![allow(clippy::excessive_precision)]
2//! Color gamut conversion matrices for BT.709, BT.2020, and Display P3.
3//!
4//! All matrices operate in **linear light** — apply EOTF first, then gamut
5//! matrix, then OETF. The matrices are computed from the CIE 1931 xy
6//! chromaticity coordinates of each primary set using D65 white point.
7//!
8//! These are exact, standards-derived values (not fitted or approximated).
9
10use crate::ColorPrimaries;
11
12/// A 3×3 row-major matrix for linear RGB ↔ linear RGB gamut conversion.
13pub type GamutMatrix = [[f32; 3]; 3];
14
15// ---------------------------------------------------------------------------
16// RGB ↔ CIE XYZ (D65 white point) matrices
17// ---------------------------------------------------------------------------
18
19/// BT.709/sRGB → CIE XYZ (D65). Derived from CIE 1931 xy chromaticities.
20pub(crate) const BT709_TO_XYZ: GamutMatrix = [
21    [0.4123907993, 0.3575843394, 0.1804807884],
22    [0.2126390059, 0.7151686788, 0.0721923154],
23    [0.0193308187, 0.1191947798, 0.9505321522],
24];
25
26/// CIE XYZ (D65) → BT.709/sRGB.
27pub(crate) const XYZ_TO_BT709: GamutMatrix = [
28    [3.2409699419, -1.5373831776, -0.4986107603],
29    [-0.9692436363, 1.8759675015, 0.0415550574],
30    [0.0556300797, -0.2039769589, 1.0569715142],
31];
32
33/// Display P3 → CIE XYZ (D65). Same white point as sRGB, wider primaries.
34pub(crate) const DISPLAY_P3_TO_XYZ: GamutMatrix = [
35    [0.4865709486, 0.2656676932, 0.1982172852],
36    [0.2289745641, 0.6917385218, 0.0792869141],
37    [0.0000000000, 0.0451133819, 1.0439443689],
38];
39
40/// CIE XYZ (D65) → Display P3.
41pub(crate) const XYZ_TO_DISPLAY_P3: GamutMatrix = [
42    [2.4934969119, -0.9313836179, -0.4027107845],
43    [-0.8294889696, 1.7626640603, 0.0236246858],
44    [0.0358458302, -0.0761723893, 0.9568845240],
45];
46
47/// BT.2020 → CIE XYZ (D65).
48pub(crate) const BT2020_TO_XYZ: GamutMatrix = [
49    [0.6369580484, 0.1446169036, 0.1688809752],
50    [0.2627002120, 0.6779980715, 0.0593017165],
51    [0.0000000000, 0.0280726930, 1.0609850578],
52];
53
54/// CIE XYZ (D65) → BT.2020.
55pub(crate) const XYZ_TO_BT2020: GamutMatrix = [
56    [1.7166511880, -0.3556707838, -0.2533662814],
57    [-0.6666843518, 1.6164812366, 0.0157685458],
58    [0.0176398574, -0.0427706133, 0.9421031212],
59];
60
61/// Multiply two 3×3 row-major matrices: `C = A × B`.
62pub fn mat3_mul(a: &GamutMatrix, b: &GamutMatrix) -> GamutMatrix {
63    let mut c = [[0.0f32; 3]; 3];
64    for i in 0..3 {
65        for j in 0..3 {
66            c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
67        }
68    }
69    c
70}
71
72/// Apply a 3×3 gamut matrix to a single linear RGB pixel (in-place).
73#[inline]
74pub fn apply_matrix_f32(rgb: &mut [f32; 3], m: &GamutMatrix) {
75    let [r, g, b] = *rgb;
76    rgb[0] = m[0][0] * r + m[0][1] * g + m[0][2] * b;
77    rgb[1] = m[1][0] * r + m[1][1] * g + m[1][2] * b;
78    rgb[2] = m[2][0] * r + m[2][1] * g + m[2][2] * b;
79}
80
81/// Apply a gamut matrix to a row of linear F32 RGB pixels.
82///
83/// `data` is `&mut [f32]` with `width * 3` elements (RGB, no alpha).
84/// For RGBA data, call [`apply_matrix_row_rgba_f32`] instead.
85pub fn apply_matrix_row_f32(data: &mut [f32], width: usize, m: &GamutMatrix) {
86    for i in 0..width {
87        let base = i * 3;
88        let r = data[base];
89        let g = data[base + 1];
90        let b = data[base + 2];
91        data[base] = m[0][0] * r + m[0][1] * g + m[0][2] * b;
92        data[base + 1] = m[1][0] * r + m[1][1] * g + m[1][2] * b;
93        data[base + 2] = m[2][0] * r + m[2][1] * g + m[2][2] * b;
94    }
95}
96
97/// Apply a gamut matrix to a row of linear F32 RGBA pixels.
98///
99/// Alpha channel is preserved unchanged.
100pub fn apply_matrix_row_rgba_f32(data: &mut [f32], width: usize, m: &GamutMatrix) {
101    for i in 0..width {
102        let base = i * 4;
103        let r = data[base];
104        let g = data[base + 1];
105        let b = data[base + 2];
106        data[base] = m[0][0] * r + m[0][1] * g + m[0][2] * b;
107        data[base + 1] = m[1][0] * r + m[1][1] * g + m[1][2] * b;
108        data[base + 2] = m[2][0] * r + m[2][1] * g + m[2][2] * b;
109        // data[base + 3] (alpha) unchanged.
110    }
111}
112
113/// Look up the gamut conversion matrix between two color primary sets.
114///
115/// Delegates to [`ColorPrimaries::gamut_matrix_to`], which computes the matrix
116/// from chromaticity coordinates with Bradford chromatic adaptation when needed.
117/// Returns `None` if the primaries are the same or if either is `Unknown`.
118pub fn conversion_matrix(from: ColorPrimaries, to: ColorPrimaries) -> Option<GamutMatrix> {
119    if from as u8 == to as u8 {
120        return None;
121    }
122    from.gamut_matrix_to(to)
123}
124
125#[cfg(test)]
126mod tests {
127    use super::*;
128
129    use crate::ColorPrimaries;
130
131    /// Verify that forward × inverse ≈ identity for BT.709 ↔ BT.2020.
132    #[test]
133    fn bt709_bt2020_roundtrip() {
134        let fwd = ColorPrimaries::Bt709
135            .gamut_matrix_to(ColorPrimaries::Bt2020)
136            .unwrap();
137        let inv = ColorPrimaries::Bt2020
138            .gamut_matrix_to(ColorPrimaries::Bt709)
139            .unwrap();
140        let test_rgb = [0.5f32, 0.3, 0.8];
141        let mut rgb = test_rgb;
142        apply_matrix_f32(&mut rgb, &fwd);
143        apply_matrix_f32(&mut rgb, &inv);
144        for c in 0..3 {
145            assert!(
146                (rgb[c] - test_rgb[c]).abs() < 1e-4,
147                "BT.709→BT.2020→BT.709 roundtrip error in ch{c}: {:.6} vs {:.6}",
148                rgb[c],
149                test_rgb[c]
150            );
151        }
152    }
153
154    /// Verify that forward × inverse ≈ identity for BT.709 ↔ Display P3.
155    #[test]
156    fn bt709_displayp3_roundtrip() {
157        let fwd = ColorPrimaries::Bt709
158            .gamut_matrix_to(ColorPrimaries::DisplayP3)
159            .unwrap();
160        let inv = ColorPrimaries::DisplayP3
161            .gamut_matrix_to(ColorPrimaries::Bt709)
162            .unwrap();
163        let test_rgb = [0.5f32, 0.3, 0.8];
164        let mut rgb = test_rgb;
165        apply_matrix_f32(&mut rgb, &fwd);
166        apply_matrix_f32(&mut rgb, &inv);
167        for c in 0..3 {
168            assert!(
169                (rgb[c] - test_rgb[c]).abs() < 1e-4,
170                "BT.709→P3→BT.709 roundtrip error in ch{c}: {:.6} vs {:.6}",
171                rgb[c],
172                test_rgb[c]
173            );
174        }
175    }
176
177    /// White point preservation: [1,1,1] in BT.709 → BT.2020 should remain ~[1,1,1].
178    #[test]
179    fn white_point_preservation() {
180        let m = ColorPrimaries::Bt709
181            .gamut_matrix_to(ColorPrimaries::Bt2020)
182            .unwrap();
183        let mut rgb = [1.0f32, 1.0, 1.0];
184        apply_matrix_f32(&mut rgb, &m);
185        for (c, &val) in rgb.iter().enumerate() {
186            assert!(
187                (val - 1.0).abs() < 1e-4,
188                "White point not preserved in ch{c}: {val:.6}",
189            );
190        }
191    }
192
193    /// RGBA row preserves alpha.
194    #[test]
195    fn rgba_alpha_preserved() {
196        let m = ColorPrimaries::Bt709
197            .gamut_matrix_to(ColorPrimaries::Bt2020)
198            .unwrap();
199        let mut row = [0.5f32, 0.3, 0.8, 0.42, 0.1, 0.9, 0.2, 0.99];
200        apply_matrix_row_rgba_f32(&mut row, 2, &m);
201        assert_eq!(row[3], 0.42);
202        assert_eq!(row[7], 0.99);
203    }
204
205    /// Verify XYZ→RGB × RGB→XYZ ≈ identity for BT.709.
206    #[test]
207    fn xyz_bt709_roundtrip() {
208        let rgb = [0.5f32, 0.3, 0.8];
209        let mut v = rgb;
210        apply_matrix_f32(&mut v, &BT709_TO_XYZ);
211        apply_matrix_f32(&mut v, &XYZ_TO_BT709);
212        for c in 0..3 {
213            assert!(
214                (v[c] - rgb[c]).abs() < 1e-4,
215                "XYZ BT.709 roundtrip ch{c}: {:.6} vs {:.6}",
216                v[c],
217                rgb[c]
218            );
219        }
220    }
221
222    /// Verify XYZ→RGB × RGB→XYZ ≈ identity for Display P3.
223    #[test]
224    fn xyz_displayp3_roundtrip() {
225        let rgb = [0.5f32, 0.3, 0.8];
226        let mut v = rgb;
227        apply_matrix_f32(&mut v, &DISPLAY_P3_TO_XYZ);
228        apply_matrix_f32(&mut v, &XYZ_TO_DISPLAY_P3);
229        for c in 0..3 {
230            assert!(
231                (v[c] - rgb[c]).abs() < 1e-4,
232                "XYZ P3 roundtrip ch{c}: {:.6} vs {:.6}",
233                v[c],
234                rgb[c]
235            );
236        }
237    }
238
239    /// Verify XYZ→RGB × RGB→XYZ ≈ identity for BT.2020.
240    #[test]
241    fn xyz_bt2020_roundtrip() {
242        let rgb = [0.5f32, 0.3, 0.8];
243        let mut v = rgb;
244        apply_matrix_f32(&mut v, &BT2020_TO_XYZ);
245        apply_matrix_f32(&mut v, &XYZ_TO_BT2020);
246        for c in 0..3 {
247            assert!(
248                (v[c] - rgb[c]).abs() < 1e-4,
249                "XYZ BT.2020 roundtrip ch{c}: {:.6} vs {:.6}",
250                v[c],
251                rgb[c]
252            );
253        }
254    }
255
256    /// XYZ white point preservation: [1,1,1] → XYZ → RGB should remain ~[1,1,1].
257    #[test]
258    fn xyz_white_point() {
259        for (name, to, from) in [
260            ("BT.709", &BT709_TO_XYZ, &XYZ_TO_BT709),
261            ("P3", &DISPLAY_P3_TO_XYZ, &XYZ_TO_DISPLAY_P3),
262            ("BT.2020", &BT2020_TO_XYZ, &XYZ_TO_BT2020),
263        ] {
264            let mut rgb = [1.0f32; 3];
265            apply_matrix_f32(&mut rgb, to);
266            apply_matrix_f32(&mut rgb, from);
267            for (c, &val) in rgb.iter().enumerate() {
268                assert!(
269                    (val - 1.0).abs() < 1e-3,
270                    "{name} XYZ white point ch{c}: {val:.6}",
271                );
272            }
273        }
274    }
275
276    /// mat3_mul: verify A × A_inv ≈ identity.
277    #[test]
278    fn mat3_mul_inverse() {
279        let identity = mat3_mul(&BT709_TO_XYZ, &XYZ_TO_BT709);
280        for (i, row) in identity.iter().enumerate() {
281            for (j, &val) in row.iter().enumerate() {
282                let expected = if i == j { 1.0 } else { 0.0 };
283                assert!(
284                    (val - expected).abs() < 1e-4,
285                    "mat3_mul identity [{i}][{j}] = {val:.6}, expected {expected:.1}",
286                );
287            }
288        }
289    }
290
291    /// Verify cross-gamut via XYZ: BT.709→XYZ→BT.2020 ≈ direct BT.709→BT.2020.
292    #[test]
293    fn xyz_cross_gamut_consistency() {
294        let via_xyz = mat3_mul(&XYZ_TO_BT2020, &BT709_TO_XYZ);
295        let direct = ColorPrimaries::Bt709
296            .gamut_matrix_to(ColorPrimaries::Bt2020)
297            .unwrap();
298        let rgb = [0.5f32, 0.3, 0.8];
299        let mut v1 = rgb;
300        apply_matrix_f32(&mut v1, &via_xyz);
301        let mut v2 = rgb;
302        apply_matrix_f32(&mut v2, &direct);
303        for c in 0..3 {
304            assert!(
305                (v1[c] - v2[c]).abs() < 1e-3,
306                "cross-gamut ch{c}: via_xyz={:.6} vs direct={:.6}",
307                v1[c],
308                v2[c]
309            );
310        }
311    }
312}
zenpixels_convert/gamut.rs

zenpixels_convert/
gamut.rs
