25 #ifndef WSYMMETRICSPHERICALHARMONIC_H
26 #define WSYMMETRICSPHERICALHARMONIC_H
33 #include <boost/math/special_functions/spherical_harmonic.hpp>
35 #include "core/common/WLogger.h"
36 #include "core/common/math/WGeometryFunctions.h"
37 #include "../exceptions/WPreconditionNotMet.h"
38 #include "linearAlgebra/WVectorFixed.h"
39 #include "WLinearAlgebraFunctions.h"
42 #include "WTensorSym.h"
43 #include "WUnitSphereCoordinates.h"
278 template<
typename T >
281 m_SHCoefficients( 0 )
285 template<
typename T >
287 m_SHCoefficients( SHCoefficients )
291 T q = std::sqrt( 0.25 + 2.0 * static_cast< T >(
m_SHCoefficients.size() ) ) - 1.5;
292 m_order =
static_cast<size_t>( q );
295 template<
typename T >
300 template<
typename T >
305 const T rootOf2 = std::sqrt( 2.0 );
306 for(
int k = 0; k <= static_cast<int>( m_order ); k+=2 )
309 for(
int m = 1; m <= k; m++ )
311 j = ( k*k+k+2 ) / 2 + m;
312 result += m_SHCoefficients[ j-1 ] * rootOf2 *
313 static_cast< T
>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
316 result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
318 for(
int m = -k; m < 0; m++ )
320 j = ( k*k+k+2 ) / 2 + m;
321 result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
327 template<
typename T >
333 template<
typename T >
336 return m_SHCoefficients;
339 template<
typename T >
344 for(
int l = 0; l <= static_cast< int >( m_order ); l += 2 )
346 for(
int m = -l; m <= l; ++m )
348 res[ k ] = m_SHCoefficients[ k ];
349 if( m < 0 && ( ( -m ) % 2 == 1 ) )
353 else if( m > 0 && ( m % 2 == 0 ) )
363 template<
typename T >
368 T r = 1.0 / sqrt( 2.0 );
369 std::complex< T > i( 0.0, -1.0 );
370 for(
int l = 0; l <= static_cast< int >( m_order ); l += 2 )
372 for(
int m = -l; m <= l; ++m )
376 res[ k ] = m_SHCoefficients[ k ];
380 res[ k ] += i * r * m_SHCoefficients[ k - 2 * m ];
381 res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
385 res[ k ] += i * ( -m % 2 == 0 ? -r : r ) * m_SHCoefficients[ k ];
386 res[ k ] += r * m_SHCoefficients[ k - 2 * m ];
394 template<
typename T >
397 T n =
static_cast< T
>( orientations.size() );
401 std::vector< T > v( orientations.size() );
403 for( std::size_t i = 0; i < orientations.size(); ++i )
405 v[ i ] = getValue( orientations[ i ] );
410 for( std::size_t i = 0; i < orientations.size(); ++i )
419 if( gfa == 0.0 || n <= 1.0 )
424 gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
426 WAssert( gfa >= -0.01 && gfa <= 1.01,
"" );
438 template<
typename T >
443 WAssert( B.
getNbCols() == m_SHCoefficients.size(),
"" );
444 for( std::size_t k = 0; k < B.
getNbCols(); ++k )
446 s( k, 0 ) = m_SHCoefficients[ k ];
449 WAssert( s.getNbRows() == B.
getNbRows(),
"" );
450 WAssert( s.getNbCols() == 1u,
"" );
452 T n =
static_cast< T
>( s.getNbRows() );
456 for( std::size_t i = 0; i < s.getNbRows(); ++i )
462 for( std::size_t i = 0; i < s.getNbRows(); ++i )
471 if( gfa == 0.0 || n <= 1.0 )
476 gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
478 WAssert( gfa >= -0.01 && gfa <= 1.01,
"" );
490 template<
typename T >
493 WAssert( frtMat.
getNbCols() == m_SHCoefficients.size(),
"" );
494 WAssert( frtMat.
getNbRows() == m_SHCoefficients.size(),
"" );
496 for(
size_t j = 0; j < m_SHCoefficients.size(); j++ )
498 m_SHCoefficients[ j ] *= frtMat( j, j );
502 template<
typename T >
508 template<
typename T >
515 std::vector< WUnitSphereCoordinates< T > > sphereCoordinates;
516 for(
typename std::vector<
WMatrixFixed< T, 3, 1 > >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
523 template<
typename T >
538 result = pseudoInverse( result )*Bt;
542 return ( P * result );
547 template<
typename T >
553 std::vector< WUnitSphereCoordinates< T > > sphereCoordinates;
554 for(
typename std::vector<
WMatrixFixed< T, 3, 1 > >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
561 template<
typename T >
578 result = pseudoInverse( result )*Bt;
593 T correctionFactor = 1.0 / ( 16.0 * std::pow( static_cast< T >( piDouble ), 2 ) );
594 result *= correctionFactor;
600 template<
typename T >
604 WMatrix< T > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
608 const T rootOf2 = std::sqrt( 2.0 );
609 for(uint32_t row = 0; row < orientations.size(); row++ )
611 const T theta = orientations[row].getTheta();
612 const T phi = orientations[row].getPhi();
613 for(
int k = 0; k <= order; k+=2 )
616 for(
int m = 1; m <= k; m++ )
618 j = ( k * k + k + 2 ) / 2 + m;
619 B( row, j-1 ) = rootOf2 *
static_cast< T
>( std::pow( static_cast< T >( -1 ), m + 1 ) )
620 * boost::math::spherical_harmonic_i( k, m, theta, phi );
623 B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
625 for(
int m = -k; m < 0; m++ )
627 j = ( k * k + k + 2 ) / 2 + m;
628 B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
635 template<
typename T >
641 for( std::size_t row = 0; row < orientations.size(); row++ )
643 const T theta = orientations[ row ].getTheta();
644 const T phi = orientations[ row ].getPhi();
647 for(
int k = 0; k <= order; k += 2 )
649 for(
int m = -k; m < 0; m++ )
651 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
654 B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
656 for(
int m = 1; m <= k; m++ )
658 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
666 template<
typename T >
669 size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
672 for(
size_t k = 0; k <= order; k += 2 )
674 for(
int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
676 result[ i ] = -
static_cast< T
> ( k * ( k + 1 ) );
683 template<
typename T >
688 for( std::size_t i = 0; i < eigenvalues.
size(); ++i )
690 L( i, i ) = eigenvalues[ i ];
695 template<
typename T >
700 for( std::size_t i = 0; i < eigenvalues.
size(); ++i )
702 L( i, i ) = std::pow( eigenvalues[ i ], 2 );
707 template<
typename T >
710 size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
713 for(
size_t k = 0; k <= order; k += 2 )
715 T h = 2.0 *
static_cast< T
>( piDouble ) * static_cast< T >( std::pow( static_cast< T >( -1 ), static_cast< T >( k / 2 ) ) ) *
716 static_cast< T >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast< T >( evenFactorial( k ) );
717 for(
int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
726 template<
typename T >
729 std::vector< WVector3d > vertices;
730 std::vector< unsigned int > triangles;
732 int level =
static_cast< int >( log( static_cast< float >( ( order + 1 ) * ( order + 2 ) ) ) / 2.0f + 1.0f );
733 tesselateIcosahedron( &vertices, &triangles, level );
734 std::vector< WUnitSphereCoordinates< T > > orientations;
735 for(
typename std::vector< WVector3d >::const_iterator cit = vertices.begin(); cit != vertices.end(); cit++ )
738 if( ( *cit )[ 0 ] >= 0.0001 )
746 template<
typename T >
750 std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
751 WPrecondEquals( order % 2, 0u );
752 WPrecondLess( numElements, orientations.size() + 1 );
755 std::vector< WMatrixFixed< T, 3, 1 > > directions( numElements );
756 for( std::size_t i = 0; i < numElements; ++i )
758 directions[ i ] = orientations[ i ].getEuclidean();
762 std::vector< std::size_t > indices( order, 0 );
766 for( std::size_t j = 0; j < numElements; ++j )
769 std::size_t amount[ 3 ] = { 0, 0, 0 };
770 for( std::size_t k = 0; k < order; ++k )
772 ++amount[ indices[ k ] ];
776 std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
777 for( std::size_t i = 0; i < numElements; ++i )
779 TEMat( i, j ) = multiplicity;
780 for( std::size_t k = 0; k < order; ++k )
782 TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
787 positionIterateSortedOneStep( order, 3, indices );
792 std::vector< WUnitSphereCoordinates< T > > ori2( orientations.begin(), orientations.begin() + numElements );
796 WAssert( ( TEMat*p ).isIdentity( 1e-3 ),
"Test of inverse matrix failed. Are the given orientations linear independent?" );
798 return p * calcBaseMatrix( ori2, order );
801 template<
typename T >
805 if( m_SHCoefficients.size() > 0 )
807 scale = std::sqrt( 4.0 * static_cast< T >( piDouble ) ) * m_SHCoefficients[0];
809 for(
size_t i = 0; i < m_SHCoefficients.size(); i++ )
811 m_SHCoefficients[ i ] /= scale;
815 #endif // WSYMMETRICSPHERICALHARMONIC_H
WSymmetricSphericalHarmonic()
Default constructor.
static WMatrix< T > calcFRTMatrix(size_t order)
Calculates the Funk-Radon-Transformation-Matrix P from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".
WValue< T > m_SHCoefficients
coefficients of the spherical harmonic
static WMatrix< T > getSHFittingMatrixForConstantSolidAngle(const std::vector< WMatrixFixed< T, 3, 1 > > &orientations, int order, T lambda)
This calculates the transformation/fitting matrix T like in the 2010 Aganj paper. ...
static WMatrix< T > calcSmoothingMatrix(size_t order)
This calcs the smoothing matrix L from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".
T getValue(T theta, T phi) const
Return the value on the sphere.
void normalize()
Normalize this SH in place.
Base class for all higher level values like tensors, vectors, matrices and so on. ...
size_t getNbCols() const
Get number of columns.
virtual ~WSymmetricSphericalHarmonic()
Destructor.
This class stores coordinates on the unit sphere.
T getPhi() const
Return the phi angle.
Matrix template class with variable number of rows and columns.
const WValue< T > & getCoefficients() const
Returns the used coefficients (stored like in the mentioned 2007 Descoteaux paper).
WValue< std::complex< T > > getCoefficientsComplex() const
Returns the coefficients for the complex base.
size_t getOrder() const
Return the order of the spherical harmonic.
WMatrix transposed() const
Returns the transposed matrix.
static WValue< T > calcEigenvalues(size_t order)
Calc eigenvalues for SH elements.
Class for symmetric spherical harmonics The index scheme of the coefficients/basis values is like in ...
T calcGFA(std::vector< WUnitSphereCoordinates< T > > const &orientations) const
Calculate the generalized fractional anisotropy for this ODF.
static WMatrix< T > calcBaseMatrix(const std::vector< WUnitSphereCoordinates< T > > &orientations, int order)
Calculates the base matrix B like in the dissertation of Descoteaux.
static WMatrix< T > calcSHToTensorSymMatrix(std::size_t order)
Calculates a matrix that converts spherical harmonics to symmetric tensors of equal or lower order...
size_t m_order
order of the spherical harmonic
size_t size() const
Get number of components the value consists of.
static WMatrix< T > calcMatrixWithEigenvalues(size_t order)
Calc matrix with the eigenvalues of the SH elements on its diagonal.
size_t getNbRows() const
Get number of rows.
T getTheta() const
Return the theta angle.
static WMatrix< std::complex< T > > calcComplexBaseMatrix(std::vector< WUnitSphereCoordinates< T > > const &orientations, int order)
Calculates the base matrix B for the complex spherical harmonics.
WValue< T > getCoefficientsSchultz() const
Returns the coefficients for Schultz' SH base.
WStreamedLogger debug(const std::string &source)
Logging a debug message.
void applyFunkRadonTransformation(WMatrix< T > const &frtMat)
Applies the Funk-Radon-Transformation.
static WMatrix< T > getSHFittingMatrix(const std::vector< WMatrixFixed< T, 3, 1 > > &orientations, int order, T lambda, bool withFRT)
This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper.