diff --git a/src/common/math/WSymmetricSphericalHarmonic.cpp b/src/common/math/WSymmetricSphericalHarmonic.cpp index 9239553a466b9712df35ab9d020f360a0620ba95..59fb0312fd10f6284bd8c22ff2da11157cc905ad 100644 --- a/src/common/math/WSymmetricSphericalHarmonic.cpp +++ b/src/common/math/WSymmetricSphericalHarmonic.cpp @@ -22,16 +22,14 @@ // //--------------------------------------------------------------------------- -#include - #include #include "stdint.h" #include "WLinearAlgebraFunctions.h" #include "WUnitSphereCoordinates.h" -#include "WValue.h" #include "WMatrix.h" +#include "WValue.h" #include "WSymmetricSphericalHarmonic.h" @@ -119,6 +117,7 @@ const wmath::WValue& WSymmetricSphericalHarmonic::getCoefficients() cons void WSymmetricSphericalHarmonic::applyFunkRadonTransformation() { // Funk-Radon-Transformation like in the 2007 Descoteaux paper + // TODO(philips): implement size_t lj = 0; // double ljDbl = 0.0; double factor = 0.0; @@ -148,9 +147,9 @@ size_t WSymmetricSphericalHarmonic::getOrder() const void WSymmetricSphericalHarmonic::calcLj( size_t order ) { m_lj.clear(); - for ( size_t k = 0; k <= order; k += 2 ) + for ( size_t k=0; k<=order; k+=2 ) { - for ( int m = -k; m <= static_cast( k ); m++ ) + for ( int m=-k; m<=static_cast( k ); m++ ) { m_lj.push_back( k ); } @@ -158,10 +157,7 @@ void WSymmetricSphericalHarmonic::calcLj( size_t order ) } -wmath::WMatrix WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< wmath::WVector3D >& orientations, - int order, - double lambda, - bool withFRT ) +wmath::WMatrix WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< wmath::WVector3D >& orientations, int order, double lambda, bool withFRT ) { // convert euclid 3d orientations/gradients to sphere coordinates std::vector< wmath::WUnitSphereCoordinates > sphereCoordinates; @@ -172,18 +168,15 @@ wmath::WMatrix WSymmetricSphericalHarmonic::getSHFittingMatrix( const st return WSymmetricSphericalHarmonic::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT ); } -wmath::WMatrix WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< wmath::WUnitSphereCoordinates >& orientations, - int order, - double lambda, - bool withFRT ) +wmath::WMatrix WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< wmath::WUnitSphereCoordinates >& orientations, int order, double lambda, bool withFRT ) { - wmath::WMatrix B( WSymmetricSphericalHarmonic::calcBMatrix( orientations, order ) ); + wmath::WMatrix B( WSymmetricSphericalHarmonic::calcBaseMatrix( orientations, order ) ); wmath::WMatrix Bt( B.transposed() ); wmath::WMatrix result( Bt*B ); if ( lambda != 0.0 ) { wmath::WMatrix L( WSymmetricSphericalHarmonic::calcSmoothingMatrix( order ) ); - result += lambda*L; + result += lambda*L; } result = wmath::pseudoInverse( result )*Bt; // result *= Bt; @@ -195,9 +188,9 @@ wmath::WMatrix WSymmetricSphericalHarmonic::getSHFittingMatrix( const st return result; } -wmath::WMatrix WSymmetricSphericalHarmonic::calcBMatrix( const std::vector< wmath::WUnitSphereCoordinates >& orientations, int order ) +wmath::WMatrix WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector< wmath::WUnitSphereCoordinates >& orientations, int order ) { - wmath::WMatrix B( orientations.size(), ( ( order+1 ) * ( order+2 ) ) / 2 ); + wmath::WMatrix 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; @@ -212,15 +205,14 @@ wmath::WMatrix WSymmetricSphericalHarmonic::calcBMatrix( const std::vect for ( int m = 1; m <= k; m++ ) { j = ( k*k+k+2 ) / 2 + m; - B( row, j-1 ) = rootOf2 * static_cast( std::pow( static_cast( -1 ), m+1 ) ) - * boost::math::spherical_harmonic_i( k, m, theta, phi ); + B( row, j-1 ) = rootOf2 * static_cast( std::pow( static_cast( -1 ), m+1 ) ) * 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 ); + 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; + j = (k*k+k+2)/2+m; B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi ); } } @@ -228,7 +220,7 @@ wmath::WMatrix WSymmetricSphericalHarmonic::calcBMatrix( const std::vect return B; // old implementation // wmath::WMatrix 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++ ) @@ -261,14 +253,14 @@ wmath::WMatrix WSymmetricSphericalHarmonic::calcSmoothingMatrix( size_t { std::vector lj; lj.clear(); - for ( size_t k = 0; k <= order; k += 2 ) + for ( size_t k=0; k<=order; k+=2 ) { - for ( int m =- k; m <= static_cast( k ); m++ ) + for ( int m=-k; m<=static_cast( k ); m++ ) { lj.push_back( k ); } } - size_t R = ( ( order+1 ) * ( order+2 ) ) / 2; + size_t R = ( (order+1) * (order+2) ) / 2; wmath::WMatrix L( R, R ); for ( size_t i = 0; i < R; i++ ) L( i, i ) = static_cast ( lj[i] * lj[i] * ( lj[i] + 1) * ( lj[i] + 1) ); return L; @@ -278,8 +270,9 @@ wmath::WMatrix WSymmetricSphericalHarmonic::calcFRTMatrix( size_t order { size_t dim = ( order + 1 ) * ( order + 2 ) / 2; if ( WSymmetricSphericalHarmonic::m_lj.size() < dim ) WSymmetricSphericalHarmonic::calcLj( order ); - + wmath::WMatrix result( dim, dim ); + double factor = 0.0; for ( size_t j = 0; j < dim; j++ ) { size_t lj = WSymmetricSphericalHarmonic::m_lj[j];