diff --git a/Cargo.lock b/Cargo.lock index a56882f..aa85c37 100644 --- a/Cargo.lock +++ b/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)", diff --git a/Cargo.toml b/Cargo.toml index 80e86cf..088dc73 100644 --- a/Cargo.toml +++ b/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" diff --git a/src/bbox.rs b/src/bbox.rs index b3e8b11..11995ca 100644 --- a/src/bbox.rs +++ b/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; diff --git a/src/color.rs b/src/color.rs index a15af12..6d2ff37 100644 --- a/src/color.rs +++ b/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 for SpectralSample { } } +//---------------------------------------------------------------- #[derive(Copy, Clone, Debug)] pub struct XYZ { @@ -257,35 +262,7 @@ impl DivAssign 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)] -fn x_1931(wavelength: f32) -> f32 { +pub 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)] -fn y_1931(wavelength: f32) -> f32 { +pub 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)] -fn z_1931(wavelength: f32) -> f32 { +pub 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); - } -} diff --git a/src/image.rs b/src/image.rs index d35dbe9..3a9844c 100644 --- a/src/image.rs +++ b/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)) } diff --git a/src/main.rs b/src/main.rs index 4b6bd8b..ac46a1e 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,3 +1,4 @@ +extern crate color as color_util; extern crate float4; extern crate halton; extern crate math3d; diff --git a/src/parse/psy.rs b/src/parse/psy.rs index 8314ed8..a376092 100644 --- a/src/parse/psy.rs +++ b/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, 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( diff --git a/src/parse/psy_assembly.rs b/src/parse/psy_assembly.rs index 3a88373..021aa23 100644 --- a/src/parse/psy_assembly.rs +++ b/src/parse/psy_assembly.rs @@ -110,7 +110,7 @@ pub fn parse_assembly<'a>(arena: &'a MemArena, tree: &'a DataTree) -> Result { - if let &DataTree::Internal { ident: Some(ident), .. } = child { + if let &DataTree::Internal { ident: Some(_), .. } = child { // TODO //unimplemented!() } else { diff --git a/src/parse/psy_light.rs b/src/parse/psy_light.rs index ad82226..dad9c5b 100644 --- a/src/parse/psy_light.rs +++ b/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)); diff --git a/src/renderer.rs b/src/renderer.rs index 0e55d29..4b4e8af 100644 --- a/src/renderer.rs +++ b/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; - Some(img_bucket.rgba_base64(xyz_to_rec709e)) + use color::xyz_to_rec709_e; + Some(img_bucket.rgba_base64(xyz_to_rec709_e)) } else { None }; diff --git a/src/shading/mod.rs b/src/shading/mod.rs index 05e6b46..db94e3e 100644 --- a/src/shading/mod.rs +++ b/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 } } diff --git a/sub_crates/color/Cargo.toml b/sub_crates/color/Cargo.toml new file mode 100644 index 0000000..a54edf6 --- /dev/null +++ b/sub_crates/color/Cargo.toml @@ -0,0 +1,10 @@ +[package] +name = "color" +version = "0.1.0" +authors = ["Nathan Vegdahl "] +license = "MIT" +build = "build.rs" + +[lib] +name = "color" +path = "src/lib.rs" diff --git a/sub_crates/color/build.rs b/sub_crates/color/build.rs new file mode 100644 index 0000000..7bc54a4 --- /dev/null +++ b/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 +} diff --git a/sub_crates/color/src/lib.rs b/sub_crates/color/src/lib.rs new file mode 100644 index 0000000..8f31cb2 --- /dev/null +++ b/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"));