colorimetry 0.0.9

// SPDX-License-Identifier: Apache-2.0 OR MIT
// Copyright (c) 2024-2025, Harbers Bik LLC

//! Standard colorimetric observers
//!
//! A *standard observer* models average human color vision as three *color-matching functions*
//! (CMFs) — x̄(λ), ȳ(λ), z̄(λ) — derived from psychophysical experiments in which observers
//! matched test colors by adjusting three narrow-band primaries.  Integrating any spectral power
//! distribution (SPD) against these curves yields the *XYZ tristimulus values* that underpin all
//! CIE color metrics: CCT, CRI, CIELAB, TM-30 Rf, and more.
//!
//! # Available observers
//!
//! | Variant | Field of view | Basis | Primary use |
//! |---------|--------------|-------|-------------|
//! | [`Cie1931`](Observer::Cie1931) | 2° (foveal) | Wright & Guild 1931 | Default for most colorimetry |
//! | [`Cie1964`](Observer::Cie1964) | 10° (peripheral) | Stiles & Burch 1964 | CIE 224:2017 / TM-30; large fields |
//! | [`Cie2015`](Observer::Cie2015) | 2° | Cone fundamentals (Stockman & Sharpe) | Higher accuracy, especially in blue |
//! | [`Cie2015_10`](Observer::Cie2015_10) | 10° | Cone fundamentals (Stockman & Sharpe) | Higher accuracy, wide fields |
//!
//! The **CIE 1931 2° observer** dominates industry practice — sRGB, ICC profiles, and CIE
//! CRI Ra all specify their reference calculations in terms of it.  The **CIE 1964 10°
//! observer** is used by CIE 224:2017 / ANSI/IES TM-30 for colour fidelity calculations,
//! and is also preferred whenever the viewed area subtends more than roughly 4° at the eye
//! (larger than a business card held at arm’s length).  The **CIE 2015** observers construct
//! the CMFs as linear transforms of the Stockman & Sharpe (2000) cone fundamentals —
//! psychophysical estimates of L, M, S cone spectral sensitivities derived from improved
//! color-matching experiments — which corrects an inaccuracy in the short-wavelength region
//! of the 1931 functions.
//!
//! # Background: color-matching experiments
//!
//! Human color vision relies on three cone types — L (“red”), M (“green”), S (“blue”) — each
//! with a distinct spectral sensitivity.  The CIE quantified these indirectly: observers
//! adjusted monochromatic primaries until the mixture appeared identical to a test stimulus
//! at each wavelength.  Because some stimuli required adding light to the *test* side (a
//! “negative” primary contribution), the 1931 committee applied a linear transformation to a
//! set of *imaginary* primaries, ensuring x̄(λ), ȳ(λ), z̄(λ) are everywhere non-negative.
//! The ȳ function was additionally chosen to equal the photopic luminous efficiency V(λ),
//! so that Y always represents luminance (CIE 15:2018, §3).
//!
//! # Example
//!
//! ```
//! use colorimetry::observer::Observer;
//!
//! // D65 white point under the CIE 1931 2° observer — Y is 100 by convention
//! let xyz = Observer::Cie1931.xyz_d65();
//! assert!((xyz.to_array()[1] - 100.0).abs() < 0.1);
//! ```
//!
//! # Reference
//!
//! CIE 15:2018 *Colorimetry*, 4th ed., §3 (color-matching functions and standard observers).

pub mod observer_data;

mod rgbxyz;

mod optimal;
pub use optimal::OptimalColors;

mod spectral_locus;
pub use spectral_locus::SpectralLocus;

use crate::{
    cam::{CieCam02, CieCam16, ViewConditions},
    error::Error,
    illuminant::{CieIlluminant, Planck},
    lab::CieLab,
    observer::observer_data::ObserverData,
    rgb::RgbSpace,
    spectrum::{to_wavelength, Spectrum, NS, SPECTRUM_WAVELENGTH_RANGE},
    traits::{Filter, Light},
    xyz::{RelXYZ, XYZ},
};
use nalgebra::{Matrix3, Vector3};
use std::{fmt, ops::RangeInclusive};

use strum::{AsRefStr, EnumCount, EnumIter};

