Commit 04c0c51f by Alexander Wiebel

[MERGE]

parents 1cf0d972 fa7054fc
......@@ -37,8 +37,6 @@
namespace wmath
{
std::vector<size_t> wmath::WSymmetricSphericalHarmonic::m_lj;
WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic() :
m_order( 0 ),
m_SHCoefficients( 0 )
......@@ -52,7 +50,6 @@ WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic( const WValue<double>&
// - this is done by reversing the R=(l+1)*(l+2)/2 formula from the Descoteaux paper
double q = std::sqrt( 0.25 + 2.0 * static_cast<double>( m_SHCoefficients.size() ) ) - 1.5;
m_order = static_cast<size_t>( q );
if ( WSymmetricSphericalHarmonic::m_lj.size() < m_SHCoefficients.size() ) WSymmetricSphericalHarmonic::calcLj( m_order );
}
WSymmetricSphericalHarmonic::~WSymmetricSphericalHarmonic()
......@@ -83,27 +80,6 @@ double WSymmetricSphericalHarmonic::getValue( double theta, double phi ) const
}
}
return result;
// old implementation
// double result = 0.0;
// int j = 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 ] * std::sqrt( 2.0 ) * 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 ] * std::sqrt( 2.0 ) * boost::math::spherical_harmonic_r( k, m, theta, phi );
// }
// }
// return result;
}
double WSymmetricSphericalHarmonic::getValue( const WUnitSphereCoordinates& coordinates ) const
......@@ -116,29 +92,15 @@ const wmath::WValue<double>& WSymmetricSphericalHarmonic::getCoefficients() cons
return m_SHCoefficients;
}
void WSymmetricSphericalHarmonic::applyFunkRadonTransformation()
void WSymmetricSphericalHarmonic::applyFunkRadonTransformation( wmath::WMatrix< double > const& frtMat )
{
// Funk-Radon-Transformation like in the 2007 Descoteaux paper
// TODO(philips): implement
size_t lj = 0;
// double ljDbl = 0.0;
double factor = 0.0;
for ( size_t j = 0; j < m_SHCoefficients.size(); j++ )
{
// lj = 2 * ( ( j == 1 ) ? 1 : ( j+2 )/4 );
lj = WSymmetricSphericalHarmonic::m_lj[j];
// ljDbl = static_cast<double>( lj );
// std::cerr << "j: " << j << " lj: " << ljDbl << std::endl;
// lj is always even!
factor = 2.0 * wmath::piDouble * static_cast<double>( std::pow( static_cast<double>( -1.0 ), static_cast<double>( lj / 2 ) ) ) *
static_cast<double>( wmath::oddFactorial( ( lj <= 1 ) ? 1 : lj-1 ) ) / static_cast<double>( wmath::evenFactorial( lj ) );
// the following line is the (wrong) first calculation
// static_cast<double>( wmath::oddFactorial( ( lj == 1 ) ? lj-1 : 1 ) ) / static_cast<double>( wmath::evenFactorial( lj ) );
// factor = (double) std::pow( -1, ( lj / 2 ) ) * (double)wmath::oddFaculty( lj ) / (double)wmath::evenFaculty( lj );
// std::cerr << "factor: " << factor << std::endl;
m_SHCoefficients[ j ] = m_SHCoefficients[ j ] * factor;
}
WAssert( frtMat.getNbCols() == m_SHCoefficients.size(), "" );
WAssert( frtMat.getNbRows() == m_SHCoefficients.size(), "" );
// Funk-Radon-Transformation as in Descoteaux's thesis
for ( size_t j = 0; j < m_SHCoefficients.size(); j++ )
{
m_SHCoefficients[ j ] *= frtMat( j, j );
}
}
size_t WSymmetricSphericalHarmonic::getOrder() const
......@@ -146,19 +108,6 @@ size_t WSymmetricSphericalHarmonic::getOrder() const
return m_order;
}
void WSymmetricSphericalHarmonic::calcLj( size_t order )
{
m_lj.clear();
for ( size_t k = 0; k <= order; k += 2 )
{
for ( int m = -k; m <= static_cast<int>( k ); ++m )
{
m_lj.push_back( k );
}
}
}
wmath::WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< wmath::WVector3D >& orientations,
int order,
double lambda,
......@@ -228,67 +177,40 @@ wmath::WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::v
}
}
return B;
// old implementation
// wmath::WMatrix<double> 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;
// for (uint row = 0; row < orientations.size(); row++ )
// {
// const double theta = orientations[row].getTheta();
// const double phi = orientations[row].getPhi();
// for ( int k = 0; k <= order; k+=2 )
// {
// // m = 1 .. k
// for ( int m = 1; m <= k; m++ )
// {
// j = (k*k+k+2)/2+m;
// B( row, j-1 ) = std::sqrt(2.0)*boost::math::spherical_harmonic_i( k, m, theta, phi );
// }
// // m = 0
// B( row, (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;
// B( row, j-1 ) = std::sqrt(2.0)*boost::math::spherical_harmonic_r( k, m, theta, phi );
// }
// }
// }
// return B;
}
wmath::WMatrix<double> WSymmetricSphericalHarmonic::calcSmoothingMatrix( size_t order )
{
std::vector<int> lj;
lj.clear();
for ( size_t k = 0; k <= order; k += 2 )
{
for ( int m = -k; m <= static_cast<int>( k ); ++m )
size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
std::size_t i = 0;
wmath::WMatrix<double> L( R, R );
for ( size_t k = 0; k <= order; k += 2 )
{
lj.push_back( k );
for ( int m = -k; m <= static_cast< int >( k ); ++m )
{
L( i, i ) = static_cast< double > ( k * k * ( k + 1 ) * ( k + 1 ) );
++i;
}
}
}
size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
wmath::WMatrix<double> L( R, R );
for ( size_t i = 0; i < R; i++ ) L( i, i ) = static_cast<double> ( lj[i] * lj[i] * ( lj[i] + 1) * ( lj[i] + 1) );
return L;
return L;
}
wmath::WMatrix<double> WSymmetricSphericalHarmonic::calcFRTMatrix( size_t order )
{
size_t dim = ( order + 1 ) * ( order + 2 ) / 2;
if ( WSymmetricSphericalHarmonic::m_lj.size() < dim ) WSymmetricSphericalHarmonic::calcLj( order );
wmath::WMatrix<double> result( dim, dim );
for ( size_t j = 0; j < dim; j++ )
{
size_t lj = WSymmetricSphericalHarmonic::m_lj[j];
// lj is always even!
result( j, j ) = 2.0 * wmath::piDouble * static_cast<double>( std::pow( static_cast<double>( -1 ), static_cast<double>( lj / 2 ) ) ) *
static_cast<double>( wmath::oddFactorial( ( lj <= 1 ) ? 1 : lj-1 ) ) / static_cast<double>( wmath::evenFactorial( lj ) );
}
return result;
size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
std::size_t i = 0;
wmath::WMatrix< double > result( R, R );
for ( size_t k = 0; k <= order; k += 2 )
{
double h = 2.0 * wmath::piDouble * static_cast< double >( std::pow( static_cast< double >( -1 ), static_cast< double >( k / 2 ) ) ) *
static_cast< double >( wmath::oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast<double>( wmath::evenFactorial( k ) );
for ( int m = -k; m <= static_cast< int >( k ); ++m )
{
result( i, i ) = h;
++i;
}
}
return result;
}
} // Namespace WMath
......@@ -77,28 +77,18 @@ public:
const wmath::WValue<double>& getCoefficients() const;
/**
* Applies the Funk-Radon-Transformation.
* Applies the Funk-Radon-Transformation. This is faster than matrix multiplication.
* ( O(n) instead of O(n²) )
*
* \param frtMat the frt matrix as calculated by calcFRTMatrix()
*/
void applyFunkRadonTransformation();
void applyFunkRadonTransformation( wmath::WMatrix< double > const& frtMat );
/**
* Return the order of the spherical harmonic.
*/
size_t getOrder() const;
/**
* This calculates the l_j array (std::vector) for the given order.
* \param order order
*/
static void calcLj( size_t order );
/**
* The l_j array stores the order to the index.
* index: j={1,2,3,4,5,6,7,8, ...} order l_j={0,2,2,2,2,2,4,4,...}
*/
static std::vector<size_t> m_lj;
/**
* This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper. The orientations are given as wmath::WVector3D.
* \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
......@@ -151,6 +141,7 @@ protected:
private:
/** order of the spherical harmonic */
size_t m_order;
/** coefficients of the spherical harmonic */
WValue<double> m_SHCoefficients;
};
......
......@@ -25,6 +25,7 @@
#ifndef WVALUE_H
#define WVALUE_H
#include <algorithm>
#include <cmath>
#include <vector>
......@@ -272,6 +273,31 @@ public:
return result;
}
/**
* Returns the mean value of all values stored in this WValue.
*/
T mean() const
{
WAssert( !m_components.empty(), "WValue has no entries." );
T sum = 0;
for ( typename std::vector< T >::const_iterator it = m_components.begin(); it != m_components.end(); it++ )
{
sum += ( *it );
}
return ( sum / static_cast< T >( m_components.size() ) );
}
/**
* Returns the median of all values stored in this WValue.
*/
T median() const
{
WAssert( !m_components.empty(), "WValue has no entries. " );
std::vector< T > components( m_components );
std::sort( components.begin(), components.end() );
return components[ components.size() / 2 ];
}
protected:
private:
/**
......
......@@ -643,6 +643,30 @@ public:
ss << val;
TS_ASSERT_EQUALS( ss.str(), expected );
}
/**
* Test the mean calculation.
*/
void testMean( void )
{
WValue< double > val( 3 );
val[0] = 1.0;
val[1] = 2.0;
val[2] = 3.0;
TS_ASSERT_EQUALS( val.mean(), 2.0 );
}
/**
* Test the median calculation.
*/
void testMedian( void )
{
WValue< double > val( 3 );
val[0] = 1.0;
val[1] = 2.0;
val[2] = 3.0;
TS_ASSERT_EQUALS( val.mean(), 2.0 );
}
};
#endif // WVALUE_TEST_H
......@@ -173,9 +173,6 @@ void WMHomeGlyphs::properties()
void WMHomeGlyphs::moduleMain()
{
// TODO(wiebel): remove this as soon as wmath::WSymmetricSphericalHarmonic is thread save
wmath::WSymmetricSphericalHarmonic::calcLj( 4 );
m_moduleState.add( m_input->getDataChangedCondition() );
m_moduleState.add( m_recompute );
......
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