ogeom_math/
transform.rs

1//! Rigid and similarity transforms, and general affine transforms.
2//!
3//! [`Transform`] is a similarity: an orthonormal linear part, a uniform scale
4//! and a translation. That covers everything a solid modeller applies to a
5//! shape (placement, rotation, mirroring, uniform scaling) while preserving
6//! the two properties the geometry depends on: angles are unchanged, and an
7//! analytic surface stays the same *kind* of analytic surface. A cylinder
8//! remains a cylinder.
9//!
10//! [`GeneralTransform`] drops both guarantees, allowing non-uniform scaling and
11//! shear. It is a separate type on purpose: applying one turns a circle into an
12//! ellipse and a cylinder into something with no analytic form at all, so it
13//! cannot be used interchangeably.
14//!
15//! # Form classification
16//!
17//! Every [`Transform`] carries a [`TransformKind`], and applying one dispatches
18//! on it: a translation adds a vector, the identity does nothing at all. That
19//! matters because transforms are applied to every control point of every
20//! curve, every vertex of every tessellation, over an entire model: the
21//! difference between a branch and nine multiplies, repeated a hundred million
22//! times, is real.
23//!
24//! The kind is *derived from* the data rather than asserted alongside it, so it
25//! cannot drift out of agreement with the matrix it describes. Every
26//! constructor routes through one private classifier. There is no way to build a
27//! transform that claims more structure than it has.
28
29use core::ops::Mul;
30
31use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
32
33use crate::{
34    Axis, Direction, Direction2, Frame, Frame2, Matrix2, Matrix3, Point, Point2, Quaternion,
35    Vector, Vector2,
36};
37
38/// How much structure a [`Transform`] has, for dispatch.
39#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
40pub enum TransformKind {
41    /// Does nothing.
42    #[default]
43    Identity,
44    /// Translation only.
45    Translation,
46    /// Rotation about an axis through the origin, possibly with a translation.
47    Rotation,
48    /// Reflection through a point (equivalently, a scale of `-1`).
49    PointMirror,
50    /// Reflection in a plane.
51    PlaneMirror,
52    /// Uniform scaling, possibly with a translation.
53    Scale,
54    /// Anything else: a combination with no simpler description.
55    Compound,
56}
57
58/// A similarity transform: orthonormal rotation or reflection, uniform scale,
59/// translation.
60///
61/// Applied as `p -> linear * (scale * p) + translation`.
62#[derive(Debug, Clone, Copy, PartialEq)]
63pub struct Transform {
64    linear: Matrix3,
65    scale: f64,
66    translation: Vector,
67    kind: TransformKind,
68}
69
70/// A similarity transform in the plane.
71#[derive(Debug, Clone, Copy, PartialEq)]
72pub struct Transform2 {
73    linear: Matrix2,
74    scale: f64,
75    translation: Vector2,
76    kind: TransformKind,
77}
78
79/// A general affine transform: any linear part, plus a translation.
80///
81/// Non-uniform scaling and shear are allowed, so angles are not preserved and
82/// analytic geometry does not survive intact. Kept distinct from [`Transform`]
83/// so that cannot happen by accident.
84#[derive(Debug, Clone, Copy, PartialEq)]
85pub struct GeneralTransform {
86    /// The linear part.
87    pub linear: Matrix3,
88    /// The translation.
89    pub translation: Vector,
90}
91
92/// How far from `1` a scale factor may sit and still count as unit, and how far
93/// a matrix may stray from a canonical form and still be classified as it.
94/// Dimensionless; classification is a fast-path hint, and a value that just
95/// misses the threshold is merely applied by the general path.
96const CLASSIFY_EPS: f64 = 1e-12;
97
98impl Default for Transform {
99    fn default() -> Self {
100        Self::IDENTITY
101    }
102}
103
104impl Transform {
105    /// The identity.
106    pub const IDENTITY: Self = Self {
107        linear: Matrix3::IDENTITY,
108        scale: 1.0,
109        translation: Vector::ZERO,
110        kind: TransformKind::Identity,
111    };
112
113    /// Derive the kind from the parts, and build the transform.
114    ///
115    /// The single constructor everything else routes through, so a transform's
116    /// kind can never disagree with what it actually does.
117    fn build(linear: Matrix3, scale: f64, translation: Vector) -> Self {
118        let kind = Self::classify(&linear, scale, translation);
119        Self {
120            linear,
121            scale,
122            translation,
123            kind,
124        }
125    }
126
127    /// A transform from the parts it is stored as.
128    ///
129    /// `linear` is the orthonormal part alone and `scale` the uniform factor
130    /// beside it, which is how a [`Transform`] holds them. The kind is
131    /// re-derived rather than taken on trust, since it is a function of the
132    /// other three.
133    ///
134    /// For reading a document back. Going the long way round (multiplying the
135    /// scale into the matrix and asking
136    /// [`GeneralTransform::to_similarity`](crate::GeneralTransform::to_similarity)
137    /// to factor it out again) recovers a transform that is *close*, not the
138    /// one that was written, and a round trip that drifts a little each time is
139    /// not a round trip.
140    ///
141    /// # Errors
142    ///
143    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `linear` is
144    /// not orthonormal within `eps`, or `scale` is not finite and non-zero; a
145    /// placement that squashes space is not a placement.
146    pub fn from_parts(
147        linear: Matrix3,
148        scale: f64,
149        translation: Vector,
150        eps: f64,
151    ) -> OgeomResult<Self> {
152        if !scale.is_finite() || scale == 0.0 {
153            ogeom_bail!(
154                Construction,
155                "a placement's scale must be finite and non-zero; got {scale}"
156            );
157        }
158        if !linear.is_orthonormal(eps) {
159            ogeom_bail!(
160                Construction,
161                "a placement's linear part must be orthonormal; this one shears                  or scales unevenly"
162            );
163        }
164        Ok(Self::build(linear, scale, translation))
165    }
166
167    /// Classify a similarity from its parts.
168    fn classify(linear: &Matrix3, scale: f64, translation: Vector) -> TransformKind {
169        let is_identity_linear = linear.is_equal(&Matrix3::IDENTITY, CLASSIFY_EPS);
170        let unit_scale = (scale - 1.0).abs() <= CLASSIFY_EPS;
171        let negative_unit_scale = (scale + 1.0).abs() <= CLASSIFY_EPS;
172        let no_translation = translation.square_magnitude() == 0.0;
173
174        if is_identity_linear && unit_scale {
175            return if no_translation {
176                TransformKind::Identity
177            } else {
178                TransformKind::Translation
179            };
180        }
181        if is_identity_linear && negative_unit_scale {
182            return TransformKind::PointMirror;
183        }
184        if is_identity_linear {
185            return TransformKind::Scale;
186        }
187        if unit_scale && linear.is_orthonormal(CLASSIFY_EPS) {
188            // Determinant separates a rotation from a reflection. Both are
189            // orthonormal, and confusing them flips the sense of every face.
190            return if linear.determinant() > 0.0 {
191                TransformKind::Rotation
192            } else {
193                TransformKind::PlaneMirror
194            };
195        }
196        TransformKind::Compound
197    }
198
199    /// A translation.
200    #[must_use]
201    pub fn translation(v: Vector) -> Self {
202        Self::build(Matrix3::IDENTITY, 1.0, v)
203    }
204
205    /// A rotation of `angle` radians about `axis`.
206    #[must_use]
207    pub fn rotation(axis: Axis, angle: f64) -> Self {
208        let linear = Matrix3::rotation(axis.direction, angle);
209        // Rotating about an axis that misses the origin: move the axis point to
210        // the origin, rotate, move back. Written as one translation so the
211        // result stays a single transform.
212        let p = axis.location.to_vector();
213        Self::build(linear, 1.0, p - linear * p)
214    }
215
216    /// A rotation given as a quaternion, about an axis through the origin.
217    #[must_use]
218    pub fn from_quaternion(q: Quaternion) -> Self {
219        Self::build(q.to_matrix(), 1.0, Vector::ZERO)
220    }
221
222    /// A uniform scaling about `centre`.
223    ///
224    /// # Errors
225    ///
226    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `factor` is
227    /// zero or non-finite. A zero scale collapses every shape to a point and is
228    /// not invertible, so it is refused rather than allowed to produce
229    /// degenerate geometry.
230    pub fn scaling(centre: Point, factor: f64, tol: Tolerances) -> OgeomResult<Self> {
231        if !factor.is_finite() || factor.abs() <= tol.confusion() {
232            ogeom_bail!(Construction, "scale factor {factor} is degenerate");
233        }
234        let c = centre.to_vector();
235        Ok(Self::build(Matrix3::IDENTITY, factor, c - c * factor))
236    }
237
238    /// Reflection through a point.
239    #[must_use]
240    pub fn point_mirror(centre: Point) -> Self {
241        let c = centre.to_vector();
242        Self::build(Matrix3::IDENTITY, -1.0, c + c)
243    }
244
245    /// Reflection in the plane through `origin` with the given `normal`.
246    #[must_use]
247    pub fn plane_mirror(origin: Point, normal: Direction) -> Self {
248        let linear = Matrix3::reflection(normal);
249        let p = origin.to_vector();
250        Self::build(linear, 1.0, p - linear * p)
251    }
252
253    /// Reflection in a line: a half turn about it.
254    #[must_use]
255    pub fn axis_mirror(axis: Axis) -> Self {
256        Self::rotation(axis, core::f64::consts::PI)
257    }
258
259    /// The transform taking world coordinates into `frame`'s local coordinates.
260    #[must_use]
261    pub fn to_frame(frame: &Frame) -> Self {
262        let linear = frame.to_matrix().transposed();
263        Self::build(linear, 1.0, -(linear * frame.origin().to_vector()))
264    }
265
266    /// The transform taking `frame`'s local coordinates into world coordinates.
267    #[must_use]
268    pub fn from_frame(frame: &Frame) -> Self {
269        Self::build(frame.to_matrix(), 1.0, frame.origin().to_vector())
270    }
271
272    /// The transform taking `from`'s local coordinates into `to`'s.
273    #[must_use]
274    pub fn between_frames(from: &Frame, to: &Frame) -> Self {
275        Self::to_frame(to) * Self::from_frame(from)
276    }
277
278    /// This transform's classification.
279    #[must_use]
280    pub const fn kind(&self) -> TransformKind {
281        self.kind
282    }
283
284    /// The orthonormal part.
285    #[must_use]
286    pub const fn linear(&self) -> Matrix3 {
287        self.linear
288    }
289
290    /// The uniform scale factor. Negative for a point mirror.
291    #[must_use]
292    pub const fn scale_factor(&self) -> f64 {
293        self.scale
294    }
295
296    /// The translation.
297    #[must_use]
298    pub const fn translation_vector(&self) -> Vector {
299        self.translation
300    }
301
302    /// Whether this transform preserves handedness.
303    ///
304    /// A shape transformed by a transform that does not must have its
305    /// orientation flipped to stay consistent. Otherwise a mirrored solid ends
306    /// up inside out.
307    #[must_use]
308    pub fn preserves_handedness(&self) -> bool {
309        self.linear.determinant() * self.scale.signum() > 0.0
310    }
311
312    /// Apply to a point.
313    #[must_use]
314    pub fn apply(&self, p: Point) -> Point {
315        match self.kind {
316            TransformKind::Identity => p,
317            TransformKind::Translation => p + self.translation,
318            TransformKind::PointMirror | TransformKind::Scale => {
319                Point::from_vector(p.to_vector() * self.scale + self.translation)
320            }
321            TransformKind::Rotation | TransformKind::PlaneMirror => {
322                Point::from_vector(self.linear * p.to_vector() + self.translation)
323            }
324            TransformKind::Compound => {
325                Point::from_vector(self.linear * (p.to_vector() * self.scale) + self.translation)
326            }
327        }
328    }
329
330    /// Apply to a free vector. Translation does not affect it.
331    #[must_use]
332    pub fn apply_vector(&self, v: Vector) -> Vector {
333        match self.kind {
334            TransformKind::Identity | TransformKind::Translation => v,
335            TransformKind::PointMirror | TransformKind::Scale => v * self.scale,
336            TransformKind::Rotation | TransformKind::PlaneMirror => self.linear * v,
337            TransformKind::Compound => self.linear * (v * self.scale),
338        }
339    }
340
341    /// Apply to a direction, renormalizing.
342    ///
343    /// A similarity maps unit vectors to vectors of length `|scale|`, so the
344    /// result is rescaled. Under a negative scale the direction reverses, which
345    /// is the correct behaviour and not a sign error.
346    ///
347    /// # Errors
348    ///
349    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the result
350    /// cannot be normalized, which a valid similarity never produces.
351    pub fn apply_direction(&self, d: Direction, tol: Tolerances) -> OgeomResult<Direction> {
352        match self.kind {
353            TransformKind::Identity | TransformKind::Translation => Ok(d),
354            // A negative scale is a point mirror scaled: it reverses.
355            TransformKind::Scale if self.scale > 0.0 => Ok(d),
356            TransformKind::PointMirror | TransformKind::Scale => Ok(d.reversed()),
357            _ => Direction::new(self.apply_vector(d.vector()), tol),
358        }
359    }
360
361    /// Apply to a frame, transforming origin and all three axes.
362    ///
363    /// # Errors
364    ///
365    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
366    /// transformed axes cannot be renormalized.
367    pub fn apply_frame(&self, f: &Frame, tol: Tolerances) -> OgeomResult<Frame> {
368        Frame::from_axes(
369            self.apply(f.origin()),
370            self.apply_direction(f.x(), tol)?,
371            self.apply_direction(f.y(), tol)?,
372            self.apply_direction(f.z(), tol)?,
373            tol,
374        )
375    }
376
377    /// The inverse.
378    ///
379    /// # Errors
380    ///
381    /// [`OgeomError::Numeric`](ogeom_core::OgeomError::Numeric) if the linear part is
382    /// singular, which a valid similarity never is.
383    pub fn inverse(&self) -> OgeomResult<Self> {
384        if self.kind == TransformKind::Identity {
385            return Ok(Self::IDENTITY);
386        }
387        if self.kind == TransformKind::Translation {
388            return Ok(Self::translation(-self.translation));
389        }
390        if self.scale == 0.0 {
391            ogeom_bail!(Numeric, "transform has a zero scale and no inverse");
392        }
393        // The linear part is orthonormal, so its inverse is its transpose. No
394        // need to go through a general inversion, and no rounding beyond the
395        // transpose itself.
396        let inv_linear = self.linear.transposed();
397        let inv_scale = 1.0 / self.scale;
398        Ok(Self::build(
399            inv_linear,
400            inv_scale,
401            -(inv_linear * self.translation) * inv_scale,
402        ))
403    }
404
405    /// Whether two transforms agree in effect.
406    #[must_use]
407    pub fn is_equal(&self, other: &Self, tol: Tolerances) -> bool {
408        (self.scale - other.scale).abs() <= CLASSIFY_EPS
409            && self.linear.is_equal(&other.linear, CLASSIFY_EPS)
410            && self.translation.is_equal(other.translation, tol)
411    }
412
413    /// This transform as a general affine one.
414    #[must_use]
415    pub fn to_general(&self) -> GeneralTransform {
416        GeneralTransform {
417            linear: self.linear * self.scale,
418            translation: self.translation,
419        }
420    }
421}
422
423impl Mul for Transform {
424    type Output = Self;
425    /// Composition. `(a * b)` applies `b` first, then `a`.
426    fn mul(self, b: Self) -> Self {
427        if self.kind == TransformKind::Identity {
428            return b;
429        }
430        if b.kind == TransformKind::Identity {
431            return self;
432        }
433        Self::build(
434            self.linear * b.linear,
435            self.scale * b.scale,
436            self.linear * (b.translation * self.scale) + self.translation,
437        )
438    }
439}
440
441impl Default for Transform2 {
442    fn default() -> Self {
443        Self::IDENTITY
444    }
445}
446
447impl Transform2 {
448    /// The identity.
449    pub const IDENTITY: Self = Self {
450        linear: Matrix2::IDENTITY,
451        scale: 1.0,
452        translation: Vector2::ZERO,
453        kind: TransformKind::Identity,
454    };
455
456    fn build(linear: Matrix2, scale: f64, translation: Vector2) -> Self {
457        let is_identity_linear = linear.is_equal(&Matrix2::IDENTITY, CLASSIFY_EPS);
458        let unit_scale = (scale - 1.0).abs() <= CLASSIFY_EPS;
459        let kind = if is_identity_linear && unit_scale {
460            if translation.square_magnitude() == 0.0 {
461                TransformKind::Identity
462            } else {
463                TransformKind::Translation
464            }
465        } else if is_identity_linear && (scale + 1.0).abs() <= CLASSIFY_EPS {
466            TransformKind::PointMirror
467        } else if is_identity_linear {
468            TransformKind::Scale
469        } else if unit_scale && (linear.determinant().abs() - 1.0).abs() <= CLASSIFY_EPS {
470            if linear.determinant() > 0.0 {
471                TransformKind::Rotation
472            } else {
473                TransformKind::PlaneMirror
474            }
475        } else {
476            TransformKind::Compound
477        };
478        Self {
479            linear,
480            scale,
481            translation,
482            kind,
483        }
484    }
485
486    /// A translation.
487    #[must_use]
488    pub fn translation(v: Vector2) -> Self {
489        Self::build(Matrix2::IDENTITY, 1.0, v)
490    }
491
492    /// A rotation about `centre`.
493    #[must_use]
494    pub fn rotation(centre: Point2, angle: f64) -> Self {
495        let linear = Matrix2::rotation(angle);
496        let c = centre.to_vector();
497        Self::build(linear, 1.0, c - linear * c)
498    }
499
500    /// A uniform scaling about `centre`.
501    ///
502    /// # Errors
503    ///
504    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `factor` is
505    /// degenerate.
506    pub fn scaling(centre: Point2, factor: f64, tol: Tolerances) -> OgeomResult<Self> {
507        if !factor.is_finite() || factor.abs() <= tol.confusion() {
508            ogeom_bail!(Construction, "scale factor {factor} is degenerate");
509        }
510        let c = centre.to_vector();
511        Ok(Self::build(Matrix2::IDENTITY, factor, c - c * factor))
512    }
513
514    /// Reflection in the line through `origin` with the given `normal`.
515    #[must_use]
516    pub fn line_mirror(origin: Point2, normal: Direction2) -> Self {
517        let (x, y) = (normal.x(), normal.y());
518        let linear = Matrix2::new([
519            [(-2.0f64).mul_add(x * x, 1.0), -2.0 * x * y],
520            [-2.0 * x * y, (-2.0f64).mul_add(y * y, 1.0)],
521        ]);
522        let p = origin.to_vector();
523        Self::build(linear, 1.0, p - linear * p)
524    }
525
526    /// The transform taking `frame`'s local coordinates into world
527    /// coordinates; a reflection where the frame is left-handed.
528    #[must_use]
529    pub fn from_frame(frame: &Frame2) -> Self {
530        let (x, y) = (frame.x(), frame.y());
531        let linear = Matrix2::new([[x.x(), y.x()], [x.y(), y.y()]]);
532        Self::build(linear, 1.0, frame.origin().to_vector())
533    }
534
535    /// This transform's classification.
536    #[must_use]
537    pub const fn kind(&self) -> TransformKind {
538        self.kind
539    }
540
541    /// The orthonormal part.
542    #[must_use]
543    pub const fn linear(&self) -> Matrix2 {
544        self.linear
545    }
546
547    /// The uniform scale factor. Negative for a point mirror.
548    #[must_use]
549    pub const fn scale_factor(&self) -> f64 {
550        self.scale
551    }
552
553    /// The translation.
554    #[must_use]
555    pub const fn translation_vector(&self) -> Vector2 {
556        self.translation
557    }
558
559    /// Whether this transform preserves handedness. In the plane a
560    /// negative scale is a half turn, which does: only the linear part's
561    /// reflection decides.
562    #[must_use]
563    pub fn preserves_handedness(&self) -> bool {
564        self.linear.determinant() > 0.0
565    }
566
567    /// Apply to a point.
568    #[must_use]
569    pub fn apply(&self, p: Point2) -> Point2 {
570        match self.kind {
571            TransformKind::Identity => p,
572            TransformKind::Translation => p + self.translation,
573            TransformKind::PointMirror | TransformKind::Scale => {
574                Point2::from_vector(p.to_vector() * self.scale + self.translation)
575            }
576            TransformKind::Rotation | TransformKind::PlaneMirror => {
577                Point2::from_vector(self.linear * p.to_vector() + self.translation)
578            }
579            TransformKind::Compound => {
580                Point2::from_vector(self.linear * (p.to_vector() * self.scale) + self.translation)
581            }
582        }
583    }
584
585    /// Apply to a direction, renormalizing.
586    ///
587    /// # Errors
588    ///
589    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the result
590    /// cannot be normalized, which a valid similarity never produces.
591    pub fn apply_direction(&self, d: Direction2, tol: Tolerances) -> OgeomResult<Direction2> {
592        match self.kind {
593            TransformKind::Identity | TransformKind::Translation => Ok(d),
594            TransformKind::Scale if self.scale > 0.0 => Ok(d),
595            TransformKind::PointMirror | TransformKind::Scale => Ok(d.reversed()),
596            _ => Direction2::new(self.apply_vector(d.vector()), tol),
597        }
598    }
599
600    /// Apply to a frame, transforming the origin and both axes.
601    ///
602    /// # Errors
603    ///
604    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
605    /// transformed axes cannot be renormalized or are no longer perpendicular.
606    pub fn apply_frame(&self, f: &Frame2, tol: Tolerances) -> OgeomResult<Frame2> {
607        Frame2::from_axes(
608            self.apply(f.origin()),
609            self.apply_direction(f.x(), tol)?,
610            self.apply_direction(f.y(), tol)?,
611            tol,
612        )
613    }
614
615    /// Apply to a free vector.
616    #[must_use]
617    pub fn apply_vector(&self, v: Vector2) -> Vector2 {
618        match self.kind {
619            TransformKind::Identity | TransformKind::Translation => v,
620            TransformKind::PointMirror | TransformKind::Scale => v * self.scale,
621            TransformKind::Rotation | TransformKind::PlaneMirror => self.linear * v,
622            TransformKind::Compound => self.linear * (v * self.scale),
623        }
624    }
625
626    /// The inverse.
627    ///
628    /// # Errors
629    ///
630    /// [`OgeomError::Numeric`](ogeom_core::OgeomError::Numeric) if the transform is
631    /// degenerate.
632    pub fn inverse(&self) -> OgeomResult<Self> {
633        if self.kind == TransformKind::Identity {
634            return Ok(Self::IDENTITY);
635        }
636        if self.scale == 0.0 {
637            ogeom_bail!(Numeric, "transform has a zero scale and no inverse");
638        }
639        let inv_linear = self.linear.transposed();
640        let inv_scale = 1.0 / self.scale;
641        Ok(Self::build(
642            inv_linear,
643            inv_scale,
644            -(inv_linear * self.translation) * inv_scale,
645        ))
646    }
647
648    /// Whether two transforms agree in effect.
649    #[must_use]
650    pub fn is_equal(&self, other: &Self, tol: Tolerances) -> bool {
651        (self.scale - other.scale).abs() <= CLASSIFY_EPS
652            && self.linear.is_equal(&other.linear, CLASSIFY_EPS)
653            && self.translation.is_equal(other.translation, tol)
654    }
655}
656
657impl Mul for Transform2 {
658    type Output = Self;
659    fn mul(self, b: Self) -> Self {
660        if self.kind == TransformKind::Identity {
661            return b;
662        }
663        if b.kind == TransformKind::Identity {
664            return self;
665        }
666        Self::build(
667            self.linear * b.linear,
668            self.scale * b.scale,
669            self.linear * (b.translation * self.scale) + self.translation,
670        )
671    }
672}
673
674impl Default for GeneralTransform {
675    fn default() -> Self {
676        Self::IDENTITY
677    }
678}
679
680impl GeneralTransform {
681    /// The identity.
682    pub const IDENTITY: Self = Self {
683        linear: Matrix3::IDENTITY,
684        translation: Vector::ZERO,
685    };
686
687    /// From a linear part and a translation.
688    #[must_use]
689    pub const fn new(linear: Matrix3, translation: Vector) -> Self {
690        Self {
691            linear,
692            translation,
693        }
694    }
695
696    /// Non-uniform scaling about the origin.
697    #[must_use]
698    pub const fn scaling_xyz(x: f64, y: f64, z: f64) -> Self {
699        Self::new(Matrix3::scaling_xyz(x, y, z), Vector::ZERO)
700    }
701
702    /// Apply to a point.
703    #[must_use]
704    pub fn apply(&self, p: Point) -> Point {
705        Point::from_vector(self.linear * p.to_vector() + self.translation)
706    }
707
708    /// Apply to a free vector.
709    #[must_use]
710    pub fn apply_vector(&self, v: Vector) -> Vector {
711        self.linear * v
712    }
713
714    /// Apply to a normal vector.
715    ///
716    /// Normals transform by the inverse transpose, not by the linear part
717    /// itself. Using the linear part directly is only correct for a similarity;
718    /// under any shear or non-uniform scale it tilts normals off the surface
719    /// they belong to, which then breaks every orientation test downstream.
720    ///
721    /// The result is not renormalized. It is a direction, not a length.
722    ///
723    /// # Errors
724    ///
725    /// [`OgeomError::Numeric`](ogeom_core::OgeomError::Numeric) if the linear part is
726    /// singular.
727    pub fn apply_normal(&self, n: Vector) -> OgeomResult<Vector> {
728        Ok(self.linear.inverse()?.transposed() * n)
729    }
730
731    /// Whether this transform preserves handedness.
732    #[must_use]
733    pub fn preserves_handedness(&self) -> bool {
734        self.linear.determinant() > 0.0
735    }
736
737    /// The factor by which volumes are multiplied. Negative if handedness flips.
738    #[must_use]
739    pub fn volume_ratio(&self) -> f64 {
740        self.linear.determinant()
741    }
742
743    /// The inverse.
744    ///
745    /// # Errors
746    ///
747    /// [`OgeomError::Numeric`](ogeom_core::OgeomError::Numeric) if the linear part is
748    /// singular.
749    pub fn inverse(&self) -> OgeomResult<Self> {
750        let inv = self.linear.inverse()?;
751        Ok(Self::new(inv, -(inv * self.translation)))
752    }
753
754    /// Whether this is a similarity, and so can be narrowed to a [`Transform`].
755    #[must_use]
756    pub fn is_similarity(&self, eps: f64) -> bool {
757        self.to_similarity(eps).is_some()
758    }
759
760    /// This transform as a [`Transform`], if it is in fact a similarity.
761    ///
762    /// Returns `None` when the linear part contains shear or non-uniform
763    /// scaling, since no similarity describes it.
764    #[must_use]
765    pub fn to_similarity(&self, eps: f64) -> Option<Transform> {
766        // A similarity's linear part is `s * R` with `R` orthonormal, so its
767        // columns are mutually orthogonal and all of length |s|.
768        let det = self.linear.determinant();
769        if det == 0.0 {
770            return None;
771        }
772        let scale = det.abs().cbrt() * det.signum();
773        let rotation = self.linear * (1.0 / scale);
774        if !rotation.is_orthonormal(eps) {
775            return None;
776        }
777        Some(Transform::build(rotation, scale, self.translation))
778    }
779}
780
781impl Mul for GeneralTransform {
782    type Output = Self;
783    /// Composition. `(a * b)` applies `b` first, then `a`.
784    fn mul(self, b: Self) -> Self {
785        Self::new(
786            self.linear * b.linear,
787            self.linear * b.translation + self.translation,
788        )
789    }
790}
791
792impl From<Transform> for GeneralTransform {
793    fn from(t: Transform) -> Self {
794        t.to_general()
795    }
796}
797
798#[cfg(test)]
799#[allow(clippy::unwrap_used)]
800mod tests {
801    use super::*;
802    use approx::assert_relative_eq;
803
804    const T: Tolerances = Tolerances::millimetres();
805
806    /// A negative scale is a point mirror with a scale: a direction comes
807    /// out as the image of a vector along it does, reversed, in space and
808    /// in the plane.
809    #[test]
810    fn a_negative_scale_reverses_directions() {
811        let s = Transform::scaling(Point::ORIGIN, -2.0, T).unwrap();
812        let d = s.apply_direction(Direction::X, T).unwrap();
813        let v = s.apply_vector(Vector::new(1.0, 0.0, 0.0));
814        assert!(d.vector().dot(v) > 0.0, "{d:?} against {v:?}");
815        assert!(!s.preserves_handedness());
816        let s2 = Transform2::scaling(Point2::ORIGIN, -3.0, T).unwrap();
817        let d2 = s2.apply_direction(Direction2::X, T).unwrap();
818        assert!(d2.vector().dot(s2.apply_vector(Direction2::X.vector())) > 0.0);
819        // A positive scale keeps them.
820        let p = Transform::scaling(Point::ORIGIN, 2.0, T).unwrap();
821        assert_eq!(p.apply_direction(Direction::X, T).unwrap(), Direction::X);
822    }
823
824    fn sample_points() -> [Point; 4] {
825        [
826            Point::ORIGIN,
827            Point::new(1.0, 0.0, 0.0),
828            Point::new(-3.0, 7.5, 2.25),
829            Point::new(1e3, -1e3, 0.5),
830        ]
831    }
832
833    #[test]
834    fn classification_matches_what_the_transform_does() {
835        assert_eq!(Transform::IDENTITY.kind(), TransformKind::Identity);
836        assert_eq!(
837            Transform::translation(Vector::X).kind(),
838            TransformKind::Translation
839        );
840        assert_eq!(
841            Transform::rotation(Axis::Z, 0.5).kind(),
842            TransformKind::Rotation
843        );
844        assert_eq!(
845            Transform::point_mirror(Point::ORIGIN).kind(),
846            TransformKind::PointMirror
847        );
848        assert_eq!(
849            Transform::plane_mirror(Point::ORIGIN, Direction::Z).kind(),
850            TransformKind::PlaneMirror
851        );
852        assert_eq!(
853            Transform::scaling(Point::ORIGIN, 3.0, T).unwrap().kind(),
854            TransformKind::Scale
855        );
856        let compound =
857            Transform::rotation(Axis::Z, 0.5) * Transform::scaling(Point::ORIGIN, 3.0, T).unwrap();
858        assert_eq!(compound.kind(), TransformKind::Compound);
859    }
860
861    #[test]
862    fn a_zero_rotation_classifies_as_identity_not_rotation() {
863        // The classification is derived from the data, so it cannot claim more
864        // structure than the transform has, or less.
865        assert_eq!(
866            Transform::rotation(Axis::Z, 0.0).kind(),
867            TransformKind::Identity
868        );
869        assert_eq!(
870            Transform::translation(Vector::ZERO).kind(),
871            TransformKind::Identity
872        );
873        assert_eq!(
874            Transform::scaling(Point::ORIGIN, 1.0, T).unwrap().kind(),
875            TransformKind::Identity
876        );
877    }
878
879    #[test]
880    fn every_dispatch_path_gives_the_same_answer_as_the_general_one() {
881        // Classification exists for speed, so each fast path must agree
882        // exactly with the general formula it replaces.
883        let cases = [
884            Transform::IDENTITY,
885            Transform::translation(Vector::new(1.0, -2.0, 3.0)),
886            Transform::rotation(Axis::new(Point::new(1.0, 0.0, 0.0), Direction::Z), 0.7),
887            Transform::point_mirror(Point::new(2.0, 0.0, -1.0)),
888            Transform::plane_mirror(Point::new(0.0, 1.0, 0.0), Direction::Y),
889            Transform::scaling(Point::new(1.0, 1.0, 1.0), 2.5, T).unwrap(),
890        ];
891        for t in cases {
892            for p in sample_points() {
893                let general = Point::from_vector(
894                    t.linear() * (p.to_vector() * t.scale_factor()) + t.translation_vector(),
895                );
896                assert!(
897                    t.apply(p).is_equal(general, T),
898                    "fast path for {:?} disagrees",
899                    t.kind()
900                );
901            }
902        }
903    }
904
905    #[test]
906    fn rotation_about_an_off_origin_axis_leaves_the_axis_fixed() {
907        let axis = Axis::new(Point::new(5.0, 3.0, 0.0), Direction::Z);
908        let t = Transform::rotation(axis, 1.234);
909        assert!(t.apply(axis.location).is_equal(axis.location, T));
910        assert!(
911            t.apply(axis.point_at(10.0))
912                .is_equal(axis.point_at(10.0), T)
913        );
914        // A point off the axis moves, staying at the same radius.
915        let p = Point::new(6.0, 3.0, 0.0);
916        assert_relative_eq!(axis.distance_to(t.apply(p)), 1.0, epsilon = 1e-14);
917    }
918
919    #[test]
920    fn scaling_about_a_centre_leaves_the_centre_fixed() {
921        let c = Point::new(3.0, -1.0, 2.0);
922        let t = Transform::scaling(c, 4.0, T).unwrap();
923        assert!(t.apply(c).is_equal(c, T));
924        let p = c + Vector::new(1.0, 0.0, 0.0);
925        assert!(t.apply(p).is_equal(c + Vector::new(4.0, 0.0, 0.0), T));
926    }
927
928    #[test]
929    fn degenerate_scales_are_refused() {
930        assert!(Transform::scaling(Point::ORIGIN, 0.0, T).is_err());
931        assert!(Transform::scaling(Point::ORIGIN, f64::NAN, T).is_err());
932        assert!(Transform::scaling(Point::ORIGIN, f64::INFINITY, T).is_err());
933        assert!(
934            Transform::scaling(Point::ORIGIN, -2.0, T).is_ok(),
935            "negative is fine"
936        );
937    }
938
939    #[test]
940    fn handedness_tracks_mirroring() {
941        assert!(Transform::rotation(Axis::Z, 1.0).preserves_handedness());
942        assert!(Transform::translation(Vector::X).preserves_handedness());
943        assert!(
944            Transform::scaling(Point::ORIGIN, 3.0, T)
945                .unwrap()
946                .preserves_handedness()
947        );
948        assert!(!Transform::plane_mirror(Point::ORIGIN, Direction::Z).preserves_handedness());
949        assert!(!Transform::point_mirror(Point::ORIGIN).preserves_handedness());
950        // Two mirrors make a rotation.
951        let twice = Transform::plane_mirror(Point::ORIGIN, Direction::Z)
952            * Transform::plane_mirror(Point::ORIGIN, Direction::X);
953        assert!(twice.preserves_handedness());
954    }
955
956    #[test]
957    fn axis_mirror_is_a_half_turn() {
958        let t = Transform::axis_mirror(Axis::Z);
959        assert!(
960            t.apply(Point::new(1.0, 0.0, 5.0))
961                .is_equal(Point::new(-1.0, 0.0, 5.0), T)
962        );
963        assert!(t.preserves_handedness(), "a half turn is a rotation");
964    }
965
966    #[test]
967    fn inverse_round_trips_for_every_kind() {
968        let cases = [
969            Transform::IDENTITY,
970            Transform::translation(Vector::new(1.0, -2.0, 3.0)),
971            Transform::rotation(Axis::new(Point::new(1.0, 2.0, 3.0), Direction::Y), 2.1),
972            Transform::point_mirror(Point::new(1.0, 1.0, 1.0)),
973            Transform::plane_mirror(Point::new(0.0, 0.0, 4.0), Direction::Z),
974            Transform::scaling(Point::new(-1.0, 0.0, 0.0), 0.25, T).unwrap(),
975        ];
976        for t in cases {
977            let inv = t.inverse().unwrap();
978            for p in sample_points() {
979                assert!(inv.apply(t.apply(p)).is_equal(p, T), "{:?}", t.kind());
980                assert!(t.apply(inv.apply(p)).is_equal(p, T), "{:?}", t.kind());
981            }
982        }
983    }
984
985    #[test]
986    fn composition_applies_right_to_left() {
987        let a = Transform::translation(Vector::new(10.0, 0.0, 0.0));
988        let b = Transform::rotation(Axis::Z, core::f64::consts::FRAC_PI_2);
989        let p = Point::new(1.0, 0.0, 0.0);
990        assert!((a * b).apply(p).is_equal(a.apply(b.apply(p)), T));
991        assert!((b * a).apply(p).is_equal(b.apply(a.apply(p)), T));
992        // And the two orders genuinely differ.
993        assert!(!(a * b).is_equal(&(b * a), T));
994    }
995
996    #[test]
997    fn composition_is_associative() {
998        let a = Transform::rotation(Axis::X, 0.3);
999        let b = Transform::scaling(Point::new(1.0, 0.0, 0.0), 2.0, T).unwrap();
1000        let c = Transform::translation(Vector::new(0.0, 5.0, 0.0));
1001        assert!(((a * b) * c).is_equal(&(a * (b * c)), T));
1002    }
1003
1004    #[test]
1005    fn vectors_ignore_translation_and_directions_stay_unit() {
1006        let t = Transform::translation(Vector::new(100.0, 0.0, 0.0))
1007            * Transform::rotation(Axis::Z, 0.9);
1008        let v = Vector::new(1.0, 2.0, 3.0);
1009        assert!(
1010            t.apply_vector(v)
1011                .is_equal(Transform::rotation(Axis::Z, 0.9).apply_vector(v), T)
1012        );
1013        let d = t.apply_direction(Direction::X, T).unwrap();
1014        assert_relative_eq!(d.vector().magnitude(), 1.0, epsilon = 1e-15);
1015    }
1016
1017    #[test]
1018    fn a_point_mirror_reverses_directions() {
1019        let t = Transform::point_mirror(Point::new(5.0, 5.0, 5.0));
1020        assert!(
1021            t.apply_direction(Direction::X, T)
1022                .unwrap()
1023                .is_equal(-Direction::X, T)
1024        );
1025        // A positive scale does not.
1026        let s = Transform::scaling(Point::ORIGIN, 3.0, T).unwrap();
1027        assert!(
1028            s.apply_direction(Direction::X, T)
1029                .unwrap()
1030                .is_equal(Direction::X, T)
1031        );
1032    }
1033
1034    #[test]
1035    fn frame_transforms_round_trip_through_world() {
1036        let f = Frame::new(
1037            Point::new(1.0, 2.0, 3.0),
1038            Direction::from_coords(1.0, 1.0, 0.0, T).unwrap(),
1039            Direction::Z,
1040            T,
1041        )
1042        .unwrap();
1043        let to = Transform::to_frame(&f);
1044        let from = Transform::from_frame(&f);
1045        for p in sample_points() {
1046            assert!(from.apply(to.apply(p)).is_equal(p, T));
1047            // And they agree with the frame's own conversion.
1048            assert!(to.apply(p).is_equal(f.to_local(p), T));
1049            assert!(from.apply(f.to_local(p)).is_equal(p, T));
1050        }
1051    }
1052
1053    #[test]
1054    fn between_frames_composes_correctly() {
1055        let a = Frame::new(Point::new(1.0, 0.0, 0.0), Direction::Z, Direction::X, T).unwrap();
1056        let b = Frame::new(Point::new(0.0, 5.0, 0.0), Direction::X, Direction::Y, T).unwrap();
1057        let t = Transform::between_frames(&a, &b);
1058        // A point at local (1,2,3) in `a` must land at the same world position
1059        // when read back out of `b`.
1060        let local = Point::new(1.0, 2.0, 3.0);
1061        assert!(b.to_world(t.apply(local)).is_equal(a.to_world(local), T));
1062    }
1063
1064    #[test]
1065    fn general_transform_normals_use_the_inverse_transpose() {
1066        // Non-uniform scaling: the plane z = x has normal (1, 0, -1) up to
1067        // scale. Scale x by 2 and the plane becomes z = x/2, whose normal is
1068        // (1, 0, -2) up to scale, not (2, 0, -1), which is what applying the
1069        // linear part directly would give.
1070        let g = GeneralTransform::scaling_xyz(2.0, 1.0, 1.0);
1071        let n = Vector::new(1.0, 0.0, -1.0);
1072        let transformed = g.apply_normal(n).unwrap();
1073
1074        let on_plane = Vector::new(1.0, 0.0, 1.0);
1075        assert_relative_eq!(n.dot(on_plane), 0.0, epsilon = 1e-15);
1076        assert_relative_eq!(
1077            transformed.dot(g.apply_vector(on_plane)),
1078            0.0,
1079            epsilon = 1e-14,
1080            max_relative = 1e-14
1081        );
1082        // The naive answer does not stay perpendicular.
1083        assert!(g.apply_vector(n).dot(g.apply_vector(on_plane)).abs() > 1e-6);
1084    }
1085
1086    #[test]
1087    fn general_transform_recognizes_similarities() {
1088        let similar: GeneralTransform = Transform::rotation(Axis::Z, 0.4).into();
1089        assert!(similar.is_similarity(1e-12));
1090        let narrowed = similar.to_similarity(1e-12).unwrap();
1091        assert_eq!(narrowed.kind(), TransformKind::Rotation);
1092
1093        let scaled: GeneralTransform = Transform::scaling(Point::ORIGIN, 3.0, T).unwrap().into();
1094        assert!(scaled.is_similarity(1e-12));
1095        assert_relative_eq!(
1096            scaled.to_similarity(1e-12).unwrap().scale_factor(),
1097            3.0,
1098            epsilon = 1e-12
1099        );
1100
1101        assert!(!GeneralTransform::scaling_xyz(1.0, 2.0, 3.0).is_similarity(1e-12));
1102        let shear = GeneralTransform::new(
1103            Matrix3::new([[1.0, 0.5, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]),
1104            Vector::ZERO,
1105        );
1106        assert!(!shear.is_similarity(1e-12));
1107    }
1108
1109    #[test]
1110    fn general_transform_volume_ratio_and_inverse() {
1111        let g = GeneralTransform::scaling_xyz(2.0, 3.0, 4.0);
1112        assert_relative_eq!(g.volume_ratio(), 24.0);
1113        assert!(g.preserves_handedness());
1114        let inv = g.inverse().unwrap();
1115        for p in sample_points() {
1116            assert!(inv.apply(g.apply(p)).is_equal(p, T));
1117        }
1118        let flip = GeneralTransform::scaling_xyz(-1.0, 1.0, 1.0);
1119        assert!(!flip.preserves_handedness());
1120        assert!(
1121            GeneralTransform::scaling_xyz(0.0, 1.0, 1.0)
1122                .inverse()
1123                .is_err()
1124        );
1125    }
1126
1127    #[test]
1128    fn transform2_behaves_like_its_3d_counterpart() {
1129        let r = Transform2::rotation(Point2::new(1.0, 1.0), core::f64::consts::FRAC_PI_2);
1130        assert_eq!(r.kind(), TransformKind::Rotation);
1131        assert!(
1132            r.apply(Point2::new(1.0, 1.0))
1133                .is_equal(Point2::new(1.0, 1.0), T)
1134        );
1135        assert!(
1136            r.apply(Point2::new(2.0, 1.0))
1137                .is_equal(Point2::new(1.0, 2.0), T)
1138        );
1139        assert!(
1140            r.inverse()
1141                .unwrap()
1142                .apply(r.apply(Point2::ORIGIN))
1143                .is_equal(Point2::ORIGIN, T)
1144        );
1145
1146        let m = Transform2::line_mirror(Point2::ORIGIN, Direction2::Y);
1147        assert_eq!(m.kind(), TransformKind::PlaneMirror);
1148        assert!(!m.preserves_handedness());
1149        assert!(
1150            m.apply(Point2::new(3.0, 2.0))
1151                .is_equal(Point2::new(3.0, -2.0), T)
1152        );
1153
1154        assert!(Transform2::scaling(Point2::ORIGIN, 0.0, T).is_err());
1155    }
1156}
ogeom_math/transform.rs

ogeom_math/
transform.rs
