Commit 40d11e56 authored by Stefan Philips's avatar Stefan Philips
Browse files

[CHANGE] Adapt WSymmetricSphericalHarmonic to Eigen vector and matrix

parent 87e22617
......@@ -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++ )
{
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++ )
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 * 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,31 @@ 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, "" );
// TODO(philips): is this a problem?
s = B * WMatrix_2( s );
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 +240,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 +256,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 +270,33 @@ WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vect
return WSymmetricSphericalHarmonic::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
}
WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
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 ( size_t 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
......@@ -27,10 +27,20 @@
#include <vector>
// the following block is only for the last termporarly test
// #include <boost/nondet_random.hpp>
// #include <boost/random.hpp>
// #include <boost/random/normal_distribution.hpp>
// #include <boost/random/uniform_real.hpp>
// #include <boost/random/linear_congruential.hpp>
// #include <boost/random/uniform_real.hpp>
// #include <boost/random/variate_generator.hpp>
#include <cxxtest/TestSuite.h>
#include "../WMatrix.h"
#include "../WValue.h"
#include "../WVector3D.h"
#include "../WGeometryFunctions.h"
#include "../WSymmetricSphericalHarmonic.h"
......@@ -51,8 +61,8 @@ public:
*/
void testCalcFRTMatrix( void )
{
WMatrix<double> result( WSymmetricSphericalHarmonic::calcFRTMatrix( 4 ) );
WMatrix<double> reference( 15, 15 );
WMatrix_2 result( WSymmetricSphericalHarmonic::calcFRTMatrix( 4 ) );
WMatrix_2 reference( 15, 15 );
// j 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
// lj 0 2 2 2 2 2 4 4 4 4 4 4 4 4 4
reference( 0, 0 ) = 2.0 * piDouble;
......@@ -73,8 +83,8 @@ public:
*/
void testCalcSmoothingMatrix( void )
{
WMatrix<double> result( WSymmetricSphericalHarmonic::calcSmoothingMatrix( 4 ) );
WMatrix<double> reference( 15, 15 );
WMatrix_2 result( WSymmetricSphericalHarmonic::calcSmoothingMatrix( 4 ) );
WMatrix_2 reference( 15, 15 );
// j 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
// lj 0 2 2 2 2 2 4 4 4 4 4 4 4 4 4
reference( 0, 0 ) = 0.0;
......@@ -96,7 +106,7 @@ public:
void testCalcSHtoTensorMatrix()
{
#ifdef OW_USE_OSSIM
WValue< double > w( 6 );
WVector_2 w( 6 );
for( int i = 0; i < 6; ++i )
{
w[ i ] = exp( i / 6.0 );
......@@ -114,11 +124,13 @@ public:
orientations.push_back( WUnitSphereCoordinates( WVector3D( 0.0, 0.0, -4.0 ).normalized() ) );
orientations.push_back( WUnitSphereCoordinates( WVector3D( 0.0, 4.0, 1.0 ).normalized() ) );
WMatrix< double > SHToTensor = WSymmetricSphericalHarmonic::calcSHToTensorSymMatrix( 2, orientations );
WTensorSym< 2, 3, double > t( SHToTensor * w );
WMatrix_2 SHToTensor = WSymmetricSphericalHarmonic::calcSHToTensorSymMatrix( 2, orientations );
// TODO(all): remove the WValue from the following line, when WTensorSym supports WVector_2
WTensorSym< 2, 3, double > t( WValue<double>( SHToTensor * w ) );
for( std::vector< WUnitSphereCoordinates >::iterator it = orientations.begin();
it != orientations.end(); ++it )
it != orientations.end();
++it )
{
TS_ASSERT_DELTA( i.getValue( *it ), evaluateSphericalFunction( t, it->getEuclidean() ), 0.001 );
}
......@@ -146,18 +158,18 @@ public:
}
grad.clear();
WValue< double > values( 15 );
WVector_2 values( 15 );
for( std::size_t i = 0; i < 15; ++i )
{
values[ i ] = i / 15.0;
}
WSymmetricSphericalHarmonic sh( values );
WValue< std::complex< double > > values2 = sh.getCoefficientsComplex();
WMatrix< std::complex< double > > complexBaseMatrix = WSymmetricSphericalHarmonic::calcComplexBaseMatrix( orientations, 4 );
WVectorComplex_2 values2 = sh.getCoefficientsComplex();
WValue< std::complex< double > > res = complexBaseMatrix * values2;
WMatrixComplex_2 complexBaseMatrix = WSymmetricSphericalHarmonic::calcComplexBaseMatrix( orientations, 4 );
WVectorComplex_2 res = complexBaseMatrix * values2;
for( std::size_t k = 0; k < orientations.size(); ++k )
{
......
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