/// Selects a CIE standard colorimetric observer.
///
/// The tag is embedded in every [`XYZ`] and [`Rgb`](crate::rgb::Rgb) value so that
/// operations across incompatible observers can be detected at runtime.  Each variant
/// is a lightweight index; the color-matching function tables are stored in
/// [`observer_data`].
#[cfg_attr(target_arch = "wasm32", wasm_bindgen::prelude::wasm_bindgen)]
#[derive(Clone, Copy, Default, PartialEq, Eq, Hash, Debug, EnumIter, EnumCount, AsRefStr)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[non_exhaustive]
pub enum Observer {
    /// CIE 1931 2° standard observer — the default for most colorimetry.
    ///
    /// Used by sRGB, ICC profiles, and CIE CRI Ra.
    #[default]
    Cie1931,
    /// CIE 1964 10° supplementary standard observer.
    ///
    /// Preferred when the viewed area subtends more than ~4° at the eye.
    /// Used by CIE 224:2017 / ANSI/IES TM-30 for colour fidelity calculations.
    Cie1964,
    /// CIE 2015 2° observer — CMFs constructed as linear transforms of the Stockman & Sharpe
    /// (2000) cone fundamentals.
    ///
    /// More accurate than `Cie1931` in the short-wavelength (blue) region.
    Cie2015,
    /// CIE 2015 10° observer — CMFs constructed as linear transforms of the Stockman & Sharpe
    /// (2000) cone fundamentals.
    ///
    /// Wide-field counterpart of [`Cie2015`](Observer::Cie2015).
    Cie2015_10,
}

