Implement Matrix3x3::InverseTranpose() InverseTransposed() Inverse() InverseFast()

include/J3ML/LinearAlgebra/Matrix3x3.h

@@ -341,8 +341,13 @@ namespace J3ML::LinearAlgebra {
         Matrix3x3 Transposed() const;
 
         /// Returns the inverse transpose of this matrix.
+        /// Use the matrix to transform covariant vectors (normal vectors).
         Matrix3x3 InverseTransposed() const;
 
+        /// Computes the inverse transpose of this matrix.
+        /// Use the matrix to transform covariant vectors (normal vectors).
+        bool InverseTranspose();
+
         /// Inverts this matrix using numerically stable Gaussian elimination.
         /// @return Returns true on success, false otherwise;
         bool Inverse(float epsilon = 1e-6f);
@@ -382,7 +387,7 @@ namespace J3ML::LinearAlgebra {
 
         void Transpose();
 
-        bool InverseTranspose();
+
 
         void RemoveScale();
 

src/J3ML/LinearAlgebra/Matrix3x3.cpp

@@ -639,6 +639,57 @@ namespace J3ML::LinearAlgebra {
         return a * (d*f - e*e) + b * (2.f * c * e - b * f) - c*c*d;
     }
 
+    Matrix3x3 Matrix3x3::InverseTransposed() const {
+        Matrix3x3 copy = *this;
+        copy.Transpose();
+        copy.Inverse();
+        return copy;
+    }
+
+    bool Matrix3x3::InverseTranspose() {
+        bool success = Inverse();
+        Transpose();
+        return success;
+    }
+
+    bool Matrix3x3::Inverse(float epsilon) {
+        /// There exists a generic matrix inverse calculator that uses Gaussian elimination
+
+        Matrix3x3 i = *this;
+        bool success = InverseMatrix(i, epsilon);
+        if (!success)
+            return false;
+
+        *this = i;
+        return true;
+    }
+
+    bool Matrix3x3::InverseFast(float epsilon) {
+        // Compute the inverse directly using Cramer's rule.
+        // Warning: This method is numerically very unstable!
+        float d = Determinant();
+        if (Math::EqualAbs(d, 0.f, epsilon))
+            return true;
+
+        d = 1.f / d;
+        Matrix3x3 i;
+        i[0][0] = d * (At(1, 1) * At(2, 2) - At(1, 2) * At(2, 1));
+        i[0][1] = d * (At(0, 2) * At(2, 1) - At(0, 1) * At(2, 2));
+        i[0][2] = d * (At(0, 1) * At(1, 2) - At(0, 2) * At(1, 1));
+
+        i[1][0] = d * (At(1, 2) * At(2, 0) - At(1, 0) * At(2, 2));
+        i[1][1] = d * (At(0, 0) * At(2, 2) - At(0, 2) * At(2, 0));
+        i[1][2] = d * (At(0, 2) * At(1, 0) - At(0, 0) * At(1, 2));
+
+        i[2][0] = d * (At(1, 0) * At(2, 1) - At(1, 1) * At(2, 0));
+        i[2][1] = d * (At(2, 0) * At(0, 1) - At(0, 0) * At(2, 1));
+        i[2][2] = d * (At(0, 0) * At(1, 1) - At(0, 1) * At(1, 0));
+
+
+        *this = i;
+        return true;
+    }
+
 
 }