indicatrix 0.1.0

A physically-based spectral gemstone renderer that turns GemCAD-style cutting schedules into rendered output.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121

//! CIE 1931 2° Standard Observer Color Matching Functions (CMFs), as an analytic fit.
//!
//! Single source of truth: [`optics::raytracer`](crate::optics::raytracer) delegates to
//! [`cie_1931_cmf`] rather than keeping its own copy of the lobe constants.

/// Piecewise-asymmetric Gaussian lobe: sigma differs on the rising (x < mu) vs falling
/// (x >= mu) side of the peak. Required by the Wyman/Sloan/Shirley fit; a symmetric
/// sigma distorts lobe shapes, worst on the z-bar 437 nm lobe's falling side.
#[inline]
fn g(x: f32, mu: f32, sigma_lo: f32, sigma_hi: f32) -> f32 {
    let sigma = if x < mu { sigma_lo } else { sigma_hi };
    let t = (x - mu) / sigma;
    (-0.5 * t * t).exp()
}

/// Wyman, Sloan, Shirley (2013) multi-lobe analytic fit to the CIE 1931 2° Standard
/// Observer Color Matching Functions.
///
/// Uses per-[`g`] asymmetric sigmas deliberately: a symmetric-sigma version skewed
/// chromaticities (equal-energy white ~= (0.33311, 0.33359) vs. (0.3413, 0.3644)).
#[must_use]
pub fn cie_1931_cmf(lambda_nm: f32) -> [f32; 3] {
    let l = lambda_nm;

    let x = 0.065f32.mul_add(
        -g(l, 501.1, 20.4, 26.2),
        1.056f32.mul_add(g(l, 599.8, 37.9, 31.0), 0.362 * g(l, 442.0, 16.0, 26.7)),
    );

    let y = 0.286f32.mul_add(g(l, 530.9, 16.3, 31.1), 0.821 * g(l, 568.8, 46.9, 40.5));

    let z = 0.681f32.mul_add(g(l, 459.0, 26.0, 13.8), 1.217 * g(l, 437.0, 11.8, 36.0));

    [x.max(0.0), y.max(0.0), z.max(0.0)]
}

/// Batched [`g`] over 8 wavelengths: same per-lane arithmetic, one
/// [`crate::simd::exp_f32x8`] call instead of 8 scalar `f32::exp` calls.
#[inline]
fn g_x8(lambdas: &[f32; 8], mu: f32, sigma_lo: f32, sigma_hi: f32) -> [f32; 8] {
    let mut args = [0f32; 8];
    for (k, &l) in lambdas.iter().enumerate() {
        let sigma = if l < mu { sigma_lo } else { sigma_hi };
        let t = (l - mu) / sigma;
        args[k] = -0.5 * t * t;
    }
    crate::simd::exp_f32x8(args)
}

/// Batched form of [`cie_1931_cmf`] over 8 wavelengths: one [`crate::simd::exp_f32x8`]
/// call per lobe (7 total) instead of 56 scalar `f32::exp` calls.
///
/// Not bit-identical to 8 scalar calls: every non-exponential step matches exactly, so
/// the only difference is `exp_f32x8` vs libm's `f32::exp`, a few ULP -- see
/// [`crate::simd::exp_f32x8`]'s docs for the determinism contract.
#[must_use]
pub fn cie_1931_cmf_x8(lambdas: &[f32; 8]) -> [[f32; 3]; 8] {
    let lobe_x0 = g_x8(lambdas, 501.1, 20.4, 26.2);
    let lobe_x1 = g_x8(lambdas, 599.8, 37.9, 31.0);
    let lobe_x2 = g_x8(lambdas, 442.0, 16.0, 26.7);
    let lobe_y0 = g_x8(lambdas, 530.9, 16.3, 31.1);
    let lobe_y1 = g_x8(lambdas, 568.8, 46.9, 40.5);
    let lobe_z0 = g_x8(lambdas, 459.0, 26.0, 13.8);
    let lobe_z1 = g_x8(lambdas, 437.0, 11.8, 36.0);

    let mut out = [[0f32; 3]; 8];
    for k in 0..8 {
        let x = 0.065f32.mul_add(
            -lobe_x0[k],
            1.056f32.mul_add(lobe_x1[k], 0.362 * lobe_x2[k]),
        );
        let y = 0.286f32.mul_add(lobe_y0[k], 0.821 * lobe_y1[k]);
        let z = 0.681f32.mul_add(lobe_z0[k], 1.217 * lobe_z1[k]);
        out[k] = [x.max(0.0), y.max(0.0), z.max(0.0)];
    }
    out
}

#[cfg(test)]
mod tests {
    use super::*;

    /// [`cie_1931_cmf_x8`] must agree with 8 separate [`cie_1931_cmf`] calls to a tight
    /// relative tolerance -- the documented ULP-level `exp_f32x8`-vs-libm gap, not more.
    #[test]
    fn cmf_x8_matches_scalar_closely() {
        let mut state = 99u64;
        let mut next = move || {
            state = state
                .wrapping_mul(6_364_136_223_846_793_005)
                .wrapping_add(1);
            (((state >> 33) as f32) / (u32::MAX as f32)).mul_add(700.0 - 380.0, 380.0)
        };
        for _round in 0..200 {
            let lambdas = [
                next(),
                next(),
                next(),
                next(),
                next(),
                next(),
                next(),
                next(),
            ];
            let batched = cie_1931_cmf_x8(&lambdas);
            for (k, &l) in lambdas.iter().enumerate() {
                let scalar = cie_1931_cmf(l);
                for c in 0..3 {
                    let a = f64::from(batched[k][c]);
                    let b = f64::from(scalar[c]);
                    let tol = 1e-5 * b.abs().max(1e-6);
                    assert!(
                        (a - b).abs() <= tol,
                        "lambda={l} channel={c}: batched={a} scalar={b} diff={}",
                        (a - b).abs()
                    );
                }
            }
        }
    }
}