From b63adf24042c620fa1dfedd48860d3b5cce7e0b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vladim=C3=ADr=20Vondru=C5=A1?= Date: Wed, 6 Feb 2013 11:56:03 +0100 Subject: [PATCH] Math: added Gram-Schmidt orthogonalization algorithm. --- src/Math/Algorithms/GramSchmidt.h | 25 ++++++++++- src/Math/Algorithms/Test/GramSchmidtTest.cpp | 46 +++++++++++++++----- 2 files changed, 60 insertions(+), 11 deletions(-) diff --git a/src/Math/Algorithms/GramSchmidt.h b/src/Math/Algorithms/GramSchmidt.h index 905641766..b27f51b3c 100644 --- a/src/Math/Algorithms/GramSchmidt.h +++ b/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 void gramSchmidtOrthogonalizeInPlace(RectangularMatrix& 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 RectangularMatrix gramSchmidtOrthogonalize(RectangularMatrix matrix) { + gramSchmidtOrthogonalizeInPlace(matrix); + return matrix; +} + /** @brief In-place Gram-Schmidt matrix orthonormalization @param[in,out] matrix Matrix to perform orthonormalization on diff --git a/src/Math/Algorithms/Test/GramSchmidtTest.cpp b/src/Math/Algorithms/Test/GramSchmidtTest.cpp index 85700deaa..67e1e7595 100644 --- a/src/Math/Algorithms/Test/GramSchmidtTest.cpp +++ b/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(normalized[1].length(), 1.0f); - CORRADE_COMPARE(normalized[2].length(), 1.0f); + CORRADE_COMPARE(orthonormalized[0].length(), 1.0f); + CORRADE_COMPARE(orthonormalized[1].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(normalized[0], normalized[2]), 0.0f); - CORRADE_COMPARE(Vector3::dot(normalized[1], normalized[2]), 0.0f); + CORRADE_COMPARE(Vector3::dot(orthonormalized[0], orthonormalized[1]), 0.0f); + CORRADE_COMPARE(Vector3::dot(orthonormalized[0], orthonormalized[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); } }}}}