srgb 0.3.4

//! Functions converting a colour into a shade of grey.  All function is this
//! module and submodules return the shade of grey as a single component in the
//! same colour space the source colour was in.  To construct actual grey
//! colour, the returned component simply needs to be repeated three times in an
//! arary.
//!
//! With the exception of [`from_linear()`] function — which operates in
//! linear sRGB space — all functions in this module (but not submodules)
//! perform proper gamma calculations.  This guarantees accurate result but may
//! be needlessly slow for some applications.  To address that, a faster,
//! approximate algorithms are are provided in [`approx_gamma`] and [`no_gamma`]
//! submodules.
//!
//! For more information and comparisons of accuracy and speed of various
//! methods see [Greyscale, you might be doing it
//! wrong](https://mina86.com/2021/rgb-to-greyscale/).

/// Calculates shade of grey for a colour given in linear sRGB space.  This is
/// morally equivalent to converting the colour to XYZ and returning the
/// Y component.  Because the colour is already in linear space, this function
/// does not need to perform any gamma calculations and thus is fast (boiling
/// down to only three floating point multiplies and additions).
///
/// # Example
/// ```
/// assert_eq!(0.0       , srgb::grey::from_linear([0.0, 0.0, 0.0]));
/// assert_eq!(0.22809356, srgb::grey::from_linear([0.4, 0.2, 0.0]));
/// assert_eq!(0.17190644, srgb::grey::from_linear([0.0, 0.2, 0.4]));
/// assert_eq!(0.19255298, srgb::grey::from_linear([0.6, 0.0, 0.9]));
/// assert_eq!(1.0       , srgb::grey::from_linear([1.0, 1.0, 1.0]));
/// ```
pub fn from_linear(linear: [f32; 3]) -> f32 { from_floats(linear) }

/// Converts an sRGB colour in normalised representation into greyscale
/// returning the shade of grey as a normalised number in the range 0–1.  This
/// perorms gamma calculation first converting the colour to linear space and
/// later gamma compressing the resulting shade of grey.
///
/// # Example
/// ```
/// assert_eq!(0.0       , srgb::grey::from_normalised([0.0, 0.0, 0.0]));
/// assert_eq!(0.2526165 , srgb::grey::from_normalised([0.4, 0.2, 0.0]));
/// assert_eq!(0.2005172 , srgb::grey::from_normalised([0.0, 0.2, 0.4]));
/// assert_eq!(0.38794443, srgb::grey::from_normalised([0.6, 0.0, 0.9]));
/// // Floating point Maths is unfortunately hard:
/// assert_eq!(0.99999994, srgb::grey::from_normalised([1.0, 1.0, 1.0]));
/// ```
pub fn from_normalised(normalised: [f32; 3]) -> f32 {
    let y = from_linear(super::linear_from_normalised(normalised));
    super::gamma::compress_normalised(y)
}

/// Converts a 24-bit sRGB colour (also known as true colour) into greyscale
/// returning the shade of grey as a single 8-bit integer.  This perorms gamma
/// calculation first converting the colour to linear space and later gamma
/// compressing the resulting shade of grey.
///
/// # Example
/// ```
/// assert_eq!(  0, srgb::grey::from_u8([  0,   0,   0]));
/// assert_eq!( 62, srgb::grey::from_u8([ 42,  69,   0]));
/// assert_eq!( 40, srgb::grey::from_u8([  0,  42,  69]));
/// assert_eq!( 27, srgb::grey::from_u8([  1,  33,   7]));
/// assert_eq!(255, srgb::grey::from_u8([255, 255, 255]));
/// ```
pub fn from_u8(encoded: [u8; 3]) -> u8 {
    super::gamma::compress_u8(from_linear(super::linear_from_u8(encoded)))
}


