Beginnings of light sources. Yay!

src/color/mod.rs

@@ -2,6 +2,7 @@ mod spectra_xyz;
 
 use std::ops::{Add, AddAssign, Mul, MulAssign, Div, DivAssign};
 
+use lerp::Lerp;
 use self::spectra_xyz::{spectrum_xyz_to_p, EQUAL_ENERGY_REFLECTANCE};
 
 // Minimum and maximum wavelength of light we care about, in nanometers
@@ -12,6 +13,26 @@ const WL_RANGE_Q: f32 = WL_RANGE / 4.0;
 
 pub trait Color {
     fn sample_spectrum(&self, wavelength: f32) -> f32;
+
+    fn get_spectral_sample(&self, hero_wavelength: f32) -> SpectralSample {
+        SpectralSample {
+            e: [self.sample_spectrum(nth_wavelength(hero_wavelength, 0)),
+                self.sample_spectrum(nth_wavelength(hero_wavelength, 1)),
+                self.sample_spectrum(nth_wavelength(hero_wavelength, 2)),
+                self.sample_spectrum(nth_wavelength(hero_wavelength, 3))],
+
+            hero_wavelength: hero_wavelength,
+        }
+    }
+}
+
+fn nth_wavelength(hero_wavelength: f32, n: usize) -> f32 {
+    let mut wl = hero_wavelength + (WL_RANGE_Q * n as f32);
+    if wl > WL_MAX {
+        wl - WL_RANGE
+    } else {
+        wl
+    }
 }
 
 
@@ -87,6 +108,12 @@ impl Color for XYZ {
     }
 }
 
