Switch to colorbox and jakob upsampling for color handling.

Cargo.lock

@@ -126,6 +126,9 @@ dependencies = [
 [[package]]
 name = "color"
 version = "0.1.0"
+dependencies = [
+ "colorbox",
+]
 
 [[package]]
 name = "colorbox"

src/color.rs

@@ -1,13 +1,8 @@
 use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign};
 
 use crate::math::Float4;
-pub use color::{
-    rec709_e_to_xyz, rec709_to_xyz, xyz_to_aces_ap0, xyz_to_aces_ap0_e, xyz_to_rec709,
-    xyz_to_rec709_e,
-};
 use compact::fluv::fluv32;
 use half::f16;
-use spectral_upsampling::meng::{spectrum_xyz_to_p_4, EQUAL_ENERGY_REFLECTANCE};
 
 use crate::{lerp::Lerp, math::fast_exp};
 
@@ -600,8 +595,14 @@ impl DivAssign<f32> for XYZ {
 /// Colours" by Meng et al.
 #[inline(always)]
 fn xyz_to_spectrum_4(xyz: (f32, f32, f32), wavelengths: Float4) -> Float4 {
-    spectrum_xyz_to_p_4(wavelengths, xyz) * Float4::splat(1.0 / EQUAL_ENERGY_REFLECTANCE)
+    use spectral_upsampling as su;
-    // aces_to_spectrum_p4(wavelengths, xyz_to_aces_ap0_e(xyz))
+
+    // su::meng::spectrum_xyz_to_p_4(wavelengths, xyz)
+    //     * Float4::splat(1.0 / su::meng::EQUAL_ENERGY_REFLECTANCE)
+
+    su::jakob::rec2020_to_spectrum_p4(wavelengths, color::xyz_to_rec2020_e(xyz))
+
+    // su::jakob::rec709_to_spectrum_p4(wavelengths, color::xyz_to_rec709_e(xyz))
 }
 
 /// Close analytic approximations of the CIE 1931 XYZ color curves.

src/image.rs

@@ -13,7 +13,9 @@ use std::{
 
 use half::f16;
 
-use crate::color::{xyz_to_rec709_e, XYZ};
+pub use color::{rec709_e_to_xyz, xyz_to_rec709_e};
+
+use crate::color::XYZ;
 
 #[derive(Debug)]
 #[allow(clippy::type_complexity)]

src/parse/psy.rs

@@ -6,14 +6,11 @@ use nom::{combinator::all_consuming, sequence::tuple, IResult};
 
 use kioku::Arena;
 
+use color::rec709_e_to_xyz;
+
 use crate::{
-    camera::Camera,
+    camera::Camera, color::Color, light::WorldLightSource, math::Xform, renderer::Renderer,
-    color::{rec709_e_to_xyz, Color},
+    scene::Scene, scene::World,
-    light::WorldLightSource,
-    math::Xform,
-    renderer::Renderer,
-    scene::Scene,
-    scene::World,
 };
 
 use super::{

src/renderer.rs

@@ -285,8 +285,7 @@ impl<'a> Renderer<'a> {
 
             // Pre-calculate base64 encoding if needed
             let base64_enc = if do_blender_output {
-                use crate::color::xyz_to_rec709_e;
+                Some(img_bucket.rgba_base64(color::xyz_to_rec709_e))
-                Some(img_bucket.rgba_base64(xyz_to_rec709_e))
             } else {
                 None
             };

src/tracer.rs

@@ -1,5 +1,5 @@
 use crate::{
-    color::{rec709_to_xyz, Color},
+    color::Color,
     lerp::lerp_slice,
     math::XformFull,
     ray::{LocalRay, Ray},
@@ -109,7 +109,7 @@ impl<'a> Tracer<'a> {
         match *obj {
             Object::Surface(surface) => {
                 let unassigned_shader = SimpleSurfaceShader::Emit {
-                    color: Color::new_xyz(rec709_to_xyz((1.0, 0.0, 1.0))),
+                    color: Color::new_xyz(color::rec709_to_xyz((1.0, 0.0, 1.0))),
                 };
                 let shader = surface_shader.unwrap_or(&unassigned_shader);
 
@@ -119,7 +119,7 @@ impl<'a> Tracer<'a> {
             Object::SurfaceLight(surface) => {
                 // Lights don't use shaders
                 let bogus_shader = SimpleSurfaceShader::Emit {
-                    color: Color::new_xyz(rec709_to_xyz((1.0, 0.0, 1.0))),
+                    color: Color::new_xyz(color::rec709_to_xyz((1.0, 0.0, 1.0))),
                 };
 
                 surface.intersect_ray(ray, local_ray, space, isect, &bogus_shader);

sub_crates/color/Cargo.toml

@@ -6,6 +6,9 @@ edition = "2018"
 license = "MIT, Apache 2.0"
 build = "build.rs"
 
+[build-dependencies]
+colorbox = "0.3"
+
 [lib]
 name = "color"
 path = "src/lib.rs"

sub_crates/color/build.rs

@@ -1,76 +1,46 @@
 use std::{env, fs::File, io::Write, path::Path};
 
-#[derive(Copy, Clone)]
+use colorbox::{
-struct Chromaticities {
+    chroma::{self, Chromaticities},
-    r: (f64, f64),
+    matrix::{invert, rgb_to_xyz_matrix, xyz_chromatic_adaptation_matrix, AdaptationMethod},
-    g: (f64, f64),
+    matrix_compose,
-    b: (f64, f64),
+};
-    w: (f64, f64),
-}
 
 fn main() {
     let out_dir = env::var("OUT_DIR").unwrap();
 
     // Rec709
     {
-        let chroma = Chromaticities {
-            r: (0.640, 0.330),
-            g: (0.300, 0.600),
-            b: (0.150, 0.060),
-            w: (0.3127, 0.3290),
-        };
-
         let dest_path = Path::new(&out_dir).join("rec709_inc.rs");
         let mut f = File::create(&dest_path).unwrap();
-        write_conversion_functions("rec709", chroma, &mut f);
+        write_conversion_functions("rec709", chroma::REC709, &mut f);
     }
 
     // Rec2020
     {
-        let chroma = Chromaticities {
-            r: (0.708, 0.292),
-            g: (0.170, 0.797),
-            b: (0.131, 0.046),
-            w: (0.3127, 0.3290),
-        };
-
         let dest_path = Path::new(&out_dir).join("rec2020_inc.rs");
         let mut f = File::create(&dest_path).unwrap();
-        write_conversion_functions("rec2020", chroma, &mut f);
+        write_conversion_functions("rec2020", chroma::REC2020, &mut f);
     }
 
     // ACES AP0
     {
-        let chroma = Chromaticities {
-            r: (0.73470, 0.26530),
-            g: (0.00000, 1.00000),
-            b: (0.00010, -0.07700),
-            w: (0.32168, 0.33767),
-        };
-
         let dest_path = Path::new(&out_dir).join("aces_ap0_inc.rs");
         let mut f = File::create(&dest_path).unwrap();
-        write_conversion_functions("aces_ap0", chroma, &mut f);
+        write_conversion_functions("aces_ap0", chroma::ACES_AP0, &mut f);
     }
 
     // ACES AP1
     {
-        let chroma = Chromaticities {
-            r: (0.713, 0.293),
-            g: (0.165, 0.830),
-            b: (0.128, 0.044),
-            w: (0.32168, 0.33767),
-        };
-
         let dest_path = Path::new(&out_dir).join("aces_ap1_inc.rs");
         let mut f = File::create(&dest_path).unwrap();
-        write_conversion_functions("aces_ap1", chroma, &mut f);
+        write_conversion_functions("aces_ap1", chroma::ACES_AP1, &mut f);
     }
 }
 
 /// Generates conversion functions for the given rgb to xyz transform matrix.
 fn write_conversion_functions(space_name: &str, chroma: Chromaticities, f: &mut File) {
-    let to_xyz = rgb_to_xyz(chroma, 1.0);
+    let to_xyz = rgb_to_xyz_matrix(chroma);
 
     f.write_all(
         format!(
@@ -99,7 +69,7 @@ pub fn {}_to_xyz(rgb: (f32, f32, f32)) -> (f32, f32, f32) {{
     )
     .unwrap();
 
-    let inv = inverse(to_xyz);
+    let inv = invert(to_xyz).unwrap();
     f.write_all(
         format!(
             r#"
@@ -127,12 +97,14 @@ pub fn xyz_to_{}(xyz: (f32, f32, f32)) -> (f32, f32, f32) {{
     )
     .unwrap();
 
-    let e_chroma = {
+    let e_to_xyz = matrix_compose!(
-        let mut e_chroma = chroma;
+        rgb_to_xyz_matrix(chroma),
-        e_chroma.w = (1.0 / 3.0, 1.0 / 3.0);
+        xyz_chromatic_adaptation_matrix(
-        e_chroma
+            chroma.w,
-    };
+            (1.0 / 3.0, 1.0 / 3.0),
-    let e_to_xyz = rgb_to_xyz(e_chroma, 1.0);
+            AdaptationMethod::Bradford,
+        ),
+    );
     f.write_all(
         format!(
             r#"
@@ -160,7 +132,7 @@ pub fn {}_e_to_xyz(rgb: (f32, f32, f32)) -> (f32, f32, f32) {{
     )
     .unwrap();
 
-    let inv_e = inverse(e_to_xyz);
+    let inv_e = invert(e_to_xyz).unwrap();
     f.write_all(
         format!(
             r#"
@@ -188,135 +160,3 @@ pub fn xyz_to_{}_e(xyz: (f32, f32, f32)) -> (f32, f32, f32) {{
     )
     .unwrap();
 }
-
-/// Port of the RGBtoXYZ function from the ACES CTL reference implementation.
-/// See lib/IlmCtlMath/CtlColorSpace.cpp in the CTL reference implementation.
-///
-/// This takes the chromaticities of an RGB colorspace and generates a
-/// transform matrix from that space to XYZ.
-///
-/// * `chroma` is the chromaticities.
-/// * `y` is the XYZ "Y" value that should map to RGB (1,1,1)
-fn rgb_to_xyz(chroma: Chromaticities, y: f64) -> [[f64; 3]; 3] {
-    // X and Z values of RGB value (1, 1, 1), or "white"
-    let x = chroma.w.0 * y / chroma.w.1;
-    let z = (1.0 - chroma.w.0 - chroma.w.1) * y / chroma.w.1;
-
-    // Scale factors for matrix rows
-    let d = chroma.r.0 * (chroma.b.1 - chroma.g.1)
-        + chroma.b.0 * (chroma.g.1 - chroma.r.1)
-        + chroma.g.0 * (chroma.r.1 - chroma.b.1);
-
-    let sr = (x * (chroma.b.1 - chroma.g.1)
-        - chroma.g.0 * (y * (chroma.b.1 - 1.0) + chroma.b.1 * (x + z))
-        + chroma.b.0 * (y * (chroma.g.1 - 1.0) + chroma.g.1 * (x + z)))
-        / d;
-
-    let sg = (x * (chroma.r.1 - chroma.b.1)
-        + chroma.r.0 * (y * (chroma.b.1 - 1.0) + chroma.b.1 * (x + z))
-        - chroma.b.0 * (y * (chroma.r.1 - 1.0) + chroma.r.1 * (x + z)))
-        / d;
-
-    let sb = (x * (chroma.g.1 - chroma.r.1)
-        - chroma.r.0 * (y * (chroma.g.1 - 1.0) + chroma.g.1 * (x + z))
-        + chroma.g.0 * (y * (chroma.r.1 - 1.0) + chroma.r.1 * (x + z)))
-        / d;
-
-    // Assemble the matrix
-    let mut mat = [[0.0; 3]; 3];
-
-    mat[0][0] = sr * chroma.r.0;
-    mat[0][1] = sg * chroma.g.0;
-    mat[0][2] = sb * chroma.b.0;
-
-    mat[1][0] = sr * chroma.r.1;
-    mat[1][1] = sg * chroma.g.1;
-    mat[1][2] = sb * chroma.b.1;
-
-    mat[2][0] = sr * (1.0 - chroma.r.0 - chroma.r.1);
-    mat[2][1] = sg * (1.0 - chroma.g.0 - chroma.g.1);
-    mat[2][2] = sb * (1.0 - chroma.b.0 - chroma.b.1);
-
-    mat
-}
-
-/// Calculates the inverse of the given 3x3 matrix.
-///
-/// Ported to Rust from `gjInverse()` in IlmBase's Imath/ImathMatrix.h
-fn inverse(m: [[f64; 3]; 3]) -> [[f64; 3]; 3] {
-    let mut s = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
-    let mut t = m;
-
-    // Forward elimination
-    for i in 0..2 {
-        let mut pivot = i;
-        let mut pivotsize = t[i][i];
-
-        if pivotsize < 0.0 {
-            pivotsize = -pivotsize;
-        }
-
-        for j in (i + 1)..3 {
-            let mut tmp = t[j][i];
-
-            if tmp < 0.0 {
-                tmp = -tmp;
-            }
-
-            if tmp > pivotsize {
-                pivot = j;
-                pivotsize = tmp;
-            }
-        }
-
-        if pivotsize == 0.0 {
-            panic!("Cannot invert singular matrix.");
-        }
-
-        if pivot != i {
-            for j in 0..3 {
-                let mut tmp = t[i][j];
-                t[i][j] = t[pivot][j];
-                t[pivot][j] = tmp;
-
-                tmp = s[i][j];
-                s[i][j] = s[pivot][j];
-                s[pivot][j] = tmp;
-            }
-        }
-
-        for j in (i + 1)..3 {
-            let f = t[j][i] / t[i][i];
-
-            for k in 0..3 {
-                t[j][k] -= f * t[i][k];
-                s[j][k] -= f * s[i][k];
-            }
-        }
-    }
-
-    // Backward substitution
-    for i in (0..3).rev() {
-        let f = t[i][i];
-
-        if t[i][i] == 0.0 {
-            panic!("Cannot invert singular matrix.");
-        }
-
-        for j in 0..3 {
-            t[i][j] /= f;
-            s[i][j] /= f;
-        }
-
-        for j in 0..i {
-            let f = t[j][i];
-
-            for k in 0..3 {
-                t[j][k] -= f * t[i][k];
-                s[j][k] -= f * s[i][k];
-            }
-        }
-    }
-
-    s
-}