running into some nan problems

assignment-1d/src/lib.rs

@@ -12,6 +12,9 @@ pub mod ray;
 pub mod scene;
 pub mod utils;
 
+// Creating a bunch of aliases here to make it more obvious which one I'm
+// expecting a variable to be
+
 pub type Point2 = Vector2<f64>;
 pub type Point = Vector3<f64>;
 pub type Vector = Vector3<f64>;

assignment-1d/src/main.rs

@@ -84,40 +84,7 @@ fn main() -> Result<()> {
 
       let ray = Ray::from_endpoints(ray_start, pixel_in_space);
 
-      let intersections = scene
-        .objects
-        .iter()
-        .enumerate()
-        .filter_map(|(i, object)| {
-          match object.kind.intersects_ray_at(&scene, &ray) {
-            Ok(Some(t)) => {
-              // Return both the t and the sphere, because we want to sort on
-              // the t but later retrieve attributes from the sphere
-              Some(Ok((i, t, object)))
-            }
-            Ok(None) => None,
-            Err(err) => {
-              error!("Error: {err}");
-              Some(Err(err))
-            }
-          }
-        })
-        .collect::<Result<Vec<_>>>()?;
-
-      // Sort the list of intersection times by the lowest one.
-      let earliest_intersection =
-        intersections.into_iter().min_by_key(|(_, t, _)| t.time);
-
-      Ok(match earliest_intersection {
-        // Take the object's material color
-        Some((obj_idx, intersection_context, object)) => {
-          scene.compute_pixel_color(obj_idx, object, intersection_context)?
-        }
-
-        // There was no intersection, so this should default to the scene's
-        // background color
-        None => scene.bkg_color,
-      })
+      scene.trace_single_ray(ray, 0)
     })
     .collect::<Result<Vec<_>>>()?;
 

assignment-1d/src/ray.rs

@@ -11,6 +11,10 @@ pub struct Ray {
 }
 
 impl Ray {
+  pub fn new(origin: Point, direction: Vector) -> Self {
+    Ray { origin, direction }
+  }
+
   /// Construct a ray from endpoints
   pub fn from_endpoints(start: Point, end: Point) -> Self {
     let delta = (end - start).normalize();

assignment-1d/src/scene/cylinder.rs

@@ -116,6 +116,7 @@ impl Cylinder {
         time,
         point: ray_point,
         normal,
+        exiting: todo!(),
       })
     });
 
@@ -165,6 +166,7 @@ impl Cylinder {
         time,
         point: ray_point,
         normal,
+        exiting: todo!(),
       })
     });
 

assignment-1d/src/scene/data.rs

@@ -23,6 +23,12 @@ pub struct Material {
   pub k_d: f64,
   pub k_s: f64,
   pub exponent: f64,
+
+  /// Opacity
+  pub alpha: f64,
+
+  /// Index of refraction
+  pub eta: f64,
 }
 
 #[derive(Debug)]

assignment-1d/src/scene/illumination.rs