/// Functions converting a colour into a shade of grey which trade accuracy of
/// the conversion for speed.  This is done by approximating the gamma
/// calculations.  Specifically, functions in this module use γ = 2 which
/// simplifies and speeds up all the formulæ.
///
/// Note: This module does not include a `from_linear` function since that
/// function does not need to worry about gamma at all.  If you have a colour
/// given in linear sRGB space simply use [`grey::from_linear()`] function.
pub mod approx_gamma {
    /// Converts an sRGB colour in normalised representation into greyscale
    /// using γ = 2 approximation for gamma correction rather than using proper
    /// sRGB gamma compression and expansion functions.
    ///
    /// # Example
    /// ```
    /// use srgb::grey::approx_gamma::from_normalised;
    ///
    /// assert_eq!(0.0       , from_normalised([0.0, 0.0, 0.0]));
    /// assert_eq!(0.2502612 , from_normalised([0.4, 0.2, 0.0]));
    /// assert_eq!(0.20038915, from_normalised([0.0, 0.2, 0.4]));
    /// assert_eq!(0.36745176, from_normalised([0.6, 0.0, 0.9]));
    /// assert_eq!(1.0       , from_normalised([1.0, 1.0, 1.0]));
    /// ```
    pub fn from_normalised(normalised: [f32; 3]) -> f32 {
        let [r, g, b] = normalised;
        super::from_floats([r * r, g * g, b * b]).sqrt()
    }

    /// Converts a 24-bit sRGB colour (also known as true colour) into greyscale
    /// returning the shade of grey as a single 8-bit integer.  Uses γ = 2
    /// approximation for gamma correction rather than using proper sRGB gamma
    /// compression and expansion functions.
    ///
    /// # Example
    /// ```
    /// use srgb::grey::approx_gamma::from_u8;
    ///
    /// assert_eq!(  0, from_u8([  0,   0,   0]));
    /// assert_eq!( 58, from_u8([ 42,  69,   0]));
    /// assert_eq!( 35, from_u8([  0,  42,  69]));
    /// assert_eq!( 24, from_u8([  1,  33,   7]));
    /// assert_eq!(255, from_u8([255, 255, 255]));
    /// ```
    pub fn from_u8(encoded: [u8; 3]) -> u8 {
        let [r, g, b] = encoded;
        (super::from_floats([r as f32, g as f32, b as f32]) + 0.5) as u8
    }

    fn isqrt(n: u32) -> u8 {
        let mut n = n;
        let mut x = 0;
        let mut b = 1 << 16;

        while b > n {
            b = b / 4;
        }

        while b != 0 {
            if n >= x + b {
                n = n - x - b;
                x = x / 2 + b;
            } else {
                x = x / 2;
            }
            b = b / 4;
        }
        x as u8
    }

    /// Converts a 24-bit sRGB colour (also known as true colour) into greyscale
    /// returning the shade of grey as a single 8-bit integer.  Uses γ = 2
    /// approximation for gamma correction (rather than using proper sRGB gamma
    /// compression and expansion functions) and furthermore performs all
    /// calculations using integer arithmetic.  This has better precision than
    /// [`from_u8()`] function (which uses FPU) but is slower (on systems with
    /// FPU).
    ///
    /// # Example
    /// ```
    /// use srgb::grey::approx_gamma::from_u8_no_fpu;
    ///
    /// assert_eq!(  0, from_u8_no_fpu([  0,   0,   0]));
    /// assert_eq!( 61, from_u8_no_fpu([ 42,  69,   0]));
    /// assert_eq!( 40, from_u8_no_fpu([  0,  42,  69]));
    /// assert_eq!( 27, from_u8_no_fpu([  1,  33,   7]));
    /// assert_eq!(255, from_u8_no_fpu([255, 255, 255]));
    /// ```
    pub fn from_u8_no_fpu(encoded: [u8; 3]) -> u8 {
        let [r, g, b] = encoded;
        isqrt(
            (r as u32 * r as u32 * 13936 +
                g as u32 * g as u32 * 46869 +
                b as u32 * b as u32 * 4731) >>
                16,
        )
    }
}

