From 0a7963097ce65c8a5a0d1047b440855591139ce4 Mon Sep 17 00:00:00 2001 From: Matthew Gordon Date: Wed, 27 Nov 2019 17:02:23 -0500 Subject: [PATCH] Add (co)tangent to IntersectionInfo and redid sphere intersection --- src/camera.rs | 2 + src/raycasting.rs | 115 +++++++++++++++++++++++++++++++++++++++++----- 2 files changed, 105 insertions(+), 12 deletions(-) diff --git a/src/camera.rs b/src/camera.rs index 5ac468b..f771cd6 100644 --- a/src/camera.rs +++ b/src/camera.rs @@ -153,6 +153,8 @@ mod tests { location, distance: _, normal: _, + tangent: _, + cotangent: _, retro: _, material: _, }) => location, diff --git a/src/raycasting.rs b/src/raycasting.rs index ef86849..b53c7ac 100644 --- a/src/raycasting.rs +++ b/src/raycasting.rs @@ -32,6 +32,8 @@ pub struct IntersectionInfo { pub distance: T, pub location: Vector3, pub normal: Vector3, + pub tangent: Vector3, + pub cotangent: Vector3, pub retro: Vector3, pub material: Rc>, } @@ -58,23 +60,30 @@ impl Sphere { impl Intersect for Sphere { fn intersect<'a>(&'a self, ray: &Ray) -> Option> { - let ray_origin_to_sphere_centre = self.centre - ray.origin; + /*let ray_origin_to_sphere_centre = self.centre - ray.origin; let radius_squared = self.radius * self.radius; let is_inside_sphere = ray_origin_to_sphere_centre.norm_squared() <= radius_squared; // t0/p0 is the point on the ray that's closest to the centre of the sphere + // ray.direction is normalized, so it's not necessary to divide by its length. let t0 = ray_origin_to_sphere_centre.dot(&ray.direction); if !is_inside_sphere && t0 < T::zero() { // Sphere is behind ray origin return None; } // Squared distance between ray origin and sphere centre - let d0_squared = (ray.origin - self.centre).norm_squared(); - // Squared distance petween p0 and sphere centre + let d0_squared = (ray_origin_to_sphere_centre).norm_squared(); + // p0, ray.origin and sphere.centre form a right triangle, with p0 at the right corner, + // Squared distance petween p0 and sphere centre, using Pythagoras let p0_dist_from_centre_squared = d0_squared - t0 * t0; if p0_dist_from_centre_squared > radius_squared { // Sphere is in front of ray but ray misses return None; } + let p0_dist_from_centre =p0_dist_from_centre_squared.sqrt(); + // Two more right triangles are formed by p0, the sphere centre, and the two places + // where the ray intersects the sphere. (Or the ray may be a tangent to the sphere + // in which case these triangles are degenerate. Here we use Pythagoras again to find + .// find the distance between p0 and the two intersection points. let delta = (radius_squared - p0_dist_from_centre_squared).sqrt(); let distance = if is_inside_sphere { // radius origin is inside sphere @@ -84,19 +93,69 @@ impl Intersect for Sphere { }; let location = ray.point_at(distance); let normal = (location - self.centre).normalize(); - let retro = -ray.direction; - Some(IntersectionInfo { - distance, - location, - normal, - retro, - material: Rc::clone(&self.material), - }) + let tangent = normal.cross(&Vector3::z_axis()); + let cotangent = normal.cross(&tangent); + let retro = -ray.direction;*/ + let two: T = convert(2.0); + let four: T = convert(4.0); + let a = ray + .direction + .component_mul(&ray.direction) + .iter() + .fold(T::zero(), |a, b| a + *b); + let b = ((ray.origin.component_mul(&ray.direction) + - self.centre.component_mul(&ray.direction)) + * two) + .iter() + .fold(T::zero(), |a, b| a + *b); + let c = (ray.origin.component_mul(&ray.origin) + self.centre.component_mul(&self.centre) + - self.centre.component_mul(&ray.origin) * two) + .iter() + .fold(T::zero(), |a, b| a + *b) + - self.radius * self.radius; + let delta_squared: T = b * b - four * a * c; + if delta_squared < T::zero() { + None + } else { + let delta = delta_squared.sqrt(); + let one_over_2_a = T::one() / (two * a); + let t1 = (-b - delta) * one_over_2_a; + let t2 = (-b + delta) * one_over_2_a; + let distance = if t1 < T::zero() { + t2 + } else if t2 < T::zero() { + t1 + } else if t1 < t2 { + t1 + } else { + t2 + }; + if distance <= T::zero() { + None + } else { + let location = ray.point_at(distance); + let normal = (location - self.centre).normalize(); + let tangent = normal.cross(&Vector3::z_axis()); + let cotangent = normal.cross(&tangent); + let retro = -ray.direction; + Some(IntersectionInfo { + distance, + location, + normal, + tangent, + cotangent, + retro, + material: Rc::clone(&self.material), + }) + } + } } } pub struct Plane { normal: Vector3, + tangent: Vector3, + cotangent: Vector3, distance_from_origin: T, material: Rc>, } @@ -108,8 +167,14 @@ impl Plane { material: Rc>, ) -> Plane { normal.normalize(); + let mut axis_closest_to_tangent = Vector3::zeros(); + axis_closest_to_tangent[normal.iamin()] = T::one(); + let cotangent = normal.cross(&axis_closest_to_tangent); + let tangent = normal.cross(&cotangent); Plane { normal, + tangent, + cotangent, distance_from_origin, material, } @@ -137,6 +202,8 @@ impl Intersect for Plane { distance: t, location: ray.point_at(t), normal: self.normal, + tangent: self.tangent, + cotangent: self.cotangent, retro: -ray.direction, material: Rc::clone(&self.material), }) @@ -159,7 +226,7 @@ mod tests { use super::*; use crate::materials::LambertianMaterial; - use quickcheck::{Arbitrary, Gen}; + use quickcheck::{Arbitrary, Gen, TestResult}; impl Arbitrary for Ray { fn arbitrary(g: &mut G) -> Ray { let origin = as Arbitrary>::arbitrary(g); @@ -240,6 +307,28 @@ mod tests { assert_matches!(s.intersect(&r), Some(_)); } + #[quickcheck] + fn ray_intersects_sphere_centre_at_correct_distance( + ray_origin: Vector3, + sphere_centre: Vector3, + radius: f64, + ) -> TestResult { + if radius <= 0.0 || radius + 0.000001 >= (ray_origin - sphere_centre).norm() { + return TestResult::discard(); + }; + let sphere = Sphere::new( + sphere_centre, + radius, + Rc::new(LambertianMaterial::new_dummy()), + ); + let ray = Ray::new(ray_origin, sphere_centre - ray_origin); + let info = sphere.intersect(&ray).unwrap(); + let distance_to_centre = (sphere_centre - ray.origin).norm(); + TestResult::from_bool( + (distance_to_centre - (info.distance + sphere.radius)).abs() < 0.00001, + ) + } + #[test] fn ray_intersects_plane() { let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(-1.0, 0.0, 1.0)); @@ -275,6 +364,8 @@ mod tests { distance: _, location, normal: _, + tangent: _, + cotangent: _, retro: _, material: _, }) => assert!((location.x - (-5.0f64)).abs() < 0.0000000001),