Commit ed2f3bb2 authored by Stefan Philips's avatar Stefan Philips

[CHANGE ] Template parameter for pseudoInverse and computeSVD

parent 13794c18
......@@ -24,8 +24,6 @@
#include <vector>
#include <Eigen/SVD>
#include "../WAssert.h"
#include "../WLimits.h"
......@@ -292,34 +290,3 @@ bool linearIndependent( const WVector3d& u, const WVector3d& v )
}
return true;
}
void computeSVD( const WMatrix<double>& A,
WMatrix<double>& U,
WMatrix<double>& V,
WValue<double>& S )
{
Eigen::MatrixXd eigenA( A );
Eigen::JacobiSVD<Eigen::MatrixXd> svd( eigenA, Eigen::ComputeFullU | Eigen::ComputeFullV );
U = WMatrix<double>( svd.matrixU() );
V = WMatrix<double>( svd.matrixV() );
S = WValue<double>( svd.singularValues() );
}
WMatrix<double> pseudoInverse( const WMatrix<double>& input )
{
// calc pseudo inverse
WMatrix<double> U( input.getNbRows(), input.getNbCols() );
WMatrix<double> V( input.getNbCols(), input.getNbCols() );
WValue<double> Svec( input.size() );
computeSVD( input, U, V, Svec );
// create diagonal matrix S
WMatrix<double> S( input.getNbCols(), input.getNbCols() );
S.setZero();
for( size_t i = 0; i < Svec.size() && i < S.getNbRows() && i < S.getNbCols(); i++ )
{
S( i, i ) = ( Svec[ i ] == 0.0 ) ? 0.0 : 1.0 / Svec[ i ];
}
return WMatrix<double>( V*S*U.transposed() );
}
......@@ -25,6 +25,7 @@
#ifndef WLINEARALGEBRAFUNCTIONS_H
#define WLINEARALGEBRAFUNCTIONS_H
#include <Eigen/SVD>
#include "WMatrix.h"
#include "linearAlgebra/WLinearAlgebra.h"
......@@ -92,22 +93,60 @@ bool linearIndependent( const WVector3d& u, const WVector3d& v );
*
* A=U*S*V^T
*
* \tparam T The data type.
* \param A Input Matrix
* \param U Output Matrix
* \param S Output of the entries in the diagonal matrix
* \param V Output Matrix
*
*/
void computeSVD( const WMatrix<double>& A, WMatrix<double>& U, WMatrix<double>& V, WValue<double>& S );
template< typename T >
void computeSVD( const WMatrix< T >& A, WMatrix< T >& U, WMatrix< T >& V, WValue< T >& S );
/**
* Calculates for a matrix the pseudo inverse.
*
* \tparam T The data type.
* \param input Matrix to invert
*
* \return Inverted Matrix
*
*/
WMatrix<double> pseudoInverse( const WMatrix<double>& input );
template< typename T >
WMatrix< T > pseudoInverse( const WMatrix< T >& input );
template< typename T >
void computeSVD( const WMatrix< T >& A,
WMatrix< T >& U,
WMatrix< T >& V,
WValue< T >& S )
{
Eigen::Matrix< T, -1, -1 > eigenA( A );
Eigen::JacobiSVD< Eigen::Matrix< T, -1, -1 > > svd( eigenA, Eigen::ComputeFullU | Eigen::ComputeFullV );
U = WMatrix< T >( svd.matrixU() );
V = WMatrix< T >( svd.matrixV() );
S = WValue< T >( svd.singularValues() );
}
template< typename T >
WMatrix< T > pseudoInverse( const WMatrix< T >& input )
{
// calc pseudo inverse
WMatrix< T > U( input.getNbRows(), input.getNbCols() );
WMatrix< T > V( input.getNbCols(), input.getNbCols() );
WValue< T > Svec( input.size() );
computeSVD( input, U, V, Svec );
// create diagonal matrix S
WMatrix< T > S( input.getNbCols(), input.getNbCols() );
S.setZero();
for( size_t i = 0; i < Svec.size() && i < S.getNbRows() && i < S.getNbCols(); i++ )
{
S( i, i ) = ( Svec[ i ] == 0.0 ) ? 0.0 : 1.0 / Svec[ i ];
}
return WMatrix< T >( V*S*U.transposed() );
}
#endif // WLINEARALGEBRAFUNCTIONS_H
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment