oklab 1.1.2

OK Lab is a perceptually uniform color space for image processing. Improvement of CIELAB.
Documentation
//! [What is Oklab color space?](https://bottosson.github.io/posts/oklab/)
#![allow(clippy::excessive_precision)]
#![allow(clippy::inline_always)]
#![allow(clippy::similar_names)]
#![allow(clippy::many_single_char_names)]
#![allow(clippy::cast_possible_truncation)]

use rgb::ComponentMap;

#[doc(hidden)]
#[deprecated(note = "Rename RGB to Rgb")]
pub type RGB<T> = rgb::Rgb<T>;

/// Interoperable with the [`rgb` crate](https://lib.rs/crates/rgb).
pub use rgb::Rgb;

/// Linearized sRGB (AKA linear light RGB, sRGB with gamma 1.0)
pub type LinearRgb = Rgb<f32>;

/// A color in Oklab is represented with three coordinates, similar to how CIELAB works, but with better perceptual properties.
///
/// Oklab uses a D65 whitepoint, since this is what sRGB and other common color spaces use.
#[derive(Copy, Clone, Debug, PartialOrd, PartialEq)]
#[repr(C)]
pub struct Oklab {
    /// L – perceived lightness
    pub l: f32,
    /// a – how green/red the color is
    pub a: f32,
    /// b – how blue/yellow the color is
    pub b: f32,
}

/// Converts from _linearized_ sRGB (in 0..1 range) to [`Oklab`]
///
/// (this is what is usually called a "linear light" or "without gamma" RGB.)
#[must_use]
pub fn linear_srgb_to_oklab(c: LinearRgb) -> Oklab {
    let l = 0.0514459929f32.mul_add(c.b, 0.4122214708f32.mul_add(c.r, 0.5363325363 * c.g));
    let m = 0.1073969566f32.mul_add(c.b, 0.2119034982f32.mul_add(c.r, 0.6806995451 * c.g));
    let s = 0.6299787005f32.mul_add(c.b, 0.0883024619f32.mul_add(c.r, 0.2817188376 * c.g));
    let l_ = cbrt(l);
    let m_ = cbrt(m);
    let s_ = cbrt(s);
    Oklab {
        l:(-0.0040720468f32).mul_add(s_, 0.2104542553f32.mul_add(l_, 0.7936177850 * m_)),
        a:  0.4505937099f32 .mul_add(s_, 1.9779984951f32.mul_add(l_,-2.4285922050 * m_)),
        b:(-0.8086757660f32).mul_add(s_, 0.0259040371f32.mul_add(l_, 0.7827717662 * m_)),
    }
}

/// (Optionally) a custom implementation of cube root
#[doc(hidden)]
#[inline(always)]
pub fn cbrt(x: f32) -> f32 {
    #[cfg(feature = "std-cbrt")]
    return x.cbrt();

    #[cfg(not(feature = "std-cbrt"))]
    {
        const B: u32 =  709957561;
        const C: f32 =  5.4285717010e-1;
        const D: f32 = -7.0530611277e-1;
        const E: f32 =  1.4142856598e+0;
        const F: f32 =  1.6071428061e+0;
        const G: f32 =  3.5714286566e-1;

        let mut t = f32::from_bits((x.to_bits() / 3).wrapping_add(B));
        let s = C + (t * t) * (t / x);
        t *= G + F / (s + E + D / s);
        t
    }
}

#[test]
fn cbrt_precision() {
    let mut max_err = 0.;
    let mut total_err = 0.;
    for i in 0..10000u16 {
        let i = (f64::from(i) / 10000.0) as f32;
        let good = (i as f64).cbrt();
        let approx = cbrt(i) as f64;
        let err = (good - approx).abs();
        max_err = err.max(max_err);
        total_err += err;
    }
    assert!(max_err < 0.000000190, "{max_err} {total_err}");
    assert!(total_err < 0.000308, "{max_err} {total_err}");
}

/// Converts [`Oklab`] to _linear_ sRGB (in 0..1 range)
///
/// (this is what is usually called a "linear light" or "without gamma" RGB.)
///
/// Outputs aren't clamped, and can be negative or slightly over 1.0 due to floating-point rounding errors.
#[must_use]
pub fn oklab_to_linear_srgb(c: Oklab) -> LinearRgb {
    let l_ =   0.2158037573f32 .mul_add(c.b,  0.3963377774f32 .mul_add(c.a, c.l));
    let m_ = (-0.0638541728f32).mul_add(c.b,(-0.1055613458f32).mul_add(c.a, c.l));
    let s_ = (-1.2914855480f32).mul_add(c.b,(-0.0894841775f32).mul_add(c.a, c.l));
    let l = l_ * l_ * l_;
    let m = m_ * m_ * m_;
    let s = s_ * s_ * s_;
    LinearRgb {
        r:  0.2309699292f32 .mul_add(s,   4.0767416621f32 .mul_add(l, -3.3077115913 * m)),
        g:(-0.3413193965f32).mul_add(s, (-1.2684380046f32).mul_add(l,  2.6097574011 * m)),
        b:  1.7076147010f32 .mul_add(s, (-0.0041960863f32).mul_add(l, -0.7034186147 * m)),
    }
}

