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);
78#[cfg(feature = "f64")]
79impl_symmetric_eigen2!(DSymmetricEigen2, glam::DMat2, glam::DVec2, f64);
80
81// f32 <-> f64 conversions
82#[cfg(feature = "f64")]
83impl From<SymmetricEigen2> for DSymmetricEigen2 {
84    #[inline]
85    fn from(e: SymmetricEigen2) -> Self {
86        Self {
87            eigenvalues: e.eigenvalues.as_dvec2(),
88            eigenvectors: e.eigenvectors.as_dmat2(),
89        }
90    }
91}
92
93#[cfg(feature = "f64")]
94impl From<DSymmetricEigen2> for SymmetricEigen2 {
95    #[inline]
96    fn from(e: DSymmetricEigen2) -> Self {
97        Self {
98            eigenvalues: e.eigenvalues.as_vec2(),
99            eigenvectors: e.eigenvectors.as_mat2(),
100        }
101    }
102}
103
104#[cfg(test)]
105mod tests {
106    use super::*;
107    use approx::assert_relative_eq;
108
109    #[test]
110    fn eigen_2x2() {
111        let mat = glam::Mat2::from_cols_array_2d(&[[6.0, 3.0], [3.0, 4.0]]);
112        let eigen = SymmetricEigen2::new(mat);
113
114        assert_relative_eq!(
115            eigen.eigenvalues,
116            glam::Vec2::new(8.16228, 1.83772),
117            epsilon = 0.001
118        );
119        assert_relative_eq!(
120            glam::Mat2::from_cols(eigen.eigenvectors.x_axis, eigen.eigenvectors.y_axis,),
121            glam::Mat2::from_cols(
122                glam::Vec2::new(0.811242, 0.58471),
123                glam::Vec2::new(0.58471, -0.811242),
124            ),
125            epsilon = 0.001
126        );
127    }
128
129    #[cfg(feature = "f64")]
130    #[test]
131    fn eigen_2x2_f64() {
132        let mat = glam::DMat2::from_cols_array_2d(&[[6.0, 3.0], [3.0, 4.0]]);
133        let eigen = DSymmetricEigen2::new(mat);
134
135        assert_relative_eq!(
136            eigen.eigenvalues,
137            glam::DVec2::new(8.16227766016838, 1.8377223398316202),
138            epsilon = 1e-10
139        );
140    }
141}
glamx/eigen2.rs

glamx/
eigen2.rs
