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/// BT.709 → BT.2020 (row-major).
16///
17/// Converts linear BT.709/sRGB RGB to linear BT.2020 RGB.
18pub(crate) const BT709_TO_BT2020: GamutMatrix = [
19    [0.6274_0389, 0.3292_8303, 0.0433_1307],
20    [0.0690_9729, 0.9195_4040, 0.0113_6232],
21    [0.0163_9170, 0.0880_1327, 0.8955_9503],
22];
23
24/// BT.2020 → BT.709 (row-major).
25///
26/// Converts linear BT.2020 RGB to linear BT.709/sRGB RGB.
27/// Values outside \[0,1\] indicate out-of-gamut colors.
28pub(crate) const BT2020_TO_BT709: GamutMatrix = [
29    [1.6604_9100, -0.5876_5614, -0.0728_3486],
30    [-0.1245_5047, 1.1328_9990, -0.0083_4942],
31    [-0.0181_5076, -0.1005_7890, 1.1187_2966],
32];
33
34/// BT.709 → Display P3 (row-major).
35pub(crate) const BT709_TO_DISPLAY_P3: GamutMatrix = [
36    [0.8224_5811, 0.1775_4189, 0.0000_0000],
37    [0.0331_9419, 0.9668_0581, 0.0000_0000],
38    [0.0170_8263, 0.0723_9744, 0.9105_3993],
39];
40
41/// Display P3 → BT.709 (row-major).
42pub(crate) const DISPLAY_P3_TO_BT709: GamutMatrix = [
43    [1.2249_4018, -0.2249_4018, 0.0000_0000],
44    [-0.0420_4986, 1.0420_4986, 0.0000_0000],
45    [-0.0196_4113, -0.0786_4905, 1.0982_5018],
46];
47
48/// BT.2020 → Display P3 (row-major).
49pub(crate) const BT2020_TO_DISPLAY_P3: GamutMatrix = [
50    [1.3434_6376, -0.2826_7869, -0.0607_8507],
51    [-0.0652_8279, 1.0764_0361, -0.0111_2082],
52    [-0.0028_8423, -0.0193_4633, 1.0222_3056],
53];
54
55/// Display P3 → BT.2020 (row-major).
56pub(crate) const DISPLAY_P3_TO_BT2020: GamutMatrix = [
57    [0.7536_7740, 0.1985_4087, 0.0477_8174],
58    [0.0457_0150, 0.9417_7793, 0.0125_2057],
59    [0.0011_7409, 0.0176_4065, 0.9811_8526],
60];
61
62// ---------------------------------------------------------------------------
63// RGB ↔ CIE XYZ (D65 white point) matrices
64// ---------------------------------------------------------------------------
65
66/// BT.709/sRGB → CIE XYZ (D65). Derived from CIE 1931 xy chromaticities.
67pub(crate) const BT709_TO_XYZ: GamutMatrix = [
68    [0.4123907993, 0.3575843394, 0.1804807884],
69    [0.2126390059, 0.7151686788, 0.0721923154],
70    [0.0193308187, 0.1191947798, 0.9505321522],
71];
72
73/// CIE XYZ (D65) → BT.709/sRGB.
74pub(crate) const XYZ_TO_BT709: GamutMatrix = [
75    [3.2409699419, -1.5373831776, -0.4986107603],
76    [-0.9692436363, 1.8759675015, 0.0415550574],
77    [0.0556300797, -0.2039769589, 1.0569715142],
78];
79
80/// Display P3 → CIE XYZ (D65). Same white point as sRGB, wider primaries.
81pub(crate) const DISPLAY_P3_TO_XYZ: GamutMatrix = [
82    [0.4865709486, 0.2656676932, 0.1982172852],
83    [0.2289745641, 0.6917385218, 0.0792869141],
84    [0.0000000000, 0.0451133819, 1.0439443689],
85];
86
87/// CIE XYZ (D65) → Display P3.
88pub(crate) const XYZ_TO_DISPLAY_P3: GamutMatrix = [
89    [2.4934969119, -0.9313836179, -0.4027107845],
90    [-0.8294889696, 1.7626640603, 0.0236246858],
91    [0.0358458302, -0.0761723893, 0.9568845240],
92];
93
94/// BT.2020 → CIE XYZ (D65).
95pub(crate) const BT2020_TO_XYZ: GamutMatrix = [
96    [0.6369580484, 0.1446169036, 0.1688809752],
97    [0.2627002120, 0.6779980715, 0.0593017165],
98    [0.0000000000, 0.0280726930, 1.0609850578],
99];
100
101/// CIE XYZ (D65) → BT.2020.
102pub(crate) const XYZ_TO_BT2020: GamutMatrix = [
103    [1.7166511880, -0.3556707838, -0.2533662814],
104    [-0.6666843518, 1.6164812366, 0.0157685458],
105    [0.0176398574, -0.0427706133, 0.9421031212],
106];
107
108/// Multiply two 3×3 row-major matrices: `C = A × B`.
109pub fn mat3_mul(a: &GamutMatrix, b: &GamutMatrix) -> GamutMatrix {
110    let mut c = [[0.0f32; 3]; 3];
111    for i in 0..3 {
112        for j in 0..3 {
113            c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
114        }
115    }
116    c
117}
118
119/// Apply a 3×3 gamut matrix to a single linear RGB pixel (in-place).
120#[inline]
121pub fn apply_matrix_f32(rgb: &mut [f32; 3], m: &GamutMatrix) {
122    let [r, g, b] = *rgb;
123    rgb[0] = m[0][0] * r + m[0][1] * g + m[0][2] * b;
124    rgb[1] = m[1][0] * r + m[1][1] * g + m[1][2] * b;
125    rgb[2] = m[2][0] * r + m[2][1] * g + m[2][2] * b;
126}
127
128/// Apply a gamut matrix to a row of linear F32 RGB pixels.
129///
130/// `data` is `&mut [f32]` with `width * 3` elements (RGB, no alpha).
131/// For RGBA data, call [`apply_matrix_row_rgba_f32`] instead.
132pub fn apply_matrix_row_f32(data: &mut [f32], width: usize, m: &GamutMatrix) {
133    for i in 0..width {
134        let base = i * 3;
135        let r = data[base];
136        let g = data[base + 1];
137        let b = data[base + 2];
138        data[base] = m[0][0] * r + m[0][1] * g + m[0][2] * b;
139        data[base + 1] = m[1][0] * r + m[1][1] * g + m[1][2] * b;
140        data[base + 2] = m[2][0] * r + m[2][1] * g + m[2][2] * b;
141    }
142}
143
144/// Apply a gamut matrix to a row of linear F32 RGBA pixels.
145///
146/// Alpha channel is preserved unchanged.
147pub fn apply_matrix_row_rgba_f32(data: &mut [f32], width: usize, m: &GamutMatrix) {
148    for i in 0..width {
149        let base = i * 4;
150        let r = data[base];
151        let g = data[base + 1];
152        let b = data[base + 2];
153        data[base] = m[0][0] * r + m[0][1] * g + m[0][2] * b;
154        data[base + 1] = m[1][0] * r + m[1][1] * g + m[1][2] * b;
155        data[base + 2] = m[2][0] * r + m[2][1] * g + m[2][2] * b;
156        // data[base + 3] (alpha) unchanged.
157    }
158}
159
160/// Look up the gamut conversion matrix between two color primary sets.
161///
162/// Returns `None` if the primaries are the same or if no matrix is available.
163pub fn conversion_matrix(from: ColorPrimaries, to: ColorPrimaries) -> Option<&'static GamutMatrix> {
164    match (from, to) {
165        (ColorPrimaries::Bt709, ColorPrimaries::Bt2020) => Some(&BT709_TO_BT2020),
166        (ColorPrimaries::Bt2020, ColorPrimaries::Bt709) => Some(&BT2020_TO_BT709),
167        (ColorPrimaries::Bt709, ColorPrimaries::DisplayP3) => Some(&BT709_TO_DISPLAY_P3),
168        (ColorPrimaries::DisplayP3, ColorPrimaries::Bt709) => Some(&DISPLAY_P3_TO_BT709),
169        (ColorPrimaries::Bt2020, ColorPrimaries::DisplayP3) => Some(&BT2020_TO_DISPLAY_P3),
170        (ColorPrimaries::DisplayP3, ColorPrimaries::Bt2020) => Some(&DISPLAY_P3_TO_BT2020),
171        _ => None,
172    }
173}
174
175#[cfg(test)]
176mod tests {
177    use super::*;
178
179    /// Verify that forward × inverse ≈ identity for BT.709 ↔ BT.2020.
180    #[test]
181    fn bt709_bt2020_roundtrip() {
182        let test_rgb = [0.5f32, 0.3, 0.8];
183        let mut rgb = test_rgb;
184        apply_matrix_f32(&mut rgb, &BT709_TO_BT2020);
185        apply_matrix_f32(&mut rgb, &BT2020_TO_BT709);
186        for c in 0..3 {
187            assert!(
188                (rgb[c] - test_rgb[c]).abs() < 1e-5,
189                "BT.709→BT.2020→BT.709 roundtrip error in ch{c}: {:.6} vs {:.6}",
190                rgb[c],
191                test_rgb[c]
192            );
193        }
194    }
195
196    /// Verify that forward × inverse ≈ identity for BT.709 ↔ Display P3.
197    #[test]
198    fn bt709_displayp3_roundtrip() {
199        let test_rgb = [0.5f32, 0.3, 0.8];
200        let mut rgb = test_rgb;
201        apply_matrix_f32(&mut rgb, &BT709_TO_DISPLAY_P3);
202        apply_matrix_f32(&mut rgb, &DISPLAY_P3_TO_BT709);
203        for c in 0..3 {
204            assert!(
205                (rgb[c] - test_rgb[c]).abs() < 1e-5,
206                "BT.709→P3→BT.709 roundtrip error in ch{c}: {:.6} vs {:.6}",
207                rgb[c],
208                test_rgb[c]
209            );
210        }
211    }
212
213    /// White point preservation: [1,1,1] in BT.709 → BT.2020 should remain ~[1,1,1].
214    #[test]
215    fn white_point_preservation() {
216        let mut rgb = [1.0f32, 1.0, 1.0];
217        apply_matrix_f32(&mut rgb, &BT709_TO_BT2020);
218        for (c, &val) in rgb.iter().enumerate() {
219            assert!(
220                (val - 1.0).abs() < 1e-4,
221                "White point not preserved in ch{c}: {val:.6}",
222            );
223        }
224    }
225
226    /// RGBA row preserves alpha.
227    #[test]
228    fn rgba_alpha_preserved() {
229        let mut row = [0.5f32, 0.3, 0.8, 0.42, 0.1, 0.9, 0.2, 0.99];
230        apply_matrix_row_rgba_f32(&mut row, 2, &BT709_TO_BT2020);
231        assert_eq!(row[3], 0.42);
232        assert_eq!(row[7], 0.99);
233    }
234
235    /// Verify XYZ→RGB × RGB→XYZ ≈ identity for BT.709.
236    #[test]
237    fn xyz_bt709_roundtrip() {
238        let rgb = [0.5f32, 0.3, 0.8];
239        let mut v = rgb;
240        apply_matrix_f32(&mut v, &BT709_TO_XYZ);
241        apply_matrix_f32(&mut v, &XYZ_TO_BT709);
242        for c in 0..3 {
243            assert!(
244                (v[c] - rgb[c]).abs() < 1e-4,
245                "XYZ BT.709 roundtrip ch{c}: {:.6} vs {:.6}",
246                v[c],
247                rgb[c]
248            );
249        }
250    }
251
252    /// Verify XYZ→RGB × RGB→XYZ ≈ identity for Display P3.
253    #[test]
254    fn xyz_displayp3_roundtrip() {
255        let rgb = [0.5f32, 0.3, 0.8];
256        let mut v = rgb;
257        apply_matrix_f32(&mut v, &DISPLAY_P3_TO_XYZ);
258        apply_matrix_f32(&mut v, &XYZ_TO_DISPLAY_P3);
259        for c in 0..3 {
260            assert!(
261                (v[c] - rgb[c]).abs() < 1e-4,
262                "XYZ P3 roundtrip ch{c}: {:.6} vs {:.6}",
263                v[c],
264                rgb[c]
265            );
266        }
267    }
268
269    /// Verify XYZ→RGB × RGB→XYZ ≈ identity for BT.2020.
270    #[test]
271    fn xyz_bt2020_roundtrip() {
272        let rgb = [0.5f32, 0.3, 0.8];
273        let mut v = rgb;
274        apply_matrix_f32(&mut v, &BT2020_TO_XYZ);
275        apply_matrix_f32(&mut v, &XYZ_TO_BT2020);
276        for c in 0..3 {
277            assert!(
278                (v[c] - rgb[c]).abs() < 1e-4,
279                "XYZ BT.2020 roundtrip ch{c}: {:.6} vs {:.6}",
280                v[c],
281                rgb[c]
282            );
283        }
284    }
285
286    /// XYZ white point preservation: [1,1,1] → XYZ → RGB should remain ~[1,1,1].
287    #[test]
288    fn xyz_white_point() {
289        for (name, to, from) in [
290            ("BT.709", &BT709_TO_XYZ, &XYZ_TO_BT709),
291            ("P3", &DISPLAY_P3_TO_XYZ, &XYZ_TO_DISPLAY_P3),
292            ("BT.2020", &BT2020_TO_XYZ, &XYZ_TO_BT2020),
293        ] {
294            let mut rgb = [1.0f32; 3];
295            apply_matrix_f32(&mut rgb, to);
296            apply_matrix_f32(&mut rgb, from);
297            for (c, &val) in rgb.iter().enumerate() {
298                assert!(
299                    (val - 1.0).abs() < 1e-3,
300                    "{name} XYZ white point ch{c}: {val:.6}",
301                );
302            }
303        }
304    }
305
306    /// mat3_mul: verify A × A_inv ≈ identity.
307    #[test]
308    fn mat3_mul_inverse() {
309        let identity = mat3_mul(&BT709_TO_XYZ, &XYZ_TO_BT709);
310        for (i, row) in identity.iter().enumerate() {
311            for (j, &val) in row.iter().enumerate() {
312                let expected = if i == j { 1.0 } else { 0.0 };
313                assert!(
314                    (val - expected).abs() < 1e-4,
315                    "mat3_mul identity [{i}][{j}] = {val:.6}, expected {expected:.1}",
316                );
317            }
318        }
319    }
320
321    /// Verify cross-gamut via XYZ: BT.709→XYZ→BT.2020 ≈ direct BT.709→BT.2020.
322    #[test]
323    fn xyz_cross_gamut_consistency() {
324        let via_xyz = mat3_mul(&XYZ_TO_BT2020, &BT709_TO_XYZ);
325        let rgb = [0.5f32, 0.3, 0.8];
326        let mut v1 = rgb;
327        apply_matrix_f32(&mut v1, &via_xyz);
328        let mut v2 = rgb;
329        apply_matrix_f32(&mut v2, &BT709_TO_BT2020);
330        for c in 0..3 {
331            assert!(
332                (v1[c] - v2[c]).abs() < 1e-3,
333                "cross-gamut ch{c}: via_xyz={:.6} vs direct={:.6}",
334                v1[c],
335                v2[c]
336            );
337        }
338    }
339}
zenpixels_convert/gamut.rs

zenpixels_convert/
gamut.rs