#[cfg(feature = "fast-srgb8")]
use fast_srgb8::f32_to_srgb8;

#[inline(always)]
#[cfg(not(feature = "fast-srgb8"))]
fn f32_to_srgb8(u: f32) -> u8 {
    (to_gamma(u) * 255.0 + 0.5) as u8
}

#[cfg(feature = "fast-srgb8")]
use fast_srgb8::srgb8_to_f32;

#[inline]
#[cfg(not(feature = "fast-srgb8"))]
fn srgb8_to_f32(c: u8) -> f32 {
    const SRGB_LUT: [f32; 256] = [0.,0.000303527,0.000607054,0.000910581,0.001214108,0.001517635,0.0018211621,0.0021246891,0.0024282159,0.002731743,0.0030352699,0.0033465356,0.0036765069,0.004024717,0.0043914421,0.0047769533,0.0051815161,0.0056053917,0.0060488326,0.0065120901,0.0069954102,0.0074990317,0.0080231922,0.0085681248,0.009134057,0.0097212177,0.010329823,0.0109600937,0.0116122449,0.012286487,0.0129830306,0.0137020806,0.0144438408,0.0152085144,0.0159962922,0.0168073755,0.0176419523,0.0185002182,0.0193823576,0.0202885587,0.0212190095,0.0221738834,0.0231533647,0.0241576303,0.0251868572,0.0262412187,0.0273208879,0.0284260381,0.0295568425,0.0307134502,0.0318960398,0.0331047736,0.0343398117,0.0356013179,0.0368894525,0.0382043757,0.0395462364,0.0409151986,0.0423114225,0.0437350422,0.0451862141,0.0466650948,0.0481718332,0.0497065745,0.0512694679,0.0528606549,0.0544802807,0.0561284944,0.0578054339,0.0595112406,0.0612460561,0.063010022,0.0648032799,0.0666259527,0.068478182,0.0703601092,0.0722718611,0.0742135793,0.0761853904,0.0781874284,0.0802198276,0.0822827145,0.0843762159,0.0865004659,0.088655591,0.0908417106,0.093058981,0.0953074843,0.0975873619,0.0998987406,0.1022417471,0.104616493,0.1070231125,0.1094617173,0.1119324341,0.1144353822,0.1169706956,0.1195384488,0.1221387982,0.1247718409,0.1274376959,0.1301365197,0.1328683645,0.1356333643,0.1384316534,0.141263321,0.1441285014,0.1470272988,0.1499598175,0.1529261768,0.1559264958,0.1589608639,0.1620294005,0.1651322246,0.1682694256,0.1714411527,0.1746474504,0.1778884679,0.1811642945,0.1844750345,0.1878208071,0.1912017167,0.1946178675,0.1980693489,0.2015562952,0.2050787657,0.2086368948,0.2122307867,0.2158605307,0.2195262313,0.2232279778,0.2269658893,0.2307400703,0.2345505953,0.2383975834,0.2422811389,0.246201396,0.2501583695,0.2541521788,0.2581829131,0.2622507215,0.2663556635,0.2704978585,0.274677366,0.278894335,0.2831487954,0.287440896,0.2917706966,0.2961383164,0.3005438447,0.304987371,0.3094689548,0.3139887452,0.3185468316,0.3231432438,0.3277781308,0.332451582,0.3371636569,0.3419144452,0.3467040956,0.3515326381,0.356400162,0.3613067865,0.3662526011,0.3712377846,0.3762622178,0.3813261092,0.3864295185,0.3915725648,0.3967553079,0.4019778669,0.4072403014,0.4125427008,0.4178851545,0.4232677519,0.4286905527,0.4341537058,0.4396572411,0.4452012479,0.4507858455,0.4564110935,0.4620770514,0.4677838385,0.4735315442,0.4793202281,0.4851499796,0.4910208881,0.496933043,0.5028864741,0.5088813305,0.5149176717,0.5209956169,0.5271152258,0.5332764983,0.5394796729,0.54572469,0.5520116091,0.5583406091,0.5647116899,0.5711250305,0.5775806308,0.5840786099,0.5906190276,0.597202003,0.6038275361,0.6104957461,0.6172067523,0.6239605546,0.6307573318,0.6375970244,0.6444798708,0.6514058113,0.6583749652,0.6653874516,0.6724433303,0.6795426011,0.6866854429,0.6938719153,0.7011020184,0.7083759308,0.7156936526,0.7230552435,0.7304610014,0.7379106879,0.7454044819,0.7529424429,0.7605247498,0.7681514025,0.7758224607,0.7835380435,0.791298151,0.799102962,0.8069524765,0.8148468137,0.8227859735,0.8307700753,0.8387992382,0.8468734622,0.8549928069,0.8631573915,0.8713673353,0.8796225786,0.8879231811,0.8962695599,0.9046612382,0.9130988121,0.921581924,0.9301110506,0.9386857748,0.9473066926,0.955973506,0.9646865726,0.9734454751,0.9822508693,0.9911022186,1.];
    SRGB_LUT[c as usize]
}

