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

matrices: move tests from doctest to unit.

Browse Source
This commit is contained in:
Bill Thiede 2021-07-29 20:33:03 -07:00
parent cd2a4770ca
commit 3aea76b35c
1 changed files with 308 additions and 422 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

730
rtchallenge/src/matrices.rs
View File

@ -4,82 +4,30 @@ use std::ops::{Index, IndexMut, Mul, Sub};
use crate::{tuples::Tuple, Float, EPSILON}; use crate::{tuples::Tuple, Float, EPSILON};
/// Short hand for creating a Matrix4x4 set to the identity matrix. /// Short hand for creating a Matrix4x4 set to the identity matrix.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::{identity, Matrix4x4};
///
/// assert_eq!(identity(), Matrix4x4::identity());
/// ```
pub fn identity() -> Matrix4x4 { pub fn identity() -> Matrix4x4 {
Matrix4x4::identity() Matrix4x4::identity()
} }
/// Short hand for creating a Matrix4x4 for rotating around the X-axis. /// Short hand for creating a Matrix4x4 for rotating around the X-axis.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::{rotation_x, Matrix4x4};
///
/// assert_eq!(rotation_x(10.), Matrix4x4::rotation_x(10.));
/// ```
pub fn rotation_x(radians: Float) -> Matrix4x4 { pub fn rotation_x(radians: Float) -> Matrix4x4 {
Matrix4x4::rotation_x(radians) Matrix4x4::rotation_x(radians)
} }
/// Short hand for creating a Matrix4x4 for rotating around the Y-axis. /// Short hand for creating a Matrix4x4 for rotating around the Y-axis.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::{rotation_y, Matrix4x4};
///
/// assert_eq!(rotation_y(10.), Matrix4x4::rotation_y(10.));
/// ```
pub fn rotation_y(radians: Float) -> Matrix4x4 { pub fn rotation_y(radians: Float) -> Matrix4x4 {
Matrix4x4::rotation_y(radians) Matrix4x4::rotation_y(radians)
} }
/// Short hand for creating a Matrix4x4 for rotating around the Z-axis. /// Short hand for creating a Matrix4x4 for rotating around the Z-axis.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::{rotation_z, Matrix4x4};
///
/// assert_eq!(rotation_z(10.), Matrix4x4::rotation_z(10.));
/// ```
pub fn rotation_z(radians: Float) -> Matrix4x4 { pub fn rotation_z(radians: Float) -> Matrix4x4 {
Matrix4x4::rotation_z(radians) Matrix4x4::rotation_z(radians)
} }
/// Short hand for creating a Matrix4x4 that scales in the given x,y,z axis. /// Short hand for creating a Matrix4x4 that scales in the given x,y,z axis.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::{scaling, Matrix4x4};
///
/// assert_eq!(scaling(1., 2., 3.), Matrix4x4::scaling(1., 2., 3.));
/// ```
pub fn scaling(x: Float, y: Float, z: Float) -> Matrix4x4 { pub fn scaling(x: Float, y: Float, z: Float) -> Matrix4x4 {
Matrix4x4::scaling(x, y, z) Matrix4x4::scaling(x, y, z)
} }
/// Short hand for creating a Matrix4x4 that shears across the given axis pairs. /// Short hand for creating a Matrix4x4 that shears across the given axis pairs.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::{shearing, Matrix4x4};
///
/// assert_eq!(
/// shearing(1., 2., 3., 4., 5., 6.),
/// Matrix4x4::shearing(1., 2., 3., 4., 5., 6.)
/// );
/// ```
pub fn shearing(xy: Float, xz: Float, yx: Float, yz: Float, zx: Float, zy: Float) -> Matrix4x4 { pub fn shearing(xy: Float, xz: Float, yx: Float, yz: Float, zx: Float, zy: Float) -> Matrix4x4 {
Matrix4x4::shearing(xy, xz, yx, yz, zx, zy) Matrix4x4::shearing(xy, xz, yx, yz, zx, zy)
} }
/// Short hand for creating a Matrix4x4 that translations along the given x,y,z axis. /// Short hand for creating a Matrix4x4 that translations along the given x,y,z axis.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::{translation, Matrix4x4};
///
/// assert_eq!(translation(1., 2., 3.), Matrix4x4::translation(1., 2., 3.));
/// ```
pub fn translation(x: Float, y: Float, z: Float) -> Matrix4x4 { pub fn translation(x: Float, y: Float, z: Float) -> Matrix4x4 {
Matrix4x4::translation(x, y, z) Matrix4x4::translation(x, y, z)
} }
@ -95,16 +43,6 @@ impl Matrix2x2 {
} }
/// Calculate the determinant of a 2x2. /// Calculate the determinant of a 2x2.
///
/// # Examples
///
/// ```
/// use rtchallenge::matrices::Matrix2x2;
///
/// let a = Matrix2x2::new([1., 5.], [-3., 2.]);
///
/// assert_eq!(a.determinant(), 17.);
/// ```
pub fn determinant(&self) -> Float { pub fn determinant(&self) -> Float {
let m = self; let m = self;
m[(0, 0)] * m[(1, 1)] - m[(0, 1)] * m[(1, 0)] m[(0, 0)] * m[(1, 1)] - m[(0, 1)] * m[(1, 0)]
@ -142,16 +80,6 @@ impl Matrix3x3 {
Matrix3x3 { m: [r0, r1, r2] } Matrix3x3 { m: [r0, r1, r2] }
} }
/// submatrix extracts a 2x2 matrix ignoring the 0-based `row` and `col` given. /// submatrix extracts a 2x2 matrix ignoring the 0-based `row` and `col` given.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::{Matrix2x2, Matrix3x3};
///
/// assert_eq!(
/// Matrix3x3::new([1., 5., 0.], [-3., 2., 7.], [0., 6., -3.],).submatrix(0, 2),
/// Matrix2x2::new([-3., 2.], [0., 6.])
/// );
/// ```
pub fn submatrix(&self, row: usize, col: usize) -> Matrix2x2 { pub fn submatrix(&self, row: usize, col: usize) -> Matrix2x2 {
assert!(row < 3); assert!(row < 3);
assert!(col < 3); assert!(col < 3);
@ -172,49 +100,17 @@ impl Matrix3x3 {
} }
/// Compute minor of a 3x3 matrix. /// Compute minor of a 3x3 matrix.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::Matrix3x3;
///
/// let a = Matrix3x3::new([3., 5., 0.], [2., -1., -7.], [6., -1., 5.]);
/// let b = a.submatrix(1, 0);
/// assert_eq!(b.determinant(), 25.0);
/// assert_eq!(b.determinant(), a.minor(1, 0));
/// ```
pub fn minor(&self, row: usize, col: usize) -> Float { pub fn minor(&self, row: usize, col: usize) -> Float {
self.submatrix(row, col).determinant() self.submatrix(row, col).determinant()
} }
/// Compute cofactor of a 3x3 matrix. /// Compute cofactor of a 3x3 matrix.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::Matrix3x3;
///
/// let a = Matrix3x3::new([3., 5., 0.], [2., -1., -7.], [6., -1., 5.]);
/// assert_eq!(a.minor(0, 0), -12.);
/// assert_eq!(a.cofactor(0, 0), -12.);
/// assert_eq!(a.minor(1, 0), 25.);
/// assert_eq!(a.cofactor(1, 0), -25.);
/// ```
pub fn cofactor(&self, row: usize, col: usize) -> Float { pub fn cofactor(&self, row: usize, col: usize) -> Float {
let negate = if (row + col) % 2 == 0 { 1. } else { -1. }; let negate = if (row + col) % 2 == 0 { 1. } else { -1. };
self.submatrix(row, col).determinant() * negate self.submatrix(row, col).determinant() * negate
} }
/// Compute determinant of a 3x3 matrix. /// Compute determinant of a 3x3 matrix.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::Matrix3x3;
///
/// let a = Matrix3x3::new([1., 2., 6.], [-5., 8., -4.], [2., 6., 4.]);
/// assert_eq!(a.cofactor(0, 0), 56.);
/// assert_eq!(a.cofactor(0, 1), 12.);
/// assert_eq!(a.cofactor(0, 2), -46.);
/// assert_eq!(a.determinant(), -196.);
/// ```
pub fn determinant(&self) -> Float { pub fn determinant(&self) -> Float {
(0..3).map(|i| self.cofactor(0, i) * self[(0, i)]).sum() (0..3).map(|i| self.cofactor(0, i) * self[(0, i)]).sum()
} }
@ -247,30 +143,34 @@ impl PartialEq for Matrix3x3 {
/// ///
/// # Examples /// # Examples
/// ``` /// ```
/// use rtchallenge::{float::consts::PI, matrices::Matrix4x4, tuples::Tuple}; /// use rtchallenge::{
/// float::consts::PI,
/// matrices::Matrix4x4,
/// tuples::{point, vector},
/// };
/// ///
/// // Individual transformations are applied in sequence. /// // Individual transformations are applied in sequence.
/// let p = Tuple::point(1., 0., 1.); /// let p = point(1., 0., 1.);
/// let a = Matrix4x4::rotation_x(PI / 2.); /// let a = Matrix4x4::rotation_x(PI / 2.);
/// let b = Matrix4x4::scaling(5., 5., 5.); /// let b = Matrix4x4::scaling(5., 5., 5.);
/// let c = Matrix4x4::translation(10., 5., 7.); /// let c = Matrix4x4::translation(10., 5., 7.);
/// // Apply rotation first. /// // Apply rotation first.
/// let p2 = a * p; /// let p2 = a * p;
/// assert_eq!(p2, Tuple::point(1., -1., 0.)); /// assert_eq!(p2, point(1., -1., 0.));
/// // Then apply scaling. /// // Then apply scaling.
/// let p3 = b * p2; /// let p3 = b * p2;
/// assert_eq!(p3, Tuple::point(5., -5., 0.)); /// assert_eq!(p3, point(5., -5., 0.));
/// // Then apply translation. /// // Then apply translation.
/// let p4 = c * p3; /// let p4 = c * p3;
/// assert_eq!(p4, Tuple::point(15., 0., 7.)); /// assert_eq!(p4, point(15., 0., 7.));
/// ///
/// // Chained transformations must be applied in reverse order. /// // Chained transformations must be applied in reverse order.
/// let p = Tuple::point(1., 0., 1.); /// let p = point(1., 0., 1.);
/// let a = Matrix4x4::rotation_x(PI / 2.); /// let a = Matrix4x4::rotation_x(PI / 2.);
/// let b = Matrix4x4::scaling(5., 5., 5.); /// let b = Matrix4x4::scaling(5., 5., 5.);
/// let c = Matrix4x4::translation(10., 5., 7.); /// let c = Matrix4x4::translation(10., 5., 7.);
/// let t = c * b * a; /// let t = c * b * a;
/// assert_eq!(t * p, Tuple::point(15., 0., 7.)); /// assert_eq!(t * p, point(15., 0., 7.));
/// ``` /// ```
#[derive(Copy, Clone, Default)] #[derive(Copy, Clone, Default)]
pub struct Matrix4x4 { pub struct Matrix4x4 {
@ -292,21 +192,6 @@ impl From<[Float; 16]> for Matrix4x4 {
impl Matrix4x4 { impl Matrix4x4 {
/// Create a `Matrix4x4` containing the identity, all zeros with ones along the diagonal. /// Create a `Matrix4x4` containing the identity, all zeros with ones along the diagonal.
/// # Examples
///
/// ```
/// use rtchallenge::matrices::Matrix4x4;
///
/// let a = Matrix4x4::new(
/// [0., 1., 2., 3.],
/// [1., 2., 4., 8.],
/// [2., 4., 8., 16.],
/// [4., 8., 16., 32.],
/// );
/// let i = Matrix4x4::identity();
///
/// assert_eq!(a * i, a);
/// ```
pub fn identity() -> Matrix4x4 { pub fn identity() -> Matrix4x4 {
Matrix4x4::new( Matrix4x4::new(
[1., 0., 0., 0.], [1., 0., 0., 0.],
@ -324,22 +209,6 @@ impl Matrix4x4 {
} }
/// Creates a 4x4 matrix representing a translation of x,y,z. /// Creates a 4x4 matrix representing a translation of x,y,z.
///
/// # Examples
///
/// ```
/// use rtchallenge::{matrices::Matrix4x4, tuples::Tuple};
///
/// let transform = Matrix4x4::translation(5., -3., 2.);
/// let p = Tuple::point(-3., 4., 5.);
/// assert_eq!(transform * p, Tuple::point(2., 1., 7.));
///
/// let inv = transform.inverse();
/// assert_eq!(inv * p, Tuple::point(-8., 7., 3.));
///
/// let v = Tuple::vector(-3., 4., 5.);
/// assert_eq!(transform * v, v);
/// ```
pub fn translation(x: Float, y: Float, z: Float) -> Matrix4x4 { pub fn translation(x: Float, y: Float, z: Float) -> Matrix4x4 {
Matrix4x4::new( Matrix4x4::new(
[1., 0., 0., x], [1., 0., 0., x],
@ -350,30 +219,6 @@ impl Matrix4x4 {
} }
/// Creates a 4x4 matrix representing a scaling of x,y,z. /// Creates a 4x4 matrix representing a scaling of x,y,z.
///
/// # Examples
///
/// ```
/// use rtchallenge::{matrices::Matrix4x4, tuples::Tuple};
///
/// // A scaling matrix applied to a point.
/// let transform = Matrix4x4::scaling(2., 3., 4.);
/// let p = Tuple::point(-4., 6., 8.);
/// assert_eq!(transform * p, Tuple::point(-8., 18., 32.));
///
/// // A scaling matrix applied to a vector.
/// let v = Tuple::vector(-4., 6., 8.);
/// assert_eq!(transform * v, Tuple::vector(-8., 18., 32.));
///
/// // Multiplying by the inverse of a scaling matrix.
/// let inv = transform.inverse();
/// assert_eq!(inv * v, Tuple::vector(-2., 2., 2.));
///
/// // Reflection is scaling by a negative value.
/// let transform = Matrix4x4::scaling(-1., 1., 1.);
/// let p = Tuple::point(2., 3., 4.);
/// assert_eq!(transform * p, Tuple::point(-2., 3., 4.));
/// ```
pub fn scaling(x: Float, y: Float, z: Float) -> Matrix4x4 { pub fn scaling(x: Float, y: Float, z: Float) -> Matrix4x4 {
Matrix4x4::new( Matrix4x4::new(
[x, 0., 0., 0.], [x, 0., 0., 0.],
@ -384,23 +229,6 @@ impl Matrix4x4 {
} }
/// Creates a 4x4 matrix representing a rotation around the x-axis. /// Creates a 4x4 matrix representing a rotation around the x-axis.
///
/// # Examples
///
/// ```
/// use rtchallenge::{float::consts::PI, matrices::Matrix4x4, tuples::Tuple, Float};
///
/// // A scaling matrix applied to a point.
/// let p = Tuple::point(0., 1., 0.);
/// let half_quarter = Matrix4x4::rotation_x(PI / 4.);
/// let full_quarter = Matrix4x4::rotation_x(PI / 2.);
///
/// assert_eq!(
/// half_quarter * p,
/// Tuple::point(0., (2.0 as Float).sqrt() / 2., (2.0 as Float).sqrt() / 2.)
/// );
/// assert_eq!(full_quarter * p, Tuple::point(0., 0., 1.),);
/// ```
pub fn rotation_x(radians: Float) -> Matrix4x4 { pub fn rotation_x(radians: Float) -> Matrix4x4 {
let r = radians; let r = radians;
Matrix4x4::new( Matrix4x4::new(
@ -412,23 +240,6 @@ impl Matrix4x4 {
} }
/// Creates a 4x4 matrix representing a rotation around the y-axis. /// Creates a 4x4 matrix representing a rotation around the y-axis.
///
/// # Examples
///
/// ```
/// use rtchallenge::{float::consts::PI, matrices::Matrix4x4, tuples::Tuple, Float};
///
/// // A scaling matrix applied to a point.
/// let p = Tuple::point(0., 0., 1.);
/// let half_quarter = Matrix4x4::rotation_y(PI / 4.);
/// let full_quarter = Matrix4x4::rotation_y(PI / 2.);
///
/// assert_eq!(
/// half_quarter * p,
/// Tuple::point((2.0 as Float).sqrt() / 2., 0., (2.0 as Float).sqrt() / 2.)
/// );
/// assert_eq!(full_quarter * p, Tuple::point(1., 0., 0.,),);
/// ```
pub fn rotation_y(radians: Float) -> Matrix4x4 { pub fn rotation_y(radians: Float) -> Matrix4x4 {
let r = radians; let r = radians;
Matrix4x4::new( Matrix4x4::new(
@ -440,23 +251,6 @@ impl Matrix4x4 {
} }
/// Creates a 4x4 matrix representing a rotation around the z-axis. /// Creates a 4x4 matrix representing a rotation around the z-axis.
///
/// # Examples
///
/// ```
/// use rtchallenge::{float::consts::PI, matrices::Matrix4x4, tuples::Tuple, Float};
///
/// // A scaling matrix applied to a point.
/// let p = Tuple::point(0., 1., 0.);
/// let half_quarter = Matrix4x4::rotation_z(PI / 4.);
/// let full_quarter = Matrix4x4::rotation_z(PI / 2.);
///
/// assert_eq!(
/// half_quarter * p,
/// Tuple::point(-(2.0 as Float).sqrt() / 2., (2.0 as Float).sqrt() / 2., 0.)
/// );
/// assert_eq!(full_quarter * p, Tuple::point(-1., 0., 0.,),);
/// ```
pub fn rotation_z(radians: Float) -> Matrix4x4 { pub fn rotation_z(radians: Float) -> Matrix4x4 {
let r = radians; let r = radians;
Matrix4x4::new( Matrix4x4::new(
@ -468,26 +262,6 @@ impl Matrix4x4 {
} }
/// Transpose self, returning a new matrix that has been reflected across the diagonal. /// 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);
///
/// assert_eq!(Matrix4x4::identity(), Matrix4x4::identity().transpose());
pub fn transpose(&self) -> Matrix4x4 { pub fn transpose(&self) -> Matrix4x4 {
let m = self.m; let m = self.m;
Matrix4x4 { Matrix4x4 {
@ -501,40 +275,6 @@ impl Matrix4x4 {
} }
/// Create a transform matrix that will shear (skew) points. /// Create a transform matrix that will shear (skew) points.
/// # Examples /// # Examples
///
/// ```
/// use rtchallenge::{matrices::Matrix4x4, tuples::Tuple};
///
/// // A shearing transform moves x in proportion to y.
/// let transform = Matrix4x4::shearing(1.,0.,0.,0.,0.,0.);
/// let p = Tuple::point(2.,3.,4.);
/// assert_eq!(transform * p, Tuple::point(5.,3.,4.));
///
/// // A shearing transform moves x in proportion to z.
/// let transform = Matrix4x4::shearing(0.,1.,0.,0.,0.,0.);
/// let p = Tuple::point(2.,3.,4.);
/// assert_eq!(transform * p, Tuple::point(6.,3.,4.));
///
/// // A shearing transform moves y in proportion to x.
/// let transform = Matrix4x4::shearing(0.,0.,1.,0.,0.,0.);
/// let p = Tuple::point(2.,3.,4.);
/// assert_eq!(transform * p, Tuple::point(2.,5.,4.));
///
/// // A shearing transform moves y in proportion to z.
/// let transform = Matrix4x4::shearing(0.,0.,0.,1.,0.,0.);
/// let p = Tuple::point(2.,3.,4.);
/// assert_eq!(transform * p, Tuple::point(2.,7.,4.));
///
/// // A shearing transform moves z in proportion to x.
/// let transform = Matrix4x4::shearing(0.,0.,0.,0.,1.,0.);
/// let p = Tuple::point(2.,3.,4.);
/// assert_eq!(transform * p, Tuple::point(2.,3.,6.));
///
/// // A shearing transform moves z in proportion to y.
/// let transform = Matrix4x4::shearing(0.,0.,0.,0.,0.,1.);
/// let p = Tuple::point(2.,3.,4.);
/// assert_eq!(transform * p, Tuple::point(2.,3.,7.));
pub fn shearing(xy: Float, xz: Float, yx: Float, yz: Float, zx: Float, zy: Float) -> Matrix4x4 { pub fn shearing(xy: Float, xz: Float, yx: Float, yz: Float, zx: Float, zy: Float) -> Matrix4x4 {
Matrix4x4::new( Matrix4x4::new(
[1., xy, xz, 0.], [1., xy, xz, 0.],
@ -547,24 +287,6 @@ impl Matrix4x4 {
/// Returns a new matrix that is the inverse of self. If self is A, inverse returns A<sup>-1</sup>, where /// 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. /// AA<sup>-1</sup> = I.
/// This implementation uses a numerically stable Gauss–Jordan elimination routine to compute the inverse. /// 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_rtiow() * 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_rtiow() * m, i);
/// assert_eq!(m * m.inverse_rtiow(), i);
/// ```
pub fn inverse_rtiow(&self) -> Matrix4x4 { pub fn inverse_rtiow(&self) -> Matrix4x4 {
// TODO(wathiede): how come the C++ version doesn't need to deal with non-invertable // TODO(wathiede): how come the C++ version doesn't need to deal with non-invertable
// matrix. // matrix.
@ -639,22 +361,6 @@ impl Matrix4x4 {
Matrix4x4 { m: minv } Matrix4x4 { m: minv }
} }
/// submatrix extracts a 3x3 matrix ignoring the 0-based `row` and `col` given. /// submatrix extracts a 3x3 matrix ignoring the 0-based `row` and `col` given.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::{Matrix3x3, Matrix4x4};
///
/// assert_eq!(
/// Matrix4x4::new(
/// [-6., 1., 1., 6.],
/// [-8., 5., 8., 6.],
/// [-1., 0., 8., 2.],
/// [-7., 1., -1., 1.],
/// )
/// .submatrix(2, 1),
/// Matrix3x3::new([-6., 1., 6.], [-8., 8., 6.], [-7., -1., 1.],)
/// );
/// ```
pub fn submatrix(&self, row: usize, col: usize) -> Matrix3x3 { pub fn submatrix(&self, row: usize, col: usize) -> Matrix3x3 {
assert!(row < 4); assert!(row < 4);
assert!(col < 4); assert!(col < 4);
@ -688,133 +394,16 @@ impl Matrix4x4 {
self.submatrix(row, col).determinant() * negate self.submatrix(row, col).determinant() * negate
} }
/// Compute determinant of a 4x4 matrix. /// Compute determinant of a 4x4 matrix.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::Matrix4x4;
///
/// let a = Matrix4x4::new(
/// [-2., -8., 3., 5.],
/// [-3., 1., 7., 3.],
/// [1., 2., -9., 6.],
/// [-6., 7., 7., -9.],
/// );
/// assert_eq!(a.cofactor(0, 0), 690.);
/// assert_eq!(a.cofactor(0, 1), 447.);
/// assert_eq!(a.cofactor(0, 2), 210.);
/// assert_eq!(a.cofactor(0, 3), 51.);
/// assert_eq!(a.determinant(), -4071.);
/// ```
pub fn determinant(&self) -> Float { pub fn determinant(&self) -> Float {
(0..4).map(|i| self.cofactor(0, i) * self[(0, i)]).sum() (0..4).map(|i| self.cofactor(0, i) * self[(0, i)]).sum()
} }
/// Compute invertibility of matrix (i.e. non-zero determinant. /// Compute invertibility of matrix (i.e. non-zero determinant.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::Matrix4x4;
///
/// let a = Matrix4x4::new(
/// [6., 4., 4., 4.],
/// [5., 5., 7., 6.],
/// [4., -9., 3., -7.],
/// [9., 1., 7., -6.],
/// );
/// assert_eq!(a.determinant(), -2120.);
/// assert_eq!(a.invertable(), true);
///
/// let a = Matrix4x4::new(
/// [-4., 2., -2., -3.],
/// [9., 6., 2., 6.],
/// [0., -5., 1., -5.],
/// [0., 0., 0., 0.],
/// );
/// assert_eq!(a.determinant(), 0.);
/// assert_eq!(a.invertable(), false);
/// ```
pub fn invertable(&self) -> bool { pub fn invertable(&self) -> bool {
self.determinant() != 0. self.determinant() != 0.
} }
/// Compute the inverse of a 4x4 matrix. /// Compute the inverse of a 4x4 matrix.
///
/// # Examples
/// ```
/// use rtchallenge::matrices::Matrix4x4;
///
/// let a = Matrix4x4::new(
/// [-5., 2., 6., -8.],
/// [1., -5., 1., 8.],
/// [7., 7., -6., -7.],
/// [1., -3., 7., 4.],
/// );
/// let b = a.inverse();
///
/// assert_eq!(a.determinant(), 532.);
/// assert_eq!(a.cofactor(2, 3), -160.);
/// assert_eq!(b[(3, 2)], -160. / 532.);
/// assert_eq!(a.cofactor(3, 2), 105.);
/// assert_eq!(b[(2, 3)], 105. / 532.);
/// assert_eq!(
/// b,
/// Matrix4x4::new(
/// [0.21804512, 0.45112783, 0.24060151, -0.04511278],
/// [-0.8082707, -1.456767, -0.44360903, 0.5206767],
/// [-0.078947365, -0.2236842, -0.05263158, 0.19736843],
/// [-0.52255636, -0.81390977, -0.30075186, 0.30639097]
/// )
/// );
///
/// // Second test case
/// assert_eq!(
/// Matrix4x4::new(
/// [8., -5., 9., 2.],
/// [7., 5., 6., 1.],
/// [-6., 0., 9., 6.],
/// [-3., 0., -9., -4.],
/// )
/// .inverse(),
/// Matrix4x4::new(
/// [-0.15384616, -0.15384616, -0.2820513, -0.53846157],
/// [-0.07692308, 0.12307692, 0.025641026, 0.03076923],
/// [0.35897437, 0.35897437, 0.43589744, 0.9230769],
/// [-0.6923077, -0.6923077, -0.7692308, -1.9230769]
/// ),
/// );
///
/// // Third test case
/// assert_eq!(
/// Matrix4x4::new(
/// [9., 3., 0., 9.],
/// [-5., -2., -6., -3.],
/// [-4., 9., 6., 4.],
/// [-7., 6., 6., 2.],
/// )
/// .inverse(),
/// Matrix4x4::new(
/// [-0.04074074, -0.07777778, 0.14444445, -0.22222222],
/// [-0.07777778, 0.033333335, 0.36666667, -0.33333334],
/// [-0.029012345, -0.14629629, -0.10925926, 0.12962963],
/// [0.17777778, 0.06666667, -0.26666668, 0.33333334]
/// ),
/// );
///
/// let a = Matrix4x4::new(
/// [3., -9., 7., 3.],
/// [3., -8., 2., -9.],
/// [-4., 4., 4., 1.],
/// [-6., 5., -1., 1.],
/// );
/// let b = Matrix4x4::new(
/// [8., 2., 2., 2.],
/// [3., -1., 7., 0.],
/// [7., 0., 5., 4.],
/// [6., -2., 0., 5.],
/// );
/// let c = a * b;
/// assert_eq!(c * b.inverse(), a);
/// ```
pub fn inverse(&self) -> Matrix4x4 { pub fn inverse(&self) -> Matrix4x4 {
self.inverse_rtc() self.inverse_rtc()
} }
@ -962,6 +551,209 @@ impl IndexMut<(usize, usize)> for Matrix4x4 {
mod tests { mod tests {
use super::*; use super::*;
use crate::{
float::consts::PI,
tuples::{point, vector},
};
#[test]
fn translation() {
let transform = Matrix4x4::translation(5., -3., 2.);
let p = point(-3., 4., 5.);
assert_eq!(transform * p, point(2., 1., 7.));
let inv = transform.inverse();
assert_eq!(inv * p, point(-8., 7., 3.));
let v = vector(-3., 4., 5.);
assert_eq!(transform * v, v);
}
#[test]
fn scaling() {
// A scaling matrix applied to a point.
let transform = Matrix4x4::scaling(2., 3., 4.);
let p = point(-4., 6., 8.);
assert_eq!(transform * p, point(-8., 18., 32.));
// A scaling matrix applied to a vector.
let v = vector(-4., 6., 8.);
assert_eq!(transform * v, vector(-8., 18., 32.));
// Multiplying by the inverse of a scaling matrix.
let inv = transform.inverse();
assert_eq!(inv * v, vector(-2., 2., 2.));
// Reflection is scaling by a negative value.
let transform = Matrix4x4::scaling(-1., 1., 1.);
let p = point(2., 3., 4.);
assert_eq!(transform * p, point(-2., 3., 4.));
}
#[test]
fn rotation_x() {
// A scaling matrix applied to a point.
let p = point(0., 1., 0.);
let half_quarter = Matrix4x4::rotation_x(PI / 4.);
let full_quarter = Matrix4x4::rotation_x(PI / 2.);
assert_eq!(
half_quarter * p,
point(0., (2.0 as Float).sqrt() / 2., (2.0 as Float).sqrt() / 2.)
);
assert_eq!(full_quarter * p, point(0., 0., 1.),);
}
#[test]
fn rotation_y() {
// A scaling matrix applied to a point.
let p = point(0., 0., 1.);
let half_quarter = Matrix4x4::rotation_y(PI / 4.);
let full_quarter = Matrix4x4::rotation_y(PI / 2.);
assert_eq!(
half_quarter * p,
point((2.0 as Float).sqrt() / 2., 0., (2.0 as Float).sqrt() / 2.)
);
assert_eq!(full_quarter * p, point(1., 0., 0.,),);
}
#[test]
fn rotation_z() {
// A scaling matrix applied to a point.
let p = point(0., 1., 0.);
let half_quarter = Matrix4x4::rotation_z(PI / 4.);
let full_quarter = Matrix4x4::rotation_z(PI / 2.);
assert_eq!(
half_quarter * p,
point(-(2.0 as Float).sqrt() / 2., (2.0 as Float).sqrt() / 2., 0.)
);
assert_eq!(full_quarter * p, point(-1., 0., 0.,),);
}
#[test]
fn transpose() {
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);
assert_eq!(Matrix4x4::identity(), Matrix4x4::identity().transpose());
}
#[test]
fn shearing() {
// A shearing transform moves x in proportion to y.
let transform = Matrix4x4::shearing(1., 0., 0., 0., 0., 0.);
let p = point(2., 3., 4.);
assert_eq!(transform * p, point(5., 3., 4.));
// A shearing transform moves x in proportion to z.
let transform = Matrix4x4::shearing(0., 1., 0., 0., 0., 0.);
let p = point(2., 3., 4.);
assert_eq!(transform * p, point(6., 3., 4.));
// A shearing transform moves y in proportion to x.
let transform = Matrix4x4::shearing(0., 0., 1., 0., 0., 0.);
let p = point(2., 3., 4.);
assert_eq!(transform * p, point(2., 5., 4.));
// A shearing transform moves y in proportion to z.
let transform = Matrix4x4::shearing(0., 0., 0., 1., 0., 0.);
let p = point(2., 3., 4.);
assert_eq!(transform * p, point(2., 7., 4.));
// A shearing transform moves z in proportion to x.
let transform = Matrix4x4::shearing(0., 0., 0., 0., 1., 0.);
let p = point(2., 3., 4.);
assert_eq!(transform * p, point(2., 3., 6.));
// A shearing transform moves z in proportion to y.
let transform = Matrix4x4::shearing(0., 0., 0., 0., 0., 1.);
let p = point(2., 3., 4.);
assert_eq!(transform * p, point(2., 3., 7.));
}
#[test]
fn inverse_rtiow() {
let i = Matrix4x4::identity();
assert_eq!(i.inverse_rtiow() * 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_rtiow() * m, i);
assert_eq!(m * m.inverse_rtiow(), i);
}
#[test]
fn determinant_2x2() {
let a = Matrix2x2::new([1., 5.], [-3., 2.]);
assert_eq!(a.determinant(), 17.);
}
#[test]
fn determinant_3x3() {
let a = Matrix3x3::new([1., 2., 6.], [-5., 8., -4.], [2., 6., 4.]);
assert_eq!(a.cofactor(0, 0), 56.);
assert_eq!(a.cofactor(0, 1), 12.);
assert_eq!(a.cofactor(0, 2), -46.);
assert_eq!(a.determinant(), -196.);
}
#[test]
fn determinant_4x4() {
let a = Matrix4x4::new(
[-2., -8., 3., 5.],
[-3., 1., 7., 3.],
[1., 2., -9., 6.],
[-6., 7., 7., -9.],
);
assert_eq!(a.cofactor(0, 0), 690.);
assert_eq!(a.cofactor(0, 1), 447.);
assert_eq!(a.cofactor(0, 2), 210.);
assert_eq!(a.cofactor(0, 3), 51.);
assert_eq!(a.determinant(), -4071.);
}
#[test]
fn submatrix_3x3() {
assert_eq!(
Matrix3x3::new([1., 5., 0.], [-3., 2., 7.], [0., 6., -3.],).submatrix(0, 2),
Matrix2x2::new([-3., 2.], [0., 6.])
);
}
#[test]
fn submatrix_4x4() {
assert_eq!(
Matrix4x4::new(
[-6., 1., 1., 6.],
[-8., 5., 8., 6.],
[-1., 0., 8., 2.],
[-7., 1., -1., 1.],
)
.submatrix(2, 1),
Matrix3x3::new([-6., 1., 6.], [-8., 8., 6.], [-7., -1., 1.],)
);
}
#[test]
fn minor_3x3() {
let a = Matrix3x3::new([3., 5., 0.], [2., -1., -7.], [6., -1., 5.]);
let b = a.submatrix(1, 0);
assert_eq!(b.determinant(), 25.0);
assert_eq!(b.determinant(), a.minor(1, 0));
}
#[test]
fn cofactor_3x3() {
let a = Matrix3x3::new([3., 5., 0.], [2., -1., -7.], [6., -1., 5.]);
assert_eq!(a.minor(0, 0), -12.);
assert_eq!(a.cofactor(0, 0), -12.);
assert_eq!(a.minor(1, 0), 25.);
assert_eq!(a.cofactor(1, 0), -25.);
}
#[test] #[test]
fn construct2x2() { fn construct2x2() {
let m = Matrix2x2::new([-3., 5.], [1., -2.]); let m = Matrix2x2::new([-3., 5.], [1., -2.]);
@ -980,6 +772,100 @@ mod tests {
assert_eq!(m[(2, 2)], 1.); assert_eq!(m[(2, 2)], 1.);
} }
#[test] #[test]
fn invertable() {
let a = Matrix4x4::new(
[6., 4., 4., 4.],
[5., 5., 7., 6.],
[4., -9., 3., -7.],
[9., 1., 7., -6.],
);
assert_eq!(a.determinant(), -2120.);
assert_eq!(a.invertable(), true);
let a = Matrix4x4::new(
[-4., 2., -2., -3.],
[9., 6., 2., 6.],
[0., -5., 1., -5.],
[0., 0., 0., 0.],
);
assert_eq!(a.determinant(), 0.);
assert_eq!(a.invertable(), false);
}
#[test]
fn inverse() {
let a = Matrix4x4::new(
[-5., 2., 6., -8.],
[1., -5., 1., 8.],
[7., 7., -6., -7.],
[1., -3., 7., 4.],
);
let b = a.inverse();
assert_eq!(a.determinant(), 532.);
assert_eq!(a.cofactor(2, 3), -160.);
assert_eq!(b[(3, 2)], -160. / 532.);
assert_eq!(a.cofactor(3, 2), 105.);
assert_eq!(b[(2, 3)], 105. / 532.);
assert_eq!(
b,
Matrix4x4::new(
[0.21804512, 0.45112783, 0.24060151, -0.04511278],
[-0.8082707, -1.456767, -0.44360903, 0.5206767],
[-0.078947365, -0.2236842, -0.05263158, 0.19736843],
[-0.52255636, -0.81390977, -0.30075186, 0.30639097]
)
);
// Second test case
assert_eq!(
Matrix4x4::new(
[8., -5., 9., 2.],
[7., 5., 6., 1.],
[-6., 0., 9., 6.],
[-3., 0., -9., -4.],
)
.inverse(),
Matrix4x4::new(
[-0.15384616, -0.15384616, -0.2820513, -0.53846157],
[-0.07692308, 0.12307692, 0.025641026, 0.03076923],
[0.35897437, 0.35897437, 0.43589744, 0.9230769],
[-0.6923077, -0.6923077, -0.7692308, -1.9230769]
),
);
// Third test case
assert_eq!(
Matrix4x4::new(
[9., 3., 0., 9.],
[-5., -2., -6., -3.],
[-4., 9., 6., 4.],
[-7., 6., 6., 2.],
)
.inverse(),
Matrix4x4::new(
[-0.04074074, -0.07777778, 0.14444445, -0.22222222],
[-0.07777778, 0.033333335, 0.36666667, -0.33333334],
[-0.029012345, -0.14629629, -0.10925926, 0.12962963],
[0.17777778, 0.06666667, -0.26666668, 0.33333334]
),
);
let a = Matrix4x4::new(
[3., -9., 7., 3.],
[3., -8., 2., -9.],
[-4., 4., 4., 1.],
[-6., 5., -1., 1.],
);
let b = Matrix4x4::new(
[8., 2., 2., 2.],
[3., -1., 7., 0.],
[7., 0., 5., 4.],
[6., -2., 0., 5.],
);
let c = a * b;
assert_eq!(c * b.inverse(), a);
}
#[test]
fn construct4x4() { fn construct4x4() {
let m = Matrix4x4::new( let m = Matrix4x4::new(
[1., 2., 3., 4.], [1., 2., 3., 4.],