qualia_core_db/render/
spectral_blend.rs

1//! P7.3 — σ spectral blend as interpolation on the spectral manifold.
2//!
3//! Blending colours in spectral space (interpolating SPDs) is physically
4//! correct: it preserves the spectral power distribution's structure,
5//! unlike RGB lerp which can produce intermediate colours that don't
6//! correspond to any physical spectrum.
7//!
8//! ## Algorithm
9//!
10//! Spectral blend: `S_blend(t) = (1-t)·S_a + t·S_b` in SPD space, then
11//! project to XYZ/RGB.
12//!
13//! RGB lerp: `RGB_blend(t) = (1-t)·RGB_a + t·RGB_b` in sRGB space.
14//!
15//! The divergence between the two is the ΔE of their XYZ projections.
16//!
17//! ## Determinism
18//!
19//! All operations are deterministic: SPD lerp is a pure function, and
20//! the CMF projection is a compile-time constant.
21
22use super::spectral_kernel::{
23    delta_e_76, emf_to_spd, linear_to_srgb_channel, spd_to_xyz, xyz_to_linear_srgb, Spd, Xyz,
24};
25
26// ───────────────────────────────────────────────────────────────────────────
27//  Spectral blend
28// ───────────────────────────────────────────────────────────────────────────
29
30/// Blend two SPDs in spectral space: `S(t) = (1-t)·a + t·b`.
31#[inline]
32pub fn spectral_blend_spd(a: &Spd, b: &Spd, t: f32) -> Spd {
33    a.lerp(b, t)
34}
35
36/// Blend two EMF payloads in spectral space and return the resulting XYZ.
37#[inline]
38pub fn spectral_blend_emf(
39    alpha_a: f32,
40    mu_a: f32,
41    sigma_a: f32,
42    alpha_b: f32,
43    mu_b: f32,
44    sigma_b: f32,
45    t: f32,
46) -> Xyz {
47    let spd_a = emf_to_spd(alpha_a, mu_a, sigma_a);
48    let spd_b = emf_to_spd(alpha_b, mu_b, sigma_b);
49    let blended = spectral_blend_spd(&spd_a, &spd_b, t);
50    spd_to_xyz(&blended)
51}
52
53/// RGB lerp in linear sRGB space (for comparison with spectral blend).
54#[inline]
55pub fn rgb_lerp(rgb_a: [f32; 3], rgb_b: [f32; 3], t: f32) -> [f32; 3] {
56    [
57        rgb_a[0] * (1.0 - t) + rgb_b[0] * t,
58        rgb_a[1] * (1.0 - t) + rgb_b[1] * t,
59        rgb_a[2] * (1.0 - t) + rgb_b[2] * t,
60    ]
61}
62
63/// Compute the ΔE divergence between spectral blend and gamma-encoded
64/// sRGB lerp at parameter `t`.
65///
66/// The CMF projection is linear, so SPD lerp = XYZ lerp. The divergence
67/// with "RGB lerp" appears because display RGB lerp happens in gamma-encoded
68/// sRGB space, which is non-linear. We decode the gamma-encoded endpoints,
69/// lerp in linear space, then re-encode — no, that would be the same.
70///
71/// The real divergence is: spectral blend produces an SPD with two peaks
72/// (for narrow-band inputs), which when projected gives a different XYZ
73/// than the gamma-encoded sRGB lerp. We simulate the common "lerp in
74/// 8-bit sRGB" workflow: encode both endpoints to 8-bit sRGB, lerp the
75/// 8-bit values, decode back to linear, convert to XYZ.
76pub fn blend_divergence(
77    alpha_a: f32,
78    mu_a: f32,
79    sigma_a: f32,
80    alpha_b: f32,
81    mu_b: f32,
82    sigma_b: f32,
83    t: f32,
84) -> f32 {
85    let spd_a = emf_to_spd(alpha_a, mu_a, sigma_a);
86    let spd_b = emf_to_spd(alpha_b, mu_b, sigma_b);
87
88    // Spectral blend → XYZ.
89    let blended_spd = spectral_blend_spd(&spd_a, &spd_b, t);
90    let xyz_spectral = spd_to_xyz(&blended_spd);
91
92    // sRGB 8-bit lerp: encode both to 8-bit sRGB, lerp, decode back.
93    let xyz_a = spd_to_xyz(&spd_a);
94    let xyz_b = spd_to_xyz(&spd_b);
95    let rgb_a = xyz_to_linear_srgb(&xyz_a);
96    let rgb_b = xyz_to_linear_srgb(&xyz_b);
97
98    // Encode to 8-bit sRGB (gamma encode + quantise to 0-255).
99    let enc_a = [
100        (linear_to_srgb_channel(rgb_a.r) * 255.0).round() / 255.0,
101        (linear_to_srgb_channel(rgb_a.g) * 255.0).round() / 255.0,
102        (linear_to_srgb_channel(rgb_a.b) * 255.0).round() / 255.0,
103    ];
104    let enc_b = [
105        (linear_to_srgb_channel(rgb_b.r) * 255.0).round() / 255.0,
106        (linear_to_srgb_channel(rgb_b.g) * 255.0).round() / 255.0,
107        (linear_to_srgb_channel(rgb_b.b) * 255.0).round() / 255.0,
108    ];
109
110    // Lerp in gamma-encoded space.
111    let enc_blend = [
112        enc_a[0] * (1.0 - t) + enc_b[0] * t,
113        enc_a[1] * (1.0 - t) + enc_b[1] * t,
114        enc_a[2] * (1.0 - t) + enc_b[2] * t,
115    ];
116
117    // Decode back to linear (inverse sRGB gamma).
118    let decode_channel = |c: f32| -> f32 {
119        if c <= 0.04045 {
120            c / 12.92
121        } else {
122            ((c + 0.055) / 1.055).powf(2.4)
123        }
124    };
125    let rgb_blend = [
126        decode_channel(enc_blend[0]),
127        decode_channel(enc_blend[1]),
128        decode_channel(enc_blend[2]),
129    ];
130
131    // Convert blended linear RGB to XYZ.
132    let xyz_srgb_lerp = Xyz::new(
133        0.4124564 * rgb_blend[0] + 0.3575761 * rgb_blend[1] + 0.1804375 * rgb_blend[2],
134        0.2126729 * rgb_blend[0] + 0.7151522 * rgb_blend[1] + 0.0721750 * rgb_blend[2],
135        0.0193339 * rgb_blend[0] + 0.1191920 * rgb_blend[1] + 0.9503041 * rgb_blend[2],
136    );
137
138    delta_e_76(&xyz_spectral, &xyz_srgb_lerp)
139}
140
141// ───────────────────────────────────────────────────────────────────────────
142//  Tests
143// ───────────────────────────────────────────────────────────────────────────
144
145#[cfg(test)]
146mod tests {
147    use super::*;
148
149    #[test]
150    fn blend_at_t0_returns_first() {
151        let spd_a = emf_to_spd(1.0, 0.1, 0.2);
152        let spd_b = emf_to_spd(1.0, 0.1, 0.8);
153        let blended = spectral_blend_spd(&spd_a, &spd_b, 0.0);
154        assert_eq!(blended, spd_a, "t=0 should return first SPD");
155    }
156
157    #[test]
158    fn blend_at_t1_returns_second() {
159        let spd_a = emf_to_spd(1.0, 0.1, 0.2);
160        let spd_b = emf_to_spd(1.0, 0.1, 0.8);
161        let blended = spectral_blend_spd(&spd_a, &spd_b, 1.0);
162        assert_eq!(blended, spd_b, "t=1 should return second SPD");
163    }
164
165    #[test]
166    fn spectral_blend_differs_from_rgb_lerp() {
167        // Blend between a blue (σ=0) and red (σ=1) — the divergence
168        // should be significant because spectral blending preserves the
169        // two peaks while RGB lerp produces a mid-colour.
170        let de = blend_divergence(1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.5);
171        assert!(
172            de > 0.5,
173            "spectral blend should differ from RGB lerp: ΔE={}",
174            de
175        );
176    }
177
178    #[test]
179    fn blend_no_nans_in_sweep() {
180        for i in 0..=100 {
181            let t = i as f32 / 100.0;
182            let xyz = spectral_blend_emf(1.0, 0.2, 0.3, 0.8, 0.5, 0.7, t);
183            assert!(xyz.x.is_finite(), "X NaN at t={}", t);
184            assert!(xyz.y.is_finite(), "Y NaN at t={}", t);
185            assert!(xyz.z.is_finite(), "Z NaN at t={}", t);
186        }
187    }
188
189    #[test]
190    fn blend_determinism() {
191        let xyz1 = spectral_blend_emf(1.0, 0.2, 0.3, 0.8, 0.5, 0.7, 0.5);
192        let xyz2 = spectral_blend_emf(1.0, 0.2, 0.3, 0.8, 0.5, 0.7, 0.5);
193        assert_eq!(xyz1, xyz2, "blend must be deterministic");
194    }
195
196    #[test]
197    fn blend_monotone_continuity() {
198        // The blended XYZ should change continuously (no jumps).
199        let mut prev = spectral_blend_emf(1.0, 0.1, 0.2, 1.0, 0.1, 0.8, 0.0);
200        for i in 1..=100 {
201            let t = i as f32 / 100.0;
202            let curr = spectral_blend_emf(1.0, 0.1, 0.2, 1.0, 0.1, 0.8, t);
203            let dx = (curr.x - prev.x).abs();
204            let dy = (curr.y - prev.y).abs();
205            let dz = (curr.z - prev.z).abs();
206            assert!(dx < 0.1, "X discontinuity at t={}: {}", t, dx);
207            assert!(dy < 0.1, "Y discontinuity at t={}: {}", t, dy);
208            assert!(dz < 0.1, "Z discontinuity at t={}: {}", t, dz);
209            prev = curr;
210        }
211    }
212}
qualia_core_db/render/spectral_blend.rs

qualia_core_db/render/
spectral_blend.rs
