WIP for cylinder

assignment-1/Cargo.lock

@@ -8,12 +8,22 @@ version = "1.0.68"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "2cb2f989d18dd141ab8ae82f64d1a8cdd37e0840f73a406896cf5e99502fab61"
 
+[[package]]
+name = "approx"
+version = "0.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6"
+dependencies = [
+ "num-traits",
+]
+
 [[package]]
 name = "assignment-1"
 version = "0.1.0"
 dependencies = [
  "anyhow",
  "clap",
+ "nalgebra",
  "num",
  "ordered-float",
  "rayon",
@@ -32,10 +42,16 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a"
 
 [[package]]
-name = "cc"
+name = "bytemuck"
-version = "1.0.78"
+version = "1.13.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a20104e2335ce8a659d6dd92a51a767a0c062599c73b343fd152cb401e828c3d"
+checksum = "c041d3eab048880cb0b86b256447da3f18859a163c3b8d8893f4e6368abe6393"
+
+[[package]]
+name = "cc"
+version = "1.0.79"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "50d30906286121d95be3d479533b458f87493b30a4b5f79a607db8f5d11aa91f"
 
 [[package]]
 name = "cfg-if"
@@ -152,9 +168,9 @@ dependencies = [
 
 [[package]]
 name = "heck"
-version = "0.4.0"
+version = "0.4.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2540771e65fc8cb83cd6e8a237f70c319bd5c29f78ed1084ba5d50eeac86f7f9"
+checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8"
 
 [[package]]
 name = "hermit-abi"
@@ -199,6 +215,15 @@ version = "0.1.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "f051f77a7c8e6957c0696eac88f26b0117e54f52d3fc682ab19397a8812846a4"
 
+[[package]]
+name = "matrixmultiply"
+version = "0.3.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "add85d4dd35074e6fedc608f8c8f513a3548619a9024b751949ef0e8e45a4d84"
+dependencies = [
+ "rawpointer",
+]
+
 [[package]]
 name = "memoffset"
 version = "0.7.1"
@@ -208,6 +233,33 @@ dependencies = [
  "autocfg",
 ]
 
+[[package]]
+name = "nalgebra"
+version = "0.32.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f6515c882ebfddccaa73ead7320ca28036c4bc84c9bcca3cc0cbba8efe89223a"
+dependencies = [
+ "approx",
+ "matrixmultiply",
+ "nalgebra-macros",
+ "num-complex",
+ "num-rational",
+ "num-traits",
+ "simba",
+ "typenum",
+]
+
+[[package]]
+name = "nalgebra-macros"
+version = "0.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d232c68884c0c99810a5a4d333ef7e47689cfd0edc85efc9e54e1e6bf5212766"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn",
+]
+
 [[package]]
 name = "num"
 version = "0.4.0"
@@ -318,6 +370,12 @@ version = "6.4.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "9b7820b9daea5457c9f21c69448905d723fbd21136ccf521748f23fd49e723ee"
 
+[[package]]
+name = "paste"
+version = "1.0.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d01a5bd0424d00070b0098dd17ebca6f961a959dead1dbcbbbc1d1cd8d3deeba"
+
 [[package]]
 name = "proc-macro-error"
 version = "1.0.4"
@@ -360,6 +418,12 @@ dependencies = [
  "proc-macro2",
 ]
 
+[[package]]
+name = "rawpointer"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
+
 [[package]]
 name = "rayon"
 version = "1.6.1"
@@ -396,6 +460,15 @@ dependencies = [
  "windows-sys",
 ]
 
+[[package]]
+name = "safe_arch"
+version = "0.6.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "794821e4ccb0d9f979512f9c1973480123f9bd62a90d74ab0f9426fcf8f4a529"
+dependencies = [
+ "bytemuck",
+]
+
 [[package]]
 name = "scopeguard"
 version = "1.1.0"
@@ -408,6 +481,19 @@ version = "1.0.152"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "bb7d1f0d3021d347a83e556fc4683dea2ea09d87bccdf88ff5c12545d89d5efb"
 
+[[package]]
+name = "simba"
+version = "0.8.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "50582927ed6f77e4ac020c057f37a268fc6aebc29225050365aacbb9deeeddc4"
+dependencies = [
+ "approx",
+ "num-complex",
+ "num-traits",
+ "paste",
+ "wide",
+]
+
 [[package]]
 name = "strsim"
 version = "0.10.0"
@@ -434,6 +520,12 @@ dependencies = [
  "winapi-util",
 ]
 
+[[package]]
+name = "typenum"
+version = "1.16.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "497961ef93d974e23eb6f433eb5fe1b7930b659f06d12dec6fc44a8f554c0bba"
+
 [[package]]
 name = "unicode-ident"
 version = "1.0.6"
@@ -446,6 +538,16 @@ version = "0.9.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "49874b5167b65d7193b8aba1567f5c7d93d001cafc34600cee003eda787e483f"
 
+[[package]]
+name = "wide"
+version = "0.7.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ae41ecad2489a1655c8ef8489444b0b113c0a0c795944a3572a0931cf7d2525c"
+dependencies = [
+ "bytemuck",
+ "safe_arch",
+]
+
 [[package]]
 name = "winapi"
 version = "0.3.9"

assignment-1/Cargo.toml

@@ -7,6 +7,7 @@ edition = "2021"
 [dependencies]
 anyhow = "1.0.68"
 clap = { version = "4.1.4", features = ["derive"] }
+nalgebra = "0.32.1"
 num = { version = "0.4.0", features = ["serde"] }
 ordered-float = "3.4.0"
 rayon = "1.6.1"

assignment-1/Makefile

@@ -40,4 +40,4 @@ writeup.pdf: writeup.md
 	$(PANDOC) -o $@ $<
 
 clean:
-	rm -f $(HANDIN) $(BINARY) $(WRITEUP)
+	rm -f $(HANDIN) $(BINARY) $(WRITEUP) $(EXAMPLES_PPM) $(EXAMPLES_PNG)

assignment-1/src/input_file.rs

@@ -1,11 +1,11 @@
 use std::{fs::File, io::Read, path::Path};
 
 use anyhow::Result;
+use nalgebra::Vector3;
 
 use crate::{
   image::Color,
   scene_data::{Cylinder, Object, ObjectKind, Scene, Sphere},
-  vec3::Vec3,
 };
 
 /// Parse the input file into a scene
@@ -50,7 +50,7 @@ pub fn parse_input_file(path: impl AsRef<Path>) -> Result<Scene> {
         if parts.len() < 3 {
           bail!("Vec3 requires 3 components.");
         }
-        Ok(Vec3::new(parts[start], parts[start + 1], parts[start + 2]))
+        Ok(Vector3::new(parts[start], parts[start + 1], parts[start + 2]))
       };
 
       let read_color = || {

assignment-1/src/main.rs

@@ -5,7 +5,6 @@ mod image;
 mod input_file;
 mod ray;
 mod scene_data;
-mod vec3;
 
 use std::fs::File;
 use std::path::PathBuf;
@@ -91,11 +90,12 @@ fn main() -> Result<()> {
       let pixel_in_space = translate_pixel(px, py);
 
       let ray_start = if scene.parallel_projection {
+        println!("Using parallel projection.");
         // For a parallel projection, we'll just take the view direction and
         // subtract it from the target point. This means every single
         // ray will be viewed from a point at infinity, rather than a single eye
         // position.
-        let n = scene.view_dir.unit();
+        let n = scene.view_dir.normalize();
         let view_dir = n * distance;
         pixel_in_space - view_dir
       } else {

assignment-1/src/ray.rs

@@ -1,5 +1,6 @@
+use nalgebra::Vector3;
+
 use crate::scene_data::{Cylinder, Sphere};
-use crate::vec3::Vec3;
 
 /// A normalized parametric Ray of the form (origin + direction * time)
 ///
@@ -7,14 +8,14 @@ use crate::vec3::Vec3;
 /// time occurs on the ray.
 #[derive(Debug)]
 pub struct Ray {
-  origin: Vec3,
+  origin: Vector3<f64>,
-  direction: Vec3,
+  direction: Vector3<f64>,
 }
 
 impl Ray {
   /// Construct a ray from endpoints
-  pub fn from_endpoints(start: Vec3, end: Vec3) -> Self {
+  pub fn from_endpoints(start: Vector3<f64>, end: Vector3<f64>) -> Self {
-    let delta = (end - start).unit();
+    let delta = (end - start).normalize();
     Ray {
       origin: start,
       direction: delta,
@@ -22,7 +23,7 @@ impl Ray {
   }
 
   /// Evaluate the ray at a certain point in time, yielding a point
-  pub fn eval(&self, time: f64) -> Vec3 {
+  pub fn eval(&self, time: f64) -> Vector3<f64> {
     self.origin + self.direction * time
   }
 
@@ -70,6 +71,16 @@ impl Ray {
   ///
   /// 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 rotation_matrix = {
+      let target_direction = Vector3::new(0.0, 0.0, 1.0);
+      let axis = target_direction.cross(&cylinder.direction).normalize();
+      let angle = (target_direction.dot(&cylinder.direction)
+        / (target_direction.norm() * cylinder.direction.norm()))
+      .acos();
+    };
+
     // TODO: Implement
     None
   }
@@ -77,36 +88,37 @@ impl Ray {
 
 #[cfg(test)]
 mod tests {
+  use nalgebra::Vector3;
+
   use crate::scene_data::Sphere;
-  use crate::vec3::Vec3;
 
   use super::Ray;
 
   #[test]
   fn practice_problem_slide_154() {
     let ray = Ray {
-      origin: Vec3::new(0.0, 0.0, 0.0),
+      origin: Vector3::new(0.0, 0.0, 0.0),
-      direction: Vec3::new(0.0, 0.0, -1.0),
+      direction: Vector3::new(0.0, 0.0, -1.0),
     };
     let sphere = Sphere {
-      center: Vec3::new(0.0, 0.0, -10.0),
+      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(Vec3::new(0.0, 0.0, -6.0)));
+    assert_eq!(point, Some(Vector3::new(0.0, 0.0, -6.0)));
   }
 
   #[test]
   fn practice_problem_slide_158() {
     let ray = Ray {
-      origin: Vec3::new(0.0, 0.0, 0.0),
+      origin: Vector3::new(0.0, 0.0, 0.0),
-      direction: Vec3::new(0.0, 0.5, -1.0),
+      direction: Vector3::new(0.0, 0.5, -1.0),
     };
     let sphere = Sphere {
-      center: Vec3::new(0.0, 0.0, -10.0),
+      center: Vector3::new(0.0, 0.0, -10.0),
       radius: 4.0,
     };
 

assignment-1/src/scene_data.rs

@@ -1,15 +1,17 @@
-use crate::{image::Color, vec3::Vec3};
+use nalgebra::Vector3;
+
+use crate::image::Color;
 
 #[derive(Debug)]
 pub struct Sphere {
-  pub center: Vec3,
+  pub center: Vector3<f64>,
   pub radius: f64,
 }
 
 #[derive(Debug)]
 pub struct Cylinder {
-  pub center: Vec3,
+  pub center: Vector3<f64>,
-  pub direction: Vec3,
+  pub direction: Vector3<f64>,
   pub radius: f64,
   pub length: f64,
 }
@@ -30,17 +32,17 @@ pub enum ObjectKind {
 
 #[derive(Debug)]
 pub struct Rect {
-  pub upper_left: Vec3,
+  pub upper_left: Vector3<f64>,
-  pub upper_right: Vec3,
+  pub upper_right: Vector3<f64>,
-  pub lower_left: Vec3,
+  pub lower_left: Vector3<f64>,
-  pub lower_right: Vec3,
+  pub lower_right: Vector3<f64>,
 }
 
 #[derive(Debug, Default)]
 pub struct Scene {
-  pub eye_pos: Vec3,
+  pub eye_pos: Vector3<f64>,
-  pub view_dir: Vec3,
+  pub view_dir: Vector3<f64>,
-  pub up_dir: Vec3,
+  pub up_dir: Vector3<f64>,
 
   /// Horizontal field of view (in degrees)
   pub hfov: f64,
@@ -60,8 +62,8 @@ impl Scene {
   /// Determine the boundaries of the viewing window in world coordinates
   pub fn compute_viewing_window(&self, distance: f64) -> Rect {
     // Compute viewing directions
-    let u = Vec3::cross(self.view_dir, self.up_dir).unit();
+    let u = self.view_dir.cross(&self.up_dir).normalize();
-    let v = Vec3::cross(u, self.view_dir).unit();
+    let v = u.cross(&self.view_dir).normalize();
 
     // Compute dimensions of viewing window based on field of view
     let viewing_width = {
@@ -80,7 +82,7 @@ impl Scene {
     let viewing_height = viewing_width / aspect_ratio;
 
     // Compute viewing window corners
-    let n = self.view_dir.unit();
+    let n = self.view_dir.normalize();
     #[rustfmt::skip] // Otherwise this line wraps over
     let view_window = Rect {
       upper_left: self.eye_pos + n * distance - u * (viewing_width / 2.0) + v * (viewing_height / 2.0),

assignment-1/src/vec3.rs

@@ -1,117 +0,0 @@
-use std::fmt;
-use std::ops::{Add, Div, Mul, Sub};
-
-use num::Float;
-
-#[derive(Copy, Clone, Default, Debug, PartialEq, Eq)]
-pub struct Vec3<T = f64> {
-  pub x: T,
-  pub y: T,
-  pub z: T,
-}
-
-impl<T> Vec3<T> {
-  pub fn new(x: T, y: T, z: T) -> Self {
-    Vec3 { x, y, z }
-  }
-}
-
-impl<T: fmt::Display> fmt::Display for Vec3<T> {
-  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-    f.write_str("(")?;
-    self.x.fmt(f)?;
-    f.write_str(", ")?;
-    self.y.fmt(f)?;
-    f.write_str(", ")?;
-    self.z.fmt(f)?;
-    f.write_str(")")?;
-
-    Ok(())
-  }
-}
-
-impl<T: Float> Vec3<T> {
-  /// Cross product on floats
-  pub fn cross(u: Self, v: Self) -> Self {
-    Vec3::new(
-      u.y * v.z - u.z * v.y,
-      u.z * v.x - u.x * v.z,
-      u.x * v.y - u.y * v.x,
-    )
-  }
-
-  /// Vector L2-norm
-  pub fn norm(&self) -> T {
-    (self.x.powi(2) + self.y.powi(2) + self.z.powi(2)).sqrt()
-  }
-
-  /// Normalize
-  pub fn unit(&self) -> Self {
-    let norm = self.norm();
-    Vec3::new(self.x / norm, self.y / norm, self.z / norm)
-  }
-}
-
-/// Vector addition
-impl<T> Add<Vec3<T>> for Vec3<T>
-where
-  T: Add<T>,
-{
-  type Output = Vec3<T::Output>;
-
-  fn add(self, rhs: Vec3<T>) -> Self::Output {
-    Vec3::new(self.x + rhs.x, self.y + rhs.y, self.z + rhs.z)
-  }
-}
-
-/// Vector subtraction
-impl<T> Sub<Vec3<T>> for Vec3<T>
-where
-  T: Sub<T>,
-{
-  type Output = Vec3<T::Output>;
-
-  fn sub(self, rhs: Vec3<T>) -> Self::Output {
-    Vec3::new(self.x - rhs.x, self.y - rhs.y, self.z - rhs.z)
-  }
-}
-
-/// Scalar multiplication
-impl<T, U> Mul<U> for Vec3<T>
-where
-  T: Mul<U, Output = U>,
-  U: Copy,
-{
-  type Output = Vec3<U>;
-
-  fn mul(self, rhs: U) -> Self::Output {
-    Vec3::new(self.x * rhs, self.y * rhs, self.z * rhs)
-  }
-}
-
-/// Scalar division
-impl<T, U> Div<U> for Vec3<T>
-where
-  T: Div<U, Output = U>,
-  U: Copy,
-{
-  type Output = Vec3<U>;
-
-  fn div(self, rhs: U) -> Self::Output {
-    Vec3::new(self.x / rhs, self.y / rhs, self.z / rhs)
-  }
-}
-
-#[cfg(test)]
-mod tests {
-  use super::Vec3;
-
-  #[test]
-  fn test_cross() {
-    let u = Vec3::new(1.0, 0.0, 0.0);
-    let v = Vec3::new(0.0, 1.0, 0.0);
-
-    assert_eq!(Vec3::cross(u, v), Vec3::new(0.0, 0.0, 1.0));
-    assert_eq!(Vec3::cross(v, u), Vec3::new(0.0, 0.0, -1.0));
-  }
-}

assignment-1/writeup.md

@@ -72,6 +72,12 @@ location for any pixel.
 (Technically really we would want the middle of the pixel, so just add
 $\frac{\Delta x + \Delta y}{2}$ to the point to get that)
 
+### Parallel Projection Notes
+
+Because of the way I implemented parallel projection, it's recommended to use
+`--distance` to force a much bigger distance from the eye for the raycaster.
+This is due to the size of the image.
+
 ### Cylinder Intersection Notes
 
 First, we will transform the current point into the vector space of the