Skip to main content

mcraw_tui/
color.rs

1use crate::agx::{AgxConfig, AgxPipeline, Gamut, OutputTransfer, Transfer};
2use crate::file::BayerPattern;
3use anyhow::Result;
4use rayon::prelude::*;
5
6pub trait Demosaic {
7    fn process(&self, bayer: &[u16], stride_width: u32, offset_x: u32, offset_y: u32, active_width: u32, active_height: u32, pattern: &BayerPattern) -> Result<Vec<f32>>;
8}
9
10pub trait ColorSpaceConverter {
11    fn process(&self, pixels: &mut [f32], ccm: &[f32; 9]);
12}
13
14pub trait TransferFunctionProcessor {
15    fn process(&self, pixels: &mut [f32]);
16}
17
18#[derive(Debug, Clone, Copy, PartialEq, Eq)]
19pub enum ColorSpace {
20    ACESAP1, AppleWideGamut, ARRIWideGamut3, ARRIWideGamut4, CanonCinemaGamut,
21    DaVinciWideGamut, DciP3, DisplayP3, FGamut, FGamutC, PanasonicVGamut, Rec2020,
22    Rec709, SGamut3, SGamut3Cine, Srgb,
23}
24
25impl ColorSpace {
26    pub fn name(&self) -> &'static str {
27        match self {
28            ColorSpace::ACESAP1 => "ACES AP1",
29            ColorSpace::AppleWideGamut => "Apple Wide Gamut",
30            ColorSpace::ARRIWideGamut3 => "ARRI Wide Gamut 3", ColorSpace::ARRIWideGamut4 => "ARRI Wide Gamut 4",
31            ColorSpace::CanonCinemaGamut => "Canon Cinema Gamut",
32            ColorSpace::DaVinciWideGamut => "DaVinci Wide Gamut",
33            ColorSpace::DciP3 => "DCI-P3", ColorSpace::DisplayP3 => "Display P3",
34            ColorSpace::FGamut => "F-Gamut", ColorSpace::FGamutC => "F-Gamut C",
35            ColorSpace::PanasonicVGamut => "Panasonic V-Gamut",
36            ColorSpace::Rec2020 => "Rec.2020", ColorSpace::Rec709 => "Rec.709",
37            ColorSpace::SGamut3 => "S-Gamut3", ColorSpace::SGamut3Cine => "S-Gamut3.Cinema",
38            ColorSpace::Srgb => "sRGB",
39        }
40    }
41
42    pub fn get_white_point_chromaticities(&self) -> (f32, f32) {
43        match self {
44            ColorSpace::DciP3 => (0.314, 0.351),
45            ColorSpace::ACESAP1 => (0.32168, 0.33767),
46            _ => (0.3127, 0.3290),
47        }
48    }
49
50    pub fn get_xyz_to_rgb_matrix(&self) -> [f32; 9] {
51        match self {
52            ColorSpace::AppleWideGamut => xyz_to_rgb_from_primaries(0.725, 0.301, 0.221, 0.814, 0.068, -0.076, 0.3127, 0.3290),
53            ColorSpace::Rec709 | ColorSpace::Srgb => xyz_to_rec709(),
54            ColorSpace::Rec2020 | ColorSpace::FGamut => xyz_to_rgb_from_primaries(0.708, 0.292, 0.170, 0.797, 0.131, 0.046, 0.3127, 0.3290),
55            ColorSpace::DciP3 => xyz_to_rgb_from_primaries(0.680, 0.320, 0.265, 0.690, 0.150, 0.060, 0.314, 0.351),
56            ColorSpace::DisplayP3 => xyz_to_rgb_from_primaries(0.680, 0.320, 0.265, 0.690, 0.150, 0.060, 0.3127, 0.3290),
57            ColorSpace::SGamut3Cine => xyz_to_rgb_from_primaries(0.76600, 0.27500, 0.22500, 0.80000, 0.08900, -0.08700, 0.3127, 0.3290),
58            ColorSpace::SGamut3 => xyz_to_rgb_from_primaries(0.7300, 0.2800, 0.1400, 0.8550, 0.1000, -0.0500, 0.3127, 0.3290),
59            ColorSpace::ARRIWideGamut3 => xyz_to_rgb_from_primaries(0.6840, 0.3130, 0.2210, 0.8480, 0.0861, -0.1020, 0.3127, 0.3290),
60            ColorSpace::ARRIWideGamut4 => xyz_to_rgb_from_primaries(0.7347, 0.2653, 0.1424, 0.8576, 0.0991, -0.0308, 0.3127, 0.3290),
61            ColorSpace::CanonCinemaGamut => xyz_to_rgb_from_primaries(0.7400, 0.2700, 0.1700, 1.1400, 0.0800, -0.1000, 0.3127, 0.3290),
62            ColorSpace::PanasonicVGamut => xyz_to_rgb_from_primaries(0.7300, 0.2800, 0.1650, 0.8400, 0.1000, -0.0300, 0.3127, 0.3290),
63            ColorSpace::FGamutC => xyz_to_rgb_from_primaries(0.7347, 0.2653, 0.0263, 0.9737, 0.1173, -0.0224, 0.3127, 0.3290),
64            ColorSpace::DaVinciWideGamut => xyz_to_rgb_from_primaries(0.8000, 0.3130, 0.1682, 0.9877, 0.0790, -0.1155, 0.3127, 0.3290),
65            ColorSpace::ACESAP1 => xyz_to_rgb_from_primaries(0.71300, 0.29300, 0.16500, 0.83000, 0.12800, 0.04400, 0.32168, 0.33767),
66        }
67    }
68
69    pub fn all() -> &'static [ColorSpace] {
70        // Alphabetical order for deterministic, pleasing cycle order.
71        &[ColorSpace::ACESAP1, ColorSpace::AppleWideGamut, ColorSpace::ARRIWideGamut3, ColorSpace::ARRIWideGamut4,
72          ColorSpace::CanonCinemaGamut, ColorSpace::DaVinciWideGamut, ColorSpace::DciP3,
73          ColorSpace::DisplayP3, ColorSpace::FGamut, ColorSpace::FGamutC,
74          ColorSpace::PanasonicVGamut, ColorSpace::Rec2020, ColorSpace::Rec709,
75          ColorSpace::SGamut3, ColorSpace::SGamut3Cine, ColorSpace::Srgb]
76    }
77    pub fn next(self) -> Self { let all = Self::all(); let pos = all.iter().position(|&x| x == self).unwrap_or(0); all[(pos + 1) % all.len()] }
78    pub fn prev(self) -> Self { let all = Self::all(); let pos = all.iter().position(|&x| x == self).unwrap_or(0); all[(pos + all.len() - 1) % all.len()] }
79}
80
81#[derive(Debug, Clone, Copy, PartialEq, Eq)]
82pub enum TransferFunction {
83    ACESCCT, ARRIlog3, ARRIlog4, AppleLog, AppleLog2, CLog3, DaVinciIntermediate,
84    FLog2, Gamma24, HLG, Linear, PQ, Rec709, SLog3, VLog,
85}
86
87impl TransferFunction {
88    pub fn name(&self) -> &'static str {
89        match self {
90            TransferFunction::ACESCCT => "ACES CCT",
91            TransferFunction::ARRIlog3 => "ARRI LogC3", TransferFunction::ARRIlog4 => "ARRI LogC4",
92            TransferFunction::AppleLog => "Apple Log", TransferFunction::AppleLog2 => "Apple Log 2",
93            TransferFunction::CLog3 => "C-Log3",
94            TransferFunction::DaVinciIntermediate => "DaVinci Intermediate",
95            TransferFunction::FLog2 => "F-Log2", TransferFunction::Gamma24 => "Gamma 2.4",
96            TransferFunction::HLG => "HLG (BT.2100)", TransferFunction::Linear => "Linear",
97            TransferFunction::PQ => "PQ (ST.2084)", TransferFunction::Rec709 => "Rec.709",
98            TransferFunction::SLog3 => "S-Log3", TransferFunction::VLog => "V-Log",
99        }
100    }
101
102    /// Apply the OETF (linear → log) for the selected transfer function.
103    ///
104    /// **Source-of-truth references** for each branch:
105    ///
106    /// | Variant | Spec / document |
107    /// |---|---|
108    /// | `Rec709`         | ITU-R BT.709-6 OETF |
109    /// | `SLog3`          | Sony "S-Log3 Technical Specification" (Sept 2014) — canonical form: code = `(420 + 261.5×log₁₀((x+0.01)/0.19)) / 1023`, knee at `0.01125`, black code `95`, 18% grey code `420` |
110    /// | `VLog`           | Panasonic "V-Log/V-Gamut Reference Manual" (2014) — `5.6x+0.125` / `0.241514*log10(x+0.00873)+0.598206`, knee at `0.01` |
111    /// | `ARRIlog3`       | ARRI "LogC-3 Logarithmic Color Space" spec (2020), EI 800 variant |
112    /// | `ARRIlog4`       | ARRI "LogC4 Encoding Function" (Cooper & Brendel, 2022; ALEV4 / Alexa 35), EI-independent |
113    /// | `CLog3`          | Canon Cinema EOS C-Log3 characteristics (2016) — three-segment with negative-side graft |
114    /// | `FLog2`          | Fujifilm "F-Log2 Data Sheet" (2021) — Fujifilm-internal anchor at `0.000889` |
115    /// | `AppleLog`/`AppleLog2` | Apple "Apple Log Profile White Paper" (Sept 2023) — `R0=-0.05641088`, `C=47.28711236` |
116    /// | `ACESCCT`        | AMPAS ACEScc specification (TB-2022-002), knee at `2^-7 = 0.0078125`, log slope `17.52` |
117    /// | `PQ`             | ITU-R BT.2100-2 ST.2084 PQ (2022) — `m1=0.1593017578125`, `m2=78.84375`, `c1=0.8359375`, `c2=18.8515625`, `c3=18.6875` |
118    /// | `HLG`            | ITU-R BT.2100-2 HLG OETF (2022) — knee at `1/12`, `a=0.17883277`, `b=0.28466892`, `c=0.55991073` |
119    /// | `DaVinciIntermediate` | Blackmagic "DaVinci YRGB Intermediate" — knee at `0.00262409`, log slope `0.07329248` |
120    /// | `Gamma24`        | Display gamma `1/2.4` (Rec.1886 EOTF approximation) |
121    /// | `Linear`         | identity |
122    pub fn process(&self, pixels: &mut [f32]) {
123        match self {
124            TransferFunction::Linear => {}
125            // Source: ITU-R BT.709-6 §3.
126            TransferFunction::Rec709 => { pixels.par_iter_mut().for_each(|v| { *v = rec709_oetf(*v).min(1.0).max(0.0); }); }
127            // Source: Sony "S-Log3 Technical Summary" (Sept 2014).
128            // Canonical form per colour-science and ACES CTL ref.
129            // Knee at 0.01125; above: log segment maps 18% grey (0.18) to
130            // code 420/1023; below: linear segment maps black (0.0) to
131            // code 95/1023.
132            TransferFunction::SLog3 => { pixels.par_iter_mut().for_each(|v| { let x = *v; *v = if x >= 0.01125_f32 { (420.0_f32 + 261.5_f32 * ((x + 0.01_f32) / 0.19_f32).log10()) / 1023.0_f32 } else { (x * (171.2102946929_f32 - 95.0_f32) / 0.01125_f32 + 95.0_f32) / 1023.0_f32 }; }); }
133            // Source: Panasonic V-Log/V-Gamut Reference Manual (2014).
134            TransferFunction::VLog => { pixels.par_iter_mut().for_each(|v| { let x = *v; *v = if x < 0.01 { 5.6_f32 * x + 0.125_f32 } else { 0.241514_f32 * (x + 0.00873_f32).log10() + 0.598206_f32 }; }); }
135            // Source: ARRI LogC-3 spec (2020), EI 800.
136            TransferFunction::ARRIlog3 => { pixels.par_iter_mut().for_each(|v| { let x = *v; *v = if x > 0.010591_f32 { 0.247190_f32 * (5.555556_f32 * x + 0.052272_f32).log10() + 0.385537_f32 } else { 5.367655_f32 * x + 0.092809_f32 }; }); }
137            // Source: ARRI "LogC4 Logarithmic Color Space SPECIFICATION"
138            // (Cooper & Brendel, 2022). EI-independent log encoding optimised
139            // for 12-bit ALEV4 sensors. Two-segment with a linear-to-log
140            // threshold at x = t ≈ -0.0180967. Constants a/b/c/s/t are
141            // defined in arri_logc4_constants() below; see also colour-
142            // science/colour (`log_encoding_ARRILogC4`).
143            TransferFunction::ARRIlog4 => {
144                let (a, b, c, s, t) = arri_logc4_constants();
145                pixels.par_iter_mut().for_each(|v| {
146                    let x = *v;
147                    *v = if x >= t {
148                        ((a * x + 64.0_f32).log2() - 6.0_f32) / 14.0_f32 * b + c
149                    } else {
150                        (x - t) / s
151                    };
152                });
153            }
154            // Source: Canon C-Log3 characteristics (2016). Three-segment
155            // with a negative-side log graft and a linear middle.
156            TransferFunction::CLog3 => {
157                let neg_graft_lin = (0.097465473_f32 - 0.12512219_f32) / 1.9754798_f32;
158                let pos_graft_lin = (0.15277891_f32 - 0.12512219_f32) / 1.9754798_f32;
159                pixels.par_iter_mut().for_each(|v| {
160                    let x = *v;
161                    *v = if x < neg_graft_lin { -0.36726845_f32 * ((-x * 14.98325_f32 + 1.0_f32).max(1e-10_f32)).log10() + 0.12783901_f32 }
162                         else if x <= pos_graft_lin { 1.9754798_f32 * x + 0.12512219_f32 }
163                         else { 0.36726845_f32 * (x * 14.98325_f32 + 1.0_f32).log10() + 0.12240537_f32 };
164                });
165            }
166            // Source: Fujifilm F-Log2 Data Sheet (2021).
167            TransferFunction::FLog2 => { pixels.par_iter_mut().for_each(|v| { let x = *v; *v = if x >= 0.000889_f32 { 0.245281_f32 * (5.555556_f32 * x + 0.064829_f32).log10() + 0.384316_f32 } else { 8.799461_f32 * x + 0.092864_f32 }; }); }
168            // Source: Apple "Apple Log Profile White Paper" (Sept 2023).
169            TransferFunction::AppleLog | TransferFunction::AppleLog2 => {
170                pixels.par_iter_mut().for_each(|v| {
171                    let x = *v;
172                    const R0: f32 = -0.05641088; const RT: f32 = 0.01; const C: f32 = 47.28711236;
173                    const BETA: f32 = 0.00964052; const GAMMA: f32 = 0.08550479; const DELTA: f32 = 0.69336945;
174                    *v = if x < R0 { 0.0 } else if x < RT { C * (x - R0) * (x - R0) } else { GAMMA * (x + BETA).log2() + DELTA };
175                });
176            }
177            // Source: AMPAS ACEScc specification (TB-2022-002).
178            TransferFunction::ACESCCT => { pixels.par_iter_mut().for_each(|v| { let x = *v; *v = if x > 0.0078125_f32 { (x.log2() + 9.72_f32) / 17.52_f32 } else { 10.5402377416545_f32 * x + 0.0729055341958355_f32 }; }); }
179            // Source: ITU-R BT.2100-2 ST.2084 PQ. Input clamped to ≥0 to
180            // prevent NaN from negative values entering the power function.
181            TransferFunction::PQ => { pixels.par_iter_mut().for_each(|v| { let x = (*v).max(0.0_f32); let x_m1 = x.powf(0.1593017578125_f32); *v = ((0.8359375_f32 + 18.8515625_f32 * x_m1) / (1.0_f32 + 18.6875_f32 * x_m1)).powf(78.84375_f32); }); }
182            // Source: ITU-R BT.2100-2 HLG OETF. Input clamped to ≥0 to
183            // prevent NaN from negative values entering sqrt/ln.
184            // Knee at L = 1/12; below the knee V = sqrt(3L), above
185            // V = a*ln(12L - b) + c with a=0.17883277, b=0.28466892, c=0.55991073.
186            TransferFunction::HLG => { pixels.par_iter_mut().for_each(|v| { let x = (*v).max(0.0_f32); *v = if x < (1.0_f32 / 12.0_f32) { (3.0_f32 * x).sqrt() } else { 0.17883277_f32 * (12.0_f32 * x - 0.28466892_f32).ln() + 0.55991073_f32 }; }); }
187            // Source: Blackmagic DaVinci YRGB Intermediate white paper.
188            TransferFunction::DaVinciIntermediate => { pixels.par_iter_mut().for_each(|v| { let x = *v; *v = if x <= 0.00262409_f32 { x * 10.44426855_f32 } else { 0.07329248_f32 * ((x + 0.0075_f32).log2() + 7.0_f32) }; }); }
189            // Display gamma 1/2.4. Not a log curve; for 8-bit preview
190            // only — production use should always pick a real OETF.
191            TransferFunction::Gamma24 => { pixels.par_iter_mut().for_each(|v| { *v = v.max(0.0).powf(1.0 / 2.4); }); }
192        }
193    }
194
195    pub fn all() -> &'static [TransferFunction] {
196        // Alphabetical order for deterministic, pleasing cycle order.
197        &[TransferFunction::ACESCCT, TransferFunction::ARRIlog3, TransferFunction::ARRIlog4,
198          TransferFunction::AppleLog, TransferFunction::AppleLog2, TransferFunction::CLog3,
199          TransferFunction::DaVinciIntermediate, TransferFunction::FLog2,
200          TransferFunction::Gamma24, TransferFunction::HLG, TransferFunction::Linear,
201          TransferFunction::PQ, TransferFunction::Rec709, TransferFunction::SLog3,
202          TransferFunction::VLog]
203    }
204    pub fn next(self) -> Self { let all = Self::all(); let pos = all.iter().position(|&x| x == self).unwrap_or(0); all[(pos + 1) % all.len()] }
205    pub fn prev(self) -> Self { let all = Self::all(); let pos = all.iter().position(|&x| x == self).unwrap_or(0); all[(pos + all.len() - 1) % all.len()] }
206    pub fn is_log_bypass(&self) -> bool { !matches!(self, TransferFunction::Linear | TransferFunction::Rec709 | TransferFunction::Gamma24) }
207    pub fn requires_10bit(&self) -> bool { !matches!(self, TransferFunction::Linear | TransferFunction::Rec709 | TransferFunction::Gamma24) }
208}
209
210#[inline] pub fn rec709_oetf(x: f32) -> f32 { if x < 0.018 { 4.5 * x } else { 1.099 * x.powf(0.45) - 0.099 } }
211#[inline] pub fn rec709_eotf(x: f32) -> f32 { if x < 0.0812429 { x / 4.5 } else { ((x + 0.099) / 1.099).powf(1.0 / 0.45) } }
212
213/// ARRI LogC4 constants (a, b, c, s, t) from the 2022 LogC4 specification.
214///
215/// Derivation (Cooper & Brendel 2022, §4.1.1):
216///   a = (2^18 - 16) / 117.45
217///   b = (1023 - 95) / 1023
218///   c = 95 / 1023
219///   s = (7 · ln 2 · 2^(7 - 14·c/b)) / (a · b)
220///   t = (2^(14·(-c/b) + 6) - 64) / a
221///
222/// Cross-checked against colour-science/colour
223/// `colour.models.rgb.transfer_functions.arri` and antlerpost.com/colour-spaces/LogC4.
224pub fn arri_logc4_constants() -> (f32, f32, f32, f32, f32) {
225    let a: f32 = ((1u32 << 18) as f32 - 16.0) / 117.45;
226    let b: f32 = (1023.0 - 95.0) / 1023.0;
227    let c: f32 = 95.0 / 1023.0;
228    let s: f32 = (7.0 * std::f32::consts::LN_2 * (7.0 - 14.0 * c / b).exp2()) / (a * b);
229    let t: f32 = ((14.0 * (-c / b) + 6.0).exp2() - 64.0) / a;
230    (a, b, c, s, t)
231}
232
233/// ARRI LogC4 scene-linear → normalized log encoding (E_scene → E').
234/// Reference: ARRI "LogC4 Encoding Function" (Cooper & Brendel, 2022).
235#[inline]
236pub fn arri_logc4_oetf(x: f32) -> f32 {
237    let (a, b, c, s, t) = arri_logc4_constants();
238    if x >= t {
239        ((a * x + 64.0).log2() - 6.0) / 14.0 * b + c
240    } else {
241        (x - t) / s
242    }
243}
244
245/// ARRI LogC4 normalized log → scene-linear decoding (E' → E_scene).
246/// Reference: ARRI "LogC4 Decoding Function" (Cooper & Brendel, 2022).
247#[inline]
248pub fn arri_logc4_eotf(y: f32) -> f32 {
249    let (a, b, c, s, t) = arri_logc4_constants();
250    if y >= 0.0 {
251        ((14.0 * ((y - c) / b) + 6.0).exp2() - 64.0) / a
252    } else {
253        y * s + t
254    }
255}
256#[inline] pub fn apply_ccm(r: f32, g: f32, b: f32, ccm: &[f32; 9]) -> [f32; 3] { [r * ccm[0] + g * ccm[1] + b * ccm[2], r * ccm[3] + g * ccm[4] + b * ccm[5], r * ccm[6] + g * ccm[7] + b * ccm[8]] }
257pub fn identity_ccm() -> [f32; 9] { [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0] }
258
259pub fn invert_3x3(m: &[f32; 9]) -> [f32; 9] {
260    let det = m[0] * (m[4] * m[8] - m[5] * m[7]) - m[1] * (m[3] * m[8] - m[5] * m[6]) + m[2] * (m[3] * m[7] - m[4] * m[6]);
261    let inv_det = 1.0 / det;
262    [
263        (m[4] * m[8] - m[5] * m[7]) * inv_det, (m[2] * m[7] - m[1] * m[8]) * inv_det, (m[1] * m[5] - m[2] * m[4]) * inv_det,
264        (m[5] * m[6] - m[3] * m[8]) * inv_det, (m[0] * m[8] - m[2] * m[6]) * inv_det, (m[2] * m[3] - m[0] * m[5]) * inv_det,
265        (m[3] * m[7] - m[4] * m[6]) * inv_det, (m[1] * m[6] - m[0] * m[7]) * inv_det, (m[0] * m[4] - m[1] * m[3]) * inv_det,
266    ]
267}
268
269pub fn mat_mul_3x3(a: &[f32; 9], b: &[f32; 9]) -> [f32; 9] {
270    let mut out = [0.0; 9];
271    for i in 0..3 { for j in 0..3 { out[i * 3 + j] = a[i * 3] * b[j] + a[i * 3 + 1] * b[3 + j] + a[i * 3 + 2] * b[6 + j]; } }
272    out
273}
274
275pub fn camera_to_rec709_matrix(color_matrix: &[f32; 9]) -> [f32; 9] {
276    let cam_to_xyz = detect_camera_to_xyz(color_matrix);
277    let d50_to_d65 = [0.9555, -0.0230, 0.0633, -0.0283, 1.0099, 0.0210, 0.0123, -0.0205, 1.3300];
278    let cam_to_xyz_d65 = mat_mul_3x3(&d50_to_d65, &cam_to_xyz);
279    mat_mul_3x3(&xyz_to_rec709(), &cam_to_xyz_d65)
280}
281
282pub fn rec709_to_xyz() -> [f32; 9] { [0.4124564, 0.3575761, 0.1804375, 0.2126729, 0.7151522, 0.0721750, 0.0193339, 0.1191920, 0.9503041] }
283
284pub(crate) const MCAT16: [f32; 9] = [0.401288, 0.650173, -0.051461, -0.250268, 1.204414, 0.045854, -0.002079, 0.048952, 0.953127];
285pub(crate) const MCAT16_INV: [f32; 9] = [1.86206786, -1.01125463, 0.14918678, 0.38752654, 0.62144744, -0.00897398, -0.01584150, -0.03412294, 1.04996444];
286pub const D50_XYZ: [f32; 3] = [0.96422, 1.0, 0.82521];
287pub const D65_XYZ: [f32; 3] = [0.95047, 1.0, 1.08883];
288
289pub fn xyz_from_chromaticities(x: f32, y: f32) -> [f32; 3] { let z = 1.0 - x - y; [x / y, 1.0, z / y] }
290
291pub fn cat16_adapt(xyz: &[f32; 3], src_white: &[f32; 3], dst_white: &[f32; 3]) -> [f32; 3] {
292    let [l_s, m_s, s_s] = mat_mul_vec3(&MCAT16, src_white);
293    let [l_d, m_d, s_d] = mat_mul_vec3(&MCAT16, dst_white);
294    let lms = mat_mul_vec3(&MCAT16, xyz);
295    let adapted = [lms[0] * (l_d / l_s), lms[1] * (m_d / m_s), lms[2] * (s_d / s_s)];
296    mat_mul_vec3(&MCAT16_INV, &adapted)
297}
298
299pub fn build_cat16_output_matrix(cam_to_xyz: &[f32; 9], scene_white_xyz: &[f32; 3], dst_white: &[f32; 3], xyz_to_output: &[f32; 9]) -> [f32; 9] {
300    let [l_s, m_s, s_s] = mat_mul_vec3(&MCAT16, scene_white_xyz);
301    let [l_d, m_d, s_d] = mat_mul_vec3(&MCAT16, dst_white);
302    let r_l = l_d / l_s; let r_m = m_d / m_s; let r_s = s_d / s_s;
303    let rgb_to_lms = mat_mul_3x3(&MCAT16, cam_to_xyz);
304    let rgb_to_adapted = [
305        rgb_to_lms[0] * r_l, rgb_to_lms[1] * r_l, rgb_to_lms[2] * r_l,
306        rgb_to_lms[3] * r_m, rgb_to_lms[4] * r_m, rgb_to_lms[5] * r_m,
307        rgb_to_lms[6] * r_s, rgb_to_lms[7] * r_s, rgb_to_lms[8] * r_s,
308    ];
309    let rgb_to_xyz = mat_mul_3x3(&MCAT16_INV, &rgb_to_adapted);
310    mat_mul_3x3(xyz_to_output, &rgb_to_xyz)
311}
312
313#[inline]
314pub fn mat_mul_vec3(m: &[f32; 9], v: &[f32; 3]) -> [f32; 3] {
315    [m[0] * v[0] + m[1] * v[1] + m[2] * v[2], m[3] * v[0] + m[4] * v[1] + m[5] * v[2], m[6] * v[0] + m[7] * v[1] + m[8] * v[2]]
316}
317
318/// Bradford cone-response matrix
319pub(crate) const BRADFORD: [f32; 9] = [
320    0.8951000,  0.2664000, -0.1614000,
321   -0.7502000,  1.7135000,  0.0367000,
322    0.0389000, -0.0685000,  1.0296000,
323];
324
325/// Inverse Bradford matrix
326pub(crate) const BRADFORD_INV: [f32; 9] = [
327    0.9869929, -0.1470543,  0.1599627,
328    0.4323053,  0.5183603,  0.0492912,
329   -0.0085287,  0.0400428,  0.9684867,
330];
331
332/// Build a fused Bradford adaptation matrix: src_white → dst_white
333pub fn build_bradford_matrix(src_white: &[f32; 3], dst_white: &[f32; 3]) -> [f32; 9] {
334    let [rho_s, gamma_s, beta_s] = mat_mul_vec3(&BRADFORD, src_white);
335    let [rho_d, gamma_d, beta_d] = mat_mul_vec3(&BRADFORD, dst_white);
336
337    let scale = [
338        rho_d / rho_s, 0.0, 0.0,
339        0.0, gamma_d / gamma_s, 0.0,
340        0.0, 0.0, beta_d / beta_s,
341    ];
342
343    let temp = mat_mul_3x3(&scale, &BRADFORD);
344    mat_mul_3x3(&BRADFORD_INV, &temp)
345}
346
347/// DNG 1.4 specification:
348///   * `ColorMatrix1` is a 3x3 matrix that maps the camera's native color
349///     values to CIE XYZ (D50 / 2° observer). It is the FORWARD matrix.
350///   * `ForwardMatrix1` is a 3x3 matrix that maps XYZ (D50) to the camera's
351///     native color values. It is the INVERSE direction relative to the
352///     camera→XYZ transform we need for the rendering pipeline.
353///   * `CalibrationMatrix1` is applied in camera-native space BEFORE
354///     `ColorMatrix1`, so the effective transform is
355///     `ColorMatrix1 * CalibrationMatrix1 * camera_native`.
356///
357/// MCRAW embeds these matrices verbatim. We don't know whether a given file
358/// stores them row-major or column-major, nor whether any tooling has
359/// pre-inverted the direction, so we evaluate the four possible orientations
360/// and pick the one whose implied scene white point best matches D50.
361pub fn detect_camera_to_xyz(m: &[f32; 9]) -> [f32; 9] {
362    let d50 = D50_XYZ;
363    let transposed = [
364        m[0], m[3], m[6],
365        m[1], m[4], m[7],
366        m[2], m[5], m[8],
367    ];
368    let inv = invert_3x3(m);
369    let inv_t = invert_3x3(&transposed);
370
371    let candidates: [[f32; 9]; 4] = [*m, transposed, inv, inv_t];
372
373    // For each candidate compute the implied white point: a forward
374    // Camera→XYZ matrix sends (1,1,1) to the white in XYZ, so the
375    // row-sums of that matrix equal that white point. An XYZ→Camera
376    // matrix has its row-sums equal to the row-basis sums (not the white
377    // point) and is rejected by the distance check below.
378    let mut best = *m;
379    let mut best_dist = f32::MAX;
380    for c in &candidates {
381        let w = [c[0] + c[1] + c[2], c[3] + c[4] + c[5], c[6] + c[7] + c[8]];
382        let dx = w[0] - d50[0];
383        let dy = w[1] - d50[1];
384        let dz = w[2] - d50[2];
385        let dist = dx * dx + dy * dy + dz * dz;
386        if dist < best_dist {
387            best_dist = dist;
388            best = *c;
389        }
390    }
391    tracing::debug!(
392        "detect_camera_to_xyz: white=[{:.3},{:.3},{:.3}] dist={:.4}",
393        best[0] + best[1] + best[2],
394        best[3] + best[4] + best[5],
395        best[6] + best[7] + best[8],
396        best_dist.sqrt()
397    );
398    best
399}
400
401/// Build a camera→XYZ matrix from a DNG `ColorMatrix1` and (optional)
402/// `CalibrationMatrix1`. Orientation is auto-detected — see
403/// [`detect_camera_to_xyz`].
404pub fn camera_to_xyz_matrix(color_matrix: &[f32; 9], calibration_matrix: Option<&[f32; 9]>) -> [f32; 9] {
405    let cam_to_xyz = detect_camera_to_xyz(color_matrix);
406    match calibration_matrix {
407        // Calibration is applied in camera-native space BEFORE the
408        // forward XYZ transform, so the product order is
409        // `ColorMatrix1 * CalibrationMatrix1`.
410        Some(cal) => mat_mul_3x3(&cam_to_xyz, cal),
411        None => cam_to_xyz,
412    }
413}
414
415/// Invert a forward matrix to recover Camera→XYZ when only a
416/// `ForwardMatrix1` (XYZ→Camera) is available.
417pub fn forward_to_camera_xyz(forward_matrix: &[f32; 9]) -> [f32; 9] {
418    detect_camera_to_xyz(forward_matrix)
419}
420
421/// Build a fused Camera→Rec709 CCM for the preview thumbnail path.
422///
423/// Mirrors the export pipeline's CCM construction (pipeline.rs:128-204):
424///   1. Prefer ForwardMatrix1+2 (averaged) when available — already D50-adapted
425///   2. Fall back to ColorMatrix1+2 + calibration + chromatic adaptation
426///   3. Detect matrix orientation automatically (D50 white-point row-sum check)
427///   4. Bradford-adapt from reference illuminant to D65
428///   5. Fuse with Rec709 primaries
429///
430/// This fixes the green/pink tint on older MOTION files where the raw
431/// `ColorMatrix1` was applied without orientation detection or D50→D65
432/// chromatic adaptation.
433pub fn build_preview_ccm(
434    color_matrix: Option<&[f64; 9]>,
435    forward_matrix1: Option<&[f64; 9]>,
436    forward_matrix2: Option<&[f64; 9]>,
437    color_matrix2: Option<&[f64; 9]>,
438    calibration_matrix1: Option<&[f64; 9]>,
439) -> [f32; 9] {
440    let cm1_f32 = color_matrix.map(|m| [m[0] as f32, m[1] as f32, m[2] as f32, m[3] as f32, m[4] as f32, m[5] as f32, m[6] as f32, m[7] as f32, m[8] as f32]);
441    let cm2_f32 = color_matrix2.map(|m| [m[0] as f32, m[1] as f32, m[2] as f32, m[3] as f32, m[4] as f32, m[5] as f32, m[6] as f32, m[7] as f32, m[8] as f32]);
442    let fm1_f32 = forward_matrix1.map(|m| [m[0] as f32, m[1] as f32, m[2] as f32, m[3] as f32, m[4] as f32, m[5] as f32, m[6] as f32, m[7] as f32, m[8] as f32]);
443    let fm2_f32 = forward_matrix2.map(|m| [m[0] as f32, m[1] as f32, m[2] as f32, m[3] as f32, m[4] as f32, m[5] as f32, m[6] as f32, m[7] as f32, m[8] as f32]);
444    let cal1_f32 = calibration_matrix1.map(|m| [m[0] as f32, m[1] as f32, m[2] as f32, m[3] as f32, m[4] as f32, m[5] as f32, m[6] as f32, m[7] as f32, m[8] as f32]);
445
446    let cam_to_xyz: [f32; 9] = if let (Some(ref fm1), Some(ref fm2)) = (fm1_f32, fm2_f32) {
447        let fm_avg = interpolate_matrix(fm1, fm2, 0.5);
448        let rs = [fm_avg[0] + fm_avg[1] + fm_avg[2], fm_avg[3] + fm_avg[4] + fm_avg[5], fm_avg[6] + fm_avg[7] + fm_avg[8]];
449        let d = (rs[0] - D50_XYZ[0]).powi(2) + (rs[1] - D50_XYZ[1]).powi(2) + (rs[2] - D50_XYZ[2]).powi(2);
450        if d < 0.05 { fm_avg } else { detect_camera_to_xyz(&fm_avg) }
451    } else if let Some(ref fm1) = fm1_f32 {
452        let rs = [fm1[0] + fm1[1] + fm1[2], fm1[3] + fm1[4] + fm1[5], fm1[6] + fm1[7] + fm1[8]];
453        let d = (rs[0] - D50_XYZ[0]).powi(2) + (rs[1] - D50_XYZ[1]).powi(2) + (rs[2] - D50_XYZ[2]).powi(2);
454        if d < 0.05 { *fm1 } else { detect_camera_to_xyz(fm1) }
455    } else if let Some(ref cm1) = cm1_f32 {
456        let cal = cal1_f32;
457        match cm2_f32 {
458            Some(ref cm2) => {
459                let cm_avg = interpolate_matrix(cm1, cm2, 0.5);
460                camera_to_xyz_matrix(&cm_avg, cal.as_ref())
461            }
462            None => camera_to_xyz_matrix(cm1, cal.as_ref()),
463        }
464    } else {
465        identity_ccm()
466    };
467
468    // Determine reference illuminant from the selected matrix
469    let rs = [cam_to_xyz[0] + cam_to_xyz[1] + cam_to_xyz[2], cam_to_xyz[3] + cam_to_xyz[4] + cam_to_xyz[5], cam_to_xyz[6] + cam_to_xyz[7] + cam_to_xyz[8]];
470    let cam_illuminant_xyz = if fm1_f32.is_some() {
471        D50_XYZ
472    } else {
473        let l = rs[0].max(rs[1]).max(rs[2]);
474        if l < 0.1 || l > 5.0 { D50_XYZ } else { rs }
475    };
476
477    let bradford_static = build_bradford_matrix(&cam_illuminant_xyz, &D65_XYZ);
478    let cam_to_xyz_d65 = mat_mul_3x3(&bradford_static, &cam_to_xyz);
479    mat_mul_3x3(&xyz_to_rec709(), &cam_to_xyz_d65)
480}
481
482pub fn interpolate_matrix(a: &[f32; 9], b: &[f32; 9], t: f32) -> [f32; 9] {
483    let s = 1.0 - t; let mut out = [0.0; 9];
484    for i in 0..9 { out[i] = a[i] * s + b[i] * t; }
485    out
486}
487
488pub fn xyz_to_rec709() -> [f32; 9] { [3.2404542, -1.5371385, -0.4985354, -0.9689294, 1.8767608, 0.0415560, 0.0556434, -0.2040259, 1.0572252] }
489
490pub fn xyz_to_rgb_from_primaries(xr: f32, yr: f32, xg: f32, yg: f32, xb: f32, yb: f32, xw: f32, yw: f32) -> [f32; 9] {
491    let xr_z = (1.0 - xr - yr) / yr; let xg_z = (1.0 - xg - yg) / yg; let xb_z = (1.0 - xb - yb) / yb;
492    let m = [xr / yr, xg / yg, xb / yb, 1.0, 1.0, 1.0, xr_z, xg_z, xb_z];
493    let wx = xw / yw; let wy = 1.0; let wz = (1.0 - xw - yw) / yw;
494    let det_m = m[0] * (m[4] * m[8] - m[5] * m[7]) - m[1] * (m[3] * m[8] - m[5] * m[6]) + m[2] * (m[3] * m[7] - m[4] * m[6]);
495    let inv_det = 1.0 / det_m;
496    let inv_m = [
497        (m[4] * m[8] - m[5] * m[7]) * inv_det, (m[2] * m[7] - m[1] * m[8]) * inv_det, (m[1] * m[5] - m[2] * m[4]) * inv_det,
498        (m[5] * m[6] - m[3] * m[8]) * inv_det, (m[0] * m[8] - m[2] * m[6]) * inv_det, (m[2] * m[3] - m[0] * m[5]) * inv_det,
499        (m[3] * m[7] - m[4] * m[6]) * inv_det, (m[1] * m[6] - m[0] * m[7]) * inv_det, (m[0] * m[4] - m[1] * m[3]) * inv_det,
500    ];
501    let sr = inv_m[0] * wx + inv_m[1] * wy + inv_m[2] * wz;
502    let sg = inv_m[3] * wx + inv_m[4] * wy + inv_m[5] * wz;
503    let sb = inv_m[6] * wx + inv_m[7] * wy + inv_m[8] * wz;
504    let rgb_to_xyz = [m[0] * sr, m[1] * sg, m[2] * sb, m[3] * sr, m[4] * sg, m[5] * sb, m[6] * sr, m[7] * sg, m[8] * sb];
505    invert_3x3(&rgb_to_xyz)
506}
507
508pub struct BilinearDemosaic { pattern: BayerPattern }
509impl BilinearDemosaic {
510    pub fn new(pattern: BayerPattern) -> Self { BilinearDemosaic { pattern } }
511    
512    fn get_pixel(&self, bayer: &[u16], stride_width: u32, x: i32, y: i32) -> f64 {
513        if x < 0 || y < 0 || x >= stride_width as i32 { return 0.0; }
514        let idx = (y as usize) * (stride_width as usize) + (x as usize);
515        if idx >= bayer.len() { return 0.0; }
516        bayer[idx] as f64
517    }
518
519    fn is_red_site(&self, x: i32, y: i32, pattern: BayerPattern) -> bool {
520        match pattern {
521            BayerPattern::RGGB => x % 2 == 0 && y % 2 == 0,
522            BayerPattern::BGGR => x % 2 == 1 && y % 2 == 1,
523            BayerPattern::GRBG => x % 2 == 1 && y % 2 == 0,
524            BayerPattern::GBRG => x % 2 == 0 && y % 2 == 1,
525            _ => false,
526        }
527    }
528
529    fn is_blue_site(&self, x: i32, y: i32, pattern: BayerPattern) -> bool {
530        match pattern {
531            BayerPattern::RGGB => x % 2 == 1 && y % 2 == 1,
532            BayerPattern::BGGR => x % 2 == 0 && y % 2 == 0,
533            BayerPattern::GRBG => x % 2 == 0 && y % 2 == 1,
534            BayerPattern::GBRG => x % 2 == 1 && y % 2 == 0,
535            _ => false,
536        }
537    }
538
539    fn interp_green_at_red(&self, bayer: &[u16], stride: u32, _height: u32, x: i32, y: i32, pattern: BayerPattern) -> f64 {
540        let mut sum = 0.0; let mut count = 0.0;
541        let positions = [(0, -1), (0, 1), (-1, 0), (1, 0)];
542        for (dx, dy) in positions.iter() {
543            let px = x + dx; let py = y + dy;
544            if self.is_green_site(px, py, pattern) { sum += self.get_pixel(bayer, stride, px, py); count += 1.0; }
545        }
546        if count > 0.0 { sum / count } else { self.get_pixel(bayer, stride, x, y) }
547    }
548
549    fn interp_green_at_blue(&self, bayer: &[u16], stride: u32, height: u32, x: i32, y: i32, pattern: BayerPattern) -> f64 {
550        self.interp_green_at_red(bayer, stride, height, x, y, pattern)
551    }
552
553    fn interp_blue_at_red(&self, bayer: &[u16], stride: u32, _height: u32, x: i32, y: i32, pattern: BayerPattern) -> f64 {
554        let mut sum = 0.0; let mut count = 0.0;
555        let positions = [(-1, -1), (1, -1), (-1, 1), (1, 1)];
556        for (dx, dy) in positions.iter() {
557            let px = x + dx; let py = y + dy;
558            if self.is_blue_site(px, py, pattern) { sum += self.get_pixel(bayer, stride, px, py); count += 1.0; }
559        }
560        if count > 0.0 { sum / count } else { self.get_pixel(bayer, stride, x, y) }
561    }
562
563    fn interp_red_at_blue(&self, bayer: &[u16], stride: u32, _height: u32, x: i32, y: i32, pattern: BayerPattern) -> f64 {
564        let mut sum = 0.0; let mut count = 0.0;
565        let positions = [(-1, -1), (1, -1), (-1, 1), (1, 1)];
566        for (dx, dy) in positions.iter() {
567            let px = x + dx; let py = y + dy;
568            if self.is_red_site(px, py, pattern) { sum += self.get_pixel(bayer, stride, px, py); count += 1.0; }
569        }
570        if count > 0.0 { sum / count } else { self.get_pixel(bayer, stride, x, y) }
571    }
572    
573    fn is_green_site(&self, x: i32, y: i32, pattern: BayerPattern) -> bool {
574        !self.is_red_site(x, y, pattern) && !self.is_blue_site(x, y, pattern)
575    }
576
577    pub fn process_par(&self, bayer: &[u16], stride_width: u32, offset_x: u32, offset_y: u32, active_width: u32, active_height: u32, pattern: &BayerPattern) -> Result<Vec<f32>> {
578        let stride = stride_width as usize; let ox = offset_x as i32; let oy = offset_y as i32;
579        let aw = active_width as usize; let ah = active_height as usize;
580        let min_len = (stride * (oy as usize + ah - 1) + ox as usize + aw - 1) + 1;
581        if bayer.len() < min_len { anyhow::bail!("Bayer data too short"); }
582        let mut rgb = vec![0.0f32; aw * ah * 3]; let pat = *pattern; let row_len = aw * 3;
583        rgb.par_chunks_exact_mut(row_len).enumerate().for_each(|(sy, row)| {
584            let y = sy as i32 + oy;
585            for sx in 0..aw {
586                let x = sx as i32 + ox;
587                let is_red = self.is_red_site(x, y, pat); let is_blue = self.is_blue_site(x, y, pat);
588                let (r, g, b) = if is_red {
589                    (self.get_pixel(bayer, stride_width, x, y), self.interp_green_at_red(bayer, stride_width, active_height, x, y, pat), self.interp_blue_at_red(bayer, stride_width, active_height, x, y, pat))
590                } else if is_blue {
591                    (self.interp_red_at_blue(bayer, stride_width, active_height, x, y, pat), self.interp_green_at_blue(bayer, stride_width, active_height, x, y, pat), self.get_pixel(bayer, stride_width, x, y))
592                } else {
593                    // FIXED: GBRG top-green logic
594                    let is_top_green = match pat {
595                        BayerPattern::RGGB | BayerPattern::BGGR => y % 2 == 0,
596                        BayerPattern::GRBG => y % 2 == 0,
597                        BayerPattern::GBRG => y % 2 == 0, 
598                        _ => y % 2 == 0,
599                    };
600                    if is_top_green {
601                        (self.interp_red_at_blue(bayer, stride_width, active_height, x + 1, y, pat), self.get_pixel(bayer, stride_width, x, y), self.interp_blue_at_red(bayer, stride_width, active_height, x - 1, y, pat))
602                    } else {
603                        (self.interp_red_at_blue(bayer, stride_width, active_height, x - 1, y, pat), self.get_pixel(bayer, stride_width, x, y), self.interp_blue_at_red(bayer, stride_width, active_height, x + 1, y, pat))
604                    }
605                };
606                let base = sx * 3; row[base] = r as f32; row[base + 1] = g as f32; row[base + 2] = b as f32;
607            }
608        });
609        Ok(rgb)
610    }
611
612    pub fn process_par_into(&self, bayer: &[u16], stride_width: u32, offset_x: u32, offset_y: u32, active_width: u32, active_height: u32, pattern: &BayerPattern, output: &mut [f32]) -> Result<()> {
613        let stride = stride_width as usize; let ox = offset_x as i32; let oy = offset_y as i32;
614        let aw = active_width as usize; let ah = active_height as usize;
615        let min_len = (stride * (oy as usize + ah - 1) + ox as usize + aw - 1) + 1;
616        if bayer.len() < min_len { anyhow::bail!("Bayer data too short"); }
617        if output.len() < aw * ah * 3 { anyhow::bail!("Output buffer too short"); }
618        let pat = *pattern; let row_len = aw * 3;
619        output.par_chunks_exact_mut(row_len).enumerate().for_each(|(sy, row)| {
620            let y = sy as i32 + oy;
621            for sx in 0..aw {
622                let x = sx as i32 + ox;
623                let is_red = self.is_red_site(x, y, pat); let is_blue = self.is_blue_site(x, y, pat);
624                let (r, g, b) = if is_red {
625                    (self.get_pixel(bayer, stride_width, x, y), self.interp_green_at_red(bayer, stride_width, active_height, x, y, pat), self.interp_blue_at_red(bayer, stride_width, active_height, x, y, pat))
626                } else if is_blue {
627                    (self.interp_red_at_blue(bayer, stride_width, active_height, x, y, pat), self.interp_green_at_blue(bayer, stride_width, active_height, x, y, pat), self.get_pixel(bayer, stride_width, x, y))
628                } else {
629                    // FIXED: GBRG top-green logic
630                    let is_top_green = match pat {
631                        BayerPattern::RGGB | BayerPattern::BGGR => y % 2 == 0,
632                        BayerPattern::GRBG => y % 2 == 0,
633                        BayerPattern::GBRG => y % 2 == 0,
634                        _ => y % 2 == 0,
635                    };
636                    if is_top_green {
637                        (self.interp_red_at_blue(bayer, stride_width, active_height, x + 1, y, pat), self.get_pixel(bayer, stride_width, x, y), self.interp_blue_at_red(bayer, stride_width, active_height, x - 1, y, pat))
638                    } else {
639                        (self.interp_red_at_blue(bayer, stride_width, active_height, x - 1, y, pat), self.get_pixel(bayer, stride_width, x, y), self.interp_blue_at_red(bayer, stride_width, active_height, x + 1, y, pat))
640                    }
641                };
642                let base = sx * 3; row[base] = r as f32; row[base + 1] = g as f32; row[base + 2] = b as f32;
643            }
644        });
645        Ok(())
646    }
647}
648
649impl Demosaic for BilinearDemosaic {
650    fn process(&self, bayer: &[u16], stride_width: u32, offset_x: u32, offset_y: u32, active_width: u32, active_height: u32, pattern: &BayerPattern) -> Result<Vec<f32>> {
651        let stride = stride_width as usize; let ox = offset_x as i32; let oy = offset_y as i32;
652        let aw = active_width as usize; let ah = active_height as usize;
653        let min_len = (stride * (oy as usize + ah - 1) + ox as usize + aw - 1) + 1;
654        if bayer.len() < min_len { anyhow::bail!("Bayer data too short"); }
655        let mut rgb = Vec::with_capacity(aw * ah * 3); let pat = *pattern;
656        for sy in 0..ah as i32 {
657            for sx in 0..aw as i32 {
658                let x = sx + ox; let y = sy + oy;
659                let is_red = self.is_red_site(x, y, pat); let is_blue = self.is_blue_site(x, y, pat);
660                let (r, g, b) = if is_red {
661                    (self.get_pixel(bayer, stride_width, x, y), self.interp_green_at_red(bayer, stride_width, active_height, x, y, pat), self.interp_blue_at_red(bayer, stride_width, active_height, x, y, pat))
662                } else if is_blue {
663                    (self.interp_red_at_blue(bayer, stride_width, active_height, x, y, pat), self.interp_green_at_blue(bayer, stride_width, active_height, x, y, pat), self.get_pixel(bayer, stride_width, x, y))
664                } else {
665                    // FIXED: GBRG top-green logic
666                    let is_top_green = match pat {
667                        BayerPattern::RGGB | BayerPattern::BGGR => y % 2 == 0,
668                        BayerPattern::GRBG => y % 2 == 0,
669                        BayerPattern::GBRG => y % 2 == 0,
670                        _ => y % 2 == 0,
671                    };
672                    if is_top_green {
673                        (self.interp_red_at_blue(bayer, stride_width, active_height, x + 1, y, pat), self.get_pixel(bayer, stride_width, x, y), self.interp_blue_at_red(bayer, stride_width, active_height, x - 1, y, pat))
674                    } else {
675                        (self.interp_red_at_blue(bayer, stride_width, active_height, x - 1, y, pat), self.get_pixel(bayer, stride_width, x, y), self.interp_blue_at_red(bayer, stride_width, active_height, x + 1, y, pat))
676                    }
677                };
678                rgb.push(r as f32); rgb.push(g as f32); rgb.push(b as f32);
679            }
680        }
681        Ok(rgb)
682    }
683}
684
685pub struct CcmColorSpaceConverter;
686impl CcmColorSpaceConverter { pub fn new() -> Self { CcmColorSpaceConverter } }
687impl Default for CcmColorSpaceConverter { fn default() -> Self { Self::new() } }
688impl ColorSpaceConverter for CcmColorSpaceConverter {
689    fn process(&self, pixels: &mut [f32], ccm: &[f32; 9]) {
690        for chunk in pixels.chunks_exact_mut(3) {
691            let [r_out, g_out, b_out] = apply_ccm(chunk[0], chunk[1], chunk[2], ccm);
692            chunk[0] = r_out.max(0.0).min(1.0); chunk[1] = g_out.max(0.0).min(1.0); chunk[2] = b_out.max(0.0).min(1.0);
693        }
694    }
695}
696
697pub struct Rec709TransferFunction;
698impl Rec709TransferFunction { pub fn new() -> Self { Rec709TransferFunction } }
699impl TransferFunctionProcessor for Rec709TransferFunction {
700    fn process(&self, pixels: &mut [f32]) { pixels.par_iter_mut().for_each(|v| { *v = rec709_oetf(*v).min(1.0).max(0.0); }); }
701}
702
703pub struct LinearTransferFunction;
704impl LinearTransferFunction { pub fn new() -> Self { LinearTransferFunction } }
705impl TransferFunctionProcessor for LinearTransferFunction { fn process(&self, _pixels: &mut [f32]) {} }
706
707pub struct AgxKrakenPipeline { demosaic: BilinearDemosaic, agx: AgxPipeline, output_gamma: f32, enable_tonemap: bool }
708impl AgxKrakenPipeline {
709    pub fn new(pattern: BayerPattern) -> Self {
710        let config = ColorPipelineConfig::broadcast(); let demosaic = BilinearDemosaic::new(pattern);
711        let agx = AgxPipeline::new(config.tonemap_config.clone()); let output_gamma = config.output_gamma.gamma();
712        let enable_tonemap = config.enable_tonemapping;
713        AgxKrakenPipeline { demosaic, agx, output_gamma, enable_tonemap }
714    }
715}
716
717#[derive(Debug, Clone, Copy, PartialEq, Eq)] pub enum OutputGamma { Srgb, Bt1886, Linear }
718impl OutputGamma { pub fn gamma(&self) -> f32 { match self { OutputGamma::Srgb => 2.2, OutputGamma::Bt1886 => 2.4, OutputGamma::Linear => 1.0 } } }
719
720pub struct ColorPipelineConfig {
721    pub input_color_space: ColorSpace, pub input_transfer: TransferFunction, pub output_color_space: ColorSpace,
722    pub output_transfer: TransferFunction, pub output_gamma: OutputGamma, pub enable_tonemapping: bool, pub tonemap_config: AgxConfig,
723}
724impl Default for ColorPipelineConfig {
725    fn default() -> Self {
726        Self { input_color_space: ColorSpace::Rec709, input_transfer: TransferFunction::Linear, output_color_space: ColorSpace::Rec709, output_transfer: TransferFunction::Rec709, output_gamma: OutputGamma::Bt1886, enable_tonemapping: true, tonemap_config: AgxConfig::default() }
727    }
728}
729impl ColorPipelineConfig {
730    pub fn broadcast() -> Self {
731        let mut config = AgxConfig::default(); config.in_gamut = Gamut::Rec709; config.in_transfer = Transfer::Linear;
732        config.working_curve = Transfer::AgxLogKraken; config.out_gamut = Gamut::Rec709; config.out_transfer = OutputTransfer::Bt1886InverseEotf;
733        config.toe_power = 3.0; config.shoulder_power = 3.25; config.slope = 2.0; config.working_mid_grey = 0.606060; config.log_output = false;
734        Self { input_color_space: ColorSpace::Rec709, input_transfer: TransferFunction::Linear, output_color_space: ColorSpace::Rec709, output_transfer: TransferFunction::Rec709, output_gamma: OutputGamma::Bt1886, enable_tonemapping: true, tonemap_config: config }
735    }
736    pub fn log_output(log_space: TransferFunction, gamut: ColorSpace) -> Self {
737        let mut config = AgxConfig::default(); config.in_gamut = Gamut::Rec709; config.in_transfer = Transfer::Linear;
738        config.working_curve = Transfer::AgxLogKraken;
739        config.out_gamut = match gamut {
740            ColorSpace::Rec709 => Gamut::Rec709, ColorSpace::Rec2020 => Gamut::Rec2020,
741            ColorSpace::DciP3 | ColorSpace::DisplayP3 => Gamut::P3D65, ColorSpace::SGamut3Cine => Gamut::SGamut3Cine,
742            ColorSpace::SGamut3 => Gamut::SGamut3, ColorSpace::ARRIWideGamut3 | ColorSpace::ARRIWideGamut4 => Gamut::Awg3,
743            ColorSpace::CanonCinemaGamut => Gamut::CanonCinema, ColorSpace::ACESAP1 => Gamut::Ap1,
744            ColorSpace::FGamut | ColorSpace::PanasonicVGamut => Gamut::Rwg, ColorSpace::FGamutC => Gamut::Ap0,
745            ColorSpace::DaVinciWideGamut => Gamut::DaVinciWg, _ => Gamut::Rec709,
746        };
747        config.out_transfer = OutputTransfer::Linear; config.log_output = true;
748        Self { input_color_space: ColorSpace::Rec709, input_transfer: TransferFunction::Linear, output_color_space: gamut, output_transfer: log_space, output_gamma: OutputGamma::Linear, enable_tonemapping: false, tonemap_config: config }
749    }
750}
751
752pub fn pipeline_convert_to_u16(pixels: &[f32]) -> Vec<u16> { pixels.iter().map(|&v| (v.clamp(0.0, 1.0) * 65535.0) as u16).collect() }
753
754pub fn highlight_clip(pixels: &mut [f32], threshold: f32) {
755    let range = 1.0 - threshold; if range <= 0.0 { return; }
756    for chunk in pixels.chunks_exact_mut(3) {
757        let r = chunk[0]; let g = chunk[1]; let b = chunk[2];
758        let max_val = r.max(g).max(b);
759        if max_val > threshold {
760            let t = ((max_val - threshold) / range).min(1.0);
761            chunk[0] = r + (max_val - r) * t; chunk[1] = g + (max_val - g) * t; chunk[2] = b + (max_val - b) * t;
762        }
763    }
764}
765
766pub fn normalize_linear(pixels: &mut [f32], black_level: f64, white_level: f64) {
767    let range = if white_level > black_level { white_level - black_level } else { 1.0 }; let inv_range = 1.0 / range;
768    for v in pixels.iter_mut() { *v = ((*v as f64 - black_level) * inv_range).clamp(0.0, 1.0) as f32; }
769}
770
771pub fn normalize_linear_f32(pixels: &mut [f32], black_level: f32, white_level: f32) {
772    let range = if white_level > black_level { white_level - black_level } else { 1.0 }; let inv_range = 1.0 / range;
773    pixels.par_iter_mut().for_each(|v| { *v = (*v - black_level) * inv_range; if *v < 0.0 { *v = 0.0; } else if *v > 1.0 { *v = 1.0; } });
774}
775
776#[cfg(test)]
777mod tests {
778    use super::*;
779
780    /// The detector should pick the identity matrix as-is when given the
781    /// identity (the row-sum white is exactly D50).
782    #[test]
783    fn detect_camera_to_xyz_picks_identity_when_input_is_identity() {
784        let id = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0];
785        let out = detect_camera_to_xyz(&id);
786        for i in 0..9 {
787            assert!((out[i] - id[i]).abs() < 1e-5, "entry {} differs: {} vs {}", i, out[i], id[i]);
788        }
789    }
790
791    /// A forward Camera→XYZ matrix that maps (1,1,1)→D50 should be picked
792    /// over its inverse / transpose. The detector should not fall back to
793    /// the inverse of a forward matrix.
794    #[test]
795    fn detect_camera_to_xyz_prefers_forward_over_inverse() {
796        // Build a known forward matrix: diag(s) with a D50 row-sum.
797        // Row-sums must equal D50_XYZ. Simplest: identity (above) is the
798        // forward direction; the inverse IS also identity — so we use a
799        // non-trivial scaling. Let s = D50_XYZ (so the matrix is
800        // diag(d50)). Forward row-sum = D50. Its inverse has row-sum =
801        // (1/d50_x, 1/d50_y, 1/d50_z), which is far from D50.
802        let m = [
803            D50_XYZ[0], 0.0, 0.0,
804            0.0, D50_XYZ[1], 0.0,
805            0.0, 0.0, D50_XYZ[2],
806        ];
807        let out = detect_camera_to_xyz(&m);
808        for i in 0..9 {
809            assert!((out[i] - m[i]).abs() < 1e-5, "entry {} differs: {} vs {}", i, out[i], m[i]);
810        }
811    }
812
813    /// HLG OETF at the knee L = 1/12 should give V = 0.5 on both sides
814    /// of the branch, and the function must be monotonic.
815    #[test]
816    fn hlg_knee_is_continuous_at_one_twelfth() {
817        let below = TransferFunction::HLG.process_apply(1.0 / 12.0);
818        let above = TransferFunction::HLG.process_apply(1.0 / 12.0 + 1e-4);
819        let mid = (below + above) * 0.5;
820        assert!((below - 0.5).abs() < 1e-4, "HLG at knee: {} (want 0.5)", below);
821        assert!((above - 0.5).abs() < 5e-3, "HLG just above knee: {} (want ~0.5)", above);
822        assert!((mid - 0.5).abs() < 5e-3, "HLG mid (knee avg): {}", mid);
823        // Monotonicity sanity at three points.
824        let a = TransferFunction::HLG.process_apply(0.001);
825        let b = TransferFunction::HLG.process_apply(0.1);
826        let c = TransferFunction::HLG.process_apply(0.8);
827        assert!(a < b && b < c, "HLG must be monotonic: a={} b={} c={}", a, b, c);
828    }
829
830    /// PQ forward then inverse (the inverse function is not exported but
831    /// we can sanity-check the forward is monotone and stays in [0,1] for
832    /// inputs in [0,1]).
833    #[test]
834    fn pq_forward_is_monotone_bounded() {
835        let pf = |x: f32| {
836            let x_m1 = x.powf(0.1593017578125_f32);
837            ((0.8359375_f32 + 18.8515625_f32 * x_m1) / (1.0_f32 + 18.6875_f32 * x_m1)).powf(78.84375_f32)
838        };
839        for s in [0.0_f32, 0.01, 0.1, 0.18, 0.5, 1.0] {
840            let v = pf(s);
841            assert!(v.is_finite() && v >= 0.0 && v <= 1.0, "PQ({}) = {}", s, v);
842        }
843        // Monotonicity
844        let a = pf(0.10);
845        let b = pf(0.18);
846        let c = pf(0.50);
847        assert!(a < b && b < c, "PQ must be monotonic: a={} b={} c={}", a, b, c);
848    }
849
850    /// `build_bradford_matrix` from D65 to D65 must be the identity.
851    #[test]
852    fn bradford_identity_for_same_white() {
853        let m = build_bradford_matrix(&D65_XYZ, &D65_XYZ);
854        for i in 0..9 {
855            let expected = if i == 0 || i == 4 || i == 8 { 1.0 } else { 0.0 };
856            assert!((m[i] - expected).abs() < 1e-4, "entry {}: {} (want {})", i, m[i], expected);
857        }
858    }
859
860    /// Rec.709 OETF spot checks. The 0.018 knee and the `1.099`/`0.099`
861    /// coefficients are the only place the linear and power segments
862    /// meet. Below the knee the slope is 4.5; above the knee the
863    /// `x^0.45` form is used.
864    #[test]
865    fn rec709_oetf_at_key_points() {
866        let v_zero = TransferFunction::Rec709.process_apply(0.0);
867        let v_low  = TransferFunction::Rec709.process_apply(0.01);
868        let v_knee = TransferFunction::Rec709.process_apply(0.018);
869        let v_high = TransferFunction::Rec709.process_apply(0.5);
870        let v_one  = TransferFunction::Rec709.process_apply(1.0);
871        assert!(v_zero.abs() < 1e-6, "Rec.709 at 0 = {}", v_zero);
872        // Linear segment: V = 4.5 * 0.01 = 0.045.
873        assert!((v_low - 0.045).abs() < 1e-4, "Rec.709 at 0.01 = {}", v_low);
874        // Power segment: V = 1.099 * 0.018^0.45 - 0.099.
875        // (Linear segment would be V = 4.5*0.018 = 0.081, so the
876        // power segment value is the more diagnostic of the two.)
877        let power_at_knee = 1.099_f32 * 0.018_f32.powf(0.45) - 0.099;
878        assert!((v_knee - power_at_knee).abs() < 1e-4, "Rec.709 at 0.018 = {}", v_knee);
879        // At x=1.0, V = 1.099 - 0.099 = 1.0.
880        assert!((v_one - 1.0).abs() < 1e-4, "Rec.709 at 1.0 = {}", v_one);
881        // Monotonicity.
882        assert!(v_zero < v_low && v_low < v_knee && v_knee < v_high && v_high < v_one,
883                "Rec.709 must be monotonic");
884        // Spot v_high should land in the power branch.
885        let power_high = 1.099_f32 * 0.5_f32.powf(0.45) - 0.099;
886        assert!((v_high - power_high).abs() < 1e-4, "Rec.709 at 0.5 = {} (power={})", v_high, power_high);
887    }
888
889    /// V-Log (Panasonic) spot checks. Knee at x=0.01; below the knee
890    /// the linear slope is 5.6 (offset 0.125), above the knee the
891    /// log10 form with offset 0.00873 is used.
892    #[test]
893    fn vlog_at_key_points() {
894        let v_knee = TransferFunction::VLog.process_apply(0.01);
895        // Below the knee: 5.6 * 0.01 + 0.125 = 0.181.
896        assert!((v_knee - 0.181).abs() < 1e-4, "V-Log at knee = {} (want 0.181)", v_knee);
897        let v_one = TransferFunction::VLog.process_apply(1.0);
898        // Log branch: 0.241514 * log10(1.00873) + 0.598206.
899        let expected = 0.241514_f32 * (1.0_f32 + 0.00873_f32).log10() + 0.598206_f32;
900        assert!((v_one - expected).abs() < 1e-3, "V-Log at 1.0 = {} (want {})", v_one, expected);
901    }
902
903    /// ARRI LogC3 (EI 800) spot check. Knee at x=0.010591; below
904    /// the knee linear with slope 5.367655, above log with the
905    /// published coefficients.
906    #[test]
907    fn arri_logc3_at_key_points() {
908        let v_one = TransferFunction::ARRIlog3.process_apply(1.0);
909        let expected = 0.247190_f32 * (5.555556_f32 + 0.052272_f32).log10() + 0.385537_f32;
910        assert!((v_one - expected).abs() < 1e-3, "ARRI LogC3 at 1.0 = {} (want {})", v_one, expected);
911        let v_low = TransferFunction::ARRIlog3.process_apply(0.0);
912        // Linear segment: 5.367655 * 0 + 0.092809 = 0.092809.
913        assert!((v_low - 0.092809).abs() < 1e-4, "ARRI LogC3 at 0 = {} (want 0.092809)", v_low);
914    }
915
916    /// ARRI LogC4 spot check. Cross-checked against colour-science/colour
917    /// `log_encoding_ARRILogC4` / `log_decoding_ARRILogC4`. Encoding of
918    /// 0.18 (18% grey) must be ≈ 0.2783958, and the round-trip must hold.
919    #[test]
920    fn arri_logc4_at_key_points() {
921        use crate::color::{arri_logc4_constants, arri_logc4_eotf, arri_logc4_oetf};
922        // Spot-check: constants from the spec (Cooper & Brendel, 2022).
923        // Reference values computed independently with Python and match
924        // colour-science/colour to 12+ decimal places.
925        let (a, b, c, s, t) = arri_logc4_constants();
926        assert!((a - 2231.8263091).abs() < 1e-3, "a = {} (want 2231.8263)", a);
927        assert!((b - 0.90713587).abs() < 1e-6, "b = {} (want 0.9071)", b);
928        assert!((c - 0.09286413).abs() < 1e-6, "c = {} (want 0.0929)", c);
929        assert!((s - 0.1135972).abs() < 1e-5, "s = {} (want 0.1135972)", s);
930        assert!((t - (-0.0180570)).abs() < 1e-5, "t = {} (want -0.0180570)", t);
931
932        // Spec: 18% grey → ≈ 0.2783958.
933        let v_18 = arri_logc4_oetf(0.18);
934        assert!((v_18 - 0.2783958).abs() < 1e-5, "LogC4 OETF(0.18) = {} (want 0.2783958)", v_18);
935
936        // Spec: scene-linear 1.0 → ≈ 0.4275194 (unbounded formula;
937        // the hardware form clamps to 1.0 for highlights).
938        let v_one = arri_logc4_oetf(1.0);
939        let expected_one = (((a * 1.0 + 64.0).log2() - 6.0) / 14.0) * b + c;
940        assert!((v_one - expected_one).abs() < 1e-5, "LogC4 OETF(1.0) = {} (want {})", v_one, expected_one);
941        assert!((v_one - 0.4275194).abs() < 1e-5, "LogC4 OETF(1.0) = {} (want 0.4275194)", v_one);
942
943        // Linear branch (x < t ≈ -0.018): pure slope.
944        let v_below = arri_logc4_oetf(t - 0.001);
945        let expected_below = (t - 0.001 - t) / s; // = -0.001 / s
946        assert!((v_below - expected_below).abs() < 1e-5, "LogC4 linear branch");
947
948        // Round-trip: decode the encoded 18% grey back to scene-linear.
949        let rt = arri_logc4_eotf(v_18);
950        assert!((rt - 0.18).abs() < 1e-4, "LogC4 round-trip: encode→decode(0.18) = {} (want 0.18)", rt);
951
952        // Round-trip for a couple more stops.
953        for x in [0.001_f32, 0.01, 0.1, 0.5, 2.0, 10.0] {
954            let enc = arri_logc4_oetf(x);
955            let dec = arri_logc4_eotf(enc);
956            assert!((dec - x).abs() < 1e-4, "LogC4 round-trip at x={}: encode→decode = {} (want {})", x, dec, x);
957        }
958
959        // Sanity-check the full TransferFunction::ARRIlog4 path agrees with
960        // the standalone helper (so the production code is correct).
961        let v_18_full = TransferFunction::ARRIlog4.process_apply(0.18);
962        assert!((v_18_full - v_18).abs() < 1e-5, "TransferFunction::ARRIlog4 disagrees with arri_logc4_oetf: {} vs {}", v_18_full, v_18);
963    }
964
965    /// S-Log3 must follow Sony's canonical form per the Sony specification
966    /// (2014), colour-science, and ACES CTL reference implementation.
967    /// Formula:
968    ///   x >= 0.01125: V = (420 + 261.5 * log10((x + 0.01) / 0.19)) / 1023
969    ///   x <  0.01125: V = (x * (knee_val - 95) / 0.01125 + 95) / 1023
970    ///                 where knee_val = 420 + 261.5 * log10((0.01125+0.01)/0.19)
971    ///
972    /// 18% grey (x=0.18) maps to code 420, normalized 420/1023 ≈ 0.4106.
973    /// Black (x=0) maps to code 95, normalized 95/1023 ≈ 0.0929.
974    /// These match the known Sony S-Log3 encoding and DaVinci Resolve.
975    #[test]
976    fn slog3_canonical_at_key_points() {
977        let v_low = TransferFunction::SLog3.process_apply(0.009);
978        let v_at = TransferFunction::SLog3.process_apply(0.01125);
979        let v_high = TransferFunction::SLog3.process_apply(0.013);
980        assert!(v_low.is_finite() && v_at.is_finite() && v_high.is_finite());
981        assert!(v_low < v_high, "S-Log3 must be monotonic across the knee: low={} high={}", v_low, v_high);
982        // Spot-check x=0.18 (18% grey). Canonical S-Log3 gives:
983        //   V(0.18) = 420/1023 ≈ 0.4106 (code 420)
984        let v_18 = TransferFunction::SLog3.process_apply(0.18);
985        assert!((v_18 - 0.4106).abs() < 0.01, "S-Log3 at 0.18 = {} (want ~0.4106)", v_18);
986        // Spot-check x=1.0 (peak white, V ≈ 0.596, code ~610).
987        let v_1 = TransferFunction::SLog3.process_apply(1.0);
988        assert!((v_1 - 0.596).abs() < 0.02, "S-Log3 at 1.0 = {} (want ~0.596)", v_1);
989        // Black (x=0) should be code 95.
990        let v_0 = TransferFunction::SLog3.process_apply(0.0);
991        assert!((v_0 - 0.0929).abs() < 0.001, "S-Log3 at 0 = {} (want ~0.0929)", v_0);
992    }
993}
994
995// Tiny helper so the unit tests can invoke TransferFunction::process on
996// single pixels without spinning up rayon. Mirrors the per-pixel math
997// in the existing match arms exactly; if a new variant is added this
998// must be updated.
999impl TransferFunction {
1000    #[cfg(test)]
1001    fn process_apply(&self, x: f32) -> f32 {
1002        match self {
1003            TransferFunction::Linear => x,
1004            TransferFunction::Rec709 => rec709_oetf(x).min(1.0).max(0.0),
1005            TransferFunction::SLog3 => if x >= 0.01125_f32 { (420.0_f32 + 261.5_f32 * ((x + 0.01_f32) / 0.19_f32).log10()) / 1023.0_f32 } else { (x * (171.2102946929_f32 - 95.0_f32) / 0.01125_f32 + 95.0_f32) / 1023.0_f32 },
1006            TransferFunction::VLog => if x < 0.01 { 5.6_f32 * x + 0.125_f32 } else { 0.241514_f32 * (x + 0.00873_f32).log10() + 0.598206_f32 },
1007            TransferFunction::ARRIlog3 => if x > 0.010591_f32 { 0.247190_f32 * (5.555556_f32 * x + 0.052272_f32).log10() + 0.385537_f32 } else { 5.367655_f32 * x + 0.092809_f32 },
1008            TransferFunction::ARRIlog4 => {
1009                let (a, b, c, s, t) = crate::color::arri_logc4_constants();
1010                if x >= t { ((a * x + 64.0_f32).log2() - 6.0_f32) / 14.0_f32 * b + c } else { (x - t) / s }
1011            },
1012            TransferFunction::CLog3 => {
1013                let neg_graft_lin = (0.097465473_f32 - 0.12512219_f32) / 1.9754798_f32;
1014                let pos_graft_lin = (0.15277891_f32 - 0.12512219_f32) / 1.9754798_f32;
1015                if x < neg_graft_lin { -0.36726845_f32 * ((-x * 14.98325_f32 + 1.0_f32).max(1e-10_f32)).log10() + 0.12783901_f32 }
1016                else if x <= pos_graft_lin { 1.9754798_f32 * x + 0.12512219_f32 }
1017                else { 0.36726845_f32 * (x * 14.98325_f32 + 1.0_f32).log10() + 0.12240537_f32 }
1018            }
1019            TransferFunction::FLog2 => if x >= 0.000889_f32 { 0.245281_f32 * (5.555556_f32 * x + 0.064829_f32).log10() + 0.384316_f32 } else { 8.799461_f32 * x + 0.092864_f32 },
1020            TransferFunction::AppleLog | TransferFunction::AppleLog2 => {
1021                const R0: f32 = -0.05641088; const RT: f32 = 0.01; const C: f32 = 47.28711236;
1022                const BETA: f32 = 0.00964052; const GAMMA: f32 = 0.08550479; const DELTA: f32 = 0.69336945;
1023                if x < R0 { 0.0 } else if x < RT { C * (x - R0) * (x - R0) } else { GAMMA * (x + BETA).log2() + DELTA }
1024            }
1025            TransferFunction::ACESCCT => if x > 0.0078125_f32 { (x.log2() + 9.72_f32) / 17.52_f32 } else { 10.5402377416545_f32 * x + 0.0729055341958355_f32 },
1026            TransferFunction::PQ => { let x_m1 = x.powf(0.1593017578125_f32); ((0.8359375_f32 + 18.8515625_f32 * x_m1) / (1.0_f32 + 18.6875_f32 * x_m1)).powf(78.84375_f32) }
1027            TransferFunction::HLG => if x < (1.0_f32 / 12.0_f32) { (3.0_f32 * x).sqrt() } else { 0.17883277_f32 * (12.0_f32 * x - 0.28466892_f32).ln() + 0.55991073_f32 },
1028            TransferFunction::DaVinciIntermediate => if x <= 0.00262409_f32 { x * 10.44426855_f32 } else { 0.07329248_f32 * ((x + 0.0075_f32).log2() + 7.0_f32) },
1029            TransferFunction::Gamma24 => x.max(0.0).powf(1.0 / 2.4),
1030        }
1031    }
1032}