From 1cbf20e67faece76c45196d89b8b9594fa216663 Mon Sep 17 00:00:00 2001 From: Nathan Vegdahl Date: Sun, 19 Jun 2016 14:32:13 -0700 Subject: [PATCH] Beginnings of light sources. Yay! --- src/color/mod.rs | 27 ++++++++ src/light/mod.rs | 64 +++++++++++++++++++ src/light/sphere_light.rs | 129 ++++++++++++++++++++++++++++++++++++++ src/main.rs | 2 + src/math/mod.rs | 17 +++++ src/sampling.rs | 71 +++++++++++++++++++++ 6 files changed, 310 insertions(+) create mode 100644 src/light/mod.rs create mode 100644 src/light/sphere_light.rs create mode 100644 src/sampling.rs diff --git a/src/color/mod.rs b/src/color/mod.rs index df1763c..5451324 100644 --- a/src/color/mod.rs +++ b/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 { diff --git a/src/light/mod.rs b/src/light/mod.rs new file mode 100644 index 0000000..c88c38c --- /dev/null +++ b/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; +} diff --git a/src/light/sphere_light.rs b/src/light/sphere_light.rs new file mode 100644 index 0000000..7ea6343 --- /dev/null +++ b/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, + radii: Vec, + colors: Vec, + bounds_: Vec, +} + +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 + } +} diff --git a/src/main.rs b/src/main.rs index 5c692d7..ad00634 100644 --- a/src/main.rs +++ b/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; diff --git a/src/math/mod.rs b/src/math/mod.rs index 6808924..8ec35a5 100644 --- a/src/math/mod.rs +++ b/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, diff --git a/src/sampling.rs b/src/sampling.rs new file mode 100644 index 0000000..dda2157 --- /dev/null +++ b/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)) +}