devela/media/color/spectral/

kubelka_munk.rs

// devela::media::color::spectral::kubelka-munk
//
//!
//!
//! Vendored from: <https://github.com/rvanwijnen/spectral.js> v2.0.2
//!
//

use devela::{ExtFloat, Float};

#[cfg(feature = "alloc")]
use crate::{format, String, ToString, Vec};

/* Format constants */

pub const RGB: u8 = 0;
pub const RGBA: u8 = 1;
pub const HEX: u8 = 2;
pub const HEXA: u8 = 3;

/* Configuration constants */

pub const SIZE: usize = 38;
pub const GAMMA: f64 = 2.4;
pub const EPSILON: f64 = 1e-8;

/* Refractive index constants */

// K1 and K2 are used for the Saunderson correction,
// which is used for how the color looks behind a refractive surface like glass.

/// Refractive Index
///
/// - <https://en.wikipedia.org/wiki/Refractive_index>
/// - <https://en.wikipedia.org/wiki/List_of_refractive_indices>
///
/// - Vacuum = 1.0 (no surface)
/// - Glass  = 1.5
pub const RI: f64 = 1.0;

// MAYBE
// /// K1: Constant used for Saunderson correction with refractive index.
// pub const K1: f64 = Float(RI - 1.0).const_powi(2).0 / Float(RI + 1.0).const_powi(2).0;

/// K2: Modifier for lightening or darkening color appearance.
/// (0 = neutral, - = lighten, + = darken)
pub const K2: f64 = 0.0;

/* Spectral data arrays */

/// Spectral Power Distribution constants for cyan.
#[allow(clippy::unreadable_literal)]
pub const SPD_C: [f64; SIZE] = [
    0.96853629, 0.96855103, 0.96859338, 0.96877345, 0.96942204, 0.97143709, 0.97541862, 0.98074186,
    0.98580992, 0.98971194, 0.99238027, 0.99409844, 0.995172, 0.99576545, 0.99593552, 0.99564041,
    0.99464769, 0.99229579, 0.98638762, 0.96829712, 0.89228016, 0.53740239, 0.15360445, 0.05705719,
    0.03126539, 0.02205445, 0.01802271, 0.0161346, 0.01520947, 0.01475977, 0.01454263, 0.01444459,
    0.01439897, 0.0143762, 0.01436343, 0.01435687, 0.0143537, 0.01435408,
];

/// Spectral Power Distribution constants for magenta.
#[allow(clippy::unreadable_literal)]
pub const SPD_M: [f64; SIZE] = [
    0.51567122, 0.5401552, 0.62645502, 0.75595012, 0.92826996, 0.97223624, 0.98616174, 0.98955255,
    0.98676237, 0.97312575, 0.91944277, 0.32564851, 0.13820628, 0.05015143, 0.02912336, 0.02421691,
    0.02660696, 0.03407586, 0.04835936, 0.0001172, 0.00008554, 0.85267882, 0.93188793, 0.94810268,
    0.94200977, 0.91478045, 0.87065445, 0.78827548, 0.65738359, 0.59909403, 0.56817268, 0.54031997,
    0.52110241, 0.51041094, 0.50526577, 0.5025508, 0.50126452, 0.50083021,
];

/// Spectral Power Distribution constants for yellow.
#[allow(clippy::unreadable_literal)]
pub const SPD_Y: [f64; SIZE] = [
    0.02055257, 0.02059936, 0.02062723, 0.02073387, 0.02114202, 0.02233154, 0.02556857, 0.03330189,
    0.05185294, 0.10087639, 0.24000413, 0.53589066, 0.79874659, 0.91186529, 0.95399623, 0.97137099,
    0.97939505, 0.98345207, 0.98553736, 0.98648905, 0.98674535, 0.98657555, 0.98611877, 0.98559942,
    0.98507063, 0.98460039, 0.98425301, 0.98403909, 0.98388535, 0.98376116, 0.98368246, 0.98365023,
    0.98361309, 0.98357259, 0.98353856, 0.98351247, 0.98350101, 0.98350852,
];

/// Spectral Power Distribution constants for red.
#[allow(clippy::unreadable_literal)]
pub const SPD_R: [f64; SIZE] = [
    0.03147571, 0.03146636, 0.03140624, 0.03119611, 0.03053888, 0.02856855, 0.02459485, 0.0192952,
    0.01423112, 0.01033111, 0.00765876, 0.00593693, 0.00485616, 0.00426186, 0.00409039, 0.00438375,
    0.00537525, 0.00772962, 0.0136612, 0.03181352, 0.10791525, 0.46249516, 0.84604333, 0.94275572,
    0.96860996, 0.97783966, 0.98187757, 0.98377315, 0.98470202, 0.98515481, 0.98537114, 0.98546685,
    0.98550011, 0.98551031, 0.98550741, 0.98551323, 0.98551563, 0.98551547,
];

