mesh_geometry/
transforms.rs

1//! Affine transforms and Jacobians for mesh-geometry.
2
3use crate::{Float, Point3, Vec3};
4
5/// 3×3 matrix + translation = affine transform in ℝ³.
6#[derive(Debug, Clone, Copy, PartialEq)]
7pub struct Transform3<T: Float> {
8    /// 3x3 matrix for the linear part of the transform
9    pub m: [[T;3];3],
10    /// Translation vector
11    pub t: Vec3<T>,
12}
13
14impl<T: Float> Transform3<T> {
15    /// Identity transform
16    pub fn identity() -> Self {
17        Transform3 {
18            m: [
19                [T::one(), T::zero(), T::zero()],
20                [T::zero(), T::one(), T::zero()],
21                [T::zero(), T::zero(), T::one()],
22            ],
23            t: Vec3::new(T::zero(), T::zero(), T::zero()),
24        }
25    }
26
27    /// Construct from 3×3 matrix and translation
28    pub fn new(m: [[T;3];3], t: Vec3<T>) -> Self {
29        Transform3 { m, t }
30    }
31
32    /// Apply to a point: x ↦ M·x + t
33    pub fn transform_point(&self, p: Point3<T>) -> Point3<T> {
34        let x = self.m[0][0]*p.x + self.m[0][1]*p.y + self.m[0][2]*p.z + self.t.x;
35        let y = self.m[1][0]*p.x + self.m[1][1]*p.y + self.m[1][2]*p.z + self.t.y;
36        let z = self.m[2][0]*p.x + self.m[2][1]*p.y + self.m[2][2]*p.z + self.t.z;
37        Point3::new(x,y,z)
38    }
39
40    /// Apply to a vector (no translation)
41    pub fn transform_vec(&self, v: Vec3<T>) -> Vec3<T> {
42        Vec3::new(
43            self.m[0][0]*v.x + self.m[0][1]*v.y + self.m[0][2]*v.z,
44            self.m[1][0]*v.x + self.m[1][1]*v.y + self.m[1][2]*v.z,
45            self.m[2][0]*v.x + self.m[2][1]*v.y + self.m[2][2]*v.z,
46        )
47    }
48
49    /// Invert (requires det(M) ≠ 0)
50    pub fn inverse(&self) -> Option<Transform3<T>> {
51        // compute adjugate/det of 3×3 matrix
52        let m = &self.m;
53        let det = 
54            m[0][0]*(m[1][1]*m[2][2] - m[1][2]*m[2][1]) -
55            m[0][1]*(m[1][0]*m[2][2] - m[1][2]*m[2][0]) +
56            m[0][2]*(m[1][0]*m[2][1] - m[1][1]*m[2][0]);
57        if det.abs() < T::epsilon() { return None; }
58        let inv_det = T::one()/det;
59        let adj = [
60            [
61                ( m[1][1]*m[2][2] - m[1][2]*m[2][1]) * inv_det,
62               -(m[0][1]*m[2][2] - m[0][2]*m[2][1]) * inv_det,
63                ( m[0][1]*m[1][2] - m[0][2]*m[1][1]) * inv_det,
64            ],
65            [
66               -(m[1][0]*m[2][2] - m[1][2]*m[2][0]) * inv_det,
67                ( m[0][0]*m[2][2] - m[0][2]*m[2][0]) * inv_det,
68               -(m[0][0]*m[1][2] - m[0][2]*m[1][0]) * inv_det,
69            ],
70            [
71                ( m[1][0]*m[2][1] - m[1][1]*m[2][0]) * inv_det,
72               -(m[0][0]*m[2][1] - m[0][1]*m[2][0]) * inv_det,
73                ( m[0][0]*m[1][1] - m[0][1]*m[1][0]) * inv_det,
74            ],
75        ];
76        // Inverse translation: -M⁻¹·t
77        let inv_t = {
78            let t = &self.t;
79            Vec3::new(
80                -(adj[0][0]*t.x + adj[0][1]*t.y + adj[0][2]*t.z),
81                -(adj[1][0]*t.x + adj[1][1]*t.y + adj[1][2]*t.z),
82                -(adj[2][0]*t.x + adj[2][1]*t.y + adj[2][2]*t.z),
83            )
84        };
85        Some(Transform3 { m: adj, t: inv_t })
86    }
87}
88
89#[cfg(test)]
90mod tests {
91    use super::*;
92    use crate::{Point3, Vec3};
93
94    #[test]
95    fn transform_point_vec_roundtrip() {
96        let m = [
97            [2.0, 0.0, 0.0],
98            [0.0, 3.0, 0.0],
99            [0.0, 0.0, 4.0],
100        ];
101        let t = Vec3::new(1.0, -1.0, 2.0);
102        let tf = Transform3::new(m, t);
103        let p = Point3::new(1.0, 1.0, 1.0);
104        let p2 = tf.transform_point(p);
105        assert_eq!(p2, Point3::new(3.0, 2.0, 6.0));
106    }
107}