From 72b15e5516a4e5832f36a5f5d9403734e5e48872 Mon Sep 17 00:00:00 2001 From: Bill Thiede Date: Tue, 29 Jun 2021 20:12:23 -0700 Subject: [PATCH] Copy Matrix4x4 impl from pbrt and start on tests --- rtchallenge/src/lib.rs | 1 + rtchallenge/src/matrices.rs | 261 ++++++++++++++++++++++++++++++++++++ 2 files changed, 262 insertions(+) create mode 100644 rtchallenge/src/matrices.rs diff --git a/rtchallenge/src/lib.rs b/rtchallenge/src/lib.rs index afabcff..58b9c78 100644 --- a/rtchallenge/src/lib.rs +++ b/rtchallenge/src/lib.rs @@ -1,2 +1,3 @@ pub mod canvas; +pub mod matrices; pub mod tuples; diff --git a/rtchallenge/src/matrices.rs b/rtchallenge/src/matrices.rs new file mode 100644 index 0000000..3d3c35f --- /dev/null +++ b/rtchallenge/src/matrices.rs @@ -0,0 +1,261 @@ +use std::fmt; +use std::ops::{Add, Div, Index, Mul, Neg, Sub}; + +#[derive(Default, Clone, Copy)] +/// Matrix4x4 represents a 4x4 matrix in row-major form. So, element `m[i][j]` corresponds to mi,j +/// where `i` is the row number and `j` is the column number. +pub struct Matrix4x4 { + m: [[f32; 4]; 4], +} + +impl From<[f32; 16]> for Matrix4x4 { + fn from(t: [f32; 16]) -> Self { + Matrix4x4 { + m: [ + [t[0], t[1], t[2], t[3]], + [t[4], t[5], t[6], t[7]], + [t[8], t[9], t[10], t[11]], + [t[12], t[13], t[14], t[15]], + ], + } + } +} + +impl Matrix4x4 { + /// Create a `Matrix4x4` containing the identity, all zeros with ones along the diagonal. + pub fn identity() -> Matrix4x4 { + Matrix4x4::new( + [1., 0., 0., 0.], + [0., 1., 0., 0.], + [0., 0., 1., 0.], + [0., 0., 0., 1.], + ) + } + + /// Create a `Matrix4x4` with each of the given rows. + pub fn new(r0: [f32; 4], r1: [f32; 4], r2: [f32; 4], r3: [f32; 4]) -> Matrix4x4 { + Matrix4x4 { + m: [r0, r1, r2, r3], + } + } + + /// Transpose self, returning a new matrix that has been reflected across the diagonal. + /// # Examples + /// + /// ``` + /// use rtchallenge::matrices::Matrix4x4; + /// + /// let m = Matrix4x4::new( + /// [2., 0., 0., 0.], + /// [3., 1., 0., 0.], + /// [4., 0., 1., 0.], + /// [5., 6., 7., 1.], + /// ); + /// let m_t = Matrix4x4::new( + /// [2., 3., 4., 5.], + /// [0., 1., 0., 6.], + /// [0., 0., 1., 7.], + /// [0., 0., 0., 1.], + /// ); + /// assert_eq!(m.transpose(), m_t); + pub fn transpose(&self) -> Matrix4x4 { + let m = self.m; + Matrix4x4 { + m: [ + [m[0][0], m[1][0], m[2][0], m[3][0]], + [m[0][1], m[1][1], m[2][1], m[3][1]], + [m[0][2], m[1][2], m[2][2], m[3][2]], + [m[0][3], m[1][3], m[2][3], m[3][3]], + ], + } + } + + /// Returns a new matrix that is the inverse of self. If self is A, inverse returns A-1, where + /// AA-1 = I. + /// This implementation uses a numerically stable Gauss–Jordan elimination routine to compute the inverse. + /// + /// # Examples + /// + /// ``` + /// use rtchallenge::matrices::Matrix4x4; + /// + /// let i = Matrix4x4::identity(); + /// assert_eq!(i.inverse() * i, i); + /// + /// let m = Matrix4x4::new( + /// [2., 0., 0., 0.], + /// [0., 3., 0., 0.], + /// [0., 0., 4., 0.], + /// [0., 0., 0., 1.], + /// ); + /// assert_eq!(m.inverse() * m, i); + /// assert_eq!(m * m.inverse(), i); + /// ``` + pub fn inverse(&self) -> Matrix4x4 { + // TODO(wathiede): how come the C++ version doesn't need to deal with non-invertable + // matrix. + let mut indxc: [usize; 4] = Default::default(); + let mut indxr: [usize; 4] = Default::default(); + let mut ipiv: [usize; 4] = Default::default(); + let mut minv = self.m; + + for i in 0..4 { + let mut irow: usize = 0; + let mut icol: usize = 0; + let mut big: f32 = 0.; + // Choose pivot + for j in 0..4 { + if ipiv[j] != 1 { + for (k, ipivk) in ipiv.iter().enumerate() { + if *ipivk == 0 { + if minv[j][k].abs() >= big { + big = minv[j][k].abs(); + irow = j; + icol = k; + } + } else if *ipivk > 1 { + eprintln!("Singular matrix in MatrixInvert"); + } + } + } + } + ipiv[icol] += 1; + // Swap rows _irow_ and _icol_ for pivot + if irow != icol { + // Can't figure out how to make swap work here. + #[allow(clippy::manual_swap)] + for k in 0..4 { + let tmp = minv[irow][k]; + minv[irow][k] = minv[icol][k]; + minv[icol][k] = tmp; + } + } + indxr[i] = irow; + indxc[i] = icol; + if minv[icol][icol] == 0. { + eprintln!("Singular matrix in MatrixInvert"); + } + + // Set $m[icol][icol]$ to one by scaling row _icol_ appropriately + let pivinv: f32 = minv[icol][icol].recip(); + minv[icol][icol] = 1.; + for j in 0..4 { + minv[icol][j] *= pivinv; + } + + // Subtract this row from others to zero out their columns + for j in 0..4 { + if j != icol { + let save = minv[j][icol]; + minv[j][icol] = 0.; + for k in 0..4 { + minv[j][k] -= minv[icol][k] * save; + } + } + } + } + // Swap columns to reflect permutation + for j in (0..4).rev() { + if indxr[j] != indxc[j] { + for mi in &mut minv { + mi.swap(indxr[j], indxc[j]) + } + } + } + Matrix4x4 { m: minv } + } +} + +impl fmt::Debug for Matrix4x4 { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + if f.alternate() { + write!( + f, + "{:?}\n {:?}\n {:?}\n {:?}", + self.m[0], self.m[1], self.m[2], self.m[3] + ) + } else { + write!( + f, + "[{:?} {:?} {:?} {:?}]", + self.m[0], self.m[1], self.m[2], self.m[3] + ) + } + } +} + +impl Mul for Matrix4x4 { + type Output = Matrix4x4; + + /// Implement matrix multiplication for `Matrix4x4`. + /// + /// # Examples + /// ``` + /// use rtchallenge::matrices::Matrix4x4; + /// + /// let i = Matrix4x4::identity(); + /// let m1 = Matrix4x4::identity(); + /// let m2 = Matrix4x4::identity(); + /// + /// assert_eq!(m1 * m2, i); + /// ``` + fn mul(self, m2: Matrix4x4) -> Matrix4x4 { + let m1 = self; + let mut r: Matrix4x4 = Default::default(); + for i in 0..4 { + for j in 0..4 { + r.m[i][j] = m1.m[i][0] * m2.m[0][j] + + m1.m[i][1] * m2.m[1][j] + + m1.m[i][2] * m2.m[2][j] + + m1.m[i][3] * m2.m[3][j]; + } + } + r + } +} + +impl PartialEq for Matrix4x4 { + fn eq(&self, rhs: &Matrix4x4) -> bool { + let l = self.m; + let r = rhs.m; + for i in 0..4 { + for j in 0..4 { + let d = (l[i][j] - r[i][j]).abs(); + if d > f32::EPSILON { + return false; + } + } + } + true + } +} + +impl Index<(usize, usize)> for Matrix4x4 { + type Output = f32; + fn index(&self, (row, col): (usize, usize)) -> &Self::Output { + &self.m[row][col] + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn construct4x4() { + let m = Matrix4x4::new( + [1., 2., 3., 4.], + [5.5, 6.5, 7.5, 8.5], + [9., 10., 11., 12.], + [13.5, 14.5, 15.5, 16.5], + ); + + assert_eq!(m[(0, 0)], 1.); + assert_eq!(m[(0, 3)], 4.); + assert_eq!(m[(1, 0)], 5.5); + assert_eq!(m[(1, 2)], 7.5); + assert_eq!(m[(2, 2)], 11.); + assert_eq!(m[(3, 0)], 13.5); + assert_eq!(m[(3, 2)], 15.5); + } +}