From a7b9ec249a340a1264756019c93ab55d9b0bdc17 Mon Sep 17 00:00:00 2001 From: Michael Zhang Date: Thu, 2 Feb 2023 02:08:17 -0600 Subject: [PATCH] Refactored large part of the code to make it more extensible --- assignment-1/src/input_file.rs | 16 +- assignment-1/src/main.rs | 20 +- assignment-1/src/math.rs | 21 ++ assignment-1/src/ray.rs | 211 +----------------- assignment-1/src/scene/cylinder.rs | 121 ++++++++++ .../src/{scene_data.rs => scene/data.rs} | 43 ++-- assignment-1/src/scene/mod.rs | 3 + assignment-1/src/scene/sphere.rs | 92 ++++++++ 8 files changed, 283 insertions(+), 244 deletions(-) create mode 100644 assignment-1/src/math.rs create mode 100644 assignment-1/src/scene/cylinder.rs rename assignment-1/src/{scene_data.rs => scene/data.rs} (72%) create mode 100644 assignment-1/src/scene/mod.rs create mode 100644 assignment-1/src/scene/sphere.rs diff --git a/assignment-1/src/input_file.rs b/assignment-1/src/input_file.rs index 419abef..7d68f19 100644 --- a/assignment-1/src/input_file.rs +++ b/assignment-1/src/input_file.rs @@ -5,7 +5,11 @@ use nalgebra::Vector3; use crate::{ image::Color, - scene_data::{Cylinder, Object, ObjectKind, Scene, Sphere}, + scene::{ + cylinder::Cylinder, + data::{Object, Scene}, + sphere::Sphere, + }, }; /// Parse the input file into a scene @@ -50,7 +54,11 @@ pub fn parse_input_file(path: impl AsRef) -> Result { if parts.len() < 3 { bail!("Vec3 requires 3 components."); } - Ok(Vector3::new(parts[start], parts[start + 1], parts[start + 2])) + Ok(Vector3::new( + parts[start], + parts[start + 1], + parts[start + 2], + )) }; let read_color = || { @@ -75,7 +83,7 @@ pub fn parse_input_file(path: impl AsRef) -> Result { } "sphere" => scene.objects.push(Object { - kind: ObjectKind::Sphere(Sphere { + kind: Box::new(Sphere { center: read_vec3(0)?, radius: parts[3], }), @@ -86,7 +94,7 @@ pub fn parse_input_file(path: impl AsRef) -> Result { }), "cylinder" => scene.objects.push(Object { - kind: ObjectKind::Cylinder(Cylinder { + kind: Box::new(Cylinder { center: read_vec3(0)?, direction: read_vec3(3)?, radius: parts[6], diff --git a/assignment-1/src/main.rs b/assignment-1/src/main.rs index d1dc9be..588a297 100644 --- a/assignment-1/src/main.rs +++ b/assignment-1/src/main.rs @@ -3,8 +3,9 @@ extern crate anyhow; mod image; mod input_file; +mod math; mod ray; -mod scene_data; +mod scene; use std::fs::File; use std::path::PathBuf; @@ -13,7 +14,6 @@ use anyhow::Result; use clap::Parser; use ordered_float::NotNan; use rayon::prelude::{IntoParallelIterator, ParallelIterator}; -use scene_data::ObjectKind; use crate::image::Image; use crate::input_file::parse_input_file; @@ -107,11 +107,7 @@ fn main() -> Result<()> { .objects .iter() .filter_map(|object| { - use ObjectKind::*; - let intersection_point_opt = match &object.kind { - Sphere(sphere) => ray.intersects_sphere_at(&sphere), - Cylinder(cylinder) => ray.intersects_cylinder_at(&cylinder), - }; + let intersection_point_opt = object.kind.intersects_ray_at(&ray); intersection_point_opt.and_then(|t| { // Unfortunately, IEEE floats in Rust don't have total ordering, @@ -128,16 +124,14 @@ fn main() -> Result<()> { // Sort the list of intersection times by the lowest one. .min_by_key(|(t, _)| *t); - let pixel_color = match earliest_intersection { + match earliest_intersection { // Take the object's material color - Some((_, object)) => scene.material_colors[object.material], + Some((_, object)) => scene.compute_pixel_color(object.material), // There was no intersection, so this should default to the scene's // background color None => scene.bkg_color, - }; - - pixel_color + } }) .collect::>(); @@ -149,7 +143,7 @@ fn main() -> Result<()> { }; { - let file = File::create(&out_file)?; + let file = File::create(out_file)?; image.write(file)?; } diff --git a/assignment-1/src/math.rs b/assignment-1/src/math.rs new file mode 100644 index 000000000..e04f80b --- /dev/null +++ b/assignment-1/src/math.rs @@ -0,0 +1,21 @@ +use nalgebra::{Matrix3, Vector3}; + +/// Calculate the rotation matrix between the 2 given vectors +/// Based on the method here: https://math.stackexchange.com/a/897677 +pub fn compute_rotation_matrix( + a: Vector3, + b: Vector3, +) -> Matrix3 { + let cos_t = a.dot(&b); + let sin_t = a.cross(&b).norm(); + let g = Matrix3::new(cos_t, -sin_t, 0.0, sin_t, cos_t, 0.0, 0.0, 0.0, 1.0); + + // New basis vectors + let u = a; + let v = (b - a.dot(&b) * a).normalize(); + let w = b.cross(&a); + // Not sure if this is required to be invertible? + let f = Matrix3::from_columns(&[u, v, w]).try_inverse().unwrap(); + + f.try_inverse().unwrap() * g * f +} diff --git a/assignment-1/src/ray.rs b/assignment-1/src/ray.rs index be6bb45..a722ddf 100644 --- a/assignment-1/src/ray.rs +++ b/assignment-1/src/ray.rs @@ -1,7 +1,4 @@ -use nalgebra::{Matrix3, Vector3}; -use ordered_float::NotNan; - -use crate::scene_data::{Cylinder, Sphere}; +use nalgebra::Vector3; /// A normalized parametric Ray of the form (origin + direction * time) /// @@ -9,8 +6,8 @@ use crate::scene_data::{Cylinder, Sphere}; /// time occurs on the ray. #[derive(Debug)] pub struct Ray { - origin: Vector3, - direction: Vector3, + pub origin: Vector3, + pub direction: Vector3, } impl Ray { @@ -27,206 +24,4 @@ impl Ray { pub fn eval(&self, time: f64) -> Vector3 { self.origin + self.direction * time } - - /// Given a sphere, returns the first time at which this ray intersects the - /// sphere. - /// - /// If there is no intersection point, returns None. - pub fn intersects_sphere_at(&self, sphere: &Sphere) -> Option { - let a = self.direction.x.powi(2) - + self.direction.y.powi(2) - + self.direction.z.powi(2); - let b = 2.0 - * (self.direction.x * (self.origin.x - sphere.center.x) - + self.direction.y * (self.origin.y - sphere.center.y) - + self.direction.z * (self.origin.z - sphere.center.z)); - let c = (self.origin.x - sphere.center.x).powi(2) - + (self.origin.y - sphere.center.y).powi(2) - + (self.origin.z - sphere.center.z).powi(2) - - sphere.radius.powi(2); - let discriminant = b * b - 4.0 * a * c; - - match discriminant { - // Discriminant < 0, means the equation has no solutions. - d if d < 0.0 => return None, - - // Discriminant == 0 - d if d == 0.0 => { - return Some(-b / (2.0 * a)); - } - - d if d > 0.0 => { - let solution_1 = (-b + discriminant.sqrt()) / (2.0 * a); - let solution_2 = (-b - discriminant.sqrt()) / (2.0 * a); - - return Some(solution_1.min(solution_2)); - } - - // Probably hit some NaN or Infinity value due to faulty inputs... - _ => unreachable!("Invalid determinant value: {discriminant}"), - } - } - - /// Given a cylinder, returns the first time at which this ray intersects the - /// cylinder. - /// - /// If there is no intersection point, returns None. - pub fn intersects_cylinder_at(&self, cylinder: &Cylinder) -> Option { - // Determine rotation matrix for turning the cylinder upright along the - // Z-axis - let target_direction = Vector3::new(0.0, 0.0, 1.0); - let rotation_matrix = - compute_rotation_matrix(cylinder.direction, target_direction); - - // Transform all parameters according to this rotation matrix - let rotated_cylinder_center = rotation_matrix * cylinder.center; - let rotated_ray_origin = rotation_matrix * self.origin; - let rotated_ray_direction = rotation_matrix * self.direction; - - // Now that we know the cylinder is upright, we can start checking against - // the formula: - // - // (ox + t*rx - cx)^2 + (oy + t*ry - cy)^2 = r^2 - // - // where o{xy} is the ray origin, r{xy} is the ray direction, and c{xy} is - // the cylinder center. The z will be taken care of after the fact. To - // solve, we must put it into the form At^2 + Bt + c = 0. The variables - // are: - // - // A: rx^2 + ry^2 - // B: 2(rx(ox - cx) + ry(oy - cy)) - // C: (cx - ox)^2 + (cy - oy)^2 - r^2 - let (a, b, c) = { - let o = rotated_ray_origin; - let r = rotated_ray_direction; - let c = rotated_cylinder_center; - - ( - r.x.powi(2) + r.y.powi(2), - 2.0 * (r.x * (o.x - c.x) + r.y * (o.y - c.y)), - (c.x - o.x).powi(2) + (c.y - o.y).powi(2) - cylinder.radius.powi(2), - ) - }; - - let discriminant = b * b - 4.0 * a * c; - - let mut solutions = match discriminant { - // Discriminant < 0, means the equation has no solutions. - d if d < 0.0 => vec![], - - // Discriminant == 0 - d if d == 0.0 => vec![-b / 2.0 * a], - - // Discriminant > 0, 2 solutions available. - d if d > 0.0 => { - vec![ - (-b + discriminant.sqrt()) / (2.0 * a), - (-b - discriminant.sqrt()) / (2.0 * a), - ] - } - - // Probably hit some NaN or Infinity value due to faulty inputs... - _ => unreachable!("Invalid determinant value: {discriminant}"), - }; - - // We also need to add solutions for the two ends of the cylinder, which - // uses a similar method except backwards: check intersection points - // with the correct z-plane and then see if the points are within the - // circle. - // - // Luckily, this means we only need to care about one dimension at first, - // and don't need to perform the quadratic equation method above. - // - // oz + t * rz = cz +- (len / 2) - // t = (oz + cz +- (len / 2)) / rz - let possible_z_intersections = { - let o = rotated_ray_origin; - let r = rotated_ray_direction; - let c = rotated_cylinder_center; - - vec![ - (o.z + c.z + cylinder.length / 2.0) / r.z, - (o.z + c.z - cylinder.length / 2.0) / r.z, - ] - }; - // Filter out all the solutions where the z does not lie in the circle - solutions.extend(possible_z_intersections.into_iter().filter(|t| { - let ray_point = self.origin + self.direction * (*t); - ray_point.x.powi(2) + ray_point.y.powi(2) <= cylinder.radius.powi(2) - })); - - // Filter out solutions that don't have a valid Z position. - let solutions = solutions - .into_iter() - .filter(|t| { - let ray_point = self.origin + self.direction * (*t); - - let rotated_ray_point = rotation_matrix * ray_point; - let z = rotated_ray_point.z - rotated_cylinder_center.z; - - // Check to see if z is between -len/2 and len/2 - z.abs() < cylinder.length / 2.0 - }) - .filter_map(|t| NotNan::new(t).ok()); - - // Return the minimum solution - solutions.min().map(|t| t.into_inner()) - } -} - -/// Calculate the rotation matrix between the 2 given vectors -/// Based on the method here: https://math.stackexchange.com/a/897677 -fn compute_rotation_matrix(a: Vector3, b: Vector3) -> Matrix3 { - let cos_t = a.dot(&b); - let sin_t = a.cross(&b).norm(); - let G = Matrix3::new(cos_t, -sin_t, 0.0, sin_t, cos_t, 0.0, 0.0, 0.0, 1.0); - - // Basis - let u = a; - let v = (b - a.dot(&b) * a).normalize(); - let w = b.cross(&a); - let F = Matrix3::from_columns(&[u, v, w]).try_inverse().unwrap(); - - F.try_inverse().unwrap() * G * F -} - -#[cfg(test)] -mod tests { - use nalgebra::Vector3; - - use crate::scene_data::Sphere; - - use super::Ray; - - #[test] - fn practice_problem_slide_154() { - let ray = Ray { - origin: Vector3::new(0.0, 0.0, 0.0), - direction: Vector3::new(0.0, 0.0, -1.0), - }; - let sphere = Sphere { - center: Vector3::new(0.0, 0.0, -10.0), - radius: 4.0, - }; - - let point = ray.intersects_sphere_at(&sphere).map(|t| ray.eval(t)); - - // the intersection point in this case is (0, 0, -6) - assert_eq!(point, Some(Vector3::new(0.0, 0.0, -6.0))); - } - - #[test] - fn practice_problem_slide_158() { - let ray = Ray { - origin: Vector3::new(0.0, 0.0, 0.0), - direction: Vector3::new(0.0, 0.5, -1.0), - }; - let sphere = Sphere { - center: Vector3::new(0.0, 0.0, -10.0), - radius: 4.0, - }; - - // oops! In this case, the ray does not intersect the sphere. - assert_eq!(ray.intersects_sphere_at(&sphere), None); - } } diff --git a/assignment-1/src/scene/cylinder.rs b/assignment-1/src/scene/cylinder.rs new file mode 100644 index 000000000..7d388b2 --- /dev/null +++ b/assignment-1/src/scene/cylinder.rs @@ -0,0 +1,121 @@ +use nalgebra::Vector3; +use ordered_float::NotNan; + +use crate::{math::compute_rotation_matrix, ray::Ray}; + +use super::data::ObjectKind; + +#[derive(Debug)] +pub struct Cylinder { + pub center: Vector3, + pub direction: Vector3, + pub radius: f64, + pub length: f64, +} + +impl ObjectKind for Cylinder { + /// Given a cylinder, returns the first time at which this ray intersects the + /// cylinder. + /// + /// If there is no intersection point, returns None. + fn intersects_ray_at(&self, ray: &Ray) -> Option { + // Determine rotation matrix for turning the cylinder upright along the + // Z-axis + let target_direction = Vector3::new(0.0, 0.0, 1.0); + let rotation_matrix = + compute_rotation_matrix(self.direction, target_direction); + + // Transform all parameters according to this rotation matrix + let rotated_cylinder_center = rotation_matrix * self.center; + let rotated_ray_origin = rotation_matrix * ray.origin; + let rotated_ray_direction = rotation_matrix * ray.direction; + + // Now that we know the cylinder is upright, we can start checking against + // the formula: + // + // (ox + t*rx - cx)^2 + (oy + t*ry - cy)^2 = r^2 + // + // where o{xy} is the ray origin, r{xy} is the ray direction, and c{xy} is + // the cylinder center. The z will be taken care of after the fact. To + // solve, we must put it into the form At^2 + Bt + c = 0. The variables + // are: + // + // A: rx^2 + ry^2 + // B: 2(rx(ox - cx) + ry(oy - cy)) + // C: (cx - ox)^2 + (cy - oy)^2 - r^2 + let (a, b, c) = { + let o = rotated_ray_origin; + let r = rotated_ray_direction; + let c = rotated_cylinder_center; + + ( + r.x.powi(2) + r.y.powi(2), + 2.0 * (r.x * (o.x - c.x) + r.y * (o.y - c.y)), + (c.x - o.x).powi(2) + (c.y - o.y).powi(2) - self.radius.powi(2), + ) + }; + + let discriminant = b * b - 4.0 * a * c; + + let mut solutions = match discriminant { + // Discriminant < 0, means the equation has no solutions. + d if d < 0.0 => vec![], + + // Discriminant == 0 + d if d == 0.0 => vec![-b / 2.0 * a], + + // Discriminant > 0, 2 solutions available. + d if d > 0.0 => { + vec![ + (-b + discriminant.sqrt()) / (2.0 * a), + (-b - discriminant.sqrt()) / (2.0 * a), + ] + } + + // Probably hit some NaN or Infinity value due to faulty inputs... + _ => unreachable!("Invalid determinant value: {discriminant}"), + }; + + // We also need to add solutions for the two ends of the cylinder, which + // uses a similar method except backwards: check intersection points + // with the correct z-plane and then see if the points are within the + // circle. + // + // Luckily, this means we only need to care about one dimension at first, + // and don't need to perform the quadratic equation method above. + // + // oz + t * rz = cz +- (len / 2) + // t = (oz + cz +- (len / 2)) / rz + let possible_z_intersections = { + let o = rotated_ray_origin; + let r = rotated_ray_direction; + let c = rotated_cylinder_center; + + vec![ + (o.z + c.z + self.length / 2.0) / r.z, + (o.z + c.z - self.length / 2.0) / r.z, + ] + }; + // Filter out all the solutions where the z does not lie in the circle + solutions.extend(possible_z_intersections.into_iter().filter(|t| { + let ray_point = ray.eval(*t); + ray_point.x.powi(2) + ray_point.y.powi(2) <= self.radius.powi(2) + })); + + // Filter out solutions that don't have a valid Z position. + let solutions = solutions + .into_iter() + .filter(|t| { + let ray_point = ray.eval(*t); + let rotated_ray_point = rotation_matrix * ray_point; + let z = rotated_ray_point.z - rotated_cylinder_center.z; + + // Check to see if z is between -len/2 and len/2 + z.abs() < self.length / 2.0 + }) + .filter_map(|t| NotNan::new(t).ok()); + + // Return the minimum solution + solutions.min().map(|t| t.into_inner()) + } +} diff --git a/assignment-1/src/scene_data.rs b/assignment-1/src/scene/data.rs similarity index 72% rename from assignment-1/src/scene_data.rs rename to assignment-1/src/scene/data.rs index e24ffac..cbf3881 100644 --- a/assignment-1/src/scene_data.rs +++ b/assignment-1/src/scene/data.rs @@ -1,35 +1,27 @@ +use std::fmt::Debug; + use nalgebra::Vector3; use crate::image::Color; +use crate::ray::Ray; -#[derive(Debug)] -pub struct Sphere { - pub center: Vector3, - pub radius: f64, -} - -#[derive(Debug)] -pub struct Cylinder { - pub center: Vector3, - pub direction: Vector3, - pub radius: f64, - pub length: f64, +pub trait ObjectKind: Debug + Send + Sync { + /// Determine where the ray intersects this object, returning the earliest + /// time this happens. Returns None if no intersection occurs. + /// + /// Also known as Trace_Ray in the slides, except not the part where it calls + /// Shade_Ray. + fn intersects_ray_at(&self, ray: &Ray) -> Option; } #[derive(Debug)] pub struct Object { - pub kind: ObjectKind, + pub kind: Box, /// Index into the scene's material color list pub material: usize, } -#[derive(Debug)] -pub enum ObjectKind { - Sphere(Sphere), - Cylinder(Cylinder), -} - #[derive(Debug)] pub struct Rect { pub upper_left: Vector3, @@ -59,6 +51,19 @@ pub struct Scene { } impl Scene { + /// Determine the color that should be used to fill this pixel + /// + /// Also known as Shade_Ray in the slides. + pub fn compute_pixel_color(&self, material_idx: usize) -> Color { + // TODO: Does it make sense to make this function fallible from an API + // design standpoint? + self + .material_colors + .get(material_idx) + .cloned() + .unwrap_or(self.bkg_color) + } + /// Determine the boundaries of the viewing window in world coordinates pub fn compute_viewing_window(&self, distance: f64) -> Rect { // Compute viewing directions diff --git a/assignment-1/src/scene/mod.rs b/assignment-1/src/scene/mod.rs new file mode 100644 index 000000000..e0465c9 --- /dev/null +++ b/assignment-1/src/scene/mod.rs @@ -0,0 +1,3 @@ +pub mod data; +pub mod sphere; +pub mod cylinder; diff --git a/assignment-1/src/scene/sphere.rs b/assignment-1/src/scene/sphere.rs new file mode 100644 index 000000000..0acd760 --- /dev/null +++ b/assignment-1/src/scene/sphere.rs @@ -0,0 +1,92 @@ +use nalgebra::Vector3; + +use crate::ray::Ray; + +use super::data::ObjectKind; + +#[derive(Debug)] +pub struct Sphere { + pub center: Vector3, + pub radius: f64, +} + +impl ObjectKind for Sphere { + /// Given a sphere, returns the first time at which this ray intersects the + /// sphere. + /// + /// If there is no intersection point, returns None. + fn intersects_ray_at(&self, ray: &Ray) -> Option { + let a = ray.direction.x.powi(2) + + ray.direction.y.powi(2) + + ray.direction.z.powi(2); + let b = 2.0 + * (ray.direction.x * (ray.origin.x - self.center.x) + + ray.direction.y * (ray.origin.y - self.center.y) + + ray.direction.z * (ray.origin.z - self.center.z)); + let c = (ray.origin.x - self.center.x).powi(2) + + (ray.origin.y - self.center.y).powi(2) + + (ray.origin.z - self.center.z).powi(2) + - self.radius.powi(2); + let discriminant = b * b - 4.0 * a * c; + + match discriminant { + // Discriminant < 0, means the equation has no solutions. + d if d < 0.0 => None, + + // Discriminant == 0 + d if d == 0.0 => Some(-b / (2.0 * a)), + + d if d > 0.0 => { + let solution_1 = (-b + discriminant.sqrt()) / (2.0 * a); + let solution_2 = (-b - discriminant.sqrt()) / (2.0 * a); + + Some(solution_1.min(solution_2)) + } + + // Probably hit some NaN or Infinity value due to faulty inputs... + _ => unreachable!("Invalid determinant value: {discriminant}"), + } + } +} + +#[cfg(test)] +mod tests { + use nalgebra::Vector3; + + use crate::ray::Ray; + use crate::scene::data::ObjectKind; + + use super::Sphere; + + #[test] + fn practice_problem_slide_154() { + let ray = Ray { + origin: Vector3::new(0.0, 0.0, 0.0), + direction: Vector3::new(0.0, 0.0, -1.0), + }; + let sphere = Sphere { + center: Vector3::new(0.0, 0.0, -10.0), + radius: 4.0, + }; + + let point = sphere.intersects_ray_at(&ray).map(|t| ray.eval(t)); + + // the intersection point in this case is (0, 0, -6) + assert_eq!(point, Some(Vector3::new(0.0, 0.0, -6.0))); + } + + #[test] + fn practice_problem_slide_158() { + let ray = Ray { + origin: Vector3::new(0.0, 0.0, 0.0), + direction: Vector3::new(0.0, 0.5, -1.0), + }; + let sphere = Sphere { + center: Vector3::new(0.0, 0.0, -10.0), + radius: 4.0, + }; + + // oops! In this case, the ray does not intersect the sphere. + assert_eq!(sphere.intersects_ray_at(&ray), None); + } +}