colorimetry/observer/
spectral_locus.rs

1// SPDX-License-Identifier: Apache-2.0 OR MIT
2// Copyright (c) 2024-2025, Harbers Bik LLC
3
4//! Spectral Locus module provides functionality to work with the spectral locus in colorimetry.
5//!
6//! The spectral locus represents the boundary of visible colors in a chromaticity diagram,
7//! formed by plotting monochromatic light sources of different wavelengths.
8//!
9//! The module provides:
10//! - `SpectralLocus`: Main struct representing the spectral locus
11//! - `SpectralLocusIterator`: Iterator for traversing spectral locus points
12//!
13//! The spectral locus is represented as a closed polygon, where each point corresponds
14//! to the chromaticity coordinates of monochromatic light at different wavelengths.
15
16use std::ops::RangeBounds;
17
18use geo::{Contains, Polygon};
19
20use crate::{
21    math::QuadraticCurve,
22    observer::Observer,
23    spectrum::{NS, SPECTRUM_WAVELENGTH_RANGE},
24    xyz::XYZ,
25};
26
27/// Represents the spectral locus, which includes:
28/// - An observer for spectral data.
29/// - A vector of tuples where each tuple contains:
30///   - A wavelength (usize).
31///   - A chromaticity coordinate ([f64; 2]).
32pub struct SpectralLocus {
33    observer: Observer,
34    polygon: Polygon,
35}
36
37impl SpectralLocus {
38    pub fn new(observer: Observer) -> Self {
39        let obs_data = observer.data().data;
40        let mut v = Vec::with_capacity(NS + 1);
41        for i in 0..NS {
42            let xyz = obs_data.column(i).into();
43            let chromaticity = XYZ::from_vec(xyz, observer).chromaticity();
44            v.push(chromaticity.to_array());
45        }
46        v.push(v[0]); // close the polygon
47        let polygon = vec_to_polygon(v);
48        SpectralLocus { observer, polygon }
49    }
50
51    pub fn contains(&self, point: [f64; 2]) -> bool {
52        self.polygon.contains(&geo::Point::new(point[0], point[1]))
53    }
54
55    pub fn observer(&self) -> Observer {
56        self.observer
57    }
58
59    /// Iterate over the spectral locus points over a given wavelength range and step size producing
60    /// the spectral locus coordinates (f64, f64) and slope angles f64 in radians.
61    /// The slope is calculated using a quadratic curve fit for interior points, and a simple
62    /// difference for the endpoints.
63    ///
64    /// This iterator is useful for applications that require not only the locus points,
65    /// but also the direction of the locus at each point, for example when plotting tick marks along the spectral locus
66    /// in plots, to find its boundaries.
67    ///
68    /// # Item
69    /// Each iteration yields a tuple:
70    /// - `((f64, f64), f64)`
71    ///   - The chromaticity coordinate (x, y) of the locus point.
72    ///   - The slope angle (in radians) at that point.
73    pub fn iter_range_with_slope<'a>(
74        &'a self,
75        range: impl RangeBounds<usize>,
76        step: usize,
77    ) -> SpectralLocusSlopeIterator<'a> {
78        let min_wavelength = *SPECTRUM_WAVELENGTH_RANGE.start();
79        let clip = |s, lim| {
80            if s > min_wavelength {
81                s - min_wavelength
82            } else {
83                lim
84            }
85        };
86        let start = match range.start_bound() {
87            std::ops::Bound::Included(&s) => clip(s, 0),
88            std::ops::Bound::Excluded(&s) => clip(s + 1, 0),
89            std::ops::Bound::Unbounded => 0,
90        };
91        let end = match range.end_bound() {
92            std::ops::Bound::Included(&e) => clip(e + 1, NS - 1),
93            std::ops::Bound::Excluded(&e) => clip(e, NS - 1),
94            std::ops::Bound::Unbounded => NS - 1,
95        };
96
97        if end <= start {
98            panic!("End of range must be greater than start");
99        }
100
101        SpectralLocusSlopeIterator {
102            locus: self,
103            index: start,
104            step,
105            end,
106        }
107    }
108
109    pub fn polygon(&self) -> &Polygon {
110        &self.polygon
111    }
112}
113
114pub struct SpectralLocusIterator<'a> {
115    locus: &'a SpectralLocus,
116    index: usize,
117}
118
119impl<'a> IntoIterator for &'a SpectralLocus {
120    type Item = (f64, f64);
121    type IntoIter = SpectralLocusIterator<'a>;
122
123    fn into_iter(self) -> Self::IntoIter {
124        SpectralLocusIterator {
125            locus: self,
126            index: 0,
127        }
128    }
129}
130
131impl<'a> Iterator for SpectralLocusIterator<'a> {
132    type Item = (f64, f64);
133
134    fn next(&mut self) -> Option<Self::Item> {
135        if self.index < self.locus.polygon.exterior().points().count() {
136            let coord = self.locus.polygon.exterior()[self.index];
137            self.index += 1;
138            Some((coord.x, coord.y))
139        } else {
140            None
141        }
142    }
143}
144
145pub struct SpectralLocusSlopeIterator<'a> {
146    locus: &'a SpectralLocus,
147    index: usize,
148    step: usize,
149    end: usize,
150}
151
152impl Iterator for SpectralLocusSlopeIterator<'_> {
153    type Item = ((f64, f64), f64);
154
155    fn next(&mut self) -> Option<Self::Item> {
156        if self.index < self.end {
157            let coord = self.locus.polygon.exterior()[self.index];
158            let angle = match self.index {
159                0 => {
160                    let coord_plus = self.locus.polygon.exterior()[1];
161                    (coord_plus.y - coord.y).atan2(coord_plus.x - coord.x)
162                }
163                i if i < NS - 1 => {
164                    let coord_min = self.locus.polygon.exterior()[i - 1];
165                    let coord_plus = self.locus.polygon.exterior()[i + 1];
166                    let qurve = QuadraticCurve::new(
167                        (coord_min.x, coord_min.y),
168                        (coord.x, coord.y),
169                        (coord_plus.x, coord_plus.y),
170                    )
171                    .unwrap();
172                    qurve.slope_angle(0.5)
173                }
174                n if n == NS - 1 => {
175                    let coord_min = self.locus.polygon.exterior()[n - 1];
176                    (coord.y - coord_min.y).atan2(coord.x - coord_min.x)
177                }
178                _ => 0.0, // Should not be reached given the if self.index < self.end condition
179            };
180            self.index += self.step;
181            Some(((coord.x, coord.y), angle))
182        } else {
183            None
184        }
185    }
186}
187
188impl<'a> DoubleEndedIterator for SpectralLocusSlopeIterator<'a> {
189    fn next_back(&mut self) -> Option<<Self as Iterator>::Item> {
190        if self.index < self.end && self.end > 0 {
191            self.end = self.end.saturating_sub(self.step);
192            let coord = self.locus.polygon.exterior()[self.end];
193            let angle = match self.end {
194                0 => {
195                    let coord_plus = self.locus.polygon.exterior()[1];
196                    (coord_plus.y - coord.y).atan2(coord_plus.x - coord.x)
197                }
198                i if i < NS - 1 => {
199                    let coord_min = self.locus.polygon.exterior()[i - 1];
200                    let coord_plus = self.locus.polygon.exterior()[i + 1];
201                    let qurve = QuadraticCurve::new(
202                        (coord_min.x, coord_min.y),
203                        (coord.x, coord.y),
204                        (coord_plus.x, coord_plus.y),
205                    )
206                    .unwrap();
207                    qurve.slope_angle(0.5)
208                }
209                n if n == NS - 1 => {
210                    let coord_min = self.locus.polygon.exterior()[n - 1];
211                    (coord.y - coord_min.y).atan2(coord.x - coord_min.x)
212                }
213                _ => 0.0,
214            };
215            Some(((coord.x, coord.y), angle))
216        } else {
217            None
218        }
219    }
220}
221
222// make sure Vec is closed
223use geo::{Coord, LineString};
224fn vec_to_polygon(coords: Vec<[f64; 2]>) -> Polygon {
225    let linestring = LineString::from(
226        coords
227            .into_iter()
228            .map(|[x, y]| Coord { x, y })
229            .collect::<Vec<_>>(),
230    );
231
232    Polygon::new(linestring, vec![]) // no interior rings
233}
234
235#[cfg(test)]
236mod test {
237    use approx::assert_abs_diff_eq;
238
239    use super::*;
240    use crate::observer::Observer::Cie1931;
241
242    #[test]
243    fn test_spectral_locus() {
244        let locus = SpectralLocus::new(Cie1931);
245        assert_eq!(locus.observer(), Cie1931);
246        assert!(locus.contains([0.3, 0.3]));
247        assert!(!locus.contains([0.05, 0.05]));
248        assert!(!locus.contains([0.7, 0.7]));
249
250        let points: Vec<(f64, f64)> = locus.into_iter().collect();
251        assert_eq!(points.len(), NS + 1); // closed polygon
252    }
253
254    #[test]
255    fn test_spectral_locus_iterator() {
256        let locus = SpectralLocus::new(Cie1931);
257        let mut iter = locus.into_iter();
258        let point = iter.next().unwrap();
259        assert_abs_diff_eq!(point.0, 0.17411, epsilon = 0.00005);
260        assert_abs_diff_eq!(point.1, 0.00496, epsilon = 0.00005);
261        let last = iter.last().unwrap();
262        assert_abs_diff_eq!(last.0, 0.17411, epsilon = 0.00005);
263        assert_abs_diff_eq!(last.1, 0.00496, epsilon = 0.00005);
264    }
265
266    #[test]
267    fn test_spectral_locus_slope_iterator() {
268        let locus = SpectralLocus::new(Cie1931);
269        let mut iter = locus.iter_range_with_slope(400..700, 10);
270
271        // Check the first point
272        let (point, angle) = iter.next().unwrap();
273        assert_abs_diff_eq!(point.0, 0.173336, epsilon = 0.00005);
274        assert_abs_diff_eq!(point.1, 0.0047967, epsilon = 0.00005);
275        assert_abs_diff_eq!(angle, -2.84106, epsilon = 0.00005); // slope at first point
276
277        // Check the last point
278        let (point, angle) = iter.next_back().unwrap();
279        assert_abs_diff_eq!(point.0, 0.734390, epsilon = 0.00005);
280        assert_abs_diff_eq!(point.1, 0.265609, epsilon = 0.00005);
281        assert_abs_diff_eq!(angle, -0.785398, epsilon = 0.00005); // slope at last point
282    }
283}
colorimetry/observer/spectral_locus.rs

colorimetry/observer/
spectral_locus.rs
