Support cylinders

assignment-1/examples/objects.txt

@@ -1,7 +1,7 @@
 imsize 640 480
-eye 0 0 3
+eye 0 0 15
 viewdir 0 0 -1
-hfov 120
+hfov 60
 updir 0 1 0
 bkgcolor 0.1 0.1 0.1
 
@@ -18,4 +18,4 @@ sphere 5 5 -1 1
 sphere -6 -4 -8 7
 
 mtlcolor 0.5 1 0.5
-cylinder 2 1 -2 1 -2 1 1 2
+cylinder 5 1 -2 1 -2 1 1 2

assignment-1/src/main.rs

@@ -90,7 +90,6 @@ 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

assignment-1/src/ray.rs

@@ -1,4 +1,5 @@
-use nalgebra::Vector3;
+use nalgebra::{Matrix3, Vector3};
+use ordered_float::NotNan;
 
 use crate::scene_data::{Cylinder, Sphere};
 
@@ -73,17 +74,120 @@ impl Ray {
   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();
+    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),
+      )
     };
 
-    // TODO: Implement
-    None
+    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)]

assignment-1/writeup.md

@@ -74,9 +74,9 @@ $\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.
+Because of the way I implemented parallel projection, it's recommended to
+either put the eye much farther back, or 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
 
@@ -85,4 +85,8 @@ cylinder, so that the cylinder location is $(0, 0, 0)$ and the direction vector
 is normalized into $(0, 0, 1)$.
 
 Then it's a matter of determining if the $x$ and $y$ coordinates fall into the
-space constrained by the equation $x^2 + y^2 = r^2$ and if $z \le L$.
+space constrained by the equation $(ox + t\times rx - cx)^2 + (oy + t\times ty -
+cy)^2 = r^2$ and if $z \le L$.
+
+See the comments in the code for a more detailed explanation of the equation
+used.