1use crate::color::TristimulusObserver;
2use crate::cri::{
3 spds_to_ciera, spds_to_ciera_result, spds_to_ciera_special, spds_to_cierf,
4 spds_to_cierf_result, spds_to_cierf_special, spds_to_cierg, spds_to_ies_tm30_result,
5 spds_to_iesrf, spds_to_iesrf_result, spds_to_iesrf_special, spds_to_iesrg, spds_to_tm30_result,
6 CieRaResult, CieRfResult, Tm30Result,
7};
8use crate::error::{LuxError, LuxResult};
9use crate::photometry::{spd_to_power, PowerType};
10use crate::spectral_mismatch::{spectral_mismatch_correction_factors, spectral_mismatch_f1primes};
11
12#[derive(Debug, Clone, Copy, PartialEq)]
13pub struct WavelengthGrid {
14 pub start: f64,
15 pub end: f64,
16 pub step: f64,
17}
18
19impl WavelengthGrid {
20 pub fn new(start: f64, end: f64, step: f64) -> LuxResult<Self> {
21 if !start.is_finite() || !end.is_finite() || !step.is_finite() || step <= 0.0 || end < start
22 {
23 return Err(LuxError::InvalidGridSpec);
24 }
25 Ok(Self { start, end, step })
26 }
27}
28
29#[derive(Debug, Clone, PartialEq)]
30pub(crate) struct SingleSpectrum {
31 wavelengths: Vec<f64>,
32 values: Vec<f64>,
33}
34
35#[derive(Debug, Clone, PartialEq)]
36pub struct Spectrum {
37 wavelengths: Vec<f64>,
38 spectra: Vec<Vec<f64>>,
39}
40
41#[derive(Debug, Clone, Copy, PartialEq)]
42pub enum SpectrumNormalization {
43 Max(f64),
44 Area(f64),
45 Lambda(f64),
46 Radiometric(f64),
47 Photometric(f64),
48 Quantal(f64),
49}
50
51#[doc(hidden)]
52pub trait IntoSpectrumRows {
53 fn into_rows(self) -> Vec<Vec<f64>>;
54}
55
56impl IntoSpectrumRows for Vec<f64> {
57 fn into_rows(self) -> Vec<Vec<f64>> {
58 vec![self]
59 }
60}
61
62impl IntoSpectrumRows for Vec<Vec<f64>> {
63 fn into_rows(self) -> Vec<Vec<f64>> {
64 self
65 }
66}
67
68impl SingleSpectrum {
69 fn as_batch(&self) -> LuxResult<Spectrum> {
70 Spectrum::new(self.wavelengths.clone(), self.values.clone())
71 }
72
73 pub fn new(wavelengths: Vec<f64>, values: Vec<f64>) -> LuxResult<Self> {
74 if wavelengths.is_empty() || values.is_empty() {
75 return Err(LuxError::EmptyInput);
76 }
77 if wavelengths.len() != values.len() {
78 return Err(LuxError::MismatchedLengths {
79 wavelengths: wavelengths.len(),
80 values: values.len(),
81 });
82 }
83 if wavelengths
84 .windows(2)
85 .any(|pair| !(pair[1].is_finite() && pair[0].is_finite() && pair[1] > pair[0]))
86 {
87 return Err(LuxError::NonMonotonicWavelengths);
88 }
89 Ok(Self {
90 wavelengths,
91 values,
92 })
93 }
94
95 pub fn values(&self) -> &[f64] {
96 &self.values
97 }
98
99 pub fn spacing(&self) -> LuxResult<Vec<f64>> {
100 getwld(&self.wavelengths)
101 }
102
103 pub fn interpolate_linear(&self, target_wavelengths: &[f64]) -> LuxResult<Self> {
104 if target_wavelengths.is_empty() {
105 return Err(LuxError::EmptyInput);
106 }
107 if target_wavelengths
108 .windows(2)
109 .any(|pair| !(pair[1].is_finite() && pair[0].is_finite() && pair[1] > pair[0]))
110 {
111 return Err(LuxError::NonMonotonicWavelengths);
112 }
113
114 let mut values = Vec::with_capacity(target_wavelengths.len());
115 for &target in target_wavelengths {
116 values.push(self.interpolate_one_linear(target));
117 }
118 SingleSpectrum::new(target_wavelengths.to_vec(), values)
119 }
120
121 pub fn cie_interp_linear(
122 &self,
123 target_wavelengths: &[f64],
124 negative_values_allowed: bool,
125 ) -> LuxResult<Self> {
126 let mut interpolated = self.interpolate_linear(target_wavelengths)?;
127 if !negative_values_allowed {
128 for value in &mut interpolated.values {
129 if *value < 0.0 {
130 *value = 0.0;
131 }
132 }
133 }
134 Ok(interpolated)
135 }
136
137 fn interpolate_one_linear(&self, target: f64) -> f64 {
138 let wavelengths = &self.wavelengths;
139 let values = &self.values;
140
141 if target <= wavelengths[0] {
142 return linear_segment(wavelengths[0], values[0], wavelengths[1], values[1], target);
143 }
144 if target >= wavelengths[wavelengths.len() - 1] {
145 let last = wavelengths.len() - 1;
146 return linear_segment(
147 wavelengths[last - 1],
148 values[last - 1],
149 wavelengths[last],
150 values[last],
151 target,
152 );
153 }
154
155 let idx = wavelengths.partition_point(|wavelength| *wavelength < target);
156 if wavelengths[idx] == target {
157 values[idx]
158 } else {
159 linear_segment(
160 wavelengths[idx - 1],
161 values[idx - 1],
162 wavelengths[idx],
163 values[idx],
164 target,
165 )
166 }
167 }
168
169 pub fn normalize(
170 &self,
171 mode: SpectrumNormalization,
172 observer: Option<&TristimulusObserver>,
173 ) -> LuxResult<Self> {
174 let scale = match mode {
175 SpectrumNormalization::Max(target_max) => {
176 target_max
177 / self
178 .values
179 .iter()
180 .copied()
181 .fold(f64::NEG_INFINITY, f64::max)
182 }
183 SpectrumNormalization::Area(target_area) => {
184 let area: f64 = self
185 .values
186 .iter()
187 .zip(self.spacing()?.iter())
188 .map(|(value, dl)| value * dl)
189 .sum();
190 target_area / area
191 }
192 SpectrumNormalization::Lambda(target_wavelength) => {
193 let (index, _) = self
194 .wavelengths
195 .iter()
196 .enumerate()
197 .map(|(index, wavelength)| (index, (wavelength - target_wavelength).abs()))
198 .min_by(|(_, lhs), (_, rhs)| lhs.partial_cmp(rhs).unwrap())
199 .ok_or(LuxError::EmptyInput)?;
200 1.0 / self.values[index]
201 }
202 SpectrumNormalization::Radiometric(target_power) => {
203 target_power / spd_to_power(&self.as_batch()?, PowerType::Radiometric, None)?
204 }
205 SpectrumNormalization::Photometric(target_power) => {
206 target_power / spd_to_power(&self.as_batch()?, PowerType::Photometric, observer)?
207 }
208 SpectrumNormalization::Quantal(target_power) => {
209 target_power / spd_to_power(&self.as_batch()?, PowerType::Quantal, None)?
210 }
211 };
212
213 SingleSpectrum::new(
214 self.wavelengths.clone(),
215 self.values
216 .iter()
217 .map(|value| value * scale)
218 .collect::<Vec<_>>(),
219 )
220 }
221}
222
223impl Spectrum {
224 pub fn new<T: IntoSpectrumRows>(wavelengths: Vec<f64>, spectra: T) -> LuxResult<Self> {
225 let spectra = spectra.into_rows();
226 if wavelengths.is_empty() || spectra.is_empty() {
227 return Err(LuxError::EmptyInput);
228 }
229 if wavelengths
230 .windows(2)
231 .any(|pair| !(pair[1].is_finite() && pair[0].is_finite() && pair[1] > pair[0]))
232 {
233 return Err(LuxError::NonMonotonicWavelengths);
234 }
235 for values in &spectra {
236 if values.len() != wavelengths.len() {
237 return Err(LuxError::MismatchedLengths {
238 wavelengths: wavelengths.len(),
239 values: values.len(),
240 });
241 }
242 }
243
244 Ok(Self {
245 wavelengths,
246 spectra,
247 })
248 }
249
250 pub fn wavelengths(&self) -> &[f64] {
251 &self.wavelengths
252 }
253
254 pub fn values(&self) -> &[f64] {
255 self.spectra.first().map(Vec::as_slice).unwrap_or(&[])
256 }
257
258 pub fn spectra(&self) -> &[Vec<f64>] {
259 &self.spectra
260 }
261
262 pub fn spectrum_count(&self) -> usize {
263 self.spectra.len()
264 }
265
266 pub fn wavelength_count(&self) -> usize {
267 self.wavelengths.len()
268 }
269
270 pub fn spacing(&self) -> LuxResult<Vec<f64>> {
271 getwld(&self.wavelengths)
272 }
273
274 pub fn into_vec(self) -> Vec<Vec<f64>> {
275 self.spectra
276 }
277
278 pub fn interpolate_linear(&self, target_wavelengths: &[f64]) -> LuxResult<Self> {
279 let mut spectra = Vec::with_capacity(self.spectra.len());
280 for values in &self.spectra {
281 let spectrum = SingleSpectrum::new(self.wavelengths.clone(), values.clone())?;
282 let interpolated = spectrum.interpolate_linear(target_wavelengths)?;
283 spectra.push(interpolated.values().to_vec());
284 }
285 Spectrum::new(target_wavelengths.to_vec(), spectra)
286 }
287
288 pub fn cie_interp_linear(
289 &self,
290 target_wavelengths: &[f64],
291 negative_values_allowed: bool,
292 ) -> LuxResult<Self> {
293 let mut spectra = Vec::with_capacity(self.spectra.len());
294 for values in &self.spectra {
295 let spectrum = SingleSpectrum::new(self.wavelengths.clone(), values.clone())?;
296 let interpolated =
297 spectrum.cie_interp_linear(target_wavelengths, negative_values_allowed)?;
298 spectra.push(interpolated.values().to_vec());
299 }
300 Spectrum::new(target_wavelengths.to_vec(), spectra)
301 }
302
303 pub fn normalize(
304 &self,
305 mode: SpectrumNormalization,
306 observer: Option<&TristimulusObserver>,
307 ) -> LuxResult<Self> {
308 self.normalize_each(&[mode], observer)
309 }
310
311 pub fn normalize_each(
312 &self,
313 modes: &[SpectrumNormalization],
314 observer: Option<&TristimulusObserver>,
315 ) -> LuxResult<Self> {
316 if modes.is_empty() {
317 return Err(LuxError::EmptyInput);
318 }
319
320 let mut spectra = Vec::with_capacity(self.spectra.len());
321 for (index, values) in self.spectra.iter().enumerate() {
322 let mode = modes.get(index).copied().unwrap_or(modes[0]);
323 let spectrum = SingleSpectrum::new(self.wavelengths.clone(), values.clone())?;
324 let normalized = spectrum.normalize(mode, observer)?;
325 spectra.push(normalized.values().to_vec());
326 }
327
328 Spectrum::new(self.wavelengths.clone(), spectra)
329 }
330
331 pub fn spd_to_xyz(
332 &self,
333 observer: &TristimulusObserver,
334 relative: bool,
335 ) -> LuxResult<Vec<[f64; 3]>> {
336 let wavelengths = self.wavelengths();
337 let x_bar = observer.x_bar_spectrum()?.interpolate_linear(wavelengths)?;
338 let y_bar = observer.vl_spectrum()?.interpolate_linear(wavelengths)?;
339 let z_bar = observer.z_bar_spectrum()?.interpolate_linear(wavelengths)?;
340
341 self.spectra
342 .iter()
343 .map(|values| {
344 let spectrum = Spectrum::new(wavelengths.to_vec(), values.clone())?;
345 crate::photometry::integrate_xyz(
346 &spectrum,
347 x_bar.values(),
348 y_bar.values(),
349 z_bar.values(),
350 observer.k,
351 relative,
352 )
353 })
354 .collect()
355 }
356
357 pub fn spd_to_ler(&self, observer: &TristimulusObserver) -> LuxResult<Vec<f64>> {
358 self.spectra
359 .iter()
360 .map(|values| {
361 let spectrum = Spectrum::new(self.wavelengths.to_vec(), values.clone())?;
362 crate::photometry::spd_to_ler(&spectrum, observer)
363 })
364 .collect()
365 }
366
367 pub fn spd_to_ciera(&self) -> LuxResult<Vec<f64>> {
368 spds_to_ciera(self)
369 }
370
371 pub fn spd_to_ciera_special(&self) -> LuxResult<Vec<Vec<f64>>> {
372 spds_to_ciera_special(self)
373 }
374
375 pub fn spd_to_ciera_result(&self) -> LuxResult<Vec<CieRaResult>> {
376 spds_to_ciera_result(self)
377 }
378
379 pub fn spd_to_cierf(&self) -> LuxResult<Vec<f64>> {
380 spds_to_cierf(self)
381 }
382
383 pub fn spd_to_iesrf(&self) -> LuxResult<Vec<f64>> {
384 spds_to_iesrf(self)
385 }
386
387 pub fn spd_to_cierg(&self) -> LuxResult<Vec<f64>> {
388 spds_to_cierg(self)
389 }
390
391 pub fn spd_to_iesrg(&self) -> LuxResult<Vec<f64>> {
392 spds_to_iesrg(self)
393 }
394
395 pub fn spd_to_cierf_special(&self) -> LuxResult<Vec<Vec<f64>>> {
396 spds_to_cierf_special(self)
397 }
398
399 pub fn spd_to_iesrf_special(&self) -> LuxResult<Vec<Vec<f64>>> {
400 spds_to_iesrf_special(self)
401 }
402
403 pub fn spd_to_cierf_result(&self) -> LuxResult<Vec<CieRfResult>> {
404 spds_to_cierf_result(self)
405 }
406
407 pub fn spd_to_iesrf_result(&self) -> LuxResult<Vec<CieRfResult>> {
408 spds_to_iesrf_result(self)
409 }
410
411 pub fn spd_to_tm30_result(&self) -> LuxResult<Vec<Tm30Result>> {
412 spds_to_tm30_result(self)
413 }
414
415 pub fn spd_to_ies_tm30_result(&self) -> LuxResult<Vec<Tm30Result>> {
416 spds_to_ies_tm30_result(self)
417 }
418
419 pub fn spectral_mismatch_f1primes(
420 &self,
421 calibration_illuminant: &Spectrum,
422 target_responsivity: &Spectrum,
423 ) -> LuxResult<Vec<f64>> {
424 spectral_mismatch_f1primes(self, calibration_illuminant, target_responsivity)
425 }
426
427 pub fn spectral_mismatch_correction_factors(
428 &self,
429 detectors: &Spectrum,
430 calibration_illuminant: &Spectrum,
431 target_responsivity: &Spectrum,
432 ) -> LuxResult<Vec<Vec<f64>>> {
433 spectral_mismatch_correction_factors(
434 self,
435 detectors,
436 calibration_illuminant,
437 target_responsivity,
438 )
439 }
440}
441
442fn linear_segment(x0: f64, y0: f64, x1: f64, y1: f64, x: f64) -> f64 {
443 y0 + (y1 - y0) * ((x - x0) / (x1 - x0))
444}
445
446pub fn getwlr(grid: WavelengthGrid) -> LuxResult<Vec<f64>> {
447 let mut wavelengths = Vec::new();
448 let mut current = grid.start;
449 let epsilon = grid.step.abs() * 1e-9;
450
451 while current <= grid.end + epsilon {
452 wavelengths.push(current);
453 current += grid.step;
454 }
455
456 if wavelengths.is_empty() {
457 return Err(LuxError::InvalidGridSpec);
458 }
459
460 Ok(wavelengths)
461}
462
463pub fn getwld(wavelengths: &[f64]) -> LuxResult<Vec<f64>> {
464 if wavelengths.is_empty() {
465 return Err(LuxError::EmptyInput);
466 }
467 if wavelengths.len() == 1 {
468 return Ok(vec![0.0]);
469 }
470 if wavelengths
471 .windows(2)
472 .any(|pair| !(pair[1].is_finite() && pair[0].is_finite() && pair[1] > pair[0]))
473 {
474 return Err(LuxError::NonMonotonicWavelengths);
475 }
476
477 let diffs: Vec<f64> = wavelengths
478 .windows(2)
479 .map(|pair| pair[1] - pair[0])
480 .collect::<Vec<_>>();
481 let mut spacing = Vec::with_capacity(wavelengths.len());
482 spacing.push(diffs[0]);
483 for idx in 1..wavelengths.len() - 1 {
484 spacing.push((diffs[idx - 1] + diffs[idx]) / 2.0);
485 }
486 spacing.push(*diffs.last().unwrap());
487 Ok(spacing)
488}