diff --git a/Cargo.lock b/Cargo.lock
index 20dada6..1199248 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -11,9 +11,9 @@ dependencies = [
 
 [[package]]
 name = "approx"
-version = "0.3.2"
+version = "0.4.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f0e60b75072ecd4168020818c0107f2857bb6c4e64252d8d3983f6263b40a5c3"
+checksum = "3f2a05fd1bd10b2527e20a2cd32d8873d115b8b39fe219ee25f42a8aca6ba278"
 dependencies = [
  "num-traits",
 ]
@@ -190,11 +190,9 @@ dependencies = [
 
 [[package]]
 name = "glam"
-version = "0.7.1"
-source = "git+https://github.com/bitshifter/glam-rs.git?rev=0f314f99#0f314f990710ff9357e5896de2b55ec82fe88e0d"
-dependencies = [
- "approx",
-]
+version = "0.15.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "411e0584defa447c328f25c756ba3d0685727ecc126b46c3c1176001141cd4b6"
 
 [[package]]
 name = "half"
@@ -614,9 +612,9 @@ checksum = "1d51f5df5af43ab3f1360b429fa5e0152ac5ce8c0bd6485cae490332e96846a8"
 
 [[package]]
 name = "sobol_burley"
-version = "0.2.0"
+version = "0.3.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b9d3ea39927d337a987e2334d9ff049617d2d132d64e3ce3cd6fb6e4aa00a864"
+checksum = "26e3528b09b1f1b1e152342a4462d1e80d568dc5623a0772252a6e584a53d550"
 
 [[package]]
 name = "spectral_upsampling"
diff --git a/Cargo.toml b/Cargo.toml
index 6cdcb23..10f1424 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -31,12 +31,12 @@ nom = "5"
 num_cpus = "1.8"
 openexr = "0.7"
 kioku = "0.3"
-sobol_burley = "0.2"
+sobol_burley = "0.3"
 png_encode_mini = "0.1.2"
 rustc-serialize = "0.3"
 scoped_threadpool = "0.1"
 time = "0.1"
-glam = {git="https://github.com/bitshifter/glam-rs.git", rev="0f314f99", default-features=false, features=["approx"]}
+glam = "0.15"
 fastapprox = "0.3"
 
 # Local crate dependencies
diff --git a/src/accel/bvh4.rs b/src/accel/bvh4.rs
index 211f476..3b78077 100644
--- a/src/accel/bvh4.rs
+++ b/src/accel/bvh4.rs
@@ -6,7 +6,7 @@
 
 use std::mem::{transmute, MaybeUninit};
 
-use glam::Vec4Mask;
+use glam::BVec4A;
 
 use kioku::Arena;
 
@@ -123,12 +123,12 @@ impl<'a> BVH4<'a> {
                     traversal_code,
                 } => {
                     node_tests += ray_stack.ray_count_in_next_task() as u64;
-                    let mut all_hits = Vec4Mask::default();
+                    let mut all_hits = BVec4A::default();
 
                     // Ray testing
                     ray_stack.pop_do_next_task_and_push_rays(children.len(), |ray_idx| {
                         if rays.is_done(ray_idx) {
-                            Vec4Mask::default()
+                            BVec4A::default()
                         } else {
                             let hits = if bounds.len() == 1 {
                                 bounds[0].intersect_ray(
diff --git a/src/bbox.rs b/src/bbox.rs
index b7706f2..f4a2ab6 100644
--- a/src/bbox.rs
+++ b/src/bbox.rs
@@ -7,7 +7,7 @@ use std::{
 
 use crate::{
     lerp::{lerp, lerp_slice, Lerp},
-    math::{fast_minf32, Matrix4x4, Point, Vector},
+    math::{fast_minf32, Point, Transform, Vector},
 };
 
 const BBOX_MAXT_ADJUST: f32 = 1.000_000_24;
@@ -41,8 +41,8 @@ impl BBox {
     // Returns whether the given ray intersects with the bbox.
     pub fn intersect_ray(&self, orig: Point, dir_inv: Vector, max_t: f32) -> bool {
         // Calculate slab intersections
-        let t1 = (self.min.co - orig.co).truncate() * dir_inv.co;
-        let t2 = (self.max.co - orig.co).truncate() * dir_inv.co;
+        let t1 = (self.min.co - orig.co) * dir_inv.co;
+        let t2 = (self.max.co - orig.co) * dir_inv.co;
 
         // Find the far and near intersection
         let far_t = t1.max(t2).extend(std::f32::INFINITY);
@@ -55,7 +55,7 @@ impl BBox {
     }
 
     // Creates a new BBox transformed into a different space.
-    pub fn transformed(&self, xform: Matrix4x4) -> BBox {
+    pub fn transformed(&self, xform: Transform) -> BBox {
         // BBox corners
         let vs = [
             Point::new(self.min.x(), self.min.y(), self.min.z()),
@@ -150,7 +150,7 @@ impl Lerp for BBox {
     }
 }
 
-pub fn transform_bbox_slice_from(bbs_in: &[BBox], xforms: &[Matrix4x4], bbs_out: &mut Vec<BBox>) {
+pub fn transform_bbox_slice_from(bbs_in: &[BBox], xforms: &[Transform], bbs_out: &mut Vec<BBox>) {
     bbs_out.clear();
 
     // Transform the bounding boxes
diff --git a/src/bbox4.rs b/src/bbox4.rs
index 7f80023..9464388 100644
--- a/src/bbox4.rs
+++ b/src/bbox4.rs
@@ -9,7 +9,7 @@ use crate::{
     math::{Point, Vector},
 };
 
-use glam::{Vec4, Vec4Mask};
+use glam::{BVec4A, Vec4};
 
 const BBOX_MAXT_ADJUST: f32 = 1.000_000_24;
 
@@ -60,14 +60,14 @@ impl BBox4 {
     }
 
     // Returns whether the given ray intersects with the bboxes.
-    pub fn intersect_ray(&self, orig: Point, dir_inv: Vector, max_t: f32) -> Vec4Mask {
+    pub fn intersect_ray(&self, orig: Point, dir_inv: Vector, max_t: f32) -> BVec4A {
         // Get the ray data into SIMD format.
-        let ro_x = Vec4::splat(orig.co.x());
-        let ro_y = Vec4::splat(orig.co.y());
-        let ro_z = Vec4::splat(orig.co.z());
-        let rdi_x = Vec4::splat(dir_inv.co.x());
-        let rdi_y = Vec4::splat(dir_inv.co.y());
-        let rdi_z = Vec4::splat(dir_inv.co.z());
+        let ro_x = Vec4::splat(orig.co[0]);
+        let ro_y = Vec4::splat(orig.co[1]);
+        let ro_z = Vec4::splat(orig.co[2]);
+        let rdi_x = Vec4::splat(dir_inv.co[0]);
+        let rdi_y = Vec4::splat(dir_inv.co[1]);
+        let rdi_z = Vec4::splat(dir_inv.co[2]);
         let max_t = Vec4::splat(max_t);
 
         // Slab tests
diff --git a/src/camera.rs b/src/camera.rs
index 296bb3d..788e207 100644
--- a/src/camera.rs
+++ b/src/camera.rs
@@ -4,14 +4,14 @@ use kioku::Arena;
 
 use crate::{
     lerp::lerp_slice,
-    math::{Matrix4x4, Point, Vector},
+    math::{Point, Transform, Vector},
     ray::Ray,
     sampling::square_to_circle,
 };
 
 #[derive(Copy, Clone, Debug)]
 pub struct Camera<'a> {
-    transforms: &'a [Matrix4x4],
+    transforms: &'a [Transform],
     fovs: &'a [f32],
     tfovs: &'a [f32],
     aperture_radii: &'a [f32],
@@ -21,7 +21,7 @@ pub struct Camera<'a> {
 impl<'a> Camera<'a> {
     pub fn new(
         arena: &'a Arena,
-        transforms: &[Matrix4x4],
+        transforms: &[Transform],
         fovs: &[f32],
         mut aperture_radii: &[f32],
         mut focus_distances: &[f32],
diff --git a/src/color.rs b/src/color.rs
index 4517df7..5fdb120 100644
--- a/src/color.rs
+++ b/src/color.rs
@@ -93,10 +93,10 @@ impl Color {
                 SpectralSample::from_parts(
                     // TODO: make this SIMD
                     Vec4::new(
-                        plancks_law(temperature, wls.x()) * factor,
-                        plancks_law(temperature, wls.y()) * factor,
-                        plancks_law(temperature, wls.z()) * factor,
-                        plancks_law(temperature, wls.w()) * factor,
+                        plancks_law(temperature, wls[0]) * factor,
+                        plancks_law(temperature, wls[1]) * factor,
+                        plancks_law(temperature, wls[2]) * factor,
+                        plancks_law(temperature, wls[3]) * factor,
                     ),
                     hero_wavelength,
                 )
@@ -108,10 +108,10 @@ impl Color {
                 SpectralSample::from_parts(
                     // TODO: make this SIMD
                     Vec4::new(
-                        plancks_law_normalized(temperature, wls.x()) * factor,
-                        plancks_law_normalized(temperature, wls.y()) * factor,
-                        plancks_law_normalized(temperature, wls.z()) * factor,
-                        plancks_law_normalized(temperature, wls.w()) * factor,
+                        plancks_law_normalized(temperature, wls[0]) * factor,
+                        plancks_law_normalized(temperature, wls[1]) * factor,
+                        plancks_law_normalized(temperature, wls[2]) * factor,
+                        plancks_law_normalized(temperature, wls[3]) * factor,
                     ),
                     hero_wavelength,
                 )
@@ -578,10 +578,10 @@ impl XYZ {
     }
 
     pub fn from_spectral_sample(ss: &SpectralSample) -> XYZ {
-        let xyz0 = XYZ::from_wavelength(ss.wl_n(0), ss.e.x());
-        let xyz1 = XYZ::from_wavelength(ss.wl_n(1), ss.e.y());
-        let xyz2 = XYZ::from_wavelength(ss.wl_n(2), ss.e.z());
-        let xyz3 = XYZ::from_wavelength(ss.wl_n(3), ss.e.w());
+        let xyz0 = XYZ::from_wavelength(ss.wl_n(0), ss.e[0]);
+        let xyz1 = XYZ::from_wavelength(ss.wl_n(1), ss.e[1]);
+        let xyz2 = XYZ::from_wavelength(ss.wl_n(2), ss.e[2]);
+        let xyz3 = XYZ::from_wavelength(ss.wl_n(3), ss.e[3]);
         (xyz0 + xyz1 + xyz2 + xyz3) * 0.75
     }
 
diff --git a/src/lerp.rs b/src/lerp.rs
index 9078b5c..65a6479 100644
--- a/src/lerp.rs
+++ b/src/lerp.rs
@@ -1,6 +1,6 @@
 #![allow(dead_code)]
 
-use math3d::{Matrix4x4, Normal, Point, Vector};
+use math3d::{Normal, Point, Transform, Vector};
 
 /// Trait for allowing a type to be linearly interpolated.
 pub trait Lerp: Copy {
@@ -106,8 +106,8 @@ impl Lerp for glam::Vec4 {
     }
 }
 
-impl Lerp for Matrix4x4 {
-    fn lerp(self, other: Matrix4x4, alpha: f32) -> Matrix4x4 {
+impl Lerp for Transform {
+    fn lerp(self, other: Transform, alpha: f32) -> Transform {
         (self * (1.0 - alpha)) + (other * alpha)
     }
 }
@@ -120,8 +120,8 @@ impl Lerp for Normal {
 
 impl Lerp for Point {
     fn lerp(self, other: Point, alpha: f32) -> Point {
-        let s = self.norm();
-        let o = other.norm();
+        let s = self;
+        let o = other;
         Point {
             co: (s.co * (1.0 - alpha)) + (o.co * alpha),
         }
@@ -215,23 +215,21 @@ mod tests {
 
     #[test]
     fn lerp_matrix() {
-        let a = Matrix4x4::new_from_values(
-            0.0, 2.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
+        let a = Transform::new_from_values(
+            0.0, 2.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0,
         );
-        let b = Matrix4x4::new_from_values(
-            -1.0, 1.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
+        let b = Transform::new_from_values(
+            -1.0, 1.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
         );
 
-        let c1 = Matrix4x4::new_from_values(
-            -0.25, 1.75, 2.25, 3.25, 4.25, 5.25, 6.25, 7.25, 8.25, 9.25, 10.25, 11.25, 12.25,
-            13.25, 14.25, 15.25,
+        let c1 = Transform::new_from_values(
+            -0.25, 1.75, 2.25, 3.25, 4.25, 5.25, 6.25, 7.25, 8.25, 9.25, 10.25, 11.25,
         );
-        let c2 = Matrix4x4::new_from_values(
-            -0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5,
+        let c2 = Transform::new_from_values(
+            -0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5,
         );
-        let c3 = Matrix4x4::new_from_values(
-            -0.75, 1.25, 2.75, 3.75, 4.75, 5.75, 6.75, 7.75, 8.75, 9.75, 10.75, 11.75, 12.75,
-            13.75, 14.75, 15.75,
+        let c3 = Transform::new_from_values(
+            -0.75, 1.25, 2.75, 3.75, 4.75, 5.75, 6.75, 7.75, 8.75, 9.75, 10.75, 11.75,
         );
 
         assert_eq!(a.lerp(b, 0.0), a);
diff --git a/src/light/mod.rs b/src/light/mod.rs
index e38f53a..a45c567 100644
--- a/src/light/mod.rs
+++ b/src/light/mod.rs
@@ -6,7 +6,7 @@ use std::fmt::Debug;
 
 use crate::{
     color::SpectralSample,
-    math::{Matrix4x4, Normal, Point, Vector},
+    math::{Normal, Point, Transform, Vector},
     surface::Surface,
 };
 
@@ -34,7 +34,7 @@ pub trait SurfaceLight: Surface {
     /// - The pdf of the sample.
     fn sample_from_point(
         &self,
-        space: &Matrix4x4,
+        space: &Transform,
         arr: Point,
         u: f32,
         v: f32,
diff --git a/src/light/rectangle_light.rs b/src/light/rectangle_light.rs
index 4d09552..2d7c9a3 100644
--- a/src/light/rectangle_light.rs
+++ b/src/light/rectangle_light.rs
@@ -5,7 +5,7 @@ use crate::{
     boundable::Boundable,
     color::{Color, SpectralSample},
     lerp::lerp_slice,
-    math::{cross, dot, Matrix4x4, Normal, Point, Vector},
+    math::{cross, dot, Normal, Point, Transform, Vector},
     ray::{RayBatch, RayStack},
     sampling::{
         spherical_triangle_solid_angle, triangle_surface_area, uniform_sample_spherical_triangle,
@@ -51,7 +51,7 @@ impl<'a> RectangleLight<'a> {
     // more efficiently by inlining it there.
     fn sample_pdf(
         &self,
-        space: &Matrix4x4,
+        space: &Transform,
         arr: Point,
         sample_dir: Vector,
         hit_point: Point,
@@ -97,7 +97,7 @@ impl<'a> RectangleLight<'a> {
 
     // fn outgoing(
     //     &self,
-    //     space: &Matrix4x4,
+    //     space: &Transform,
     //     dir: Vector,
     //     u: f32,
     //     v: f32,
@@ -120,7 +120,7 @@ impl<'a> RectangleLight<'a> {
 impl<'a> SurfaceLight for RectangleLight<'a> {
     fn sample_from_point(
         &self,
-        space: &Matrix4x4,
+        space: &Transform,
         arr: Point,
         u: f32,
         v: f32,
@@ -261,7 +261,7 @@ impl<'a> Surface for RectangleLight<'a> {
         ray_stack: &mut RayStack,
         isects: &mut [SurfaceIntersection],
         shader: &dyn SurfaceShader,
-        space: &[Matrix4x4],
+        space: &[Transform],
     ) {
         let _ = shader; // Silence 'unused' warning
 
diff --git a/src/light/sphere_light.rs b/src/light/sphere_light.rs
index cfa1635..079ff0a 100644
--- a/src/light/sphere_light.rs
+++ b/src/light/sphere_light.rs
@@ -7,7 +7,7 @@ use crate::{
     boundable::Boundable,
     color::{Color, SpectralSample},
     lerp::lerp_slice,
-    math::{coordinate_system_from_vector, dot, Matrix4x4, Normal, Point, Vector},
+    math::{coordinate_system_from_vector, dot, Normal, Point, Transform, Vector},
     ray::{RayBatch, RayStack},
     sampling::{uniform_sample_cone, uniform_sample_cone_pdf, uniform_sample_sphere},
     shading::surface_closure::SurfaceClosure,
@@ -50,7 +50,7 @@ impl<'a> SphereLight<'a> {
     // more efficiently by inlining it there.
     fn sample_pdf(
         &self,
-        space: &Matrix4x4,
+        space: &Transform,
         arr: Point,
         sample_dir: Vector,
         sample_u: f32,
@@ -84,7 +84,7 @@ impl<'a> SphereLight<'a> {
 impl<'a> SurfaceLight for SphereLight<'a> {
     fn sample_from_point(
         &self,
-        space: &Matrix4x4,
+        space: &Transform,
         arr: Point,
         u: f32,
         v: f32,
@@ -210,7 +210,7 @@ impl<'a> Surface for SphereLight<'a> {
         ray_stack: &mut RayStack,
         isects: &mut [SurfaceIntersection],
         shader: &dyn SurfaceShader,
-        space: &[Matrix4x4],
+        space: &[Transform],
     ) {
         let _ = shader; // Silence 'unused' warning
 
diff --git a/src/math.rs b/src/math.rs
index 6df51e4..fec2f06 100644
--- a/src/math.rs
+++ b/src/math.rs
@@ -2,7 +2,7 @@
 
 use std::f32;
 
-pub use math3d::{cross, dot, CrossProduct, DotProduct, Matrix4x4, Normal, Point, Vector};
+pub use math3d::{cross, dot, CrossProduct, DotProduct, Normal, Point, Transform, Vector};
 
 /// Clamps a value between a min and max.
 pub fn clamp<T: PartialOrd>(v: T, lower: T, upper: T) -> T {
diff --git a/src/parse/psy.rs b/src/parse/psy.rs
index 6b2bd0f..32c30a8 100644
--- a/src/parse/psy.rs
+++ b/src/parse/psy.rs
@@ -11,7 +11,7 @@ use crate::{
     camera::Camera,
     color::{rec709_e_to_xyz, Color},
     light::WorldLightSource,
-    math::Matrix4x4,
+    math::Transform,
     // renderer::Renderer,
     scene::Scene,
     scene::World,
@@ -650,16 +650,17 @@ fn parse_shaders<'a>(
     return Ok(shaders);
 }
 
-pub fn parse_matrix(contents: &str) -> PsyResult<Matrix4x4> {
+pub fn parse_matrix(contents: &str) -> PsyResult<Transform> {
     if let IResult::Ok((leftover, ns)) = all_consuming(tuple((
         ws_f32, ws_f32, ws_f32, ws_f32, ws_f32, ws_f32, ws_f32, ws_f32, ws_f32, ws_f32, ws_f32,
         ws_f32, ws_f32, ws_f32, ws_f32, ws_f32,
     )))(contents)
     {
         if leftover.is_empty() {
-            return Ok(Matrix4x4::new_from_values(
-                ns.0, ns.4, ns.8, ns.12, ns.1, ns.5, ns.9, ns.13, ns.2, ns.6, ns.10, ns.14, ns.3,
-                ns.7, ns.11, ns.15,
+            return Ok(Transform::new_from_values(
+                // We throw away the last row, since it's not necessarily affine.
+                // TODO: is there a more correct way to handle this?
+                ns.0, ns.4, ns.8, ns.12, ns.1, ns.5, ns.9, ns.13, ns.2, ns.6, ns.10, ns.14,
             ));
         }
     }
diff --git a/src/parse/psy_mesh_surface.rs b/src/parse/psy_mesh_surface.rs
index 069df75..47e9f5c 100644
--- a/src/parse/psy_mesh_surface.rs
+++ b/src/parse/psy_mesh_surface.rs
@@ -9,7 +9,7 @@ use kioku::Arena;
 use data_tree::{DataTreeReader, Event};
 
 use crate::{
-    math::{Matrix4x4, Normal, Point},
+    math::{Normal, Point, Transform},
     surface::triangle_mesh::TriangleMesh,
 };
 
@@ -21,7 +21,7 @@ use super::{
 pub fn parse_mesh_surface<'a>(
     arena: &'a Arena,
     instance_xform_idxs: &[std::ops::Range<usize>],
-    xforms: &[Matrix4x4],
+    xforms: &[Transform],
     events: &mut DataTreeReader<impl BufRead>,
 ) -> PsyResult<TriangleMesh<'a>> {
     let mut verts = Vec::new(); // Vec of vecs, one for each time sample
diff --git a/src/ray.rs b/src/ray.rs
index 2ad37b8..daf29ab 100644
--- a/src/ray.rs
+++ b/src/ray.rs
@@ -1,8 +1,8 @@
 #![allow(dead_code)]
 
-use glam::Vec4Mask;
+use glam::BVec4A;
 
-use crate::math::{Matrix4x4, Point, Vector};
+use crate::math::{Point, Transform, Vector};
 
 type RayIndexType = u16;
 type FlagType = u8;
@@ -86,7 +86,7 @@ impl RayBatch {
     pub fn set_from_ray(&mut self, ray: &Ray, is_occlusion: bool, idx: usize) {
         self.hot[idx].orig_local = ray.orig;
         self.hot[idx].dir_inv_local = Vector {
-            co: ray.dir.co.reciprocal(),
+            co: ray.dir.co.recip(),
         };
         self.hot[idx].max_t = ray.max_t;
         self.hot[idx].time = ray.time;
@@ -119,10 +119,10 @@ impl RayBatch {
     ///
     /// This should be called when entering (and exiting) traversal of a
     /// new transform space.
-    pub fn update_local(&mut self, idx: usize, xform: &Matrix4x4) {
+    pub fn update_local(&mut self, idx: usize, xform: &Transform) {
         self.hot[idx].orig_local = self.cold[idx].orig * *xform;
         self.hot[idx].dir_inv_local = Vector {
-            co: (self.cold[idx].dir * *xform).co.reciprocal(),
+            co: (self.cold[idx].dir * *xform).co.recip(),
         };
     }
 
@@ -349,7 +349,7 @@ impl RayStack {
     /// indicated lanes.
     pub fn pop_do_next_task_and_push_rays<F>(&mut self, output_lane_count: usize, mut handle_ray: F)
     where
-        F: FnMut(usize) -> Vec4Mask,
+        F: FnMut(usize) -> BVec4A,
     {
         // Pop the task and do necessary bookkeeping.
         let task = self.tasks.pop().unwrap();
diff --git a/src/renderer.rs b/src/renderer.rs
index e2d8483..9cbd1e6 100644
--- a/src/renderer.rs
+++ b/src/renderer.rs
@@ -706,7 +706,7 @@ fn get_sample_4d(
     let seed = pixel_co.0 ^ (pixel_co.1 << 16) ^ seed.wrapping_mul(0x736caf6f);
 
     match dimension_set {
-        ds if ds < sobol_burley::NUM_DIMENSION_SETS as u32 => {
+        ds if ds < sobol_burley::NUM_DIMENSION_SETS_4D as u32 => {
             // Sobol sampling.
             let n4 = sobol_burley::sample_4d(i, ds, seed);
             (n4[0], n4[1], n4[2], n4[3])
diff --git a/src/scene/assembly.rs b/src/scene/assembly.rs
index 87c177d..77facd5 100644
--- a/src/scene/assembly.rs
+++ b/src/scene/assembly.rs
@@ -1,12 +1,12 @@
 use std::{collections::HashMap, ops::Range};
 
-use crate::{light::SurfaceLight, math::Matrix4x4, surface::Surface};
+use crate::{light::SurfaceLight, math::Transform, surface::Surface};
 
 /// Stores the objects of a scene and its acceleration structures.
 #[derive(Debug)]
 pub struct Assembly<'a> {
     pub objects: HashMap<String, Object<'a>>, // Name, Object.
-    pub xforms: Vec<Matrix4x4>,
+    pub xforms: Vec<Transform>,
 }
 
 impl<'a> Assembly<'a> {
diff --git a/src/shading/surface_closure.rs b/src/shading/surface_closure.rs
index a769e27..50d7e5f 100644
--- a/src/shading/surface_closure.rs
+++ b/src/shading/surface_closure.rs
@@ -545,23 +545,23 @@ mod ggx_closure {
             let spectrum_sample = col.to_spectral_sample(wavelength);
             let rev_fresnel = 1.0 - fresnel;
             let c0 = lerp(
-                schlick_fresnel_from_fac(spectrum_sample.e.x(), hb),
-                spectrum_sample.e.x(),
+                schlick_fresnel_from_fac(spectrum_sample.e[0], hb),
+                spectrum_sample.e[0],
                 rev_fresnel,
             );
             let c1 = lerp(
-                schlick_fresnel_from_fac(spectrum_sample.e.y(), hb),
-                spectrum_sample.e.y(),
+                schlick_fresnel_from_fac(spectrum_sample.e[1], hb),
+                spectrum_sample.e[1],
                 rev_fresnel,
             );
             let c2 = lerp(
-                schlick_fresnel_from_fac(spectrum_sample.e.z(), hb),
-                spectrum_sample.e.z(),
+                schlick_fresnel_from_fac(spectrum_sample.e[2], hb),
+                spectrum_sample.e[2],
                 rev_fresnel,
             );
             let c3 = lerp(
-                schlick_fresnel_from_fac(spectrum_sample.e.w(), hb),
-                spectrum_sample.e.w(),
+                schlick_fresnel_from_fac(spectrum_sample.e[3], hb),
+                spectrum_sample.e[3],
                 rev_fresnel,
             );
 
diff --git a/src/surface/micropoly_batch.rs b/src/surface/micropoly_batch.rs
index 3dbb75d..67504fd 100644
--- a/src/surface/micropoly_batch.rs
+++ b/src/surface/micropoly_batch.rs
@@ -9,7 +9,7 @@ use crate::{
     bbox::BBox,
     boundable::Boundable,
     lerp::lerp_slice,
-    math::{cross, dot, Matrix4x4, Normal, Point},
+    math::{cross, dot, Normal, Point, Transform},
     ray::{RayBatch, RayStack},
     shading::SurfaceClosure,
 };
@@ -149,13 +149,13 @@ impl<'a> MicropolyBatch<'a> {
         rays: &mut RayBatch,
         ray_stack: &mut RayStack,
         isects: &mut [SurfaceIntersection],
-        space: &[Matrix4x4],
+        space: &[Transform],
     ) {
         // Precalculate transform for non-motion blur cases
         let static_mat_space = if space.len() == 1 {
             lerp_slice(space, 0.0).inverse()
         } else {
-            Matrix4x4::new()
+            Transform::new()
         };
 
         self.accel
diff --git a/src/surface/mod.rs b/src/surface/mod.rs
index 7122261..a718f4d 100644
--- a/src/surface/mod.rs
+++ b/src/surface/mod.rs
@@ -10,7 +10,7 @@ use std::fmt::Debug;
 
 use crate::{
     boundable::Boundable,
-    math::{Matrix4x4, Normal, Point, Vector},
+    math::{Normal, Point, Transform, Vector},
     ray::{RayBatch, RayStack},
     shading::surface_closure::SurfaceClosure,
     shading::SurfaceShader,
@@ -25,7 +25,7 @@ pub trait Surface: Boundable + Debug + Sync {
         ray_stack: &mut RayStack,
         isects: &mut [SurfaceIntersection],
         shader: &dyn SurfaceShader,
-        space: &[Matrix4x4],
+        space: &[Transform],
     );
 }
 
@@ -86,7 +86,7 @@ pub struct SurfaceIntersectionData {
     // a cube centered around `pos` with dimensions of `2 * pos_err`.
     pub nor: Normal,            // Shading normal
     pub nor_g: Normal,          // True geometric normal
-    pub local_space: Matrix4x4, // Matrix from global space to local space
+    pub local_space: Transform, // Matrix from global space to local space
     pub t: f32,                 // Ray t-value at the intersection point
     pub sample_pdf: f32,        // The PDF of getting this point by explicitly sampling the surface
 }
diff --git a/src/surface/triangle_mesh.rs b/src/surface/triangle_mesh.rs
index be1ea8f..6b16ab6 100644
--- a/src/surface/triangle_mesh.rs
+++ b/src/surface/triangle_mesh.rs
@@ -7,7 +7,7 @@ use crate::{
     bbox::BBox,
     boundable::Boundable,
     lerp::lerp_slice,
-    math::{cross, dot, Matrix4x4, Normal, Point},
+    math::{cross, dot, Normal, Point, Transform},
     ray::{RayBatch, RayStack},
     shading::SurfaceShader,
 };
@@ -128,13 +128,13 @@ impl<'a> Surface for TriangleMesh<'a> {
         ray_stack: &mut RayStack,
         isects: &mut [SurfaceIntersection],
         shader: &dyn SurfaceShader,
-        space: &[Matrix4x4],
+        space: &[Transform],
     ) {
         // Precalculate transform for non-motion blur cases
         let static_mat_space = if space.len() == 1 {
             lerp_slice(space, 0.0).inverse()
         } else {
-            Matrix4x4::new()
+            Transform::new()
         };
 
         self.accel
diff --git a/src/tracer.rs b/src/tracer.rs
index 95cf8cb..8e85db6 100644
--- a/src/tracer.rs
+++ b/src/tracer.rs
@@ -4,7 +4,7 @@ use crate::{
     accel::ray_code,
     color::{rec709_to_xyz, Color},
     lerp::lerp_slice,
-    math::Matrix4x4,
+    math::Transform,
     ray::{RayBatch, RayStack},
     scene::{Assembly, InstanceType, Object},
     shading::{SimpleSurfaceShader, SurfaceShader},
@@ -63,7 +63,7 @@ impl<'a> TracerInner<'a> {
 
         // Prep the accel part of the rays.
         {
-            let ident = Matrix4x4::new();
+            let ident = Transform::new();
             for i in 0..rays.len() {
                 rays.update_local(i, &ident);
             }
@@ -140,7 +140,7 @@ impl<'a> TracerInner<'a> {
                             rays.update_local(ray_idx, &lerp_slice(xforms, t));
                         });
                     } else {
-                        let ident = Matrix4x4::new();
+                        let ident = Transform::new();
                         ray_stack.pop_do_next_task(|ray_idx| {
                             rays.update_local(ray_idx, &ident);
                         });
diff --git a/src/transform_stack.rs b/src/transform_stack.rs
index 2612ce1..8219799 100644
--- a/src/transform_stack.rs
+++ b/src/transform_stack.rs
@@ -3,10 +3,10 @@ use std::{
     mem::{transmute, MaybeUninit},
 };
 
-use crate::{algorithm::merge_slices_to, math::Matrix4x4};
+use crate::{algorithm::merge_slices_to, math::Transform};
 
 pub struct TransformStack {
-    stack: Vec<MaybeUninit<Matrix4x4>>,
+    stack: Vec<MaybeUninit<Transform>>,
     stack_indices: Vec<usize>,
 }
 
@@ -30,11 +30,11 @@ impl TransformStack {
         self.stack_indices.push(0);
     }
 
-    pub fn push(&mut self, xforms: &[Matrix4x4]) {
+    pub fn push(&mut self, xforms: &[Transform]) {
         assert!(!xforms.is_empty());
 
         if self.stack.is_empty() {
-            let xforms: &[MaybeUninit<Matrix4x4>] = unsafe { transmute(xforms) };
+            let xforms: &[MaybeUninit<Transform>] = unsafe { transmute(xforms) };
             self.stack.extend(xforms);
         } else {
             let sil = self.stack_indices.len();
@@ -73,7 +73,7 @@ impl TransformStack {
         self.stack_indices.pop();
     }
 
-    pub fn top(&self) -> &[Matrix4x4] {
+    pub fn top(&self) -> &[Transform] {
         let sil = self.stack_indices.len();
         let i1 = self.stack_indices[sil - 2];
         let i2 = self.stack_indices[sil - 1];
diff --git a/sub_crates/math3d/Cargo.toml b/sub_crates/math3d/Cargo.toml
index 0f69425..dba6e93 100644
--- a/sub_crates/math3d/Cargo.toml
+++ b/sub_crates/math3d/Cargo.toml
@@ -11,5 +11,5 @@ path = "src/lib.rs"
 
 # Local crate dependencies
 [dependencies]
-glam = {git="https://github.com/bitshifter/glam-rs.git", rev="0f314f99", default-features=false, features=["approx"]}
-approx = "0.3"
+glam = "0.15"
+approx = "0.4"
diff --git a/sub_crates/math3d/src/lib.rs b/sub_crates/math3d/src/lib.rs
index 6795dee..3c51eed 100644
--- a/sub_crates/math3d/src/lib.rs
+++ b/sub_crates/math3d/src/lib.rs
@@ -1,11 +1,11 @@
 #![allow(dead_code)]
 
-mod matrix;
 mod normal;
 mod point;
+mod transform;
 mod vector;
 
-pub use self::{matrix::Matrix4x4, normal::Normal, point::Point, vector::Vector};
+pub use self::{normal::Normal, point::Point, transform::Transform, vector::Vector};
 
 /// Trait for calculating dot products.
 pub trait DotProduct {
diff --git a/sub_crates/math3d/src/matrix.rs b/sub_crates/math3d/src/matrix.rs
deleted file mode 100644
index e804064..0000000
--- a/sub_crates/math3d/src/matrix.rs
+++ /dev/null
@@ -1,193 +0,0 @@
-#![allow(dead_code)]
-
-use std::ops::{Add, Mul};
-
-use approx::RelativeEq;
-use glam::{Mat4, Vec4};
-
-use super::Point;
-
-/// A 4x4 matrix, used for transforms
-#[derive(Debug, Copy, Clone, PartialEq)]
-pub struct Matrix4x4(pub Mat4);
-
-impl Matrix4x4 {
-    /// Creates a new identity matrix
-    #[inline]
-    pub fn new() -> Matrix4x4 {
-        Matrix4x4(Mat4::identity())
-    }
-
-    /// Creates a new matrix with the specified values:
-    /// a b c d
-    /// e f g h
-    /// i j k l
-    /// m n o p
-    #[inline]
-    #[allow(clippy::many_single_char_names)]
-    #[allow(clippy::too_many_arguments)]
-    pub fn new_from_values(
-        a: f32,
-        b: f32,
-        c: f32,
-        d: f32,
-        e: f32,
-        f: f32,
-        g: f32,
-        h: f32,
-        i: f32,
-        j: f32,
-        k: f32,
-        l: f32,
-        m: f32,
-        n: f32,
-        o: f32,
-        p: f32,
-    ) -> Matrix4x4 {
-        Matrix4x4(Mat4::new(
-            Vec4::new(a, e, i, m),
-            Vec4::new(b, f, j, n),
-            Vec4::new(c, g, k, o),
-            Vec4::new(d, h, l, p),
-        ))
-    }
-
-    #[inline]
-    pub fn from_location(loc: Point) -> Matrix4x4 {
-        Matrix4x4(Mat4::from_translation(loc.co.truncate()))
-    }
-
-    /// Returns whether the matrices are approximately equal to each other.
-    /// Each corresponding element in the matrices cannot have a relative
-    /// error exceeding epsilon.
-    #[inline]
-    pub fn aprx_eq(&self, other: Matrix4x4, epsilon: f32) -> bool {
-        self.0.relative_eq(&other.0, std::f32::EPSILON, epsilon)
-    }
-
-    /// Returns the transpose of the matrix
-    #[inline]
-    pub fn transposed(&self) -> Matrix4x4 {
-        Matrix4x4(self.0.transpose())
-    }
-
-    /// Returns the inverse of the Matrix
-    #[inline]
-    pub fn inverse(&self) -> Matrix4x4 {
-        Matrix4x4(self.0.inverse())
-    }
-}
-
-impl Default for Matrix4x4 {
-    fn default() -> Self {
-        Self::new()
-    }
-}
-
-/// Multiply two matrices together
-impl Mul for Matrix4x4 {
-    type Output = Self;
-
-    #[inline]
-    fn mul(self, other: Self) -> Self {
-        Self(other.0.mul_mat4(&self.0))
-    }
-}
-
-/// Multiply a matrix by a f32
-impl Mul<f32> for Matrix4x4 {
-    type Output = Self;
-
-    #[inline]
-    fn mul(self, other: f32) -> Self {
-        Self(self.0 * other)
-    }
-}
-
-/// Add two matrices together
-impl Add for Matrix4x4 {
-    type Output = Self;
-
-    #[inline]
-    fn add(self, other: Self) -> Self {
-        Self(self.0 + other.0)
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    #[test]
-    fn equality_test() {
-        let a = Matrix4x4::new();
-        let b = Matrix4x4::new();
-        let c = Matrix4x4::new_from_values(
-            1.1, 0.0, 0.0, 0.0, 0.0, 1.1, 0.0, 0.0, 0.0, 0.0, 1.1, 0.0, 0.0, 0.0, 0.0, 1.1,
-        );
-
-        assert_eq!(a, b);
-        assert!(a != c);
-    }
-
-    #[test]
-    fn approximate_equality_test() {
-        let a = Matrix4x4::new();
-        let b = Matrix4x4::new_from_values(
-            1.000001, 0.0, 0.0, 0.0, 0.0, 1.000001, 0.0, 0.0, 0.0, 0.0, 1.000001, 0.0, 0.0, 0.0,
-            0.0, 1.000001,
-        );
-        let c = Matrix4x4::new_from_values(
-            1.000003, 0.0, 0.0, 0.0, 0.0, 1.000003, 0.0, 0.0, 0.0, 0.0, 1.000003, 0.0, 0.0, 0.0,
-            0.0, 1.000003,
-        );
-        let d = Matrix4x4::new_from_values(
-            -1.000001, 0.0, 0.0, 0.0, 0.0, -1.000001, 0.0, 0.0, 0.0, 0.0, -1.000001, 0.0, 0.0, 0.0,
-            0.0, -1.000001,
-        );
-
-        assert!(a.aprx_eq(b, 0.000001));
-        assert!(!a.aprx_eq(c, 0.000001));
-        assert!(!a.aprx_eq(d, 0.000001));
-    }
-
-    #[test]
-    fn multiply_test() {
-        let a = Matrix4x4::new_from_values(
-            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0, 13.0, 7.0, 15.0, 3.0,
-        );
-        let b = Matrix4x4::new_from_values(
-            1.0, 5.0, 9.0, 13.0, 2.0, 6.0, 10.0, 14.0, 3.0, 7.0, 11.0, 15.0, 4.0, 8.0, 12.0, 16.0,
-        );
-        let c = Matrix4x4::new_from_values(
-            266.0, 141.0, 331.0, 188.5, 292.0, 158.0, 366.0, 213.0, 318.0, 175.0, 401.0, 237.5,
-            344.0, 192.0, 436.0, 262.0,
-        );
-
-        assert_eq!(a * b, c);
-    }
-
-    #[test]
-    fn inverse_test() {
-        let a = Matrix4x4::new_from_values(
-            1.0, 0.33, 0.0, -2.0, 0.0, 1.0, 0.0, 0.0, 2.1, 0.7, 1.3, 0.0, 0.0, 0.0, 0.0, -1.0,
-        );
-        let b = a.inverse();
-        let c = Matrix4x4::new();
-
-        assert!((dbg!(a * b)).aprx_eq(dbg!(c), 0.0000001));
-    }
-
-    #[test]
-    fn transpose_test() {
-        let a = Matrix4x4::new_from_values(
-            1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
-        );
-        let b = Matrix4x4::new_from_values(
-            1.0, 5.0, 9.0, 13.0, 2.0, 6.0, 10.0, 14.0, 3.0, 7.0, 11.0, 15.0, 4.0, 8.0, 12.0, 16.0,
-        );
-        let c = a.transposed();
-
-        assert_eq!(b, c);
-    }
-}
diff --git a/sub_crates/math3d/src/normal.rs b/sub_crates/math3d/src/normal.rs
index c02bac9..c5d55a7 100644
--- a/sub_crates/math3d/src/normal.rs
+++ b/sub_crates/math3d/src/normal.rs
@@ -5,21 +5,21 @@ use std::{
     ops::{Add, Div, Mul, Neg, Sub},
 };
 
-use glam::Vec3;
+use glam::Vec3A;
 
-use super::{CrossProduct, DotProduct, Matrix4x4, Vector};
+use super::{CrossProduct, DotProduct, Transform, Vector};
 
 /// A surface normal in 3d homogeneous space.
 #[derive(Debug, Copy, Clone)]
 pub struct Normal {
-    pub co: Vec3,
+    pub co: Vec3A,
 }
 
 impl Normal {
     #[inline(always)]
     pub fn new(x: f32, y: f32, z: f32) -> Normal {
         Normal {
-            co: Vec3::new(x, y, z),
+            co: Vec3A::new(x, y, z),
         }
     }
 
@@ -57,32 +57,32 @@ impl Normal {
 
     #[inline(always)]
     pub fn x(&self) -> f32 {
-        self.co.x()
+        self.co[0]
     }
 
     #[inline(always)]
     pub fn y(&self) -> f32 {
-        self.co.y()
+        self.co[1]
     }
 
     #[inline(always)]
     pub fn z(&self) -> f32 {
-        self.co.z()
+        self.co[2]
     }
 
     #[inline(always)]
     pub fn set_x(&mut self, x: f32) {
-        self.co.set_x(x);
+        self.co[0] = x;
     }
 
     #[inline(always)]
     pub fn set_y(&mut self, y: f32) {
-        self.co.set_y(y);
+        self.co[1] = y;
     }
 
     #[inline(always)]
     pub fn set_z(&mut self, z: f32) {
-        self.co.set_z(z);
+        self.co[2] = z;
     }
 }
 
@@ -126,14 +126,13 @@ impl Mul<f32> for Normal {
     }
 }
 
-impl Mul<Matrix4x4> for Normal {
+impl Mul<Transform> for Normal {
     type Output = Normal;
 
     #[inline]
-    fn mul(self, other: Matrix4x4) -> Normal {
-        let mat = other.0.inverse().transpose();
+    fn mul(self, other: Transform) -> Normal {
         Normal {
-            co: mat.transform_vector3(self.co),
+            co: other.0.matrix3.inverse().transpose().mul_vec3a(self.co),
         }
     }
 }
@@ -176,9 +175,9 @@ impl CrossProduct for Normal {
 
 #[cfg(test)]
 mod tests {
-    use super::super::{CrossProduct, DotProduct, Matrix4x4};
+    use super::super::{CrossProduct, DotProduct, Transform};
     use super::*;
-    use approx::UlpsEq;
+    use approx::assert_ulps_eq;
 
     #[test]
     fn add() {
@@ -210,12 +209,14 @@ mod tests {
     #[test]
     fn mul_matrix_1() {
         let n = Normal::new(1.0, 2.5, 4.0);
-        let m = Matrix4x4::new_from_values(
-            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0, 13.0, 7.0, 15.0, 3.0,
+        let m = Transform::new_from_values(
+            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0,
         );
         let nm = n * m;
-        let nm2 = Normal::new(-19.258825, 5.717648, -1.770588);
-        assert!(nm.co.ulps_eq(&nm2.co, 0.0, 4));
+        let nm2 = Normal::new(-4.0625, 1.78125, -0.03125);
+        for i in 0..3 {
+            assert_ulps_eq!(nm.co[i], nm2.co[i], max_ulps = 4);
+        }
     }
 
     #[test]
diff --git a/sub_crates/math3d/src/point.rs b/sub_crates/math3d/src/point.rs
index 33d993c..19deb56 100644
--- a/sub_crates/math3d/src/point.rs
+++ b/sub_crates/math3d/src/point.rs
@@ -5,37 +5,28 @@ use std::{
     ops::{Add, Mul, Sub},
 };
 
-use glam::Vec4;
+use glam::Vec3A;
 
-use super::{Matrix4x4, Vector};
+use super::{Transform, Vector};
 
 /// A position in 3d homogeneous space.
 #[derive(Debug, Copy, Clone)]
 pub struct Point {
-    pub co: Vec4,
+    pub co: Vec3A,
 }
 
 impl Point {
     #[inline(always)]
     pub fn new(x: f32, y: f32, z: f32) -> Point {
         Point {
-            co: Vec4::new(x, y, z, 1.0),
-        }
-    }
-
-    /// Returns the point in standardized coordinates, where the
-    /// fourth homogeneous component has been normalized to 1.0.
-    #[inline(always)]
-    pub fn norm(&self) -> Point {
-        Point {
-            co: self.co / self.co.w(),
+            co: Vec3A::new(x, y, z),
         }
     }
 
     #[inline(always)]
     pub fn min(&self, other: Point) -> Point {
-        let n1 = self.norm();
-        let n2 = other.norm();
+        let n1 = self;
+        let n2 = other;
 
         Point {
             co: n1.co.min(n2.co),
@@ -44,8 +35,8 @@ impl Point {
 
     #[inline(always)]
     pub fn max(&self, other: Point) -> Point {
-        let n1 = self.norm();
-        let n2 = other.norm();
+        let n1 = self;
+        let n2 = other;
 
         Point {
             co: n1.co.max(n2.co),
@@ -54,9 +45,7 @@ impl Point {
 
     #[inline(always)]
     pub fn into_vector(self) -> Vector {
-        Vector {
-            co: self.co.truncate(),
-        }
+        Vector { co: self.co }
     }
 
     #[inline(always)]
@@ -71,32 +60,32 @@ impl Point {
 
     #[inline(always)]
     pub fn x(&self) -> f32 {
-        self.co.x()
+        self.co[0]
     }
 
     #[inline(always)]
     pub fn y(&self) -> f32 {
-        self.co.y()
+        self.co[1]
     }
 
     #[inline(always)]
     pub fn z(&self) -> f32 {
-        self.co.z()
+        self.co[2]
     }
 
     #[inline(always)]
     pub fn set_x(&mut self, x: f32) {
-        self.co.set_x(x);
+        self.co[0] = x;
     }
 
     #[inline(always)]
     pub fn set_y(&mut self, y: f32) {
-        self.co.set_y(y);
+        self.co[1] = y;
     }
 
     #[inline(always)]
     pub fn set_z(&mut self, z: f32) {
-        self.co.set_z(z);
+        self.co[2] = z;
     }
 }
 
@@ -113,7 +102,7 @@ impl Add<Vector> for Point {
     #[inline(always)]
     fn add(self, other: Vector) -> Point {
         Point {
-            co: self.co + other.co.extend(0.0),
+            co: self.co + other.co,
         }
     }
 }
@@ -124,7 +113,7 @@ impl Sub for Point {
     #[inline(always)]
     fn sub(self, other: Point) -> Vector {
         Vector {
-            co: (self.norm().co - other.norm().co).truncate(),
+            co: self.co - other.co,
         }
     }
 }
@@ -135,36 +124,27 @@ impl Sub<Vector> for Point {
     #[inline(always)]
     fn sub(self, other: Vector) -> Point {
         Point {
-            co: self.co - other.co.extend(0.0),
+            co: self.co - other.co,
         }
     }
 }
 
-impl Mul<Matrix4x4> for Point {
+impl Mul<Transform> for Point {
     type Output = Point;
 
     #[inline]
-    fn mul(self, other: Matrix4x4) -> Point {
+    fn mul(self, other: Transform) -> Point {
         Point {
-            co: other.0.mul_vec4(self.co),
+            co: other.0.transform_point3a(self.co),
         }
     }
 }
 
 #[cfg(test)]
 mod tests {
-    use super::super::{Matrix4x4, Vector};
+    use super::super::{Transform, Vector};
     use super::*;
 
-    #[test]
-    fn norm() {
-        let mut p1 = Point::new(1.0, 2.0, 3.0);
-        let p2 = Point::new(2.0, 4.0, 6.0);
-        p1.co.set_w(0.5);
-
-        assert_eq!(p2, p1.norm());
-    }
-
     #[test]
     fn add() {
         let p1 = Point::new(1.0, 2.0, 3.0);
@@ -186,8 +166,8 @@ mod tests {
     #[test]
     fn mul_matrix_1() {
         let p = Point::new(1.0, 2.5, 4.0);
-        let m = Matrix4x4::new_from_values(
-            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0, 0.0, 0.0, 0.0, 1.0,
+        let m = Transform::new_from_values(
+            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0,
         );
         let pm = Point::new(15.5, 54.0, 70.0);
         assert_eq!(p * m, pm);
@@ -196,11 +176,10 @@ mod tests {
     #[test]
     fn mul_matrix_2() {
         let p = Point::new(1.0, 2.5, 4.0);
-        let m = Matrix4x4::new_from_values(
-            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0, 2.0, 3.0, 1.0, 5.0,
+        let m = Transform::new_from_values(
+            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0,
         );
-        let mut pm = Point::new(15.5, 54.0, 70.0);
-        pm.co.set_w(18.5);
+        let pm = Point::new(15.5, 54.0, 70.0);
         assert_eq!(p * m, pm);
     }
 
@@ -208,12 +187,11 @@ mod tests {
     fn mul_matrix_3() {
         // Make sure matrix multiplication composes the way one would expect
         let p = Point::new(1.0, 2.5, 4.0);
-        let m1 = Matrix4x4::new_from_values(
-            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0, 13.0, 7.0, 15.0, 3.0,
-        );
-        let m2 = Matrix4x4::new_from_values(
-            4.0, 1.0, 2.0, 3.5, 3.0, 6.0, 5.0, 2.0, 2.0, 2.0, 4.0, 12.0, 5.0, 7.0, 8.0, 11.0,
+        let m1 = Transform::new_from_values(
+            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0,
         );
+        let m2 =
+            Transform::new_from_values(4.0, 1.0, 2.0, 3.5, 3.0, 6.0, 5.0, 2.0, 2.0, 2.0, 4.0, 12.0);
         println!("{:?}", m1 * m2);
 
         let pmm1 = p * (m1 * m2);
diff --git a/sub_crates/math3d/src/transform.rs b/sub_crates/math3d/src/transform.rs
new file mode 100644
index 0000000..29cd069
--- /dev/null
+++ b/sub_crates/math3d/src/transform.rs
@@ -0,0 +1,178 @@
+#![allow(dead_code)]
+
+use std::ops::{Add, Mul};
+
+use approx::relative_eq;
+use glam::{Affine3A, Mat3, Mat4, Vec3};
+
+use super::Point;
+
+/// A 4x3 affine transform matrix, used for transforms.
+#[derive(Debug, Copy, Clone, PartialEq)]
+pub struct Transform(pub Affine3A);
+
+impl Transform {
+    /// Creates a new identity matrix
+    #[inline]
+    pub fn new() -> Transform {
+        Transform(Affine3A::IDENTITY)
+    }
+
+    /// Creates a new matrix with the specified values:
+    /// a b c d
+    /// e f g h
+    /// i j k l
+    /// m n o p
+    #[inline]
+    #[allow(clippy::many_single_char_names)]
+    #[allow(clippy::too_many_arguments)]
+    pub fn new_from_values(
+        a: f32,
+        b: f32,
+        c: f32,
+        d: f32,
+        e: f32,
+        f: f32,
+        g: f32,
+        h: f32,
+        i: f32,
+        j: f32,
+        k: f32,
+        l: f32,
+    ) -> Transform {
+        Transform(Affine3A::from_mat3_translation(
+            Mat3::from_cols(Vec3::new(a, e, i), Vec3::new(b, f, j), Vec3::new(c, g, k)),
+            Vec3::new(d, h, l),
+        ))
+    }
+
+    #[inline]
+    pub fn from_location(loc: Point) -> Transform {
+        Transform(Affine3A::from_translation(loc.co.into()))
+    }
+
+    /// Returns whether the matrices are approximately equal to each other.
+    /// Each corresponding element in the matrices cannot have a relative
+    /// error exceeding epsilon.
+    #[inline]
+    pub fn aprx_eq(&self, other: Transform, epsilon: f32) -> bool {
+        let mut eq = true;
+        for c in 0..3 {
+            for r in 0..3 {
+                let a = self.0.matrix3.col(c)[r];
+                let b = other.0.matrix3.col(c)[r];
+                eq &= relative_eq!(a, b, epsilon = epsilon);
+            }
+        }
+        for i in 0..3 {
+            let a = self.0.translation[i];
+            let b = other.0.translation[i];
+            eq &= relative_eq!(a, b, epsilon = epsilon);
+        }
+        eq
+    }
+
+    /// Returns the inverse of the Matrix
+    #[inline]
+    pub fn inverse(&self) -> Transform {
+        Transform(self.0.inverse())
+    }
+}
+
+impl Default for Transform {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+/// Multiply two matrices together
+impl Mul for Transform {
+    type Output = Self;
+
+    #[inline]
+    fn mul(self, other: Self) -> Self {
+        Self(other.0 * self.0)
+    }
+}
+
+/// Multiply a matrix by a f32
+impl Mul<f32> for Transform {
+    type Output = Self;
+
+    #[inline]
+    fn mul(self, other: f32) -> Self {
+        Self(Affine3A::from_mat4(Mat4::from(self.0) * other))
+    }
+}
+
+/// Add two matrices together
+impl Add for Transform {
+    type Output = Self;
+
+    #[inline]
+    fn add(self, other: Self) -> Self {
+        Self(Affine3A::from_mat4(
+            Mat4::from(self.0) + Mat4::from(other.0),
+        ))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn equality_test() {
+        let a = Transform::new();
+        let b = Transform::new();
+        let c =
+            Transform::new_from_values(1.1, 0.0, 0.0, 0.0, 0.0, 1.1, 0.0, 0.0, 0.0, 0.0, 1.1, 0.0);
+
+        assert_eq!(a, b);
+        assert!(a != c);
+    }
+
+    #[test]
+    fn approximate_equality_test() {
+        let a = Transform::new();
+        let b = Transform::new_from_values(
+            1.000001, 0.0, 0.0, 0.0, 0.0, 1.000001, 0.0, 0.0, 0.0, 0.0, 1.000001, 0.0,
+        );
+        let c = Transform::new_from_values(
+            1.000003, 0.0, 0.0, 0.0, 0.0, 1.000003, 0.0, 0.0, 0.0, 0.0, 1.000003, 0.0,
+        );
+        let d = Transform::new_from_values(
+            -1.000001, 0.0, 0.0, 0.0, 0.0, -1.000001, 0.0, 0.0, 0.0, 0.0, -1.000001, 0.0,
+        );
+
+        assert!(a.aprx_eq(b, 0.000001));
+        assert!(!a.aprx_eq(c, 0.000001));
+        assert!(!a.aprx_eq(d, 0.000001));
+    }
+
+    #[test]
+    fn multiply_test() {
+        let a = Transform::new_from_values(
+            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0,
+        );
+        let b = Transform::new_from_values(
+            1.0, 5.0, 9.0, 13.0, 2.0, 6.0, 10.0, 14.0, 3.0, 7.0, 11.0, 15.0,
+        );
+        let c = Transform::new_from_values(
+            97.0, 50.0, 136.0, 162.5, 110.0, 60.0, 156.0, 185.0, 123.0, 70.0, 176.0, 207.5,
+        );
+
+        assert_eq!(a * b, c);
+    }
+
+    #[test]
+    fn inverse_test() {
+        let a = Transform::new_from_values(
+            1.0, 0.33, 0.0, -2.0, 0.0, 1.0, 0.0, 0.0, 2.1, 0.7, 1.3, 0.0,
+        );
+        let b = a.inverse();
+        let c = Transform::new();
+
+        assert!((dbg!(a * b)).aprx_eq(dbg!(c), 0.0000001));
+    }
+}
diff --git a/sub_crates/math3d/src/vector.rs b/sub_crates/math3d/src/vector.rs
index 3e409f0..365e892 100644
--- a/sub_crates/math3d/src/vector.rs
+++ b/sub_crates/math3d/src/vector.rs
@@ -5,21 +5,21 @@ use std::{
     ops::{Add, Div, Mul, Neg, Sub},
 };
 
-use glam::Vec3;
+use glam::Vec3A;
 
-use super::{CrossProduct, DotProduct, Matrix4x4, Normal, Point};
+use super::{CrossProduct, DotProduct, Normal, Point, Transform};
 
 /// A direction vector in 3d homogeneous space.
 #[derive(Debug, Copy, Clone)]
 pub struct Vector {
-    pub co: Vec3,
+    pub co: Vec3A,
 }
 
 impl Vector {
     #[inline(always)]
     pub fn new(x: f32, y: f32, z: f32) -> Vector {
         Vector {
-            co: Vec3::new(x, y, z),
+            co: Vec3A::new(x, y, z),
         }
     }
 
@@ -43,15 +43,13 @@ impl Vector {
     #[inline(always)]
     pub fn abs(&self) -> Vector {
         Vector {
-            co: self.co * self.co.sign(),
+            co: self.co * self.co.signum(),
         }
     }
 
     #[inline(always)]
     pub fn into_point(self) -> Point {
-        Point {
-            co: self.co.extend(1.0),
-        }
+        Point { co: self.co }
     }
 
     #[inline(always)]
@@ -71,32 +69,32 @@ impl Vector {
 
     #[inline(always)]
     pub fn x(&self) -> f32 {
-        self.co.x()
+        self.co[0]
     }
 
     #[inline(always)]
     pub fn y(&self) -> f32 {
-        self.co.y()
+        self.co[1]
     }
 
     #[inline(always)]
     pub fn z(&self) -> f32 {
-        self.co.z()
+        self.co[2]
     }
 
     #[inline(always)]
     pub fn set_x(&mut self, x: f32) {
-        self.co.set_x(x);
+        self.co[0] = x;
     }
 
     #[inline(always)]
     pub fn set_y(&mut self, y: f32) {
-        self.co.set_y(y);
+        self.co[1] = y;
     }
 
     #[inline(always)]
     pub fn set_z(&mut self, z: f32) {
-        self.co.set_z(z);
+        self.co[2] = z;
     }
 }
 
@@ -140,13 +138,13 @@ impl Mul<f32> for Vector {
     }
 }
 
-impl Mul<Matrix4x4> for Vector {
+impl Mul<Transform> for Vector {
     type Output = Vector;
 
     #[inline]
-    fn mul(self, other: Matrix4x4) -> Vector {
+    fn mul(self, other: Transform) -> Vector {
         Vector {
-            co: other.0.transform_vector3(self.co),
+            co: other.0.transform_vector3a(self.co),
         }
     }
 }
@@ -189,7 +187,7 @@ impl CrossProduct for Vector {
 
 #[cfg(test)]
 mod tests {
-    use super::super::{CrossProduct, DotProduct, Matrix4x4};
+    use super::super::{CrossProduct, DotProduct, Transform};
     use super::*;
 
     #[test]
@@ -222,8 +220,8 @@ mod tests {
     #[test]
     fn mul_matrix_1() {
         let v = Vector::new(1.0, 2.5, 4.0);
-        let m = Matrix4x4::new_from_values(
-            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0, 13.0, 7.0, 15.0, 3.0,
+        let m = Transform::new_from_values(
+            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0,
         );
         assert_eq!(v * m, Vector::new(14.0, 46.0, 58.0));
     }
@@ -231,8 +229,8 @@ mod tests {
     #[test]
     fn mul_matrix_2() {
         let v = Vector::new(1.0, 2.5, 4.0);
-        let m = Matrix4x4::new_from_values(
-            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0, 0.0, 0.0, 0.0, 1.0,
+        let m = Transform::new_from_values(
+            1.0, 2.0, 2.0, 1.5, 3.0, 6.0, 7.0, 8.0, 9.0, 2.0, 11.0, 12.0,
         );
         assert_eq!(v * m, Vector::new(14.0, 46.0, 58.0));
     }
diff --git a/sub_crates/spectral_upsampling/Cargo.toml b/sub_crates/spectral_upsampling/Cargo.toml
index 7d22204..88be212 100644
--- a/sub_crates/spectral_upsampling/Cargo.toml
+++ b/sub_crates/spectral_upsampling/Cargo.toml
@@ -10,4 +10,4 @@ name = "spectral_upsampling"
 path = "src/lib.rs"
 
 [dependencies]
-glam = {git="https://github.com/bitshifter/glam-rs.git", rev="0f314f99", default-features=false, features=["approx"]}
\ No newline at end of file
+glam = "0.15"
\ No newline at end of file
diff --git a/sub_crates/spectral_upsampling/src/jakob.rs b/sub_crates/spectral_upsampling/src/jakob.rs
index b33f619..9a8e07c 100644
--- a/sub_crates/spectral_upsampling/src/jakob.rs
+++ b/sub_crates/spectral_upsampling/src/jakob.rs
@@ -118,14 +118,14 @@ fn small_rgb_to_spectrum_p4(
 
     // Evaluate the spectral function and return the result.
     if max_val <= table_mid_value {
-        rgb2spec_eval_4([c[0].x(), c[0].y(), c[0].z()], lambdas) * (1.0 / table_mid_value) * max_val
+        rgb2spec_eval_4([c[0][0], c[0][1], c[0][2]], lambdas) * (1.0 / table_mid_value) * max_val
     } else if max_val < 1.0 {
         let n = (max_val - table_mid_value) / (1.0 - table_mid_value);
-        let s0 = rgb2spec_eval_4([c[0].x(), c[0].y(), c[0].z()], lambdas);
-        let s1 = rgb2spec_eval_4([c[1].x(), c[1].y(), c[1].z()], lambdas);
+        let s0 = rgb2spec_eval_4([c[0][0], c[0][1], c[0][2]], lambdas);
+        let s1 = rgb2spec_eval_4([c[1][0], c[1][1], c[1][2]], lambdas);
         (s0 * (1.0 - n)) + (s1 * n)
     } else {
-        rgb2spec_eval_4([c[1].x(), c[1].y(), c[1].z()], lambdas) * max_val
+        rgb2spec_eval_4([c[1][0], c[1][1], c[1][2]], lambdas) * max_val
     }
 }
 
@@ -147,7 +147,7 @@ fn rgb2spec_eval_4(coeff: [f32; RGB2SPEC_N_COEFFS], lambda: Vec4) -> Vec4 {
     let y = {
         // TODO: replace this with a SIMD sqrt op.
         let (x, y, z, w) = rgb2spec_fma_4(x, x, Vec4::splat(1.0)).into();
-        Vec4::new(x.sqrt(), y.sqrt(), z.sqrt(), w.sqrt()).reciprocal()
+        Vec4::new(x.sqrt(), y.sqrt(), z.sqrt(), w.sqrt()).recip()
     };
 
     rgb2spec_fma_4(Vec4::splat(0.5) * x, y, Vec4::splat(0.5))
diff --git a/sub_crates/spectral_upsampling/src/meng.rs b/sub_crates/spectral_upsampling/src/meng.rs
index 953940e..f3ad3e7 100644
--- a/sub_crates/spectral_upsampling/src/meng.rs
+++ b/sub_crates/spectral_upsampling/src/meng.rs
@@ -227,13 +227,13 @@ pub fn spectrum_xyz_to_p_4(lambdas: Vec4, xyz: (f32, f32, f32)) -> Vec4 {
     // Get the spectral values for the vertices of the grid cell.
     // TODO: use integer SIMD intrinsics to make this part faster.
     let mut p = [Vec4::splat(0.0); 6];
-    let sb0: [i32; 4] = [sb.x() as i32, sb.y() as i32, sb.z() as i32, sb.w() as i32];
+    let sb0: [i32; 4] = [sb[0] as i32, sb[1] as i32, sb[2] as i32, sb[3] as i32];
     assert!(sb0[0].max(sb0[1]).max(sb0[2].max(sb0[3])) < SPECTRUM_NUM_SAMPLES);
     let sb1: [i32; 4] = [
-        (sb.x() as i32 + 1).min(SPECTRUM_NUM_SAMPLES - 1),
-        (sb.y() as i32 + 1).min(SPECTRUM_NUM_SAMPLES - 1),
-        (sb.z() as i32 + 1).min(SPECTRUM_NUM_SAMPLES - 1),
-        (sb.w() as i32 + 1).min(SPECTRUM_NUM_SAMPLES - 1),
+        (sb[0] as i32 + 1).min(SPECTRUM_NUM_SAMPLES - 1),
+        (sb[1] as i32 + 1).min(SPECTRUM_NUM_SAMPLES - 1),
+        (sb[2] as i32 + 1).min(SPECTRUM_NUM_SAMPLES - 1),
+        (sb[3] as i32 + 1).min(SPECTRUM_NUM_SAMPLES - 1),
     ];
     let sbf = sb - Vec4::new(sb0[0] as f32, sb0[1] as f32, sb0[2] as f32, sb0[3] as f32);
     for i in 0..(num as usize) {