Commit d531d348 authored by Sebastian Eichelbaum's avatar Sebastian Eichelbaum
Browse files

[MERGE]

parents e2a3af19 dcdcf30d
......@@ -24,17 +24,16 @@
#include <vector>
#include "../../ext/Eigen/SVD"
#include "../WAssert.h"
#include "WLinearAlgebraFunctions.h"
#include "WMatrix.h"
#include "WMatrix4x4.h"
#include "WValue.h"
#include "WVector3D.h"
#ifdef OW_USE_OSSIM
#include "WOSSIMHelper.h"
#endif
WVector3D multMatrixWithVector3D( WMatrix<double> mat, WVector3D vec )
{
WVector3D result;
......@@ -301,32 +300,32 @@ 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 )
void computeSVD( const WMatrix_2& A,
WMatrix_2& U,
WMatrix_2& V,
WVector_2& S )
{
#ifdef OW_USE_OSSIM
WOSSIMHelper::computeSVD( A, U, V, S );
#else
(void) A; (void) U; (void) V; (void) S; // NOLINT to prevent "unused variable" warnings
WAssert( false, "OpenWalnut must be compiled with OSSIM to support SVD." );
#endif
Eigen::JacobiSVD<WMatrix_2> svd( A, Eigen::ComputeFullU | Eigen::ComputeFullV );
U = svd.matrixU();
V = svd.matrixV();
S = svd.singularValues();
}
WMatrix<double> pseudoInverse( const WMatrix<double>& input )
WMatrix_2 pseudoInverse( const WMatrix_2& input )
{
// calc pseudo inverse
WMatrix< double > U( input.getNbRows(), input.getNbCols() );
WMatrix< double > V( input.getNbCols(), input.getNbCols() );
WValue< double > Svec( input.getNbCols() );
WMatrix_2 U( input.rows(), input.cols() );
WMatrix_2 V( input.cols(), input.cols() );
WVector_2 Svec( input.cols() );
computeSVD( input, U, V, Svec );
// create diagonal matrix S
WMatrix< double > S( input.getNbCols(), input.getNbCols() );
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 ];
WMatrix_2 S( input.cols(), input.cols() );
S.setZero();
for ( int i = 0; i < Svec.size() && i < S.rows() && i < S.cols(); i++ )
{
S( i, i ) = ( Svec[ i ] == 0.0 ) ? 0.0 : 1.0 / Svec[ i ];
}
return WMatrix< double >( V*S*U.transposed() );
return WMatrix_2( V*S*U.transpose() );
}
......@@ -26,6 +26,7 @@
#define WLINEARALGEBRAFUNCTIONS_H
#include "../WExportCommon.h"
#include "WMatrix.h"
namespace osg
{
......@@ -115,8 +116,7 @@ bool OWCOMMON_EXPORT linearIndependent( const WVector3D& u, const WVector3D& v )
*
* \note This function needs the OSSIM library.
*/
void OWCOMMON_EXPORT computeSVD( const WMatrix< double >& A, WMatrix< double >& U,
WMatrix< double >& V, WValue< double >& S );
void OWCOMMON_EXPORT computeSVD( const WMatrix_2& A, WMatrix_2& U, WMatrix_2& V, WVector_2& S );
/**
* Calculates for a matrix the pseudo inverse.
......@@ -126,6 +126,6 @@ void OWCOMMON_EXPORT computeSVD( const WMatrix< double >& A, WMatrix< double >&
* \return Inverted Matrix
*
*/
WMatrix< double > OWCOMMON_EXPORT pseudoInverse( const WMatrix<double>& input );
WMatrix_2 OWCOMMON_EXPORT pseudoInverse( const WMatrix_2& input );
#endif // WLINEARALGEBRAFUNCTIONS_H
......@@ -72,6 +72,15 @@ typedef Eigen::Matrix< double, 4, 4 > WMatrix4x4_2;
*/
typedef Eigen::MatrixXd WMatrix_2;
/**
* A complex double matrix of dynamic size. Heap-allocated.
* If you want to access coefficients using the operator( size_t, size_t ), the first parameter is still the row index, starting with 0.
*
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html
* \see http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html
*/
typedef Eigen::MatrixXcd WMatrixComplex_2;
/**
* Matrix template class with variable number of rows and columns.
* The access function are row-major, which means that the rows
......
......@@ -34,15 +34,16 @@
#include "WSymmetricSphericalHarmonic.h"
#include "WTensorSym.h"
#include "WUnitSphereCoordinates.h"
#include "WValue.h"
#include "WVector3D.h"
WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic() :
m_order( 0 ),
m_SHCoefficients( 0 )
m_order( 0 )
// ,
// m_SHCoefficients( 0 )
{
}
WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic( const WValue<double>& SHCoefficients ) :
WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic( const WVector_2& SHCoefficients ) :
m_SHCoefficients( SHCoefficients )
{
// calculate order from SHCoefficients.size:
......@@ -57,28 +58,28 @@ WSymmetricSphericalHarmonic::~WSymmetricSphericalHarmonic()
double WSymmetricSphericalHarmonic::getValue( double theta, double phi ) const
{
double result = 0.0;
int j = 0;
const double rootOf2 = std::sqrt( 2.0 );
for ( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
{
// m = 1 .. k
for ( int m = 1; m <= k; m++ )
double result = 0.0;
int j = 0;
const double rootOf2 = std::sqrt( 2.0 );
for ( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
{
j = ( k*k+k+2 ) / 2 + m;
result += m_SHCoefficients[ j-1 ] * rootOf2 *
static_cast<double>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
}
// m = 0
result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
// m = -k .. -1
for ( int m = -k; m < 0; m++ )
{
j = ( k*k+k+2 ) / 2 + m;
result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
// m = 1 .. k
for ( int m = 1; m <= k; m++ )
{
j = ( k*k+k+2 ) / 2 + m;
result += m_SHCoefficients[ j-1 ] * rootOf2 *
static_cast<double>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
}
// m = 0
result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
// m = -k .. -1
for ( int m = -k; m < 0; m++ )
{
j = ( k*k+k+2 ) / 2 + m;
result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
}
}
}
return result;
return result;
}
double WSymmetricSphericalHarmonic::getValue( const WUnitSphereCoordinates& coordinates ) const
......@@ -86,14 +87,14 @@ double WSymmetricSphericalHarmonic::getValue( const WUnitSphereCoordinates& coor
return getValue( coordinates.getTheta(), coordinates.getPhi() );
}
const WValue<double>& WSymmetricSphericalHarmonic::getCoefficients() const
const WVector_2& WSymmetricSphericalHarmonic::getCoefficients() const
{
return m_SHCoefficients;
}
WValue< double > WSymmetricSphericalHarmonic::getCoefficientsSchultz() const
WVector_2 WSymmetricSphericalHarmonic::getCoefficientsSchultz() const
{
WValue< double > res( m_SHCoefficients.size() );
WVector_2 res( m_SHCoefficients.size() );
size_t k = 0;
for ( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
{
......@@ -114,9 +115,9 @@ WValue< double > WSymmetricSphericalHarmonic::getCoefficientsSchultz() const
return res;
}
WValue< std::complex< double > > WSymmetricSphericalHarmonic::getCoefficientsComplex() const
WVectorComplex_2 WSymmetricSphericalHarmonic::getCoefficientsComplex() const
{
WValue< std::complex< double > > res( m_SHCoefficients.size() );
WVectorComplex_2 res( m_SHCoefficients.size() );
size_t k = 0;
double r = 1.0 / sqrt( 2.0 );
std::complex< double > i( 0.0, -1.0 );
......@@ -187,30 +188,30 @@ double WSymmetricSphericalHarmonic::calcGFA( std::vector< WUnitSphereCoordinates
return gfa;
}
double WSymmetricSphericalHarmonic::calcGFA( WMatrix< double > const& B ) const
double WSymmetricSphericalHarmonic::calcGFA( const WMatrix_2& B ) const
{
// sh coeffs
WMatrix< double > s( B.getNbCols(), 1 );
WAssert( B.getNbCols() == m_SHCoefficients.size(), "" );
for( std::size_t k = 0; k < B.getNbCols(); ++k )
WMatrix_2 s( B.cols(), 1 );
WAssert( B.cols() == m_SHCoefficients.size(), "" );
for( int k = 0; k < B.cols(); ++k )
{
s( k, 0 ) = m_SHCoefficients[ k ];
}
s = B * s;
WAssert( s.getNbRows() == B.getNbRows(), "" );
WAssert( s.getNbCols() == 1u, "" );
WAssert( s.rows() == B.rows(), "" );
WAssert( s.cols() == 1u, "" );
double n = static_cast< double >( s.getNbRows() );
double n = static_cast< double >( s.rows() );
double d = 0.0;
double gfa = 0.0;
double mean = 0.0;
for( std::size_t i = 0; i < s.getNbRows(); ++i )
for( int i = 0; i < s.rows(); ++i )
{
mean += s( i, 0 );
}
mean /= n;
for( std::size_t i = 0; i < s.getNbRows(); ++i )
for( int i = 0; i < s.rows(); ++i )
{
double f = s( i, 0 );
// we use gfa as a temporary here
......@@ -238,12 +239,12 @@ double WSymmetricSphericalHarmonic::calcGFA( WMatrix< double > const& B ) const
return gfa;
}
void WSymmetricSphericalHarmonic::applyFunkRadonTransformation( WMatrix< double > const& frtMat )
void WSymmetricSphericalHarmonic::applyFunkRadonTransformation( const WMatrix_2& frtMat )
{
WAssert( frtMat.getNbCols() == m_SHCoefficients.size(), "" );
WAssert( frtMat.getNbRows() == m_SHCoefficients.size(), "" );
WAssert( frtMat.cols() == m_SHCoefficients.size(), "" );
WAssert( frtMat.rows() == m_SHCoefficients.size(), "" );
// Funk-Radon-Transformation as in Descoteaux's thesis
for ( size_t j = 0; j < m_SHCoefficients.size(); j++ )
for ( int j = 0; j < m_SHCoefficients.size(); j++ )
{
m_SHCoefficients[ j ] *= frtMat( j, j );
}
......@@ -254,7 +255,7 @@ size_t WSymmetricSphericalHarmonic::getOrder() const
return m_order;
}
WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WVector3D >& orientations,
WMatrix_2 WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WVector3D >& orientations,
int order,
double lambda,
bool withFRT )
......@@ -268,33 +269,34 @@ WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vect
return WSymmetricSphericalHarmonic::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
}
WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
int order,
double lambda,
bool withFRT )
WMatrix_2 WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
int order,
double lambda,
bool withFRT )
{
WMatrix<double> B( WSymmetricSphericalHarmonic::calcBaseMatrix( orientations, order ) );
WMatrix<double> Bt( B.transposed() );
WMatrix<double> result( Bt*B );
WMatrix_2 B( WSymmetricSphericalHarmonic::calcBaseMatrix( orientations, order ) );
WMatrix_2 Bt( B.transpose() );
WMatrix_2 result( Bt*B );
if ( lambda != 0.0 )
{
WMatrix<double> L( WSymmetricSphericalHarmonic::calcSmoothingMatrix( order ) );
WMatrix_2 L( WSymmetricSphericalHarmonic::calcSmoothingMatrix( order ) );
result += lambda*L;
}
result = pseudoInverse( result )*Bt;
// result *= Bt;
if ( withFRT )
{
WMatrix<double> P( WSymmetricSphericalHarmonic::calcFRTMatrix( order ) );
WMatrix_2 P( WSymmetricSphericalHarmonic::calcFRTMatrix( order ) );
return ( P * result );
}
return result;
}
WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
WMatrix_2 WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
int order )
{
WMatrix<double> B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
WMatrix_2 B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
// calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
int j = 0;
......@@ -325,10 +327,10 @@ WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector<
return B;
}
WMatrix< std::complex< double > >
WMatrixComplex_2
WSymmetricSphericalHarmonic::calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations, int order )
{
WMatrix< std::complex< double > > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
WMatrixComplex_2 B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
for( std::size_t row = 0; row < orientations.size(); row++ )
{
......@@ -355,11 +357,11 @@ WSymmetricSphericalHarmonic::calcComplexBaseMatrix( std::vector< WUnitSphereCoor
return B;
}
WMatrix<double> WSymmetricSphericalHarmonic::calcSmoothingMatrix( size_t order )
WMatrix_2 WSymmetricSphericalHarmonic::calcSmoothingMatrix( size_t order )
{
size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
std::size_t i = 0;
WMatrix<double> L( R, R );
WMatrix_2 L( R, R );
for ( size_t k = 0; k <= order; k += 2 )
{
for ( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
......@@ -371,11 +373,11 @@ WMatrix<double> WSymmetricSphericalHarmonic::calcSmoothingMatrix( size_t order )
return L;
}
WMatrix<double> WSymmetricSphericalHarmonic::calcFRTMatrix( size_t order )
WMatrix_2 WSymmetricSphericalHarmonic::calcFRTMatrix( size_t order )
{
size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
std::size_t i = 0;
WMatrix< double > result( R, R );
WMatrix_2 result( R, R );
for ( size_t k = 0; k <= order; k += 2 )
{
double h = 2.0 * piDouble * static_cast< double >( std::pow( static_cast< double >( -1 ), static_cast< double >( k / 2 ) ) ) *
......@@ -391,8 +393,8 @@ WMatrix<double> WSymmetricSphericalHarmonic::calcFRTMatrix( size_t order )
#ifdef OW_USE_OSSIM
WMatrix< double > WSymmetricSphericalHarmonic::calcSHToTensorSymMatrix( std::size_t order,
const std::vector< WUnitSphereCoordinates >& orientations )
WMatrix_2 WSymmetricSphericalHarmonic::calcSHToTensorSymMatrix( std::size_t order,
const std::vector< WUnitSphereCoordinates >& orientations )
{
std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
WPrecondEquals( order % 2, 0u );
......@@ -409,7 +411,7 @@ WMatrix< double > WSymmetricSphericalHarmonic::calcSHToTensorSymMatrix( std::siz
std::vector< std::size_t > indices( order, 0 );
// calc tensor evaluation matrix
WMatrix< double > TEMat( numElements, numElements );
WMatrix_2 TEMat( numElements, numElements );
for( std::size_t j = 0; j < numElements; ++j ) // j is the 'permutation index'
{
// stores how often each value is represented in the index array
......@@ -438,9 +440,20 @@ WMatrix< double > WSymmetricSphericalHarmonic::calcSHToTensorSymMatrix( std::siz
// we do not want more orientations than nessessary
std::vector< WUnitSphereCoordinates > ori2( orientations.begin(), orientations.begin() + numElements );
WMatrix< double > p = pseudoInverse( TEMat );
WMatrix_2 p = pseudoInverse( TEMat );
return p * calcBaseMatrix( ori2, order );
}
#endif // OW_USE_OSSIM
void WSymmetricSphericalHarmonic::normalize()
{
double scale = 0.0;
if ( m_SHCoefficients.size() > 0 )
scale = std::sqrt( 4.0 * piDouble ) * m_SHCoefficients[0];
for ( int i = 0; i < m_SHCoefficients.size(); i++ )
{
m_SHCoefficients[ i ] /= scale;
}
}
......@@ -31,7 +31,7 @@
#include "WMath.h"
#include "WMatrix.h"
#include "WUnitSphereCoordinates.h"
#include "WValue.h"
#include "WVector3D.h"
/**
* Class for symmetric spherical harmonics
......@@ -50,7 +50,7 @@ public:
* Constructor.
* \param SHCoefficients the initial coefficients (stored like in the mentioned Descoteaux paper).
*/
explicit WSymmetricSphericalHarmonic( const WValue<double>& SHCoefficients );
explicit WSymmetricSphericalHarmonic( const WVector_2& SHCoefficients );
/**
* Destructor.
......@@ -73,17 +73,17 @@ public:
/**
* Returns the used coefficients (stored like in the mentioned 2007 Descoteaux paper).
*/
const WValue<double>& getCoefficients() const;
const WVector_2& getCoefficients() const;
/**
* Returns the coefficients for Schultz' SH base.
*/
WValue< double > getCoefficientsSchultz() const;
WVector_2 getCoefficientsSchultz() const;
/**
* Returns the coefficients for the complex base.
*/
WValue< std::complex< double > > getCoefficientsComplex() const;
WVectorComplex_2 getCoefficientsComplex() const;
/**
* Applies the Funk-Radon-Transformation. This is faster than matrix multiplication.
......@@ -91,7 +91,7 @@ public:
*
* \param frtMat the frt matrix as calculated by calcFRTMatrix()
*/
void applyFunkRadonTransformation( WMatrix< double > const& frtMat );
void applyFunkRadonTransformation( const WMatrix_2& frtMat );
/**
* Return the order of the spherical harmonic.
......@@ -125,7 +125,7 @@ public:
*
* \return The generalized fractional anisotropy.
*/
double calcGFA( WMatrix< double > const& B ) const;
double calcGFA( const WMatrix_2& B ) const;
/**
* This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper. The orientations are given as WVector3D.
......@@ -135,10 +135,10 @@ public:
* \param withFRT include the Funk-Radon-Transformation?
* \return Transformation matrix
*/
static WMatrix<double> getSHFittingMatrix( const std::vector< WVector3D >& orientations,
int order,
double lambda,
bool withFRT );
static WMatrix_2 getSHFittingMatrix( const std::vector< WVector3D >& orientations,
int order,
double lambda,
bool withFRT );
/**
* This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper. The orientations are given as WUnitSphereCoordinates .
......@@ -148,7 +148,7 @@ public:
* \param withFRT include the Funk-Radon-Transformation?
* \return Transformation matrix
*/
static WMatrix<double> getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
static WMatrix_2 getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
int order,
double lambda,
bool withFRT );
......@@ -159,7 +159,7 @@ public:
* \param order The order of the spherical harmonics intended to create
* \return The base Matrix B
*/
static WMatrix<double> calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations, int order );
static WMatrix_2 calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations, int order );
/**
* Calculates the base matrix B for the complex spherical harmonics.
......@@ -167,7 +167,7 @@ public:
* \param order The order of the spherical harmonics intended to create
* \return The base Matrix B
*/
static WMatrix< std::complex< double > > calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations,
static WMatrixComplex_2 calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations,
int order );
/**
......@@ -175,14 +175,14 @@ public:
* \param order The order of the spherical harmonic
* \return The smoothing matrix L
*/
static WMatrix<double> calcSmoothingMatrix( size_t order );
static WMatrix_2 calcSmoothingMatrix( size_t order );
/**
* Calculates the Funk-Radon-Transformation-Matrix P from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging"
* \param order The order of the spherical harmonic
* \return The Funk-Radon-Matrix P
*/
static WMatrix<double> calcFRTMatrix( size_t order );
static WMatrix_2 calcFRTMatrix( size_t order );
#ifdef OW_USE_OSSIM
/**
......@@ -191,8 +191,9 @@ public:
* \param order The order of the symmetric tensor.
* \param orientations A vector of at least (orderTensor+1) * (orderTensor+2) / 2 orientations.
*/
static WMatrix< double > calcSHToTensorSymMatrix( std::size_t order, const std::vector< WUnitSphereCoordinates >& orientations );
static WMatrix_2 calcSHToTensorSymMatrix( std::size_t order, const std::vector< WUnitSphereCoordinates >& orientations );
#endif // OW_USE_OSSIM
void normalize();
protected:
......@@ -201,7 +202,7 @@ private:
size_t m_order;
/** coefficients of the spherical harmonic */
WValue<double> m_SHCoefficients;
WVector_2 m_SHCoefficients;
};
#endif // WSYMMETRICSPHERICALHARMONIC_H
......@@ -31,6 +31,7 @@
#include "../WAssert.h"
#include "../WStringUtils.h"
#include "WVector3D.h"
/**
* Base class for all higher level values like tensors, vectors, matrices and so on.
......@@ -77,6 +78,19 @@ public:
}
}
/**
* Create a WValue from the given WVector_2.
* \param newValues The WVector_2 with the values..
*/
explicit WValue( const WVector_2& newValues )
: m_components( static_cast< std::size_t >( newValues.size() ) )
{
for ( std::size_t i = 0; i < m_components.size(); ++i )
{
m_components[ i ] = static_cast< T >( newValues( i ) );
}
}
/**
* Get number of components the value consists of.
*/
......@@ -320,6 +334,19 @@ public:
m_components.resize( size );
}
/**
* Returns this WValue as WVector_2.
*/
WVector_2 toWVector()
{
WVector_2 result( m_components.size() );
for ( size_t i = 0; i < m_components.size(); ++i )
{
result( i ) = static_cast<double>( m_components[ i ] );
}
return result;
}
protected:
private:
/**
......
......@@ -70,6 +70,14 @@ typedef WVector3D_2 WPosition_2;
*/
typedef Eigen::Matrix< double, Eigen::Dynamic, 1 > WVector_2;