retrofire_core/math/
mat.rs

1#![allow(clippy::needless_range_loop)]
2
3//! Matrices and linear transforms.
4
5use core::array;
6use core::fmt::{self, Debug, Formatter};
7use core::marker::PhantomData;
8use core::ops::Range;
9
10use crate::render::{NdcToScreen, ViewToProj};
11
12use super::space::{Linear, Proj4, Real};
13use super::vec::{ProjVec4, Vec2, Vec2u, Vec3, Vector};
14
15/// A linear transform from one space (or basis) to another.
16///
17/// This is a tag trait with no functionality in itself. It is used to
18/// statically ensure that only compatible maps can be composed, and that
19/// only compatible vectors can be transformed.
20pub trait LinearMap {
21    /// The source space, or domain, of `Self`.
22    type Source;
23    /// The destination space, or range, of `Self`.
24    type Dest;
25}
26
27/// Composition of two `LinearMap`s, `Self` ∘ `Inner`.
28/// If `Self` maps from `B` to `C`, and `Inner` maps from `A` to `B`,
29/// `Self::Result` maps from `A` to `C`.
30pub trait Compose<Inner: LinearMap>: LinearMap<Source = Inner::Dest> {
31    /// The result of composing `Self` with `Inner`.
32    type Result: LinearMap<Source = Inner::Source, Dest = Self::Dest>;
33}
34
35/// A change of basis in real vector space of dimension `DIM`.
36#[derive(Copy, Clone, Default, Eq, PartialEq)]
37pub struct RealToReal<const DIM: usize, SrcBasis = (), DstBasis = ()>(
38    PhantomData<(SrcBasis, DstBasis)>,
39);
40
41/// Mapping from real to projective space.
42#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)]
43pub struct RealToProj<SrcBasis>(PhantomData<SrcBasis>);
44
45/// A generic matrix type.
46#[repr(transparent)]
47#[derive(Copy, Eq, PartialEq)]
48pub struct Matrix<Repr, Map>(pub Repr, PhantomData<Map>);
49
50/// Type alias for a 3x3 float matrix.
51pub type Mat3x3<Map = ()> = Matrix<[[f32; 3]; 3], Map>;
52/// Type alias for a 4x4 float matrix.
53pub type Mat4x4<Map = ()> = Matrix<[[f32; 4]; 4], Map>;
54
55//
56// Inherent impls
57//
58
59impl<Repr, Map> Matrix<Repr, Map> {
60    /// Returns a matrix with the given elements.
61    #[inline]
62    pub const fn new(els: Repr) -> Self {
63        Self(els, PhantomData)
64    }
65
66    /// Returns a matrix equal to `self` but with mapping `M`.
67    ///
68    /// This method can be used to coerce a matrix to a different
69    /// mapping in case it is needed to make types match.
70    #[inline]
71    pub fn to<M>(&self) -> Matrix<Repr, M>
72    where
73        Repr: Clone,
74    {
75        Matrix::new(self.0.clone())
76    }
77}
78
79impl<Sc, const N: usize, const M: usize, Map> Matrix<[[Sc; N]; M], Map>
80where
81    Sc: Copy,
82    Map: LinearMap,
83{
84    /// Returns the row vector of `self` with index `i`.
85    /// The returned vector is in space `Map::Source`.
86    ///
87    /// # Panics
88    /// If `i >= M`.
89    #[inline]
90    pub fn row_vec(&self, i: usize) -> Vector<[Sc; N], Map::Source> {
91        Vector::new(self.0[i])
92    }
93    /// Returns the column vector of `self` with index `i`.
94    ///
95    /// The returned vector is in space `Map::Dest`.
96    ///
97    /// # Panics
98    /// If `i >= N`.
99    #[inline]
100    pub fn col_vec(&self, i: usize) -> Vector<[Sc; M], Map::Dest> {
101        Vector::new(self.0.map(|row| row[i]))
102    }
103}
104impl<Sc: Copy, const N: usize, const DIM: usize, S, D>
105    Matrix<[[Sc; N]; N], RealToReal<DIM, S, D>>
106{
107    /// Returns `self` with its rows and columns swapped.
108    pub fn transpose(self) -> Matrix<[[Sc; N]; N], RealToReal<DIM, D, S>> {
109        const { assert!(N >= DIM, "map dimension >= matrix dimension") }
110        array::from_fn(|j| array::from_fn(|i| self.0[i][j])).into()
111    }
112}
113
114impl<const N: usize, Map> Matrix<[[f32; N]; N], Map> {
115    /// Returns the `N`×`N` identity matrix.
116    pub const fn identity() -> Self {
117        let mut els = [[0.0; N]; N];
118        let mut i = 0;
119        while i < N {
120            els[i][i] = 1.0;
121            i += 1;
122        }
123        Self::new(els)
124    }
125
126    /// Returns whether every element of `self` is finite
127    /// (ie. neither `Inf`, `-Inf`, or `NaN`).
128    fn is_finite(&self) -> bool {
129        self.0.iter().flatten().all(|e| e.is_finite())
130    }
131}
132
133impl<M: LinearMap> Mat4x4<M> {
134    /// Constructs a matrix from a set of basis vectors.
135    pub const fn from_basis(i: Vec3, j: Vec3, k: Vec3) -> Self {
136        Self::new([
137            [i.0[0], i.0[1], i.0[2], 0.0],
138            [j.0[0], j.0[1], j.0[2], 0.0],
139            [k.0[0], k.0[1], k.0[2], 0.0],
140            [0.0, 0.0, 0.0, 1.0],
141        ])
142    }
143}
144
145impl<Sc, const N: usize, Map> Matrix<[[Sc; N]; N], Map>
146where
147    Sc: Linear<Scalar = Sc> + Copy,
148    Map: LinearMap,
149{
150    /// Returns the composite transform of `self` and `other`.
151    ///
152    /// Computes the matrix product of `self` and `other` such that applying
153    /// the resulting transformation is equivalent to first applying `other`
154    /// and then `self`. More succinctly,
155    /// ```text
156    /// (𝗠 ∘ 𝗡)𝘃 = 𝗠(𝗡𝘃)
157    /// ```
158    /// for some matrices 𝗠 and 𝗡 and a vector 𝘃.
159    #[must_use]
160    pub fn compose<Inner: LinearMap>(
161        &self,
162        other: &Matrix<[[Sc; N]; N], Inner>,
163    ) -> Matrix<[[Sc; N]; N], <Map as Compose<Inner>>::Result>
164    where
165        Map: Compose<Inner>,
166    {
167        let cols: [_; N] = array::from_fn(|i| other.col_vec(i));
168        array::from_fn(|j| {
169            let row = self.row_vec(j);
170            array::from_fn(|i| row.dot(&cols[i]))
171        })
172        .into()
173    }
174    /// Returns the composite transform of `other` and `self`.
175    ///
176    /// Computes the matrix product of `self` and `other` such that applying
177    /// the resulting matrix is equivalent to first applying `self` and then
178    /// `other`. The call `self.then(other)` is thus equivalent to
179    /// `other.compose(self)`.
180    #[must_use]
181    pub fn then<Outer: Compose<Map>>(
182        &self,
183        other: &Matrix<[[Sc; N]; N], Outer>,
184    ) -> Matrix<[[Sc; N]; N], <Outer as Compose<Map>>::Result> {
185        other.compose(self)
186    }
187}
188
189impl<Src, Dst> Mat3x3<RealToReal<2, Src, Dst>> {
190    /// Maps the real 2-vector 𝘃 from basis `Src` to basis `Dst`.
191    ///
192    /// Computes the matrix–vector multiplication 𝝡𝘃 where 𝘃 is interpreted as
193    /// a column vector with an implicit 𝘃<sub>2</sub> component with value 1:
194    ///
195    /// ```text
196    ///         / M00 ·  ·  \ / v0 \
197    ///  Mv  =  |  ·  ·  ·  | | v1 |  =  ( v0' v1' 1 )
198    ///         \  ·  · M22 / \  1 /
199    /// ```
200    #[must_use]
201    pub fn apply(&self, v: &Vec2<Src>) -> Vec2<Dst> {
202        let v = [v.x(), v.y(), 1.0].into();
203        array::from_fn(|i| self.row_vec(i).dot(&v)).into()
204    }
205}
206
207impl<Src, Dst> Mat4x4<RealToReal<3, Src, Dst>> {
208    /// Maps the real 3-vector 𝘃 from basis `Src` to basis `Dst`.
209    ///
210    /// Computes the matrix–vector multiplication 𝝡𝘃 where 𝘃 is interpreted as
211    /// a column vector with an implicit 𝘃<sub>3</sub> component with value 1:
212    ///
213    /// ```text
214    ///         / M00 ·  ·  ·  \ / v0 \
215    ///  Mv  =  |  ·  ·  ·  ·  | | v1 |  =  ( v0' v1' v2' 1 )
216    ///         |  ·  ·  ·  ·  | | v2 |
217    ///         \  ·  ·  · M33 / \  1 /
218    /// ```
219    #[must_use]
220    pub fn apply(&self, v: &Vec3<Src>) -> Vec3<Dst> {
221        let v = [v.x(), v.y(), v.z(), 1.0].into();
222        array::from_fn(|i| self.row_vec(i).dot(&v)).into()
223    }
224
225    /// Returns the determinant of `self`.
226    ///
227    /// Given a matrix M,
228    /// ```text
229    ///         / a  b  c  d \
230    ///  M  =   | e  f  g  h |
231    ///         | i  j  k  l |
232    ///         \ m  n  o  p /
233    /// ```
234    /// its determinant can be computed by recursively computing
235    /// the determinants of sub-matrices on rows 1.. and multiplying
236    /// them by the elements on row 0:
237    /// ```text
238    ///              | f g h |       | e g h |
239    /// det(M) = a · | j k l | - b · | i k l |  + - ···
240    ///              | n o p |       | m o p |
241    /// ```
242    pub fn determinant(&self) -> f32 {
243        let [[a, b, c, d], r, s, t] = self.0;
244        let det2 = |m, n| s[m] * t[n] - s[n] * t[m];
245        let det3 =
246            |j, k, l| r[j] * det2(k, l) - r[k] * det2(j, l) + r[l] * det2(j, k);
247
248        a * det3(1, 2, 3) - b * det3(0, 2, 3) + c * det3(0, 1, 3)
249            - d * det3(0, 1, 2)
250    }
251
252    /// Returns the inverse matrix of `self`.
253    ///
254    /// The inverse 𝝡<sup>-1</sup> of matrix 𝝡 is a matrix that, when
255    /// composed with 𝝡, results in the identity matrix:
256    ///
257    /// 𝝡 ∘ 𝝡<sup>-1</sup> = 𝝡<sup>-1</sup> ∘ 𝝡 = 𝐈
258    ///
259    /// In other words, it applies the transform of 𝝡 in reverse.
260    /// Given vectors 𝘃 and 𝘂,
261    ///
262    /// 𝝡𝘃 = 𝘂 ⇔ 𝝡<sup>-1</sup> 𝘂 = 𝘃.
263    ///
264    /// Only matrices with a nonzero determinant have a defined inverse.
265    /// A matrix without an inverse is said to be singular.
266    ///
267    /// Note: This method uses naive Gauss–Jordan elimination and may
268    /// suffer from imprecision or numerical instability in certain cases.
269    ///
270    /// # Panics
271    /// If `self` is singular or near-singular:
272    /// * Panics in debug mode.
273    /// * Does not panic in release mode, but the result may be inaccurate
274    /// or contain `Inf`s or `NaN`s.
275    #[must_use]
276    pub fn inverse(&self) -> Mat4x4<RealToReal<3, Dst, Src>> {
277        use super::float::f32;
278        if cfg!(debug_assertions) {
279            let det = self.determinant();
280            assert!(
281                f32::abs(det) > f32::EPSILON,
282                "a singular, near-singular, or non-finite matrix does not \
283                 have a well-defined inverse (determinant = {det})"
284            );
285        }
286
287        // Elementary row operation subtracting one row,
288        // multiplied by a scalar, from another
289        fn sub_row(m: &mut Mat4x4, from: usize, to: usize, mul: f32) {
290            m.0[to] = (m.row_vec(to) - m.row_vec(from) * mul).0;
291        }
292
293        // Elementary row operation multiplying one row with a scalar
294        fn mul_row(m: &mut Mat4x4, row: usize, mul: f32) {
295            m.0[row] = (m.row_vec(row) * mul).0;
296        }
297
298        // Elementary row operation swapping two rows
299        fn swap_rows(m: &mut Mat4x4, r: usize, s: usize) {
300            m.0.swap(r, s);
301        }
302
303        // This algorithm attempts to reduce `this` to the identity matrix
304        // by simultaneously applying elementary row operations to it and
305        // another matrix `inv` which starts as the identity matrix. Once
306        // `this` is reduced, the value of `inv` has become the inverse of
307        // `this` and thus of `self`.
308
309        let inv = &mut Mat4x4::identity();
310        let this = &mut self.to();
311
312        // Apply row operations to reduce the matrix to an upper echelon form
313        for idx in 0..4 {
314            let pivot = (idx..4)
315                .max_by(|&r1, &r2| {
316                    let v1 = f32::abs(this.0[r1][idx]);
317                    let v2 = f32::abs(this.0[r2][idx]);
318                    v1.partial_cmp(&v2).unwrap()
319                })
320                .unwrap();
321
322            if this.0[pivot][idx] != 0.0 {
323                swap_rows(this, idx, pivot);
324                swap_rows(inv, idx, pivot);
325
326                let div = 1.0 / this.0[idx][idx];
327                for r in (idx + 1)..4 {
328                    let x = this.0[r][idx] * div;
329                    sub_row(this, idx, r, x);
330                    sub_row(inv, idx, r, x);
331                }
332            }
333        }
334        // now in upper echelon form, back-substitute variables
335        for &idx in &[3, 2, 1] {
336            let diag = this.0[idx][idx];
337            for r in 0..idx {
338                let x = this.0[r][idx] / diag;
339
340                sub_row(this, idx, r, x);
341                sub_row(inv, idx, r, x);
342            }
343        }
344        // normalize
345        for r in 0..4 {
346            let x = 1.0 / this.0[r][r];
347            mul_row(this, r, x);
348            mul_row(inv, r, x);
349        }
350        debug_assert!(inv.is_finite());
351        inv.to()
352    }
353}
354
355impl<Src> Mat4x4<RealToProj<Src>> {
356    /// Maps the real 3-vector 𝘃 from basis 𝖡 to the projective 4-space.
357    ///
358    /// Computes the matrix–vector multiplication 𝝡𝘃 where 𝘃 is interpreted as
359    /// a column vector with an implicit 𝘃<sub>3</sub> component with value 1:
360    ///
361    /// ```text
362    ///         / M00  ·  · \ / v0 \
363    ///  Mv  =  |    ·      | | v1 |  =  ( v0' v1' v2' v3' )
364    ///         |      ·    | | v2 |
365    ///         \ ·  ·  M33 / \  1 /
366    /// ```
367    #[must_use]
368    pub fn apply(&self, v: &Vec3<Src>) -> ProjVec4 {
369        let v = Vector::from([v.x(), v.y(), v.z(), 1.0]);
370        [
371            self.row_vec(0).dot(&v),
372            self.row_vec(1).dot(&v),
373            self.row_vec(2).dot(&v),
374            self.row_vec(3).dot(&v),
375        ]
376        .into()
377    }
378}
379
380//
381// Local trait impls
382//
383
384impl<const DIM: usize, S, D> LinearMap for RealToReal<DIM, S, D> {
385    type Source = Real<DIM, S>;
386    type Dest = Real<DIM, D>;
387}
388
389impl<const DIM: usize, S, I, D> Compose<RealToReal<DIM, S, I>>
390    for RealToReal<DIM, I, D>
391{
392    type Result = RealToReal<DIM, S, D>;
393}
394
395impl<S> LinearMap for RealToProj<S> {
396    type Source = Real<3, S>;
397    type Dest = Proj4;
398}
399
400impl<S, I> Compose<RealToReal<3, S, I>> for RealToProj<I> {
401    type Result = RealToProj<S>;
402}
403
404/// Dummy `LinearMap` to help with generic code.
405impl LinearMap for () {
406    type Source = ();
407    type Dest = ();
408}
409
410//
411// Foreign trait impls
412//
413
414impl<R: Clone, M> Clone for Matrix<R, M> {
415    fn clone(&self) -> Self {
416        self.to()
417    }
418}
419
420impl<const N: usize, Map> Default for Matrix<[[f32; N]; N], Map> {
421    /// Returns the `N`×`N` identity matrix.
422    fn default() -> Self {
423        Self::identity()
424    }
425}
426
427impl<S: Debug, Map: Debug + Default, const N: usize, const M: usize> Debug
428    for Matrix<[[S; N]; M], Map>
429{
430    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
431        writeln!(f, "Matrix<{:?}>[", Map::default())?;
432        for i in 0..M {
433            writeln!(f, "    {:6.2?}", self.0[i])?;
434        }
435        write!(f, "]")
436    }
437}
438
439impl<const DIM: usize, F, T> Debug for RealToReal<DIM, F, T>
440where
441    F: Debug + Default,
442    T: Debug + Default,
443{
444    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
445        write!(f, "{:?}→{:?}", F::default(), T::default())
446    }
447}
448
449impl<Repr, M> From<Repr> for Matrix<Repr, M> {
450    fn from(repr: Repr) -> Self {
451        Self(repr, PhantomData)
452    }
453}
454
455//
456// Free functions
457//
458
459/// Returns a matrix applying a scaling by `s`.
460///
461/// Tip: use [`splat`][super::vec::splat] to scale uniformly:
462/// ```
463/// # use retrofire_core::math::{mat::scale, vec::splat};
464/// let m = scale(splat(2.0));
465/// assert_eq!(m.0[0][0], 2.0);
466/// assert_eq!(m.0[1][1], 2.0);
467/// assert_eq!(m.0[2][2], 2.0);
468/// ```
469pub const fn scale(s: Vec3) -> Mat4x4<RealToReal<3>> {
470    let [x, y, z] = s.0;
471    Matrix::new([
472        [x, 0.0, 0.0, 0.0],
473        [0.0, y, 0.0, 0.0],
474        [0.0, 0.0, z, 0.0],
475        [0.0, 0.0, 0.0, 1.0],
476    ])
477}
478
479/// Returns a matrix applying a translation by `t`.
480pub const fn translate(t: Vec3) -> Mat4x4<RealToReal<3>> {
481    let [x, y, z] = t.0;
482    Matrix::new([
483        [1.0, 0.0, 0.0, x],
484        [0.0, 1.0, 0.0, y],
485        [0.0, 0.0, 1.0, z],
486        [0.0, 0.0, 0.0, 1.0],
487    ])
488}
489
490/// Returns a matrix applying a rotation such that the original y axis
491/// is now parallel with `new_y` and the new z axis is orthogonal to
492/// both `x` and `new_y`.
493///
494/// Returns an orthogonal basis. If `new_y` and `x` are unit vectors,
495/// the result is orthonormal.
496#[cfg(feature = "fp")]
497pub fn orient_y(new_y: Vec3, x: Vec3) -> Mat4x4<RealToReal<3>> {
498    orient(new_y, x.cross(&new_y).normalize())
499}
500/// Returns a matrix applying a rotation such that the original z axis
501/// is now parallel with `new_z` and the new y axis is orthogonal to
502/// both `new_z` and `x`.
503///
504/// Returns an orthogonal basis. If `new_z` and `x` are unit vectors,
505/// the result is orthonormal.
506#[cfg(feature = "fp")]
507pub fn orient_z(new_z: Vec3, x: Vec3) -> Mat4x4<RealToReal<3>> {
508    orient(new_z.cross(&x).normalize(), new_z)
509}
510
511#[cfg(feature = "fp")]
512fn orient(new_y: Vec3, new_z: Vec3) -> Mat4x4<RealToReal<3>> {
513    use crate::math::{space::Linear, ApproxEq};
514
515    assert!(!new_y.approx_eq(&Vec3::zero()));
516    assert!(!new_z.approx_eq(&Vec3::zero()));
517
518    let new_x = new_y.cross(&new_z);
519    Mat4x4::from_basis(new_x, new_y, new_z)
520}
521
522// TODO constify rotate_* functions once we have const trig functions
523
524/// Returns a matrix applying a rotation by `a` about the x axis.
525#[cfg(feature = "fp")]
526pub fn rotate_x(a: super::angle::Angle) -> Mat4x4<RealToReal<3>> {
527    let (sin, cos) = a.sin_cos();
528    [
529        [1.0, 0.0, 0.0, 0.0],
530        [0.0, cos, sin, 0.0],
531        [0.0, -sin, cos, 0.0],
532        [0.0, 0.0, 0.0, 1.0],
533    ]
534    .into()
535}
536/// Returns a matrix applying a rotation by `a` about the y axis.
537#[cfg(feature = "fp")]
538pub fn rotate_y(a: super::angle::Angle) -> Mat4x4<RealToReal<3>> {
539    let (sin, cos) = a.sin_cos();
540    [
541        [cos, 0.0, -sin, 0.0],
542        [0.0, 1.0, 0.0, 0.0],
543        [sin, 0.0, cos, 0.0],
544        [0.0, 0.0, 0.0, 1.0],
545    ]
546    .into()
547}
548/// Returns a matrix applying a rotation of angle `a` about the z axis.
549#[cfg(feature = "fp")]
550pub fn rotate_z(a: super::angle::Angle) -> Mat4x4<RealToReal<3>> {
551    let (sin, cos) = a.sin_cos();
552    [
553        [cos, sin, 0.0, 0.0],
554        [-sin, cos, 0.0, 0.0],
555        [0.0, 0.0, 1.0, 0.0],
556        [0.0, 0.0, 0.0, 1.0],
557    ]
558    .into()
559}
560
561/// Creates a perspective projection matrix.
562///
563/// # Parameters
564/// * `focal_ratio`: Focal length/aperture ratio. Larger values mean
565/// a smaller angle of view, with 1.0 corresponding to a horizontal
566/// field of view of 90 degrees.
567/// * `aspect_ratio`: Viewport width/height ratio. Larger values mean
568/// a wider field of view.
569/// * `near_far`: Depth range between the near and far clipping planes.
570/// Objects outside this range are clipped or culled.
571///
572/// # Panics
573/// * If any parameter value is nonpositive.
574/// * If `near_far` is an empty range.
575pub fn perspective(
576    focal_ratio: f32,
577    aspect_ratio: f32,
578    near_far: Range<f32>,
579) -> Mat4x4<ViewToProj> {
580    let (near, far) = (near_far.start, near_far.end);
581
582    assert!(focal_ratio > 0.0, "focal ratio must be positive");
583    assert!(aspect_ratio > 0.0, "aspect ratio must be positive");
584    assert!(near > 0.0, "near must be positive");
585    assert!(!near_far.is_empty(), "far must be greater than near");
586
587    let e00 = focal_ratio;
588    let e11 = e00 * aspect_ratio;
589    let e22 = (far + near) / (far - near);
590    let e23 = 2.0 * far * near / (near - far);
591    [
592        [e00, 0.0, 0.0, 0.0],
593        [0.0, e11, 0.0, 0.0],
594        [0.0, 0.0, e22, e23],
595        [0.0, 0.0, 1.0, 0.0],
596    ]
597    .into()
598}
599
600/// Creates an orthographic projection matrix.
601///
602/// # Parameters
603/// * `lbn`: The left-bottom-near corner of the projection box.
604/// * `rtf`: The right-bottom-far corner of the projection box.
605pub fn orthographic(lbn: Vec3, rtf: Vec3) -> Mat4x4<ViewToProj> {
606    let [dx, dy, dz] = (rtf - lbn).0;
607    let [sx, sy, sz] = (rtf + lbn).0;
608    [
609        [2.0 / dx, 0.0, 0.0, -sx / dx],
610        [0.0, 2.0 / dy, 0.0, -sy / dy],
611        [0.0, 0.0, 2.0 / dz, -sz / dz],
612        [0.0, 0.0, 0.0, 1.0],
613    ]
614    .into()
615}
616
617/// Creates a viewport transform matrix. A viewport matrix is used to
618/// transform points from the NDC space to screen space for rasterization.
619///
620/// # Parameters
621/// * `bounds`: the left-top and right-bottom coordinates of the viewport.
622pub fn viewport(bounds: Range<Vec2u>) -> Mat4x4<NdcToScreen> {
623    let Range { start, end } = bounds;
624    let h = (end.x() - start.x()) as f32 / 2.0;
625    let v = (end.y() - start.y()) as f32 / 2.0;
626    [
627        [h, 0.0, 0.0, h + start.x() as f32],
628        [0.0, v, 0.0, v + start.y() as f32],
629        [0.0, 0.0, 1.0, 0.0],
630        [0.0, 0.0, 0.0, 1.0],
631    ]
632    .into()
633}
634
635#[cfg(test)]
636mod tests {
637    use crate::assert_approx_eq;
638    use crate::math::vec::{splat, vec2, vec3};
639
640    #[cfg(feature = "fp")]
641    use crate::math::angle::degs;
642
643    use super::*;
644
645    #[derive(Debug, Default, Eq, PartialEq)]
646    struct Basis1;
647    #[derive(Debug, Default, Eq, PartialEq)]
648    struct Basis2;
649
650    type Map<const N: usize = 3> = RealToReal<N, Basis1, Basis2>;
651    type InvMap<const N: usize = 3> = RealToReal<N, Basis2, Basis1>;
652
653    mod mat3x3 {
654        use super::*;
655
656        const MAT: Mat3x3<Map> = Matrix::new([
657            [0.0, 1.0, 2.0], //
658            [10.0, 11.0, 12.0],
659            [20.0, 21.0, 22.0],
660        ]);
661
662        #[test]
663        fn row_col_vecs() {
664            assert_eq!(MAT.row_vec(2), vec3::<_, Basis1>(20.0, 21.0, 22.0));
665            assert_eq!(MAT.col_vec(2), vec3::<_, Basis2>(2.0, 12.0, 22.0));
666        }
667
668        #[test]
669        fn composition() {
670            let t = Mat3x3::<Map<2>>::new([
671                [1.0, 0.0, 2.0], //
672                [0.0, 1.0, -3.0],
673                [0.0, 0.0, 1.0],
674            ]);
675            let s = Mat3x3::<InvMap<2>>::new([
676                [-1.0, 0.0, 0.0], //
677                [0.0, 2.0, 0.0],
678                [0.0, 0.0, 1.0],
679            ]);
680
681            let ts = t.then(&s);
682            let st = s.then(&t);
683
684            assert_eq!(ts, s.compose(&t));
685            assert_eq!(st, t.compose(&s));
686
687            assert_eq!(ts.apply(&splat(0.0)), vec2(-2.0, -6.0));
688            assert_eq!(st.apply(&splat(0.0)), vec2(2.0, -3.0));
689        }
690
691        #[test]
692        fn scaling() {
693            let m = Mat3x3::<Map<2>>::new([
694                [2.0, 0.0, 0.0], //
695                [0.0, -3.0, 0.0],
696                [0.0, 0.0, 1.0],
697            ]);
698            assert_eq!(m.apply(&vec2(1.0, 2.0)), vec2(2.0, -6.0));
699        }
700
701        #[test]
702        fn translation() {
703            let m = Mat3x3::<Map<2>>::new([
704                [1.0, 0.0, 2.0], //
705                [0.0, 1.0, -3.0],
706                [0.0, 0.0, 1.0],
707            ]);
708            assert_eq!(m.apply(&vec2(1.0, 2.0)), vec2(3.0, -1.0));
709        }
710
711        #[test]
712        fn matrix_debug() {
713            assert_eq!(
714                alloc::format!("{MAT:?}"),
715                r#"Matrix<Basis1→Basis2>[
716    [  0.00,   1.00,   2.00]
717    [ 10.00,  11.00,  12.00]
718    [ 20.00,  21.00,  22.00]
719]"#
720            );
721        }
722    }
723
724    mod mat4x4 {
725        use super::*;
726
727        const MAT: Mat4x4<Map> = Matrix::new([
728            [0.0, 1.0, 2.0, 3.0],
729            [10.0, 11.0, 12.0, 13.0],
730            [20.0, 21.0, 22.0, 23.0],
731            [30.0, 31.0, 32.0, 33.0],
732        ]);
733
734        #[test]
735        fn row_col_vecs() {
736            assert_eq!(MAT.row_vec(1), [10.0, 11.0, 12.0, 13.0].into());
737            assert_eq!(MAT.col_vec(3), [3.0, 13.0, 23.0, 33.0].into());
738        }
739
740        #[test]
741        fn composition() {
742            let t = translate(vec3(1.0, 2.0, 3.0)).to::<Map>();
743            let s = scale(vec3(3.0, 2.0, 1.0)).to::<InvMap>();
744
745            let ts = t.then(&s);
746            let st = s.then(&t);
747
748            assert_eq!(ts, s.compose(&t));
749            assert_eq!(st, t.compose(&s));
750
751            assert_eq!(ts.apply(&splat(0.0)), vec3::<_, Basis1>(3.0, 4.0, 3.0));
752            assert_eq!(st.apply(&splat(0.0)), vec3::<_, Basis2>(1.0, 2.0, 3.0));
753        }
754
755        #[test]
756        fn scaling_vec3() {
757            let m = scale(vec3(1.0, -2.0, 3.0));
758            let v = vec3(0.0, 4.0, -3.0);
759            assert_eq!(m.apply(&v), vec3(0.0, -8.0, -9.0));
760        }
761
762        #[test]
763        fn translation_vec3() {
764            let m = translate(vec3(1.0, 2.0, 3.0));
765            let v = vec3(0.0, 5.0, -3.0);
766            assert_eq!(m.apply(&v), vec3(1.0, 7.0, 0.0));
767        }
768
769        #[cfg(feature = "fp")]
770        #[test]
771        fn rotation_x() {
772            let m = rotate_x(degs(90.0));
773            assert_eq!(m.apply(&splat(0.0)), splat(0.0));
774            assert_approx_eq!(
775                m.apply(&vec3(0.0, 0.0, 1.0)),
776                vec3(0.0, 1.0, 0.0)
777            );
778        }
779
780        #[cfg(feature = "fp")]
781        #[test]
782        fn rotation_y() {
783            let m = rotate_y(degs(90.0));
784            assert_eq!(m.apply(&splat(0.0)), splat(0.0));
785            assert_approx_eq!(
786                m.apply(&vec3(1.0, 0.0, 0.0)),
787                vec3(0.0, 0.0, 1.0)
788            );
789        }
790
791        #[cfg(feature = "fp")]
792        #[test]
793        fn rotation_z() {
794            let m = rotate_z(degs(90.0));
795            assert_eq!(m.apply(&splat(0.0)), splat(0.0));
796            assert_approx_eq!(
797                m.apply(&vec3(0.0, 1.0, 0.0)),
798                vec3(1.0, 0.0, 0.0)
799            );
800        }
801
802        #[test]
803        fn matrix_debug() {
804            assert_eq!(
805                alloc::format!("{MAT:?}"),
806                r#"Matrix<Basis1→Basis2>[
807    [  0.00,   1.00,   2.00,   3.00]
808    [ 10.00,  11.00,  12.00,  13.00]
809    [ 20.00,  21.00,  22.00,  23.00]
810    [ 30.00,  31.00,  32.00,  33.00]
811]"#
812            );
813        }
814    }
815
816    #[test]
817    fn transposition() {
818        let m = Matrix::<_, Map>::new([
819            [0.0, 1.0, 2.0], //
820            [10.0, 11.0, 12.0],
821            [20.0, 21.0, 22.0],
822        ]);
823        assert_eq!(
824            m.transpose(),
825            Matrix::<_, InvMap>::new([
826                [0.0, 10.0, 20.0], //
827                [1.0, 11.0, 21.0],
828                [2.0, 12.0, 22.0],
829            ])
830        );
831    }
832
833    #[test]
834    fn determinant_of_identity_is_one() {
835        let id: Mat4x4<Map> = Mat4x4::identity();
836        assert_eq!(id.determinant(), 1.0);
837    }
838
839    #[test]
840    fn determinant_of_scaling_is_product_of_diagonal() {
841        let scale: Mat4x4<_> = scale(vec3(2.0, 3.0, 4.0));
842        assert_eq!(scale.determinant(), 24.0);
843    }
844
845    #[cfg(feature = "fp")]
846    #[test]
847    fn determinant_of_rotation_is_one() {
848        let rot = rotate_x(degs(73.0)).then(&rotate_y(degs(-106.0)));
849        assert_approx_eq!(rot.determinant(), 1.0);
850    }
851
852    #[cfg(feature = "fp")]
853    #[test]
854    fn mat_times_mat_inverse_is_identity() {
855        let m = translate(vec3(1.0e3, -2.0e2, 0.0))
856            .then(&scale(vec3(0.5, 100.0, 42.0)))
857            .to::<Map>();
858
859        let m_inv: Mat4x4<InvMap> = m.inverse();
860
861        assert_eq!(
862            m.compose(&m_inv),
863            Mat4x4::<RealToReal<3, Basis2, Basis2>>::identity()
864        );
865        assert_eq!(
866            m_inv.compose(&m),
867            Mat4x4::<RealToReal<3, Basis1, Basis1>>::identity()
868        );
869    }
870
871    #[test]
872    fn inverse_reverts_transform() {
873        let m: Mat4x4<Map> = scale(vec3(1.0, 2.0, 0.5))
874            .then(&translate(vec3(-2.0, 3.0, 0.0)))
875            .to();
876        let m_inv: Mat4x4<InvMap> = m.inverse();
877
878        let v1: Vec3<Basis1> = vec3(1.0, -2.0, 3.0);
879        let v2: Vec3<Basis2> = vec3(2.0, 0.0, -2.0);
880
881        assert_eq!(m_inv.apply(&m.apply(&v1)), v1);
882        assert_eq!(m.apply(&m_inv.apply(&v2)), v2);
883    }
884
885    #[test]
886    fn orthographic_box_maps_to_unit_cube() {
887        let lbn = vec3(-20.0, 0.0, 0.01);
888        let rtf = vec3(100.0, 50.0, 100.0);
889
890        let m = orthographic(lbn, rtf);
891
892        assert_approx_eq!(m.apply(&lbn.to()), [-1.0, -1.0, -1.0, 1.0].into());
893        assert_approx_eq!(m.apply(&rtf.to()), splat(1.0));
894    }
895
896    #[test]
897    fn perspective_frustum_maps_to_unit_cube() {
898        let left_bot_near = vec3(-0.125, -0.0625, 0.1);
899        let right_top_far = vec3(125.0, 62.5, 100.0);
900
901        let m = perspective(0.8, 2.0, 0.1..100.0);
902
903        let lbn = m.apply(&left_bot_near);
904        assert_approx_eq!(lbn / lbn.w(), [-1.0, -1.0, -1.0, 1.0].into());
905
906        let rtf = m.apply(&right_top_far);
907        assert_approx_eq!(rtf / rtf.w(), splat(1.0));
908    }
909}
retrofire_core/math/mat.rs

retrofire_core/math/
mat.rs