/// Functions converting a colour into a shade of grey which trade accuracy of
/// the conversion for speed.  This is done by forgoing any gamma calculations.
/// I generally don’t recommend using those functions unless you’ve determined
/// that greyscale conversion is the bottleneck of your application and you are
/// happy with severe precision degradation.
///
/// Note: This module does not include a `from_linear` function since that
/// function does not need to worry about gamma at all.  If you have a colour
/// given in linear sRGB space simply use [`grey::from_linear()`] function.
pub mod no_gamma {
    /// Calculates shade of grey for a colour given in linear RGB space but
    /// without performing any gamma correction.  This makes the function
    /// especially fast but at the same time very inaccurate.
    ///
    /// # Example
    /// ```
    /// use srgb::grey::no_gamma::from_normalised;
    ///
    /// assert_eq!(0.0       , from_normalised([0.0, 0.0, 0.0]));
    /// assert_eq!(0.22809356, from_normalised([0.4, 0.2, 0.0]));
    /// assert_eq!(0.17190644, from_normalised([0.0, 0.2, 0.4]));
    /// assert_eq!(0.19255298, from_normalised([0.6, 0.0, 0.9]));
    /// assert_eq!(1.0       , from_normalised([1.0, 1.0, 1.0]));
    /// ```
    pub fn from_normalised(normalised: [f32; 3]) -> f32 {
        super::from_floats(normalised)
    }

    /// Converts a 24-bit sRGB colour (also known as true colour) into greyscale
    /// returning the shade of grey as a single 8-bit integer.  Does not perform
    /// any gamma correction nor any floating point arithmetic.  This makes the
    /// function especially fast but at the same time very inaccurate.
    ///
    /// # Example
    /// ```
    /// use srgb::grey::no_gamma::from_u8;
    ///
    /// assert_eq!(  0, from_u8([  0,   0,   0]));
    /// assert_eq!( 58, from_u8([ 42,  69,   0]));
    /// assert_eq!( 35, from_u8([  0,  42,  69]));
    /// assert_eq!( 24, from_u8([  1,  33,   7]));
    /// assert_eq!(255, from_u8([255, 255, 255]));
    /// ```
    pub fn from_u8(encoded: [u8; 3]) -> u8 {
        let [r, g, b] = encoded;
        let y = 3567454 * r as u32 + 11998779 * g as u32 + 1210983 * b as u32;
        ((y + (1 << 23)) >> 24) as u8
    }
}


fn from_floats(arr: [f32; 3]) -> f32 {
    super::dot_product(&super::xyz::XYZ_FROM_SRGB_MATRIX[1], &arr)
}


#[cfg(test)]
mod test {
    use approx::assert_ulps_eq;

    #[test]
    fn test_is_grey() {
        for i in 0..=255 {
            assert_eq!(i, super::from_u8([i, i, i]));
            assert_eq!(i, super::approx_gamma::from_u8([i, i, i]));
            assert_eq!(i, super::approx_gamma::from_u8_no_fpu([i, i, i]));
            assert_eq!(i, super::no_gamma::from_u8([i, i, i]));

            let i = i as f32 / 255.0;
            assert_ulps_eq!(i, super::from_normalised([i, i, i]), max_ulps = 1);
            assert_ulps_eq!(i, super::from_linear([i, i, i]), max_ulps = 1);
            assert_ulps_eq!(
                i,
                super::approx_gamma::from_normalised([i, i, i]),
                max_ulps = 1
            );
            assert_ulps_eq!(
                i,
                super::no_gamma::from_normalised([i, i, i]),
                max_ulps = 1
            );
        }
    }