impl Observer {
    /// Get a reference to the data for the specified `Observer`.
    fn data(&self) -> &'static ObserverData {
        match self {
            Observer::Cie1931 => &observer_data::CIE1931,
            Observer::Cie1964 => &observer_data::CIE1964,
            Observer::Cie2015 => &observer_data::CIE2015,
            Observer::Cie2015_10 => &observer_data::CIE2015_10,
        }
    }

    /// Returns the name of the observer.
    pub fn name(&self) -> &'static str {
        self.data().name
    }

    /// Returns the spectral locus for this observer.
    ///
    /// The spectral locus is the boundary of the area of all physical colors in a chromaticity diagram.
    pub fn spectral_locus(&self) -> &SpectralLocus {
        self.data()
            .spectral_locus
            .get_or_init(|| SpectralLocus::new(*self))
    }

    /// Computes relative XYZ tristimulus values for a filtered light source.
    ///
    /// The unfiltered light spectrum sets the white point (normalised to Y = 100);
    /// the product of the light and filter spectra is integrated to give the stimulus XYZ.
    /// The result is a [`RelXYZ`] that carries both the stimulus and the white point,
    /// as required by CIELAB and CIECAM calculations.
    pub fn rel_xyz(&self, light: &dyn Light, filter: &dyn Filter) -> RelXYZ {
        let xyzn = self.xyz_from_spectrum(&light.spectrum());
        let s = *light.spectrum() * *filter.spectrum();
        let xyz = self.xyz_from_spectrum(&s);
        let scale = 100.0 / xyzn.xyz.y;
        // unwrap OK as we are using only one observer (self) here
        RelXYZ::from_xyz(xyz * scale, xyzn * scale).unwrap()
    }

    /// Calulates Tristimulus values for an object implementing the [Light] trait, and an optional [Filter],
    /// filtering the light.
    ///
    /// The Light trait is implemented by [`CieIlluminant`] and [Illuminant](crate::illuminant::Illuminant).
    ///
    /// [`Colorant`](crate::colorant::Colorant) implments the [`Filter`] trait.
    /// [`Rgb`](crate::rgb::Rgb), which represents a display pixel, implements both in this library.
    /// As a light, it is the light emitted from the pixel, as a filter it is the RGB-composite
    /// filter which is applied to the underlying standard illuminant of color space.
    pub fn xyz(&self, light: &dyn Light, filter: Option<&dyn Filter>) -> XYZ {
        let xyzn = light.xyzn(self.data().tag, None);
        if let Some(flt) = filter {
            let s = *light.spectrum() * *flt.spectrum();
            let xyz = self.xyz_from_spectrum(&s);
            let scale = 100.0 / xyzn.xyz.y;
            xyz * scale
        } else {
            xyzn
        }
    }

    /// Computes CIELAB L\*a\*b\* values for a light source and filter combination,
    /// normalised so the unfiltered source has Y = 100.
    ///
    /// # Arguments
    /// * `light` - A light source implementing [`Light`], such as [`CieIlluminant`].
    /// * `filter` - A filter implementing [`Filter`], such as `Colorant`.
    ///
    /// # Returns
    /// [`CieLab`] — the CIELAB representation of the filtered stimulus relative to the white point.
    pub fn lab(&self, light: &dyn Light, filter: &dyn Filter) -> CieLab {
        let rxyz = self.rel_xyz(light, filter);
        // unwrap OK as we are using only one observer (self) here
        CieLab::from_rxyz(rxyz)
    }

    /// Calculates the L*a*b* CIELAB D65 values of a Colorant, using D65 as an illuminant.
    /// Convenience method for `lab` with D65 illuminant, and using an illuminance value of 100.0.
    pub fn lab_d65(&self, filter: &dyn Filter) -> CieLab {
        self.lab(&crate::illuminant::D65, filter)
    }

    /// Calculates the L*a*b* CIELAB D50 values of a Colorant, using D50 as an illuminant.
    /// Convenience method for `lab` with D50 illuminant, and using an illuminance value of 100.0.
    pub fn lab_d50(&self, filter: &dyn Filter) -> CieLab {
        self.lab(&crate::illuminant::D50, filter)
    }

    /// Computes CIECAM16 color appearance model values for a light source and filter combination,
    /// normalised so the unfiltered source has Y = 100.
    ///
    /// # Arguments
    /// * `light` - A light source implementing [`Light`], such as [`CieIlluminant`].
    /// * `filter` - A filter implementing [`Filter`], such as `Colorant`.
    /// * `vc` - Viewing conditions for the CIECAM16 model.
    ///
    /// # Returns
    /// [`CieCam16`] — CIECAM16 correlates for the filtered stimulus under the given viewing conditions.
    pub fn ciecam16(&self, light: &dyn Light, filter: &dyn Filter, vc: ViewConditions) -> CieCam16 {
        let rxyz = self.rel_xyz(light, filter);
        // unwrap OK as we are using only one observer (self) here
        CieCam16::from_xyz(rxyz, vc)
    }

    /// Computes CIECAM02 color appearance model values for a light source and filter combination,
    /// normalised so the unfiltered source has Y = 100.
    ///
    /// # Arguments
    /// * `light` - A light source implementing [`Light`], such as [`CieIlluminant`].
    /// * `filter` - A filter implementing [`Filter`], such as `Colorant`.
    /// * `vc` - Viewing conditions for the CIECAM02 model.
    ///
    /// # Returns
    /// [`CieCam02`] — CIECAM02 correlates for the filtered stimulus under the given viewing conditions.
    pub fn ciecam02(&self, light: &dyn Light, filter: &dyn Filter, vc: ViewConditions) -> CieCam02 {
        let rxyz = self.rel_xyz(light, filter);
        // unwrap OK as we are using only one observer (self) here
        CieCam02::from_xyz(rxyz, vc)
    }

    /// Calculates tristimulus values for a general spectrum, returned as an [`XYZ`] object.
    ///
    /// The spectrum is interpreted as an illuminant SPD.  This method produces raw (unnormalised)
    /// XYZ data; the Y value is not scaled to 100.
    pub fn xyz_from_spectrum(&self, spectrum: &Spectrum) -> XYZ {
        let xyz = self.data().data * spectrum.0 * self.data().lumconst;
        XYZ::from_vec(xyz, self.data().tag)
    }

    /// Calculates the luminous value (Y tristimulus value) for a general spectrum.
    pub fn y_from_spectrum(&self, spectrum: &Spectrum) -> f64 {
        (self.data().data.row(1) * spectrum.0 * self.data().lumconst)[(0, 0)]
    }

    /// Returns the observer's color matching function (CMF) data as an [XYZ] tristimulus
    /// value for the given wavelength.
    ///
    /// This allows access to the underlying data for the entire range of wavelengths, 380-780nm.
    /// However, please read the documentation for the
    /// [`spectral_locus_wavelength_range`](Self::spectral_locus_wavelength_range) method on
    /// situations where you might not want to sample the full range.
    pub fn xyz_at_wavelength(&self, wavelength: usize) -> Result<XYZ, Error> {
        if !SPECTRUM_WAVELENGTH_RANGE.contains(&wavelength) {
            return Err(Error::WavelengthOutOfRange);
        };
        let &[x, y, z] = self
            .data()
            .data
            .column(wavelength - SPECTRUM_WAVELENGTH_RANGE.start())
            .as_ref();
        Ok(XYZ::from_vec(Vector3::new(x, y, z), self.data().tag))
    }

    /// Calculates the relative XYZ tristimulus values of monochromatic stimuli.
    ///
    /// Monochromatic stimuli are pure spectral colors, like those seen when white light
    /// is dispersed by a prism. This method computes their XYZ values under a given illuminant.
    ///
    /// # Details
    /// - Computes XYZ values for wavelengths from 380nm to 780nm in 1nm steps
    /// - Multiplies observer color matching functions by illuminant spectrum
    /// - Normalizes results relative to the illuminant's white point
    /// - Scales output to 100 lux illuminance
    ///
    /// # Difference from Spectral Locus
    /// While the spectral locus shows only chromaticity coordinates (x,y) of pure spectral
    /// colors, this method provides full XYZ values including luminance information.
    ///
    /// # Implementation Notes
    /// - Each wavelength is treated as a monochromatic stimulus (delta function)
    /// - Results are typically low in magnitude due to the narrow bandwidth
    /// - Values are normalized relative to the illuminant's total energy
    ///
    /// # Arguments
    /// - `ref_white`: Reference illuminant (e.g., D65, D50) for normalisation.
    ///
    /// # Returns
    /// A vector of `(wavelength_nm, RelXYZ)` pairs (380–780 nm), where each `RelXYZ`
    /// is scaled to 100 lux illuminance.
    pub fn monochromes(&self, ref_white: CieIlluminant) -> Vec<(usize, RelXYZ)> {
        let mut obs = self.data().data;
        let white = &ref_white.illuminant().as_ref().0;
        for r in 0..3 {
            for c in 0..NS {
                obs[(r, c)] *= white[c];
            }
        }
        let xyzn_vec = obs.column_sum();
        let xyzn = XYZ::from_vec(xyzn_vec, self.data().tag);
        let mut v = Vec::with_capacity(NS);
        for w in SPECTRUM_WAVELENGTH_RANGE {
            let xyz = obs.column(w - SPECTRUM_WAVELENGTH_RANGE.start()).into();
            let rxyz = RelXYZ::from_vec(xyz, xyzn);
            v.push((w, rxyz.set_illuminance(100.0)));
        }
        v
    }

    /// Returns relative XYZ values for the spectral locus, restricted to the unique wavelength range.
    ///
    /// This is a convenience wrapper around [`monochromes`](Self::monochromes) that discards
    /// the endpoints where the locus folds back on itself — see
    /// [`spectral_locus_wavelength_range`](Self::spectral_locus_wavelength_range).
    pub fn trimmed_spectral_locus(&self, ref_white: CieIlluminant) -> Vec<(usize, RelXYZ)> {
        let sl_full = self.monochromes(ref_white);
        let valid_range = self.spectral_locus_wavelength_range();
        sl_full
            .into_iter()
            .filter(|(w, _)| valid_range.contains(w))
            .collect()
    }

    /// Tristimulus values for the Standard Illuminants in this library.
    ///
    /// Values are not normalized by default, unless an illuminance value is provided.
    ///
    /// Results are not cached. For the common cases of D65 and D50, use the dedicated
    /// [`Observer::xyz_d65`] and [`Observer::xyz_d50`] methods, which cache their results.
    pub fn xyz_cie_table(&self, std_illuminant: &CieIlluminant, illuminance: Option<f64>) -> XYZ {
        let xyz = self.xyz_from_spectrum(std_illuminant.illuminant().as_ref());
        if let Some(l) = illuminance {
            xyz.set_illuminance(l)
        } else {
            xyz
        }
    }

    /// XYZ tristimulus values for the CIE standard daylight illuminant D65 (buffered).
    pub fn xyz_d65(&self) -> XYZ {
        *self.data().d65.get_or_init(|| {
            self.xyz_from_spectrum(CieIlluminant::D65.illuminant().as_ref())
                .set_illuminance(100.0)
        })
    }

    /// XYZ tristimulus values for the CIE standard daylight illuminant D50 (buffered).
    pub fn xyz_d50(&self) -> XYZ {
        *self.data().d50.get_or_init(|| {
            self.xyz_from_spectrum(CieIlluminant::D50.illuminant().as_ref())
                .set_illuminance(100.0)
        })
    }

    /// Computes XYZ tristimulus values for a Planckian (blackbody) radiator at the given CCT.
    ///
    /// The spectral radiance is evaluated from Planck's law and integrated against
    /// the observer's color-matching functions.  Combine with
    /// [`xyz_planckian_locus_slope`](Self::xyz_planckian_locus_slope) to follow
    /// the Planckian locus in XYZ space.
    pub fn xyz_planckian_locus(&self, cct: f64) -> XYZ {
        let p = Planck::new(cct);
        self.xyz_from_fn(|l| p.at_wavelength(to_wavelength(l, 0.0, 1.0)))
    }

    /// Calculates the chromaticity coordinates (x, y) of the Planckian locus
    /// over a range from 1000K to 100_000K.
    pub fn planckian_locus(&self) -> Vec<(f64, f64)> {
        let mut v = Vec::with_capacity(100);
        for i in 1..100 {
            let cct = 1E6 / (i as f64 * 10.0);
            let xy = self.xyz_planckian_locus(cct).chromaticity();
            v.push((xy.x(), xy.y()));
        }
        v
    }

    /// The slope of the Plancking locus as a (dX/dT, dY/dT, dZ/dT) contained in
    /// a XYZ object.
    pub fn xyz_planckian_locus_slope(&self, cct: f64) -> XYZ {
        let p = Planck::new(cct);
        self.xyz_from_fn(|l| p.slope_at_wavelength(to_wavelength(l, 0.0, 1.0)))
    }

    /// Computes the 3×3 RGB-to-XYZ matrix for a color space with an optional alternate white point.
    ///
    /// When `opt_white` is `None` the color space's own white point is used, producing the
    /// same matrix as [`rgb2xyz_matrix`](Self::rgb2xyz_matrix).  Supply a D50 [`XYZ`] value
    /// to obtain a matrix adapted to D50 — this is required for ICC profile `rXYZ`/`gXYZ`/`bXYZ`
    /// tags when the color space is natively D65 (e.g. sRGB).
    ///
    /// # Notes
    ///
    /// * XYZ is used for the white point (rather than a spectrum) to match the ICC profile
    ///   standard, which specifies primaries as XYZ coordinates.
    ///
    pub fn calc_rgb2xyz_matrix_with_alt_white(
        &self,
        rgbspace: RgbSpace,
        opt_white: Option<XYZ>,
    ) -> Matrix3<f64> {
        let mut rgb2xyz: nalgebra::Matrix<
            f64,
            nalgebra::Const<3>,
            nalgebra::Const<3>,
            nalgebra::ArrayStorage<f64, 3, 3>,
        > = Matrix3::from_iterator(
            rgbspace
                .primaries()
                .iter()
                .flat_map(|s| self.xyz_from_spectrum(s).set_illuminance(1.0).to_array()),
        );
        let xyzw = opt_white
            .unwrap_or(self.xyz(&rgbspace.white(), None))
            .set_illuminance(1.0);
        // check if the observers match
        assert_eq!(&xyzw.observer(), self);
        let decomp = rgb2xyz.lu();
        // unwrap: only used with library color spaces
        let rgbw = decomp.solve(&xyzw.xyz).unwrap();
        for (i, mut col) in rgb2xyz.column_iter_mut().enumerate() {
            col *= rgbw[i];
        }
        rgb2xyz
    }

    /// Computes the 3×3 RGB-to-XYZ matrix for a color space from first principles.
    ///
    /// Pre-computed, cached matrices for all built-in observer / color space combinations are
    /// returned by [`rgb2xyz_matrix`](Self::rgb2xyz_matrix).  Use this method only when
    /// adding a new color space or regenerating built-in tables via `cargo xtask gen`.
    pub fn calc_rgb2xyz_matrix(&self, rgbspace: RgbSpace) -> Matrix3<f64> {
        //   let space = rgbspace.data();
        let mut rgb2xyz: nalgebra::Matrix<
            f64,
            nalgebra::Const<3>,
            nalgebra::Const<3>,
            nalgebra::ArrayStorage<f64, 3, 3>,
        > = Matrix3::from_iterator(
            rgbspace
                .primaries()
                .iter()
                .flat_map(|s| self.xyz_from_spectrum(s).set_illuminance(1.0).to_array()),
        );
        let xyzw = self.xyz(&rgbspace.white(), None).set_illuminance(1.0);
        let decomp = rgb2xyz.lu();
        // unwrap: only used with library color spaces
        let rgbw = decomp.solve(&xyzw.xyz).unwrap();
        for (i, mut col) in rgb2xyz.column_iter_mut().enumerate() {
            col *= rgbw[i];
        }
        rgb2xyz
    }

    /// Computes the 3×3 XYZ-to-RGB matrix for a color space (the inverse of
    /// [`calc_rgb2xyz_matrix`](Self::calc_rgb2xyz_matrix)).
    ///
    /// Pre-computed, cached matrices for all built-in combinations are returned by
    /// [`xyz2rgb_matrix`](Self::xyz2rgb_matrix).  Use this method only when regenerating
    /// built-in tables via `cargo xtask gen`.
    pub fn calc_xyz2rgb_matrix(&self, rgbspace: RgbSpace) -> Matrix3<f64> {
        // unwrap: only used with library color spaces
        self.calc_rgb2xyz_matrix(rgbspace).try_inverse().unwrap()
    }

    //  pub fn rgb(&self, rgbspace: RgbSpace) -> Matrix3<f64> {
    //      self.rgb2xyz(rgbspace)
    //  }

    /// Returns the wavelength range (in nanometres) over which the spectral locus is monotone.
    ///
    /// The spectral locus (the horseshoe boundary of all physical colors in a chromaticity
    /// diagram) folds back on itself near the blue and red extremes: consecutive wavelengths
    /// can map to the same or a reversed chromaticity coordinate.  This causes problems for
    /// dominant-wavelength calculations and locus plots.
    ///
    /// This method returns the sub-range for which each wavelength maps to a unique
    /// chromaticity point distinct from its neighbours.
    /// See Wikipedia's [CIE 1931 Color Space](https://en.wikipedia.org/wiki/CIE_1931_color_space).
    ///
    /// To get the tristimulus values of the spectral locus, use
    /// [`xyz_at_wavelength`](Self::xyz_at_wavelength).
    pub fn spectral_locus_wavelength_range(&self) -> RangeInclusive<usize> {
        self.data().spectral_locus_range.clone()
    }

    /// Calculates the XYZ tristimulus values for a spectrum defined by a function.
    ///
    /// The input function `f` should accept a floating-point value in the range `[0.0, 1.0]`,
    /// where `0.0` corresponds to a wavelength of 380 nm and `1.0` to 780 nm.
    /// The function will be called once for each wavelength step (401 times at 1 nm intervals).
    ///
    /// # Arguments
    /// * `f` - A function that takes a floating-point value in the range `[0.0, 1.0]` and returns
    ///   the spectral value at that wavelength, in units of watts per square meter per nanometer (W/m²/nm) for
    ///   illuminants, or Watts per square meter per steradian per nanometer (W/m²/sr/nm) for stimuli.
    ///
    /// # Notes
    /// - This method is used in the library to compute the Planckian locus (the color of blackbody
    ///   radiators), as described by Planck's law.
    /// - For colorants, use [`xyz_from_colorant_fn`](Self::xyz_from_colorant_fn).
    pub fn xyz_from_fn(&self, f: impl Fn(f64) -> f64) -> XYZ {
        let xyz = self
            .data()
            .data
            .column_iter()
            .enumerate()
            .fold(Vector3::zeros(), |acc, (i, cmf)| {
                acc + cmf * f(i as f64 / (NS - 1) as f64)
            });
        XYZ {
            xyz,
            observer: self.data().tag,
        }
    }

    /// Calculates XYZ tristimulus values for an analytic representation of a spectral distribution of
    /// a filter or a color patch, using a normalized wavelength domain ranging from a value of 0.0 to 1.0,
    /// illuminated with a standard illuminant.
    ///
    /// The spectral values should be defined within a range from 0.0 to 1.0, and are clamped otherwise.
    /// The resulting XYZ value will have relative Y values in the range from 0 to 100.0.
    ///
    /// # Examples
    /// A linear high pass filter, with a value of 0.0 for a wavelength of 380nm, and a value of 1.0 for 780nm,
    /// and converting the resulting value to RGB values.
    /// ```
    /// use colorimetry::{rgb::RgbSpace::SRGB, observer::Observer::Cie1931, illuminant::CieIlluminant};
    /// let rgb: [u8;3] = Cie1931.xyz_from_colorant_fn(&CieIlluminant::D65, |x|x).rgb(SRGB).clamp().into();
    /// assert_eq!(rgb, [212, 171, 109]);
    /// ```
    /// Linear low pass filter, with a value of 1.0 for a wavelength of 380nm, and a value of 0.0 for 780nm,
    /// and converting the resulting value to RGB values.
    /// ```
    /// use colorimetry::{rgb::RgbSpace::SRGB, observer::Observer::Cie1931, illuminant::CieIlluminant};
    /// let rgb: [u8;3] = Cie1931.xyz_from_colorant_fn(&CieIlluminant::D65, |x|1.0-x).rgb(SRGB).clamp().into();
    /// assert_eq!(rgb, [158, 202, 237]);
    /// ```
    pub fn xyz_from_colorant_fn(&self, illuminant: &CieIlluminant, f: impl Fn(f64) -> f64) -> XYZ {
        let ill = illuminant.illuminant();
        let xyzn = self.xyz_cie_table(illuminant, None);
        let xyz = self
            .data()
            .data
            .column_iter()
            .zip(ill.0 .0.iter())
            .enumerate()
            .fold(Vector3::zeros(), |acc, (i, (cmfi, sv))| {
                acc + cmfi * f(i as f64 / (NS - 1) as f64).clamp(0.0, 1.0) * *sv
            });
        let scale = 100.0 / xyzn.xyz.y;
        XYZ {
            xyz: xyz * self.data().lumconst * scale,
            observer: self.data().tag,
        }
    }
}

