From adf8fa87faec0a5c8d0c28280287d424e2f0ed06 Mon Sep 17 00:00:00 2001 From: Michael Zhang Date: Thu, 2 Feb 2023 01:43:54 -0600 Subject: [PATCH] Support cylinders --- assignment-1/examples/objects.txt | 6 +- assignment-1/src/main.rs | 1 - assignment-1/src/ray.rs | 122 +++++++++++++++++++++++++++--- assignment-1/writeup.md | 12 ++- 4 files changed, 124 insertions(+), 17 deletions(-) diff --git a/assignment-1/examples/objects.txt b/assignment-1/examples/objects.txt index 38fa760..0dd8545 100644 --- a/assignment-1/examples/objects.txt +++ b/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 diff --git a/assignment-1/src/main.rs b/assignment-1/src/main.rs index cd928e7..d1dc9be 100644 --- a/assignment-1/src/main.rs +++ b/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 diff --git a/assignment-1/src/ray.rs b/assignment-1/src/ray.rs index ebfd40e..be6bb45 100644 --- a/assignment-1/src/ray.rs +++ b/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,19 +74,122 @@ impl Ray { pub fn intersects_cylinder_at(&self, cylinder: &Cylinder) -> Option { // 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 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), + ) }; - // 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, 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; diff --git a/assignment-1/writeup.md b/assignment-1/writeup.md index 3adbbf5..f9edae1 100644 --- a/assignment-1/writeup.md +++ b/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.