use nalgebra::{convert, RealField, Vector3}; use super::materials::Material; use std::rc::Rc; #[derive(Clone, Debug)] pub struct Ray { origin: Vector3, direction: Vector3, } impl Ray { pub fn new(origin: Vector3, direction: Vector3) -> Ray { Ray { origin, direction: direction.normalize(), } } pub fn point_at(&self, t: T) -> Vector3 { return self.origin + self.direction * t; } pub fn bias(&self, amount: T) -> Ray { Ray::new(self.origin + self.direction * amount, self.direction) } } #[derive(Debug)] pub struct IntersectionInfo { pub distance: T, pub location: Vector3, pub normal: Vector3, pub tangent: Vector3, pub cotangent: Vector3, pub retro: Vector3, pub material: Rc>, } pub trait Intersect { fn intersect<'a>(&'a self, ray: &Ray) -> Option>; } pub struct Sphere { centre: Vector3, radius: T, material: Rc>, } impl Sphere { pub fn new(centre: Vector3, radius: T, material: Rc>) -> Sphere { Sphere { centre, radius, material, } } } impl Intersect for Sphere { fn intersect<'a>(&'a self, ray: &Ray) -> Option> { /*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_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 t0 + delta } else { t0 - delta }; 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;*/ 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>, } impl Plane { pub fn new( normal: Vector3, distance_from_origin: T, 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, } } } impl Intersect for Plane { fn intersect<'a>(&'a self, ray: &Ray) -> Option> { let ray_direction_dot_plane_normal = ray.direction.dot(&self.normal); let point_on_plane = self.normal * self.distance_from_origin; let point_on_plane_minus_ray_origin_dot_normal = (point_on_plane - ray.origin).dot(&self.normal); if ray_direction_dot_plane_normal == convert(0.0) { //Ray is parallel to plane if point_on_plane_minus_ray_origin_dot_normal != convert(0.0) { //Ray is not in plane return None; } } let t = point_on_plane_minus_ray_origin_dot_normal / ray_direction_dot_plane_normal; if t < convert(0.0) { return None; } Some(IntersectionInfo { distance: t, location: ray.point_at(t), normal: self.normal, tangent: self.tangent, cotangent: self.cotangent, retro: -ray.direction, material: Rc::clone(&self.material), }) } } #[cfg(test)] mod tests { use quickcheck_macros::quickcheck; macro_rules! assert_matches { ($expression:expr, $($pattern:tt)+) => { match $expression { $($pattern)+ => (), ref e => panic!("assertion failed: `{:?}` does not match `{}`", e, stringify!($($pattern)+)), } } } use super::*; use crate::materials::LambertianMaterial; use quickcheck::{Arbitrary, Gen, TestResult}; impl Arbitrary for Ray { fn arbitrary(g: &mut G) -> Ray { let origin = as Arbitrary>::arbitrary(g); let direction = as Arbitrary>::arbitrary(g); return Ray::new(origin, direction); } } #[quickcheck] fn t0_is_origin(ray: Ray) -> bool { ray.point_at(0.0) == ray.origin } #[quickcheck] fn t1_is_origin_plus_direction(ray: Ray) -> bool { ray.point_at(1.0) == ray.origin + ray.direction } #[quickcheck] fn points_are_colinear(ray: Ray, t1: f64, t2: f64, t3: f64) -> bool { let p1 = ray.point_at(t1); let p2 = ray.point_at(t2); let p3 = ray.point_at(t3); let epsilon = [t1, t2, t3, ray.origin[0], ray.origin[1], ray.origin[2]] .iter() .fold(0.0, |a, &b| a.max(b.abs())) * std::f64::EPSILON * 256.0; (p2 - p1).cross(&(p3 - p2)).norm() < epsilon } #[quickcheck] fn t_is_distance(ray: Ray, t: f64) -> bool { (ray.point_at(t) - ray.origin).norm() - t.abs() < 0.0000000001 } #[test] fn ray_intersects_sphere() { let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0)); let s = Sphere::new( Vector3::new(1.5, 1.5, 15.0), 5.0, Rc::new(LambertianMaterial::new_dummy()), ); assert_matches!(s.intersect(&r), Some(_)); } #[test] fn ray_does_not_intersect_sphere_when_sphere_is_in_front() { let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0)); let s = Sphere::new( Vector3::new(-5.0, 1.5, 15.0), 5.0, Rc::new(LambertianMaterial::new_dummy()), ); assert_matches!(s.intersect(&r), None); } #[test] fn ray_does_not_intersect_sphere_when_sphere_is_behind() { let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0)); let s = Sphere::new( Vector3::new(1.5, 1.5, -15.0), 5.0, Rc::new(LambertianMaterial::new_dummy()), ); assert_matches!(s.intersect(&r), None); } #[test] fn ray_intersects_sphere_when_origin_is_inside() { let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0)); let s = Sphere::new( Vector3::new(1.5, 1.5, 2.0), 5.0, Rc::new(LambertianMaterial::new_dummy()), ); 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)); let p = Plane::new( Vector3::new(1.0, 0.0, 0.0), -5.0, Rc::new(LambertianMaterial::new_dummy()), ); assert_matches!(p.intersect(&r), Some(_)); } #[test] fn ray_does_not_intersect_plane() { let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(1.0, 0.0, 1.0)); let p = Plane::new( Vector3::new(1.0, 0.0, 0.0), -5.0, Rc::new(LambertianMaterial::new_dummy()), ); assert_matches!(p.intersect(&r), None); } #[test] fn intersection_point_is_on_plane() { let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(-1.0, 0.0, 1.0)); let p = Plane::new( Vector3::new(1.0, 0.0, 0.0), -5.0, Rc::new(LambertianMaterial::new_dummy()), ); match p.intersect(&r) { Some(IntersectionInfo { distance: _, location, normal: _, tangent: _, cotangent: _, retro: _, material: _, }) => assert!((location.x - (-5.0f64)).abs() < 0.0000000001), None => panic!(), } } }