impl fmt::Display for Observer {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        self.name().fmt(f)
    }
}

#[cfg(test)]
mod obs_test {

    use super::Observer::Cie1964;
    use super::{Observer, Observer::Cie1931};
    use crate::{
        illuminant::CieIlluminant, rgb::RgbSpace, spectrum::SPECTRUM_WAVELENGTH_RANGE, xyz::XYZ,
    };
    use approx::{assert_abs_diff_eq, assert_ulps_eq};
    use strum::IntoEnumIterator as _;

    #[test]
    fn test_rel_xyz() {
        use crate::colorant::Colorant;

        let obs = Cie1931;
        let light = CieIlluminant::D65;
        let filter = Colorant::gray(0.5);

        let rxyz = obs.rel_xyz(&light, &filter);
        assert_ulps_eq!(rxyz.xyz() * 2.0, rxyz.white_point(), epsilon = 1E-5);
        assert_ulps_eq!(rxyz.white_point(), obs.xyz_d65(), epsilon = 1E-5);
    }

    #[test]
    fn test_cie1931_spectral_locus_min_max() {
        let wavelength_lange = Cie1931.spectral_locus_wavelength_range();
        let chromaticity = Cie1931
            .xyz_at_wavelength(*wavelength_lange.start())
            .unwrap()
            .chromaticity();
        assert_ulps_eq!(chromaticity.x(), 0.17411, epsilon = 1E-5);
        assert_ulps_eq!(chromaticity.y(), 0.00496, epsilon = 1E-5);

        let chromaticity = Cie1931
            .xyz_at_wavelength(*wavelength_lange.end())
            .unwrap()
            .chromaticity();
        assert_ulps_eq!(chromaticity.x(), 0.73469, epsilon = 1E-5);
        assert_ulps_eq!(chromaticity.y(), 0.26531, epsilon = 1E-5);
    }

