csci5607/assignment-1b/src/scene/cylinder.rs

use anyhow::Result;
use nalgebra::Vector3;
use ordered_float::NotNan;

use crate::ray::Ray;
use crate::utils::compute_rotation_matrix;

use super::data::{IntersectionContext, 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,
  ) -> Result<Option<IntersectionContext>> {
    // 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 {
        // No solutions here
        vec![]
      } 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))
  }
}