Refactored large part of the code to make it more extensible

This commit is contained in:
Michael Zhang 2023-02-02 02:08:17 -06:00
parent adf8fa87fa
commit a7b9ec249a
8 changed files with 283 additions and 244 deletions

View file

@ -5,7 +5,11 @@ use nalgebra::Vector3;
use crate::{ use crate::{
image::Color, image::Color,
scene_data::{Cylinder, Object, ObjectKind, Scene, Sphere}, scene::{
cylinder::Cylinder,
data::{Object, Scene},
sphere::Sphere,
},
}; };
/// Parse the input file into a scene /// Parse the input file into a scene
@ -50,7 +54,11 @@ pub fn parse_input_file(path: impl AsRef<Path>) -> Result<Scene> {
if parts.len() < 3 { if parts.len() < 3 {
bail!("Vec3 requires 3 components."); 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 = || { let read_color = || {
@ -75,7 +83,7 @@ pub fn parse_input_file(path: impl AsRef<Path>) -> Result<Scene> {
} }
"sphere" => scene.objects.push(Object { "sphere" => scene.objects.push(Object {
kind: ObjectKind::Sphere(Sphere { kind: Box::new(Sphere {
center: read_vec3(0)?, center: read_vec3(0)?,
radius: parts[3], radius: parts[3],
}), }),
@ -86,7 +94,7 @@ pub fn parse_input_file(path: impl AsRef<Path>) -> Result<Scene> {
}), }),
"cylinder" => scene.objects.push(Object { "cylinder" => scene.objects.push(Object {
kind: ObjectKind::Cylinder(Cylinder { kind: Box::new(Cylinder {
center: read_vec3(0)?, center: read_vec3(0)?,
direction: read_vec3(3)?, direction: read_vec3(3)?,
radius: parts[6], radius: parts[6],

View file

@ -3,8 +3,9 @@ extern crate anyhow;
mod image; mod image;
mod input_file; mod input_file;
mod math;
mod ray; mod ray;
mod scene_data; mod scene;
use std::fs::File; use std::fs::File;
use std::path::PathBuf; use std::path::PathBuf;
@ -13,7 +14,6 @@ use anyhow::Result;
use clap::Parser; use clap::Parser;
use ordered_float::NotNan; use ordered_float::NotNan;
use rayon::prelude::{IntoParallelIterator, ParallelIterator}; use rayon::prelude::{IntoParallelIterator, ParallelIterator};
use scene_data::ObjectKind;
use crate::image::Image; use crate::image::Image;
use crate::input_file::parse_input_file; use crate::input_file::parse_input_file;
@ -107,11 +107,7 @@ fn main() -> Result<()> {
.objects .objects
.iter() .iter()
.filter_map(|object| { .filter_map(|object| {
use ObjectKind::*; let intersection_point_opt = object.kind.intersects_ray_at(&ray);
let intersection_point_opt = match &object.kind {
Sphere(sphere) => ray.intersects_sphere_at(&sphere),
Cylinder(cylinder) => ray.intersects_cylinder_at(&cylinder),
};
intersection_point_opt.and_then(|t| { intersection_point_opt.and_then(|t| {
// Unfortunately, IEEE floats in Rust don't have total ordering, // 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. // Sort the list of intersection times by the lowest one.
.min_by_key(|(t, _)| *t); .min_by_key(|(t, _)| *t);
let pixel_color = match earliest_intersection { match earliest_intersection {
// Take the object's material color // 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 // There was no intersection, so this should default to the scene's
// background color // background color
None => scene.bkg_color, None => scene.bkg_color,
}; }
pixel_color
}) })
.collect::<Vec<_>>(); .collect::<Vec<_>>();
@ -149,7 +143,7 @@ fn main() -> Result<()> {
}; };
{ {
let file = File::create(&out_file)?; let file = File::create(out_file)?;
image.write(file)?; image.write(file)?;
} }

21
assignment-1/src/math.rs Normal file
View file

@ -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<f64>,
b: Vector3<f64>,
) -> Matrix3<f64> {
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
}

View file

@ -1,7 +1,4 @@
use nalgebra::{Matrix3, Vector3}; use nalgebra::Vector3;
use ordered_float::NotNan;
use crate::scene_data::{Cylinder, Sphere};
/// A normalized parametric Ray of the form (origin + direction * time) /// 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. /// time occurs on the ray.
#[derive(Debug)] #[derive(Debug)]
pub struct Ray { pub struct Ray {
origin: Vector3<f64>, pub origin: Vector3<f64>,
direction: Vector3<f64>, pub direction: Vector3<f64>,
} }
impl Ray { impl Ray {
@ -27,206 +24,4 @@ impl Ray {
pub fn eval(&self, time: f64) -> Vector3<f64> { pub fn eval(&self, time: f64) -> Vector3<f64> {
self.origin + self.direction * time 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<f64> {
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<f64> {
// 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<f64>, b: Vector3<f64>) -> Matrix3<f64> {
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);
}
} }

View file

@ -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<f64>,
pub direction: Vector3<f64>,
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<f64> {
// 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())
}
}

View file

@ -1,35 +1,27 @@
use std::fmt::Debug;
use nalgebra::Vector3; use nalgebra::Vector3;
use crate::image::Color; use crate::image::Color;
use crate::ray::Ray;
#[derive(Debug)] pub trait ObjectKind: Debug + Send + Sync {
pub struct Sphere { /// Determine where the ray intersects this object, returning the earliest
pub center: Vector3<f64>, /// time this happens. Returns None if no intersection occurs.
pub radius: f64, ///
} /// Also known as Trace_Ray in the slides, except not the part where it calls
/// Shade_Ray.
#[derive(Debug)] fn intersects_ray_at(&self, ray: &Ray) -> Option<f64>;
pub struct Cylinder {
pub center: Vector3<f64>,
pub direction: Vector3<f64>,
pub radius: f64,
pub length: f64,
} }
#[derive(Debug)] #[derive(Debug)]
pub struct Object { pub struct Object {
pub kind: ObjectKind, pub kind: Box<dyn ObjectKind>,
/// Index into the scene's material color list /// Index into the scene's material color list
pub material: usize, pub material: usize,
} }
#[derive(Debug)]
pub enum ObjectKind {
Sphere(Sphere),
Cylinder(Cylinder),
}
#[derive(Debug)] #[derive(Debug)]
pub struct Rect { pub struct Rect {
pub upper_left: Vector3<f64>, pub upper_left: Vector3<f64>,
@ -59,6 +51,19 @@ pub struct Scene {
} }
impl 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 /// Determine the boundaries of the viewing window in world coordinates
pub fn compute_viewing_window(&self, distance: f64) -> Rect { pub fn compute_viewing_window(&self, distance: f64) -> Rect {
// Compute viewing directions // Compute viewing directions

View file

@ -0,0 +1,3 @@
pub mod data;
pub mod sphere;
pub mod cylinder;

View file

@ -0,0 +1,92 @@
use nalgebra::Vector3;
use crate::ray::Ray;
use super::data::ObjectKind;
#[derive(Debug)]
pub struct Sphere {
pub center: Vector3<f64>,
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<f64> {
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);
}
}