    /// Ensure all observers have sane spectral locus values
    #[test]
    fn test_spectral_locus_full() {
        for observer in Observer::iter() {
            let wavelength_range = observer.spectral_locus_wavelength_range();

            // Basic sanity checking of the spectral locus wavelength range values
            assert!(wavelength_range.start() >= SPECTRUM_WAVELENGTH_RANGE.start());
            assert!(wavelength_range.end() <= SPECTRUM_WAVELENGTH_RANGE.end());
            assert!(wavelength_range.start() < wavelength_range.end());

            // Check that xyz_at_wavelength returns sane values in the allowed range
            for wavelength in wavelength_range {
                let xyz = observer.xyz_at_wavelength(wavelength).unwrap();
                let chromaticity = xyz.chromaticity();
                assert!((0.0..=1.0).contains(&chromaticity.x()));
                assert!((0.0..=1.0).contains(&chromaticity.y()));

                // Check that all chromaticity coordinates on the spectral locus can be converted
                // to XYZ.
                let _xyz2 = XYZ::from_chromaticity(chromaticity, None, Some(observer)).unwrap();
            }
        }
    }

    /// Ensure XYZ values around all supported observer spectral locuses' can be converted to RGB
    #[test]
    fn test_spectral_locus_to_rgb() {
        for observer in Observer::iter() {
            eprintln!("Testing observer {observer:?}");
            for wavelength in observer.spectral_locus_wavelength_range() {
                let xyz = observer.xyz_at_wavelength(wavelength).unwrap();
                for rgbspace in RgbSpace::iter() {
                    let _rgb = xyz.rgb(rgbspace);
                }
            }
        }
    }

