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

mcraw_tui/
color.rs
