all_is_cubes_base/math/
matrix.rs

1//! Integer-coordinate matrices.
2//! This module is private but reexported by its parent.
3
4use core::cmp::Ordering;
5use core::ops;
6
7use euclid::Vector3D;
8use num_traits::One;
9
10/// Acts as polyfill for float methods
11#[cfg(not(feature = "std"))]
12#[allow(unused_imports)]
13use num_traits::float::FloatCore as _;
14
15use crate::math::{
16    Axis, Cube, Face6, Face7, FreeCoordinate, GridCoordinate, GridPoint, GridRotation, GridVector,
17    Gridgid,
18};
19
20/// A 4×3 affine transformation matrix in [`GridCoordinate`]s.
21//---
22// Design note: It would be nice to just use `euclid` types, but that doesn't get us the
23// exact semantics we want, particularly checked arithmetic.
24#[expect(clippy::exhaustive_structs)]
25#[derive(Clone, Copy, Debug, Eq, Hash, PartialEq)]
26pub struct GridMatrix {
27    /// First column
28    pub x: GridVector,
29    /// Second column
30    pub y: GridVector,
31    /// Third column
32    pub z: GridVector,
33    /// Fourth column (translation)
34    pub w: GridVector,
35}
36
37impl GridMatrix {
38    /// The zero matrix, which transforms all points to zero.
39    pub const ZERO: Self = Self {
40        x: Vector3D::new(0, 0, 0),
41        y: Vector3D::new(0, 0, 0),
42        z: Vector3D::new(0, 0, 0),
43        w: Vector3D::new(0, 0, 0),
44    };
45
46    /// For Y-down drawing
47    #[doc(hidden)] // used by all-is-cubes-content - TODO: public?
48    pub const FLIP_Y: Self = Self {
49        x: Vector3D::new(1, 0, 0),
50        y: Vector3D::new(0, -1, 0),
51        z: Vector3D::new(0, 0, 1),
52        w: Vector3D::new(0, 0, 0),
53    };
54
55    /// Note: This takes the elements in a column-major ordering, so the argument order
56    /// is transposed relative to a conventional textual display of a matrix.
57    #[expect(clippy::too_many_arguments)]
58    #[inline]
59    pub const fn new(
60        x0: GridCoordinate,
61        x1: GridCoordinate,
62        x2: GridCoordinate,
63        y0: GridCoordinate,
64        y1: GridCoordinate,
65        y2: GridCoordinate,
66        z0: GridCoordinate,
67        z1: GridCoordinate,
68        z2: GridCoordinate,
69        w0: GridCoordinate,
70        w1: GridCoordinate,
71        w2: GridCoordinate,
72    ) -> Self {
73        Self {
74            x: Vector3D::new(x0, x1, x2),
75            y: Vector3D::new(y0, y1, y2),
76            z: Vector3D::new(z0, z1, z2),
77            w: Vector3D::new(w0, w1, w2),
78        }
79    }
80
81    /// Construct a translation matrix.
82    #[inline]
83    pub fn from_translation(offset: impl Into<GridVector>) -> Self {
84        Self {
85            w: offset.into(),
86            ..Self::one()
87        }
88    }
89
90    /// Construct a uniform scaling matrix.
91    ///
92    /// Note that since this is an integer matrix, there is no possibility  of scaling less
93    /// than 1 other than 0!
94    #[inline]
95    pub fn from_scale(scale: GridCoordinate) -> Self {
96        Self {
97            x: Vector3D::new(scale, 0, 0),
98            y: Vector3D::new(0, scale, 0),
99            z: Vector3D::new(0, 0, scale),
100            w: Vector3D::new(0, 0, 0),
101        }
102    }
103
104    /// Construct a transformation to a translated and rotated coordinate system from
105    /// an origin in the target coordinate system and basis vectors expressed as [`Face7`]s.
106    ///
107    /// Skews or scaling cannot be performed using this constructor.
108    ///
109    /// ```
110    /// # extern crate all_is_cubes_base as all_is_cubes;
111    /// use all_is_cubes::math::{Face7::*, GridMatrix, GridPoint};
112    ///
113    /// let transform = GridMatrix::from_origin([10, 10, 10], PX, PZ, NY);
114    /// assert_eq!(
115    ///     transform.transform_point(GridPoint::new(1, 2, 3)),
116    ///     GridPoint::new(11, 7, 12),
117    /// );
118    /// ```
119    #[inline]
120    pub fn from_origin(origin: impl Into<GridPoint>, x: Face7, y: Face7, z: Face7) -> Self {
121        Self {
122            x: x.normal_vector(),
123            y: y.normal_vector(),
124            z: z.normal_vector(),
125            w: origin.into().to_vector(),
126        }
127    }
128
129    /// Convert this integer-valued matrix to an equivalent float-valued matrix.
130    #[inline]
131    #[rustfmt::skip]
132    pub fn to_free(self) -> euclid::Transform3D<FreeCoordinate, Cube, Cube> {
133        euclid::Transform3D::new(
134            self.x.x, self.x.y, self.x.z, 0,
135            self.y.x, self.y.y, self.y.z, 0,
136            self.z.x, self.z.y, self.z.z, 0,
137            self.w.x, self.w.y, self.w.z, 1,
138        ).cast()
139    }
140
141    /// Returns row `r` of the matrix.
142    #[inline(always)]
143    fn row(&self, r: Axis) -> MVector4<GridCoordinate> {
144        MVector4::new(self.x[r], self.y[r], self.z[r], self.w[r])
145    }
146
147    /// Equivalent to temporarily applying an offset of `[0.5, 0.5, 0.5]` while
148    /// transforming `cube` as per [`GridMatrix::transform_point`], despite the fact that
149    /// integer arithmetic is being used.
150    ///
151    /// This operation thus transforms the standard positive-octant unit cube identified
152    /// by its most negative corner the same way as the [`GridAab::single_cube`] containing
153    /// that cube.
154    ///
155    /// ```
156    /// # extern crate all_is_cubes_base as all_is_cubes;
157    /// use all_is_cubes::math::{Cube, Face7::*, GridMatrix, GridPoint};
158    ///
159    /// // Translation without rotation has the usual definition.
160    /// let matrix = GridMatrix::from_translation([10, 0, 0]);
161    /// assert_eq!(matrix.transform_cube(Cube::new(1, 1, 1)), Cube::new(11, 1, 1));
162    ///
163    /// // With a rotation or reflection, the results are different.
164    /// // TODO: Come up with a better example and explanation.
165    /// let reflected = GridMatrix::from_origin([10, 0, 0], NX, PY, PZ);
166    /// assert_eq!(reflected.transform_point(GridPoint::new(1, 5, 5)), GridPoint::new(9, 5, 5));
167    /// assert_eq!(reflected.transform_cube(Cube::new(1, 5, 5)), Cube::new(8, 5, 5));
168    /// ```
169    ///
170    /// [`GridAab::single_cube`]: crate::math::GridAab::single_cube
171    #[inline]
172    pub fn transform_cube(&self, cube: Cube) -> Cube {
173        Cube::from(
174            self.transform_point(cube.lower_bounds())
175                .min(self.transform_point(cube.upper_bounds())),
176        )
177    }
178
179    /// Decomposes a matrix into its rotation and translation components, stored in a
180    /// [`Gridgid`]. Returns `None` if the matrix has any scaling or skew.
181    ///
182    /// ```
183    /// # extern crate all_is_cubes_base as all_is_cubes;
184    /// use all_is_cubes::math::{Face6::*, Gridgid, GridMatrix, GridRotation, GridVector};
185    ///
186    /// assert_eq!(
187    ///     GridMatrix::new(
188    ///         0, -1,  0,
189    ///         1,  0,  0,
190    ///         0,  0,  1,
191    ///         7,  3, -8,
192    ///     ).decompose(),
193    ///     Some(Gridgid {
194    ///         rotation: GridRotation::from_basis([NY, PX, PZ]),
195    ///         translation: GridVector::new(7, 3, -8),
196    ///     }),
197    /// );
198    /// ```
199    #[allow(clippy::missing_inline_in_public_items)]
200    pub fn decompose(self) -> Option<Gridgid> {
201        Some(Gridgid {
202            rotation: GridRotation::from_basis([
203                Face6::try_from(self.x).ok()?,
204                Face6::try_from(self.y).ok()?,
205                Face6::try_from(self.z).ok()?,
206            ]),
207            translation: self.w,
208        })
209    }
210
211    /// Transform (rotate and scale) the given vector.
212    /// The translation part of this matrix is ignored.
213    #[inline]
214    #[must_use]
215    pub fn transform_vector(&self, vec: GridVector) -> GridVector {
216        GridVector::new(
217            self.row(Axis::X).truncate().dot(vec),
218            self.row(Axis::Y).truncate().dot(vec),
219            self.row(Axis::Z).truncate().dot(vec),
220        )
221    }
222
223    /// Transform the given point by this matrix.
224    #[inline]
225    #[must_use]
226    pub fn transform_point(&self, point: GridPoint) -> GridPoint {
227        let homogeneous = MVector4::homogeneous(point.to_vector());
228        GridPoint::new(
229            self.row(Axis::X).dot(homogeneous),
230            self.row(Axis::Y).dot(homogeneous),
231            self.row(Axis::Z).dot(homogeneous),
232        )
233    }
234
235    /// ```
236    /// # extern crate all_is_cubes_base as all_is_cubes;
237    /// use all_is_cubes::math::{GridMatrix, GridPoint};
238    ///
239    /// let transform_1 = GridMatrix::new(
240    ///     0, -1, 0,
241    ///     1, 0, 0,
242    ///     0, 0, 1,
243    ///     0, 0, 0,
244    /// );
245    /// let transform_2 = GridMatrix::from_translation([10, 20, 30]);
246    ///
247    /// // Demonstrate the directionality of concatenation.
248    /// assert_eq!(
249    ///     transform_1.concat(&transform_2).transform_point(GridPoint::new(0, 3, 0)),
250    ///     transform_1.transform_point(transform_2.transform_point(GridPoint::new(0, 3, 0))),
251    /// );
252    /// ```
253    #[must_use]
254    #[inline]
255    pub fn concat(&self, other: &Self) -> Self {
256        GridMatrix {
257            x: self.transform_vector(other.x),
258            y: self.transform_vector(other.y),
259            z: self.transform_vector(other.z),
260            // Shenanigan because we're working in 4x3 rather than 4x4
261            // with homogeneous coordinate vectors.
262            w: self.transform_point(other.w.to_point()).to_vector(),
263        }
264    }
265
266    /// Invert this matrix. Returns [`None`] if it is not invertible.
267    #[allow(clippy::missing_inline_in_public_items)]
268    pub fn inverse_transform(&self) -> Option<Self> {
269        // For now, implement this the expensive but simple way of borrowing float matrix ops.
270
271        const INVERSE_EPSILON: FreeCoordinate = 0.25 / (GridCoordinate::MAX as FreeCoordinate);
272        fn try_round(v: [FreeCoordinate; 4], expected_w: FreeCoordinate) -> Option<GridVector> {
273            let mut result = Vector3D::zero();
274            #[expect(clippy::needless_range_loop)]
275            for axis in 0..4 {
276                let rounded = v[axis].round();
277                let remainder = v[axis] - rounded;
278                if remainder.abs().partial_cmp(&INVERSE_EPSILON) != Some(Ordering::Less) {
279                    // The inverse matrix has a non-integer element.
280                    return None;
281                }
282                match axis {
283                    // TODO: check for overflow
284                    0 => result.x = rounded as GridCoordinate,
285                    1 => result.y = rounded as GridCoordinate,
286                    2 => result.z = rounded as GridCoordinate,
287                    3 =>
288                    {
289                        #[expect(clippy::float_cmp)]
290                        if rounded != expected_w {
291                            return None;
292                        }
293                    }
294                    _ => unreachable!(),
295                }
296            }
297            Some(result)
298        }
299
300        let fi = self.to_free().inverse()?.to_arrays();
301        Some(GridMatrix {
302            x: try_round(fi[0], 0.0)?,
303            y: try_round(fi[1], 0.0)?,
304            z: try_round(fi[2], 0.0)?,
305            w: try_round(fi[3], 1.0)?,
306        })
307    }
308}
309
310impl ops::Mul<Self> for GridMatrix {
311    type Output = Self;
312    #[inline]
313    fn mul(self, rhs: Self) -> Self::Output {
314        // Delegate to Transform implementation
315        self.concat(&rhs)
316    }
317}
318
319impl One for GridMatrix {
320    #[inline]
321    #[rustfmt::skip]
322    fn one() -> Self {
323        Self::new(
324            1, 0, 0,
325            0, 1, 0,
326            0, 0, 1,
327            0, 0, 0,
328        )
329    }
330}
331
332/// Mini 4D vector for matrix rows.
333#[derive(Clone, Copy, Debug, Eq, Hash, PartialEq)]
334struct MVector4<T> {
335    x: T,
336    y: T,
337    z: T,
338    w: T,
339}
340
341impl<T: Copy> MVector4<T> {
342    fn new(x: T, y: T, z: T, w: T) -> Self {
343        Self { x, y, z, w }
344    }
345
346    fn homogeneous<U>(v: Vector3D<T, U>) -> MVector4<T>
347    where
348        T: One,
349    {
350        Self {
351            x: v.x,
352            y: v.y,
353            z: v.z,
354            w: T::one(),
355        }
356    }
357
358    fn dot(self, other: Self) -> T
359    where
360        T: ops::Mul<Output = T>,
361        T: ops::Add<Output = T>,
362    {
363        self.x * other.x + self.y * other.y + self.z * other.z + self.w * other.w
364    }
365
366    fn truncate<U>(self) -> Vector3D<T, U> {
367        Vector3D::new(self.x, self.y, self.z)
368    }
369}
370
371#[cfg(test)]
372mod tests {
373    use super::*;
374    use euclid::{Transform3D, point3, vec3};
375    use rand::{Rng, SeedableRng as _};
376    use rand_xoshiro::Xoshiro256Plus;
377
378    fn random_grid_matrix(mut rng: impl Rng) -> GridMatrix {
379        let mut r = || rng.random_range(-100..=100);
380        GridMatrix::new(r(), r(), r(), r(), r(), r(), r(), r(), r(), r(), r(), r())
381    }
382
383    fn random_possibly_invertible_matrix(mut rng: impl Rng) -> GridMatrix {
384        let mut r = |n: GridCoordinate| rng.random_range(-n..=n);
385        GridMatrix::new(
386            r(1),
387            r(1),
388            r(1),
389            r(1),
390            r(1),
391            r(1),
392            r(1),
393            r(1),
394            r(1),
395            r(GridCoordinate::MAX / 10),
396            r(GridCoordinate::MAX / 10),
397            r(GridCoordinate::MAX / 10),
398        )
399    }
400
401    #[test]
402    #[rustfmt::skip]
403    fn equivalent_constructor() {
404        let m = GridMatrix::new(
405            1, 2, 3,
406            5, 6, 7,
407            9, 10, 11,
408            13, 14, 15,
409        );
410        assert_eq!(m.to_free(), Transform3D::new(
411            1., 2., 3., 0.,
412            5., 6., 7., 0.,
413            9., 10., 11., 0.,
414            13., 14., 15., 1.,
415        ));
416    }
417
418    #[test]
419    fn equivalent_transform() {
420        let mut rng = Xoshiro256Plus::seed_from_u64(2897358920346590823);
421        for _ in 1..100 {
422            let m = random_grid_matrix(&mut rng);
423            dbg!(m, m.to_free());
424            assert_eq!(
425                m.transform_point(GridPoint::new(2, 300, 40000)).map(FreeCoordinate::from),
426                m.to_free().transform_point3d(point3(2., 300., 40000.)).unwrap(),
427            );
428            assert_eq!(
429                m.transform_vector(GridVector::new(10, 20, 30)).map(FreeCoordinate::from),
430                m.to_free().transform_vector3d(vec3(10., 20., 30.)),
431            );
432        }
433    }
434
435    #[test]
436    fn equivalent_concat() {
437        let mut rng = Xoshiro256Plus::seed_from_u64(5933089223468901296);
438        for _ in 1..100 {
439            let m1 = random_grid_matrix(&mut rng);
440            let m2 = random_grid_matrix(&mut rng);
441            // our concat() orderingis inherited from cgmath which is opposite of euclid then()
442            assert_eq!(m1.concat(&m2).to_free(), m2.to_free().then(&m1.to_free()));
443        }
444    }
445
446    #[test]
447    fn equivalent_inverse() {
448        let mut rng = Xoshiro256Plus::seed_from_u64(0xca9bd0d289b4700e);
449        let mut nontrivial = 0;
450        for _ in 1..500 {
451            let m = random_possibly_invertible_matrix(&mut rng);
452            let inv = m.inverse_transform();
453            if let Some(inv) = inv {
454                assert_eq!(
455                    Some(inv.to_free()),
456                    m.to_free().inverse(),
457                    "inverse of {m:?}",
458                );
459                nontrivial += 1;
460            }
461        }
462        assert!(nontrivial > 100, "got {nontrivial} inverses");
463    }
464}
all_is_cubes_base/math/matrix.rs

all_is_cubes_base/math/
matrix.rs
