From 8159718fafc57cfbc9da7f706e599a91caccfa42 Mon Sep 17 00:00:00 2001 From: Bastien Montagne Date: Fri, 9 Oct 2015 20:55:15 +0200 Subject: [PATCH] BLI: add SVD solver for mat3 (using eigen3). --- extern/Eigen3/CMakeLists.txt | 2 + extern/Eigen3/eigen3_capi.h | 1 + extern/Eigen3/intern/svd.cc | 72 ++++++++++++++++++++ extern/Eigen3/intern/svd.h | 40 +++++++++++ source/blender/blenlib/BLI_math_solvers.h | 1 + source/blender/blenlib/intern/math_solvers.c | 13 ++++ 6 files changed, 129 insertions(+) create mode 100644 extern/Eigen3/intern/svd.cc create mode 100644 extern/Eigen3/intern/svd.h diff --git a/extern/Eigen3/CMakeLists.txt b/extern/Eigen3/CMakeLists.txt index 9bbfc9a3670..e3b63881aca 100644 --- a/extern/Eigen3/CMakeLists.txt +++ b/extern/Eigen3/CMakeLists.txt @@ -34,8 +34,10 @@ set(SRC eigen3_capi.h intern/eigenvalues.cc + intern/svd.cc intern/eigenvalues.h + intern/svd.h ) blender_add_lib(extern_eigen3 "${SRC}" "${INC}" "${INC_SYS}") diff --git a/extern/Eigen3/eigen3_capi.h b/extern/Eigen3/eigen3_capi.h index 16f223793a9..f8a7b3cbb77 100644 --- a/extern/Eigen3/eigen3_capi.h +++ b/extern/Eigen3/eigen3_capi.h @@ -28,5 +28,6 @@ #define __EIGEN3_C_API_H__ #include "intern/eigenvalues.h" +#include "intern/svd.h" #endif /* __EIGEN3_C_API_H__ */ diff --git a/extern/Eigen3/intern/svd.cc b/extern/Eigen3/intern/svd.cc new file mode 100644 index 00000000000..e39a8261edb --- /dev/null +++ b/extern/Eigen3/intern/svd.cc @@ -0,0 +1,72 @@ +/* + * ***** BEGIN GPL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * The Original Code is Copyright (C) 2015 Blender Foundation. + * All rights reserved. + * + * Contributor(s): Blender Foundation, + * Bastien Montagne + * + * ***** END GPL LICENSE BLOCK ***** + */ + +#ifndef __EIGEN3_SVD_C_API_CC__ +#define __EIGEN3_SVD_C_API_CC__ + +/* Eigen gives annoying huge amount of warnings here, silence them! */ +#ifdef __GNUC__ +# pragma GCC diagnostic ignored "-Wlogical-op" +#endif + +#include +#include + +#include "svd.h" + +using Eigen::JacobiSVD; + +using Eigen::NoQRPreconditioner; + +using Eigen::ComputeThinU; +using Eigen::ComputeThinV; + +using Eigen::MatrixXf; +using Eigen::VectorXf; +using Eigen::Map; + +void EG3_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V) +{ + /* Since our matrix is squared, we can use thinU/V. */ + unsigned int flags = (r_U ? ComputeThinU : 0) | (r_V ? ComputeThinV : 0); + + /* Blender and Eigen matrices are both column-major. */ + JacobiSVD svd(Map((float *)matrix, size, size), flags); + + if (r_U) { + Map(r_U, size, size) = svd.matrixU(); + } + + if (r_S) { + Map(r_S, size) = svd.singularValues(); + } + + if (r_V) { + Map(r_V, size, size) = svd.matrixV(); + } +} + +#endif /* __EIGEN3_SVD_C_API_CC__ */ diff --git a/extern/Eigen3/intern/svd.h b/extern/Eigen3/intern/svd.h new file mode 100644 index 00000000000..0ac51108977 --- /dev/null +++ b/extern/Eigen3/intern/svd.h @@ -0,0 +1,40 @@ +/* + * ***** BEGIN GPL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * The Original Code is Copyright (C) 2015 Blender Foundation. + * All rights reserved. + * + * Contributor(s): Blender Foundation, + * Bastien Montagne + * + * ***** END GPL LICENSE BLOCK ***** + */ + +#ifndef __EIGEN3_SVD_C_API_H__ +#define __EIGEN3_SVD_C_API_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +void EG3_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V); + +#ifdef __cplusplus +} +#endif + +#endif /* __EIGEN3_SVD_C_API_H__ */ diff --git a/source/blender/blenlib/BLI_math_solvers.h b/source/blender/blenlib/BLI_math_solvers.h index ec9ba5538e2..810c84cc830 100644 --- a/source/blender/blenlib/BLI_math_solvers.h +++ b/source/blender/blenlib/BLI_math_solvers.h @@ -46,6 +46,7 @@ extern "C" { bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], float r_eigen_values[3], float r_eigen_vectors[3][3]); +void BLI_svd_m3(const float m3[3][3], float r_U[3][3], float r_S[], float r_V[3][3]); /**************************** Inline Definitions ******************************/ #if 0 /* None so far. */ diff --git a/source/blender/blenlib/intern/math_solvers.c b/source/blender/blenlib/intern/math_solvers.c index 2f962714c8c..d1dad9a6269 100644 --- a/source/blender/blenlib/intern/math_solvers.c +++ b/source/blender/blenlib/intern/math_solvers.c @@ -59,3 +59,16 @@ bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], float r_eigen_values[3 return EG3_self_adjoint_eigen_solve(3, (const float *)m3, r_eigen_values, (float *)r_eigen_vectors); } + +/** + * \brief Compute the SVD (Singular Values Decomposition) of given 3D matrix (m3 = USV*). + * + * \param m3 the matrix to decompose. + * \return r_U the computed left singular vector of \a m3 (NULL if not needed). + * \return r_S the computed singular values of \a m3 (NULL if not needed). + * \return r_V the computed right singular vector of \a m3 (NULL if not needed). + */ +void BLI_svd_m3(const float m3[3][3], float r_U[3][3], float r_S[3], float r_V[3][3]) +{ + EG3_svd_square_matrix(3, (const float *)m3, (float *)r_U, (float *)r_S, (float *)r_V); +}