glam_matrix_extras/eigen/
symmetric_eigen2.rs

1use crate::{
2    SymmetricMat2,
3    ops::{self, FloatPow},
4};
5use glam::{Mat2, Vec2, Vec2Swizzles};
6
7/// The [eigen decomposition] of a [`SymmetricMat2`].
8///
9/// [eigen decomposition]: https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix
10#[derive(Clone, Copy, Debug, PartialEq)]
11#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
12pub struct SymmetricEigen2 {
13    /// The eigenvalues of the [`SymmetricMat2`].
14    ///
15    /// These should be in ascending order `eigen1 <= eigen2`.
16    pub eigenvalues: Vec2,
17    /// The eigenvectors of the [`SymmetricMat2`].
18    /// They should be unit length and orthogonal to each other.
19    ///
20    /// The eigenvectors are ordered to correspond to the eigenvalues. For example,
21    /// `eigenvectors.x_axis` corresponds to `eigenvalues.x`.
22    pub eigenvectors: Mat2,
23}
24
25impl SymmetricEigen2 {
26    /// Computes the eigen decomposition of the given [`SymmetricMat2`].
27    ///
28    /// The eigenvalues are returned in ascending order `eigen1 <= eigen2`.
29    /// This can be reversed with the [`reverse`](Self::reverse) method.
30    // TODO: Verify that the eigenvalues really are always returned in ascending order.
31    pub fn new(mat: SymmetricMat2) -> Self {
32        let eigenvalues = Self::eigenvalues(mat);
33        let eigenvector1 = Self::eigenvector(mat, eigenvalues.x);
34        let eigenvector2 = Self::eigenvector(mat, eigenvalues.y);
35
36        Self {
37            eigenvalues,
38            eigenvectors: Mat2::from_cols(eigenvector1, eigenvector2),
39        }
40    }
41
42    /// Reverses the order of the eigenvalues and their corresponding eigenvectors.
43    pub fn reverse(&self) -> Self {
44        Self {
45            eigenvalues: self.eigenvalues.yx(),
46            eigenvectors: Mat2::from_cols(self.eigenvectors.y_axis, self.eigenvectors.x_axis),
47        }
48    }
49
50    /// Computes the eigenvalues of a [`SymmetricMat2`].
51    ///
52    /// Reference: <https://croninprojects.org/Vince/Geodesy/FindingEigenvectors.pdf>
53    pub fn eigenvalues(mat: SymmetricMat2) -> Vec2 {
54        let [a, b, c] = [
55            1.0,
56            -(mat.m00 + mat.m11),
57            mat.m00 * mat.m11 - mat.m01 * mat.m01,
58        ];
59        // The eigenvalues are the roots of the quadratic equation:
60        // ax^2 + bx + c = 0
61        // x = (-b ± sqrt(b^2 - 4ac)) / 2a
62        let sqrt_part = ops::sqrt(b.squared() - 4.0 * a * c);
63        let eigen1 = (-b + sqrt_part) / (2.0 * a);
64        let eigen2 = (-b - sqrt_part) / (2.0 * a);
65        Vec2::new(eigen1, eigen2)
66    }
67
68    /// Computes the unit-length eigenvector corresponding to the given `eigenvalue`
69    /// of the symmetric 2x2 `mat`.
70    ///
71    /// Reference: <https://croninprojects.org/Vince/Geodesy/FindingEigenvectors.pdf>
72    pub fn eigenvector(mat: SymmetricMat2, eigenvalue: f32) -> Vec2 {
73        Vec2::new(1.0, (eigenvalue - mat.m00) / mat.m01).normalize()
74    }
75}
76
77#[cfg(test)]
78mod test {
79    use approx::assert_relative_eq;
80    use glam::{Mat2, Vec2};
81
82    use crate::SymmetricMat2;
83
84    use super::SymmetricEigen2;
85
86    #[test]
87    fn eigen_2x2() {
88        let mat = SymmetricMat2::new(6.0, 3.0, 4.0);
89        let eigen = SymmetricEigen2::new(mat);
90
91        assert_relative_eq!(
92            eigen.eigenvalues,
93            Vec2::new(8.16228, 1.83772),
94            epsilon = 0.001
95        );
96        assert_relative_eq!(
97            Mat2::from_cols(eigen.eigenvectors.x_axis, eigen.eigenvectors.y_axis,),
98            Mat2::from_cols(Vec2::new(0.811242, 0.58471), Vec2::new(0.58471, -0.811242),),
99            epsilon = 0.001
100        );
101    }
102}
glam_matrix_extras/eigen/symmetric_eigen2.rs

glam_matrix_extras/eigen/
symmetric_eigen2.rs
