From a46afd6ed54ce917f60a52fdc1fd35128cf1d693 Mon Sep 17 00:00:00 2001 From: Sam Hocevar Date: Wed, 30 Jul 2014 18:15:08 +0000 Subject: [PATCH] math: make cofactor computation and matrix inversion simpler and more generic. --- src/lol/math/matrix.h | 84 ++++++++++++++++++++++++----- src/math/matrix.cpp | 120 ------------------------------------------ 2 files changed, 71 insertions(+), 133 deletions(-) diff --git a/src/lol/math/matrix.h b/src/lol/math/matrix.h index 13cf9c97..e58fe235 100644 --- a/src/lol/math/matrix.h +++ b/src/lol/math/matrix.h @@ -126,8 +126,6 @@ private: }; static_assert(sizeof(imat2) == 16, "sizeof(imat2) == 16"); -static_assert(sizeof(umat2) == 16, "sizeof(umat2) == 16"); - #if LOL_FEATURE_CXX11_UNRESTRICTED_UNIONS static_assert(sizeof(f16mat2) == 8, "sizeof(f16mat2) == 8"); #endif @@ -249,8 +247,6 @@ private: }; static_assert(sizeof(imat3) == 36, "sizeof(imat3) == 36"); -static_assert(sizeof(umat3) == 36, "sizeof(umat3) == 36"); - #if LOL_FEATURE_CXX11_UNRESTRICTED_UNIONS static_assert(sizeof(f16mat3) == 18, "sizeof(f16mat3) == 18"); #endif @@ -419,21 +415,15 @@ private: }; static_assert(sizeof(imat4) == 64, "sizeof(imat4) == 64"); -static_assert(sizeof(umat4) == 64, "sizeof(umat4) == 64"); - #if LOL_FEATURE_CXX11_UNRESTRICTED_UNIONS static_assert(sizeof(f16mat4) == 32, "sizeof(f16mat4) == 32"); #endif static_assert(sizeof(mat4) == 64, "sizeof(mat4) == 64"); static_assert(sizeof(dmat4) == 128, "sizeof(dmat4) == 128"); -template T determinant(mat_t const &); -template T determinant(mat_t const &); -template T determinant(mat_t const &); - -template mat_t inverse(mat_t const &); -template mat_t inverse(mat_t const &); -template mat_t inverse(mat_t const &); +/* + * Transpose any matrix + */ template static inline mat_t transpose(mat_t const &m) @@ -445,6 +435,74 @@ static inline mat_t transpose(mat_t const &m) return ret; } +/* + * Compute a square submatrix, useful for minors and cofactor matrices + */ + +template +mat_t submatrix(mat_t const &m, int i, int j) +{ + ASSERT(i >= 0); ASSERT(j >= 0); ASSERT(i < N); ASSERT(j < N); + + mat_t ret; + for (int i2 = 0; i2 < N - 1; ++i2) + for (int j2 = 0; j2 < N - 1; ++j2) + ret[i2][j2] = m[i2 + (i2 >= i)][j2 + (j2 >= j)]; + + return ret; +} + +/* + * Compute square matrix cofactor + */ + +template +T cofactor(mat_t const &m, int i, int j) +{ + ASSERT(i >= 0); ASSERT(j >= 0); ASSERT(i < N); ASSERT(j < N); + T tmp = determinant(submatrix(m, i, j)); + return ((i + j) & 1) ? -tmp : tmp; +} + +/* + * Compute square matrix determinant, with a specialisation for 1×1 matrices + */ + +template +T determinant(mat_t const &m) +{ + T ret = (T)0; + for (int i = 0; i < N; ++i) + ret += m[i][0] * cofactor(m, i, 0); + return ret; +} + +template +T const & determinant(mat_t const &m) +{ + return m[0][0]; +} + +/* + * Compute square matrix inverse + */ + +template +mat_t inverse(mat_t const &m) +{ + mat_t ret; + + T d = determinant(m); + if (d) + { + /* Transpose result matrix on the fly */ + for (int i = 0; i < N; ++i) + for (int j = 0; j < N; ++j) + ret[i][j] = cofactor(m, j, i) / d; + } + return ret; +} + /* * Matrix-vector and vector-matrix multiplication */ diff --git a/src/math/matrix.cpp b/src/math/matrix.cpp index ef2f1bbf..9135c136 100644 --- a/src/math/matrix.cpp +++ b/src/math/matrix.cpp @@ -17,126 +17,6 @@ namespace lol { -/* - * Return the determinant of a 2×2 matrix. - */ -static inline float det2(float a, float b, - float c, float d) -{ - return a * d - b * c; -} - -/* - * Return the determinant of a 3×3 matrix. - */ -static inline float det3(float a, float b, float c, - float d, float e, float f, - float g, float h, float i) -{ - return a * (e * i - h * f) - + b * (f * g - i * d) - + c * (d * h - g * e); -} - -/* - * Return the cofactor of the (i,j) entry in a 2×2 matrix. - */ -static inline float cofact(mat2 const &m, int i, int j) -{ - float tmp = m[(i + 1) & 1][(j + 1) & 1]; - return ((i + j) & 1) ? -tmp : tmp; -} - -/* - * Return the cofactor of the (i,j) entry in a 3×3 matrix. - */ -static inline float cofact(mat3 const &m, int i, int j) -{ - return det2(m[(i + 1) % 3][(j + 1) % 3], - m[(i + 2) % 3][(j + 1) % 3], - m[(i + 1) % 3][(j + 2) % 3], - m[(i + 2) % 3][(j + 2) % 3]); -} - -/* - * Return the cofactor of the (i,j) entry in a 4×4 matrix. - */ -static inline float cofact(mat4 const &m, int i, int j) -{ - return det3(m[(i + 1) & 3][(j + 1) & 3], - m[(i + 2) & 3][(j + 1) & 3], - m[(i + 3) & 3][(j + 1) & 3], - m[(i + 1) & 3][(j + 2) & 3], - m[(i + 2) & 3][(j + 2) & 3], - m[(i + 3) & 3][(j + 2) & 3], - m[(i + 1) & 3][(j + 3) & 3], - m[(i + 2) & 3][(j + 3) & 3], - m[(i + 3) & 3][(j + 3) & 3]) * (((i + j) & 1) ? -1.0f : 1.0f); -} - -template<> float determinant(mat2 const &m) -{ - return det2(m[0][0], m[0][1], - m[1][0], m[1][1]); -} - -template<> mat2 inverse(mat2 const &m) -{ - mat2 ret; - float d = determinant(m); - if (d) - { - d = 1.0f / d; - for (int j = 0; j < 2; j++) - for (int i = 0; i < 2; i++) - ret[j][i] = cofact(m, i, j) * d; - } - return ret; -} - -template<> float determinant(mat3 const &m) -{ - return det3(m[0][0], m[0][1], m[0][2], - m[1][0], m[1][1], m[1][2], - m[2][0], m[2][1], m[2][2]); -} - -template<> mat3 inverse(mat3 const &m) -{ - mat3 ret; - float d = determinant(m); - if (d) - { - d = 1.0f / d; - for (int j = 0; j < 3; j++) - for (int i = 0; i < 3; i++) - ret[j][i] = cofact(m, i, j) * d; - } - return ret; -} - -template<> float determinant(mat4 const &m) -{ - float ret = 0; - for (int n = 0; n < 4; n++) - ret += m[n][0] * cofact(m, n, 0); - return ret; -} - -template<> mat4 inverse(mat4 const &m) -{ - mat4 ret; - float d = determinant(m); - if (d) - { - d = 1.0f / d; - for (int j = 0; j < 4; j++) - for (int i = 0; i < 4; i++) - ret[j][i] = cofact(m, i, j) * d; - } - return ret; -} - template<> mat3 mat3::scale(float x) { mat3 ret(1.0f);