+impl Lerp for XYZ {
+    fn lerp(self, other: XYZ, alpha: f32) -> XYZ {
+        (self * (1.0 - alpha)) + (other * alpha)
+    }
+}
+
 impl Add for XYZ {
     type Output = XYZ;
     fn add(self, rhs: XYZ) -> Self::Output {

src/light/mod.rs

@@ -0,0 +1,64 @@
+mod sphere_light;
+
+// pub use sphere_light::{SphereLight};
+
+use math::{Vector, Point};
+use color::SpectralSample;
+
+pub trait LightSource {
+    /// Samples the light source for a given point to be illuminated.
+    ///
+    ///     - arr: The point to be illuminated.
+    ///     - u: Random parameter U.
+    ///     - v: Random parameter V.
+    ///     - wavelength: The wavelength of light to sample at.
+    ///     - time: The time to sample at.
+    ///
+    /// Returns: The light arriving at the point arr, the vector to use for
+    /// shadow testing, and the pdf of the sample.
+    fn sample(&self,
+              arr: Point,
+              u: f32,
+              v: f32,
+              wavelength: f32,
+              time: f32)
+              -> (SpectralSample, Vector, f32);
+
+
+    /// Calculates the pdf of sampling the given
+    /// sample_dir/sample_u/sample_v from the given point arr.  This is used
+    /// primarily to calculate probabilities for multiple importance sampling.
+    ///
+    /// NOTE: this function CAN assume that sample_dir, sample_u, and sample_v
+    /// are a valid sample for the light source (i.e. hits/lies on the light
+    /// source).  No guarantees are made about the correctness of the return
+    /// value if they are not valid.
+    fn sample_pdf(&self,
+                  arr: Point,
+                  sample_dir: Vector,
+                  sample_u: f32,
+                  sample_v: f32,
+                  wavelength: f32,
+                  time: f32)
+                  -> f32;
+
+
+    /// Returns the color emitted in the given direction from the
+    /// given parameters on the light.
+    ///
+    ///     - dir: The direction of the outgoing light.
+    ///     - u: Random parameter U.
+    ///     - v: Random parameter V.
+    ///     - wavelength: The wavelength of light to sample at.
+    ///     - time: The time to sample at.
+    fn outgoing(&self, dir: Vector, u: f32, v: f32, wavelength: f32, time: f32) -> SpectralSample;
+
+
+
+    /// Returns whether the light has a delta distribution.
+    ///
+    /// If a light has no chance of a ray hitting it through random process
+    /// then it is a delta light source.  For example, point light sources,
+    /// lights that only emit in a single direction, etc.
+    fn is_delta(&self) -> bool;
+}

src/light/sphere_light.rs

@@ -0,0 +1,129 @@
+use math::{Vector, Point, coordinate_system_from_vector};
+use bbox::BBox;
+use color::{XYZ, SpectralSample, Color};
+use super::LightSource;
+use lerp::lerp_slice;
+use sampling::{uniform_sample_cone, uniform_sample_cone_pdf, uniform_sample_sphere};
+use std::f64::consts::PI as PI_64;
+use std::f32::consts::PI as PI_32;
+
+pub struct SphereLight {
+    positions: Vec<Point>,
+    radii: Vec<f32>,
+    colors: Vec<XYZ>,
+    bounds_: Vec<BBox>,
+}
+
+impl SphereLight {
+    fn new() -> SphereLight {
+        unimplemented!()
+    }
+}
+
+impl LightSource for SphereLight {
+    fn sample(&self,
+              arr: Point,
+              u: f32,
+              v: f32,
+              wavelength: f32,
+              time: f32)
+              -> (SpectralSample, Vector, f32) {
+        // Calculate time interpolated values
+        let pos = lerp_slice(&self.positions, time);
+        let radius: f64 = lerp_slice(&self.radii, time) as f64;
+        let col = lerp_slice(&self.colors, time);
+        let surface_area_inv: f64 = 1.0 / (4.0 * PI_64 * radius * radius);
+
+
+        // Create a coordinate system from the vector between the
+        // point and the center of the light
+        let (x, y, z): (Vector, Vector, Vector);
+        z = pos - arr;
+        let d2: f64 = z.length2() as f64;  // Distance from center of sphere squared
+        let d = d2.sqrt(); // Distance from center of sphere
+        let (z, x, y) = coordinate_system_from_vector(z);
+        let (x, y, z) = (x.normalized(), y.normalized(), z.normalized());
+
+        // If we're outside the sphere, sample the surface based on
+        // the angle it subtends from the point being lit.
+        if d > radius {
+            // Calculate the portion of the sphere visible from the point
+            let sin_theta_max2: f64 = ((radius * radius) / d2).min(1.0);
+            let cos_theta_max2: f64 = 1.0 - sin_theta_max2;
+            let sin_theta_max: f64 = sin_theta_max2.sqrt();
+            let cos_theta_max: f64 = cos_theta_max2.sqrt();
+
+            // Sample the cone subtended by the sphere and calculate
+            // useful data from that.
+            let sample = uniform_sample_cone(u, v, cos_theta_max).normalized();
+            let cos_theta: f64 = sample[2] as f64;
+            let cos_theta2: f64 = cos_theta * cos_theta;
+            let sin_theta2: f64 = (1.0 - cos_theta2).max(0.0);
+            let sin_theta: f64 = sin_theta2.sqrt();
+
+            // Convert to a point on the sphere.
+            // The technique for this is from "Akalin" on ompf2.com:
+            // http://ompf2.com/viewtopic.php?f=3&t=1914#p4414
+            let D = 1.0 - (d2 * sin_theta * sin_theta / (radius * radius));
+            let cos_a = if D <= 0.0 {
+                sin_theta_max
+            } else {
+                ((d / radius) * sin_theta2) + (cos_theta * D.sqrt())
+            };
+            let sin_a = ((1.0 - (cos_a * cos_a)).max(0.0)).sqrt();
+            let phi = v as f64 * 2.0 * PI_64;
+            let sample = Vector::new((phi.cos() * sin_a * radius) as f32,
+                                     (phi.sin() * sin_a * radius) as f32,
+                                     (d - (cos_a * radius)) as f32);
+
+            // Calculate the final values and return everything.
+            let shadow_vec = (x * sample[0]) + (y * sample[1]) + (z * sample[2]);
+            let pdf = uniform_sample_cone_pdf(cos_theta_max);
+            let spectral_sample = (col * surface_area_inv as f32).get_spectral_sample(wavelength);
+            return (spectral_sample, shadow_vec, pdf as f32);
+        } else {
+            // If we're inside the sphere, there's light from every direction.
+            let shadow_vec = uniform_sample_sphere(u, v);
+            let pdf = 1.0 / (4.0 * PI_64);
+            let spectral_sample = (col * surface_area_inv as f32).get_spectral_sample(wavelength);
+            return (spectral_sample, shadow_vec, pdf as f32);
+        }
+    }
+
+    fn sample_pdf(&self,
+                  arr: Point,
+                  sample_dir: Vector,
+                  sample_u: f32,
+                  sample_v: f32,
+                  wavelength: f32,
+                  time: f32)
+                  -> f32 {
+        let pos = lerp_slice(&self.positions, time);
+        let radius: f64 = lerp_slice(&self.radii, time) as f64;
+
+        let d2: f64 = (pos - arr).length2() as f64;  // Distance from center of sphere squared
+        let d: f64 = d2.sqrt(); // Distance from center of sphere
+
+        if d > radius {
+            // Calculate the portion of the sphere visible from the point
+            let sin_theta_max2: f64 = ((radius * radius) / d2).min(1.0);
+            let cos_theta_max2: f64 = 1.0 - sin_theta_max2;
+            let cos_theta_max: f64 = cos_theta_max2.sqrt();
+
+            return uniform_sample_cone_pdf(cos_theta_max) as f32;
+        } else {
+            return (1.0 / (4.0 * PI_64)) as f32;
+        }
+    }
+
+    fn outgoing(&self, dir: Vector, u: f32, v: f32, wavelength: f32, time: f32) -> SpectralSample {
+        let radius = lerp_slice(&self.radii, time) as f64;
+        let col = lerp_slice(&self.colors, time);
+        let surface_area = 4.0 * PI_64 * radius * radius;
+        (col / surface_area as f32).get_spectral_sample(wavelength)
+    }
+
+    fn is_delta(&self) -> bool {
+        false
+    }
+}

src/main.rs

@@ -20,10 +20,12 @@ mod image;
 mod boundable;
 mod triangle;
 mod surface;
+mod light;
 mod bvh;
 mod scene;
 mod assembly;
 mod halton;
+mod sampling;
 mod color;
 
 use std::mem;

src/math/mod.rs

@@ -40,6 +40,23 @@ pub fn fast_ln(x: f32) -> f32 {
 
 
 
+/// Creates a coordinate system from a single vector.
+pub fn coordinate_system_from_vector(v: Vector) -> (Vector, Vector, Vector) {
+    let v2 = if v[0].abs() > v[1].abs() {
+        let invlen = 1.0 / ((v[0] * v[0]) + (v[2] * v[2])).sqrt();
+        Vector::new(-v[2] * invlen, 0.0, v[0] * invlen)
+    } else {
+        let invlen = 1.0 / ((v[1] * v[1]) + (v[2] * v[2])).sqrt();
+        Vector::new(0.0, v[2] * invlen, -v[1] * invlen)
+    };
+
+    let v3 = cross(v, v2);
+
+    (v, v2, v3)
+}
+
+
+
 /// The logit function, scaled to approximate the probit function.
 ///
 /// We use this as a close approximation to the gaussian inverse CDF,

src/sampling.rs

@@ -0,0 +1,71 @@
+use math::Vector;
+
+use std::f64::consts::PI as PI_64;
+use std::f32::consts::PI as PI_32;
+use std::f32::consts::FRAC_PI_4 as QPI_32;
+
+/// Maps the unit square to the unit circle.
+/// Modifies x and y in place.
+/// NOTE: x and y should be distributed within [-1, 1],
+/// not [0, 1].
+pub fn square_to_circle(x: f32, y: f32) -> (f32, f32) {
+    debug_assert!(x >= -1.0 && x <= 1.0 && y >= -1.0 && y <= 1.0);
+
+    if x == 0.0 && y == 0.0 {
+        return (0.0, 0.0);
+    }
+
+    let (radius, angle) = if x > y.abs() {
+        // Quadrant 1
+        (x, QPI_32 * (y / x))
+    } else if y > x.abs() {
+        // Quadrant 2
+        (y, QPI_32 * (2.0 - (x / y)))
+    } else if x < -(y.abs()) {
+        // Quadrant 3
+        (-x, QPI_32 * (4.0 + (y / x)))
+    } else {
+        // Quadrant 4
+        (-y, QPI_32 * (6.0 - (x / y)))
+    };
+
+    (radius * angle.cos(), radius * angle.sin())
+}
+
+pub fn cosine_sample_hemisphere(u: f32, v: f32) -> Vector {
+    let (u, v) = square_to_circle((u * 2.0) - 1.0, (v * 2.0) - 1.0);
+    let z = (1.0 - ((u * u) + (v * v))).max(0.0).sqrt();
+    return Vector::new(u, v, z);
+}
+
+pub fn uniform_sample_hemisphere(u: f32, v: f32) -> Vector {
+    let z = u;
+    let r = (1.0 - (z * z)).max(0.0).sqrt();
+    let phi = 2.0 * PI_32 * v;
+    let x = r * phi.cos();
+    let y = r * phi.sin();
+    Vector::new(x, y, z)
+}
+
+pub fn uniform_sample_sphere(u: f32, v: f32) -> Vector {
+    let z = 1.0 - (2.0 * u);
+    let r = (1.0 - (z * z)).max(0.0).sqrt();
+    let phi = 2.0 * PI_32 * v;
+    let x = r * phi.cos();
+    let y = r * phi.sin();
+    Vector::new(x, y, z)
+}
+
+pub fn uniform_sample_cone(u: f32, v: f32, cos_theta_max: f64) -> Vector {
+    let cos_theta = (1.0 - u as f64) + (u as f64 * cos_theta_max);
+    let sin_theta = (1.0 - (cos_theta * cos_theta)).sqrt();
+    let phi = v as f64 * 2.0 * PI_64;
+    Vector::new((phi.cos() * sin_theta) as f32,
+                (phi.sin() * sin_theta) as f32,
+                cos_theta as f32)
+}
+
+pub fn uniform_sample_cone_pdf(cos_theta_max: f64) -> f64 {
+    // 1.0 / solid angle
+    1.0 / (2.0 * PI_64 * (1.0 - cos_theta_max))
+}