Replace nalgebra matrix and vector classes with own classes

Nalgebra only used to find matrix inverse now.

benches/simple_scene.rs

@@ -2,14 +2,13 @@ use criterion::{criterion_group, criterion_main, Criterion};
 
 use vanrijn::colour::{ColourRgbF, NamedColour};
 use vanrijn::materials::ReflectiveMaterial;
+use vanrijn::math::Vec3;
 use vanrijn::mesh::load_obj;
 use vanrijn::partial_render_scene;
 use vanrijn::raycasting::BoundingVolumeHierarchy;
 use vanrijn::scene::Scene;
 use vanrijn::util::Tile;
 
-use nalgebra::Point3;
-
 use std::path::Path;
 use std::sync::Arc;
 
@@ -22,7 +21,7 @@ fn simple_scene(bencher: &mut Criterion) {
 
     bencher.bench_function("simple_scene", |b| {
         let scene = Scene {
-            camera_location: Point3::new(-2.0, 1.0, -5.0),
+            camera_location: Vec3::new(-2.0, 1.0, -5.0),
             objects: vec![Box::new(BoundingVolumeHierarchy::build(
                 load_obj(
                     &model_file_path,

src/camera.rs

@@ -1,4 +1,4 @@
-use nalgebra::{convert, Point3, Vector3};
+use crate::math::Vec3;
 
 use super::colour::{ColourRgbF, NamedColour};
 use super::image::ImageRgbF;
@@ -14,13 +14,12 @@ struct ImageSampler {
 
     film_width: f64,
     film_height: f64,
-
+    camera_location: Vec3,
-    camera_location: Point3<f64>,
     film_distance: f64,
 }
 
 impl ImageSampler {
-    pub fn new(width: usize, height: usize, camera_location: Point3<f64>) -> ImageSampler {
+    pub fn new(width: usize, height: usize, camera_location: Vec3) -> ImageSampler {
         let (film_width, film_height) = {
             let width = width as f64;
             let height = height as f64;
@@ -34,7 +33,7 @@ impl ImageSampler {
         ImageSampler {
             image_height_pixels: height,
             image_width_pixels: width,
-            film_distance: convert(1.0),
+            film_distance: 1.0,
             film_width,
             film_height,
             camera_location,
@@ -51,7 +50,7 @@ impl ImageSampler {
     fn ray_for_pixel(&self, row: usize, column: usize) -> Ray {
         Ray::new(
             self.camera_location,
-            Vector3::new(
+            Vec3::new(
                 Self::scale(column, self.image_width_pixels, self.film_width)
                     - self.film_width * 0.5,
                 Self::scale(row, self.image_height_pixels, self.film_height)
@@ -75,11 +74,11 @@ const RECURSION_LIMIT: u16 = 32;
 /// # Examples
 //
 /// ```
-/// # use nalgebra::Point3;
+/// # use vanrijn::math::Vec3;
 /// # use vanrijn::scene::Scene;
 /// # use vanrijn::util::TileIterator;
 /// # use vanrijn::partial_render_scene;
-/// # let scene = Scene { camera_location: Point3::new(0.0, 0.0, 0.0), objects: vec![] };
+/// # let scene = Scene { camera_location: Vec3::new(0.0, 0.0, 0.0), objects: vec![] };
 /// let image_width = 640;
 /// let image_height = 480;
 /// let time_size = 32;
@@ -99,15 +98,15 @@ pub fn partial_render_scene(scene: &Scene, tile: Tile, height: usize, width: usi
         ambient_light: ColourRgbF::from_named(NamedColour::White) * ambient_intensity,
         lights: vec![
             DirectionalLight {
-                direction: Vector3::new(convert(1.0), convert(1.0), convert(-1.0)).normalize(),
+                direction: Vec3::new(1.0, 1.0, -1.0).normalize(),
                 colour: ColourRgbF::from_named(NamedColour::White) * directional_intensity1,
             },
             DirectionalLight {
-                direction: Vector3::new(convert(-0.5), convert(2.0), convert(-0.5)).normalize(),
+                direction: Vec3::new(-0.5, 2.0, -0.5).normalize(),
                 colour: ColourRgbF::from_named(NamedColour::White) * directional_intensity2,
             },
             DirectionalLight {
-                direction: Vector3::new(convert(-3.0), convert(0.1), convert(-0.5)).normalize(),
+                direction: Vec3::new(-3.0, 0.1, -0.5).normalize(),
                 colour: ColourRgbF::from_named(NamedColour::White) * directional_intensity3,
             },
         ],
@@ -154,10 +153,10 @@ mod tests {
 
         #[test]
         fn ray_for_pixel_returns_value_that_intersects_film_plane_at_expected_location() {
-            let target = ImageSampler::new(800, 600, Point3::new(0.0, 0.0, 0.0));
+            let target = ImageSampler::new(800, 600, Vec3::new(0.0, 0.0, 0.0));
             let ray = target.ray_for_pixel(100, 200);
             let film_plane = Plane::new(
-                Vector3::new(0.0, 0.0, 1.0),
+                Vec3::new(0.0, 0.0, 1.0),
                 target.film_distance,
                 Arc::new(LambertianMaterial::new_dummy()),
             );
@@ -175,11 +174,10 @@ mod tests {
             };
             let expected_x: f64 =
                 ImageSampler::scale(200, 800, target.film_width) - target.film_width * 0.5;
-            print!("{}, {}", expected_x, point_on_film_plane);
+            assert!((point_on_film_plane.x() - expected_x).abs() < 0.0000000001);
-            assert!((point_on_film_plane.x - expected_x).abs() < 0.0000000001);
             let expected_y =
                 ImageSampler::scale(100, 600, target.film_height) - target.film_height * 0.5;
-            assert!((point_on_film_plane.y - expected_y).abs() < 0.0000000001);
+            assert!((point_on_film_plane.y() - expected_y).abs() < 0.0000000001);
         }
     }
 }

src/colour.rs

@@ -1,16 +1,16 @@
-use nalgebra::Vector3;
+use crate::math::Vec3;
 
 use std::ops::{Add, Mul};
 
 #[derive(Copy, Clone, Debug)]
 pub struct ColourRgbF {
-    values: Vector3<f64>,
+    values: Vec3,
 }
 
 impl ColourRgbF {
     pub fn new(red: f64, green: f64, blue: f64) -> ColourRgbF {
         ColourRgbF {
-            values: Vector3::new(red, green, blue),
+            values: Vec3::new(red, green, blue),
         }
     }
 
@@ -34,23 +34,23 @@ impl ColourRgbF {
         }
     }
 
-    pub fn from_vector3(v: &Vector3<f64>) -> ColourRgbF {
+    pub fn from_vec3(v: &Vec3) -> ColourRgbF {
         ColourRgbF { values: *v }
     }
 
     pub fn red(&self) -> f64 {
-        self.values[0]
+        self.values.x()
     }
 
     pub fn green(&self) -> f64 {
-        self.values[1]
+        self.values.y()
     }
 
     pub fn blue(&self) -> f64 {
-        self.values[2]
+        self.values.z()
     }
 
-    pub fn as_vector3(&self) -> &Vector3<f64> {
+    pub fn as_vec3(&self) -> &Vec3 {
         &self.values
     }
 }
@@ -114,7 +114,7 @@ mod tests {
         use quickcheck_macros::quickcheck;
         impl Arbitrary for ColourRgbF {
             fn arbitrary<G: Gen>(g: &mut G) -> ColourRgbF {
-                let values = <Vector3<f64> as Arbitrary>::arbitrary(g);
+                let values = <Vec3 as Arbitrary>::arbitrary(g);
                 ColourRgbF { values }
             }
         }
@@ -130,8 +130,8 @@ mod tests {
         #[test]
         fn as_vector3_returns_expected_vector() {
             let target = ColourRgbF::new(1.0, 2.0, 3.0);
-            let result = target.as_vector3();
+            let result = target.as_vec3();
-            assert!(result.x == 1.0);
+            assert!(result.x() == 1.0);
         }
 
         #[quickcheck]

src/image.rs

@@ -3,9 +3,8 @@ use std::fs::File;
 use std::io::BufWriter;
 use std::path::Path;
 
-use nalgebra::{clamp, convert, Vector3};
+use crate::colour::{ColourRgbF, ColourRgbU8};
-
+use crate::math::Vec3;
-use super::colour::{ColourRgbF, ColourRgbU8};
 
 pub struct ImageRgbU8 {
     pixel_data: Vec<u8>,
@@ -104,7 +103,7 @@ impl ImageRgbF {
         ImageRgbF {
             width,
             height,
-            pixel_data: vec![convert(0.0); (width * height * 3) as usize],
+            pixel_data: vec![0.0; width * height * 3 as usize],
         }
     }
 
@@ -118,13 +117,13 @@ impl ImageRgbF {
     pub fn get_colour(&self, row: usize, column: usize) -> ColourRgbF {
         assert!(row < self.height && column < self.width);
         let index = self.calculate_index(row, column);
-        ColourRgbF::from_vector3(&Vector3::from_row_slice(&self.pixel_data[index..index + 3]))
+        ColourRgbF::from_vec3(&Vec3::from_slice(&self.pixel_data[index..index + 3]))
     }
 
     pub fn set_colour(&mut self, row: usize, column: usize, colour: ColourRgbF) {
         assert!(row < self.height && column < self.width);
         let index = self.calculate_index(row, column);
-        self.pixel_data[index..index + 3].copy_from_slice(&colour.as_vector3().as_slice());
+        self.pixel_data[index..index + 3].copy_from_slice(&colour.as_vec3().as_slice());
     }
 
     pub fn get_pixel_data(&self) -> &Vec<f64> {
@@ -183,7 +182,7 @@ pub struct ClampingToneMapper {}
 
 impl ClampingToneMapper {
     fn clamp(v: &f64) -> u8 {
-        clamp(v, &0.0, &1.0).normalized_to_byte()
+        v.clamp(0.0, 1.0).normalized_to_byte()
     }
 }
 

src/integrators.rs

@@ -1,4 +1,4 @@
-use nalgebra::{convert, Vector3};
+use crate::math::Vec3;
 
 use super::colour::ColourRgbF;
 use super::raycasting::{IntersectionInfo, Ray};
@@ -15,7 +15,7 @@ pub trait Integrator {
 }
 
 pub struct DirectionalLight {
-    pub direction: Vector3<f64>,
+    pub direction: Vec3,
     pub colour: ColourRgbF,
 }
 
@@ -42,9 +42,7 @@ impl Integrator for WhittedIntegrator {
         self.lights
             .iter()
             .map(|light| {
-                match sampler
+                match sampler.sample(&Ray::new(info.location, light.direction).bias(0.000_000_1)) {
-                    .sample(&Ray::new(info.location, light.direction).bias(convert(0.000_000_1)))
-                {
                     Some(_) => self.ambient_light,
                     None => {
                         info.material.bsdf()(
@@ -62,8 +60,7 @@ impl Integrator for WhittedIntegrator {
                     .map(|direction| {
                         let world_space_direction = bsdf_to_world_space * direction;
                         match sampler.sample(
-                            &Ray::new(info.location, world_space_direction)
+                            &Ray::new(info.location, world_space_direction).bias(0.000_000_1),
-                                .bias(convert(0.000_000_1)),
                         ) {
                             Some(recursive_hit) => {
                                 if recursion_limit > 0 {

src/lib.rs

@@ -1,4 +1,4 @@
-#![feature(external_doc)]
+#![feature(external_doc, clamp)]
 #![doc(include = "../README.md")]
 
 mod camera;

src/main.rs

@@ -7,8 +7,6 @@ use sdl2::rect::Rect;
 use sdl2::render::{Canvas, Texture};
 use sdl2::Sdl;
 
-use nalgebra::{Point3, Vector3};
-
 use clap::Arg;
 
 use std::path::{Path, PathBuf};
@@ -18,6 +16,7 @@ use std::time::Duration;
 use vanrijn::colour::{ColourRgbF, NamedColour};
 use vanrijn::image::{ClampingToneMapper, ImageRgbU8, ToneMapper};
 use vanrijn::materials::{LambertianMaterial, PhongMaterial, ReflectiveMaterial};
+use vanrijn::math::Vec3;
 use vanrijn::mesh::load_obj;
 use vanrijn::partial_render_scene;
 use vanrijn::raycasting::{Aggregate, BoundingVolumeHierarchy, Plane, Primitive, Sphere};
@@ -143,11 +142,11 @@ pub fn main() -> Result<(), Box<dyn std::error::Error>> {
     println!("Constructing Scene...");
 
     let scene = Scene {
-        camera_location: Point3::new(-2.0, 1.0, -5.0),
+        camera_location: Vec3::new(-2.0, 1.0, -5.0),
         objects: vec![
             Box::new(vec![
                 Box::new(Plane::new(
-                    Vector3::new(0.0, 1.0, 0.0),
+                    Vec3::new(0.0, 1.0, 0.0),
                     -2.0,
                     Arc::new(LambertianMaterial {
                         colour: ColourRgbF::new(0.55, 0.27, 0.04),
@@ -155,7 +154,7 @@ pub fn main() -> Result<(), Box<dyn std::error::Error>> {
                     }),
                 )) as Box<dyn Primitive>,
                 Box::new(Sphere::new(
-                    Point3::new(-6.25, -0.5, 1.0),
+                    Vec3::new(-6.25, -0.5, 1.0),
                     1.0,
                     Arc::new(LambertianMaterial {
                         colour: ColourRgbF::from_named(NamedColour::Green),
@@ -163,7 +162,7 @@ pub fn main() -> Result<(), Box<dyn std::error::Error>> {
                     }),
                 )),
                 Box::new(Sphere::new(
-                    Point3::new(-4.25, -0.5, 2.0),
+                    Vec3::new(-4.25, -0.5, 2.0),
                     1.0,
                     Arc::new(ReflectiveMaterial {
                         colour: ColourRgbF::from_named(NamedColour::Blue),
@@ -172,7 +171,7 @@ pub fn main() -> Result<(), Box<dyn std::error::Error>> {
                     }),
                 )),
                 Box::new(Sphere::new(
-                    Point3::new(-5.0, 1.5, 1.0),
+                    Vec3::new(-5.0, 1.5, 1.0),
                     1.0,
                     Arc::new(PhongMaterial {
                         colour: ColourRgbF::from_named(NamedColour::Red),

src/materials/lambertian_material.rs

@@ -1,6 +1,5 @@
-use nalgebra::Vector3;
-
 use crate::colour::ColourRgbF;
+use crate::math::Vec3;
 
 use super::{Bsdf, Material};
 
@@ -24,8 +23,6 @@ impl LambertianMaterial {
 impl Material for LambertianMaterial {
     fn bsdf(&self) -> Bsdf {
         let colour = self.colour * self.diffuse_strength;
-        Box::new(
+        Box::new(move |_w_o: Vec3, _w_i: Vec3, colour_in: ColourRgbF| colour * colour_in)
-            move |_w_o: Vector3<f64>, _w_i: Vector3<f64>, colour_in: ColourRgbF| colour * colour_in,
-        )
     }
 }

src/materials/mod.rs

@@ -1,10 +1,10 @@
-use nalgebra::Vector3;
+use crate::math::Vec3;
 
 use super::colour::ColourRgbF;
 
 use std::fmt::Debug;
 
-type Bsdf = Box<dyn Fn(Vector3<f64>, Vector3<f64>, ColourRgbF) -> ColourRgbF>;
+type Bsdf = Box<dyn Fn(Vec3, Vec3, ColourRgbF) -> ColourRgbF>;
 
 pub mod lambertian_material;
 pub use lambertian_material::LambertianMaterial;
@@ -21,7 +21,7 @@ pub use rgb_sampled_bsdf_material::RgbSampledBsdfMaterial;
 pub trait Material: Debug + Sync + Send {
     fn bsdf(&self) -> Bsdf;
 
-    fn sample(&self, _w_o: &Vector3<f64>) -> Vec<Vector3<f64>> {
+    fn sample(&self, _w_o: &Vec3) -> Vec<Vec3> {
         vec![]
     }
 }

src/materials/phong_material.rs

@@ -1,6 +1,5 @@
-use nalgebra::Vector3;
-
 use crate::colour::{ColourRgbF, NamedColour};
+use crate::math::Vec3;
 
 use std::fmt::Debug;
 
@@ -19,18 +18,16 @@ impl Material for PhongMaterial {
         let smoothness = self.smoothness;
         let specular_strength = self.specular_strength;
         let colour = self.colour * self.diffuse_strength;
-        Box::new(
+        Box::new(move |w_o: Vec3, w_i: Vec3, colour_in: ColourRgbF| {
-            move |w_o: Vector3<f64>, w_i: Vector3<f64>, colour_in: ColourRgbF| {
+            if w_i.z() < 0.0 || w_o.z() < 0.0 {
-                if w_i.z < 0.0 || w_o.z < 0.0 {
+                ColourRgbF::from_vec3(&Vec3::zeros())
-                    ColourRgbF::from_vector3(&Vector3::zeros())
+            } else {
-                } else {
+                let reflection_vector = Vec3::new(-w_i.x(), -w_i.y(), w_i.z());
-                    let reflection_vector = Vector3::new(-w_i.x, -w_i.y, w_i.z);
+                colour * colour_in
-                    colour * colour_in
+                    + ColourRgbF::from_named(NamedColour::White)
-                        + ColourRgbF::from_named(NamedColour::White)
+                        * w_o.dot(&reflection_vector).abs().powf(smoothness)
-                            * w_o.dot(&reflection_vector).abs().powf(smoothness)
+                        * (specular_strength / w_i.dot(&Vec3::unit_z()))
-                            * (specular_strength / w_i.dot(&Vector3::z_axis()))
+            }
-                }
+        })
-            },
-        )
     }
 }

src/materials/reflective_material.rs

@@ -1,10 +1,10 @@
-use nalgebra::{clamp, Vector3};
-
 use crate::colour::ColourRgbF;
+use crate::math::Vec3;
 
 use std::fmt::Debug;
 
 use super::{Bsdf, Material};
+
 #[derive(Debug)]
 pub struct ReflectiveMaterial {
     pub colour: ColourRgbF,
@@ -16,29 +16,26 @@ impl Material for ReflectiveMaterial {
     fn bsdf(&self) -> Bsdf {
         let diffuse_colour_factor = self.colour * self.diffuse_strength;
         let reflection_strength = self.reflection_strength;
-        Box::new(
+        Box::new(move |w_o: Vec3, w_i: Vec3, colour_in: ColourRgbF| {
-            move |w_o: Vector3<f64>, w_i: Vector3<f64>, colour_in: ColourRgbF| {
+            if w_i.z() <= 0.0 || w_o.z() <= 0.0 {
-                if w_i.z <= 0.0 || w_o.z <= 0.0 {
+                ColourRgbF::new(0.0, 0.0, 0.0)
-                    ColourRgbF::new(0.0, 0.0, 0.0)
+            } else {
-                } else {
+                let reflection_vector = Vec3::new(-w_o.x(), -w_o.y(), w_o.z());
-                    let reflection_vector = Vector3::new(-w_o.x, -w_o.y, w_o.z);
+                let reflection_colour = colour_in * reflection_strength;
-                    let reflection_colour = colour_in * reflection_strength;
+                let diffuse_colour = diffuse_colour_factor * colour_in;
-                    let diffuse_colour = diffuse_colour_factor * colour_in;
+                let sigma = 0.05;
-                    let sigma = 0.05;
+                let two = 2.0;
-                    let two = 2.0;
+                // These are normalized vectors, but sometimes rounding errors cause the
-                    // These are normalized vectors, but sometimes rounding errors cause the
+                // dot product to be slightly above 1 or below 0. The call to clamp
-                    // dot product to be slightly above 1 or below 0. The call to clamp
+                // ensures the values stay within the domain of acos,
-                    // ensures the values stay within the domain of acos,
+                let theta = w_i.dot(&reflection_vector).clamp(0.0, 1.0).abs().acos();
-                    let theta = clamp(w_i.dot(&reflection_vector), 0.0, 1.0).abs().acos();
+                let reflection_factor = (-(theta.powf(two)) / (two * sigma * sigma)).exp();
-                    let reflection_factor = (-(theta.powf(two)) / (two * sigma * sigma)).exp();
+                reflection_colour * reflection_factor + diffuse_colour * (1.0 - reflection_factor)
-                    reflection_colour * reflection_factor
+            }
-                        + diffuse_colour * (1.0 - reflection_factor)
+        })
-                }
-            },
-        )
     }
 
-    fn sample(&self, w_o: &Vector3<f64>) -> Vec<Vector3<f64>> {
+    fn sample(&self, w_o: &Vec3) -> Vec<Vec3> {
-        vec![Vector3::new(-w_o.x, -w_o.y, w_o.z)]
+        vec![Vec3::new(-w_o.x(), -w_o.y(), w_o.z())]
     }
 }

src/materials/rgb_sampled_bsdf_material.rs

@@ -1,7 +1,6 @@
-use nalgebra::Vector3;
-
 use super::{Bsdf, Material};
 use crate::colour::ColourRgbF;
+use crate::math::Vec3;
 use crate::realtype::NormalizedToU32;
 
 use std::error::Error;
@@ -13,7 +12,7 @@ use std::f64::consts::{FRAC_PI_2, PI};
 
 #[derive(Debug)]
 pub struct RgbSampledBsdfMaterial {
-    lut: Arc<Vec<Vec<Vec<Vec<Vector3<f64>>>>>>,
+    lut: Arc<Vec<Vec<Vec<Vec<Vec3>>>>>,
 }
 
 fn expand_and_index<T: Clone>(v: &mut Vec<T>, i: usize, default: T) -> &mut T {
@@ -48,8 +47,8 @@ impl RgbSampledBsdfMaterial {
                     Vec::new(),
                 ),
                 phi_out_index,
-                Vector3::zeros(),
+                Vec3::zeros(),
-            ) = Vector3::new(red, green, blue);
+            ) = Vec3::new(red, green, blue);
         }
         let lut = Arc::new(lut);
         Ok(RgbSampledBsdfMaterial { lut })
@@ -60,26 +59,26 @@ impl<'a> Material for RgbSampledBsdfMaterial {
     fn bsdf(&self) -> Bsdf {
         let lut = Arc::clone(&self.lut);
         Box::new(move |w_in, w_out, colour_in| {
-            if w_in.z < 0.0 || w_out.z < 0.0 {
+            if w_in.z() < 0.0 || w_out.z() < 0.0 {
                 return ColourRgbF::new(0.0, 0.0, 0.0);
             }
-            let theta_in = w_in.z.acos();
+            let theta_in = w_in.z().acos();
             let theta_in_index = (theta_in / FRAC_PI_2).normalized_to_u32(4) as usize;
-            let phi_in = w_in.y.atan2(w_in.x) + PI;
+            let phi_in = w_in.y().atan2(w_in.x()) + PI;
             let phi_in_index = (phi_in / (2.0 * PI)).normalized_to_u32(6) as usize;
-            let theta_out = w_out.z.acos();
+            let theta_out = w_out.z().acos();
             let theta_out_index = (theta_out / FRAC_PI_2).normalized_to_u32(4) as usize;
-            let phi_out = w_out.y.atan2(w_out.x) + PI;
+            let phi_out = w_out.y().atan2(w_out.x()) + PI;
             let phi_out_index = (phi_out / (2.0 * PI)).normalized_to_u32(6) as usize;
-            ColourRgbF::from_vector3(
+            ColourRgbF::from_vec3(
-                &colour_in.as_vector3().component_mul(
+                &colour_in.as_vec3().component_mul(
                     &lut[theta_in_index][phi_in_index][theta_out_index][phi_out_index],
                 ),
             )
         })
     }
 
-    fn sample(&self, w_o: &Vector3<f64>) -> Vec<Vector3<f64>> {
+    fn sample(&self, w_o: &Vec3) -> Vec<Vec3> {
-        vec![Vector3::new(-w_o.x, -w_o.y, w_o.z)]
+        vec![Vec3::new(-w_o.x(), -w_o.y(), w_o.z())]
     }
 }

src/math/mat3.rs

@@ -2,7 +2,9 @@ use super::Vec3;
 
 use std::ops::{Mul, MulAssign};
 
-#[derive(PartialEq, Debug)]
+use nalgebra::Matrix3;
+
+#[derive(PartialEq, Debug, Copy, Clone)]
 pub struct Mat3 {
     elements: [[f64; 3]; 3],
 }
@@ -24,6 +26,12 @@ impl Mat3 {
         }
     }
 
+    pub fn identity() -> Mat3 {
+        Mat3 {
+            elements: [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
+        }
+    }
+
     pub fn from_rows(r0: &Vec3, r1: &Vec3, r2: &Vec3) -> Mat3 {
         let mut elements = [[0.0; 3]; 3];
         for (row, v) in elements.iter_mut().zip([r0, r1, r2].iter()) {
@@ -51,6 +59,40 @@ impl Mat3 {
         }
         Vec3 { coords }
     }
+
+    fn from_nalgebra(m: &Matrix3<f64>) -> Mat3 {
+        Mat3::new(
+            *m.get((0, 0)).unwrap(),
+            *m.get((0, 1)).unwrap(),
+            *m.get((0, 2)).unwrap(),
+            *m.get((1, 0)).unwrap(),
+            *m.get((1, 1)).unwrap(),
+            *m.get((1, 2)).unwrap(),
+            *m.get((2, 0)).unwrap(),
+            *m.get((2, 1)).unwrap(),
+            *m.get((2, 2)).unwrap(),
+        )
+    }
+
+    fn to_nalgebra(&self) -> Matrix3<f64> {
+        Matrix3::new(
+            self.elements[0][0],
+            self.elements[0][1],
+            self.elements[0][2],
+            self.elements[1][0],
+            self.elements[1][1],
+            self.elements[1][2],
+            self.elements[2][0],
+            self.elements[2][1],
+            self.elements[2][2],
+        )
+    }
+
+    pub fn try_inverse(&self) -> Option<Self> {
+        self.to_nalgebra()
+            .try_inverse()
+            .map(|elem| Self::from_nalgebra(&elem))
+    }
 }
 
 impl Mul<Mat3> for Mat3 {
@@ -91,6 +133,18 @@ impl Mul<Vec3> for Mat3 {
     }
 }
 
+impl Mul<&Vec3> for Mat3 {
+    type Output = Vec3;
+
+    fn mul(self, rhs: &Vec3) -> Vec3 {
+        let mut coords = [0.0; 3];
+        for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
+            *coord = Vec3 { coords: *row }.dot(&rhs);
+        }
+        Vec3 { coords }
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;

src/math/mat4.rs

@@ -2,6 +2,8 @@ use super::Vec4;
 
 use std::ops::{Mul, MulAssign};
 
+use nalgebra::Matrix4;
+
 #[derive(PartialEq, Debug)]
 pub struct Mat4 {
     elements: [[f64; 4]; 4],
@@ -63,6 +65,54 @@ impl Mat4 {
         }
         Vec4 { coords }
     }
+
+    fn from_nalgebra(m: &Matrix4<f64>) -> Mat4 {
+        Mat4::new(
+            *m.get((0, 0)).unwrap(),
+            *m.get((0, 1)).unwrap(),
+            *m.get((0, 2)).unwrap(),
+            *m.get((0, 3)).unwrap(),
+            *m.get((1, 0)).unwrap(),
+            *m.get((1, 1)).unwrap(),
+            *m.get((1, 2)).unwrap(),
+            *m.get((1, 3)).unwrap(),
+            *m.get((2, 0)).unwrap(),
+            *m.get((2, 1)).unwrap(),
+            *m.get((2, 2)).unwrap(),
+            *m.get((2, 3)).unwrap(),
+            *m.get((3, 0)).unwrap(),
+            *m.get((3, 1)).unwrap(),
+            *m.get((3, 2)).unwrap(),
+            *m.get((3, 3)).unwrap(),
+        )
+    }
+
+    fn to_nalgebra(&self) -> Matrix4<f64> {
+        Matrix4::new(
+            self.elements[0][0],
+            self.elements[0][1],
+            self.elements[0][2],
+            self.elements[0][3],
+            self.elements[1][0],
+            self.elements[1][1],
+            self.elements[1][2],
+            self.elements[1][3],
+            self.elements[2][0],
+            self.elements[2][1],
+            self.elements[2][2],
+            self.elements[2][3],
+            self.elements[3][0],
+            self.elements[3][1],
+            self.elements[3][2],
+            self.elements[3][3],
+        )
+    }
+
+    pub fn try_inverse(&self) -> Option<Self> {
+        self.to_nalgebra()
+            .try_inverse()
+            .map(|mat| Self::from_nalgebra(&mat))
+    }
 }
 
 impl Mul<Mat4> for Mat4 {

src/math/vec2.rs

@@ -2,7 +2,7 @@ use itertools::izip;
 
 use std::ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign};
 
-#[derive(PartialEq, Debug)]
+#[derive(PartialEq, Debug, Copy, Clone)]
 pub struct Vec2 {
     coords: [f64; 2],
 }
@@ -27,6 +27,10 @@ impl Vec2 {
             .map(|(a_elem, b_elem)| a_elem * b_elem)
             .sum()
     }
+
+    pub fn perp(&self, rhs: &Vec2) -> f64 {
+        self.x() * rhs.y() - self.y() * rhs.x()
+    }
 }
 
 impl Add for Vec2 {
@@ -92,6 +96,13 @@ impl MulAssign<f64> for Vec2 {
 #[cfg(test)]
 mod tests {
     use super::*;
+    use quickcheck::{Arbitrary, Gen};
+
+    impl Arbitrary for Vec2 {
+        fn arbitrary<G: Gen>(g: &mut G) -> Vec2 {
+            Vec2::new(f64::arbitrary(g), f64::arbitrary(g))
+        }
+    }
 
     #[test]
     fn x_returns_first_element() {

src/math/vec3.rs

@@ -1,8 +1,10 @@
+use super::Mat3;
+
 use itertools::izip;
 
-use std::ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign};
+use std::ops::{Add, AddAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign};
 
-#[derive(PartialEq, Debug)]
+#[derive(Copy, Clone, PartialEq, Debug)]
 pub struct Vec3 {
     pub coords: [f64; 3],
 }
@@ -12,6 +14,49 @@ impl Vec3 {
         Vec3 { coords: [x, y, z] }
     }
 
+    pub fn from_slice(v: &[f64]) -> Self {
+        let mut coords = [0.0; 3];
+        coords.clone_from_slice(v);
+        Vec3 { coords }
+    }
+
+    /*pub fn from_iterator<I>(values: I) -> Vec3
+    where
+        I: Iterator<Item = f64>,
+    {
+        Vec3 {
+            coords: [
+                values.next().unwrap(),
+                values.next().unwrap(),
+                values.next().unwrap(),
+            ],
+        }
+    }*/
+
+    pub fn zeros() -> Vec3 {
+        Vec3 {
+            coords: [0.0, 0.0, 0.0],
+        }
+    }
+
+    pub fn unit_x() -> Vec3 {
+        Vec3 {
+            coords: [1.0, 0.0, 0.0],
+        }
+    }
+
+    pub fn unit_y() -> Vec3 {
+        Vec3 {
+            coords: [0.0, 1.0, 0.0],
+        }
+    }
+
+    pub fn unit_z() -> Vec3 {
+        Vec3 {
+            coords: [0.0, 0.0, 1.0],
+        }
+    }
+
     pub fn x(&self) -> f64 {
         self.coords[0]
     }
@@ -24,6 +69,10 @@ impl Vec3 {
         self.coords[2]
     }
 
+    pub fn as_slice(&self) -> &[f64] {
+        &self.coords
+    }
+
     pub fn dot(&self, rhs: &Vec3) -> f64 {
         self.coords
             .iter()
@@ -38,6 +87,94 @@ impl Vec3 {
         let z = self.x() * rhs.y() - self.y() * rhs.x();
         Vec3 { coords: [x, y, z] }
     }
+
+    pub fn abs(&self) -> Self {
+        Vec3::new(self.x().abs(), self.y().abs(), self.z().abs())
+    }
+
+    pub fn norm_squared(&self) -> f64 {
+        self.dot(&self)
+    }
+
+    pub fn norm(&self) -> f64 {
+        self.norm_squared().sqrt()
+    }
+
+    pub fn normalize(&self) -> Self {
+        let mut coords = [0.0; 3];
+        let inverse_norm = 1.0 / self.norm();
+        for (r, a) in coords.iter_mut().zip(self.coords.iter()) {
+            *r = a * inverse_norm;
+        }
+        Vec3 { coords }
+    }
+
+    pub fn smallest_coord(&self) -> usize {
+        let x = self.x().abs();
+        let y = self.y().abs();
+        let z = self.z().abs();
+        if x < y {
+            if x < z {
+                0
+            } else {
+                2
+            }
+        } else {
+            if y < z {
+                1
+            } else {
+                2
+            }
+        }
+    }
+
+    pub fn component_mul(&self, rhs: &Self) -> Self {
+        let mut coords = [0.0; 3];
+        for (elem, lhs_elem, rhs_elem) in
+            izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter())
+        {
+            *elem = lhs_elem * rhs_elem;
+        }
+        Vec3 { coords }
+    }
+}
+
+impl Index<usize> for Vec3 {
+    type Output = f64;
+
+    fn index(&self, i: usize) -> &f64 {
+        &self.coords[i]
+    }
+}
+
+impl IndexMut<usize> for Vec3 {
+    fn index_mut(&mut self, i: usize) -> &mut f64 {
+        &mut self.coords[i]
+    }
+}
+
+impl Add<Vec3> for &Vec3 {
+    type Output = Vec3;
+
+    fn add(self, rhs: Vec3) -> Vec3 {
+        let mut coords = [0.0; 3];
+        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+            *r = a + b;
+        }
+        Vec3 { coords }
+    }
+}
+
+impl Add<&Vec3> for &Vec3 {
+    type Output = Vec3;
+
+    fn add(self, rhs: &Vec3) -> Vec3 {
+        let mut coords = [0.0; 3];
+        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+            *r = a + b;
+        }
+        Vec3 { coords }
+    }
 }
 
 impl Add for Vec3 {
@@ -60,10 +197,29 @@ impl AddAssign for Vec3 {
     }
 }
 
-impl Sub for Vec3 {
+impl Neg for Vec3 {
-    type Output = Self;
+    type Output = Vec3;
+    fn neg(self) -> Vec3 {
+        Vec3::new(-self.x(), -self.y(), -self.z())
+    }
+}
 
-    fn sub(self, rhs: Self) -> Self {
+impl Sub for &Vec3 {
+    type Output = Vec3;
+
+    fn sub(self, rhs: Self) -> Vec3 {
+        let mut coords = [0.0; 3];
+        for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
+            *r = a - b;
+        }
+        Vec3 { coords }
+    }
+}
+
+impl Sub for Vec3 {
+    type Output = Vec3;
+
+    fn sub(self, rhs: Self) -> Vec3 {
         let mut coords = [0.0; 3];
         for (r, a, b) in izip!(coords.iter_mut(), self.coords.iter(), rhs.coords.iter()) {
             *r = a - b;
@@ -80,8 +236,20 @@ impl SubAssign for Vec3 {
     }
 }
 
+impl Mul<f64> for &Vec3 {
+    type Output = Vec3;
+
+    fn mul(self, rhs: f64) -> Vec3 {
+        let mut coords = [0.0; 3];
+        for (r, a) in coords.iter_mut().zip(self.coords.iter()) {
+            *r = a * rhs;
+        }
+        Vec3 { coords }
+    }
+}
+
 impl Mul<f64> for Vec3 {
-    type Output = Self;
+    type Output = Vec3;
 
     fn mul(self, rhs: f64) -> Vec3 {
         let mut coords = [0.0; 3];
@@ -100,9 +268,50 @@ impl MulAssign<f64> for Vec3 {
     }
 }
 
+impl Mul<Mat3> for &Vec3 {
+    type Output = Vec3;
+
+    fn mul(self, rhs: Mat3) -> Vec3 {
+        let mut coords = [0.0; 3];
+        for i in 0..3 {
+            coords[i] = self.dot(&rhs.get_column(i));
+        }
+        Vec3 { coords }
+    }
+}
+
+impl Mul<Mat3> for Vec3 {
+    type Output = Self;
+
+    fn mul(self, rhs: Mat3) -> Self {
+        let mut coords = [0.0; 3];
+        for i in 0..3 {
+            coords[i] = self.dot(&rhs.get_column(i));
+        }
+        Vec3 { coords }
+    }
+}
+
+impl MulAssign<Mat3> for Vec3 {
+    fn mul_assign(&mut self, rhs: Mat3) {
+        let mut coords = [0.0; 3];
+        for i in 0..3 {
+            coords[i] = self.dot(&rhs.get_column(i));
+        }
+        self.coords = coords;
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
+    use quickcheck::{Arbitrary, Gen};
+
+    impl Arbitrary for Vec3 {
+        fn arbitrary<G: Gen>(g: &mut G) -> Vec3 {
+            Vec3::new(f64::arbitrary(g), f64::arbitrary(g), f64::arbitrary(g))
+        }
+    }
 
     #[test]
     fn x_returns_first_element() {
@@ -122,6 +331,12 @@ mod tests {
         assert!(target.z() == 3.0);
     }
 
+    /*#[test]
+    fn from_iterator_takes_first_three_elements() {
+        let target = Vec3::from_iterator([1.0, 2.0, 3.0].iter());
+        assert!(target = Vec3::new(1.0, 2.0, 3.0));
+    }*/
+
     #[test]
     fn dot_product_returns_correct_result() {
         let a = Vec3::new(1.0, 2.0, 3.0);
@@ -137,6 +352,57 @@ mod tests {
         assert!(a.cross(&b) == c);
     }
 
+    #[test]
+    fn norm_returns_expected_value() {
+        let target = Vec3::new(2.0, 3.0, 6.0);
+        assert!(target.norm() == 7.0);
+    }
+
+    #[test]
+    fn normalized_vector_times_norm_yields_original() {
+        let mut target = Vec3::new(2.0, 3.0, 6.0);
+        let norm = target.norm();
+        target = target.normalize();
+        target *= norm;
+        assert!(target == Vec3::new(2.0, 3.0, 6.0));
+    }
+
+    #[test]
+    fn smallest_coord_works_for_x_when_positive() {
+        let target = Vec3::new(1.0, 2.0, 3.0);
+        assert!(target.smallest_coord() == 0);
+    }
+
+    #[test]
+    fn smallest_coord_works_for_x_when_negative() {
+        let target = Vec3::new(-2.0, -3.0, 3.0);
+        assert!(target.smallest_coord() == 0);
+    }
+
+    #[test]
+    fn smallest_coord_works_for_y_when_positive() {
+        let target = Vec3::new(2.0, 1.0, 3.0);
+        assert!(target.smallest_coord() == 1);
+    }
+
+    #[test]
+    fn smallest_coord_works_for_y_when_negative() {
+        let target = Vec3::new(-3.0, -2.0, 3.0);
+        assert!(target.smallest_coord() == 1);
+    }
+
+    #[test]
+    fn smallest_coord_works_for_z_when_positive() {
+        let target = Vec3::new(3.0, 2.0, 1.0);
+        assert!(target.smallest_coord() == 2);
+    }
+
+    #[test]
+    fn smallest_coord_works_for_z_when_negative() {
+        let target = Vec3::new(3.0, -3.0, -2.0);
+        assert!(target.smallest_coord() == 2);
+    }
+
     #[test]
     fn add_returns_correct_result() {
         let a = Vec3::new(1.0, 2.0, 3.0);
@@ -187,4 +453,29 @@ mod tests {
         a *= b;
         assert!(a == c);
     }
+
+    #[test]
+    fn mul_with_mat3_returns_expected_result() {
+        let a = Mat3::from_rows(
+            &Vec3::new(1.0, 2.0, 3.0),
+            &Vec3::new(4.0, 5.0, 6.0),
+            &Vec3::new(7.0, 8.0, 9.0),
+        );
+        let b = Vec3::new(10.0, 11.0, 12.0);
+        let c = Vec3::new(138.0, 171.0, 204.0);
+        assert!(b * a == c);
+    }
+
+    #[test]
+    fn mul_assign_with_mat3_returns_expected_result() {
+        let a = Mat3::from_rows(
+            &Vec3::new(1.0, 2.0, 3.0),
+            &Vec3::new(4.0, 5.0, 6.0),
+            &Vec3::new(7.0, 8.0, 9.0),
+        );
+        let mut b = Vec3::new(10.0, 11.0, 12.0);
+        let c = Vec3::new(138.0, 171.0, 204.0);
+        b *= a;
+        assert!(b == c);
+    }
 }

src/mesh.rs

@@ -1,10 +1,9 @@
 /// Load a model from a Wavefront .obj file
 mod wavefront_obj {
     use crate::materials::Material;
-
+    use crate::math::Vec3;
     use crate::raycasting::{Primitive, Triangle};
 
-    use nalgebra::{convert, Point3, Vector3};
     use obj::{IndexTuple, Obj, SimplePolygon};
 
     use std::io::Result;
@@ -15,14 +14,29 @@ mod wavefront_obj {
         index_tuple: &IndexTuple,
         vertex_positions: &[[f32; 3]],
         normal_positions: &[[f32; 3]],
-    ) -> (Point3<f64>, Vector3<f64>) {
+    ) -> (Vec3, Vec3) {
         let &IndexTuple(vertex_index, _, maybe_normal_index) = index_tuple;
-        let vertex = convert(Point3::from_slice(&vertex_positions[vertex_index]));
+        (
-        let normal = match maybe_normal_index {
+            {
-            Some(normal_index) => convert(Vector3::from_row_slice(&normal_positions[normal_index])),
+                let vertex_coords = &vertex_positions[vertex_index];
-            None => Vector3::zeros(),
+                Vec3::new(
-        };
+                    vertex_coords[0] as f64,
-        (vertex, normal.normalize())
+                    vertex_coords[1] as f64,
+                    vertex_coords[2] as f64,
+                )
+            },
+            match maybe_normal_index {
+                Some(normal_index) => {
+                    let normal_coords = &normal_positions[normal_index];
+                    Vec3::new(
+                        normal_coords[0] as f64,
+                        normal_coords[1] as f64,
+                        normal_coords[2] as f64,
+                    )
+                }
+                None => Vec3::zeros(),
+            },
+        )
     }
 
     fn get_triangles(

src/raycasting/axis_aligned_bounding_box.rs

@@ -9,9 +9,11 @@ pub use crate::util::axis_aligned_bounding_box::BoundingBox;
 impl IntersectP for BoundingBox {
     fn intersect(&self, ray: &Ray) -> bool {
         let mut t_interval_in_bounds = Interval::infinite();
-        for (&ray_origin, &ray_direction, bounds) in
+        for (&ray_origin, &ray_direction, bounds) in izip!(
-            izip!(ray.origin.iter(), ray.direction.iter(), self.bounds.iter())
+            ray.origin.coords.iter(),
-        {
+            ray.direction.coords.iter(),
+            self.bounds.iter()
+        ) {
             t_interval_in_bounds = t_interval_in_bounds.intersection(Interval::new(
                 (bounds.get_min() - ray_origin) / ray_direction,
                 (bounds.get_max() - ray_origin) / ray_direction,
@@ -28,11 +30,11 @@ impl IntersectP for BoundingBox {
 mod tests {
     use super::*;
 
+    use crate::math::Vec3;
+
     use quickcheck::TestResult;
     use quickcheck_macros::quickcheck;
 
-    use nalgebra::{Point3, Vector3};
-
     fn wrap_value_in_interval(value: f64, interval: Interval) -> f64 {
         let distance_from_start = (value - interval.get_min()).abs();
         let range = interval.get_max() - interval.get_min();
@@ -46,21 +48,20 @@ mod tests {
         interval.contains_value(wrap_value_in_interval(v, interval))
     }
 
-    fn wrap_point_into_bounding_box(point: Point3<f64>, bounds: &BoundingBox) -> Point3<f64> {
+    fn wrap_point_into_bounding_box(point: Vec3, bounds: &BoundingBox) -> Vec3 {
-        Point3::from(Vector3::from_iterator(
+        let mut coords = [0.0; 3];
-            point
+        for i in 0..3 {
-                .iter()
+            coords[i] = wrap_value_in_interval(point[i], bounds.bounds[i]);
-                .zip(bounds.bounds.iter())
+        }
-                .map(|(&value, &interval)| wrap_value_in_interval(value, interval)),
+        Vec3 { coords }
-        ))
     }
 
     #[quickcheck]
     fn correctly_detects_intersections(
-        ray_origin: Point3<f64>,
+        ray_origin: Vec3,
-        corner1: Point3<f64>,
+        corner1: Vec3,
-        corner2: Point3<f64>,
+        corner2: Vec3,
-        random_point: Point3<f64>,
+        random_point: Vec3,
     ) -> bool {
         let bounds = BoundingBox::from_corners(corner1, corner2);
         let point_in_bounds = wrap_point_into_bounding_box(random_point, &bounds);
@@ -70,10 +71,10 @@ mod tests {
 
     #[quickcheck]
     fn intersect_always_true_when_ray_origin_is_inside_bounds(
-        ray_origin: Point3<f64>,
+        ray_origin: Vec3,
-        corner1: Point3<f64>,
+        corner1: Vec3,
-        corner2: Point3<f64>,
+        corner2: Vec3,
-        random_point: Point3<f64>,
+        random_point: Vec3,
     ) -> TestResult {
         let bounds = BoundingBox::from_corners(corner1, corner2);
         let ray_origin = wrap_point_into_bounding_box(ray_origin, &bounds);
@@ -83,10 +84,10 @@ mod tests {
 
     #[quickcheck]
     fn no_intersection_when_behind_ray(
-        ray_origin: Point3<f64>,
+        ray_origin: Vec3,
-        corner1: Point3<f64>,
+        corner1: Vec3,
-        corner2: Point3<f64>,
+        corner2: Vec3,
-        random_point: Point3<f64>,
+        random_point: Vec3,
     ) -> TestResult {
         let bounds = BoundingBox::from_corners(corner1, corner2);
         if bounds.contains_point(ray_origin) {
@@ -100,24 +101,24 @@ mod tests {
     #[test]
     fn intersection_detected_when_ray_parallel_to_axis() {
         let target =
-            BoundingBox::from_corners(Point3::new(1.0f64, 2.0, 3.0), Point3::new(4.0, 5.0, 6.0));
+            BoundingBox::from_corners(Vec3::new(1.0f64, 2.0, 3.0), Vec3::new(4.0, 5.0, 6.0));
-        let x_ray = Ray::new(Point3::new(0.0, 3.0, 4.0), Vector3::new(1.0, 0.0, 0.0));
+        let x_ray = Ray::new(Vec3::new(0.0, 3.0, 4.0), Vec3::new(1.0, 0.0, 0.0));
         assert!(target.intersect(&x_ray));
-        let y_ray = Ray::new(Point3::new(2.0, 0.0, 4.0), Vector3::new(0.0, 1.0, 0.0));
+        let y_ray = Ray::new(Vec3::new(2.0, 0.0, 4.0), Vec3::new(0.0, 1.0, 0.0));
         assert!(target.intersect(&y_ray));
-        let z_ray = Ray::new(Point3::new(2.0, 3.0, 0.0), Vector3::new(0.0, 0.0, 1.0));
+        let z_ray = Ray::new(Vec3::new(2.0, 3.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
         assert!(target.intersect(&z_ray));
     }
 
     #[test]
     fn intersection_missed_when_ray_parallel_to_axis() {
         let target =
-            BoundingBox::from_corners(Point3::new(1.0f64, 2.0, 3.0), Point3::new(4.0, 5.0, 6.0));
+            BoundingBox::from_corners(Vec3::new(1.0f64, 2.0, 3.0), Vec3::new(4.0, 5.0, 6.0));
-        let x_ray = Ray::new(Point3::new(0.0, 0.0, 0.0), Vector3::new(1.0, 0.0, 0.0));
+        let x_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0));
         assert!(!target.intersect(&x_ray));
-        let y_ray = Ray::new(Point3::new(0.0, 0.0, 0.0), Vector3::new(0.0, 1.0, 0.0));
+        let y_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0));
         assert!(!target.intersect(&y_ray));
-        let z_ray = Ray::new(Point3::new(0.0, 0.0, 0.0), Vector3::new(0.0, 0.0, 1.0));
+        let z_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
         assert!(!target.intersect(&z_ray));
     }
 }

src/raycasting/bounding_volume_hierarchy.rs

@@ -1,9 +1,9 @@
+use crate::math::Vec3;
+
 use super::{
     Aggregate, BoundingBox, HasBoundingBox, Intersect, IntersectP, IntersectionInfo, Primitive, Ray,
 };
 
-use nalgebra::Point3;
-
 use std::cmp::Ordering;
 use std::sync::Arc;
 
@@ -27,8 +27,8 @@ pub enum BoundingVolumeHierarchy {
     },
 }
 
-fn centre(bounds: &BoundingBox) -> Point3<f64> {
+fn centre(bounds: &BoundingBox) -> Vec3 {
-    Point3::new(
+    Vec3::new(
         (bounds.bounds[0].get_min() + bounds.bounds[0].get_max()) / 2.00,
         (bounds.bounds[1].get_min() + bounds.bounds[1].get_max()) / 2.0,
         (bounds.bounds[2].get_min() + bounds.bounds[2].get_max()) / 2.0,

src/raycasting/mod.rs

@@ -1,4 +1,4 @@
-use nalgebra::{Affine3, Point3, Vector3};
+use crate::math::Vec3;
 
 use super::materials::Material;
 
@@ -28,17 +28,17 @@ pub mod vec_aggregate;
 #[derive(Clone, Debug)]
 pub struct Ray {
     /// The start point of the ray
-    pub origin: Point3<f64>,
+    pub origin: Vec3,
 
     /// The direction the ray goes in.
     ///
     /// This vector should always be kept normalized
-    pub direction: Vector3<f64>,
+    pub direction: Vec3,
 }
 
 impl Ray {
     /// Create a new ray
-    pub fn new(origin: Point3<f64>, direction: Vector3<f64>) -> Ray {
+    pub fn new(origin: Vec3, direction: Vec3) -> Ray {
         Ray {
             origin,
             direction: direction.normalize(),
@@ -46,7 +46,7 @@ impl Ray {
     }
 
     /// Return the point on the ray that is `t` units from the start
-    pub fn point_at(&self, t: f64) -> Point3<f64> {
+    pub fn point_at(&self, t: f64) -> Vec3 {
         self.origin + self.direction * t
     }
 
@@ -70,26 +70,26 @@ pub struct IntersectionInfo {
     pub distance: f64,
 
     /// The intersection point
-    pub location: Point3<f64>,
+    pub location: Vec3,
 
     /// The surface normal at the intersection point
-    pub normal: Vector3<f64>,
+    pub normal: Vec3,
 
     /// The surface tangent at the intersection point
     ///
     /// Which surface tangent direction returned is dependent on the [Primitive](Primitive)
     /// but should generally be smooth over any given surface
-    pub tangent: Vector3<f64>,
+    pub tangent: Vec3,
 
     /// Another surface tangent, perpendicular to `tangent`
     ///
     /// The cross product or `normal` and `tangent`
-    pub cotangent: Vector3<f64>,
+    pub cotangent: Vec3,
 
     /// The direction from the intersection point back towards the ray
     ///
     /// Equal to `-ray.direction`
-    pub retro: Vector3<f64>,
+    pub retro: Vec3,
 
     /// The [Material](crate::materials::Material) which describes the optical
     /// properties of the intersected surface
@@ -124,10 +124,10 @@ pub trait HasBoundingBox: Send + Sync {
 /// Any geometric object which can have an affine transformation applied to it
 ///
 /// Used for moving, rotating or scaling primitives
-pub trait Transform {
+/*pub trait Transform {
     /// Create a new object by applying the transformation to this object.
     fn transform(&self, transformation: &Affine3<f64>) -> Self;
-}
+}*/
 
 /// A basic geometric primitive such as a sphere or a triangle
 pub trait Primitive: Intersect + HasBoundingBox {
@@ -146,8 +146,8 @@ mod tests {
     use quickcheck::{Arbitrary, Gen};
     impl Arbitrary for Ray {
         fn arbitrary<G: Gen>(g: &mut G) -> Ray {
-            let origin = <Point3<f64> as Arbitrary>::arbitrary(g);
+            let origin = <Vec3 as Arbitrary>::arbitrary(g);
-            let direction = <Vector3<f64> as Arbitrary>::arbitrary(g);
+            let direction = <Vec3 as Arbitrary>::arbitrary(g);
             return Ray::new(origin, direction);
         }
     }

src/raycasting/plane.rs

@@ -1,29 +1,24 @@
-use nalgebra::{Affine3, Point3, Vector3};
-
 use crate::materials::Material;
+use crate::math::Vec3;
 
-use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray, Transform};
+use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray};
 
 use std::sync::Arc;
 
 #[derive(Clone)]
 pub struct Plane {
-    normal: Vector3<f64>,
+    normal: Vec3,
-    tangent: Vector3<f64>,
+    tangent: Vec3,
-    cotangent: Vector3<f64>,
+    cotangent: Vec3,
     distance_from_origin: f64,
     material: Arc<dyn Material>,
 }
 
 impl Plane {
-    pub fn new(
+    pub fn new(normal: Vec3, distance_from_origin: f64, material: Arc<dyn Material>) -> Plane {
-        normal: Vector3<f64>,
-        distance_from_origin: f64,
-        material: Arc<dyn Material>,
-    ) -> Plane {
         let normal = normal.normalize();
-        let mut axis_closest_to_tangent = Vector3::zeros();
+        let mut axis_closest_to_tangent = Vec3::zeros();
-        axis_closest_to_tangent[normal.iamin()] = 1.0;
+        axis_closest_to_tangent[normal.smallest_coord()] = 1.0;
         let cotangent = normal.cross(&axis_closest_to_tangent).normalize();
         let tangent = normal.cross(&cotangent);
         Plane {
@@ -36,7 +31,7 @@ impl Plane {
     }
 }
 
-impl Transform for Plane {
+/*impl Transform for Plane {
     fn transform(&self, transformation: &Affine3<f64>) -> Self {
         Plane {
             normal: transformation.transform_vector(&self.normal).normalize(),
@@ -48,14 +43,14 @@ impl Transform for Plane {
             material: Arc::clone(&self.material),
         }
     }
-}
+}*/
 
 impl Intersect for Plane {
     fn intersect<'a>(&'a self, ray: &Ray) -> Option<IntersectionInfo> {
         let ray_direction_dot_plane_normal = ray.direction.dot(&self.normal);
         let point_on_plane = self.normal * self.distance_from_origin;
         let point_on_plane_minus_ray_origin_dot_normal =
-            (point_on_plane - ray.origin.coords).dot(&self.normal);
+            (point_on_plane - ray.origin).dot(&self.normal);
         if ray_direction_dot_plane_normal == 0.0 {
             //Ray is parallel to plane
             if point_on_plane_minus_ray_origin_dot_normal != 0.0 {
@@ -81,11 +76,24 @@ impl Intersect for Plane {
 
 impl HasBoundingBox for Plane {
     fn bounding_box(&self) -> BoundingBox {
-        let p0 = Point3::from(self.normal * self.distance_from_origin);
+        let p0 = self.normal * self.distance_from_origin;
-        let f = |v: Vector3<f64>| {
+        let f = |v: Vec3| {
-            Vector3::from_iterator(
+            Vec3::new(
-                v.iter()
+                if v.x() == 0.0 {
-                    .map(|&elem| if elem == 0.0 { 0.0 } else { std::f64::INFINITY }),
+                    0.0
+                } else {
+                    std::f64::INFINITY
+                },
+                if v.y() == 0.0 {
+                    0.0
+                } else {
+                    std::f64::INFINITY
+                },
+                if v.z() == 0.0 {
+                    0.0
+                } else {
+                    std::f64::INFINITY
+                },
             )
         };
         let tangent = f(self.tangent);
@@ -102,16 +110,16 @@ impl Primitive for Plane {}
 
 #[cfg(test)]
 mod tests {
-    use nalgebra::Point3;
 
     use super::*;
     use crate::materials::LambertianMaterial;
+    use crate::math::Vec3;
 
     #[test]
     fn ray_intersects_plane() {
-        let r = Ray::new(Point3::new(1.0, 2.0, 3.0), Vector3::new(-1.0, 0.0, 1.0));
+        let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(-1.0, 0.0, 1.0));
         let p = Plane::new(
-            Vector3::new(1.0, 0.0, 0.0),
+            Vec3::new(1.0, 0.0, 0.0),
             -5.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
@@ -122,9 +130,9 @@ mod tests {
 
     #[test]
     fn ray_does_not_intersect_plane() {
-        let r = Ray::new(Point3::new(1.0, 2.0, 3.0), Vector3::new(1.0, 0.0, 1.0));
+        let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(1.0, 0.0, 1.0));
         let p = Plane::new(
-            Vector3::new(1.0, 0.0, 0.0),
+            Vec3::new(1.0, 0.0, 0.0),
             -5.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
@@ -135,9 +143,9 @@ mod tests {
 
     #[test]
     fn intersection_point_is_on_plane() {
-        let r = Ray::new(Point3::new(1.0, 2.0, 3.0), Vector3::new(-1.0, 0.0, 1.0));
+        let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(-1.0, 0.0, 1.0));
         let p = Plane::new(
-            Vector3::new(1.0, 0.0, 0.0),
+            Vec3::new(1.0, 0.0, 0.0),
             -5.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
@@ -150,7 +158,7 @@ mod tests {
                 cotangent: _,
                 retro: _,
                 material: _,
-            }) => assert!((location.x - (-5.0f64)).abs() < 0.0000000001),
+            }) => assert!((location.x() - (-5.0f64)).abs() < 0.0000000001),
             None => panic!(),
         }
     }
@@ -158,133 +166,133 @@ mod tests {
     #[test]
     fn bounding_box_is_correct_for_yz_plane() {
         let target = Plane::new(
-            Vector3::new(1.0, 0.0, 0.0),
+            Vec3::new(1.0, 0.0, 0.0),
             2.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
         let bb = target.bounding_box();
-        assert!(!bb.contains_point(Point3::new(1.0, 2.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(1.0, 2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(2.0, 2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(2.0, 2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(2.0, 2000.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(2.0, 2000.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(2.0, 0.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(2.0, 0.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(2.0, -2000.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(2.0, -2000.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(2.0, 2.0, 3000.0)));
+        assert!(bb.contains_point(Vec3::new(2.0, 2.0, 3000.0)));
-        assert!(bb.contains_point(Point3::new(2.0, 2.0, 0.0)));
+        assert!(bb.contains_point(Vec3::new(2.0, 2.0, 0.0)));
-        assert!(bb.contains_point(Point3::new(2.0, 2.0, -3000.0)));
+        assert!(bb.contains_point(Vec3::new(2.0, 2.0, -3000.0)));
-        assert!(!bb.contains_point(Point3::new(3.0, 2.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(3.0, 2.0, 3.0)));
     }
 
     #[test]
     fn bounding_box_is_correct_for_yz_plane_with_negative_normal() {
         let target = Plane::new(
-            Vector3::new(-1.0, 0.0, 0.0),
+            Vec3::new(-1.0, 0.0, 0.0),
             2.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
         let bb = target.bounding_box();
-        assert!(!bb.contains_point(Point3::new(1.0, 2.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(1.0, 2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(-2.0, 2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(-2.0, 2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(-2.0, 2000.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(-2.0, 2000.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(-2.0, 0.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(-2.0, 0.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(-2.0, -2000.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(-2.0, -2000.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(-2.0, 2.0, 3000.0)));
+        assert!(bb.contains_point(Vec3::new(-2.0, 2.0, 3000.0)));
-        assert!(bb.contains_point(Point3::new(-2.0, 2.0, 0.0)));
+        assert!(bb.contains_point(Vec3::new(-2.0, 2.0, 0.0)));
-        assert!(bb.contains_point(Point3::new(-2.0, 2.0, -3000.0)));
+        assert!(bb.contains_point(Vec3::new(-2.0, 2.0, -3000.0)));
-        assert!(!bb.contains_point(Point3::new(-3.0, 2.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(-3.0, 2.0, 3.0)));
     }
 
     #[test]
     fn bounding_box_is_correct_for_xz_plane() {
         let target = Plane::new(
-            Vector3::new(0.0, 1.0, 0.0),
+            Vec3::new(0.0, 1.0, 0.0),
             2.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
         let bb = target.bounding_box();
-        assert!(!bb.contains_point(Point3::new(1.0, 1.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(1.0, 1.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(1000.0, 2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(1000.0, 2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(0.0, 2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(0.0, 2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(-1000.0, 2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(-1000.0, 2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2.0, 3000.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2.0, 3000.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2.0, 0.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2.0, 0.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2.0, -3000.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2.0, -3000.0)));
-        assert!(!bb.contains_point(Point3::new(1.0, 3.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(1.0, 3.0, 3.0)));
     }
 
     #[test]
     fn bounding_box_is_correct_for_xz_plane_with_negative_normal() {
         let target = Plane::new(
-            Vector3::new(0.0, -1.0, 0.0),
+            Vec3::new(0.0, -1.0, 0.0),
             2.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
         let bb = target.bounding_box();
-        assert!(!bb.contains_point(Point3::new(1.0, -1.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(1.0, -1.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(1.0, -2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, -2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(1000.0, -2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(1000.0, -2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(0.0, -2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(0.0, -2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(-1000.0, -2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(-1000.0, -2.0, 3.0)));
-        assert!(bb.contains_point(Point3::new(1.0, -2.0, 3000.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, -2.0, 3000.0)));
-        assert!(bb.contains_point(Point3::new(1.0, -2.0, 0.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, -2.0, 0.0)));
-        assert!(bb.contains_point(Point3::new(1.0, -2.0, -3000.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, -2.0, -3000.0)));
-        assert!(!bb.contains_point(Point3::new(1.0, 3.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(1.0, 3.0, 3.0)));
     }
 
     #[test]
     fn bounding_box_is_correct_for_xy_plane() {
         let target = Plane::new(
-            Vector3::new(0.0, 0.0, 1.0),
+            Vec3::new(0.0, 0.0, 1.0),
             2.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
         let bb = target.bounding_box();
-        assert!(!bb.contains_point(Point3::new(1.0, 2.0, 1.0)));
+        assert!(!bb.contains_point(Vec3::new(1.0, 2.0, 1.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2000.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2000.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 0.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 0.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(1.0, -2000.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, -2000.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(2000.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(2000.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(0.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(0.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(-2000.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(-2000.0, 2.0, 2.0)));
-        assert!(!bb.contains_point(Point3::new(3.0, 2.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(3.0, 2.0, 3.0)));
     }
 
     #[test]
     fn bounding_box_is_correct_for_xy_plane_with_negative_normal() {
         let target = Plane::new(
-            Vector3::new(0.0, 0.0, -1.0),
+            Vec3::new(0.0, 0.0, -1.0),
             -2.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
         let bb = target.bounding_box();
-        assert!(!bb.contains_point(Point3::new(1.0, 2.0, 1.0)));
+        assert!(!bb.contains_point(Vec3::new(1.0, 2.0, 1.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2000.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2000.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 0.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 0.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(1.0, -2000.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, -2000.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(2000.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(2000.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(0.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(0.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(-2000.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(-2000.0, 2.0, 2.0)));
-        assert!(!bb.contains_point(Point3::new(3.0, 2.0, 3.0)));
+        assert!(!bb.contains_point(Vec3::new(3.0, 2.0, 3.0)));
     }
 
     #[test]
     fn bounding_box_is_infinite_when_normal_is_not_aligned_with_axis() {
         let target = Plane::new(
-            Vector3::new(0.1, 0.0, -1.0),
+            Vec3::new(0.1, 0.0, -1.0),
             -2.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
         let bb = target.bounding_box();
-        assert!(bb.contains_point(Point3::new(1.0, 2.0, 1.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2.0, 1.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 2000.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 2000.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(1.0, 0.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, 0.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(1.0, -2000.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(1.0, -2000.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(2000.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(2000.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(0.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(0.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(-2000.0, 2.0, 2.0)));
+        assert!(bb.contains_point(Vec3::new(-2000.0, 2.0, 2.0)));
-        assert!(bb.contains_point(Point3::new(3.0, 2.0, 3.0)));
+        assert!(bb.contains_point(Vec3::new(3.0, 2.0, 3.0)));
     }
 }

src/raycasting/sphere.rs

@@ -1,20 +1,19 @@
-use nalgebra::{Affine3, Point3, Vector3};
-
 use crate::materials::Material;
+use crate::math::Vec3;
 
-use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray, Transform};
+use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray};
 
 use std::sync::Arc;
 
 #[derive(Clone, Debug)]
 pub struct Sphere {
-    centre: Point3<f64>,
+    centre: Vec3,
     radius: f64,
     material: Arc<dyn Material>,
 }
 
 impl Sphere {
-    pub fn new(centre: Point3<f64>, radius: f64, material: Arc<dyn Material>) -> Sphere {
+    pub fn new(centre: Vec3, radius: f64, material: Arc<dyn Material>) -> Sphere {
         Sphere {
             centre,
             radius,
@@ -23,7 +22,7 @@ impl Sphere {
     }
 }
 
-impl Transform for Sphere {
+/*impl Transform for Sphere {
     fn transform(&self, transformation: &Affine3<f64>) -> Self {
         Sphere {
             centre: transformation.transform_point(&self.centre),
@@ -35,23 +34,26 @@ impl Transform for Sphere {
             material: Arc::clone(&self.material),
         }
     }
-}
+}*/
 
 impl Intersect for Sphere {
     fn intersect<'a>(&'a self, ray: &Ray) -> Option<IntersectionInfo> {
-        let r_o = ray.origin.coords;
+        let r_o = ray.origin;
-        let centre_coords = self.centre.coords;
+        let centre_coords = self.centre;
         let a = ray
             .direction
             .component_mul(&ray.direction)
+            .coords
             .iter()
             .fold(0.0, |a, b| a + *b);
         let b = ((r_o.component_mul(&ray.direction) - centre_coords.component_mul(&ray.direction))
             * 2.0)
+            .coords
             .iter()
             .fold(0.0, |a, b| a + *b);
         let c = (r_o.component_mul(&r_o) + centre_coords.component_mul(&centre_coords)
             - centre_coords.component_mul(&r_o) * 2.0)
+            .coords
             .iter()
             .fold(0.0, |a, b| a + *b)
             - self.radius * self.radius;
@@ -73,7 +75,7 @@ impl Intersect for Sphere {
             } else {
                 let location = ray.point_at(distance);
                 let normal = (location - self.centre).normalize();
-                let tangent = normal.cross(&Vector3::z_axis()).normalize();
+                let tangent = normal.cross(&Vec3::unit_z()).normalize();
                 let cotangent = normal.cross(&tangent);
                 let retro = -ray.direction;
                 Some(IntersectionInfo {
@@ -92,7 +94,7 @@ impl Intersect for Sphere {
 
 impl HasBoundingBox for Sphere {
     fn bounding_box(&self) -> BoundingBox {
-        let radius_xyz = Vector3::new(self.radius, self.radius, self.radius);
+        let radius_xyz = Vec3::new(self.radius, self.radius, self.radius);
         BoundingBox::from_corners(self.centre + radius_xyz, self.centre - radius_xyz)
     }
 }
@@ -104,16 +106,14 @@ mod tests {
     use quickcheck::TestResult;
     use quickcheck_macros::quickcheck;
 
-    use nalgebra::{Rotation3, Translation3};
-
     use super::*;
     use crate::materials::LambertianMaterial;
 
     #[test]
     fn ray_intersects_sphere() {
-        let r = Ray::new(Point3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0));
+        let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 0.0, 1.0));
         let s = Sphere::new(
-            Point3::new(1.5, 1.5, 15.0),
+            Vec3::new(1.5, 1.5, 15.0),
             5.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
@@ -124,9 +124,9 @@ mod tests {
 
     #[test]
     fn ray_does_not_intersect_sphere_when_sphere_is_in_front() {
-        let r = Ray::new(Point3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0));
+        let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 0.0, 1.0));
         let s = Sphere::new(
-            Point3::new(-5.0, 1.5, 15.0),
+            Vec3::new(-5.0, 1.5, 15.0),
             5.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
@@ -137,9 +137,9 @@ mod tests {
 
     #[test]
     fn ray_does_not_intersect_sphere_when_sphere_is_behind() {
-        let r = Ray::new(Point3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0));
+        let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 0.0, 1.0));
         let s = Sphere::new(
-            Point3::new(1.5, 1.5, -15.0),
+            Vec3::new(1.5, 1.5, -15.0),
             5.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
@@ -150,9 +150,9 @@ mod tests {
 
     #[test]
     fn ray_intersects_sphere_when_origin_is_inside() {
-        let r = Ray::new(Point3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0));
+        let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 0.0, 1.0));
         let s = Sphere::new(
-            Point3::new(1.5, 1.5, 2.0),
+            Vec3::new(1.5, 1.5, 2.0),
             5.0,
             Arc::new(LambertianMaterial::new_dummy()),
         );
@@ -163,8 +163,8 @@ mod tests {
 
     #[quickcheck]
     fn ray_intersects_sphere_centre_at_correct_distance(
-        ray_origin: Point3<f64>,
+        ray_origin: Vec3,
-        sphere_centre: Point3<f64>,
+        sphere_centre: Vec3,
         radius: f64,
     ) -> TestResult {
         if radius <= 0.0 || radius + 0.000001 >= (ray_origin - sphere_centre).norm() {
@@ -184,10 +184,7 @@ mod tests {
     }
 
     #[quickcheck]
-    fn all_points_on_sphere_are_in_bounding_box(
+    fn all_points_on_sphere_are_in_bounding_box(sphere_centre: Vec3, radius_vector: Vec3) -> bool {
-        sphere_centre: Point3<f64>,
-        radius_vector: Vector3<f64>,
-    ) -> bool {
         let target_sphere = Sphere::new(
             sphere_centre,
             radius_vector.norm(),
@@ -197,11 +194,11 @@ mod tests {
         bounding_box.contains_point(sphere_centre + radius_vector)
     }
 
-    #[quickcheck]
+    /*#[quickcheck]
     fn translation_moves_centre(
-        sphere_centre: Point3<f64>,
+        sphere_centre: Vec3,
         radius: f64,
-        translation_vector: Vector3<f64>,
+        translation_vector: Vec3,
     ) -> TestResult {
         if radius <= 0.0 {
             return TestResult::discard();
@@ -220,9 +217,9 @@ mod tests {
 
     #[quickcheck]
     fn translation_does_not_change_radius(
-        sphere_centre: Point3<f64>,
+        sphere_centre: Vec3,
         radius: f64,
-        translation_vector: Vector3<f64>,
+        translation_vector: Vec3,
     ) -> TestResult {
         if radius <= 0.0 {
             return TestResult::discard();
@@ -241,9 +238,9 @@ mod tests {
 
     #[quickcheck]
     fn rotation_about_centre_does_not_move_centre(
-        sphere_centre: Point3<f64>,
+        sphere_centre: Vec3,
         radius: f64,
-        rotation_vector: Vector3<f64>,
+        rotation_vector: Vec3,
     ) -> TestResult {
         if radius <= 0.0 {
             return TestResult::discard();
@@ -259,6 +256,6 @@ mod tests {
             * Rotation3::new(rotation_vector)
             * Translation3::from(-sphere.centre.coords);
         let sphere = sphere.transform(&transformation);
-        TestResult::from_bool(dbg!((expected_centre - sphere.centre).norm() < 0.000000001))
+        TestResult::from_bool((expected_centre - sphere.centre).norm() < 0.000000001)
-    }
+    }*/
 }

src/raycasting/triangle.rs

@@ -1,18 +1,18 @@
 use crate::materials::Material;
+use crate::math::{Vec2, Vec3};
 
-use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray, Transform};
+use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray};
-use nalgebra::{Affine3, Point3, Vector2, Vector3};
 
 use std::sync::Arc;
 
 #[derive(Debug, Clone)]
 pub struct Triangle {
-    pub vertices: [Point3<f64>; 3],
+    pub vertices: [Vec3; 3],
-    pub normals: [Vector3<f64>; 3],
+    pub normals: [Vec3; 3],
     pub material: Arc<dyn Material>,
 }
 
-impl Transform for Triangle {
+/*impl Transform for Triangle {
     fn transform(&self, transformation: &Affine3<f64>) -> Self {
         let normal_transformation =
             Affine3::from_matrix_unchecked(transformation.inverse().matrix().transpose());
@@ -30,51 +30,51 @@ impl Transform for Triangle {
             material: Arc::clone(&self.material),
         }
     }
-}
+}*/
 
 impl Intersect for Triangle {
     fn intersect<'a>(&'a self, ray: &Ray) -> Option<IntersectionInfo> {
-        let translation = -ray.origin.coords;
+        let translation = -ray.origin;
         let indices = indices_with_index_of_largest_element_last(&ray.direction);
         let permuted_ray_direction = permute_vector_elements(&ray.direction, &indices);
         let shear_slopes = calculate_shear_to_z_axis(&permuted_ray_direction);
-        let transformed_vertices: Vec<Vector3<f64>> = self
+        let transformed_vertices: Vec<Vec3> = self
             .vertices
             .iter()
             .map(|elem| {
                 apply_shear_to_z_axis(
-                    &permute_vector_elements(&(elem.coords + translation), &indices),
+                    &permute_vector_elements(&(elem + translation), &indices),
                     &shear_slopes,
                 )
             })
             .collect();
         let edge_functions = signed_edge_functions(&transformed_vertices);
-        if edge_functions.iter().all(|e| e.is_sign_positive())
+        if edge_functions.coords.iter().all(|e| e.is_sign_positive())
-            || edge_functions.iter().all(|e| e.is_sign_negative())
+            || edge_functions.coords.iter().all(|e| e.is_sign_negative())
         {
-            let barycentric_coordinates = barycentric_coordinates_from_signed_edge_functions(
+            let barycentric_coordinates =
-                Vector3::from_iterator(edge_functions.iter().map(|e| e.abs())),
+                barycentric_coordinates_from_signed_edge_functions(edge_functions.abs());
-            );
             let transformed_z = barycentric_coordinates
+                .coords
                 .iter()
                 .zip(transformed_vertices.iter())
-                .map(|(&coord, vertex)| vertex.z * coord)
+                .map(|(&coord, vertex)| vertex.z() * coord)
                 .fold(0.0, |acc, z| acc + z);
-            if transformed_z.is_sign_positive() != permuted_ray_direction.z.is_sign_positive() {
+            if transformed_z.is_sign_positive() != permuted_ray_direction.z().is_sign_positive() {
                 return None;
             }
             let location = barycentric_coordinates
+                .coords
                 .iter()
                 .zip(self.vertices.iter())
-                .map(|(&barycentric_coord, vertex)| vertex.coords * barycentric_coord)
+                .map(|(&barycentric_coord, vertex)| vertex * barycentric_coord)
-                .fold(Point3::new(0.0, 0.0, 0.0), |a, e| a + e);
+                .fold(Vec3::zeros(), |a, e| a + e);
             let distance = (ray.origin - location).norm();
-            let normal: Vector3<f64> = barycentric_coordinates
+            let normal: Vec3 = barycentric_coordinates
+                .coords
                 .iter()
                 .zip(self.normals.iter())
-                .fold(Vector3::zeros(), |acc, (&coord, vertex)| {
+                .fold(Vec3::zeros(), |acc, (&coord, vertex)| acc + vertex * coord)
-                    acc + vertex * coord
-                })
                 .normalize();
             let cotangent = (self.vertices[0] - self.vertices[1])
                 .cross(&normal)
@@ -105,15 +105,15 @@ impl HasBoundingBox for Triangle {
 
 impl Primitive for Triangle {}
 
-fn indices_with_index_of_largest_element_last(v: &Vector3<f64>) -> [usize; 3] {
+fn indices_with_index_of_largest_element_last(v: &Vec3) -> [usize; 3] {
-    if v.x > v.y {
+    if v.x() > v.y() {
-        if v.z > v.x {
+        if v.z() > v.x() {
             [0, 1, 2]
         } else {
             [1, 2, 0]
         }
     } else {
-        if v.z > v.y {
+        if v.z() > v.y() {
             [0, 1, 2]
         } else {
             [2, 0, 1]
@@ -125,61 +125,60 @@ fn is_valid_permutation(indices: &[usize; 3]) -> bool {
     (0..2).all(|i: usize| indices.iter().any(|&j| j == i))
 }
 
-fn permute_vector_elements(v: &Vector3<f64>, indices: &[usize; 3]) -> Vector3<f64> {
+fn permute_vector_elements(v: &Vec3, indices: &[usize; 3]) -> Vec3 {
     debug_assert!(is_valid_permutation(&indices));
-    Vector3::new(v[indices[0]], v[indices[1]], v[indices[2]])
+    Vec3::new(v[indices[0]], v[indices[1]], v[indices[2]])
 }
 
-fn calculate_shear_to_z_axis(v: &Vector3<f64>) -> Vector2<f64> {
+fn calculate_shear_to_z_axis(v: &Vec3) -> Vec2 {
-    Vector2::new(-v.x / v.z, -v.y / v.z)
+    Vec2::new(-v.x() / v.z(), -v.y() / v.z())
 }
 
-fn apply_shear_to_z_axis(v: &Vector3<f64>, s: &Vector2<f64>) -> Vector3<f64> {
+fn apply_shear_to_z_axis(v: &Vec3, s: &Vec2) -> Vec3 {
-    Vector3::new(v.x + s.x * v.z, v.y + s.y * v.z, v.z)
+    Vec3::new(v.x() + s.x() * v.z(), v.y() + s.y() * v.z(), v.z())
 }
 
-fn signed_edge_function(a: &Vector3<f64>, b: &Vector3<f64>) -> f64 {
+fn signed_edge_function(a: &Vec3, b: &Vec3) -> f64 {
-    a.x * b.y - b.x * a.y
+    a.x() * b.y() - b.x() * a.y()
 }
 
-fn signed_edge_functions(vertices: &[Vector3<f64>]) -> Vector3<f64> {
+fn signed_edge_functions(vertices: &[Vec3]) -> Vec3 {
     // Iterate over the inputs in such a way that each output element is calculated
     // from the twoother elements of the input. ( (y,z) -> x, (z,x) -> y, (x,y) -> z )
-    Vector3::from_iterator(
+    let coords: Vec<_> = vertices
-        vertices
+        .iter()
-            .iter()
+        .cycle()
-            .cycle()
+        .skip(1)
-            .skip(1)
+        .zip(vertices.iter().cycle().skip(2))
-            .zip(vertices.iter().cycle().skip(2))
+        .take(vertices.len())
-            .take(vertices.len())
+        .map(|(v1, v2)| signed_edge_function(v1, v2))
-            .map(|(v1, v2)| signed_edge_function(v1, v2)),
+        .collect();
-    )
+
+    Vec3::new(coords[0], coords[1], coords[2])
 }
 
-fn barycentric_coordinates_from_signed_edge_functions(e: Vector3<f64>) -> Vector3<f64> {
+fn barycentric_coordinates_from_signed_edge_functions(e: Vec3) -> Vec3 {
-    e * (1.0 / e.iter().fold(0.0, |a, &b| a + b))
+    e * (1.0 / e.coords.iter().fold(0.0, |a, &b| a + b))
 }
 
 #[cfg(test)]
 mod tests {
     use super::*;
 
-    mod triangle_transform {
+    /*mod triangle_transform {
         use super::*;
         use quickcheck_macros::quickcheck;
 
-        use nalgebra::Translation3;
-
         use crate::materials::LambertianMaterial;
 
         #[quickcheck]
         fn transform_by_identity_does_not_change_values(
-            v0: Point3<f64>,
+            v0: Vec3,
-            v1: Point3<f64>,
+            v1: Vec3,
-            v2: Point3<f64>,
+            v2: Vec3,
-            n0: Vector3<f64>,
+            n0: Vec3,
-            n1: Vector3<f64>,
+            n1: Vec3,
-            n2: Vector3<f64>,
+            n2: Vec3,
         ) -> bool {
             let n0 = n0.normalize();
             let n1 = n1.normalize();
@@ -200,13 +199,13 @@ mod tests {
 
         #[quickcheck]
         fn translate_does_not_change_normals(
-            v0: Point3<f64>,
+            v0: Vec3,
-            v1: Point3<f64>,
+            v1: Vec3,
-            v2: Point3<f64>,
+            v2: Vec3,
-            n0: Vector3<f64>,
+            n0: Vec3,
-            n1: Vector3<f64>,
+            n1: Vec3,
-            n2: Vector3<f64>,
+            n2: Vec3,
-            translation: Vector3<f64>,
+            translation: Vec3,
         ) -> bool {
             let n0 = n0.normalize();
             let n1 = n1.normalize();
@@ -223,13 +222,13 @@ mod tests {
 
         #[quickcheck]
         fn translate_translates_vertices(
-            v0: Point3<f64>,
+            v0: Vec3,
-            v1: Point3<f64>,
+            v1: Vec3,
-            v2: Point3<f64>,
+            v2: Vec3,
-            n0: Vector3<f64>,
+            n0: Vec3,
-            n1: Vector3<f64>,
+            n1: Vec3,
-            n2: Vector3<f64>,
+            n2: Vec3,
-            translation: Vector3<f64>,
+            translation: Vec3,
         ) -> bool {
             let n0 = n0.normalize();
             let n1 = n1.normalize();
@@ -245,52 +244,52 @@ mod tests {
                 && target.vertices[1] == v1 + translation
                 && target.vertices[2] == v2 + translation
         }
-    }
+    }*/
 
     mod index_of_largest_element {
         use super::*;
         use quickcheck_macros::quickcheck;
 
         #[quickcheck]
-        fn result_is_valid_permutation(v: Vector3<f64>) -> bool {
+        fn result_is_valid_permutation(v: Vec3) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             is_valid_permutation(&indices)
         }
 
         #[quickcheck]
-        fn result_includes_x(v: Vector3<f64>) -> bool {
+        fn result_includes_x(v: Vec3) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             indices.iter().any(|&i| i == 0)
         }
 
         #[quickcheck]
-        fn result_includes_y(v: Vector3<f64>) -> bool {
+        fn result_includes_y(v: Vec3) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             indices.iter().any(|&i| i == 1)
         }
 
         #[quickcheck]
-        fn result_includes_z(v: Vector3<f64>) -> bool {
+        fn result_includes_z(v: Vec3) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
             indices.iter().any(|&i| i == 2)
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_x(v: Vector3<f64>) -> bool {
+        fn last_index_is_greater_than_or_equal_to_x(v: Vec3) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
-            v[indices[2]] >= v.x
+            v[indices[2]] >= v.x()
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_y(v: Vector3<f64>) -> bool {
+        fn last_index_is_greater_than_or_equal_to_y(v: Vec3) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
-            v[indices[2]] >= v.y
+            v[indices[2]] >= v.y()
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_z(v: Vector3<f64>) -> bool {
+        fn last_index_is_greater_than_or_equal_to_z(v: Vec3) -> bool {
             let indices = indices_with_index_of_largest_element_last(&v);
-            v[indices[2]] >= v.z
+            v[indices[2]] >= v.z()
         }
     }
 
@@ -299,21 +298,21 @@ mod tests {
         use quickcheck_macros::quickcheck;
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_x(v: Vector3<f64>) -> bool {
+        fn last_index_is_greater_than_or_equal_to_x(v: Vec3) -> bool {
             let p = permute_vector_elements(&v, &indices_with_index_of_largest_element_last(&v));
-            p.z >= v.x
+            p.z() >= v.x()
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_y(v: Vector3<f64>) -> bool {
+        fn last_index_is_greater_than_or_equal_to_y(v: Vec3) -> bool {
             let p = permute_vector_elements(&v, &indices_with_index_of_largest_element_last(&v));
-            p.z >= v.y
+            p.z() >= v.y()
         }
 
         #[quickcheck]
-        fn last_index_is_greater_than_or_equal_to_z(v: Vector3<f64>) -> bool {
+        fn last_index_is_greater_than_or_equal_to_z(v: Vec3) -> bool {
             let p = permute_vector_elements(&v, &indices_with_index_of_largest_element_last(&v));
-            p.z >= v.z
+            p.z() >= v.z()
         }
     }
 
@@ -322,21 +321,21 @@ mod tests {
         use quickcheck_macros::quickcheck;
 
         #[quickcheck]
-        fn shear_to_z_axis_makes_x_zero(v: Vector3<f64>) -> bool {
+        fn shear_to_z_axis_makes_x_zero(v: Vec3) -> bool {
             let s = calculate_shear_to_z_axis(&v);
-            apply_shear_to_z_axis(&v, &s).x.abs() < 0.00001
+            apply_shear_to_z_axis(&v, &s).x().abs() < 0.00001
         }
 
         #[quickcheck]
-        fn shear_to_z_axis_makes_y_zero(v: Vector3<f64>) -> bool {
+        fn shear_to_z_axis_makes_y_zero(v: Vec3) -> bool {
             let s = calculate_shear_to_z_axis(&v);
-            apply_shear_to_z_axis(&v, &s).y.abs() < 0.00001
+            apply_shear_to_z_axis(&v, &s).y().abs() < 0.00001
         }
 
         #[quickcheck]
-        fn shear_to_z_axis_leaves_z_unchanged(v: Vector3<f64>) -> bool {
+        fn shear_to_z_axis_leaves_z_unchanged(v: Vec3) -> bool {
             let s = calculate_shear_to_z_axis(&v);
-            apply_shear_to_z_axis(&v, &s).z == v.z
+            apply_shear_to_z_axis(&v, &s).z() == v.z()
         }
     }
 
@@ -346,13 +345,10 @@ mod tests {
         use quickcheck_macros::quickcheck;
 
         #[quickcheck]
-        fn sign_of_signed_edge_function_matches_winding(
+        fn sign_of_signed_edge_function_matches_winding(a: Vec3, b: Vec3) -> TestResult {
-            a: Vector3<f64>,
+            let a_2d = Vec2::new(a.x(), a.y());
-            b: Vector3<f64>,
+            let b_2d = Vec2::new(b.x(), b.y());
-        ) -> TestResult {
+            let c_2d = Vec2::new(0.0, 0.0);
-            let a_2d = Vector2::new(a.x, a.y);
-            let b_2d = Vector2::new(b.x, b.y);
-            let c_2d = Vector2::new(0.0, 0.0);
             let winding = (b_2d - a_2d).perp(&(c_2d - b_2d));
             if winding.abs() < 0.00001 {
                 TestResult::discard()
@@ -365,9 +361,9 @@ mod tests {
 
         #[quickcheck]
         fn signed_edge_functions_has_same_result_as_signed_edge_function(
-            a: Vector3<f64>,
+            a: Vec3,
-            b: Vector3<f64>,
+            b: Vec3,
-            c: Vector3<f64>,
+            c: Vec3,
         ) -> bool {
             let es = signed_edge_functions(&vec![a, b, c]);
             es[0] == signed_edge_function(&b, &c)
@@ -376,16 +372,18 @@ mod tests {
         }
 
         #[quickcheck]
-        fn barycentric_coordinates_sum_to_one(
+        fn barycentric_coordinates_sum_to_one(a: Vec3, b: Vec3, c: Vec3) -> bool {
-            a: Vector3<f64>,
-            b: Vector3<f64>,
-            c: Vector3<f64>,
-        ) -> bool {
             let barycentric_coordinates =
                 barycentric_coordinates_from_signed_edge_functions(signed_edge_functions(&vec![
                     a, b, c,
                 ]));
-            (barycentric_coordinates.iter().fold(0.0, |a, b| a + b) - 1.0).abs() < 0.00000001
+            (barycentric_coordinates
+                .coords
+                .iter()
+                .fold(0.0, |a, b| a + b)
+                - 1.0)
+                .abs()
+                < 0.00000001
         }
     }
 
@@ -399,14 +397,14 @@ mod tests {
         fn intersection_passes_with_ray_along_z_axis_ccw_winding() {
             let target_triangle = Triangle {
                 vertices: [
-                    Point3::new(0.0, 1.0, 1.0),
+                    Vec3::new(0.0, 1.0, 1.0),
-                    Point3::new(1.0, -1.0, 1.0),
+                    Vec3::new(1.0, -1.0, 1.0),
-                    Point3::new(-1.0, -1.0, 1.0),
+                    Vec3::new(-1.0, -1.0, 1.0),
                 ],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
-            let target_ray = Ray::new(Point3::new(0.0, 0.0, 0.0), Vector3::new(0.0, 0.0, 1.0));
+            let target_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
             if let None = target_triangle.intersect(&target_ray) {
                 panic!()
             }
@@ -416,14 +414,14 @@ mod tests {
         fn intersection_passes_with_ray_along_z_axis_cw_winding() {
             let target_triangle = Triangle {
                 vertices: [
-                    Point3::new(0.0, 1.0, 1.0),
+                    Vec3::new(0.0, 1.0, 1.0),
-                    Point3::new(-1.0, -1.0, 1.0),
+                    Vec3::new(-1.0, -1.0, 1.0),
-                    Point3::new(1.0, -1.0, 1.0),
+                    Vec3::new(1.0, -1.0, 1.0),
                 ],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
-            let target_ray = Ray::new(Point3::new(0.0, 0.0, 0.0), Vector3::new(0.0, 0.0, 1.0));
+            let target_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
             if let None = target_triangle.intersect(&target_ray) {
                 panic!()
             }
@@ -433,14 +431,14 @@ mod tests {
         fn intersection_passes_with_ray_along_nagative_z_axis_ccw_winding() {
             let target_triangle = Triangle {
                 vertices: [
-                    Point3::new(0.0, 1.0, -1.0),
+                    Vec3::new(0.0, 1.0, -1.0),
-                    Point3::new(1.0, -1.0, -1.0),
+                    Vec3::new(1.0, -1.0, -1.0),
-                    Point3::new(-1.0, -1.0, -1.0),
+                    Vec3::new(-1.0, -1.0, -1.0),
                 ],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
-            let target_ray = Ray::new(Point3::new(0.0, 0.0, 0.0), Vector3::new(0.0, 0.0, -1.0));
+            let target_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, -1.0));
             if let None = target_triangle.intersect(&target_ray) {
                 panic!()
             }
@@ -450,14 +448,14 @@ mod tests {
         fn intersection_passes_with_ray_along_negativez_axis_cw_winding() {
             let target_triangle = Triangle {
                 vertices: [
-                    Point3::new(0.0, 1.0, -1.0),
+                    Vec3::new(0.0, 1.0, -1.0),
-                    Point3::new(-1.0, -1.0, -1.0),
+                    Vec3::new(-1.0, -1.0, -1.0),
-                    Point3::new(1.0, -1.0, -1.0),
+                    Vec3::new(1.0, -1.0, -1.0),
                 ],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
-            let target_ray = Ray::new(Point3::new(0.0, 0.0, 0.0), Vector3::new(0.0, 0.0, -1.0));
+            let target_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, -1.0));
             if let None = target_triangle.intersect(&target_ray) {
                 panic!()
             }
@@ -467,14 +465,14 @@ mod tests {
         fn intersection_passes_with_ray_along_z_axis_but_translated_ccw_winding() {
             let target_triangle = Triangle {
                 vertices: [
-                    Point3::new(5.0, 6.0, 6.0),
+                    Vec3::new(5.0, 6.0, 6.0),
-                    Point3::new(6.0, 4.0, 6.0),
+                    Vec3::new(6.0, 4.0, 6.0),
-                    Point3::new(4.0, 4.0, 6.0),
+                    Vec3::new(4.0, 4.0, 6.0),
                 ],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
-            let target_ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Vector3::new(0.0, 0.0, 1.0));
+            let target_ray = Ray::new(Vec3::new(5.0, 5.0, 5.0), Vec3::new(0.0, 0.0, 1.0));
             if let None = target_triangle.intersect(&target_ray) {
                 panic!()
             }
@@ -484,32 +482,32 @@ mod tests {
         fn intersection_passes_with_ray_at_angle_to_z_axisand_translated_ccw_winding() {
             let target_triangle = Triangle {
                 vertices: [
-                    Point3::new(6.0, 6.5, 6.0),
+                    Vec3::new(6.0, 6.5, 6.0),
-                    Point3::new(7.0, 4.5, 6.0),
+                    Vec3::new(7.0, 4.5, 6.0),
-                    Point3::new(5.0, 4.5, 6.0),
+                    Vec3::new(5.0, 4.5, 6.0),
                 ],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
-            let target_ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Vector3::new(1.0, 0.5, 1.0));
+            let target_ray = Ray::new(Vec3::new(5.0, 5.0, 5.0), Vec3::new(1.0, 0.5, 1.0));
             if let None = target_triangle.intersect(&target_ray) {
                 panic!()
             }
         }
 
         fn intersect_with_centroid_and_test_result<
-            F: Fn(Option<IntersectionInfo>, Point3<f64>) -> bool,
+            F: Fn(Option<IntersectionInfo>, Vec3) -> bool,
         >(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
             test: F,
         ) -> TestResult {
-            let centroid: Point3<f64> = [vertex0.coords, vertex1.coords, vertex2.coords]
+            let centroid: Vec3 = [vertex0, vertex1, vertex2]
                 .iter()
-                .fold(Point3::new(0.0, 0.0, 0.0), |acc, &elem| acc + elem)
+                .fold(Vec3::new(0.0, 0.0, 0.0), |acc, &elem| acc + elem)
-                / 3.0;
+                * (1.0 / 3.0);
             let ray_direction = (centroid - ray_origin).normalize();
             let normal = (vertex1 - vertex0).cross(&(vertex2 - vertex0)).normalize();
             if normal.dot(&ray_direction).abs() < 0.000_000_1 {
@@ -518,9 +516,9 @@ mod tests {
             }
             let target_triangle = Triangle {
                 vertices: [
-                    Point3::from(vertex0),
+                    Vec3::from(vertex0),
-                    Point3::from(vertex1),
+                    Vec3::from(vertex1),
-                    Point3::from(vertex2),
+                    Vec3::from(vertex2),
                 ],
                 normals: [normal; 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
@@ -532,15 +530,15 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_hits(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
         ) -> TestResult {
-            let centroid: Point3<f64> = [vertex0.coords, vertex1.coords, vertex2.coords]
+            let centroid: Vec3 = [vertex0, vertex1, vertex2]
                 .iter()
-                .fold(Point3::new(0.0, 0.0, 0.0), |acc, &elem| acc + elem)
+                .fold(Vec3::new(0.0, 0.0, 0.0), |acc, &elem| acc + elem)
-                / 3.0;
+                * (1.0 / 3.0);
             let ray_direction = (centroid - ray_origin).normalize();
             let normal = (vertex1 - vertex0).cross(&(vertex2 - vertex0)).normalize();
             if normal.dot(&ray_direction).abs() < 0.000_000_1 {
@@ -549,9 +547,9 @@ mod tests {
             }
             let target_triangle = Triangle {
                 vertices: [
-                    Point3::from(vertex0),
+                    Vec3::from(vertex0),
-                    Point3::from(vertex1),
+                    Vec3::from(vertex1),
-                    Point3::from(vertex2),
+                    Vec3::from(vertex2),
                 ],
                 normals: [normal; 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
@@ -567,10 +565,10 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_hits_centroid(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
         ) -> TestResult {
             intersect_with_centroid_and_test_result(
                 vertex0,
@@ -589,10 +587,10 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_hits_at_expected_distance(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
         ) -> TestResult {
             intersect_with_centroid_and_test_result(
                 vertex0,
@@ -611,10 +609,10 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_has_expected_normal(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
         ) -> TestResult {
             intersect_with_centroid_and_test_result(
                 vertex0,
@@ -635,10 +633,10 @@ mod tests {
 
         #[quickcheck]
         fn intersection_with_centroid_has_expected_retro(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
         ) -> TestResult {
             intersect_with_centroid_and_test_result(
                 vertex0,
@@ -674,18 +672,18 @@ mod tests {
         }
 
         fn intersect_with_barycentric_and_test_result<
-            F: Fn(Option<IntersectionInfo>, Point3<f64>) -> bool,
+            F: Fn(Option<IntersectionInfo>, Vec3) -> bool,
         >(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
             barycentric_coords: BarycentricCoords,
             test: F,
         ) -> TestResult {
             let point = vertex0 * barycentric_coords.alpha
-                + vertex1.coords * barycentric_coords.beta
+                + vertex1 * barycentric_coords.beta
-                + vertex2.coords * barycentric_coords.gamma;
+                + vertex2 * barycentric_coords.gamma;
             let ray_direction = (point - ray_origin).normalize();
             let normal = (vertex1 - vertex0).cross(&(vertex2 - vertex0)).normalize();
             if normal.dot(&ray_direction).abs() < 0.000_000_1 {
@@ -694,9 +692,9 @@ mod tests {
             }
             let target_triangle = Triangle {
                 vertices: [
-                    Point3::from(vertex0),
+                    Vec3::from(vertex0),
-                    Point3::from(vertex1),
+                    Vec3::from(vertex1),
-                    Point3::from(vertex2),
+                    Vec3::from(vertex2),
                 ],
                 normals: [normal; 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
@@ -708,10 +706,10 @@ mod tests {
 
         #[quickcheck]
         fn point_with_arbitrary_barycentric_coords_hits(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
             barycentric_coords: BarycentricCoords,
         ) -> TestResult {
             intersect_with_barycentric_and_test_result(
@@ -732,10 +730,10 @@ mod tests {
 
         #[quickcheck]
         fn point_with_arbitrary_barycentric_coords_has_expected_normal(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
             barycentric_coords: BarycentricCoords,
         ) -> TestResult {
             intersect_with_barycentric_and_test_result(
@@ -758,10 +756,10 @@ mod tests {
 
         #[quickcheck]
         fn point_with_arbitrary_barycentric_coords_has_expected_distance(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
             barycentric_coords: BarycentricCoords,
         ) -> TestResult {
             intersect_with_barycentric_and_test_result(
@@ -783,10 +781,10 @@ mod tests {
 
         #[quickcheck]
         fn point_with_arbitrary_barycentric_coords_has_expected_retro(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
             barycentric_coords: BarycentricCoords,
         ) -> TestResult {
             intersect_with_barycentric_and_test_result(
@@ -808,24 +806,24 @@ mod tests {
 
         #[quickcheck]
         fn intersection_fails_when_ray_outside_first_edge(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
-            uv: Vector2<f64>,
+            uv: Vec2,
         ) -> bool {
-            let uv_origin = Point3::from(vertex0);
+            let uv_origin = Vec3::from(vertex0);
             let u_axis = (vertex1 - vertex0).normalize();
             let w_axis = (vertex2 - vertex0).cross(&u_axis).normalize();
             let v_axis = w_axis.cross(&u_axis);
-            let target_point = uv_origin + u_axis * uv.x + v_axis * uv.y.abs();
+            let target_point = uv_origin + u_axis * uv.x() + v_axis * uv.y().abs();
             let ray = Ray {
                 origin: ray_origin,
                 direction: (target_point - ray_origin).normalize(),
             };
             let triangle = Triangle {
                 vertices: [vertex0, vertex1, vertex2],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
             match triangle.intersect(&ray) {
@@ -836,24 +834,24 @@ mod tests {
 
         #[quickcheck]
         fn intersection_fails_when_ray_outside_second_edge(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
-            uv: Vector2<f64>,
+            uv: Vec2,
         ) -> bool {
-            let uv_origin = Point3::from(vertex0);
+            let uv_origin = Vec3::from(vertex0);
             let u_axis = (vertex2 - vertex1).normalize();
             let w_axis = (vertex1 - vertex0).cross(&u_axis).normalize();
             let v_axis = w_axis.cross(&u_axis);
-            let target_point = uv_origin + u_axis * uv.x + v_axis * uv.y.abs();
+            let target_point = uv_origin + u_axis * uv.x() + v_axis * uv.y().abs();
             let ray = Ray {
                 origin: ray_origin,
                 direction: (target_point - ray_origin).normalize(),
             };
             let triangle = Triangle {
                 vertices: [vertex0, vertex1, vertex2],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
             match triangle.intersect(&ray) {
@@ -864,24 +862,24 @@ mod tests {
 
         #[quickcheck]
         fn intersection_fails_when_ray_outside_third_edge(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
-            uv: Vector2<f64>,
+            uv: Vec2,
         ) -> bool {
-            let uv_origin = Point3::from(vertex0);
+            let uv_origin = Vec3::from(vertex0);
             let u_axis = (vertex0 - vertex2).normalize();
             let w_axis = (vertex1 - vertex2).cross(&u_axis).normalize();
             let v_axis = w_axis.cross(&u_axis);
-            let target_point = uv_origin + u_axis * uv.x + v_axis * uv.y.abs();
+            let target_point = uv_origin + u_axis * uv.x() + v_axis * uv.y().abs();
             let ray = Ray {
                 origin: ray_origin,
                 direction: (target_point - ray_origin).normalize(),
             };
             let triangle = Triangle {
                 vertices: [vertex0, vertex1, vertex2],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
             match triangle.intersect(&ray) {
@@ -892,22 +890,22 @@ mod tests {
 
         #[quickcheck]
         fn intersection_fails_when_triangle_is_behind_ray(
-            vertex0: Point3<f64>,
+            vertex0: Vec3,
-            vertex1: Point3<f64>,
+            vertex1: Vec3,
-            vertex2: Point3<f64>,
+            vertex2: Vec3,
-            ray_origin: Point3<f64>,
+            ray_origin: Vec3,
             barycentric_coords: BarycentricCoords,
         ) -> bool {
-            let point_behind_ray = vertex0.coords * barycentric_coords.alpha
+            let point_behind_ray = vertex0 * barycentric_coords.alpha
-                + vertex1.coords * barycentric_coords.beta
+                + vertex1 * barycentric_coords.beta
-                + vertex2.coords * barycentric_coords.gamma;
+                + vertex2 * barycentric_coords.gamma;
             let ray = Ray {
                 origin: ray_origin,
-                direction: (ray_origin.coords - point_behind_ray).normalize(),
+                direction: (ray_origin - point_behind_ray).normalize(),
             };
             let triangle = Triangle {
                 vertices: [vertex0, vertex1, vertex2],
-                normals: [Vector3::zeros(); 3],
+                normals: [Vec3::zeros(); 3],
                 material: Arc::new(LambertianMaterial::new_dummy()),
             };
             match triangle.intersect(&ray) {

src/scene.rs

@@ -1,8 +1,8 @@
-use nalgebra::Point3;
+use crate::math::Vec3;
 
 use crate::raycasting::Aggregate;
 
 pub struct Scene {
-    pub camera_location: Point3<f64>,
+    pub camera_location: Vec3,
     pub objects: Vec<Box<dyn Aggregate>>,
 }

src/util/algebra_utils.rs

@@ -1,15 +1,7 @@
-use nalgebra::{Matrix3, Vector3};
+use crate::math::{Mat3, Vec3};
 
-pub fn try_change_of_basis_matrix(
+pub fn try_change_of_basis_matrix(x: &Vec3, y: &Vec3, z: &Vec3) -> Option<Mat3> {
-    x: &Vector3<f64>,
+    Some(Mat3::from_rows(x, y, z))
-    y: &Vector3<f64>,
-    z: &Vector3<f64>,
-) -> Option<Matrix3<f64>> {
-    Some(Matrix3::from_rows(&[
-        x.transpose(),
-        y.transpose(),
-        z.transpose(),
-    ]))
 }
 
 #[cfg(test)]
@@ -23,49 +15,37 @@ mod tests {
 
         #[test]
         fn produces_isentity_when_passed_axes() {
-            let target: Matrix3<f64> = try_change_of_basis_matrix(
+            let target: Mat3 =
-                &Vector3::x_axis(),
+                try_change_of_basis_matrix(&Vec3::unit_x(), &Vec3::unit_y(), &Vec3::unit_z())
-                &Vector3::y_axis(),
+                    .unwrap();
-                &Vector3::z_axis(),
+            assert!(target == Mat3::identity())
-            )
-            .unwrap();
-            assert!(target == Matrix3::identity())
         }
 
         #[quickcheck]
-        fn swap_xy_does_not_change_z(v: Vector3<f64>) {
+        fn swap_xy_does_not_change_z(v: Vec3) {
-            let target: Matrix3<f64> = try_change_of_basis_matrix(
+            let target: Mat3 =
-                &Vector3::y_axis(),
+                try_change_of_basis_matrix(&Vec3::unit_y(), &Vec3::unit_x(), &Vec3::unit_z())
-                &Vector3::x_axis(),
+                    .unwrap();
-                &Vector3::z_axis(),
-            )
-            .unwrap();
             let v2 = target * v;
-            assert!(v2.z == v.z)
+            assert!(v2.z() == v.z())
         }
 
         #[quickcheck]
-        fn swap_xy_copies_y_to_x(v: Vector3<f64>) {
+        fn swap_xy_copies_y_to_x(v: Vec3) {
-            let target: Matrix3<f64> = try_change_of_basis_matrix(
+            let target: Mat3 =
-                &Vector3::y_axis(),
+                try_change_of_basis_matrix(&Vec3::unit_y(), &Vec3::unit_x(), &Vec3::unit_z())
-                &Vector3::x_axis(),
+                    .unwrap();
-                &Vector3::z_axis(),
-            )
-            .unwrap();
             let v2 = target * v;
-            assert!(v2.x == v.y)
+            assert!(v2.x() == v.y())
         }
 
         #[quickcheck]
-        fn swap_xy_copies_x_to_y(v: Vector3<f64>) {
+        fn swap_xy_copies_x_to_y(v: Vec3) {
-            let target: Matrix3<f64> = try_change_of_basis_matrix(
+            let target: Mat3 =
-                &Vector3::y_axis(),
+                try_change_of_basis_matrix(&Vec3::unit_y(), &Vec3::unit_x(), &Vec3::unit_z())
-                &Vector3::x_axis(),
+                    .unwrap();
-                &Vector3::z_axis(),
-            )
-            .unwrap();
             let v2 = target * v;
-            assert!(v2.y == v.x)
+            assert!(v2.y() == v.x())
         }
     }
 }

src/util/axis_aligned_bounding_box.rs

@@ -1,5 +1,4 @@
-use nalgebra::Point3;
+use crate::math::Vec3;
-
 use crate::util::Interval;
 
 use itertools::izip;
@@ -10,11 +9,13 @@ pub struct BoundingBox {
 }
 
 impl BoundingBox {
-    pub fn from_corners(a: Point3<f64>, b: Point3<f64>) -> Self {
+    pub fn from_corners(a: Vec3, b: Vec3) -> Self {
         let mut result = BoundingBox {
             bounds: [Interval::infinite(); 3],
         };
-        for (bounds_elem, a_elem, b_elem) in izip!(result.bounds.iter_mut(), a.iter(), b.iter()) {
+        for (bounds_elem, a_elem, b_elem) in
+            izip!(result.bounds.iter_mut(), a.coords.iter(), b.coords.iter())
+        {
             *bounds_elem = Interval::new(*a_elem, *b_elem);
         }
         result
@@ -26,39 +27,39 @@ impl BoundingBox {
         }
     }
 
-    pub fn from_point(p: Point3<f64>) -> Self {
+    pub fn from_point(p: Vec3) -> Self {
         BoundingBox {
             bounds: [
-                Interval::degenerate(p.x),
+                Interval::degenerate(p.x()),
-                Interval::degenerate(p.y),
+                Interval::degenerate(p.y()),
-                Interval::degenerate(p.z),
+                Interval::degenerate(p.z()),
             ],
         }
     }
 
     pub fn from_points<'a, I>(points: I) -> Self
     where
-        I: IntoIterator<Item = &'a Point3<f64>>,
+        I: IntoIterator<Item = &'a Vec3>,
     {
         points
             .into_iter()
             .fold(BoundingBox::empty(), |acc, p| acc.expand_to_point(p))
     }
 
-    pub fn expand_to_point(&self, p: &Point3<f64>) -> Self {
+    pub fn expand_to_point(&self, p: &Vec3) -> Self {
         BoundingBox {
             bounds: [
-                self.bounds[0].expand_to_value(p.x),
+                self.bounds[0].expand_to_value(p.x()),
-                self.bounds[1].expand_to_value(p.y),
+                self.bounds[1].expand_to_value(p.y()),
-                self.bounds[2].expand_to_value(p.z),
+                self.bounds[2].expand_to_value(p.z()),
             ],
         }
     }
 
-    pub fn contains_point(&self, p: Point3<f64>) -> bool {
+    pub fn contains_point(&self, p: Vec3) -> bool {
         self.bounds
             .iter()
-            .zip(p.iter())
+            .zip(p.coords.iter())
             .all(|(interval, &value)| interval.contains_value(value))
     }
 
@@ -106,23 +107,23 @@ mod tests {
 
     #[test]
     fn from_corners_with_same_point_yields_degenerate_intervals() {
-        let test_point = Point3::new(0f64, 1f64, 2f64);
+        let test_point = Vec3::new(0f64, 1f64, 2f64);
         let target = BoundingBox::from_corners(test_point, test_point);
         assert!(target.bounds.iter().all(|e| e.is_degenerate()));
     }
 
     #[test]
     fn from_corners_yields_same_result_with_any_oposite_corners() {
-        let corner_000 = Point3::new(0.0, 0.0, 0.0);
+        let corner_000 = Vec3::new(0.0, 0.0, 0.0);
-        let corner_001 = Point3::new(0.0, 0.0, 1.0);
+        let corner_001 = Vec3::new(0.0, 0.0, 1.0);
-        let corner_010 = Point3::new(0.0, 1.0, 0.0);
+        let corner_010 = Vec3::new(0.0, 1.0, 0.0);
-        let corner_011 = Point3::new(0.0, 1.0, 1.0);
+        let corner_011 = Vec3::new(0.0, 1.0, 1.0);
-        let corner_100 = Point3::new(1.0, 0.0, 0.0);
+        let corner_100 = Vec3::new(1.0, 0.0, 0.0);
-        let corner_101 = Point3::new(1.0, 0.0, 1.0);
+        let corner_101 = Vec3::new(1.0, 0.0, 1.0);
-        let corner_110 = Point3::new(1.0, 1.0, 0.0);
+        let corner_110 = Vec3::new(1.0, 1.0, 0.0);
-        let corner_111 = Point3::new(1.0, 1.0, 1.0);
+        let corner_111 = Vec3::new(1.0, 1.0, 1.0);
 
-        let test_inputs: Vec<(Point3<f64>, Point3<f64>)> = vec![
+        let test_inputs: Vec<(Vec3, Vec3)> = vec![
             (corner_000, corner_111),
             (corner_001, corner_110),
             (corner_010, corner_101),
@@ -142,7 +143,7 @@ mod tests {
     }
 
     #[quickcheck]
-    fn union_with_self_yields_self(a: Point3<f64>, b: Point3<f64>) -> bool {
+    fn union_with_self_yields_self(a: Vec3, b: Vec3) -> bool {
         let target = BoundingBox::from_corners(a, b);
         let result = target.union(&target);
         target
@@ -153,12 +154,7 @@ mod tests {
     }
 
     #[quickcheck]
-    fn union_yields_full_ranges(
+    fn union_yields_full_ranges(a: Vec3, b: Vec3, c: Vec3, d: Vec3) -> bool {
-        a: Point3<f64>,
-        b: Point3<f64>,
-        c: Point3<f64>,
-        d: Point3<f64>,
-    ) -> bool {
         let target1 = BoundingBox::from_corners(a, b);
         let target2 = BoundingBox::from_corners(c, d);
         let result = target1.union(&target2);
@@ -176,28 +172,28 @@ mod tests {
     }
 
     #[quickcheck]
-    fn empty_box_contains_no_points(p: Point3<f64>) -> bool {
+    fn empty_box_contains_no_points(p: Vec3) -> bool {
         let target = BoundingBox::empty();
         !target.contains_point(p)
     }
 
     #[quickcheck]
     fn from_points_produces_box_that_contains_only_points_bounded_by_inputs_on_all_axes(
-        p: Point3<f64>,
+        p: Vec3,
-        a: Point3<f64>,
+        a: Vec3,
-        b: Point3<f64>,
+        b: Vec3,
-        c: Point3<f64>,
+        c: Vec3,
-        d: Point3<f64>,
+        d: Vec3,
-        e: Point3<f64>,
+        e: Vec3,
     ) -> bool {
         let points = vec![a, b, c, d, e];
         let target = BoundingBox::from_points(&points);
-        let is_in_bounds = points.iter().any(|elem| elem.x >= p.x)
+        let is_in_bounds = points.iter().any(|elem| elem.x() >= p.x())
-            && points.iter().any(|elem| elem.x <= p.x)
+            && points.iter().any(|elem| elem.x() <= p.x())
-            && points.iter().any(|elem| elem.y >= p.y)
+            && points.iter().any(|elem| elem.y() >= p.y())
-            && points.iter().any(|elem| elem.y <= p.y)
+            && points.iter().any(|elem| elem.y() <= p.y())
-            && points.iter().any(|elem| elem.z >= p.z)
+            && points.iter().any(|elem| elem.z() >= p.z())
-            && points.iter().any(|elem| elem.z <= p.z);
+            && points.iter().any(|elem| elem.z() <= p.z());
         target.contains_point(p) == is_in_bounds
     }
 

src/util/morton.rs

@@ -1,5 +1,5 @@
+use crate::math::Vec3;
 use crate::realtype::NormalizedToU32;
-use nalgebra::Point3;
 
 fn spread_bits(v: u32) -> u32 {
     let mut result = 0;
@@ -9,10 +9,10 @@ fn spread_bits(v: u32) -> u32 {
     result
 }
 
-pub fn morton_order_value_3d(p: Point3<f64>) -> u32 {
+pub fn morton_order_value_3d(p: Vec3) -> u32 {
-    let x = p.x.normalized_to_u32(10);
+    let x = p.x().normalized_to_u32(10);
-    let y = p.y.normalized_to_u32(10);
+    let y = p.y().normalized_to_u32(10);
-    let z = p.z.normalized_to_u32(10);
+    let z = p.z().normalized_to_u32(10);
     (spread_bits(x) << 2) | (spread_bits(y) << 1) | spread_bits(z)
 }
 

src/util/normalizer.rs

@@ -1,7 +1,7 @@
 use super::axis_aligned_bounding_box::BoundingBox;
 use super::Interval;
 
-use nalgebra::{clamp, Point3};
+use crate::math::Vec3;
 
 use itertools::izip;
 
@@ -23,7 +23,7 @@ impl RealNormalizer {
     }
 
     pub fn normalize_and_clamp(&self, value: f64) -> f64 {
-        clamp((value - self.min) / self.range, 0.0, 1.0)
+        ((value - self.min) / self.range).clamp(0.0, 1.0)
     }
 }
 
@@ -47,11 +47,11 @@ impl Point3Normalizer {
         normalizer
     }
 
-    pub fn normalize(&self, point: Point3<f64>) -> Point3<f64> {
+    pub fn normalize(&self, point: Vec3) -> Vec3 {
-        let mut result = Point3::new(0.0, 0.0, 0.0);
+        let mut result = Vec3::new(0.0, 0.0, 0.0);
         for (value_out, &value_in, normalizer) in izip!(
-            result.iter_mut(),
+            result.coords.iter_mut(),
-            point.iter(),
+            point.coords.iter(),
             self.dimension_normalizers.iter()
         ) {
             *value_out = normalizer.normalize(value_in);
@@ -59,11 +59,11 @@ impl Point3Normalizer {
         result
     }
 
-    pub fn normalize_and_clamp(&self, point: Point3<f64>) -> Point3<f64> {
+    pub fn normalize_and_clamp(&self, point: Vec3) -> Vec3 {
-        let mut result = Point3::new(0.0, 0.0, 0.0);
+        let mut result = Vec3::new(0.0, 0.0, 0.0);
         for (value_out, &value_in, normalizer) in izip!(
-            result.iter_mut(),
+            result.coords.iter_mut(),
-            point.iter(),
+            point.coords.iter(),
             self.dimension_normalizers.iter()
         ) {
             *value_out = normalizer.normalize_and_clamp(value_in);
@@ -145,34 +145,30 @@ mod test {
     }
 
     #[quickcheck]
-    fn normalize_point3_is_the_same_as_normalize_each_dimension(
+    fn normalize_point3_is_the_same_as_normalize_each_dimension(a: Vec3, b: Vec3, c: Vec3) -> bool {
-        a: Point3<f64>,
+        let x_normalizer = RealNormalizer::new(Interval::new(a.x().min(b.x()), a.x().max(b.x())));
-        b: Point3<f64>,
+        let y_normalizer = RealNormalizer::new(Interval::new(a.y().min(b.y()), a.y().max(b.y())));
-        c: Point3<f64>,
+        let z_normalizer = RealNormalizer::new(Interval::new(a.z().min(b.z()), a.z().max(b.z())));
-    ) -> bool {
-        let x_normalizer = RealNormalizer::new(Interval::new(a.x.min(b.x), a.x.max(b.x)));
-        let y_normalizer = RealNormalizer::new(Interval::new(a.y.min(b.y), a.y.max(b.y)));
-        let z_normalizer = RealNormalizer::new(Interval::new(a.z.min(b.z), a.z.max(b.z)));
         let xyz_normalizer = Point3Normalizer::new(BoundingBox::from_corners(a, b));
         let normalized_point = xyz_normalizer.normalize(c);
-        x_normalizer.normalize(c.x) == normalized_point.x
+        x_normalizer.normalize(c.x()) == normalized_point.x()
-            && y_normalizer.normalize(c.y) == normalized_point.y
+            && y_normalizer.normalize(c.y()) == normalized_point.y()
-            && z_normalizer.normalize(c.z) == normalized_point.z
+            && z_normalizer.normalize(c.z()) == normalized_point.z()
     }
 
     #[quickcheck]
     fn normalize_and_clamp_point3_is_the_same_as_normalize_and_clamp_each_dimension(
-        a: Point3<f64>,
+        a: Vec3,
-        b: Point3<f64>,
+        b: Vec3,
-        c: Point3<f64>,
+        c: Vec3,
     ) -> bool {
-        let x_normalizer = RealNormalizer::new(Interval::new(a.x.min(b.x), a.x.max(b.x)));
+        let x_normalizer = RealNormalizer::new(Interval::new(a.x().min(b.x()), a.x().max(b.x())));
-        let y_normalizer = RealNormalizer::new(Interval::new(a.y.min(b.y), a.y.max(b.y)));
+        let y_normalizer = RealNormalizer::new(Interval::new(a.y().min(b.y()), a.y().max(b.y())));
-        let z_normalizer = RealNormalizer::new(Interval::new(a.z.min(b.z), a.z.max(b.z)));
+        let z_normalizer = RealNormalizer::new(Interval::new(a.z().min(b.z()), a.z().max(b.z())));
-        let xyz_normalizer = dbg!(Point3Normalizer::new(BoundingBox::from_corners(a, b)));
+        let xyz_normalizer = Point3Normalizer::new(BoundingBox::from_corners(a, b));
         let normalized_point = xyz_normalizer.normalize_and_clamp(c);
-        x_normalizer.normalize_and_clamp(c.x) == normalized_point.x
+        x_normalizer.normalize_and_clamp(c.x()) == normalized_point.x()
-            && y_normalizer.normalize_and_clamp(c.y) == normalized_point.y
+            && y_normalizer.normalize_and_clamp(c.y()) == normalized_point.y()
-            && z_normalizer.normalize_and_clamp(c.z) == normalized_point.z
+            && z_normalizer.normalize_and_clamp(c.z()) == normalized_point.z()
     }
 }

src/util/polyhedra.rs

@@ -1,14 +1,14 @@
 use itertools::izip;
-use nalgebra::{convert, Point3, Vector3};
 
 use crate::materials::Material;
+use crate::math::Vec3;
 use crate::raycasting::Triangle;
 
 use std::sync::Arc;
 
 pub fn triangulate_polygon(
-    vertices: &Vec<Point3<f64>>,
+    vertices: &Vec<Vec3>,
-    normal: &Vector3<f64>,
+    normal: &Vec3,
     material: Arc<dyn Material>,
 ) -> Vec<Triangle> {
     assert!(vertices.len() >= 3);
@@ -23,97 +23,97 @@ pub fn triangulate_polygon(
 }
 
 pub fn generate_dodecahedron(
-    centre: Point3<f64>,
+    centre: Vec3,
     size: f64,
     material: Arc<dyn Material>,
 ) -> Vec<Triangle> {
-    let phi = convert((1.0 + (5.0_f64).sqrt()) / 2.0);
+    let phi = (1.0 + (5.0_f64).sqrt()) / 2.0;
     let phi_inv = 1.0 / phi;
 
     let faces = vec![
         vec![
-            Vector3::new(phi_inv, 0.0, phi),
+            Vec3::new(phi_inv, 0.0, phi),
-            Vector3::new(-phi_inv, 0.0, phi),
+            Vec3::new(-phi_inv, 0.0, phi),
-            Vector3::new(-1.0, -1.0, 1.0),
+            Vec3::new(-1.0, -1.0, 1.0),
-            Vector3::new(0.0, -phi, phi_inv),
+            Vec3::new(0.0, -phi, phi_inv),
-            Vector3::new(1.0, -1.0, 1.0),
+            Vec3::new(1.0, -1.0, 1.0),
         ],
         vec![
-            Vector3::new(phi_inv, 0.0, phi),
+            Vec3::new(phi_inv, 0.0, phi),
-            Vector3::new(-phi_inv, 0.0, phi),
+            Vec3::new(-phi_inv, 0.0, phi),
-            Vector3::new(-1.0, 1.0, 1.0),
+            Vec3::new(-1.0, 1.0, 1.0),
-            Vector3::new(0.0, phi, phi_inv),
+            Vec3::new(0.0, phi, phi_inv),
-            Vector3::new(1.0, 1.0, 1.0),
+            Vec3::new(1.0, 1.0, 1.0),
         ],
         vec![
-            Vector3::new(phi_inv, 0.0, phi),
+            Vec3::new(phi_inv, 0.0, phi),
-            Vector3::new(1.0, -1.0, 1.0),
+            Vec3::new(1.0, -1.0, 1.0),
-            Vector3::new(phi, -phi_inv, 0.0),
+            Vec3::new(phi, -phi_inv, 0.0),
-            Vector3::new(phi, phi_inv, 0.0),
+            Vec3::new(phi, phi_inv, 0.0),
-            Vector3::new(1.0, 1.0, 1.0),
+            Vec3::new(1.0, 1.0, 1.0),
         ],
         vec![
-            Vector3::new(-phi_inv, 0.0, phi),
+            Vec3::new(-phi_inv, 0.0, phi),
-            Vector3::new(-1.0, -1.0, 1.0),
+            Vec3::new(-1.0, -1.0, 1.0),
-            Vector3::new(-phi, -phi_inv, 0.0),
+            Vec3::new(-phi, -phi_inv, 0.0),
-            Vector3::new(-phi, phi_inv, 0.0),
+            Vec3::new(-phi, phi_inv, 0.0),
-            Vector3::new(-1.0, 1.0, 1.0),
+            Vec3::new(-1.0, 1.0, 1.0),
         ],
         vec![
-            Vector3::new(-1.0, -1.0, 1.0),
+            Vec3::new(-1.0, -1.0, 1.0),
-            Vector3::new(-phi, -phi_inv, 0.0),
+            Vec3::new(-phi, -phi_inv, 0.0),
-            Vector3::new(-1.0, -1.0, -1.0),
+            Vec3::new(-1.0, -1.0, -1.0),
-            Vector3::new(0.0, -phi, -phi_inv),
+            Vec3::new(0.0, -phi, -phi_inv),
-            Vector3::new(0.0, -phi, phi_inv),
+            Vec3::new(0.0, -phi, phi_inv),
         ],
         vec![
-            Vector3::new(0.0, -phi, phi_inv),
+            Vec3::new(0.0, -phi, phi_inv),
-            Vector3::new(0.0, -phi, -phi_inv),
+            Vec3::new(0.0, -phi, -phi_inv),
-            Vector3::new(1.0, -1.0, -1.0),
+            Vec3::new(1.0, -1.0, -1.0),
-            Vector3::new(phi, -phi_inv, 0.0),
+            Vec3::new(phi, -phi_inv, 0.0),
-            Vector3::new(1.0, -1.0, 1.0),
+            Vec3::new(1.0, -1.0, 1.0),
         ],
         vec![
-            Vector3::new(0.0, phi, phi_inv),
+            Vec3::new(0.0, phi, phi_inv),
-            Vector3::new(0.0, phi, -phi_inv),
+            Vec3::new(0.0, phi, -phi_inv),
-            Vector3::new(-1.0, 1.0, -1.0),
+            Vec3::new(-1.0, 1.0, -1.0),
-            Vector3::new(-phi, phi_inv, 0.0),
+            Vec3::new(-phi, phi_inv, 0.0),
-            Vector3::new(-1.0, 1.0, 1.0),
+            Vec3::new(-1.0, 1.0, 1.0),
         ],
         vec![
-            Vector3::new(1.0, 1.0, 1.0),
+            Vec3::new(1.0, 1.0, 1.0),
-            Vector3::new(phi, phi_inv, 0.0),
+            Vec3::new(phi, phi_inv, 0.0),
-            Vector3::new(1.0, 1.0, -1.0),
+            Vec3::new(1.0, 1.0, -1.0),
-            Vector3::new(0.0, phi, -phi_inv),
+            Vec3::new(0.0, phi, -phi_inv),
-            Vector3::new(0.0, phi, phi_inv),
+            Vec3::new(0.0, phi, phi_inv),
         ],
         vec![
-            Vector3::new(1.0, -1.0, -1.0),
+            Vec3::new(1.0, -1.0, -1.0),
-            Vector3::new(0.0, -phi, -phi_inv),
+            Vec3::new(0.0, -phi, -phi_inv),
-            Vector3::new(-1.0, -1.0, -1.0),
+            Vec3::new(-1.0, -1.0, -1.0),
-            Vector3::new(-phi_inv, 0.0, -phi),
+            Vec3::new(-phi_inv, 0.0, -phi),
-            Vector3::new(phi_inv, 0.0, -phi),
+            Vec3::new(phi_inv, 0.0, -phi),
         ],
         vec![
-            Vector3::new(1.0, 1.0, -1.0),
+            Vec3::new(1.0, 1.0, -1.0),
-            Vector3::new(0.0, phi, -phi_inv),
+            Vec3::new(0.0, phi, -phi_inv),
-            Vector3::new(-1.0, 1.0, -1.0),
+            Vec3::new(-1.0, 1.0, -1.0),
-            Vector3::new(-phi_inv, 0.0, -phi),
+            Vec3::new(-phi_inv, 0.0, -phi),
-            Vector3::new(phi_inv, 0.0, -phi),
+            Vec3::new(phi_inv, 0.0, -phi),
         ],
         vec![
-            Vector3::new(1.0, 1.0, -1.0),
+            Vec3::new(1.0, 1.0, -1.0),
-            Vector3::new(phi, phi_inv, 0.0),
+            Vec3::new(phi, phi_inv, 0.0),
-            Vector3::new(phi, -phi_inv, 0.0),
+            Vec3::new(phi, -phi_inv, 0.0),
-            Vector3::new(1.0, -1.0, -1.0),
+            Vec3::new(1.0, -1.0, -1.0),
-            Vector3::new(phi_inv, 0.0, -phi),
+            Vec3::new(phi_inv, 0.0, -phi),
         ],
         vec![
-            Vector3::new(-1.0, 1.0, -1.0),
+            Vec3::new(-1.0, 1.0, -1.0),
-            Vector3::new(-phi, phi_inv, 0.0),
+            Vec3::new(-phi, phi_inv, 0.0),
-            Vector3::new(-phi, -phi_inv, 0.0),
+            Vec3::new(-phi, -phi_inv, 0.0),
-            Vector3::new(-1.0, -1.0, -1.0),
+            Vec3::new(-1.0, -1.0, -1.0),
-            Vector3::new(-phi_inv, 0.0, -phi),
+            Vec3::new(-phi_inv, 0.0, -phi),
         ],
     ];