    fn distance(want_grey: f32, got_grey: f32) -> f64 {
        fn l_from_y(y: f64) -> f64 {
            // http://www.brucelindbloom.com/index.html?Eqn_XYZ_to_Lab.html
            if y > 216.0 / 24389.0 {
                116.0 * y.powf(1.0 / 3.0) - 16.0
            } else {
                24389.0 / 27.0 * y
            }
        }

        // ΔE₀₀ implementation with assumption that a* and b* for both colours
        // are zero.  I.e. compares two shades of grey.
        fn de(l1: f64, l2: f64) -> f64 {
            let v = ((l1 + l2) / 2.0 - 50.0).powi(2);
            let div = 1.0 + ((0.015 * v) / (20.0 + v).sqrt());
            (l2 - l1).abs() / div
        }

        de(l_from_y(want_grey as f64), l_from_y(got_grey as f64))
    }

    #[derive(Debug, PartialEq)]
    struct Bench(f64, f64);

    impl approx::AbsDiffEq for Bench {
        type Epsilon = <f64 as approx::AbsDiffEq>::Epsilon;

        fn default_epsilon() -> Self::Epsilon { 0.000001 }
        fn abs_diff_eq(&self, rhs: &Self, epsilon: Self::Epsilon) -> bool {
            self.0.abs_diff_eq(&rhs.0, epsilon) &&
                self.1.abs_diff_eq(&rhs.1, epsilon)
        }
    }

    impl approx::UlpsEq for Bench {
        fn default_max_ulps() -> u32 { f64::default_max_ulps() }
        fn ulps_eq(
            &self,
            rhs: &Self,
            epsilon: Self::Epsilon,
            max_ulps: u32,
        ) -> bool {
            self.0.ulps_eq(&rhs.0, epsilon, max_ulps) &&
                self.1.ulps_eq(&rhs.1, epsilon, max_ulps)
        }
    }

    fn accuracy(f: &dyn Fn([u8; 3]) -> f64) -> Bench {
        // Calculate statistics
        let mut sum_err = 0.0;
        let mut max_err = 0.0;
        for i in 0..(1 << 24) {
            let rgb = [(i >> 16) as u8, (i >> 8) as u8, (i) as u8];
            let err = f(rgb);
            max_err = err.max(max_err);
            sum_err += err;
        }
        Bench(sum_err, max_err)
    }

    #[test]
    fn test_from_normalised_accuracy() {
        let cases: [(f64, f64, &dyn Fn([f32; 3]) -> f32); 3] = [
            (0.0, 0.0, &crate::grey::from_normalised),
            (
                11278361.255885385,
                3.860478819921189,
                &crate::grey::approx_gamma::from_normalised,
            ),
            (
                69142133.11431149,
                28.150428835752255,
                &crate::grey::no_gamma::from_normalised,
            ),
        ];
        for (want_sum, want_max, func) in cases.iter().copied() {
            let got = accuracy(&|rgb| {
                let rgb = crate::normalised_from_u8(rgb);
                let want = crate::grey::from_normalised(rgb);
                let got = func(rgb);
                distance(
                    crate::gamma::expand_normalised(want),
                    crate::gamma::expand_normalised(got),
                )
            });
            assert_ulps_eq!(Bench(want_sum, want_max), got);
        }
    }

    #[test]
    fn test_from_u8_accuracy() {
        #[rustfmt::skip]
        let cases: [(u64, u64, &dyn Fn([u8; 3]) -> u8); 4] = [
            (        0,  0, &crate::grey::from_u8),
            (207766613, 73, &crate::grey::approx_gamma::from_u8),
            ( 41742745, 11, &crate::grey::approx_gamma::from_u8_no_fpu),
            (207766715, 73, &crate::grey::no_gamma::from_u8),
        ];
        for (i, case) in cases.iter().copied().enumerate() {
            let (want_sum, want_max, func) = case;
            let got = accuracy(&|rgb| {
                let lin = crate::grey::from_linear(crate::linear_from_u8(rgb));
                let want = crate::gamma::compress_u8(lin);
                let got = func(rgb);
                (want as i32 - got as i32).abs() as f64
            });
            assert_eq!(
                (want_sum, want_max),
                (got.0 as u64, got.1 as u64),
                " [Test Case #{}]",
                i
            );
        }
    }
}