mirror of
https://github.com/g-truc/glm.git
synced 2024-11-25 18:04:34 +00:00
Implements QR and RQ matrix decomposition functions.
This commit is contained in:
Vincent Aymong
parent
2dc6196467
commit
d6abdb7935
2
.gitignore
vendored
2
.gitignore
vendored
@ -52,3 +52,5 @@ Makefile
|
||||
# local build(s)
|
||||
build*
|
||||
|
||||
/.vs
|
||||
/CMakeSettings.json
|
||||
|
||||
65
glm/gtx/matrix_factorisation.hpp
Normal file
65
glm/gtx/matrix_factorisation.hpp
Normal file
@ -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 <algorithm>
|
||||
#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 flip lr to another file: They may be helpful in more general circumstances.
|
||||
- When rq_decompose is fed a matrix that has more rows than columns, the resulting r matrix is NOT upper triangular. Is that a bug?
|
||||
- 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(const matType<C, R, T, P>& 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(const matType<C, R, T, P>& in);
|
||||
|
||||
/// Performs QR factorisation of a matrix.
|
||||
/// Returns 2 matrices, q and r, such that q columns are orthonormal, r is an upper triangular matrix, and q*r=in.
|
||||
/// r is a square matrix whose dimensions are the same than the width of the input matrix, and q has the same dimensions than the input matrix.
|
||||
/// 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<std::min(C, R), R, T, P>& q, matType<C, std::min(C, R), T, P>& r, const matType<C, R, T, P>& in);
|
||||
|
||||
/// Performs RQ factorisation of a matrix.
|
||||
/// Returns 2 matrices, r and q, such that r is an upper triangular matrix, q rows are orthonormal, and r*q=in.
|
||||
/// q has the same dimensions than the input matrix, and r is a square matrix whose dimensions are the same than the height of the input matrix.
|
||||
/// 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<std::min(C, R), R, T, P>& r, matType<C, std::min(C, R), T, P>& q, const matType<C, R, T, P>& in);
|
||||
|
||||
/// @}
|
||||
}
|
||||
|
||||
#include "matrix_factorisation.inl"
|
||||
74
glm/gtx/matrix_factorisation.inl
Normal file
74
glm/gtx/matrix_factorisation.inl
Normal file
@ -0,0 +1,74 @@
|
||||
/// @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(const matType<C, R, T, P>& 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(const matType<C, R, T, P>& in) {
|
||||
constexpr length_t num_cols = C;
|
||||
|
||||
matType<C, R, T, P> out;
|
||||
for (length_t i = 0; i < num_cols; i++) {
|
||||
out[i] = in[(num_cols - 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<std::min(C, R), R, T, P>& q, matType<C, std::min(C, R), T, P>& r, const matType<C, R, T, P>& in) {
|
||||
// Uses modified Gram-Schmidt method
|
||||
// Source: https://en.wikipedia.org/wiki/Gram–Schmidt_process
|
||||
// And https://en.wikipedia.org/wiki/QR_decomposition
|
||||
|
||||
for (length_t i = 0; i < std::min(R, C); i++) {
|
||||
q[i] = in[i];
|
||||
|
||||
for (length_t j = 0; j < i; j++) {
|
||||
q[i] -= dot(q[i], q[j])*q[j];
|
||||
}
|
||||
|
||||
q[i] = normalize(q[i]);
|
||||
}
|
||||
|
||||
for (length_t i = 0; i < std::min(R, C); i++) {
|
||||
for (length_t j = 0; j < i; j++) {
|
||||
r[j][i] = 0;
|
||||
}
|
||||
|
||||
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<std::min(C, R), R, T, P>& r, matType<C, std::min(C, R), T, P>& q, const matType<C, R, T, P>& in) {
|
||||
// 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, std::min(C, R), T, P> tr;
|
||||
matType<std::min(C, R), C, T, P> tq;
|
||||
qr_decompose(tq, tr, tin);
|
||||
|
||||
tr = fliplr(tr);
|
||||
r = transpose(tr);
|
||||
r = fliplr(r);
|
||||
|
||||
tq = fliplr(tq);
|
||||
q = transpose(tq);
|
||||
}
|
||||
} //namespace glm
|
||||
@ -21,6 +21,7 @@ glmCreateTestGTC(gtx_io)
|
||||
glmCreateTestGTC(gtx_log_base)
|
||||
glmCreateTestGTC(gtx_matrix_cross_product)
|
||||
glmCreateTestGTC(gtx_matrix_decompose)
|
||||
glmCreateTestGTC(gtx_matrix_factorisation)
|
||||
glmCreateTestGTC(gtx_matrix_interpolation)
|
||||
glmCreateTestGTC(gtx_matrix_major_storage)
|
||||
glmCreateTestGTC(gtx_matrix_operation)
|
||||
|
||||
103
test/gtx/gtx_matrix_factorisation.cpp
Normal file
103
test/gtx/gtx_matrix_factorisation.cpp
Normal file
@ -0,0 +1,103 @@
|
||||
#define GLM_ENABLE_EXPERIMENTAL
|
||||
#include <glm/gtx/matrix_factorisation.hpp>
|
||||
|
||||
const double 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) {
|
||||
matType<std::min(C, R), R, T, P> q(-999);
|
||||
matType<C, std::min(C, R), T, P> r(-999);
|
||||
|
||||
glm::qr_decompose(q, r, m);
|
||||
|
||||
//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++) {
|
||||
if (abs(err[i][j]) > epsilon) return 1;
|
||||
}
|
||||
}
|
||||
|
||||
//Test if the columns of q are orthonormal
|
||||
for (glm::length_t i = 0; i < std::min(C, R); i++) {
|
||||
if ((length(q[i]) - 1) > epsilon) return 1;
|
||||
|
||||
for (glm::length_t j = 0; j<i; j++) {
|
||||
if (abs(dot(q[i], q[j])) > epsilon) return 1;
|
||||
}
|
||||
}
|
||||
|
||||
//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 < std::min(C, R); j++) {
|
||||
if (r[i][j] != 0) return 1;
|
||||
}
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
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) {
|
||||
matType<C, std::min(C, R), T, P> q(-999);
|
||||
matType<std::min(C, R), R, T, P> r(-999);
|
||||
|
||||
glm::rq_decompose(r, q, m);
|
||||
|
||||
//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++) {
|
||||
if (abs(err[i][j]) > epsilon) return 1;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
//Test if the rows of q are orthonormal
|
||||
matType<std::min(C, R), C, T, P> tq = transpose(q);
|
||||
|
||||
for (glm::length_t i = 0; i < std::min(C, R); i++) {
|
||||
if ((length(tq[i]) - 1) > epsilon) return 1;
|
||||
|
||||
for (glm::length_t j = 0; j<i; j++) {
|
||||
if (abs(dot(tq[i], tq[j])) > epsilon) return 1;
|
||||
}
|
||||
}
|
||||
|
||||
//Test if the matrix r is upper triangular
|
||||
for (glm::length_t i = 0; i < std::min(C, R); i++) {
|
||||
for (glm::length_t j = i + 1; j < R; j++) {
|
||||
if (r[i][j] != 0) return 1;
|
||||
}
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
|
||||
//Test QR square
|
||||
if(test_qr(glm::dmat3(12, 6, -4, -51, 167, 24, 4, -68, -41))) return 1;
|
||||
|
||||
//Test RQ square
|
||||
if (test_rq(glm::dmat3(12, 6, -4, -51, 167, 24, 4, -68, -41))) return 1;
|
||||
|
||||
//Test QR triangular 1
|
||||
if (test_qr(glm::dmat3x4(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15))) return 1;
|
||||
|
||||
//Test QR triangular 2
|
||||
if (test_qr(glm::dmat4x3(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15))) return 1;
|
||||
|
||||
//Test RQ triangular 1 : Fails at the triangular test
|
||||
//if (test_rq(glm::dmat3x4(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15))) return 1;
|
||||
|
||||
//Test QR triangular 2
|
||||
if (test_rq(glm::dmat4x3(12, 6, -4, -51, 167, 24, 4, -68, -41, 7, 2, 15))) return 1;
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user