dimitri
ago%!(EXTRA string=8 years)
commit
5d9bdbaded
60 changed files with 1074 additions and 288 deletions
@ -0,0 +1,65 @@ |
|||||||
|
/// @ref gtx_matrix_factorisation
|
||||||
|
/// @file glm/gtx/matrix_factorisation.hpp
|
||||||
|
///
|
||||||
|
/// @see core (dependence)
|
||||||
|
///
|
||||||
|
/// @defgroup gtx_matrix_factorisation GLM_GTX_matrix_factorisation
|
||||||
|
/// @ingroup gtx
|
||||||
|
///
|
||||||
|
/// @brief Functions to factor matrices in various forms
|
||||||
|
///
|
||||||
|
/// <glm/gtx/matrix_factorisation.hpp> need to be included to use these functionalities.
|
||||||
|
|
||||||
|
#pragma once |
||||||
|
|
||||||
|
// Dependency:
|
||||||
|
#include "../glm.hpp" |
||||||
|
|
||||||
|
#ifndef GLM_ENABLE_EXPERIMENTAL |
||||||
|
# error "GLM: GLM_GTX_matrix_factorisation is an experimental extension and may change in the future. Use #define GLM_ENABLE_EXPERIMENTAL before including it, if you really want to use it." |
||||||
|
#endif |
||||||
|
|
||||||
|
#if GLM_MESSAGES == GLM_MESSAGES_ENABLED && !defined(GLM_EXT_INCLUDED) |
||||||
|
# pragma message("GLM: GLM_GTX_matrix_factorisation extension included") |
||||||
|
#endif |
||||||
|
|
||||||
|
/*
|
||||||
|
Suggestions: |
||||||
|
- Move helper functions flipud and fliplr to another file: They may be helpful in more general circumstances. |
||||||
|
- Implement other types of matrix factorisation, such as: QL and LQ, L(D)U, eigendecompositions, etc... |
||||||
|
*/ |
||||||
|
|
||||||
|
namespace glm |
||||||
|
{ |
||||||
|
/// @addtogroup gtx_matrix_factorisation
|
||||||
|
/// @{
|
||||||
|
|
||||||
|
/// Flips the matrix rows up and down.
|
||||||
|
/// From GLM_GTX_matrix_factorisation extension.
|
||||||
|
template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType> |
||||||
|
GLM_FUNC_DECL matType<C, R, T, P> flipud(matType<C, R, T, P> const& in); |
||||||
|
|
||||||
|
/// Flips the matrix columns right and left.
|
||||||
|
/// From GLM_GTX_matrix_factorisation extension.
|
||||||
|
template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType> |
||||||
|
GLM_FUNC_DECL matType<C, R, T, P> fliplr(matType<C, R, T, P> const& in); |
||||||
|
|
||||||
|
/// Performs QR factorisation of a matrix.
|
||||||
|
/// Returns 2 matrices, q and r, such that the columns of q are orthonormal and span the same subspace than those of the input matrix, r is an upper triangular matrix, and q*r=in.
|
||||||
|
/// Given an n-by-m input matrix, q has dimensions min(n,m)-by-m, and r has dimensions n-by-min(n,m).
|
||||||
|
/// From GLM_GTX_matrix_factorisation extension.
|
||||||
|
template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType> |
||||||
|
GLM_FUNC_DECL void qr_decompose(matType<C, R, T, P> const& in, matType<(C < R ? C : R), R, T, P>& q, matType<C, (C < R ? C : R), T, P>& r); |
||||||
|
|
||||||
|
/// Performs RQ factorisation of a matrix.
|
||||||
|
/// Returns 2 matrices, r and q, such that r is an upper triangular matrix, the rows of q are orthonormal and span the same subspace than those of the input matrix, and r*q=in.
|
||||||
|
/// Note that in the context of RQ factorisation, the diagonal is seen as starting in the lower-right corner of the matrix, instead of the usual upper-left.
|
||||||
|
/// Given an n-by-m input matrix, r has dimensions min(n,m)-by-m, and q has dimensions n-by-min(n,m).
|
||||||
|
/// From GLM_GTX_matrix_factorisation extension.
|
||||||
|
template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType> |
||||||
|
GLM_FUNC_DECL void rq_decompose(matType<C, R, T, P> const& in, matType<(C < R ? C : R), R, T, P>& r, matType<C, (C < R ? C : R), T, P>& q); |
||||||
|
|
||||||
|
/// @}
|
||||||
|
} |
||||||
|
|
||||||
|
#include "matrix_factorisation.inl" |
||||||
@ -0,0 +1,85 @@ |
|||||||
|
/// @ref gtx_matrix_factorisation |
||||||
|
/// @file glm/gtx/matrix_factorisation.inl |
||||||
|
|
||||||
|
namespace glm |
||||||
|
{ |
||||||
|
template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType> |
||||||
|
GLM_FUNC_QUALIFIER matType<C, R, T, P> flipud(matType<C, R, T, P> const& in) |
||||||
|
{ |
||||||
|
matType<R, C, T, P> tin = transpose(in); |
||||||
|
tin = fliplr(tin); |
||||||
|
matType<C, R, T, P> out = transpose(tin); |
||||||
|
|
||||||
|
return out; |
||||||
|
} |
||||||
|
|
||||||
|
template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType> |
||||||
|
GLM_FUNC_QUALIFIER matType<C, R, T, P> fliplr(matType<C, R, T, P> const& in) |
||||||
|
{ |
||||||
|
matType<C, R, T, P> out(uninitialize); |
||||||
|
for (length_t i = 0; i < C; i++) |
||||||
|
{ |
||||||
|
out[i] = in[(C - i) - 1]; |
||||||
|
} |
||||||
|
|
||||||
|
return out; |
||||||
|
} |
||||||
|
|
||||||
|
template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType> |
||||||
|
GLM_FUNC_QUALIFIER void qr_decompose(matType<C, R, T, P> const& in, matType<(C < R ? C : R), R, T, P>& q, matType<C, (C < R ? C : R), T, P>& r) |
||||||
|
{ |
||||||
|
// Uses modified Gram-Schmidt method |
||||||
|
// Source: https://en.wikipedia.org/wiki/Gram–Schmidt_process |
||||||
|
// And https://en.wikipedia.org/wiki/QR_decomposition |
||||||
|
|
||||||
|
//For all the linearly independs columns of the input... |
||||||
|
// (there can be no more linearly independents columns than there are rows.) |
||||||
|
for (length_t i = 0; i < (C < R ? C : R); i++) |
||||||
|
{ |
||||||
|
//Copy in Q the input's i-th column. |
||||||
|
q[i] = in[i]; |
||||||
|
|
||||||
|
//j = [0,i[ |
||||||
|
// Make that column orthogonal to all the previous ones by substracting to it the non-orthogonal projection of all the previous columns. |
||||||
|
// Also: Fill the zero elements of R |
||||||
|
for (length_t j = 0; j < i; j++) |
||||||
|
{ |
||||||
|
q[i] -= dot(q[i], q[j])*q[j]; |
||||||
|
r[j][i] = 0; |
||||||
|
} |
||||||
|
|
||||||
|
//Now, Q i-th column is orthogonal to all the previous columns. Normalize it. |
||||||
|
q[i] = normalize(q[i]); |
||||||
|
|
||||||
|
//j = [i,C[ |
||||||
|
//Finally, compute the corresponding coefficients of R by computing the projection of the resulting column on the other columns of the input. |
||||||
|
for (length_t j = i; j < C; j++) |
||||||
|
{ |
||||||
|
r[j][i] = dot(in[j], q[i]); |
||||||
|
} |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
template <length_t C, length_t R, typename T, precision P, template<length_t, length_t, typename, precision> class matType> |
||||||
|
GLM_FUNC_QUALIFIER void rq_decompose(matType<C, R, T, P> const& in, matType<(C < R ? C : R), R, T, P>& r, matType<C, (C < R ? C : R), T, P>& q) |
||||||
|
{ |
||||||
|
// From https://en.wikipedia.org/wiki/QR_decomposition: |
||||||
|
// The RQ decomposition transforms a matrix A into the product of an upper triangular matrix R (also known as right-triangular) and an orthogonal matrix Q. The only difference from QR decomposition is the order of these matrices. |
||||||
|
// QR decomposition is Gram–Schmidt orthogonalization of columns of A, started from the first column. |
||||||
|
// RQ decomposition is Gram–Schmidt orthogonalization of rows of A, started from the last row. |
||||||
|
|
||||||
|
matType<R, C, T, P> tin = transpose(in); |
||||||
|
tin = fliplr(tin); |
||||||
|
|
||||||
|
matType<R, (C < R ? C : R), T, P> tr; |
||||||
|
matType<(C < R ? C : R), C, T, P> tq; |
||||||
|
qr_decompose(tin, tq, tr); |
||||||
|
|
||||||
|
tr = fliplr(tr); |
||||||
|
r = transpose(tr); |
||||||
|
r = fliplr(r); |
||||||
|
|
||||||
|
tq = fliplr(tq); |
||||||
|
q = transpose(tq); |
||||||
|
} |
||||||
|
} //namespace glm |
||||||
@ -0,0 +1,101 @@ |
|||||||
|
#define GLM_ENABLE_EXPERIMENTAL |
||||||
|
#include <glm/gtx/matrix_factorisation.hpp> |
||||||
|
|
||||||
|
float const epsilon = 1e-10f; |
||||||
|
|
||||||
|
template <glm::length_t C, glm::length_t R, typename T, glm::precision P, template<glm::length_t, glm::length_t, typename, glm::precision> class matType> |
||||||
|
int test_qr(matType<C, R, T, P> m) |
||||||
|
{ |
||||||
|
int Error = 0; |
||||||
|
|
||||||
|
matType<(C < R ? C : R), R, T, P> q(-999); |
||||||
|
matType<C, (C < R ? C : R), T, P> r(-999); |
||||||
|
|
||||||
|
glm::qr_decompose(m, q, r); |
||||||
|
|
||||||
|
//Test if q*r really equals the input matrix
|
||||||
|
matType<C, R, T, P> tm = q*r; |
||||||
|
matType<C, R, T, P> err = tm - m; |
||||||
|
|
||||||
|
for (glm::length_t i = 0; i < C; i++) |
||||||
|
for (glm::length_t j = 0; j < R; j++) |
||||||
|
Error += std::abs(err[i][j]) > epsilon ? 1 : 0; |
||||||
|
|
||||||
|
//Test if the columns of q are orthonormal
|
||||||
|
for (glm::length_t i = 0; i < (C < R ? C : R); i++) |
||||||
|
{ |
||||||
|
Error += (length(q[i]) - 1) > epsilon ? 1 : 0; |
||||||
|
|
||||||
|
for (glm::length_t j = 0; j<i; j++) |
||||||
|
Error += std::abs(dot(q[i], q[j])) > epsilon ? 1 : 0; |
||||||
|
} |
||||||
|
|
||||||
|
//Test if the matrix r is upper triangular
|
||||||
|
for (glm::length_t i = 0; i < C; i++) |
||||||
|
for (glm::length_t j = i + 1; j < (C < R ? C : R); j++) |
||||||
|
Error += r[i][j] != 0 ? 1 : 0; |
||||||
|
|
||||||
|
return Error; |
||||||
|
} |
||||||
|
|
||||||
|
template <glm::length_t C, glm::length_t R, typename T, glm::precision P, template<glm::length_t, glm::length_t, typename, glm::precision> class matType> |
||||||
|
int test_rq(matType<C, R, T, P> m) |
||||||
|
{ |
||||||
|
int Error = 0; |
||||||
|
|
||||||
|
matType<C, (C < R ? C : R), T, P> q(-999); |
||||||
|
matType<(C < R ? C : R), R, T, P> r(-999); |
||||||
|
|
||||||
|
glm::rq_decompose(m, r, q); |
||||||
|
|
||||||
|
//Test if q*r really equals the input matrix
|
||||||
|
matType<C, R, T, P> tm = r*q; |
||||||
|
matType<C, R, T, P> err = tm - m; |
||||||
|
|
||||||
|
for (glm::length_t i = 0; i < C; i++) |
||||||
|
for (glm::length_t j = 0; j < R; j++) |
||||||
|
Error += std::abs(err[i][j]) > epsilon ? 1 : 0; |
||||||
|
|
||||||
|
//Test if the rows of q are orthonormal
|
||||||
|
matType<(C < R ? C : R), C, T, P> tq = transpose(q); |
||||||
|
|
||||||
|
for (glm::length_t i = 0; i < (C < R ? C : R); i++) |
||||||
|
{ |
||||||
|
Error += (length(tq[i]) - 1) > epsilon ? 1 : 0; |
||||||
|
|
||||||
|
for (glm::length_t j = 0; j<i; j++) |
||||||
|
Error += std::abs(dot(tq[i], tq[j])) > epsilon ? 1 : 0; |
||||||
|
} |
||||||
|
|
||||||
|
//Test if the matrix r is upper triangular
|
||||||
|
for (glm::length_t i = 0; i < (C < R ? C : R); i++) |
||||||
|
for (glm::length_t j = R - (C < R ? C : R) + i + 1; j < R; j++) |
||||||
|
Error += r[i][j] != 0 ? 1 : 0; |
||||||
|
|
||||||
|
return Error; |
||||||
|
} |
||||||
|
|
||||||
|
int main() |
||||||
|
{ |
||||||
|
int Error = 0; |
||||||
|
|
||||||
|
//Test QR square
|
||||||
|
Error += test_qr(glm::dmat3(12, 6, -4, -51, 167, 24, 4, -68, -41)) ? 1 : 0; |
||||||
|
|
||||||
|
//Test RQ square
|
||||||
|
Error += test_rq(glm::dmat3(12, 6, -4, -51, 167, 24, 4, -68, -41)) ? 1 : 0; |
||||||
|
|
||||||
|
//Test QR triangular 1
|
||||||
|
Error += test_qr(glm::dmat3x4(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15)) ? 1 : 0; |
||||||
|
|
||||||
|
//Test QR triangular 2
|
||||||
|
Error += test_qr(glm::dmat4x3(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15)) ? 1 : 0; |
||||||
|
|
||||||
|
//Test RQ triangular 1 : Fails at the triangular test
|
||||||
|
Error += test_rq(glm::dmat3x4(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15)) ? 1 : 0; |
||||||
|
|
||||||
|
//Test QR triangular 2
|
||||||
|
Error += test_rq(glm::dmat4x3(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15)) ? 1 : 0; |
||||||
|
|
||||||
|
return Error; |
||||||
|
} |
||||||
Loading…
Reference in New Issue