    #[test]
    // Data from http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html Lindbloom's values
    // are reproduced with an accuracy of 3E-4, which is a small, but significant difference.  This
    // difference is caused by a difference in the display's white point, due to wavelength domain
    // differences.  Here we use a domain from 380 to 780 with a step of 1 nanometer for the
    // spectra, and in specific for the color matching functions, as recommended by CIE15:2004, and
    // the color matching functions provided by the CIE in their online dataset section. The spectra
    // for the primaries are chosen to match the RGB primary values as given by the Color Space
    // specifications.  The white point uses the standard's spectral distribution, as provided by
    // the CIE, clipped to a domain from 380 to 780 nanometer.  See `colorimetry::data::D65`.
    fn test_rgb2xyz_cie1931() {
        let want = nalgebra::Matrix3::new(
            0.4124564, 0.3575761, 0.1804375, 0.2126729, 0.7151522, 0.0721750, 0.0193339, 0.1191920,
            0.9503041,
        );
        let got = Cie1931.rgb2xyz_matrix(crate::rgb::RgbSpace::SRGB);
        approx::assert_ulps_eq!(want, got, epsilon = 3E-4);
    }

    #[test]
    // Check the inverse transformation.
    // See comments at `test_rgb2xyz_cie1931`.
    fn test_xyz2rgb_cie1931() {
        let want = nalgebra::Matrix3::new(
            3.2404542, -1.5371385, -0.4985314, -0.9692660, 1.8760108, 0.0415560, 0.0556434,
            -0.2040259, 1.0572252,
        );
        let got = Cie1931.xyz2rgb_matrix(crate::rgb::RgbSpace::SRGB);
        approx::assert_ulps_eq!(want, got, epsilon = 3E-4);
    }

