wathiede/raytracers
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Copy Matrix4x4 impl from pbrt and start on tests

Browse Source
This commit is contained in:
Bill Thiede 2021-06-29 20:12:23 -07:00
parent 3cf580f607
commit 72b15e5516
2 changed files with 262 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
rtchallenge/src/lib.rs
View File

@ -1,2 +1,3 @@
pub mod canvas; pub mod canvas;
pub mod matrices;
pub mod tuples; pub mod tuples;

261
rtchallenge/src/matrices.rs Normal file
View File

@ -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 m<sub>i,j</sub>
/// 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<sup>-1</sup>, where
/// AA<sup>-1</sup> = 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<Matrix4x4> 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);
}
}