qualia_core_db/render/
gpu_colour_kernel.rs

1//! P7.4 — GPU colour-projection / gamut batch kernel + CPU oracle.
2//!
3//! The GPU kernel processes batches of EMF payloads `[α, μ, σ]` through
4//! the full colour pipeline (SPD → XYZ → linear sRGB → gamut map →
5//! display sRGB) in parallel. This module provides:
6//!
7//! 1. **CPU oracle**: the reference implementation for correctness checks.
8//! 2. **GPU kernel spec**: the WGSL shader source and buffer layout
9//!    specification for the GPU implementation.
10//! 3. **Differential check**: compares CPU vs GPU output for a given batch.
11//!
12//! ## Buffer layout (GPU)
13//!
14//! ```text
15//! Bind group 0:
16//!   - storage buffer: EMF input  [f32×3 × N]  (α, μ, σ per element)
17//!   - storage buffer: RGB output [u8×3 × N]   (r, g, b per element)
18//!   - uniform buffer: params { sample_count: u32 }
19//! ```
20//!
21//! ## Determinism
22//!
23//! The CPU oracle is fully deterministic. The GPU kernel must produce
24//! bit-identical output to the CPU oracle for the same input (within
25//! f32 rounding tolerance).
26
27use crate::render::spectral_kernel::{emf_to_spd, spd_to_xyz, xyz_to_linear_srgb};
28
29// ───────────────────────────────────────────────────────────────────────────
30//  CPU oracle — batch EMF → gamut-mapped display RGB
31// ───────────────────────────────────────────────────────────────────────────
32
33/// CPU oracle: process a batch of EMF payloads through the full colour
34/// pipeline and produce gamut-mapped 8-bit display RGB.
35///
36/// `emf_in` is `[α, μ, σ]` triples (3*N floats).
37/// `rgb_out` is N*3 u8.
38/// Returns the number of elements processed.
39pub fn cpu_batch_emf_to_display_gamut_mapped(emf_in: &[f32], rgb_out: &mut [u8]) -> usize {
40    let n = emf_in.len() / 3;
41    let n = n.min(rgb_out.len() / 3);
42
43    for i in 0..n {
44        let alpha = emf_in[i * 3];
45        let mu = emf_in[i * 3 + 1];
46        let sigma = emf_in[i * 3 + 2];
47
48        // EMF → SPD → XYZ → linear sRGB.
49        let spd = emf_to_spd(alpha, mu, sigma);
50        let xyz = spd_to_xyz(&spd);
51        let rgb = xyz_to_linear_srgb(&xyz);
52
53        // Gamut map: clamp to [0,1].
54        let r = rgb.r.clamp(0.0, 1.0);
55        let g = rgb.g.clamp(0.0, 1.0);
56        let b = rgb.b.clamp(0.0, 1.0);
57
58        // sRGB gamma encode.
59        let enc = |c: f32| -> u8 {
60            let encoded = if c <= 0.0031308 {
61                12.92 * c
62            } else {
63                1.055 * c.powf(1.0 / 2.4) - 0.055
64            };
65            (encoded * 255.0).round().clamp(0.0, 255.0) as u8
66        };
67
68        rgb_out[i * 3] = enc(r);
69        rgb_out[i * 3 + 1] = enc(g);
70        rgb_out[i * 3 + 2] = enc(b);
71    }
72
73    n
74}
75
76// ───────────────────────────────────────────────────────────────────────────
77//  CPU oracle — batch EMF → XYZ (for differential testing)
78// ───────────────────────────────────────────────────────────────────────────
79
80/// CPU oracle: process a batch of EMF payloads and produce XYZ values.
81///
82/// `emf_in` is `[α, μ, σ]` triples (3*N floats).
83/// `xyz_out` is N*3 f32.
84pub fn cpu_batch_emf_to_xyz(emf_in: &[f32], xyz_out: &mut [f32]) -> usize {
85    let n = emf_in.len() / 3;
86    let n = n.min(xyz_out.len() / 3);
87
88    for i in 0..n {
89        let spd = emf_to_spd(emf_in[i * 3], emf_in[i * 3 + 1], emf_in[i * 3 + 2]);
90        let xyz = spd_to_xyz(&spd);
91        xyz_out[i * 3] = xyz.x;
92        xyz_out[i * 3 + 1] = xyz.y;
93        xyz_out[i * 3 + 2] = xyz.z;
94    }
95
96    n
97}
98
99// ───────────────────────────────────────────────────────────────────────────
100//  GPU kernel specification (WGSL source + buffer layout)
101// ───────────────────────────────────────────────────────────────────────────
102
103/// WGSL shader source for the GPU colour-projection batch kernel.
104///
105/// This is the GPU twin of `cpu_batch_emf_to_display_gamut_mapped`.
106/// It processes one EMF element per workgroup invocation.
107pub const GPU_COLOUR_KERNEL_WGSL: &str = r#"
108// P7.4 — GPU colour-projection / gamut batch kernel.
109// Processes one EMF [α, μ, σ] element per invocation.
110
111const SPD_SAMPLES: u32 = 41u;
112const LAMBDA_MIN: f32 = 380.0;
113const LAMBDA_STEP: f32 = 10.0;
114
115// CIE 1931 2-degree observer CMFs (41 samples, 10nm steps).
116const CMF_X = array<f32, 41>(
117    0.001368, 0.004243, 0.014310, 0.043510, 0.134380, 0.283900, 0.348280,
118    0.336200, 0.290800, 0.195360, 0.095640, 0.032010, 0.004900, 0.009300,
119    0.063270, 0.165500, 0.290400, 0.433450, 0.594500, 0.762100, 0.916300,
120    1.026300, 1.062200, 1.002600, 0.854450, 0.642400, 0.447900, 0.283500,
121    0.164900, 0.087400, 0.046770, 0.022700, 0.011359, 0.005790, 0.002899,
122    0.001440, 0.000690, 0.000332, 0.000166, 0.000083, 0.000042,
123);
124
125const CMF_Y = array<f32, 41>(
126    0.000039, 0.000120, 0.000396, 0.001210, 0.004000, 0.011600, 0.023000,
127    0.038000, 0.060000, 0.090980, 0.139020, 0.208020, 0.323000, 0.503000,
128    0.710000, 0.862000, 0.954000, 0.994950, 0.995000, 0.952000, 0.870000,
129    0.757000, 0.631000, 0.503000, 0.381000, 0.265000, 0.175000, 0.107000,
130    0.061000, 0.032000, 0.017000, 0.008210, 0.004102, 0.002091, 0.001047,
131    0.000520, 0.000249, 0.000120, 0.000060, 0.000030, 0.000015,
132);
133
134const CMF_Z = array<f32, 41>(
135    0.006450, 0.020050, 0.067850, 0.207400, 0.645600, 1.282500, 1.453000,
136    1.562100, 1.562700, 1.385600, 1.114600, 0.777500, 0.445600, 0.198700,
137    0.068100, 0.019800, 0.004100, 0.000500, 0.000200, 0.000010, 0.000000,
138    0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
139    0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
140    0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
141);
142
143// Precomputed Y normalisation = 1 / Σȳ.
144const Y_NORM: f32 = 0.01068128;
145
146// XYZ → linear sRGB matrix.
147const M_R: vec3<f32> = vec3(3.2404542, -1.5371385, -0.4985314);
148const M_G: vec3<f32> = vec3(-0.9692660, 1.8760108, 0.0415560);
149const M_B: vec3<f32> = vec3(0.0556434, -0.2040259, 1.0572252);
150
151struct EmfInput {
152    data: array<f32>,
153};
154
155struct RgbOutput {
156    data: array<u32>,
157};
158
159@group(0) @binding(0) var<storage, read> emf_in: EmfInput;
160@group(0) @binding(1) var<storage, read_write> rgb_out: RgbOutput;
161@group(0) @binding(2) var<uniform> params: Params;
162
163struct Params {
164    count: u32,
165};
166
167@compute @workgroup_size(64)
168fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
169    let idx = gid.x;
170    if (idx >= params.count) {
171        return;
172    }
173
174    let alpha = emf_in.data[idx * 3u];
175    let mu = emf_in.data[idx * 3u + 1u];
176    let sigma = emf_in.data[idx * 3u + 2u];
177
178    // EMF → SPD: Gaussian peak.
179    let lambda = 400.0 + clamp(sigma, 0.0, 1.0) * 300.0;
180    let width = 10.0 + clamp(mu, 0.0, 1.0) * 140.0;
181    let amplitude = max(alpha, 0.0);
182
183    // SPD → XYZ.
184    var x = 0.0;
185    var y = 0.0;
186    var z = 0.0;
187    for (var i = 0u; i < SPD_SAMPLES; i = i + 1u) {
188        let l = LAMBDA_MIN + f32(i) * LAMBDA_STEP;
189        let d = (l - lambda) / width;
190        let s = amplitude * exp(-0.5 * d * d);
191        x = x + s * CMF_X[i];
192        y = y + s * CMF_Y[i];
193        z = z + s * CMF_Z[i];
194    }
195    x = x * Y_NORM;
196    y = y * Y_NORM;
197    z = z * Y_NORM;
198
199    // XYZ → linear sRGB.
200    let xyz = vec3<f32>(x, y, z);
201    var r = dot(M_R, xyz);
202    var g = dot(M_G, xyz);
203    var b = dot(M_B, xyz);
204
205    // Gamut map: clamp to [0, 1].
206    r = clamp(r, 0.0, 1.0);
207    g = clamp(g, 0.0, 1.0);
208    b = clamp(b, 0.0, 1.0);
209
210    // sRGB gamma encode.
211    let encode = fn(c: f32) -> f32 {
212        if (c <= 0.0031308) {
213            return 12.92 * c;
214        }
215        return 1.055 * pow(c, 1.0 / 2.4) - 0.055;
216    };
217
218    let er = encode(r);
219    let eg = encode(g);
220    let eb = encode(b);
221
222    // Pack into u32 (RGB24).
223    let packed = u32(round(er * 255.0)) |
224                 (u32(round(eg * 255.0)) << 8u) |
225                 (u32(round(eb * 255.0)) << 16u);
226
227    rgb_out.data[idx] = packed;
228}
229"#;
230
231/// Buffer layout specification for the GPU kernel.
232#[derive(Debug, Clone, Copy, PartialEq, Eq)]
233pub struct GpuBufferLayout {
234    pub emf_stride: usize,
235    pub rgb_stride: usize,
236    pub params_size: usize,
237}
238
239impl Default for GpuBufferLayout {
240    fn default() -> Self {
241        Self {
242            emf_stride: 12, // 3 × f32
243            rgb_stride: 4,  // 1 × u32 (packed RGB24)
244            params_size: 4, // 1 × u32 (count)
245        }
246    }
247}
248
249// ───────────────────────────────────────────────────────────────────────────
250//  Differential check (CPU vs GPU)
251// ───────────────────────────────────────────────────────────────────────────
252
253/// Tolerance for CPU vs GPU differential comparison.
254/// f32 rounding in the GPU kernel may produce ±1 LSB differences
255/// in the 8-bit output.
256pub const DIFF_TOLERANCE_U8: u8 = 2;
257
258/// Compare CPU oracle output with GPU output.
259///
260/// `cpu_out` and `gpu_out` are N*3 u8 arrays. Returns the number of
261/// elements that differ by more than `DIFF_TOLERANCE_U8` in any channel.
262pub fn diff_cpu_gpu(cpu_out: &[u8], gpu_out: &[u8]) -> usize {
263    let n = cpu_out.len().min(gpu_out.len()) / 3;
264    let mut mismatches = 0;
265
266    for i in 0..n {
267        let dr = (cpu_out[i * 3] as i16 - gpu_out[i * 3] as i16).unsigned_abs() as u8;
268        let dg = (cpu_out[i * 3 + 1] as i16 - gpu_out[i * 3 + 1] as i16).unsigned_abs() as u8;
269        let db = (cpu_out[i * 3 + 2] as i16 - gpu_out[i * 3 + 2] as i16).unsigned_abs() as u8;
270        if dr > DIFF_TOLERANCE_U8 || dg > DIFF_TOLERANCE_U8 || db > DIFF_TOLERANCE_U8 {
271            mismatches += 1;
272        }
273    }
274
275    mismatches
276}
277
278// ───────────────────────────────────────────────────────────────────────────
279//  Tests
280// ───────────────────────────────────────────────────────────────────────────
281
282#[cfg(test)]
283mod tests {
284    use super::*;
285
286    #[test]
287    fn cpu_batch_deterministic() {
288        let emf = [1.0, 0.3, 0.5, 0.8, 0.2, 0.7, 1.0, 0.0, 0.0];
289        let mut out1 = [0u8; 9];
290        let mut out2 = [0u8; 9];
291        cpu_batch_emf_to_display_gamut_mapped(&emf, &mut out1);
292        cpu_batch_emf_to_display_gamut_mapped(&emf, &mut out2);
293        assert_eq!(out1, out2, "CPU batch must be deterministic");
294    }
295
296    #[test]
297    fn cpu_batch_xyz_deterministic() {
298        let emf = [1.0, 0.3, 0.5, 0.8, 0.2, 0.7];
299        let mut out1 = [0.0f32; 6];
300        let mut out2 = [0.0f32; 6];
301        cpu_batch_emf_to_xyz(&emf, &mut out1);
302        cpu_batch_emf_to_xyz(&emf, &mut out2);
303        assert_eq!(out1, out2, "CPU batch XYZ must be deterministic");
304    }
305
306    #[test]
307    fn cpu_batch_valid_rgb() {
308        let emf = [1.0, 0.0, 0.0, 1.0, 0.5, 0.5, 1.0, 0.0, 1.0];
309        let mut out = [0u8; 9];
310        let n = cpu_batch_emf_to_display_gamut_mapped(&emf, &mut out);
311        assert_eq!(n, 3);
312        // RGB is u8, so it's always in [0, 255] by type limits.
313    }
314
315    #[test]
316    fn cpu_batch_partial_output() {
317        let emf = [1.0, 0.3, 0.5, 0.8, 0.2, 0.7, 1.0, 0.0, 0.0];
318        let mut out = [0u8; 6]; // only room for 2
319        let n = cpu_batch_emf_to_display_gamut_mapped(&emf, &mut out);
320        assert_eq!(n, 2, "should process only 2 elements");
321    }
322
323    #[test]
324    fn gpu_kernel_source_contains_cmf_data() {
325        assert!(
326            GPU_COLOUR_KERNEL_WGSL.contains("CMF_X"),
327            "WGSL must contain CMF X data"
328        );
329        assert!(
330            GPU_COLOUR_KERNEL_WGSL.contains("CMF_Y"),
331            "WGSL must contain CMF Y data"
332        );
333        assert!(
334            GPU_COLOUR_KERNEL_WGSL.contains("CMF_Z"),
335            "WGSL must contain CMF Z data"
336        );
337    }
338
339    #[test]
340    fn gpu_kernel_source_contains_pipeline() {
341        assert!(
342            GPU_COLOUR_KERNEL_WGSL.contains("SPD"),
343            "WGSL must mention SPD"
344        );
345        assert!(
346            GPU_COLOUR_KERNEL_WGSL.contains("XYZ"),
347            "WGSL must mention XYZ"
348        );
349        assert!(
350            GPU_COLOUR_KERNEL_WGSL.contains("sRGB"),
351            "WGSL must mention sRGB"
352        );
353        assert!(
354            GPU_COLOUR_KERNEL_WGSL.contains("clamp"),
355            "WGSL must gamut-map with clamp"
356        );
357    }
358
359    #[test]
360    fn gpu_kernel_source_has_workgroup_size() {
361        assert!(
362            GPU_COLOUR_KERNEL_WGSL.contains("workgroup_size(64)"),
363            "WGSL must specify workgroup size"
364        );
365    }
366
367    #[test]
368    fn buffer_layout_defaults() {
369        let layout = GpuBufferLayout::default();
370        assert_eq!(layout.emf_stride, 12);
371        assert_eq!(layout.rgb_stride, 4);
372        assert_eq!(layout.params_size, 4);
373    }
374
375    #[test]
376    fn diff_cpu_gpu_self_zero() {
377        let cpu = [100u8, 150, 200, 50, 60, 70];
378        let mismatches = diff_cpu_gpu(&cpu, &cpu);
379        assert_eq!(mismatches, 0, "self-comparison should have zero mismatches");
380    }
381
382    #[test]
383    fn diff_cpu_gpu_within_tolerance() {
384        let cpu = [100u8, 150, 200, 50, 60, 70];
385        let gpu = [101u8, 149, 202, 51, 58, 72]; // all within ±2
386        let mismatches = diff_cpu_gpu(&cpu, &gpu);
387        assert_eq!(
388            mismatches, 0,
389            "within-tolerance differences should not count"
390        );
391    }
392
393    #[test]
394    fn diff_cpu_gpu_outside_tolerance() {
395        let cpu = [100u8, 150, 200, 50, 60, 70];
396        let gpu = [110u8, 150, 200, 50, 60, 70]; // dr=10 > 2
397        let mismatches = diff_cpu_gpu(&cpu, &gpu);
398        assert_eq!(mismatches, 1, "out-of-tolerance should count as mismatch");
399    }
400
401    #[test]
402    fn cpu_batch_blue_dominates_at_low_sigma() {
403        let emf = [1.0, 0.0, 0.0]; // σ=0 → blue
404        let mut out = [0u8; 3];
405        cpu_batch_emf_to_display_gamut_mapped(&emf, &mut out);
406        assert!(out[2] >= out[0], "B should dominate R at σ=0");
407        assert!(out[2] >= out[1], "B should dominate G at σ=0");
408    }
409
410    #[test]
411    fn cpu_batch_red_dominates_at_high_sigma() {
412        let emf = [1.0, 0.0, 1.0]; // σ=1 → red
413        let mut out = [0u8; 3];
414        cpu_batch_emf_to_display_gamut_mapped(&emf, &mut out);
415        assert!(out[0] >= out[1], "R should dominate G at σ=1");
416        assert!(out[0] >= out[2], "R should dominate B at σ=1");
417    }
418}
qualia_core_db/render/gpu_colour_kernel.rs

qualia_core_db/render/
gpu_colour_kernel.rs