    #[test]
    fn test_planckian_locus() {
        // see https://www.waveformlighting.com/tech/calculate-cie-1931-xy-coordinates-from-cct
        // for test data (not clear what CMF domain they use)
        let xy = Cie1931
            .xyz_planckian_locus(3000.0)
            .chromaticity()
            .to_array();
        approx::assert_abs_diff_eq!(&xy.as_ref(), &[0.43693, 0.40407].as_ref(), epsilon = 2E-5);

        let xy = Cie1931
            .xyz_planckian_locus(6500.0)
            .chromaticity()
            .to_array();
        approx::assert_abs_diff_eq!(&xy.as_ref(), &[0.31352, 0.32363].as_ref(), epsilon = 6E-5);
    }

    #[test]
    fn test_xyz_from_illuminant_x_fn() {
        let xyz = Cie1931.xyz_from_colorant_fn(&CieIlluminant::D65, |_v| 1.0);
        let d65xyz = Cie1931.xyz_d65().xyz;
        approx::assert_ulps_eq!(
            xyz,
            crate::xyz::XYZ::from_vec(d65xyz, crate::observer::Observer::Cie1931)
        );
    }

    #[test]
    fn test_xyz_of_sample_with_standard_illuminant() {
        use crate::illuminant::CieIlluminant::D65 as d65;
        let xyz = Cie1931
            .xyz(&d65, Some(&crate::colorant::Colorant::white()))
            .set_illuminance(100.0);
        approx::assert_ulps_eq!(xyz, Cie1931.xyz_from_colorant_fn(&d65, |_| 1.0));

        let xyz = Cie1931.xyz(&d65, Some(&crate::colorant::Colorant::black()));
        approx::assert_ulps_eq!(xyz, Cie1931.xyz_from_colorant_fn(&d65, |_| 0.0));
    }

