Skip to main content

fast_mesh_fields/
color.rs

1//! Colormaps and colour rules: turning a scalar field value into an RGBA8
2//! colour.
3//!
4//! A [`ColorRule`] maps values in `range` onto a [`Colormap`]. Without a
5//! `center` the mapping is linear, `t = (v - lo) / (hi - lo)`. With a
6//! `center` `c` it is piecewise linear and pins `c` to the middle of the
7//! colormap: `t = 0.5 (v - lo) / (c - lo)` on `[lo, c]` and
8//! `t = 0.5 + 0.5 (v - c) / (hi - c)` on `[c, hi]`.
9//!
10//! The colormap is sampled at `t * 255` **with linear interpolation** between
11//! the two neighbouring table entries (per channel, rounded to nearest).
12//!
13//! Outside the range, [`ColorRule::below`] decides for `v < lo` and
14//! [`ColorRule::above`] for `v > hi`: [`OutOfRange::Clamp`] gives the end
15//! colour of the colormap, [`OutOfRange::Base`] gives [`ColorRule::base`].
16//! `-∞` follows `below` and `+∞` follows `above`; because the "far" value of
17//! a distance kernel is `+∞`, it follows `above` (`Base` shows the base
18//! colour, `Clamp` the end colour). NaN gives [`ColorRule::nan`].
19//!
20//! Built-in tables have 256 RGBA8 entries with alpha 255: Turbo is the
21//! published table, the others are generated on first use:
22//!
23//! - **Turbo**: Anton Mikhailov's published 256-entry look-up table of the
24//!   Turbo colormap (Google, 2019; Apache-2.0).
25//! - **Viridis**: Matt Zucker's published degree-6 polynomial fit of
26//!   matplotlib's viridis, evaluated at 256 points (within a few 8-bit
27//!   levels of matplotlib's table).
28//! - **Coolwarm**: Kenneth Moreland's diverging map ("Diverging Color Maps
29//!   for Scientific Visualization", 2009), interpolated in Msh space between
30//!   sRGB (59, 76, 192) and (180, 4, 38) through a white-point of M = 88.
31//! - **Grayscale**: linear, entry `i` is `(i, i, i, 255)`.
32
33use std::sync::OnceLock;
34
35use crate::{FieldsError, Result};
36
37/// Number of entries in each built-in colormap table.
38const TABLE_LEN: usize = 256;
39
40/// A colormap from `t` in `[0, 1]` to RGBA8.
41#[derive(Debug, Clone, PartialEq)]
42pub enum Colormap {
43    /// Google's Turbo rainbow (published table).
44    Turbo,
45    /// Matplotlib's perceptually uniform Viridis (polynomial fit).
46    Viridis,
47    /// Moreland's blue-white-red diverging map.
48    Coolwarm,
49    /// Linear black to white.
50    Grayscale,
51    /// Colour stops `(position, rgba)`: at least two, positions in `[0, 1]`
52    /// and strictly increasing. Colours are interpolated linearly per channel
53    /// and rounded; `t` before the first or after the last stop takes that
54    /// stop's colour.
55    Custom(Vec<(f32, [u8; 4])>),
56}
57
58/// What to show for a value outside the colour range.
59#[derive(Debug, Clone, Copy, PartialEq, Eq)]
60pub enum OutOfRange {
61    /// The colormap's end colour (`t = 0` below, `t = 1` above).
62    Clamp,
63    /// The rule's [`ColorRule::base`] colour.
64    Base,
65}
66
67/// How a scalar field is turned into colours.
68#[derive(Debug, Clone, PartialEq)]
69pub struct ColorRule {
70    /// The colormap.
71    pub map: Colormap,
72    /// `(lo, hi)`: values mapped onto the colormap; `lo < hi`, both finite.
73    pub range: (f32, f32),
74    /// Colour rule for values below `lo` (and `-∞`).
75    pub below: OutOfRange,
76    /// Colour rule for values above `hi` (and `+∞`).
77    pub above: OutOfRange,
78    /// A value pinned to the middle of the colormap, within `range`.
79    pub center: Option<f32>,
80    /// The colour of [`OutOfRange::Base`].
81    pub base: [u8; 4],
82    /// The colour of NaN values.
83    pub nan: [u8; 4],
84}
85
86impl ColorRule {
87    /// Checks the rule.
88    ///
89    /// # Errors
90    ///
91    /// [`FieldsError::InvalidInput`] if `range` has a non-finite end or
92    /// `lo >= hi`, if `center` is non-finite or outside `range`, or if a
93    /// custom colormap has fewer than two stops or stop positions that are
94    /// not finite, outside `[0, 1]` or not strictly increasing.
95    pub fn validate(&self) -> Result<()> {
96        let (lo, hi) = self.range;
97        if !lo.is_finite() || !hi.is_finite() {
98            return Err(FieldsError::InvalidInput("color range must be finite"));
99        }
100        if lo >= hi {
101            return Err(FieldsError::InvalidInput("color range needs lo < hi"));
102        }
103        if let Some(c) = self.center
104            && (!c.is_finite() || c < lo || c > hi)
105        {
106            return Err(FieldsError::InvalidInput(
107                "color center must lie within the range",
108            ));
109        }
110        if let Colormap::Custom(stops) = &self.map {
111            if stops.len() < 2 {
112                return Err(FieldsError::InvalidInput(
113                    "custom colormap needs at least two stops",
114                ));
115            }
116            let mut prev = -1.0_f32;
117            for &(p, _) in stops {
118                if !(0.0..=1.0).contains(&p) || p <= prev {
119                    return Err(FieldsError::InvalidInput(
120                        "custom stops must be in [0, 1] and strictly increasing",
121                    ));
122                }
123                prev = p;
124            }
125        }
126        Ok(())
127    }
128}
129
130/// The colour of value `v` under `rule`.
131///
132/// `rule` should have passed [`ColorRule::validate`]; an invalid rule never
133/// panics but gives unspecified colours.
134pub fn color_of(rule: &ColorRule, v: f32) -> [u8; 4] {
135    if v.is_nan() {
136        return rule.nan;
137    }
138    let (lo, hi) = rule.range;
139    if v < lo {
140        return match rule.below {
141            OutOfRange::Clamp => sample(&rule.map, 0.0),
142            OutOfRange::Base => rule.base,
143        };
144    }
145    if v > hi {
146        return match rule.above {
147            OutOfRange::Clamp => sample(&rule.map, 1.0),
148            OutOfRange::Base => rule.base,
149        };
150    }
151    let t = match rule.center {
152        None => (v - lo) / (hi - lo),
153        Some(c) => {
154            if v == c {
155                0.5
156            } else if v < c {
157                0.5 * (v - lo) / (c - lo)
158            } else {
159                0.5 + 0.5 * (v - c) / (hi - c)
160            }
161        }
162    };
163    sample(&rule.map, if t.is_nan() { 0.0 } else { t })
164}
165
166/// Samples a colormap at `t` (clamped to `[0, 1]`).
167fn sample(map: &Colormap, t: f32) -> [u8; 4] {
168    let t = t.clamp(0.0, 1.0);
169    match map {
170        Colormap::Custom(stops) => sample_stops(stops, t),
171        builtin => {
172            let table = table_of(builtin);
173            let x = t * (TABLE_LEN - 1) as f32;
174            let i = (x.floor() as usize).min(TABLE_LEN - 2);
175            lerp_rgba(table[i], table[i + 1], x - i as f32)
176        }
177    }
178}
179
180fn sample_stops(stops: &[(f32, [u8; 4])], t: f32) -> [u8; 4] {
181    let (Some(first), Some(last)) = (stops.first(), stops.last()) else {
182        return [0, 0, 0, 0];
183    };
184    if t <= first.0 {
185        return first.1;
186    }
187    if t >= last.0 {
188        return last.1;
189    }
190    for w in stops.windows(2) {
191        let (p0, c0) = w[0];
192        let (p1, c1) = w[1];
193        if t <= p1 {
194            let span = p1 - p0;
195            let f = if span > 0.0 { (t - p0) / span } else { 1.0 };
196            return lerp_rgba(c0, c1, f);
197        }
198    }
199    last.1
200}
201
202fn lerp_rgba(a: [u8; 4], b: [u8; 4], f: f32) -> [u8; 4] {
203    let mut out = [0u8; 4];
204    for k in 0..4 {
205        let v = f32::from(a[k]) + (f32::from(b[k]) - f32::from(a[k])) * f;
206        out[k] = v.round().clamp(0.0, 255.0) as u8;
207    }
208    out
209}
210
211/// Google's published 256-entry Turbo look-up table (Anton Mikhailov,
212/// 2019; Apache-2.0), as 8-bit sRGB triplets with opaque alpha.
213const TURBO_TABLE: [[u8; 4]; TABLE_LEN] = [
214    [48, 18, 59, 255],
215    [50, 21, 67, 255],
216    [51, 24, 74, 255],
217    [52, 27, 81, 255],
218    [53, 30, 88, 255],
219    [54, 33, 95, 255],
220    [55, 36, 102, 255],
221    [56, 39, 109, 255],
222    [57, 42, 115, 255],
223    [58, 45, 121, 255],
224    [59, 47, 128, 255],
225    [60, 50, 134, 255],
226    [61, 53, 139, 255],
227    [62, 56, 145, 255],
228    [63, 59, 151, 255],
229    [63, 62, 156, 255],
230    [64, 64, 162, 255],
231    [65, 67, 167, 255],
232    [65, 70, 172, 255],
233    [66, 73, 177, 255],
234    [66, 75, 181, 255],
235    [67, 78, 186, 255],
236    [68, 81, 191, 255],
237    [68, 84, 195, 255],
238    [68, 86, 199, 255],
239    [69, 89, 203, 255],
240    [69, 92, 207, 255],
241    [69, 94, 211, 255],
242    [70, 97, 214, 255],
243    [70, 100, 218, 255],
244    [70, 102, 221, 255],
245    [70, 105, 224, 255],
246    [70, 107, 227, 255],
247    [71, 110, 230, 255],
248    [71, 113, 233, 255],
249    [71, 115, 235, 255],
250    [71, 118, 238, 255],
251    [71, 120, 240, 255],
252    [71, 123, 242, 255],
253    [70, 125, 244, 255],
254    [70, 128, 246, 255],
255    [70, 130, 248, 255],
256    [70, 133, 250, 255],
257    [70, 135, 251, 255],
258    [69, 138, 252, 255],
259    [69, 140, 253, 255],
260    [68, 143, 254, 255],
261    [67, 145, 254, 255],
262    [66, 148, 255, 255],
263    [65, 150, 255, 255],
264    [64, 153, 255, 255],
265    [62, 155, 254, 255],
266    [61, 158, 254, 255],
267    [59, 160, 253, 255],
268    [58, 163, 252, 255],
269    [56, 165, 251, 255],
270    [55, 168, 250, 255],
271    [53, 171, 248, 255],
272    [51, 173, 247, 255],
273    [49, 175, 245, 255],
274    [47, 178, 244, 255],
275    [46, 180, 242, 255],
276    [44, 183, 240, 255],
277    [42, 185, 238, 255],
278    [40, 188, 235, 255],
279    [39, 190, 233, 255],
280    [37, 192, 231, 255],
281    [35, 195, 228, 255],
282    [34, 197, 226, 255],
283    [32, 199, 223, 255],
284    [31, 201, 221, 255],
285    [30, 203, 218, 255],
286    [28, 205, 216, 255],
287    [27, 208, 213, 255],
288    [26, 210, 210, 255],
289    [26, 212, 208, 255],
290    [25, 213, 205, 255],
291    [24, 215, 202, 255],
292    [24, 217, 200, 255],
293    [24, 219, 197, 255],
294    [24, 221, 194, 255],
295    [24, 222, 192, 255],
296    [24, 224, 189, 255],
297    [25, 226, 187, 255],
298    [25, 227, 185, 255],
299    [26, 228, 182, 255],
300    [28, 230, 180, 255],
301    [29, 231, 178, 255],
302    [31, 233, 175, 255],
303    [32, 234, 172, 255],
304    [34, 235, 170, 255],
305    [37, 236, 167, 255],
306    [39, 238, 164, 255],
307    [42, 239, 161, 255],
308    [44, 240, 158, 255],
309    [47, 241, 155, 255],
310    [50, 242, 152, 255],
311    [53, 243, 148, 255],
312    [56, 244, 145, 255],
313    [60, 245, 142, 255],
314    [63, 246, 138, 255],
315    [67, 247, 135, 255],
316    [70, 248, 132, 255],
317    [74, 248, 128, 255],
318    [78, 249, 125, 255],
319    [82, 250, 122, 255],
320    [85, 250, 118, 255],
321    [89, 251, 115, 255],
322    [93, 252, 111, 255],
323    [97, 252, 108, 255],
324    [101, 253, 105, 255],
325    [105, 253, 102, 255],
326    [109, 254, 98, 255],
327    [113, 254, 95, 255],
328    [117, 254, 92, 255],
329    [121, 254, 89, 255],
330    [125, 255, 86, 255],
331    [128, 255, 83, 255],
332    [132, 255, 81, 255],
333    [136, 255, 78, 255],
334    [139, 255, 75, 255],
335    [143, 255, 73, 255],
336    [146, 255, 71, 255],
337    [150, 254, 68, 255],
338    [153, 254, 66, 255],
339    [156, 254, 64, 255],
340    [159, 253, 63, 255],
341    [161, 253, 61, 255],
342    [164, 252, 60, 255],
343    [167, 252, 58, 255],
344    [169, 251, 57, 255],
345    [172, 251, 56, 255],
346    [175, 250, 55, 255],
347    [177, 249, 54, 255],
348    [180, 248, 54, 255],
349    [183, 247, 53, 255],
350    [185, 246, 53, 255],
351    [188, 245, 52, 255],
352    [190, 244, 52, 255],
353    [193, 243, 52, 255],
354    [195, 241, 52, 255],
355    [198, 240, 52, 255],
356    [200, 239, 52, 255],
357    [203, 237, 52, 255],
358    [205, 236, 52, 255],
359    [208, 234, 52, 255],
360    [210, 233, 53, 255],
361    [212, 231, 53, 255],
362    [215, 229, 53, 255],
363    [217, 228, 54, 255],
364    [219, 226, 54, 255],
365    [221, 224, 55, 255],
366    [223, 223, 55, 255],
367    [225, 221, 55, 255],
368    [227, 219, 56, 255],
369    [229, 217, 56, 255],
370    [231, 215, 57, 255],
371    [233, 213, 57, 255],
372    [235, 211, 57, 255],
373    [236, 209, 58, 255],
374    [238, 207, 58, 255],
375    [239, 205, 58, 255],
376    [241, 203, 58, 255],
377    [242, 201, 58, 255],
378    [244, 199, 58, 255],
379    [245, 197, 58, 255],
380    [246, 195, 58, 255],
381    [247, 193, 58, 255],
382    [248, 190, 57, 255],
383    [249, 188, 57, 255],
384    [250, 186, 57, 255],
385    [251, 184, 56, 255],
386    [251, 182, 55, 255],
387    [252, 179, 54, 255],
388    [252, 177, 54, 255],
389    [253, 174, 53, 255],
390    [253, 172, 52, 255],
391    [254, 169, 51, 255],
392    [254, 167, 50, 255],
393    [254, 164, 49, 255],
394    [254, 161, 48, 255],
395    [254, 158, 47, 255],
396    [254, 155, 45, 255],
397    [254, 153, 44, 255],
398    [254, 150, 43, 255],
399    [254, 147, 42, 255],
400    [254, 144, 41, 255],
401    [253, 141, 39, 255],
402    [253, 138, 38, 255],
403    [252, 135, 37, 255],
404    [252, 132, 35, 255],
405    [251, 129, 34, 255],
406    [251, 126, 33, 255],
407    [250, 123, 31, 255],
408    [249, 120, 30, 255],
409    [249, 117, 29, 255],
410    [248, 114, 28, 255],
411    [247, 111, 26, 255],
412    [246, 108, 25, 255],
413    [245, 105, 24, 255],
414    [244, 102, 23, 255],
415    [243, 99, 21, 255],
416    [242, 96, 20, 255],
417    [241, 93, 19, 255],
418    [240, 91, 18, 255],
419    [239, 88, 17, 255],
420    [237, 85, 16, 255],
421    [236, 83, 15, 255],
422    [235, 80, 14, 255],
423    [234, 78, 13, 255],
424    [232, 75, 12, 255],
425    [231, 73, 12, 255],
426    [229, 71, 11, 255],
427    [228, 69, 10, 255],
428    [226, 67, 10, 255],
429    [225, 65, 9, 255],
430    [223, 63, 8, 255],
431    [221, 61, 8, 255],
432    [220, 59, 7, 255],
433    [218, 57, 7, 255],
434    [216, 55, 6, 255],
435    [214, 53, 6, 255],
436    [212, 51, 5, 255],
437    [210, 49, 5, 255],
438    [208, 47, 5, 255],
439    [206, 45, 4, 255],
440    [204, 43, 4, 255],
441    [202, 42, 4, 255],
442    [200, 40, 3, 255],
443    [197, 38, 3, 255],
444    [195, 37, 3, 255],
445    [193, 35, 2, 255],
446    [190, 33, 2, 255],
447    [188, 32, 2, 255],
448    [185, 30, 2, 255],
449    [183, 29, 2, 255],
450    [180, 27, 1, 255],
451    [178, 26, 1, 255],
452    [175, 24, 1, 255],
453    [172, 23, 1, 255],
454    [169, 22, 1, 255],
455    [167, 20, 1, 255],
456    [164, 19, 1, 255],
457    [161, 18, 1, 255],
458    [158, 16, 1, 255],
459    [155, 15, 1, 255],
460    [152, 14, 1, 255],
461    [149, 13, 1, 255],
462    [146, 11, 1, 255],
463    [142, 10, 1, 255],
464    [139, 9, 2, 255],
465    [136, 8, 2, 255],
466    [133, 7, 2, 255],
467    [129, 6, 2, 255],
468    [126, 5, 2, 255],
469    [122, 4, 3, 255],
470];
471
472/// The 256-entry table of a built-in colormap.
473fn table_of(map: &Colormap) -> &'static [[u8; 4]; TABLE_LEN] {
474    static VIRIDIS: OnceLock<[[u8; 4]; TABLE_LEN]> = OnceLock::new();
475    static COOLWARM: OnceLock<[[u8; 4]; TABLE_LEN]> = OnceLock::new();
476    static GRAY: OnceLock<[[u8; 4]; TABLE_LEN]> = OnceLock::new();
477    match map {
478        Colormap::Turbo => &TURBO_TABLE,
479        Colormap::Viridis => VIRIDIS.get_or_init(|| build_table(viridis)),
480        Colormap::Coolwarm => COOLWARM.get_or_init(|| build_table(coolwarm)),
481        Colormap::Grayscale => GRAY.get_or_init(|| build_table(grayscale)),
482        // Custom colormaps sample their stops directly and never reach here.
483        Colormap::Custom(_) => unreachable!("custom colormaps do not use tables"),
484    }
485}
486
487fn build_table(f: fn(f64) -> [f64; 3]) -> [[u8; 4]; TABLE_LEN] {
488    let mut table = [[0u8; 4]; TABLE_LEN];
489    for (i, entry) in table.iter_mut().enumerate() {
490        let c = f(i as f64 / (TABLE_LEN - 1) as f64);
491        *entry = [to_u8(c[0]), to_u8(c[1]), to_u8(c[2]), 255];
492    }
493    table
494}
495
496fn to_u8(x: f64) -> u8 {
497    (x * 255.0).round().clamp(0.0, 255.0) as u8
498}
499
500fn grayscale(t: f64) -> [f64; 3] {
501    [t, t, t]
502}
503
504/// Matt Zucker's degree-6 polynomial fit of matplotlib's viridis.
505fn viridis(t: f64) -> [f64; 3] {
506    const C: [[f64; 3]; 7] = [
507        [
508            0.277_727_327_223_417_7,
509            0.005_407_344_544_966_578,
510            0.334_099_805_335_306_1,
511        ],
512        [
513            0.105_093_043_108_577_4,
514            1.404_613_529_898_575,
515            1.384_590_162_594_685,
516        ],
517        [
518            -0.330_861_828_725_556_3,
519            0.214_847_559_468_213,
520            0.095_095_163_028_236_59,
521        ],
522        [
523            -4.634_230_498_983_486,
524            -5.799_100_973_351_585,
525            -19.332_440_956_279_87,
526        ],
527        [
528            6.228_269_936_347_081,
529            14.179_933_366_805_09,
530            56.690_552_600_681_05,
531        ],
532        [
533            4.776_384_997_670_288,
534            -13.745_145_377_746_01,
535            -65.353_032_633_372_34,
536        ],
537        [
538            -5.435_455_855_934_631,
539            4.645_852_612_178_535,
540            26.312_435_249_583_2,
541        ],
542    ];
543    let mut out = [0.0; 3];
544    for (ch, o) in out.iter_mut().enumerate() {
545        let mut acc = 0.0;
546        for c in C.iter().rev() {
547            acc = acc * t + c[ch];
548        }
549        *o = acc.clamp(0.0, 1.0);
550    }
551    out
552}
553
554/// Moreland's diverging map, interpolated in Msh space.
555fn coolwarm(t: f64) -> [f64; 3] {
556    static ENDS: OnceLock<(Msh, Msh)> = OnceLock::new();
557    let (cool, warm) = ENDS.get_or_init(|| {
558        (
559            srgb_to_msh([59.0 / 255.0, 76.0 / 255.0, 192.0 / 255.0]),
560            srgb_to_msh([180.0 / 255.0, 4.0 / 255.0, 38.0 / 255.0]),
561        )
562    });
563    msh_to_srgb(msh_interpolate(*cool, *warm, t))
564}
565
566type Msh = [f64; 3];
567
568fn srgb_to_msh(c: [f64; 3]) -> Msh {
569    let lin = c.map(|u| {
570        if u <= 0.040_45 {
571            u / 12.92
572        } else {
573            ((u + 0.055) / 1.055).powf(2.4)
574        }
575    });
576    let x = 0.412_456_4 * lin[0] + 0.357_576_1 * lin[1] + 0.180_437_5 * lin[2];
577    let y = 0.212_672_9 * lin[0] + 0.715_152_2 * lin[1] + 0.072_175_0 * lin[2];
578    let z = 0.019_333_9 * lin[0] + 0.119_192_0 * lin[1] + 0.950_304_1 * lin[2];
579    let f = |u: f64| {
580        if u > 216.0 / 24389.0 {
581            u.cbrt()
582        } else {
583            (24389.0 / 27.0 * u + 16.0) / 116.0
584        }
585    };
586    let (fx, fy, fz) = (f(x / 0.950_47), f(y), f(z / 1.088_83));
587    let l = 116.0 * fy - 16.0;
588    let a = 500.0 * (fx - fy);
589    let b = 200.0 * (fy - fz);
590    let m = (l * l + a * a + b * b).sqrt();
591    let s = if m > 0.0 {
592        (l / m).clamp(-1.0, 1.0).acos()
593    } else {
594        0.0
595    };
596    let h = b.atan2(a);
597    [m, s, h]
598}
599
600fn msh_to_srgb(msh: Msh) -> [f64; 3] {
601    let [m, s, h] = msh;
602    let l = m * s.cos();
603    let a = m * s.sin() * h.cos();
604    let b = m * s.sin() * h.sin();
605    let fy = (l + 16.0) / 116.0;
606    let fx = fy + a / 500.0;
607    let fz = fy - b / 200.0;
608    let finv = |u: f64| {
609        if u * u * u > 216.0 / 24389.0 {
610            u * u * u
611        } else {
612            (116.0 * u - 16.0) * 27.0 / 24389.0
613        }
614    };
615    let (x, y, z) = (finv(fx) * 0.950_47, finv(fy), finv(fz) * 1.088_83);
616    let lin = [
617        3.240_454_2 * x - 1.537_138_5 * y - 0.498_531_4 * z,
618        -0.969_266_0 * x + 1.876_010_8 * y + 0.041_556_0 * z,
619        0.055_643_4 * x - 0.204_025_9 * y + 1.057_225_2 * z,
620    ];
621    lin.map(|u| {
622        let u = u.clamp(0.0, 1.0);
623        if u <= 0.003_130_8 {
624            12.92 * u
625        } else {
626            1.055 * u.powf(1.0 / 2.4) - 0.055
627        }
628    })
629}
630
631/// Moreland's hue adjustment for the unsaturated end of an interpolation.
632fn adjust_hue(sat: Msh, m_unsat: f64) -> f64 {
633    if sat[0] >= m_unsat {
634        return sat[2];
635    }
636    let spin = sat[1] * (m_unsat * m_unsat - sat[0] * sat[0]).sqrt() / (sat[0] * sat[1].sin());
637    if sat[2] > -std::f64::consts::FRAC_PI_3 {
638        sat[2] + spin
639    } else {
640        sat[2] - spin
641    }
642}
643
644fn msh_interpolate(mut a: Msh, mut b: Msh, mut t: f64) -> Msh {
645    if a[1] > 0.05 && b[1] > 0.05 && (a[2] - b[2]).abs() > std::f64::consts::FRAC_PI_3 {
646        let m_mid = a[0].max(b[0]).max(88.0);
647        if t < 0.5 {
648            b = [m_mid, 0.0, 0.0];
649            t *= 2.0;
650        } else {
651            a = [m_mid, 0.0, 0.0];
652            t = 2.0 * t - 1.0;
653        }
654    }
655    if a[1] < 0.05 && b[1] > 0.05 {
656        a[2] = adjust_hue(b, a[0]);
657    } else if b[1] < 0.05 && a[1] > 0.05 {
658        b[2] = adjust_hue(a, b[0]);
659    }
660    [
661        (1.0 - t) * a[0] + t * b[0],
662        (1.0 - t) * a[1] + t * b[1],
663        (1.0 - t) * a[2] + t * b[2],
664    ]
665}
666
667#[cfg(test)]
668mod tests {
669    use super::*;
670
671    fn rule(map: Colormap) -> ColorRule {
672        ColorRule {
673            map,
674            range: (0.0, 10.0),
675            below: OutOfRange::Clamp,
676            above: OutOfRange::Clamp,
677            center: None,
678            base: [9, 8, 7, 255],
679            nan: [1, 2, 3, 255],
680        }
681    }
682
683    fn near(a: [u8; 4], b: [u8; 4], tol: i32) -> bool {
684        a.iter()
685            .zip(&b)
686            .all(|(&x, &y)| (i32::from(x) - i32::from(y)).abs() <= tol)
687    }
688
689    #[test]
690    fn range_ends_and_center_map_to_table_positions() {
691        let r = rule(Colormap::Grayscale);
692        assert_eq!(color_of(&r, 0.0), [0, 0, 0, 255]);
693        assert_eq!(color_of(&r, 10.0), [255, 255, 255, 255]);
694        assert_eq!(color_of(&r, 5.0), [128, 128, 128, 255]);
695
696        // Center pinned off-middle: 2 maps to t = 0.5, 6 to t = 0.75.
697        let mut c = rule(Colormap::Grayscale);
698        c.center = Some(2.0);
699        assert_eq!(color_of(&c, 0.0), [0, 0, 0, 255]);
700        assert_eq!(color_of(&c, 2.0), [128, 128, 128, 255]);
701        assert_eq!(color_of(&c, 6.0), [191, 191, 191, 255]);
702        assert_eq!(color_of(&c, 10.0), [255, 255, 255, 255]);
703
704        // Built-in tables: opaque, distinct ends, expected anchors.
705        for map in [Colormap::Turbo, Colormap::Viridis, Colormap::Coolwarm] {
706            let r = rule(map);
707            let (a, b) = (color_of(&r, 0.0), color_of(&r, 10.0));
708            assert_eq!((a[3], b[3]), (255, 255));
709            assert_ne!(a, b);
710        }
711        assert!(near(
712            color_of(&rule(Colormap::Viridis), 0.0),
713            [68, 1, 84, 255],
714            6
715        ));
716        assert!(near(
717            color_of(&rule(Colormap::Viridis), 10.0),
718            [253, 231, 37, 255],
719            6
720        ));
721        assert!(near(
722            color_of(&rule(Colormap::Coolwarm), 0.0),
723            [59, 76, 192, 255],
724            3
725        ));
726        assert!(near(
727            color_of(&rule(Colormap::Coolwarm), 10.0),
728            [180, 4, 38, 255],
729            3
730        ));
731        let mid = color_of(&rule(Colormap::Coolwarm), 5.0);
732        assert!(near(mid, [221, 221, 221, 255], 4), "{mid:?}");
733    }
734
735    #[test]
736    fn turbo_matches_published_table() {
737        // The embedded table is Google's published 256-entry Turbo LUT
738        // (Mikhailov, 2019): dark end (48, 18, 59), bright end (122, 4, 3).
739        let r = rule(Colormap::Turbo);
740        assert_eq!(color_of(&r, 0.0), [48, 18, 59, 255]);
741        assert_eq!(color_of(&r, 10.0), [122, 4, 3, 255]);
742    }
743
744    #[test]
745    fn viridis_midpoints_match_reference() {
746        // Mid-points against matplotlib's viridis table.
747        let r = rule(Colormap::Viridis);
748        assert!(
749            near(color_of(&r, 2.5), [59, 82, 139, 255], 6),
750            "{:?}",
751            color_of(&r, 2.5)
752        );
753        assert!(
754            near(color_of(&r, 5.0), [33, 145, 140, 255], 6),
755            "{:?}",
756            color_of(&r, 5.0)
757        );
758        assert!(
759            near(color_of(&r, 7.5), [94, 201, 98, 255], 6),
760            "{:?}",
761            color_of(&r, 7.5)
762        );
763    }
764
765    #[test]
766    fn below_and_above_are_independent() {
767        let mut r = rule(Colormap::Grayscale);
768        r.below = OutOfRange::Clamp;
769        r.above = OutOfRange::Base;
770        assert_eq!(color_of(&r, -5.0), [0, 0, 0, 255]);
771        assert_eq!(color_of(&r, 15.0), r.base);
772        r.below = OutOfRange::Base;
773        r.above = OutOfRange::Clamp;
774        assert_eq!(color_of(&r, -5.0), r.base);
775        assert_eq!(color_of(&r, 15.0), [255, 255, 255, 255]);
776    }
777
778    #[test]
779    fn infinities_follow_below_and_above() {
780        // +inf follows `above`, -inf follows `below`.
781        let mut r = rule(Colormap::Grayscale);
782        r.below = OutOfRange::Base;
783        r.above = OutOfRange::Base;
784        assert_eq!(color_of(&r, f32::INFINITY), r.base);
785        assert_eq!(color_of(&r, f32::NEG_INFINITY), r.base);
786        r.above = OutOfRange::Clamp;
787        assert_eq!(color_of(&r, f32::INFINITY), [255, 255, 255, 255]);
788        assert_eq!(color_of(&r, f32::NEG_INFINITY), r.base);
789    }
790
791    #[test]
792    fn nan_is_nan_color() {
793        let r = rule(Colormap::Viridis);
794        assert_eq!(color_of(&r, f32::NAN), r.nan);
795    }
796
797    #[test]
798    fn custom_stops_interpolate() {
799        let r = rule(Colormap::Custom(vec![
800            (0.0, [0, 0, 0, 255]),
801            (0.5, [100, 200, 50, 255]),
802            (1.0, [200, 200, 250, 255]),
803        ]));
804        assert_eq!(color_of(&r, 0.0), [0, 0, 0, 255]);
805        assert_eq!(color_of(&r, 2.5), [50, 100, 25, 255]);
806        assert_eq!(color_of(&r, 5.0), [100, 200, 50, 255]);
807        assert_eq!(color_of(&r, 7.5), [150, 200, 150, 255]);
808        assert_eq!(color_of(&r, 10.0), [200, 200, 250, 255]);
809        assert!(r.validate().is_ok());
810    }
811
812    #[test]
813    fn invalid_rules_rejected() {
814        let ok = rule(Colormap::Turbo);
815        assert!(ok.validate().is_ok());
816        let mut r = ok.clone();
817        r.range = (1.0, 1.0);
818        assert!(r.validate().is_err());
819        r.range = (2.0, 1.0);
820        assert!(r.validate().is_err());
821        r.range = (f32::NEG_INFINITY, 1.0);
822        assert!(r.validate().is_err());
823        r.range = (0.0, f32::NAN);
824        assert!(r.validate().is_err());
825        let mut r = ok.clone();
826        r.center = Some(11.0);
827        assert!(r.validate().is_err());
828        r.center = Some(-1.0);
829        assert!(r.validate().is_err());
830        r.center = Some(f32::NAN);
831        assert!(r.validate().is_err());
832        r.center = Some(10.0);
833        assert!(r.validate().is_ok());
834        let bad = |stops: Vec<(f32, [u8; 4])>| rule(Colormap::Custom(stops)).validate().is_err();
835        assert!(bad(vec![(0.0, [0; 4])]));
836        assert!(bad(vec![(0.5, [0; 4]), (0.2, [0; 4])]));
837        assert!(bad(vec![(0.5, [0; 4]), (0.5, [0; 4])]));
838        assert!(bad(vec![(-0.1, [0; 4]), (0.5, [0; 4])]));
839        assert!(bad(vec![(0.0, [0; 4]), (1.5, [0; 4])]));
840    }
841}