/// Spectral Power Distribution constants for green.
#[allow(clippy::unreadable_literal)]
pub const SPD_G: [f64; SIZE] = [
    0.49108579, 0.46944057, 0.4016578, 0.2449042, 0.0682688, 0.02732883, 0.013606, 0.01000187,
    0.01284127, 0.02636635, 0.07058713, 0.70421692, 0.85473994, 0.95081565, 0.9717037, 0.97651888,
    0.97429245, 0.97012917, 0.9425863, 0.99989207, 0.99989891, 0.13823139, 0.06968113, 0.05628787,
    0.06111561, 0.08987709, 0.13656016, 0.22169624, 0.32176956, 0.36157329, 0.4836192, 0.46488579,
    0.47440306, 0.4857699, 0.49267971, 0.49625685, 0.49807754, 0.49889859,
];

/// Spectral Power Distribution constants for blue.
#[allow(clippy::unreadable_literal)]
pub const SPD_B: [f64; SIZE] = [
    0.97901834, 0.97901649, 0.97901118, 0.97892146, 0.97858555, 0.97743705, 0.97428075, 0.96663223,
    0.94822893, 0.89937713, 0.76070164, 0.4642044, 0.20123039, 0.08808402, 0.04592894, 0.02860373,
    0.02060067, 0.01656701, 0.01451549, 0.01357964, 0.01331243, 0.01347661, 0.01387181, 0.01435472,
    0.01479836, 0.0151525, 0.01540513, 0.01557233, 0.0156571, 0.01571025, 0.01571916, 0.01572133,
    0.01572502, 0.01571717, 0.01571905, 0.01571059, 0.01569728, 0.0157002,
];

/// CIE Color Matching Function for X component.
#[allow(clippy::unreadable_literal)]
pub const CIE_CMF_X: [f64; SIZE] = [
    0.00006469, 0.00021941, 0.00112057, 0.00376661, 0.01188055, 0.02328644, 0.03455942, 0.03722379,
    0.03241838, 0.02123321, 0.01049099, 0.00329584, 0.00050704, 0.00094867, 0.00627372, 0.01686462,
    0.02868965, 0.04267481, 0.05625475, 0.0694704, 0.08305315, 0.0861261, 0.09046614, 0.08500387,
    0.07090667, 0.05062889, 0.03547396, 0.02146821, 0.01251646, 0.00680458, 0.00346457, 0.00149761,
    0.0007697, 0.00040737, 0.00016901, 0.00009522, 0.00004903, 0.00002,
];

/// CIE Color Matching Function for Y component.
#[allow(clippy::unreadable_literal)]
pub const CIE_CMF_Y: [f64; SIZE] = [
    0.00000184, 0.00000621, 0.00003101, 0.00010475, 0.00035364, 0.00095147, 0.00228226, 0.00420733,
    0.00668880, 0.00988840, 0.01524945, 0.02141831, 0.03342293, 0.05131001, 0.07040208, 0.08783871,
    0.09424905, 0.09795667, 0.09415219, 0.08678102, 0.07885653, 0.06352670, 0.05374142, 0.04264606,
    0.03161735, 0.02088521, 0.01386011, 0.00810264, 0.00463010, 0.00249138, 0.00125930, 0.00054165,
    0.00027795, 0.00014711, 0.00006103, 0.00003439, 0.00001771, 0.00000722,
];

/// CIE Color Matching Function for Z component.
#[allow(clippy::unreadable_literal)]
pub const CIE_CMF_Z: [f64; SIZE] = [
    0.00030502, 0.00103681, 0.00531314, 0.01795439, 0.05707758, 0.11365162, 0.17335873, 0.19620658,
    0.18608237, 0.13995048, 0.08917453, 0.04789621, 0.02814563, 0.01613766, 0.00775910, 0.00429615,
    0.00200551, 0.00086147, 0.00036904, 0.00019143, 0.00014956, 0.00009231, 0.00006813, 0.00002883,
    0.00001577, 0.00000394, 0.00000158, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
];

/// Transformation matrix to convert from the XYZ color space to the sRGB color space.
#[rustfmt::skip]
pub const XYZ_RGB: [[f64; 3]; 3] = [
    [ 3.243_063_33, -1.538_376_19, -0.498_932_82],
    [-0.968_963_09,  1.875_424_51,  0.041_543_03],
    [ 0.055_683_92, -0.204_174_38,  1.057_994_54],
];

/* Type definitions */

