Split colorspace transform functions out into their own crate.

They are now generated by a build.rs script from nothing but the
colorspace's primaries, which makes it super easy to add more
colorspaces.  So easy that I added three more: ACES AP0, ACES AP1
and Rec.2020.

This lays the foundation for supporting output to different
colorspaces.

Cargo.lock

@@ -50,6 +50,10 @@ dependencies = [
  "vec_map 0.8.0 (registry+https://github.com/rust-lang/crates.io-index)",
 ]
 
+[[package]]
+name = "color"
+version = "0.1.0"
+
 [[package]]
 name = "crossbeam"
 version = "0.2.10"
@@ -155,6 +159,7 @@ version = "0.1.0"
 dependencies = [
  "base64 0.5.2 (registry+https://github.com/rust-lang/crates.io-index)",
  "clap 2.24.2 (registry+https://github.com/rust-lang/crates.io-index)",
+ "color 0.1.0",
  "crossbeam 0.2.10 (registry+https://github.com/rust-lang/crates.io-index)",
  "float4 0.1.0",
  "half 1.0.0 (registry+https://github.com/rust-lang/crates.io-index)",

Cargo.toml

@@ -1,5 +1,6 @@
 [workspace]
 members = [
+    "sub_crates/color",
     "sub_crates/float4",
     "sub_crates/halton",
     "sub_crates/math3d",
@@ -35,6 +36,9 @@ scoped_threadpool = "0.1"
 time = "0.1"
 
 # Local crate dependencies
+[dependencies.color]
+path = "sub_crates/color"
+
 [dependencies.float4]
 path = "sub_crates/float4"
 

src/bbox.rs

@@ -5,7 +5,7 @@ use std::iter::Iterator;
 use std::ops::{BitOr, BitOrAssign};
 
 use lerp::{lerp, lerp_slice, Lerp};
-use math::{Point, Matrix4x4, fast_minf32, fast_maxf32};
+use math::{Point, Matrix4x4, fast_minf32};
 use ray::AccelRay;
 
 

src/color.rs

@@ -6,6 +6,8 @@ use float4::Float4;
 use lerp::Lerp;
 use math::faster_exp;
 
+pub use color_util::{xyz_to_rec709, xyz_to_rec709_e, rec709_to_xyz, rec709_e_to_xyz};
+
 
 // Minimum and maximum wavelength of light we care about, in nanometers
 const WL_MIN: f32 = 380.0;
@@ -39,6 +41,7 @@ fn nth_wavelength(hero_wavelength: f32, n: usize) -> f32 {
     if wl > WL_MAX { wl - WL_RANGE } else { wl }
 }
 
+//----------------------------------------------------------------
 
 #[derive(Copy, Clone, Debug)]
 pub struct SpectralSample {
@@ -55,6 +58,7 @@ impl SpectralSample {
         }
     }
 
+    #[allow(dead_code)]
     pub fn from_value(value: f32, wavelength: f32) -> SpectralSample {
         debug_assert!(wavelength >= WL_MIN && wavelength <= WL_MAX);
         SpectralSample {
@@ -146,6 +150,7 @@ impl DivAssign<f32> for SpectralSample {
     }
 }
 
+//----------------------------------------------------------------
 
 #[derive(Copy, Clone, Debug)]
 pub struct XYZ {
@@ -257,35 +262,7 @@ impl DivAssign<f32> for XYZ {
     }
 }
 
-/// Converts a color in XYZ colorspace to Rec.709 colorspace.
+//----------------------------------------------------------------
-/// Note: this can result in negative values, because the positive Rec.709
-/// colorspace cannot represent all colors in the XYZ colorspace.
-#[allow(dead_code)]
-pub fn xyz_to_rec709(xyz: (f32, f32, f32)) -> (f32, f32, f32) {
-    ((xyz.0 * 3.2404542) + (xyz.1 * -1.5371385) + (xyz.2 * -0.4985314), (xyz.0 * -0.9692660) + (xyz.1 * 1.8760108) + (xyz.2 * 0.0415560), (xyz.0 * 0.0556434) + (xyz.1 * -0.2040259) + (xyz.2 * 1.0572252))
-}
-
-/// Converts a color in Rec.709 colorspace to XYZ colorspace.
-#[allow(dead_code)]
-pub fn rec709_to_xyz(rec: (f32, f32, f32)) -> (f32, f32, f32) {
-    ((rec.0 * 0.4124564) + (rec.1 * 0.3575761) + (rec.2 * 0.1804375), (rec.0 * 0.2126729) + (rec.1 * 0.7151522) + (rec.2 * 0.0721750), (rec.0 * 0.0193339) + (rec.1 * 0.1191920) + (rec.2 * 0.9503041))
-}
-
-/// Converts a color in XYZ colorspace to an adjusted Rec.709 colorspace
-/// with whitepoint E.
-/// Note: this is lossy, as negative resulting values are clamped to zero.
-#[allow(dead_code)]
-pub fn xyz_to_rec709e(xyz: (f32, f32, f32)) -> (f32, f32, f32) {
-    ((xyz.0 * 3.0799600) + (xyz.1 * -1.5371400) + (xyz.2 * -0.5428160), (xyz.0 * -0.9212590) + (xyz.1 * 1.8760100) + (xyz.2 * 0.0452475), (xyz.0 * 0.0528874) + (xyz.1 * -0.2040260) + (xyz.2 * 1.1511400))
-}
-
-/// Converts a color in an adjusted Rec.709 colorspace with whitepoint E to
-/// XYZ colorspace.
-#[allow(dead_code)]
-pub fn rec709e_to_xyz(rec: (f32, f32, f32)) -> (f32, f32, f32) {
-    ((rec.0 * 0.4339499) + (rec.1 * 0.3762098) + (rec.2 * 0.1898403), (rec.0 * 0.2126729) + (rec.1 * 0.7151522) + (rec.2 * 0.0721750), (rec.0 * 0.0177566) + (rec.1 * 0.1094680) + (rec.2 * 0.8727755))
-}
-
 
 /// Samples an CIE 1931 XYZ color at a particular wavelength, according to
 /// the method in the paper "Physically Meaningful Rendering using Tristimulus
@@ -294,165 +271,24 @@ fn xyz_to_spectrum(xyz: (f32, f32, f32), wavelength: f32) -> f32 {
     spectrum_xyz_to_p(wavelength, xyz) * (1.0 / EQUAL_ENERGY_REFLECTANCE)
 }
 
-
 /// Close analytic approximations of the CIE 1931 XYZ color curves.
 /// From the paper "Simple Analytic Approximations to the CIE XYZ Color Matching
 /// Functions" by Wyman et al.
-#[allow(dead_code)]
+pub fn x_1931(wavelength: f32) -> f32 {
-fn x_1931(wavelength: f32) -> f32 {
     let t1 = (wavelength - 442.0) * (if wavelength < 442.0 { 0.0624 } else { 0.0374 });
     let t2 = (wavelength - 599.8) * (if wavelength < 599.8 { 0.0264 } else { 0.0323 });
     let t3 = (wavelength - 501.1) * (if wavelength < 501.1 { 0.0490 } else { 0.0382 });
     (0.362 * faster_exp(-0.5 * t1 * t1)) + (1.056 * faster_exp(-0.5 * t2 * t2)) - (0.065 * faster_exp(-0.5 * t3 * t3))
 }
 
-#[allow(dead_code)]
+pub fn y_1931(wavelength: f32) -> f32 {
-fn y_1931(wavelength: f32) -> f32 {
     let t1 = (wavelength - 568.8) * (if wavelength < 568.8 { 0.0213 } else { 0.0247 });
     let t2 = (wavelength - 530.9) * (if wavelength < 530.9 { 0.0613 } else { 0.0322 });
     (0.821 * faster_exp(-0.5 * t1 * t1)) + (0.286 * faster_exp(-0.5 * t2 * t2))
 }
 
-#[allow(dead_code)]
+pub fn z_1931(wavelength: f32) -> f32 {
-fn z_1931(wavelength: f32) -> f32 {
     let t1 = (wavelength - 437.0) * (if wavelength < 437.0 { 0.0845 } else { 0.0278 });
     let t2 = (wavelength - 459.0) * (if wavelength < 459.0 { 0.0385 } else { 0.0725 });
     (1.217 * faster_exp(-0.5 * t1 * t1)) + (0.681 * faster_exp(-0.5 * t2 * t2))
 }
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    fn abs_diff_tri(a: (f32, f32, f32), b: (f32, f32, f32)) -> (f32, f32, f32) {
-        ((a.0 - b.0).abs(), (a.1 - b.1).abs(), (a.2 - b.2).abs())
-    }
-
-    #[test]
-    fn rec709_xyz_01() {
-        let c1 = (1.0, 1.0, 1.0);
-        let c2 = rec709_to_xyz(c1);
-        let c3 = xyz_to_rec709(c2);
-
-        let diff = abs_diff_tri(c1, c3);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-
-    #[test]
-    fn rec709_xyz_02() {
-        let c1 = (1.0, 1.0, 1.0);
-        let c2 = xyz_to_rec709(c1);
-        let c3 = rec709_to_xyz(c2);
-
-        let diff = abs_diff_tri(c1, c3);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-
-    #[test]
-    fn rec709_xyz_03() {
-        let c1 = (0.9, 0.05, 0.8);
-        let c2 = rec709_to_xyz(c1);
-        let c3 = xyz_to_rec709(c2);
-
-        let diff = abs_diff_tri(c1, c3);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-
-    #[test]
-    fn rec709_xyz_04() {
-        let c1 = (0.9, 0.05, 0.8);
-        let c2 = xyz_to_rec709(c1);
-        let c3 = rec709_to_xyz(c2);
-
-        let diff = abs_diff_tri(c1, c3);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-
-    #[test]
-    fn rec709e_xyz_01() {
-        let c1 = (1.0, 1.0, 1.0);
-        let c2 = rec709e_to_xyz(c1);
-        let c3 = xyz_to_rec709e(c2);
-
-        let diff = abs_diff_tri(c1, c3);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-
-    #[test]
-    fn rec709e_xyz_02() {
-        let c1 = (1.0, 1.0, 1.0);
-        let c2 = xyz_to_rec709e(c1);
-        let c3 = rec709e_to_xyz(c2);
-
-        let diff = abs_diff_tri(c1, c3);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-
-    #[test]
-    fn rec709e_xyz_03() {
-        let c1 = (1.0, 1.0, 1.0);
-        let c2 = rec709e_to_xyz(c1);
-
-        let diff = abs_diff_tri(c1, c2);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-
-    #[test]
-    fn rec709e_xyz_04() {
-        let c1 = (1.0, 1.0, 1.0);
-        let c2 = xyz_to_rec709e(c1);
-
-        let diff = abs_diff_tri(c1, c2);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-
-    #[test]
-    fn rec709e_xyz_05() {
-        let c1 = (0.9, 0.05, 0.8);
-        let c2 = rec709e_to_xyz(c1);
-        let c3 = xyz_to_rec709e(c2);
-
-        let diff = abs_diff_tri(c1, c3);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-
-    #[test]
-    fn rec709e_xyz_06() {
-        let c1 = (0.9, 0.05, 0.8);
-        let c2 = xyz_to_rec709e(c1);
-        let c3 = rec709e_to_xyz(c2);
-
-        let diff = abs_diff_tri(c1, c3);
-
-        assert!(diff.0 < 0.00001);
-        assert!(diff.1 < 0.00001);
-        assert!(diff.2 < 0.00001);
-    }
-}

src/image.rs

@@ -14,7 +14,7 @@ use half::f16;
 use png_encode_mini;
 use openexr;
 
-use color::{XYZ, xyz_to_rec709e};
+use color::{XYZ, xyz_to_rec709_e};
 
 
 #[derive(Debug)]
@@ -164,7 +164,7 @@ impl Image {
         // Convert pixels
         for y in 0..self.res.1 {
             for x in 0..self.res.0 {
-                let (r, g, b) = xyz_to_rec709e(self.get(x, y).to_tuple());
+                let (r, g, b) = xyz_to_rec709_e(self.get(x, y).to_tuple());
                 image.push((f16::from_f32(r), f16::from_f32(g), f16::from_f32(b)));
             }
         }
@@ -280,7 +280,7 @@ fn srgb_inv_gamma(n: f32) -> f32 {
 }
 
 fn xyz_to_srgbe(xyz: (f32, f32, f32)) -> (f32, f32, f32) {
-    let rgb = xyz_to_rec709e(xyz);
+    let rgb = xyz_to_rec709_e(xyz);
     (srgb_gamma(rgb.0), srgb_gamma(rgb.1), srgb_gamma(rgb.2))
 }
 

src/main.rs

@@ -1,3 +1,4 @@
+extern crate color as color_util;
 extern crate float4;
 extern crate halton;
 extern crate math3d;

src/parse/psy.rs

@@ -8,7 +8,7 @@ use nom::IResult;
 use mem_arena::MemArena;
 
 use camera::Camera;
-use color::{XYZ, rec709e_to_xyz};
+use color::{XYZ, rec709_e_to_xyz};
 use light::WorldLightSource;
 use math::Matrix4x4;
 use renderer::Renderer;
@@ -533,7 +533,7 @@ fn parse_world<'a>(arena: &'a MemArena, tree: &'a DataTree) -> Result<World<'a>,
                     if let IResult::Done(_, color) = closure!(tuple!(ws_f32, ws_f32, ws_f32))(contents.trim().as_bytes()) {
                         // TODO: proper color space management, not just assuming
                         // rec.709.
-                        background_color = XYZ::from_tuple(rec709e_to_xyz(color));
+                        background_color = XYZ::from_tuple(rec709_e_to_xyz(color));
                     } else {
                         return Err(
                             PsyParseError::IncorrectLeafData(

src/parse/psy_assembly.rs

@@ -110,7 +110,7 @@ pub fn parse_assembly<'a>(arena: &'a MemArena, tree: &'a DataTree) -> Result<Ass
 
                 // Surface shader
                 "SurfaceShader" => {
-                    if let &DataTree::Internal { ident: Some(ident), .. } = child {
+                    if let &DataTree::Internal { ident: Some(_), .. } = child {
                         // TODO
                         //unimplemented!()
                     } else {

src/parse/psy_light.rs

@@ -7,7 +7,7 @@ use nom::IResult;
 use mem_arena::MemArena;
 
 use math::Vector;
-use color::{XYZ, rec709e_to_xyz};
+use color::{XYZ, rec709_e_to_xyz};
 use light::{DistantDiskLight, SphereLight, RectangleLight};
 
 use super::basics::ws_f32;
@@ -62,7 +62,7 @@ pub fn parse_distant_disk_light<'a>(arena: &'a MemArena, tree: &'a DataTree) ->
                         // TODO: handle color space conversions properly.
                         // Probably will need a special color type with its
                         // own parser...?
-                        colors.push(XYZ::from_tuple(rec709e_to_xyz(color)));
+                        colors.push(XYZ::from_tuple(rec709_e_to_xyz(color)));
                     } else {
                         // Found color, but its contents is not in the right format
                         return Err(PsyParseError::UnknownError(byte_offset));
@@ -112,7 +112,7 @@ pub fn parse_sphere_light<'a>(arena: &'a MemArena, tree: &'a DataTree) -> Result
                         // TODO: handle color space conversions properly.
                         // Probably will need a special color type with its
                         // own parser...?
-                        colors.push(XYZ::from_tuple(rec709e_to_xyz(color)));
+                        colors.push(XYZ::from_tuple(rec709_e_to_xyz(color)));
                     } else {
                         // Found color, but its contents is not in the right format
                         return Err(PsyParseError::UnknownError(byte_offset));
@@ -161,7 +161,7 @@ pub fn parse_rectangle_light<'a>(arena: &'a MemArena, tree: &'a DataTree) -> Res
                         // TODO: handle color space conversions properly.
                         // Probably will need a special color type with its
                         // own parser...?
-                        colors.push(XYZ::from_tuple(rec709e_to_xyz(color)));
+                        colors.push(XYZ::from_tuple(rec709_e_to_xyz(color)));
                     } else {
                         // Found color, but its contents is not in the right format
                         return Err(PsyParseError::UnknownError(byte_offset));

src/renderer.rs

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

src/shading/mod.rs

@@ -36,6 +36,7 @@ impl SimpleSurfaceShader {
 
 impl SurfaceShader for SimpleSurfaceShader {
     fn shade(&self, data: &SurfaceIntersectionData) -> SurfaceClosureUnion {
+        let _ = data; // Silence "unused" compiler warning
         self.closure
     }
 }

sub_crates/color/Cargo.toml

@@ -0,0 +1,10 @@
+[package]
+name = "color"
+version = "0.1.0"
+authors = ["Nathan Vegdahl <cessen@cessen.com>"]
+license = "MIT"
+build = "build.rs"
+
+[lib]
+name = "color"
+path = "src/lib.rs"

sub_crates/color/build.rs

@@ -0,0 +1,343 @@
+use std::env;
+use std::fs::File;
+use std::io::Write;
+use std::path::Path;
+
+
+#[derive(Copy, Clone)]
+struct Chromaticities {
+    r: (f64, f64),
+    g: (f64, f64),
+    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 to_xyz = rgb_to_xyz(chroma, 1.0);
+        let dest_path = Path::new(&out_dir).join("rec709_inc.rs");
+        let mut f = File::create(&dest_path).unwrap();
+        write_conversion_functions("rec709", to_xyz, &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 to_xyz = rgb_to_xyz(chroma, 1.0);
+        let dest_path = Path::new(&out_dir).join("rec2020_inc.rs");
+        let mut f = File::create(&dest_path).unwrap();
+        write_conversion_functions("rec2020", to_xyz, &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 to_xyz = rgb_to_xyz(chroma, 1.0);
+        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", to_xyz, &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 to_xyz = rgb_to_xyz(chroma, 1.0);
+        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", to_xyz, &mut f);
+    }
+}
+
+
+/// Generates conversion functions for the given rgb to xyz transform matrix.
+fn write_conversion_functions(space_name: &str, to_xyz: [[f64; 3]; 3], f: &mut File) {
+
+    f.write_all(
+            format!(
+                r#"
+#[inline]
+pub fn {}_to_xyz(rgb: (f32, f32, f32)) -> (f32, f32, f32) {{
+    (
+        (rgb.0 * {:.10}) + (rgb.1 * {:.10}) + (rgb.2 * {:.10}),
+        (rgb.0 * {:.10}) + (rgb.1 * {:.10}) + (rgb.2 * {:.10}),
+        (rgb.0 * {:.10}) + (rgb.1 * {:.10}) + (rgb.2 * {:.10}),
+    )
+}}
+        "#,
+                space_name,
+                to_xyz[0][0],
+                to_xyz[0][1],
+                to_xyz[0][2],
+                to_xyz[1][0],
+                to_xyz[1][1],
+                to_xyz[1][2],
+                to_xyz[2][0],
+                to_xyz[2][1],
+                to_xyz[2][2]
+            )
+                    .as_bytes()
+        )
+        .unwrap();
+
+    let inv = inverse(to_xyz);
+    f.write_all(
+            format!(
+                r#"
+#[inline]
+pub fn xyz_to_{}(xyz: (f32, f32, f32)) -> (f32, f32, f32) {{
+    (
+        (xyz.0 * {:.10}) + (xyz.1 * {:.10}) + (xyz.2 * {:.10}),
+        (xyz.0 * {:.10}) + (xyz.1 * {:.10}) + (xyz.2 * {:.10}),
+        (xyz.0 * {:.10}) + (xyz.1 * {:.10}) + (xyz.2 * {:.10}),
+    )
+}}
+        "#,
+                space_name,
+                inv[0][0],
+                inv[0][1],
+                inv[0][2],
+                inv[1][0],
+                inv[1][1],
+                inv[1][2],
+                inv[2][0],
+                inv[2][1],
+                inv[2][2]
+            )
+                    .as_bytes()
+        )
+        .unwrap();
+
+    let e_to_xyz = adapt_to_e(to_xyz, 1.0);
+    f.write_all(
+            format!(
+                r#"
+#[inline]
+pub fn {}_e_to_xyz(rgb: (f32, f32, f32)) -> (f32, f32, f32) {{
+    (
+        (rgb.0 * {:.10}) + (rgb.1 * {:.10}) + (rgb.2 * {:.10}),
+        (rgb.0 * {:.10}) + (rgb.1 * {:.10}) + (rgb.2 * {:.10}),
+        (rgb.0 * {:.10}) + (rgb.1 * {:.10}) + (rgb.2 * {:.10}),
+    )
+}}
+        "#,
+                space_name,
+                e_to_xyz[0][0],
+                e_to_xyz[0][1],
+                e_to_xyz[0][2],
+                e_to_xyz[1][0],
+                e_to_xyz[1][1],
+                e_to_xyz[1][2],
+                e_to_xyz[2][0],
+                e_to_xyz[2][1],
+                e_to_xyz[2][2]
+            )
+                    .as_bytes()
+        )
+        .unwrap();
+
+    let inv_e = inverse(e_to_xyz);
+    f.write_all(
+            format!(
+                r#"
+#[inline]
+pub fn xyz_to_{}_e(xyz: (f32, f32, f32)) -> (f32, f32, f32) {{
+    (
+        (xyz.0 * {:.10}) + (xyz.1 * {:.10}) + (xyz.2 * {:.10}),
+        (xyz.0 * {:.10}) + (xyz.1 * {:.10}) + (xyz.2 * {:.10}),
+        (xyz.0 * {:.10}) + (xyz.1 * {:.10}) + (xyz.2 * {:.10}),
+    )
+}}
+        "#,
+                space_name,
+                inv_e[0][0],
+                inv_e[0][1],
+                inv_e[0][2],
+                inv_e[1][0],
+                inv_e[1][1],
+                inv_e[1][2],
+                inv_e[2][0],
+                inv_e[2][1],
+                inv_e[2][2]
+            )
+                    .as_bytes()
+        )
+        .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
+}
+
+
+/// Chromatically adapts a matrix from `rgb_to_xyz` to a whitepoint of E.
+///
+/// In other words, makes it so that RGB (1,1,1) maps to XYZ (1,1,1).
+fn adapt_to_e(mat: [[f64; 3]; 3], y: f64) -> [[f64; 3]; 3] {
+    let r_fac = y / (mat[0][0] + mat[0][1] + mat[0][2]);
+    let g_fac = y / (mat[1][0] + mat[1][1] + mat[1][2]);
+    let b_fac = y / (mat[2][0] + mat[2][1] + mat[2][2]);
+
+    let mut mat2 = [[0.0; 3]; 3];
+
+    mat2[0][0] = mat[0][0] * r_fac;
+    mat2[0][1] = mat[0][1] * r_fac;
+    mat2[0][2] = mat[0][2] * r_fac;
+
+    mat2[1][0] = mat[1][0] * g_fac;
+    mat2[1][1] = mat[1][1] * g_fac;
+    mat2[1][2] = mat[1][2] * g_fac;
+
+    mat2[2][0] = mat[2][0] * b_fac;
+    mat2[2][1] = mat[2][1] * b_fac;
+    mat2[2][2] = mat[2][2] * b_fac;
+
+    mat2
+}
+
+
+/// 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
+}

sub_crates/color/src/lib.rs

@@ -0,0 +1,15 @@
+#[allow(non_camel_case_types)]
+#[derive(Copy, Clone)]
+pub enum Space {
+    XYZ,
+    ACES_AP0,
+    ACES_AP1,
+    Rec709,
+    Rec2020,
+}
+
+// Generated conversion functions between XYZ and various RGB colorspaces
+include!(concat!(env!("OUT_DIR"), "/rec709_inc.rs"));
+include!(concat!(env!("OUT_DIR"), "/rec2020_inc.rs"));
+include!(concat!(env!("OUT_DIR"), "/aces_ap0_inc.rs"));
+include!(concat!(env!("OUT_DIR"), "/aces_ap1_inc.rs"));