glamx/
eigen2.rs

1//! Eigendecomposition for symmetric 2x2 matrices.
2//!
3//! Adapted from <https://github.com/Jondolf/barry/tree/main/src/math/eigen>
4
5use glam::Vec2Swizzles;
6
7#[cfg(not(feature = "std"))]
8use simba::scalar::ComplexField;
9
10/// Macro to generate SymmetricEigen2 for a specific scalar type.
11macro_rules! impl_symmetric_eigen2 {
12    ($SymmetricEigen2:ident, $Mat2:ty, $Vec2:ty, $Real:ty $(, #[$attr:meta])*) => {
13        #[doc = concat!("The eigen decomposition of a symmetric 2x2 matrix (", stringify!($Real), " precision).")]
14        #[derive(Clone, Copy, Debug, PartialEq)]
15        #[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
16        $(#[$attr])*
17        pub struct $SymmetricEigen2 {
18            /// The eigenvalues of the symmetric 2x2 matrix.
19            pub eigenvalues: $Vec2,
20            /// The two eigenvectors of the symmetric 2x2 matrix (as columns).
21            pub eigenvectors: $Mat2,
22        }
23
24        impl $SymmetricEigen2 {
25            /// Computes the eigen decomposition of the given symmetric 2x2 matrix.
26            ///
27            /// The eigenvalues are returned in descending order `eigen1 > eigen2`.
28            /// This can be reversed with the [`reverse`](Self::reverse) method.
29            pub fn new(mat: $Mat2) -> Self {
30                let eigenvalues = Self::eigenvalues(mat);
31                let eigenvector1 = Self::eigenvector(mat, eigenvalues.x);
32                let eigenvector2 = Self::eigenvector(mat, eigenvalues.y);
33
34                Self {
35                    eigenvalues,
36                    eigenvectors: <$Mat2>::from_cols(eigenvector1, eigenvector2),
37                }
38            }
39
40            /// Reverses the order of the eigenvalues and their corresponding eigenvectors.
41            pub fn reverse(&self) -> Self {
42                Self {
43                    eigenvalues: self.eigenvalues.yx(),
44                    eigenvectors: <$Mat2>::from_cols(self.eigenvectors.y_axis, self.eigenvectors.x_axis),
45                }
46            }
47
48            /// Computes the eigenvalues of a symmetric 2x2 matrix.
49            ///
50            /// Reference: <https://croninprojects.org/Vince/Geodesy/FindingEigenvectors.pdf>
51            pub fn eigenvalues(mat: $Mat2) -> $Vec2 {
52                let [a, b, c] = [
53                    1.0,
54                    -(mat.x_axis.x + mat.y_axis.y),
55                    mat.x_axis.x * mat.y_axis.y - mat.x_axis.y * mat.y_axis.x,
56                ];
57                // The eigenvalues are the roots of the quadratic equation:
58                // ax^2 + bx + c = 0
59                // x = (-b +/- sqrt(b^2 - 4ac)) / 2a
60                let sqrt_part = (b.powi(2) - 4.0 * a * c).sqrt();
61                let eigen1 = (-b + sqrt_part) / (2.0 * a);
62                let eigen2 = (-b - sqrt_part) / (2.0 * a);
63                <$Vec2>::new(eigen1, eigen2)
64            }
65
66            /// Computes the unit-length eigenvector corresponding to the given `eigenvalue`
67            /// of the symmetric 2x2 `mat`.
68            ///
69            /// Reference: <https://croninprojects.org/Vince/Geodesy/FindingEigenvectors.pdf>
70            pub fn eigenvector(mat: $Mat2, eigenvalue: $Real) -> $Vec2 {
71                <$Vec2>::new(1.0, (eigenvalue - mat.x_axis.x) / mat.y_axis.x).normalize()
72            }
73        }
74    };
75}
76
77impl_symmetric_eigen2!(SymmetricEigen2, glam::Mat2, glam::Vec2, f32);
78impl_symmetric_eigen2!(DSymmetricEigen2, glam::DMat2, glam::DVec2, f64);
79
80// f32 <-> f64 conversions
81impl From<SymmetricEigen2> for DSymmetricEigen2 {
82    #[inline]
83    fn from(e: SymmetricEigen2) -> Self {
84        Self {
85            eigenvalues: e.eigenvalues.as_dvec2(),
86            eigenvectors: e.eigenvectors.as_dmat2(),
87        }
88    }
89}
90
91impl From<DSymmetricEigen2> for SymmetricEigen2 {
92    #[inline]
93    fn from(e: DSymmetricEigen2) -> Self {
94        Self {
95            eigenvalues: e.eigenvalues.as_vec2(),
96            eigenvectors: e.eigenvectors.as_mat2(),
97        }
98    }
99}
100
101#[cfg(test)]
102mod tests {
103    use super::*;
104    use approx::assert_relative_eq;
105
106    #[test]
107    fn eigen_2x2() {
108        let mat = glam::Mat2::from_cols_array_2d(&[[6.0, 3.0], [3.0, 4.0]]);
109        let eigen = SymmetricEigen2::new(mat);
110
111        assert_relative_eq!(
112            eigen.eigenvalues,
113            glam::Vec2::new(8.16228, 1.83772),
114            epsilon = 0.001
115        );
116        assert_relative_eq!(
117            glam::Mat2::from_cols(eigen.eigenvectors.x_axis, eigen.eigenvectors.y_axis,),
118            glam::Mat2::from_cols(
119                glam::Vec2::new(0.811242, 0.58471),
120                glam::Vec2::new(0.58471, -0.811242),
121            ),
122            epsilon = 0.001
123        );
124    }
125
126    #[test]
127    fn eigen_2x2_f64() {
128        let mat = glam::DMat2::from_cols_array_2d(&[[6.0, 3.0], [3.0, 4.0]]);
129        let eigen = DSymmetricEigen2::new(mat);
130
131        assert_relative_eq!(
132            eigen.eigenvalues,
133            glam::DVec2::new(8.16227766016838, 1.8377223398316202),
134            epsilon = 1e-10
135        );
136    }
137}
glamx/eigen2.rs

glamx/
eigen2.rs