/// Represents the color format, matching the predefined constants RGB, RGBA, HEX, and HEXA.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum ColorFormat {
    Rgb,
    Rgba,
    Hex,
    Hexa,
}

/* Core color mixing functions */

/// Computes the linear concentration between two luminances.
pub fn linear_to_concentration(l1: f64, l2: f64, t: f64) -> f64 {
    let t1 = l1 * (1.0 - t).powi(2);
    let t2 = l2 * t.powi(2);
    t2 / (t1 + t2)
}

/// Mixes two colors using spectral upsampling and outputs the result in the specified format.
#[cfg(feature = "alloc")]
#[cfg_attr(feature = "nightly_doc", doc(cfg(feature = "alloc")))]
pub fn mix(color1: &str, color2: &str, t: f64, return_format: ColorFormat) -> String {
    let color1 = unpack(color1);
    let color2 = unpack(color2);

    let lrgb1 = srgb_to_linear(&color1[..3]);
    let lrgb2 = srgb_to_linear(&color2[..3]);

    let r1 = spectral_linear_to_spd(&lrgb1);
    let r2 = spectral_linear_to_spd(&lrgb2);

    let l1 = dot_product(&r1, &CIE_CMF_Y);
    let l2 = dot_product(&r2, &CIE_CMF_Y);

    let t = linear_to_concentration(l1, l2, t);

    let mut r = [0.0; SIZE];

    for i in 0..SIZE {
        let ks = (1.0 - t) * ((1.0 - r1[i]).powi(2) / (2.0 * r1[i]))
            + t * ((1.0 - r2[i]).powi(2) / (2.0 * r2[i]));
        let km = 1.0 + ks - (ks * ks + 2.0 * ks).sqrt();

        // Saunderson correction:
        // let S = ((1.0 - K1) * (1.0 - K2) * KM) / (1.0 - K2 * KM);

        r[i] = km;
    }

    let xyz = reflectance_to_xyz(&r);
    let rgb = xyz_to_srgb(&xyz);
    let rgba = [rgb[0], rgb[1], rgb[2], lerp(color1[3], color2[3], t)];
    pack(&rgba, return_format)
}

/// Generates a palette by interpolating between two colors.
#[cfg(feature = "alloc")]
#[cfg_attr(feature = "nightly_doc", doc(cfg(feature = "alloc")))]
pub fn palette(color1: &str, color2: &str, size: usize, return_format: ColorFormat) -> Vec<String> {
    (0..size)
        .map(|i| mix(color1, color2, i as f64 / (size as f64 - 1.0), return_format))
        .collect()
}

/* Color space conversion functions */

/// Applies the sRGB uncompanding function.
pub fn uncompand(x: f64) -> f64 {
    if x < 0.04045 {
        x / 12.92
    } else {
        ((x + 0.055) / 1.055).powf(GAMMA)
    }
}

/// Applies the sRGB companding function.
pub fn compand(x: f64) -> f64 {
    if x < 0.003_130_8 {
        x * 12.92
    } else {
        1.055 * x.powf(1.0 / GAMMA) - 0.055
    }
}

/// Converts sRGB to linear RGB values.
pub fn srgb_to_linear(srgb: &[f64]) -> [f64; 3] {
    [uncompand(srgb[0] / 255.0), uncompand(srgb[1] / 255.0), uncompand(srgb[2] / 255.0)]
}

/// Converts linear RGB values to clamped and companded sRGB values.
pub fn linear_to_srgb(lrgb: &[f64]) -> [f64; 3] {
    [
        (compand(lrgb[0]).clamp(0.0, 1.0) * 255.0).round(),
        (compand(lrgb[1]).clamp(0.0, 1.0) * 255.0).round(),
        (compand(lrgb[2]).clamp(0.0, 1.0) * 255.0).round(),
    ]
}

/// Converts XYZ to sRGB.
pub fn xyz_to_srgb(xyz: &[f64]) -> [f64; 3] {
    let r = dot_product(&XYZ_RGB[0], xyz);
    let g = dot_product(&XYZ_RGB[1], xyz);
    let b = dot_product(&XYZ_RGB[2], xyz);
    linear_to_srgb(&[r, g, b])
}

/// Converts reflectance to XYZ values.
pub fn reflectance_to_xyz(r: &[f64]) -> [f64; 3] {
    [dot_product(r, &CIE_CMF_X), dot_product(r, &CIE_CMF_Y), dot_product(r, &CIE_CMF_Z)]
}

/* Spectral conversion functions */

