use anyhow::Result; use nalgebra::Vector3; use ordered_float::NotNan; use crate::ray::Ray; use crate::utils::compute_rotation_matrix; use super::{data::ObjectKind, illumination::IntersectionContext}; #[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, ) -> Result> { // 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)?; let inverse_rotation_matrix = rotation_matrix.try_inverse().ok_or_else(|| { anyhow!("Rotation matrix for some reason does not have an inverse?") })?; // 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 possible_side_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... _ => bail!("Invalid determinant value: {discriminant}"), }; // Filter out solutions that don't have a valid Z position. let side_solutions = possible_side_solutions.into_iter().filter_map(|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 if z.abs() > self.length / 2.0 { return None; } let time = NotNan::new(t).ok()?; // The point on the center of the cylinder that corresponds to the z-axis // point of the intersection let center_at_z = { let mut center_point = rotation_matrix * ray_point; center_point.x = rotated_cylinder_center.x; center_point.y = rotated_cylinder_center.y; inverse_rotation_matrix * center_point }; let normal = (ray_point - center_at_z).normalize(); Some(IntersectionContext { time, point: ray_point, normal, }) }); // 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; if r.z == 0.0 { Vec::new() // No solutions here } else { vec![ (-o.z + c.z + self.length / 2.0) / r.z, (-o.z + c.z - self.length / 2.0) / r.z, ] } }; let end_solutions = possible_z_intersections.into_iter().filter_map(|t| { let ray_point = ray.eval(t); let rotated_point = rotation_matrix * ray_point; // Filter out all the solutions where the intersection point does not lie // in the circle if rotated_point.x.powi(2) + rotated_point.y.powi(2) > self.radius.powi(2) { return None; } let normal_rotated = Vector3::new(0.0, 0.0, rotated_point.z - rotated_cylinder_center.z) .normalize(); let normal = inverse_rotation_matrix * normal_rotated; let time = NotNan::new(t).ok()?; Some(IntersectionContext { time, point: ray_point, normal, }) }); let solutions = side_solutions .into_iter() .chain(end_solutions.into_iter()) // Remove any t < 0, since that means it's behind the viewer and we // can't see it. .filter(|ctx| *ctx.time >= 0.0); // Return the minimum solution Ok(solutions.min_by_key(|ctx| ctx.time)) } } #[cfg(test)] mod tests { use nalgebra::Vector3; use crate::{ray::Ray, scene::data::ObjectKind}; use super::Cylinder; #[test] fn test_cylinder() { let cylinder = Cylinder { center: Vector3::new(0.0, 0.0, 0.0), direction: Vector3::new(0.0, 1.0, 0.0), radius: 3.0, length: 4.0, }; let eye = Vector3::new(0.0, 3.0, 3.0); let end = Vector3::new(0.0, 2.0, 2.0); let ray = Ray::from_endpoints(eye, end); let res = cylinder.intersects_ray_at(&ray); panic!("Result: {res:?}"); } }