#[inline]
fn to_gamma(u: f32) -> f32 {
    if u >= 0.0031308 {
        1.055_f32.mul_add(u.min(1.).powf(1. / 2.4), -0.055)
    } else {
        12.92 * u.max(0.)
    }
}

#[inline]
fn to_linear(u: f32) -> f32 {
    if u >= 0.04045 {
        ((u + 0.055) * (1. / (1. + 0.055))).min(1.).powf(2.4)
    } else {
        u.max(0.) * (1. / 12.92)
    }
}

/// Converts regular 8-bit s[RGB][rgb] color to [`Oklab`]
/// [rgb][rgb::RGB8]
#[inline]
#[must_use]
pub fn srgb_to_oklab(c: Rgb<u8>) -> Oklab {
    linear_srgb_to_oklab(c.map(srgb8_to_f32))
}

/// Converts sRGB as a float in range 0..=1 to [`Oklab`] (in usual Oklab range where a,b can be negative and below 1)
#[inline]
#[must_use]
pub fn srgb_f32_to_oklab(c: Rgb<f32>) -> Oklab {
    linear_srgb_to_oklab(c.map(to_linear))
}

/// Converts [`Oklab`] to regular 8-bit s[RGB][rgb] color
/// [rgb][rgb::RGB8]
#[inline]
#[must_use]
pub fn oklab_to_srgb(c: Oklab) -> Rgb<u8> {
    oklab_to_linear_srgb(c).map(f32_to_srgb8)
}

/// Converts [`Oklab`] to standard sRGB color as a float in range 0..=1
#[inline]
#[must_use]
pub fn oklab_to_srgb_f32(c: Oklab) -> Rgb<f32> {
    oklab_to_linear_srgb(c).map(to_gamma)
}

impl Oklab {
    /// Converts [`Oklab`] to regular 8-bit s[RGB][rgb] color
    /// [rgb][rgb::RGB8]
    #[inline(always)]
    #[must_use]
    pub fn to_srgb(&self) -> Rgb<u8> {
        oklab_to_srgb(*self)
    }

    /// Converts [`Oklab`] to sRGB color in range 0..=1
    #[inline(always)]
    #[must_use]
    pub fn to_srgb_f32(&self) -> Rgb<f32> {
        oklab_to_srgb_f32(*self)
    }

    /// Converts [`Oklab`] to _linear_ sRGB (in 0..1 range)
    ///
    /// (this is what is usually called a "linear light" or "without gamma" RGB.)
    ///
    /// Outputs aren't clamped, and can be negative or slightly over 1.0 due to floating-point rounding errors.
    #[inline(always)]
    #[must_use]
    pub fn to_linear_rgb(&self) -> LinearRgb {
        oklab_to_linear_srgb(*self)
    }

    /// Converts _linear_ sRGB (in 0..1 range) to [`Oklab`]
    ///
    /// (this is what is usually called a "linear light" or "without gamma" RGB.)
    #[inline(always)]
    #[must_use]
    pub fn from_linear_rgb(rgb: LinearRgb) -> Self {
        linear_srgb_to_oklab(rgb)
    }
}

/// Converts regular 8-bit s[RGB][rgb] color to [`Oklab`]
/// [rgb][rgb::RGB8]
impl From<Rgb<u8>> for Oklab {
    #[inline(always)]
    fn from(rgb: Rgb<u8>) -> Self {
        srgb_to_oklab(rgb)
    }
}

impl From<Oklab> for Rgb<u8> {
    #[inline(always)]
    fn from(lab: Oklab) -> Self {
        oklab_to_srgb(lab)
    }
}

#[test]
fn linear_int() {
    for i in 0..=255 {
        assert_eq!(i, f32_to_srgb8(srgb8_to_f32(i)));
    }
}


#[test]
fn linear_f() {
    for i in 0..4096u16 {
        let f = f32::from(i) / 4095.0;
        let lin = to_linear(f);
        let res = to_gamma(lin);
        let err = (res - f).abs();
        assert!(err < 0.001, "{f:0.2} => {lin:0.3} => {res:0.2}; diff {err:0.4}");
    }
}

#[test]
fn linear_mix() {
    for i in 0..4096u16 {
        let f = f32::from(i) / 4095.0;
        let lin = to_linear(f);
        let res = f32::from(f32_to_srgb8(lin)) / 255.0;
        let err = (res - f).abs();
        assert!(err < 1./255., "{f:0.3} ({:0.1}) => {lin:0.4} => {res:0.3}; diff {err:0.5}", (f * 255.));
    }
}

#[test]
fn roundtrip_half1() {
    for r in 128..=255 {
        for g in 0..=255 {
            for b in 0..=255 {
                let c = Rgb::new(r, g, b);
                assert_eq!(c, oklab_to_srgb(srgb_to_oklab(c)));
            }
        }
    }
}

#[test]
fn roundtrip_half2() {
    for r in 0..128 {
        for g in 0..=255 {
            for b in 0..=255 {
                let c = Rgb::new(r, g, b);
                assert_eq!(c, oklab_to_srgb(srgb_to_oklab(c)));
            }
        }
    }
}