/// Computes spectral upsampling weights.
pub fn spectral_upsampling(lrgb: &[f64; 3]) -> [f64; 7] {
    let w = lrgb.iter().copied().fold(f64::INFINITY, f64::min);
    let lrgb = [lrgb[0] - w, lrgb[1] - w, lrgb[2] - w];

    let c = lrgb[1].min(lrgb[2]);
    let m = lrgb[0].min(lrgb[2]);
    let y = lrgb[0].min(lrgb[1]);
    let r = 0.0f64.max((lrgb[0] - lrgb[2]).min(lrgb[0] - lrgb[1]));
    let g = 0.0f64.max((lrgb[1] - lrgb[2]).min(lrgb[1] - lrgb[0]));
    let b = 0.0f64.max((lrgb[2] - lrgb[1]).min(lrgb[2] - lrgb[0]));

    [w, c, m, y, r, g, b]
}

/// Converts linear RGB values to reflectance values based on spectral upsampling weights.
// NOTE: https://github.com/rvanwijnen/spectral.js/issues/13#issuecomment-2144502129
pub fn spectral_linear_to_spd(lrgb: &[f64; 3]) -> [f64; SIZE] {
    let weights = spectral_upsampling(lrgb);

    ::core::array::from_fn(|i| {
        (weights[0]
            + weights[1] * SPD_C[i]
            + weights[2] * SPD_M[i]
            + weights[3] * SPD_Y[i]
            + weights[4] * SPD_C[i]
            + weights[5] * SPD_M[i]
            + weights[6] * SPD_Y[i])
            .max(EPSILON)
    })
}

/* Utility functions */

/// Computes the dot product of two vectors.
pub fn dot_product(a: &[f64], b: &[f64]) -> f64 {
    a.iter().zip(b.iter()).map(|(x, y)| x * y).sum()
}

/// Converts a color to GLSL-friendly RGBA format.
pub const fn glsl_color(c: &[f64; 4]) -> [f64; 4] {
    [c[0] / 255.0, c[1] / 255.0, c[2] / 255.0, if c[3] > 1.0 { c[3] / 255.0 } else { c[3] }]
}

/// Unpacks a color string (e.g., "#rrggbb" or "rgba(...)" format) into an RGBA array.
#[cfg(feature = "alloc")]
#[cfg_attr(feature = "nightly_doc", doc(cfg(feature = "alloc")))]
pub fn unpack(color: &str) -> [f64; 4] {
    if color.starts_with("rgb") {
        let values: Vec<f64> = color
            .split(['(', ',', ')'].as_ref())
            .filter_map(|s| {
                let s = s.trim();
                if let Some(stripped) = s.strip_suffix('%') {
                    stripped.parse::<f64>().ok().map(|v| v * 2.55)
                } else {
                    s.parse::<f64>().ok()
                }
            })
            .collect();

        let (r, g, b, a) = (values[0], values[1], values[2], *values.get(3).unwrap_or(&1.0));
        return [r, g, b, a];
    }

    if let Some(stripped) = color.strip_prefix('#') {
        let hex = if color.len() == 4 || color.len() == 5 {
            color.chars().skip(1).flat_map(|c| core::iter::repeat(c).take(2)).collect::<String>()
        } else {
            stripped.to_string()
        };

        let values: Vec<f64> = hex
            .as_bytes()
            .chunks(2)
            .map(|chunk| {
                let hex = ::core::str::from_utf8(chunk).unwrap();
                u8::from_str_radix(hex, 16).unwrap() as f64
            })
            .collect();

        let (r, g, b, a) = (values[0], values[1], values[2], *values.get(3).unwrap_or(&255.0));
        return [r, g, b, a / 255.0];
    }

    [0.0, 0.0, 0.0, 1.0]
}

/// Packs a color into a specific format (e.g., HEX, RGBA).
#[cfg(feature = "alloc")]
#[cfg_attr(feature = "nightly_doc", doc(cfg(feature = "alloc")))]
pub fn pack(srgb: &[f64; 4], return_format: ColorFormat) -> String {
    let (r, g, b, a) = (srgb[0], srgb[1], srgb[2], 1.0);
    use ColorFormat::{Hex, Hexa, Rgb, Rgba};

    match return_format {
        Rgb => format!("rgba({}, {}, {}, {:.3})", r, g, b, a),
        Rgba => format!("rgb({}, {}, {})", r, g, b),
        Hex | Hexa => {
            let r = format!("{:02x}", r as u8);
            let g = format!("{:02x}", g as u8);
            let b = format!("{:02x}", b as u8);
            let a = format!("{:02x}", (a * 255.0_f64).round() as u8);

            if return_format == Hexa {
                format!("#{}{}{}{}", r, g, b, a)
            } else {
                format!("#{}{}{}", r, g, b)
            }
        }
    }
}

/// Linear interpolation.
fn lerp(a: f64, b: f64, t: f64) -> f64 {
    a * (1.0 - t) + b * t
}