@@ -11,11 +11,14 @@ use rayon::prelude::{
 use crate::{image::Color, ray::Ray, utils::dot, Point, Vector};
 
 use super::{
-  data::{DepthCueing, Light, LightKind},
+  data::{DepthCueing, Light, LightKind, Material},
   object::Object,
   Scene,
 };
 
+// TODO: Is this a good constant?
+const JITTER_CONST: f64 = 0.05;
+
 impl Scene {
   /// Determine the color that should be used to fill this pixel.
   ///
@@ -28,7 +31,9 @@ impl Scene {
     &self,
     obj_idx: usize,
     object: &Object,
+    incident_ray: Ray,
     intersection_context: IntersectionContext,
+    depth: usize,
   ) -> Result<Color> {
     let material = match self.materials.get(object.material_idx) {
       Some(v) => v,
@@ -110,7 +115,73 @@ impl Scene {
       })
       .sum();
 
-    let color = ambient_component + diffuse_and_specular;
+    let specular_reflection: Color = {
+      let reflection_ray = self.compute_reflection_ray(
+        incident_ray.direction,
+        intersection_context.normal,
+      );
+
+      let fresnel_coefficient = self.compute_fresnel_coefficient(
+        &material,
+        incident_ray.direction,
+        intersection_context.normal,
+      );
+
+      // Jitter a bit to reduce acne
+      let origin = intersection_context.point;
+      let origin = origin + JITTER_CONST * reflection_ray;
+
+      let ray = Ray::new(origin, reflection_ray);
+      let r_lambda = self.trace_single_ray(ray, depth + 1)?;
+
+      fresnel_coefficient * r_lambda
+    };
+
+    let transparency = {
+      let n = intersection_context.normal;
+      let i = incident_ray.direction;
+
+      let (eta_i, eta_t) = match intersection_context.exiting {
+        true => (material.eta, 1.0),
+        false => (1.0, material.eta),
+      };
+
+      // Slide 69
+      let cos_theta_i = dot(i, n);
+      assert!(cos_theta_i != std::f64::NAN);
+      let sin_theta_i = (1.0 - cos_theta_i.powi(2)).sqrt();
+      assert!(sin_theta_i != std::f64::NAN);
+      let sin_theta_t = (eta_i / eta_t) * sin_theta_i;
+      assert!(sin_theta_t != std::f64::NAN);
+      let cos_theta_t = (1.0 - sin_theta_t.powi(2)).sqrt();
+      assert!(cos_theta_t != std::f64::NAN);
+
+      let fresnel_coefficient =
+        self.compute_fresnel_coefficient(&material, i, n);
+
+      // Calculate refraction direction
+      let a = (-n).normalize() * cos_theta_t;
+      let s_direction = cos_theta_i * n - i;
+      let m_unit = s_direction.normalize();
+      let b = m_unit * sin_theta_t;
+      let t = a + b;
+
+      // Jitter a bit to reduce acne
+      // TODO: Is this a good constant?
+      let origin = intersection_context.point;
+      let origin = origin + JITTER_CONST * t;
+
+      let ray = Ray::new(origin, t);
+      let t_lambda = self.trace_single_ray(ray, depth + 1)?;
+
+      (1.0 - fresnel_coefficient) * (1.0 - material.alpha) * t_lambda
+    };
+
+    // This is the result of the Phong illumination equation.
+    let color = ambient_component
+      + diffuse_and_specular
+      + specular_reflection
+      + transparency;
 
     // Apply depth cueing to the result
     let a_dc = {
@@ -265,6 +336,34 @@ impl Scene {
 
     (JITTER_RAYS - num_obstructed_rays) as f64 / JITTER_RAYS as f64
   }
+
+  fn compute_fresnel_coefficient(
+    &self,
+    material: &Material,
+    incident_ray: Vector,
+    normal: Vector,
+  ) -> f64 {
+    let cos_theta_i = dot(incident_ray, normal);
+
+    let f0 = ((material.eta - 1.0) / (material.eta + 1.0)).powi(2);
+    let fr = f0 * 1.0 + (1.0 - f0) * (1.0 - cos_theta_i).powi(5);
+
+    fr
+  }
+
+  fn compute_reflection_ray(
+    &self,
+    incident_ray: Vector,
+    normal: Vector,
+  ) -> Vector {
+    let opposite_incident_ray = (-incident_ray).normalize();
+    let unit_normal = normal.normalize();
+
+    let a = dot(unit_normal, opposite_incident_ray);
+    let r = 2.0 * (a * unit_normal) - opposite_incident_ray;
+
+    r
+  }
 }
 
 /// Information about an intersection
@@ -287,6 +386,9 @@ pub struct IntersectionContext {
   /// intersection point
   #[derivative(PartialEq = "ignore", Ord = "ignore")]
   pub normal: Vector,
+
+  /// Is this ray exiting the material at the intersection point?
+  pub exiting: bool,
 }
 
 impl Eq for IntersectionContext {}

assignment-1d/src/scene/input_file/mod.rs

@@ -1,3 +1,5 @@
+pub mod triangle_vertex;
+
 use std::fs::File;
 use std::io::Read;
 use std::path::Path;
@@ -11,9 +13,10 @@ use crate::{
   scene::{
     cylinder::Cylinder,
     data::{Attenuation, Light, LightKind, Material},
+    input_file::triangle_vertex::TriangleVertex,
     object::{Object, ObjectKind},
     sphere::Sphere,
-    texture::{Texture, NormalMap},
+    texture::{NormalMap, Texture},
     triangle::Triangle,
     Scene,
   },
@@ -149,7 +152,7 @@ impl Scene {
           };
         }
 
-        // mtlcolor Odr Odg Odb Osr Osg Osb ka kd ks n
+        // mtlcolor Odr Odg Odb Osr Osg Osb ka kd ks n alpha eta
         "mtlcolor" => {
           let diffuse_color = r!(Color);
           let specular_color = r!(Color);
@@ -157,6 +160,8 @@ impl Scene {
           let k_d = r!(f64);
           let k_s = r!(f64);
           let exponent = r!(f64);
+          let alpha = r!(f64);
+          let eta = r!(f64);
 
           let material = Material {
             diffuse_color,
@@ -165,6 +170,8 @@ impl Scene {
             k_d,
             k_s,
             exponent,
+            alpha,
+            eta,
           };
 
           let idx = scene.materials.len();
@@ -351,48 +358,3 @@ impl Construct for Vector3<f64> {
     Ok(Vector3::new(x, y, z))
   }
 }
-
-#[derive(Debug, Copy, Clone, PartialEq)]
-struct TriangleVertex {
-  vertex_idx: usize,
-  normal_idx: Option<usize>,
-  texture_idx: Option<usize>,
-}
-
-impl Construct for TriangleVertex {
-  type Args = ();
-
-  fn construct<'a, I>(it: &mut I, _: Self::Args) -> Result<Self>
-  where
-    I: Iterator<Item = &'a str>,
-  {
-    let s = match it.next() {
-      Some(v) => v,
-      None => bail!("Waiting on another"),
-    };
-
-    // Note: indexed by 1 not 0, so we will just do the subtraction
-    // here to avoid having to deal with it later
-    let parts = s.split("/").collect_vec();
-    ensure!(parts.len() >= 1 && parts.len() <= 3);
-    let vertex_idx: usize = parts[0].parse::<usize>()? - 1;
-
-    let texture_idx =
-      match parts.get(1).and_then(|s| (!s.is_empty()).then(|| *s)) {
-        Some(s) => Some(s.parse::<usize>()? - 1),
-        None => None,
-      };
-
-    let normal_idx =
-      match parts.get(2).and_then(|s| (!s.is_empty()).then(|| *s)) {
-        Some(s) => Some(s.parse::<usize>()? - 1),
-        None => None,
-      };
-
-    Ok(TriangleVertex {
-      vertex_idx,
-      texture_idx,
-      normal_idx,
-    })
-  }
-}

assignment-1d/src/scene/input_file/triangle_vertex.rs

@@ -0,0 +1,49 @@
+use anyhow::Result;
+use itertools::Itertools;
+
+use super::Construct;
+
+#[derive(Debug, Copy, Clone, PartialEq)]
+pub struct TriangleVertex {
+  pub vertex_idx: usize,
+  pub normal_idx: Option<usize>,
+  pub texture_idx: Option<usize>,
+}
+
+impl Construct for TriangleVertex {
+  type Args = ();
+
+  fn construct<'a, I>(it: &mut I, _: Self::Args) -> Result<Self>
+  where
+    I: Iterator<Item = &'a str>,
+  {
+    let s = match it.next() {
+      Some(v) => v,
+      None => bail!("Waiting on another"),
+    };
+
+    // Note: indexed by 1 not 0, so we will just do the subtraction
+    // here to avoid having to deal with it later
+    let parts = s.split("/").collect_vec();
+    ensure!(parts.len() >= 1 && parts.len() <= 3);
+    let vertex_idx: usize = parts[0].parse::<usize>()? - 1;
+
+    let texture_idx =
+      match parts.get(1).and_then(|s| (!s.is_empty()).then(|| *s)) {
+        Some(s) => Some(s.parse::<usize>()? - 1),
+        None => None,
+      };
+
+    let normal_idx =
+      match parts.get(2).and_then(|s| (!s.is_empty()).then(|| *s)) {
+        Some(s) => Some(s.parse::<usize>()? - 1),
+        None => None,
+      };
+
+    Ok(TriangleVertex {
+      vertex_idx,
+      texture_idx,
+      normal_idx,
+    })
+  }
+}

assignment-1d/src/scene/materials.rs

@@ -0,0 +1,6 @@
+use super::data::Material;
+
+pub struct MaterialStack<'a> {
+  stack: Vec<&'a f64>,
+}
+

assignment-1d/src/scene/mod.rs

@@ -6,6 +6,8 @@ pub mod object;
 pub mod sphere;
 pub mod texture;
 pub mod triangle;
+pub mod tracing;
+pub mod materials;
 
 use crate::image::Color;
 use crate::{Point, Point2, Vector};

assignment-1d/src/scene/sphere.rs

@@ -33,6 +33,7 @@ impl Sphere {
       + (ray.origin.y - self.center.y).powi(2)
       + (ray.origin.z - self.center.z).powi(2)
       - self.radius.powi(2);
+
     let discriminant = b * b - 4.0 * a * c;
 
     let time = match discriminant {
@@ -68,10 +69,22 @@ impl Sphere {
     let point = ray.eval(*time);
     let normal = (point - self.center).normalize();
 
+    let exiting = {
+      // To figure out if we're exiting, just test if the origin is inside the
+      // sphere
+
+      let dx = ray.origin.x - self.center.x;
+      let dy = ray.origin.y - self.center.y;
+      let dz = ray.origin.z - self.center.z;
+
+      dx.powi(2) + dy.powi(2) + dz.powi(2) < self.radius.powi(2)
+    };
+
     Ok(Some(IntersectionContext {
       time,
       point,
       normal,
+      exiting,
     }))
   }
 

assignment-1d/src/scene/tracing.rs

@@ -0,0 +1,55 @@
+use anyhow::Result;
+
+use crate::{image::Color, ray::Ray};
+
+use super::Scene;
+
+const MAX_RECURSION_DEPTH: usize = 10_usize;
+
+impl Scene {
+  pub fn trace_single_ray(&self, ray: Ray, depth: usize) -> Result<Color> {
+    if depth > MAX_RECURSION_DEPTH {
+      return Ok(Color::new(0.0, 0.0, 0.0));
+    }
+
+    let intersections = self
+      .objects
+      .iter()
+      .enumerate()
+      .filter_map(|(i, object)| {
+        match object.kind.intersects_ray_at(&self, &ray) {
+          Ok(Some(t)) => {
+            // Return both the t and the sphere, because we want to sort on
+            // the t but later retrieve attributes from the sphere
+            Some(Ok((i, t, object)))
+          }
+          Ok(None) => None,
+          Err(err) => {
+            error!("Error: {err}");
+            Some(Err(err))
+          }
+        }
+      })
+      .collect::<Result<Vec<_>>>()?;
+
+    // Sort the list of intersection times by the lowest one.
+    let earliest_intersection =
+      intersections.into_iter().min_by_key(|(_, t, _)| t.time);
+
+    Ok(match earliest_intersection {
+      // Take the object's material color
+      Some((obj_idx, intersection_context, object)) => self
+        .compute_pixel_color(
+          obj_idx,
+          object,
+          ray,
+          intersection_context,
+          depth,
+        )?,
+
+      // There was no intersection, so this should default to the scene's
+      // background color
+      None => self.bkg_color,
+    })
+  }
+}

assignment-1d/src/scene/triangle.rs

@@ -100,6 +100,7 @@ impl Triangle {
       time,
       point,
       normal,
+      exiting: todo!(),
     }))
   }