Skip to main content

dotloom_geometry/
affine.rs

1use serde::{Deserialize, Serialize};
2
3use crate::{GeoResult, GeometryError, Point, Vector};
4
5/// A 2D affine transform `[a b c d e f]` with the SVG convention:
6///
7/// ```text
8/// x' = a·x + c·y + e
9/// y' = b·x + d·y + f
10/// ```
11///
12/// Serialized as a 6-element array.
13#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
14#[serde(from = "[f64; 6]", into = "[f64; 6]")]
15pub struct Affine {
16    /// Coefficients `[a, b, c, d, e, f]`.
17    pub m: [f64; 6],
18}
19
20impl From<[f64; 6]> for Affine {
21    fn from(m: [f64; 6]) -> Self {
22        Self { m }
23    }
24}
25impl From<Affine> for [f64; 6] {
26    fn from(a: Affine) -> Self {
27        a.m
28    }
29}
30
31impl Default for Affine {
32    fn default() -> Self {
33        Self::IDENTITY
34    }
35}
36
37/// Classification of the linear part of a transform.
38#[derive(Debug, Clone, Copy, PartialEq)]
39pub enum LinearKind {
40    /// Rotation + uniform scale `s` (`reflected` when the determinant is negative).
41    Similarity {
42        /// Uniform scale factor (> 0).
43        scale: f64,
44        /// Whether orientation is reversed.
45        reflected: bool,
46    },
47    /// Anything else (non-uniform scale or shear).
48    General,
49    /// Determinant is (near) zero.
50    Singular,
51}
52
53impl Affine {
54    /// Determinants with magnitude below this are treated as singular.
55    pub const SINGULAR_EPS: f64 = 1e-14;
56
57    /// The identity transform.
58    pub const IDENTITY: Self = Self { m: [1.0, 0.0, 0.0, 1.0, 0.0, 0.0] };
59
60    /// Translation by `v`.
61    #[must_use]
62    pub const fn translate(v: Vector) -> Self {
63        Self { m: [1.0, 0.0, 0.0, 1.0, v.x, v.y] }
64    }
65
66    /// Counter-clockwise rotation by `angle` radians about the origin.
67    #[must_use]
68    pub fn rotate(angle: f64) -> Self {
69        let (s, c) = crate::math::sin_cos(angle);
70        Self { m: [c, s, -s, c, 0.0, 0.0] }
71    }
72
73    /// Rotation about `center`.
74    #[must_use]
75    pub fn rotate_about(angle: f64, center: Point) -> Self {
76        Self::translate(-center.to_vector()).then(Self::rotate(angle)).then(Self::translate(center.to_vector()))
77    }
78
79    /// Scale about the origin.
80    #[must_use]
81    pub const fn scale(sx: f64, sy: f64) -> Self {
82        Self { m: [sx, 0.0, 0.0, sy, 0.0, 0.0] }
83    }
84
85    /// Scale about `center`.
86    #[must_use]
87    pub fn scale_about(sx: f64, sy: f64, center: Point) -> Self {
88        Self::translate(-center.to_vector()).then(Self::scale(sx, sy)).then(Self::translate(center.to_vector()))
89    }
90
91    /// Reflection across the line through `p` with direction `dir`.
92    pub fn mirror(p: Point, dir: Vector) -> GeoResult<Self> {
93        let d = dir.normalize().ok_or(GeometryError::Degenerate("mirror axis direction"))?;
94        let (x, y) = (d.x, d.y);
95        let lin = Self { m: [x * x - y * y, 2.0 * x * y, 2.0 * x * y, y * y - x * x, 0.0, 0.0] };
96        Ok(Self::translate(-p.to_vector()).then(lin).then(Self::translate(p.to_vector())))
97    }
98
99    /// Composition: apply `self` first, then `next`.
100    #[must_use]
101    pub fn then(self, next: Self) -> Self {
102        let [a1, b1, c1, d1, e1, f1] = self.m;
103        let [a2, b2, c2, d2, e2, f2] = next.m;
104        Self {
105            m: [
106                a2 * a1 + c2 * b1,
107                b2 * a1 + d2 * b1,
108                a2 * c1 + c2 * d1,
109                b2 * c1 + d2 * d1,
110                a2 * e1 + c2 * f1 + e2,
111                b2 * e1 + d2 * f1 + f2,
112            ],
113        }
114    }
115
116    /// Determinant of the linear part.
117    #[must_use]
118    pub fn determinant(self) -> f64 {
119        self.m[0] * self.m[3] - self.m[1] * self.m[2]
120    }
121
122    /// All coefficients are finite.
123    #[must_use]
124    pub fn is_finite(self) -> bool {
125        self.m.iter().all(|v| v.is_finite())
126    }
127
128    /// Whether this is exactly the identity.
129    #[must_use]
130    pub fn is_identity(self) -> bool {
131        self == Self::IDENTITY
132    }
133
134    /// Inverse transform. Fails for singular or non-finite transforms.
135    pub fn inverse(self) -> GeoResult<Self> {
136        if !self.is_finite() {
137            return Err(GeometryError::NonFinite("transform"));
138        }
139        let det = self.determinant();
140        let [a, b, c, d, e, f] = self.m;
141        // Relative test: |det| compared with the squared column norms, so tiny but
142        // well-conditioned scales stay invertible.
143        let scale = (a * a + b * b).max(c * c + d * d);
144        if det == 0.0 || !det.is_finite() || det.abs() <= Self::SINGULAR_EPS * scale {
145            return Err(GeometryError::SingularTransform { determinant: det });
146        }
147        let inv = 1.0 / det;
148        let na = d * inv;
149        let nb = -b * inv;
150        let nc = -c * inv;
151        let nd = a * inv;
152        Ok(Self { m: [na, nb, nc, nd, -(na * e + nc * f), -(nb * e + nd * f)] })
153    }
154
155    /// Transform a point.
156    #[must_use]
157    pub fn apply(self, p: Point) -> Point {
158        let [a, b, c, d, e, f] = self.m;
159        Point::new(a * p.x + c * p.y + e, b * p.x + d * p.y + f)
160    }
161
162    /// Transform a vector (ignores translation).
163    #[must_use]
164    pub fn apply_vector(self, v: Vector) -> Vector {
165        let [a, b, c, d, _, _] = self.m;
166        Vector::new(a * v.x + c * v.y, b * v.x + d * v.y)
167    }
168
169    /// Translation part.
170    #[must_use]
171    pub fn translation(self) -> Vector {
172        Vector::new(self.m[4], self.m[5])
173    }
174
175    /// Classify the linear part with relative tolerance `rel`.
176    #[must_use]
177    pub fn linear_kind(self, rel: f64) -> LinearKind {
178        let [a, b, c, d, _, _] = self.m;
179        let det = a * d - b * c;
180        let n1 = a * a + b * b;
181        let n2 = c * c + d * d;
182        let scale = n1.max(n2);
183        if !det.is_finite() || det.abs() <= Self::SINGULAR_EPS * scale.max(f64::MIN_POSITIVE) {
184            return LinearKind::Singular;
185        }
186        let ortho = (a * c + b * d).abs();
187        if (n1 - n2).abs() <= rel * scale && ortho <= rel * scale {
188            LinearKind::Similarity { scale: ((n1 + n2) * 0.5).sqrt(), reflected: det < 0.0 }
189        } else {
190            LinearKind::General
191        }
192    }
193
194    /// Rotation angle of the linear part (meaningful for similarities).
195    #[must_use]
196    pub fn rotation_angle(self) -> f64 {
197        crate::math::atan2(self.m[1], self.m[0])
198    }
199
200    /// Whether the linear part is a diagonal (axis aligned) scale with positive or
201    /// negative factors: no rotation or shear.
202    #[must_use]
203    pub fn is_axis_aligned(self, rel: f64) -> bool {
204        let [a, b, c, d, _, _] = self.m;
205        let scale = a.abs().max(d.abs()).max(f64::MIN_POSITIVE);
206        b.abs() <= rel * scale && c.abs() <= rel * scale
207    }
208}
209
210#[cfg(test)]
211mod tests {
212    use super::*;
213    use core::f64::consts::FRAC_PI_2;
214
215    fn close(a: Point, b: Point) -> bool {
216        a.distance(b) < 1e-9
217    }
218
219    #[test]
220    fn compose_order() {
221        let t = Affine::translate(Vector::new(10.0, 0.0)).then(Affine::rotate(FRAC_PI_2));
222        // translate first, then rotate: (1,0) -> (11,0) -> (0,11)
223        assert!(close(t.apply(Point::new(1.0, 0.0)), Point::new(0.0, 11.0)));
224    }
225
226    #[test]
227    fn inverse_roundtrip() {
228        let t = Affine::rotate_about(0.3, Point::new(5.0, -2.0))
229            .then(Affine::scale(2.0, 0.5))
230            .then(Affine::translate(Vector::new(3.0, 4.0)));
231        let inv = t.inverse().unwrap();
232        let p = Point::new(7.25, -1.5);
233        assert!(close(inv.apply(t.apply(p)), p));
234    }
235
236    #[test]
237    fn singular_is_rejected() {
238        let t = Affine::scale(1.0, 0.0);
239        assert!(matches!(t.inverse(), Err(GeometryError::SingularTransform { .. })));
240        let nan = Affine { m: [f64::NAN, 0.0, 0.0, 1.0, 0.0, 0.0] };
241        assert!(matches!(nan.inverse(), Err(GeometryError::NonFinite(_))));
242    }
243
244    #[test]
245    fn classify() {
246        assert!(matches!(
247            Affine::rotate(0.7).then(Affine::scale(3.0, 3.0)).linear_kind(1e-12),
248            LinearKind::Similarity { scale, reflected: false } if (scale - 3.0).abs() < 1e-12
249        ));
250        assert_eq!(Affine::scale(2.0, 1.0).linear_kind(1e-12), LinearKind::General);
251        assert!(matches!(
252            Affine::mirror(Point::ORIGIN, Vector::new(1.0, 1.0)).unwrap().linear_kind(1e-12),
253            LinearKind::Similarity { reflected: true, .. }
254        ));
255        assert_eq!(Affine::scale(0.0, 1.0).linear_kind(1e-12), LinearKind::Singular);
256    }
257
258    #[test]
259    fn mirror_reflects() {
260        let m = Affine::mirror(Point::new(0.0, 1.0), Vector::new(1.0, 0.0)).unwrap();
261        assert!(close(m.apply(Point::new(3.0, 3.0)), Point::new(3.0, -1.0)));
262    }
263}