    #[test]
    fn test_xyz_d65_d50() {
        let cie1931_d65_xyz = Cie1931.xyz_d65();
        approx::assert_ulps_eq!(
            cie1931_d65_xyz.to_array().as_ref(),
            [95.047, 100.0, 108.883].as_ref(),
            epsilon = 5E-2
        );

        let cie1931_d50_xyz = Cie1931.xyz_d50();
        approx::assert_ulps_eq!(
            cie1931_d50_xyz.to_array().as_ref(),
            [96.421, 100.0, 82.519].as_ref(),
            epsilon = 5E-2
        );
    }

    #[test]
    fn test_xyz_d65_d50_cie1964() {
        let cie1964_d50_xyz = Cie1964.xyz_d50();
        approx::assert_ulps_eq!(
            cie1964_d50_xyz.to_array().as_ref(),
            [96.720, 100.0, 81.427].as_ref(),
            epsilon = 5E-2
        );

        let cie1964_d65_xyz = Cie1964.xyz_d65();
        approx::assert_ulps_eq!(
            cie1964_d65_xyz.to_array().as_ref(),
            [94.811, 100.0, 107.304].as_ref(),
            epsilon = 5E-2
        );
    }

    #[test]
    fn test_spectral_locus() {
        use crate::illuminant::CieIlluminant;
        use crate::observer::Observer::Cie1931;

        let sl = Cie1931.monochromes(CieIlluminant::D65);
        // Check the first and last points of the spectral locus
        // Data obtained from spreadsheet using data directly downloaded from cie.co.at
        assert_eq!(sl.len(), 401);
        let xyz_first = sl[0].1.xyz().to_array();
        assert_abs_diff_eq!(
            xyz_first.as_ref(),
            [6.46976E-04, 1.84445E-05, 3.05044E-03].as_ref(),
            epsilon = 1E-5
        );
        let xyz_last = sl[sl.len() - 1].1.xyz().to_array();
        assert_abs_diff_eq!(
            xyz_last.as_ref(),
            [2.48982E-05, 8.99121E-06, 0.0].as_ref(),
            epsilon = 1E-5
        );
        // 550 nm
        let xyz_550 = sl[170].1.xyz().to_array();
        assert_abs_diff_eq!(
            xyz_550.as_ref(),
            [0.4268018, 0.9796899, 0.0086158].as_ref(),
            epsilon = 1E-4
        );
    }

    #[test]
    fn test_as_ref_str() {
        // Ensure that the observer names are correct
        assert_eq!(Cie1931.as_ref(), "Cie1931");
        assert_eq!(Cie1964.as_ref(), "Cie1964");
        assert_eq!(Observer::Cie2015.as_ref(), "Cie2015");
        assert_eq!(Observer::Cie2015_10.as_ref(), "Cie2015_10");
    }

    #[test]
    // Test white point, should be
    fn test_rgb_xyz_white() {
        for obs in Observer::iter() {
            for space in RgbSpace::iter() {
                let d65 = obs.xyz_d65().set_illuminance(1.0);
                let xyz2rgb = obs.xyz2rgb_matrix(space);
                let rgb = xyz2rgb * d65.xyz;
                rgb.iter().for_each(|&v| {
                    assert_abs_diff_eq!(v, 1.0, epsilon = 2E-6);
                });

                let xyz_round_trip = XYZ::from_vec(obs.rgb2xyz_matrix(space) * rgb, obs);
                assert_abs_diff_eq!(d65, xyz_round_trip, epsilon = 2E-6);
            }
        }
    }

    /*
    #[test]
    // Test CIE 1931 primaries
    fn test_rgb_primaries() {
        use nalgebra::Vector3;
        let obs = Observer::Cie1931;
        for space in RgbSpace::iter() {
            for (i, &primary) in [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
                .iter()
                .enumerate()
            {
                let want = space.primaries_chromaticity()[i].to_array();
                let prim_vec = Vector3::from(primary);
                let xy = XYZ::from_vec(obs.rgb2xyz_matrix(space) * prim_vec, obs).chromaticity();
                println!("{obs}, {space} {want:?}");
                assert_abs_diff_eq!(xy.x(), want[0], epsilon = 1E-5);
                assert_abs_diff_eq!(xy.y(), want[1], epsilon = 1E-5);
            }
        }
    }
     */
}