Converted Image over to store renders in XYZ colorspace.

src/color/mod.rs

@@ -1,5 +1,7 @@
 mod spectra_xyz;
 
+use std::ops::{Add, AddAssign, Mul, MulAssign, Div, DivAssign};
+
 use self::spectra_xyz::{spectrum_xyz_to_p, EQUAL_ENERGY_REFLECTANCE};
 
 // Minimum and maximum wavelength of light we care about, in nanometers
@@ -42,6 +44,7 @@ impl SpectralSample {
         self.e[3] *= color.sample_spectrum(self.wl_n(3));
     }
 
+    /// Returns the nth wavelength
     fn wl_n(&self, n: usize) -> f32 {
         let mut wl = self.hero_wavelength + (WL_RANGE_Q * n as f32);
         if wl > WL_MAX {
@@ -61,17 +64,21 @@ pub struct XYZ {
 }
 
 impl XYZ {
-    fn new(x: f32, y: f32, z: f32) -> XYZ {
+    pub fn new(x: f32, y: f32, z: f32) -> XYZ {
         XYZ { x: x, y: y, z: z }
     }
 
-    fn from_wavelength(wavelength: f32, intensity: f32) -> XYZ {
+    pub fn from_wavelength(wavelength: f32, intensity: f32) -> XYZ {
         XYZ {
             x: x_1931(wavelength) * intensity,
             y: y_1931(wavelength) * intensity,
             z: z_1931(wavelength) * intensity,
         }
     }
+
+    pub fn to_tuple(&self) -> (f32, f32, f32) {
+        (self.x, self.y, self.z)
+    }
 }
 
 impl Color for XYZ {
@@ -80,6 +87,100 @@ impl Color for XYZ {
     }
 }
 
+impl Add for XYZ {
+    type Output = XYZ;
+    fn add(self, rhs: XYZ) -> Self::Output {
+        XYZ {
+            x: self.x + rhs.x,
+            y: self.y + rhs.y,
+            z: self.z + rhs.z,
+        }
+    }
+}
+
+impl AddAssign for XYZ {
+    fn add_assign(&mut self, rhs: XYZ) {
+        self.x += rhs.x;
+        self.y += rhs.y;
+        self.z += rhs.z;
+    }
+}
+
+impl Mul<f32> for XYZ {
+    type Output = XYZ;
+    fn mul(self, rhs: f32) -> Self::Output {
+        XYZ {
+            x: self.x * rhs,
+            y: self.y * rhs,
+            z: self.z * rhs,
+        }
+    }
+}
+
+impl MulAssign<f32> for XYZ {
+    fn mul_assign(&mut self, rhs: f32) {
+        self.x *= rhs;
+        self.y *= rhs;
+        self.z *= rhs;
+    }
+}
+
+impl Div<f32> for XYZ {
+    type Output = XYZ;
+    fn div(self, rhs: f32) -> Self::Output {
+        XYZ {
+            x: self.x / rhs,
+            y: self.y / rhs,
+            z: self.z / rhs,
+        }
+    }
+}
+
+impl DivAssign<f32> for XYZ {
+    fn div_assign(&mut self, rhs: f32) {
+        self.x /= rhs;
+        self.y /= rhs;
+        self.z /= rhs;
+    }
+}
+
+/// 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
@@ -92,6 +193,7 @@ fn xyz_to_spectrum(xyz: (f32, f32, f32), wavelength: f32) -> f32 {
 /// 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 {
     let t1 = (wavelength - 442.0) *
              (if wavelength < 442.0 {
@@ -115,6 +217,7 @@ fn x_1931(wavelength: f32) -> f32 {
     (0.065 * (-0.5 * t3 * t3).exp())
 }
 
+#[allow(dead_code)]
 fn y_1931(wavelength: f32) -> f32 {
     let t1 = (wavelength - 568.8) *
              (if wavelength < 568.8 {
@@ -131,6 +234,7 @@ fn y_1931(wavelength: f32) -> f32 {
     (0.821 * (-0.5 * t1 * t1).exp()) + (0.286 * (-0.5 * t2 * t2).exp())
 }
 
+#[allow(dead_code)]
 fn z_1931(wavelength: f32) -> f32 {
     let t1 = (wavelength - 437.0) *
              (if wavelength < 437.0 {
@@ -146,3 +250,140 @@ fn z_1931(wavelength: f32) -> f32 {
     });
     (1.217 * (-0.5 * t1 * t1).exp()) + (0.681 * (-0.5 * t2 * t2).exp())
 }
+
+#[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

@@ -5,16 +5,18 @@ use std::io::Write;
 use std::path::Path;
 use std::fs::File;
 
+use color::{XYZ, xyz_to_rec709e};
+
 #[derive(Debug, Clone)]
 pub struct Image {
-    data: Vec<(f32, f32, f32)>,
+    data: Vec<XYZ>,
     res: (usize, usize),
 }
 
 impl Image {
     pub fn new(width: usize, height: usize) -> Image {
         Image {
-            data: vec![(0.0, 0.0, 0.0); width * height],
+            data: vec![XYZ::new(0.0, 0.0, 0.0); width * height],
             res: (width, height),
         }
     }
@@ -27,14 +29,14 @@ impl Image {
         self.res.1
     }
 
-    pub fn get(&self, x: usize, y: usize) -> (f32, f32, f32) {
+    pub fn get(&self, x: usize, y: usize) -> XYZ {
         assert!(x < self.res.0);
         assert!(y < self.res.1);
 
         self.data[self.res.0 * y + x]
     }
 
-    pub fn set(&mut self, x: usize, y: usize, value: (f32, f32, f32)) {
+    pub fn set(&mut self, x: usize, y: usize, value: XYZ) {
         assert!(x < self.res.0);
         assert!(y < self.res.1);
 
@@ -51,10 +53,7 @@ impl Image {
         // Write pixels
         for y in 0..self.res.1 {
             for x in 0..self.res.0 {
-                let (r, g, b) = self.get(x, y);
+                let (r, g, b) = quantize_tri_255(xyz_to_srgbe(self.get(x, y).to_tuple()));
-                let r = quantize_255(srgb_gamma(r));
-                let g = quantize_255(srgb_gamma(g));
-                let b = quantize_255(srgb_gamma(b));
                 try!(write!(f, "{} {} {} ", r, g, b));
             }
             try!(write!(f, "\n"));
@@ -74,10 +73,7 @@ impl Image {
         // Write pixels
         for y in 0..self.res.1 {
             for x in 0..self.res.0 {
-                let (r, g, b) = self.get(x, y);
+                let (r, g, b) = quantize_tri_255(xyz_to_srgbe(self.get(x, y).to_tuple()));
-                let r = quantize_255(srgb_gamma(r));
-                let g = quantize_255(srgb_gamma(g));
-                let b = quantize_255(srgb_gamma(b));
                 let d = [r, g, b];
                 try!(f.write_all(&d));
             }
@@ -104,7 +100,16 @@ fn srgb_inv_gamma(n: f32) -> f32 {
     }
 }
 
-fn quantize_255(n: f32) -> u8 {
+fn xyz_to_srgbe(xyz: (f32, f32, f32)) -> (f32, f32, f32) {
+    let rgb = xyz_to_rec709e(xyz);
+    (srgb_gamma(rgb.0), srgb_gamma(rgb.1), srgb_gamma(rgb.2))
+}
+
+fn quantize_tri_255(tri: (f32, f32, f32)) -> (u8, u8, u8) {
+    fn quantize(n: f32) -> u8 {
         let n = 1.0f32.min(0.0f32.max(n)) * 255.0;
         n as u8
     }
+
+    (quantize(tri.0), quantize(tri.1), quantize(tri.2))
+}

src/renderer.rs

@@ -14,6 +14,7 @@ use math::fast_logit;
 use image::Image;
 use surface;
 use scene::Scene;
+use color::{XYZ, rec709e_to_xyz};
 
 #[derive(Debug)]
 pub struct Renderer {
@@ -121,15 +122,12 @@ impl Renderer {
                                                                         nor: _,
                                                                         local_space: _,
                                                                         uv } = isect {
-
+                                let rgbcol = (uv.0, uv.1, (1.0 - uv.0 - uv.1).max(0.0));
-                                col.0 += uv.0 / self.spp as f32;
+                                let xyzcol = rec709e_to_xyz(rgbcol);
-                                col.1 += uv.1 / self.spp as f32;
+                                col += XYZ::new(xyzcol.0, xyzcol.1, xyzcol.2) / self.spp as f32;
-                                col.2 += (1.0 - uv.0 - uv.1).max(0.0) / self.spp as f32;
 
                             } else {
-                                col.0 += 0.02 / self.spp as f32;
+                                col += XYZ::new(0.02, 0.02, 0.02) / self.spp as f32;
-                                col.1 += 0.02 / self.spp as f32;
-                                col.2 += 0.02 / self.spp as f32;
                             }
                             img.set(co.0 as usize, co.1 as usize, col);
                         }

todo.txt

@@ -1,3 +1,3 @@
-//- Basic scene parsing with triangle meshes.
+- Move to spectral rendering.
-//- Implement bucketed rendering.
+- Implement basic direct lighting w/ spherical lights.
 - Unit tests for scene parsing.