Implement sRGB transfer functions

jxl/src/color/mod.rs

@@ -0,0 +1,6 @@
+// Copyright (c) the JPEG XL Project Authors. All rights reserved.
+//
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+pub mod tf;

jxl/src/color/tf.rs

@@ -0,0 +1,149 @@
+// Copyright (c) the JPEG XL Project Authors. All rights reserved.
+//
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+const SRGB_POWTABLE_UPPER: [u8; 16] = [
+    0x00, 0x0a, 0x19, 0x26, 0x32, 0x41, 0x4d, 0x5c, 0x68, 0x75, 0x83, 0x8f, 0xa0, 0xaa, 0xb9, 0xc6,
+];
+
+const SRGB_POWTABLE_LOWER: [u8; 16] = [
+    0x00, 0xb7, 0x04, 0x0d, 0xcb, 0xe7, 0x41, 0x68, 0x51, 0xd1, 0xeb, 0xf2, 0x00, 0xb7, 0x04, 0x0d,
+];
+
+/// Converts the linear samples with the sRGB transfer curve.
+// Fast linear to sRGB conversion, ported from libjxl. Max error ~1.2e-4
+pub fn linear_to_srgb_fast(samples: &mut [f32]) {
+    for s in samples {
+        let v = s.to_bits() & 0x7fff_ffff;
+        let v_adj = f32::from_bits((v | 0x3e80_0000) & 0x3eff_ffff);
+        let pow = 0.059914046f32;
+        let pow = pow * v_adj - 0.10889456;
+        let pow = pow * v_adj + 0.107963754;
+        let pow = pow * v_adj + 0.018092343;
+
+        // `mul` won't be used when `v` is small.
+        let idx = (v >> 23).wrapping_sub(118) as usize & 0xf;
+        let mul = 0x4000_0000
+            | (u32::from(SRGB_POWTABLE_UPPER[idx]) << 18)
+            | (u32::from(SRGB_POWTABLE_LOWER[idx]) << 10);
+
+        let v = f32::from_bits(v);
+        let small = v * 12.92;
+        let acc = pow * f32::from_bits(mul) - 0.055;
+
+        *s = if v <= 0.0031308 { small } else { acc }.copysign(*s);
+    }
+}
+
+/// Converts the linear samples with the sRGB transfer curve.
+// Max error ~5e-7
+pub fn linear_to_srgb(samples: &mut [f32]) {
+    #[allow(clippy::excessive_precision)]
+    const P: [f32; 5] = [
+        -5.135152395e-4,
+        5.287254571e-3,
+        3.903842876e-1,
+        1.474205315,
+        7.352629620e-1,
+    ];
+
+    #[allow(clippy::excessive_precision)]
+    const Q: [f32; 5] = [
+        1.004519624e-2,
+        3.036675394e-1,
+        1.340816930,
+        9.258482155e-1,
+        2.424867759e-2,
+    ];
+
+    for x in samples {
+        let a = x.abs();
+        *x = if a <= 0.0031308 {
+            a * 12.92
+        } else {
+            crate::util::eval_rational_poly(a.sqrt(), P, Q)
+        }
+        .copysign(*x);
+    }
+}
+
+/// Converts samples in sRGB transfer curve to linear. Inverse of `linear_to_srgb`.
+pub fn srgb_to_linear(samples: &mut [f32]) {
+    #[allow(clippy::excessive_precision)]
+    const P: [f32; 5] = [
+        2.200248328e-4,
+        1.043637593e-2,
+        1.624820318e-1,
+        7.961564959e-1,
+        8.210152774e-1,
+    ];
+
+    #[allow(clippy::excessive_precision)]
+    const Q: [f32; 5] = [
+        2.631846970e-1,
+        1.076976492,
+        4.987528350e-1,
+        -5.512498495e-2,
+        6.521209011e-3,
+    ];
+
+    for x in samples {
+        let a = x.abs();
+        *x = if a <= 0.04045 {
+            a / 12.92
+        } else {
+            crate::util::eval_rational_poly(a, P, Q)
+        }
+        .copysign(*x);
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use test_log::test;
+
+    use crate::util::test::assert_all_almost_eq;
+    use super::*;
+
+    fn arb_samples(
+        u: &mut arbtest::arbitrary::Unstructured,
+    ) -> arbtest::arbitrary::Result<Vec<f32>> {
+        let len = u.arbitrary_len::<u32>()?;
+        let mut samples = Vec::with_capacity(len);
+        // uniform distribution in [-1.0, 1.0]
+        for _ in 0..len {
+            let a: u32 = u.int_in_range(0..=(1 << 24))?;
+            let signed: bool = u.arbitrary()?;
+            let x = (a + 1) as f32 / (1 << 24) as f32;
+            samples.push(if signed { -x } else { x });
+        }
+        Ok(samples)
+    }
+
+    #[test]
+    fn srgb_roundtrip_arb() {
+        arbtest::arbtest(|u| {
+            let samples: Vec<f32> = arb_samples(u)?;
+            let mut output = samples.clone();
+
+            linear_to_srgb(&mut output);
+            srgb_to_linear(&mut output);
+            assert_all_almost_eq!(&output, &samples, 2e-6);
+            Ok(())
+        });
+    }
+
+    #[test]
+    fn linear_to_srgb_fast_arb() {
+        arbtest::arbtest(|u| {
+            let mut samples: Vec<f32> = arb_samples(u)?;
+            let mut fast = samples.clone();
+
+            linear_to_srgb(&mut samples);
+            linear_to_srgb_fast(&mut fast);
+            assert_all_almost_eq!(&samples, &fast, 1.2e-4);
+            Ok(())
+        });
+    }
+}

jxl/src/lib.rs

@@ -5,6 +5,7 @@
 
 #![deny(unsafe_code)]
 pub mod bit_reader;
+pub mod color;
 pub mod container;
 pub mod entropy_coding;
 pub mod error;

jxl/src/util.rs

@@ -10,6 +10,7 @@ mod bits;
 mod concat_slice;
 mod linalg;
 mod log2;
+mod rational_poly;
 mod shift_right_ceil;
 pub mod tracing_wrappers;
 mod vec_helpers;
@@ -18,5 +19,6 @@ pub use bits::*;
 pub use concat_slice::*;
 pub use linalg::*;
 pub use log2::*;
+pub use rational_poly::*;
 pub use shift_right_ceil::*;
 pub use vec_helpers::*;

jxl/src/util/rational_poly.rs

@@ -0,0 +1,10 @@
+/// Computes `(p0 + p1 x + p2 x^2 + ...) / (q0 + q1 x + q2 x^2 + ...)`.
+///
+/// # Panics
+/// Panics if either `P` or `Q` is zero.
+#[inline]
+pub fn eval_rational_poly<const P: usize, const Q: usize>(x: f32, p: [f32; P], q: [f32; Q]) -> f32 {
+    let yp = p.into_iter().rev().reduce(|yp, p| yp * x + p).unwrap();
+    let yq = q.into_iter().rev().reduce(|yq, q| yq * x + q).unwrap();
+    yp / yq
+}