float2x2: eigen decomposition

1 files changed, 59 insertions, 0 deletions

diff --git a/source/blender/blenlib/BLI_float2x2.hh b/source/blender/blenlib/BLI_float2x2.hh
index b64dc86a1ee..4d4f9ae522a 100644
--- a/source/blender/blenlib/BLI_float2x2.hh
+++ b/source/blender/blenlib/BLI_float2x2.hh
@@ -20,6 +20,9 @@
 #include "BLI_float2.hh"
 #include "BLI_math_matrix.h"
 
+#include <cmath>
+#include <tuple>
+
 namespace blender {
 
 /**
@@ -133,6 +136,62 @@ struct float2x2 {
     return float2x2(res);
   }
 
+  /**
+   * Computes the eigen decomposition of the 2x2 matrix and returns Q
+   * and Lambda
+   *
+   * A = Q * Lambda * Q.inverse()
+   *
+   * A is this matrix.
+   *
+   * Q is a 2x2 matrix whose ith column is the eigen vector q_i of the
+   * matrix A.
+   *
+   * Lambda is the diagonal matrix whose diagonal elements are the
+   * corresponding eigenvalues.
+   *
+   * Reference https://en.wikipedia.org/wiki/Eigenvalue_algorithm
+   * http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/index.html
+   */
+  std::tuple<float2x2, float2> eigen_decomposition() const
+  {
+    const auto a = this->ptr()[0][0];
+    const auto b = this->ptr()[1][0];
+    const auto c = this->ptr()[0][1];
+    const auto d = this->ptr()[1][1];
+
+    const auto trace = a + d;
+    const auto det = a * d - b * c;
+
+    const auto l1 = (trace + std::sqrt(trace * trace - 4.0 * det)) / 2.0;
+    const auto l2 = (trace - std::sqrt(trace * trace - 4.0 * det)) / 2.0;
+    const auto lambda = float2(l1, l2);
+
+    if (c != 0.0) {
+      float2x2 q;
+      q.ptr()[0][0] = l1 - d;
+      q.ptr()[1][0] = l2 - d;
+      q.ptr()[0][1] = c;
+      q.ptr()[1][1] = c;
+      return {q, lambda};
+    }
+    if (b != 0.0) {
+      float2x2 q;
+      q.ptr()[0][0] = b;
+      q.ptr()[1][0] = b;
+      q.ptr()[0][1] = l1 - a;
+      q.ptr()[1][1] = l2 - a;
+      return {q, lambda};
+    }
+
+    float2x2 q;
+    q.ptr()[0][0] = 1;
+    q.ptr()[1][0] = 0;
+    q.ptr()[0][1] = 0;
+    q.ptr()[1][1] = 1;
+    return {q, lambda};
+  }
+
   uint64_t hash() const
   {
     uint64_t h = 435109;