From 52acee33af929e110fd64bd2612ba0871c982fa8 Mon Sep 17 00:00:00 2001 From: Nathan Vegdahl Date: Wed, 29 Jun 2016 14:30:42 -0700 Subject: [PATCH] Implemented rectangular area lights. Also added Cornell Box test scene. --- .gitignore | 3 +- example_scenes/cornell_box.psy | 114 ++++++++++++++++++++++++ src/lerp.rs | 6 ++ src/light/mod.rs | 2 + src/light/rectangle_light.rs | 157 +++++++++++++++++++++++++++++++++ src/math/normal.rs | 11 ++- src/math/point.rs | 4 + src/math/vector.rs | 11 ++- src/parse/psy_assembly.rs | 15 +++- src/parse/psy_light.rs | 46 +++++++++- src/renderer.rs | 8 +- src/sampling.rs | 95 +++++++++++++++++++- 12 files changed, 464 insertions(+), 8 deletions(-) create mode 100644 example_scenes/cornell_box.psy create mode 100644 src/light/rectangle_light.rs diff --git a/.gitignore b/.gitignore index e04fc8d..1734597 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,5 @@ target *.rs.bk .zedstate -test_renders \ No newline at end of file +test_renders +perf.data* diff --git a/example_scenes/cornell_box.psy b/example_scenes/cornell_box.psy new file mode 100644 index 0000000..54cb937 --- /dev/null +++ b/example_scenes/cornell_box.psy @@ -0,0 +1,114 @@ +Scene $Scene_fr1 { + Output { + Path ["test_renders/cornell_box.ppm"] + } + RenderSettings { + Resolution [512 512] + SamplesPerPixel [16] + Seed [1] + } + Camera { + Fov [39.449188] + FocalDistance [10.620000] + ApertureRadius [0.000000] + Transform [1.000000 -0.000000 0.000000 0.000000 -0.000000 0.000000 1.000000 0.000000 0.000000 1.000000 -0.000000 0.000000 -2.779998 -8.000000 2.730010 1.000000] + } + World { + BackgroundShader { + Type [Color] + Color [0.000000 0.000000 0.000000] + } + } + Assembly { + SurfaceShader $Green { + Type [Lambert] + Color [0.117000 0.412500 0.115000] + } + SurfaceShader $Red { + Type [Lambert] + Color [0.611000 0.055500 0.062000] + } + SurfaceShader $White { + Type [Lambert] + Color [0.729500 0.735500 0.729000] + } + RectangleLight $__Area { + Color [84.300003 53.800003 18.500000] + Dimensions [1.350000 1.100000] + } + Instance { + Data [$__Area] + Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 2.779475 -2.794788 -5.498045 1.000000] + } + MeshSurface $__Plane.010_ { + Vertices [-2.649998 2.959996 3.299997 -4.229996 2.469997 3.299997 -3.139998 4.559995 3.299997 -4.719996 4.059995 3.299997 -4.719996 4.059996 0.000000 -3.139998 4.559995 0.000000 -4.229996 2.469997 0.000000 -2.649998 2.959997 0.000000 ] + FaceVertCounts [4 4 4 4 4 ] + FaceVertIndices [0 1 3 2 1 0 7 6 3 1 6 4 2 3 4 5 0 2 5 7 ] + } + Instance { + Data [$__Plane.010_] + SurfaceShaderBind [$White] + Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000] + } + MeshSurface $__Plane.008_ { + Vertices [-1.299999 0.649999 1.649998 -0.820000 2.249998 1.649999 -2.899997 1.139998 1.649999 -2.399998 2.719997 1.649999 -1.299999 0.649999 0.000000 -0.820000 2.249998 0.000000 -2.899997 1.139998 0.000000 -2.399998 2.719997 0.000000 ] + FaceVertCounts [4 4 4 4 4 ] + FaceVertIndices [0 2 3 1 3 2 6 7 1 3 7 5 0 1 5 4 2 0 4 6 ] + } + Instance { + Data [$__Plane.008_] + SurfaceShaderBind [$White] + Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000] + } + MeshSurface $__Plane.006_ { + Vertices [-5.495996 5.591994 0.000000 -5.527995 -0.000001 -0.000000 -5.559996 5.591993 5.487995 -5.559995 -0.000001 5.487995 ] + FaceVertCounts [4 ] + FaceVertIndices [0 1 3 2 ] + } + Instance { + Data [$__Plane.006_] + SurfaceShaderBind [$Red] + Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000] + } + MeshSurface $__Plane.004_ { + Vertices [-0.000001 5.591995 0.000000 0.000000 0.000000 0.000000 -0.000001 5.591994 5.487995 0.000000 -0.000000 5.487995 ] + FaceVertCounts [4 ] + FaceVertIndices [1 0 2 3 ] + } + Instance { + Data [$__Plane.004_] + SurfaceShaderBind [$Green] + Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000] + } + MeshSurface $__Plane.002_ { + Vertices [-5.495996 5.591994 0.000000 -0.000001 5.591995 0.000000 -5.559996 5.591993 5.487995 -0.000001 5.591994 5.487995 ] + FaceVertCounts [4 ] + FaceVertIndices [0 1 3 2 ] + } + Instance { + Data [$__Plane.002_] + SurfaceShaderBind [$White] + Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000] + } + MeshSurface $__Plane.001_ { + Vertices [-5.559996 5.591993 5.487995 -0.000001 5.591994 5.487995 -5.559995 -0.000001 5.487995 0.000000 -0.000000 5.487995 -3.429997 3.319996 5.487995 -2.129998 3.319996 5.487995 -3.429997 2.269997 5.487995 -2.129998 2.269997 5.487995 ] + FaceVertCounts [4 4 4 4 ] + FaceVertIndices [1 5 4 0 0 4 6 2 2 6 7 3 7 5 1 3 ] + } + Instance { + Data [$__Plane.001_] + SurfaceShaderBind [$White] + Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000] + } + MeshSurface $__Plane_ { + Vertices [-5.495996 5.591994 0.000000 -0.000001 5.591995 0.000000 -5.527995 -0.000001 -0.000000 0.000000 0.000000 0.000000 ] + FaceVertCounts [4 ] + FaceVertIndices [0 1 3 2 ] + } + Instance { + Data [$__Plane_] + SurfaceShaderBind [$White] + Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000] + } + } +} diff --git a/src/lerp.rs b/src/lerp.rs index a939b66..6793ee9 100644 --- a/src/lerp.rs +++ b/src/lerp.rs @@ -67,6 +67,12 @@ impl Lerp for f64 { } } +impl Lerp for (T, T) { + fn lerp(self, other: (T, T), alpha: f32) -> (T, T) { + (self.0.lerp(other.0, alpha), self.1.lerp(other.1, alpha)) + } +} + #[cfg(test)] mod tests { diff --git a/src/light/mod.rs b/src/light/mod.rs index 465b603..0011ef4 100644 --- a/src/light/mod.rs +++ b/src/light/mod.rs @@ -1,8 +1,10 @@ mod sphere_light; +mod rectangle_light; use std::fmt::Debug; pub use self::sphere_light::SphereLight; +pub use self::rectangle_light::RectangleLight; use math::{Vector, Point, Matrix4x4}; use color::SpectralSample; diff --git a/src/light/rectangle_light.rs b/src/light/rectangle_light.rs new file mode 100644 index 0000000..3242ed2 --- /dev/null +++ b/src/light/rectangle_light.rs @@ -0,0 +1,157 @@ +use math::{Vector, Point, Matrix4x4, cross}; +use bbox::BBox; +use boundable::Boundable; +use color::{XYZ, SpectralSample, Color}; +use super::LightSource; +use lerp::lerp_slice; +use sampling::{spherical_triangle_solid_angle, uniform_sample_spherical_triangle}; +use std::f64::consts::PI as PI_64; +use std::f32::consts::PI as PI_32; + +#[derive(Debug)] +pub struct RectangleLight { + dimensions: Vec<(f32, f32)>, + colors: Vec, + bounds_: Vec, +} + +impl RectangleLight { + pub fn new(dimensions: Vec<(f32, f32)>, colors: Vec) -> RectangleLight { + let bbs = dimensions.iter() + .map(|d| { + BBox { + min: Point::new(d.0 * -0.5, d.1 * -0.5, 0.0), + max: Point::new(d.0 * 0.5, d.1 * 0.5, 0.0), + } + }) + .collect(); + RectangleLight { + dimensions: dimensions, + colors: colors, + bounds_: bbs, + } + } +} + +impl LightSource for RectangleLight { + fn sample(&self, + space: &Matrix4x4, + arr: Point, + u: f32, + v: f32, + wavelength: f32, + time: f32) + -> (SpectralSample, Vector, f32) { + // Calculate time interpolated values + let dim = lerp_slice(&self.dimensions, time); + let col = lerp_slice(&self.colors, time); + + // TODO: Is this right? Do we need to get the surface area post-transform? + let surface_area_inv: f64 = 1.0 / (dim.0 as f64 * dim.1 as f64); + + // Get the four corners of the rectangle, transformed into world space + let space_inv = space.inverse(); + let p1 = Point::new(dim.0 * 0.5, dim.1 * 0.5, 0.0) * space_inv; + let p2 = Point::new(dim.0 * -0.5, dim.1 * 0.5, 0.0) * space_inv; + let p3 = Point::new(dim.0 * -0.5, dim.1 * -0.5, 0.0) * space_inv; + let p4 = Point::new(dim.0 * 0.5, dim.1 * -0.5, 0.0) * space_inv; + + // Get the four corners of the rectangle, projected on to the unit + // sphere centered around arr. + let sp1 = (p1 - arr).normalized(); + let sp2 = (p2 - arr).normalized(); + let sp3 = (p3 - arr).normalized(); + let sp4 = (p4 - arr).normalized(); + + // Get the solid angles of the rectangle split into two triangles + let area_1 = spherical_triangle_solid_angle(sp2, sp1, sp3); + let area_2 = spherical_triangle_solid_angle(sp4, sp1, sp3); + + // Normalize the solid angles for selection purposes + let prob_1 = area_1 / (area_1 + area_2); + let prob_2 = 1.0 - prob_1; + + // Select one of the triangles and sample it + let shadow_vec = if u < prob_1 { + uniform_sample_spherical_triangle(sp2, sp1, sp3, v, u / prob_1) + } else { + uniform_sample_spherical_triangle(sp4, sp1, sp3, v, 1.0 - ((u - prob_1) / prob_2)) + }; + + // Project shadow_vec back onto the light's surface + let arr_local = arr * *space; + let shadow_vec_local = shadow_vec * *space; + let shadow_vec_local = shadow_vec_local * (-arr_local[2] / shadow_vec_local[2]); + let mut sample_point_local = arr_local + shadow_vec_local; + sample_point_local[0] = sample_point_local[0].max(dim.0 * -0.5).min(dim.0 * 0.5); + sample_point_local[1] = sample_point_local[1].max(dim.1 * -0.5).min(dim.1 * 0.5); + sample_point_local[2] = 0.0; + let sample_point = sample_point_local * space_inv; + let shadow_vec = sample_point - arr; + + // Calculate pdf and light energy + let pdf = 1.0 / (area_1 + area_2); // PDF of the ray direction being sampled + let spectral_sample = (col * surface_area_inv as f32 * 0.5).to_spectral_sample(wavelength); + + return (spectral_sample, shadow_vec, pdf as f32); + } + + fn sample_pdf(&self, + space: &Matrix4x4, + arr: Point, + sample_dir: Vector, + sample_u: f32, + sample_v: f32, + wavelength: f32, + time: f32) + -> f32 { + let dim = lerp_slice(&self.dimensions, time); + + // Get the four corners of the rectangle, transformed into world space + let space_inv = space.inverse(); + let p1 = Point::new(dim.0 * 0.5, dim.1 * 0.5, 0.0) * space_inv; + let p2 = Point::new(dim.0 * -0.5, dim.1 * 0.5, 0.0) * space_inv; + let p3 = Point::new(dim.0 * -0.5, dim.1 * -0.5, 0.0) * space_inv; + let p4 = Point::new(dim.0 * 0.5, dim.1 * -0.5, 0.0) * space_inv; + + // Get the four corners of the rectangle, projected on to the unit + // sphere centered around arr. + let sp1 = (p1 - arr).normalized(); + let sp2 = (p2 - arr).normalized(); + let sp3 = (p3 - arr).normalized(); + let sp4 = (p4 - arr).normalized(); + + // Get the solid angles of the rectangle split into two triangles + let area_1 = spherical_triangle_solid_angle(sp2, sp1, sp3); + let area_2 = spherical_triangle_solid_angle(sp4, sp1, sp3); + + 1.0 / (area_1 + area_2) + } + + fn outgoing(&self, + space: &Matrix4x4, + dir: Vector, + u: f32, + v: f32, + wavelength: f32, + time: f32) + -> SpectralSample { + let dim = lerp_slice(&self.dimensions, time); + let col = lerp_slice(&self.colors, time); + + // TODO: Is this right? Do we need to get the surface area post-transform? + let surface_area_inv: f64 = 1.0 / (dim.0 as f64 * dim.1 as f64); + + (col * surface_area_inv as f32 * 0.5).to_spectral_sample(wavelength) + } + + fn is_delta(&self) -> bool { + false + } +} + +impl Boundable for RectangleLight { + fn bounds<'a>(&'a self) -> &'a [BBox] { + &self.bounds_ + } +} diff --git a/src/math/normal.rs b/src/math/normal.rs index 006970d..1d3a537 100644 --- a/src/math/normal.rs +++ b/src/math/normal.rs @@ -1,6 +1,6 @@ #![allow(dead_code)] -use std::ops::{Index, IndexMut, Add, Sub, Mul, Div}; +use std::ops::{Index, IndexMut, Add, Sub, Mul, Div, Neg}; use std::cmp::PartialEq; use lerp::Lerp; @@ -114,6 +114,15 @@ impl Div for Normal { } +impl Neg for Normal { + type Output = Normal; + + fn neg(self) -> Normal { + Normal { co: self.co * -1.0 } + } +} + + impl Lerp for Normal { fn lerp(self, other: Normal, alpha: f32) -> Normal { (self * (1.0 - alpha)) + (other * alpha) diff --git a/src/math/point.rs b/src/math/point.rs index 42098c6..05f7574 100644 --- a/src/math/point.rs +++ b/src/math/point.rs @@ -39,6 +39,10 @@ impl Point { Point { co: n1.co.v_max(n2.co) } } + + pub fn into_vector(self) -> Vector { + Vector::new(self[0], self[1], self[2]) + } } diff --git a/src/math/vector.rs b/src/math/vector.rs index 669f578..177bea0 100644 --- a/src/math/vector.rs +++ b/src/math/vector.rs @@ -1,6 +1,6 @@ #![allow(dead_code)] -use std::ops::{Index, IndexMut, Add, Sub, Mul, Div}; +use std::ops::{Index, IndexMut, Add, Sub, Mul, Div, Neg}; use std::cmp::PartialEq; use lerp::Lerp; @@ -114,6 +114,15 @@ impl Div for Vector { } +impl Neg for Vector { + type Output = Vector; + + fn neg(self) -> Vector { + Vector { co: self.co * -1.0 } + } +} + + impl Lerp for Vector { fn lerp(self, other: Vector, alpha: f32) -> Vector { (self * (1.0 - alpha)) + (other * alpha) diff --git a/src/parse/psy_assembly.rs b/src/parse/psy_assembly.rs index 6080702..c9fbc43 100644 --- a/src/parse/psy_assembly.rs +++ b/src/parse/psy_assembly.rs @@ -5,7 +5,7 @@ use std::result::Result; use super::DataTree; use super::psy::{parse_matrix, PsyParseError}; use super::psy_mesh_surface::parse_mesh_surface; -use super::psy_light::parse_sphere_light; +use super::psy_light::{parse_sphere_light, parse_rectangle_light}; use assembly::{Assembly, AssemblyBuilder, Object}; @@ -82,6 +82,19 @@ pub fn parse_assembly(tree: &DataTree) -> Result { } } + // Sphere Light + "RectangleLight" => { + if let &DataTree::Internal { ident: Some(ident), .. } = child { + builder.add_object(ident, + Object::Light(Box::new( + try!(parse_rectangle_light(&child)) + ))); + } else { + // TODO: error condition of some kind, because no ident + panic!(); + } + } + _ => { // TODO: some kind of error, because not a known type name } diff --git a/src/parse/psy_light.rs b/src/parse/psy_light.rs index bb685f9..d5de6bb 100644 --- a/src/parse/psy_light.rs +++ b/src/parse/psy_light.rs @@ -8,7 +8,7 @@ use super::DataTree; use super::basics::{ws_usize, ws_f32}; use super::psy::PsyParseError; -use light::SphereLight; +use light::{SphereLight, RectangleLight}; use math::Point; use color::XYZ; @@ -54,3 +54,47 @@ pub fn parse_sphere_light(tree: &DataTree) -> Result return Err(PsyParseError::UnknownError); } } + +pub fn parse_rectangle_light(tree: &DataTree) -> Result { + if let &DataTree::Internal { ref children, .. } = tree { + let mut dimensions = Vec::new(); + let mut colors = Vec::new(); + + // Parse + for child in children.iter() { + match child { + // Dimensions + &DataTree::Leaf { type_name, contents } if type_name == "Dimensions" => { + if let IResult::Done(_, radius) = + closure!(tuple!(ws_f32, ws_f32))(contents.as_bytes()) { + dimensions.push(radius); + } else { + // Found dimensions, but its contents is not in the right format + return Err(PsyParseError::UnknownError); + } + } + + // Color + &DataTree::Leaf { type_name, contents } if type_name == "Color" => { + if let IResult::Done(_, color) = closure!(tuple!(ws_f32, + ws_f32, + ws_f32))(contents.as_bytes()) { + // TODO: handle color space conversions properly. + // Probably will need a special color type with its + // own parser...? + colors.push(XYZ::new(color.0, color.1, color.2)); + } else { + // Found color, but its contents is not in the right format + return Err(PsyParseError::UnknownError); + } + } + + _ => {} + } + } + + return Ok(RectangleLight::new(dimensions, colors)); + } else { + return Err(PsyParseError::UnknownError); + } +} diff --git a/src/renderer.rs b/src/renderer.rs index 40819b3..34017fc 100644 --- a/src/renderer.rs +++ b/src/renderer.rs @@ -260,8 +260,12 @@ impl LightPath { let (_, shadow_vec, _) = light.sample(&space, pos, lu, lv, self.wavelength, self.time); - let la = dot(nor.normalized().into_vector(), shadow_vec.normalized()) - .max(0.0); + let rnor = if dot(nor.into_vector(), ray.dir) > 0.0 { + -nor.into_vector().normalized() + } else { + nor.into_vector().normalized() + }; + let la = dot(rnor, shadow_vec.normalized()).max(0.0); self.light_attenuation = XYZ::from_spectral_sample(&XYZ::new(la, la, la) .to_spectral_sample(self.wavelength)); *ray = Ray::new(pos + shadow_vec.normalized() * 0.0001, diff --git a/src/sampling.rs b/src/sampling.rs index dda2157..3d9e419 100644 --- a/src/sampling.rs +++ b/src/sampling.rs @@ -1,4 +1,4 @@ -use math::Vector; +use math::{Vector, dot}; use std::f64::consts::PI as PI_64; use std::f32::consts::PI as PI_32; @@ -69,3 +69,96 @@ 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)) } + +/// Calculates the projected solid angle of a spherical triangle. +/// +/// A, B, and C are the points of the triangle on a unit sphere. +pub fn spherical_triangle_solid_angle(va: Vector, vb: Vector, vc: Vector) -> f32 { + // Calculate sines and cosines of the spherical triangle's edge lengths + let cos_a: f64 = dot(vb, vc).max(-1.0).min(1.0) as f64; + let cos_b: f64 = dot(vc, va).max(-1.0).min(1.0) as f64; + let cos_c: f64 = dot(va, vb).max(-1.0).min(1.0) as f64; + let sin_a: f64 = (1.0 - (cos_a * cos_a)).sqrt(); + let sin_b: f64 = (1.0 - (cos_b * cos_b)).sqrt(); + let sin_c: f64 = (1.0 - (cos_c * cos_c)).sqrt(); + + // If two of the vertices are coincident, area is zero. + // Return early to avoid a divide by zero below. + if cos_a == 1.0 || cos_b == 1.0 || cos_c == 1.0 { + return 0.0; + } + + // Calculate the cosine of the angles at the vertices + let cos_va = ((cos_a - (cos_b * cos_c)) / (sin_b * sin_c)).max(-1.0).min(1.0); + let cos_vb = ((cos_b - (cos_c * cos_a)) / (sin_c * sin_a)).max(-1.0).min(1.0); + let cos_vc = ((cos_c - (cos_a * cos_b)) / (sin_a * sin_b)).max(-1.0).min(1.0); + + // Calculate the angles themselves, in radians + let ang_va = cos_va.acos(); + let ang_vb = cos_vb.acos(); + let ang_vc = cos_vc.acos(); + + // Calculate and return the solid angle of the triangle + (ang_va + ang_vb + ang_vc - PI_64) as f32 +} + +/// Generates a uniform sample on a spherical triangle given two uniform +/// random variables i and j in [0, 1]. +pub fn uniform_sample_spherical_triangle(va: Vector, + vb: Vector, + vc: Vector, + i: f32, + j: f32) + -> Vector { + // Calculate sines and cosines of the spherical triangle's edge lengths + let cos_a: f64 = dot(vb, vc).max(-1.0).min(1.0) as f64; + let cos_b: f64 = dot(vc, va).max(-1.0).min(1.0) as f64; + let cos_c: f64 = dot(va, vb).max(-1.0).min(1.0) as f64; + let sin_a: f64 = (1.0 - (cos_a * cos_a)).sqrt(); + let sin_b: f64 = (1.0 - (cos_b * cos_b)).sqrt(); + let sin_c: f64 = (1.0 - (cos_c * cos_c)).sqrt(); + + // If two of the vertices are coincident, area is zero. + // Return early to avoid a divide by zero below. + if cos_a == 1.0 || cos_b == 1.0 || cos_c == 1.0 { + // TODO: do something more intelligent here, in the case that it's + // an infinitely thin line. + return va; + } + + // Calculate the cosine of the angles at the vertices + let cos_va = ((cos_a - (cos_b * cos_c)) / (sin_b * sin_c)).max(-1.0).min(1.0); + let cos_vb = ((cos_b - (cos_c * cos_a)) / (sin_c * sin_a)).max(-1.0).min(1.0); + let cos_vc = ((cos_c - (cos_a * cos_b)) / (sin_a * sin_b)).max(-1.0).min(1.0); + + // Calculate sine for A + let sin_va = (1.0 - (cos_va * cos_va)).sqrt(); + + // Calculate the angles themselves, in radians + let ang_va = cos_va.acos(); + let ang_vb = cos_vb.acos(); + let ang_vc = cos_vc.acos(); + + // Calculate the area of the spherical triangle + let area = ang_va + ang_vb + ang_vc - PI_64; + + // The rest of this is from the paper "Stratified Sampling of Spherical + // Triangles" by James Arvo. + let area_2 = area * i as f64; + + let s = (area_2 - ang_va).sin(); + let t = (area_2 - ang_va).cos(); + let u = t - cos_va; + let v = s + (sin_va * cos_c); + + let q_top = (((v * t) - (u * s)) * cos_va) - v; + let q_bottom = ((v * s) + (u * t)) * sin_va; + let q = q_top / q_bottom; + + let vc_2 = (va * q as f32) + + ((vc - (va * dot(vc, va))).normalized() * (1.0 - (q * q)).sqrt() as f32); + + let z = 1.0 - (j * (1.0 - dot(vc_2, vb))); + + (vb * z) + ((vc_2 - (vb * dot(vc_2, vb))).normalized() * (1.0 - (z * z)).sqrt()) +}