Math: added Gram-Schmidt orthogonalization algorithm.

src/Math/Algorithms/GramSchmidt.h

@@ -16,13 +16,36 @@
 */
 
 /** @file
- * @brief Function Magnum::Math::Algorithms::gramSchmidtOrthonormalizeInPlace(), Magnum::Math::Algorithms::gramSchmidtOrthonormalize()
+ * @brief Function Magnum::Math::Algorithms::gramSchmidtOrthogonalizeInPlace(), Magnum::Math::Algorithms::gramSchmidtOrthogonalize(), Magnum::Math::Algorithms::gramSchmidtOrthonormalizeInPlace(), Magnum::Math::Algorithms::gramSchmidtOrthonormalize()
  */
 
 #include "Math/RectangularMatrix.h"
 
 namespace Magnum { namespace Math { namespace Algorithms {
 
+/**
+@brief In-place Gram-Schmidt matrix orthogonalization
+@param[in,out] matrix   Matrix to perform orthogonalization on
+*/
+template<std::size_t cols, std::size_t rows, class T> void gramSchmidtOrthogonalizeInPlace(RectangularMatrix<cols, rows, T>& matrix) {
+    static_assert(cols <= rows, "Unsupported matrix aspect ratio");
+    for(std::size_t i = 0; i != cols; ++i) {
+        for(std::size_t j = i+1; j != cols; ++j)
+            matrix[j] -= matrix[j].projected(matrix[i]);
+    }
+}
+
+/**
+@brief Gram-Schmidt matrix orthogonalization
+
+Unlike gramSchmidtOrthogonalizeInPlace() returns the modified matrix instead
+of performing the orthogonalization in-place.
+*/
+template<std::size_t cols, std::size_t rows, class T> RectangularMatrix<cols, rows, T> gramSchmidtOrthogonalize(RectangularMatrix<cols, rows, T> matrix) {
+    gramSchmidtOrthogonalizeInPlace(matrix);
+    return matrix;
+}
+
 /**
 @brief In-place Gram-Schmidt matrix orthonormalization
 @param[in,out] matrix   Matrix to perform orthonormalization on

src/Math/Algorithms/Test/GramSchmidtTest.cpp

@@ -23,6 +23,7 @@ class GramSchmidtTest: public Corrade::TestSuite::Tester {
     public:
         GramSchmidtTest();
 
+        void orthogonalize();
         void orthonormalize();
 };
 
@@ -30,7 +31,32 @@ typedef RectangularMatrix<3, 3, float> Matrix3;
 typedef Vector<3, float> Vector3;
 
 GramSchmidtTest::GramSchmidtTest() {
-    addTests(&GramSchmidtTest::orthonormalize);
+    addTests(&GramSchmidtTest::orthogonalize,
+             &GramSchmidtTest::orthonormalize);
+}
+
+void GramSchmidtTest::orthogonalize() {
+    Matrix3 m(Vector3(3.0f,  5.0f, 1.0f),
+              Vector3(4.0f,  4.0f, 7.0f),
+              Vector3(7.0f, -1.0f, -3.0f));
+
+    Matrix3 orthogonalized = Algorithms::gramSchmidtOrthogonalize(m);
+
+    /* Verify the first vector is in direction of first original */
+    CORRADE_COMPARE(orthogonalized[0], m[0]);
+
+    /* (The vectors don't need to unit length) */
+
+    /* Verify the vectors are orthogonal */
+    CORRADE_COMPARE(Vector3::dot(orthogonalized[0], orthogonalized[1]), 0.0f);
+    CORRADE_COMPARE(Vector3::dot(orthogonalized[0], orthogonalized[2]), 0.0f);
+    CORRADE_COMPARE(Vector3::dot(orthogonalized[1], orthogonalized[2]), 0.0f);
+
+    /* Just to be sure */
+    Matrix3 expected(Vector3(     3.0f,       5.0f,       1.0f),
+                     Vector3(0.657143f, -1.571429f,  5.885714f),
+                     Vector3(6.086759f,   -3.3379f, -1.570777f));
+    CORRADE_COMPARE(orthogonalized, expected);
 }
 
 void GramSchmidtTest::orthonormalize() {
@@ -38,26 +64,26 @@ void GramSchmidtTest::orthonormalize() {
               Vector3(4.0f,  4.0f, 7.0f),
               Vector3(7.0f, -1.0f, 8.0f));
 
-    Matrix3 normalized = Algorithms::gramSchmidtOrthonormalize(m);
+    Matrix3 orthonormalized = Algorithms::gramSchmidtOrthonormalize(m);
 
     /* Verify the first vector is in direction of first original */
-    CORRADE_COMPARE(normalized[0], m[0].normalized());
+    CORRADE_COMPARE(orthonormalized[0], m[0].normalized());
 
     /* Verify the vectors have unit length */
-    CORRADE_COMPARE(normalized[0].length(), 1.0f);
+    CORRADE_COMPARE(orthonormalized[0].length(), 1.0f);
-    CORRADE_COMPARE(normalized[1].length(), 1.0f);
+    CORRADE_COMPARE(orthonormalized[1].length(), 1.0f);
-    CORRADE_COMPARE(normalized[2].length(), 1.0f);
+    CORRADE_COMPARE(orthonormalized[2].length(), 1.0f);
 
     /* Verify the vectors are orthogonal */
-    CORRADE_COMPARE(Vector3::dot(normalized[0], normalized[1]), 0.0f);
+    CORRADE_COMPARE(Vector3::dot(orthonormalized[0], orthonormalized[1]), 0.0f);
-    CORRADE_COMPARE(Vector3::dot(normalized[0], normalized[2]), 0.0f);
+    CORRADE_COMPARE(Vector3::dot(orthonormalized[0], orthonormalized[2]), 0.0f);
-    CORRADE_COMPARE(Vector3::dot(normalized[1], normalized[2]), 0.0f);
+    CORRADE_COMPARE(Vector3::dot(orthonormalized[1], orthonormalized[2]), 0.0f);
 
     /* Just to be sure */
     Matrix3 expected(Vector3( 0.303046f,  0.505076f,  0.808122f),
                      Vector3( 0.928316f, -0.348119f, -0.130544f),
                      Vector3(-0.215388f, -0.789754f,  0.574367f));
-    CORRADE_COMPARE(normalized, expected);
+    CORRADE_COMPARE(orthonormalized, expected);
 }
 
 }}}}