29 #include <boost/math/special_functions/spherical_harmonic.hpp>
31 #include "core/common/WLogger.h"
32 #include "../exceptions/WPreconditionNotMet.h"
33 #include "linearAlgebra/WLinearAlgebra.h"
34 #include "WLinearAlgebraFunctions.h"
36 #include "WTensorSym.h"
37 #include "WUnitSphereCoordinates.h"
39 #include "WSymmetricSphericalHarmonic.h"
48 m_SHCoefficients( SHCoefficients )
53 m_order =
static_cast<size_t>( q );
64 const double rootOf2 = std::sqrt( 2.0 );
65 for(
int k = 0; k <= static_cast<int>(
m_order ); k+=2 )
68 for(
int m = 1; m <= k; m++ )
70 j = ( k*k+k+2 ) / 2 + m;
72 static_cast<double>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
75 result +=
m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
77 for(
int m = -k; m < 0; m++ )
79 j = ( k*k+k+2 ) / 2 + m;
80 result +=
m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
100 for(
int l = 0; l <= static_cast< int >(
m_order ); l += 2 )
102 for(
int m = -l; m <= l; ++m )
105 if( m < 0 && ( ( -m ) % 2 == 1 ) )
109 else if( m > 0 && ( m % 2 == 0 ) )
123 double r = 1.0 / sqrt( 2.0 );
124 std::complex< double > i( 0.0, -1.0 );
125 for(
int l = 0; l <= static_cast< int >(
m_order ); l += 2 )
127 for(
int m = -l; m <= l; ++m )
136 res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
151 double n =
static_cast< double >( orientations.size() );
157 for( std::size_t i = 0; i < orientations.size(); ++i )
159 v[ i ] =
getValue( orientations[ i ] );
164 for( std::size_t i = 0; i < orientations.size(); ++i )
173 if( gfa == 0.0 || n <= 1.0 )
178 gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
180 WAssert( gfa >= -0.01 && gfa <= 1.01,
"" );
197 for( std::size_t k = 0; k < B.
getNbCols(); ++k )
202 WAssert( s.getNbRows() == B.
getNbRows(),
"" );
203 WAssert( s.getNbCols() == 1u,
"" );
205 double n =
static_cast< double >( s.getNbRows() );
209 for( std::size_t i = 0; i < s.getNbRows(); ++i )
215 for( std::size_t i = 0; i < s.getNbRows(); ++i )
217 double f = s( i, 0 );
224 if( gfa == 0.0 || n <= 1.0 )
229 gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
231 WAssert( gfa >= -0.01 && gfa <= 1.01,
"" );
265 std::vector< WUnitSphereCoordinates > sphereCoordinates;
266 for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
287 result = pseudoInverse( result )*Bt;
291 return ( P * result );
301 std::vector< WUnitSphereCoordinates > sphereCoordinates;
302 for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
310 const std::vector< WUnitSphereCoordinates >& orientations,
325 result = pseudoInverse( result )*Bt;
340 double correctionFactor = 1.0 / ( 16.0 * std::pow( piDouble, 2 ) );
341 result *= correctionFactor;
350 WMatrix<double> B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
354 const double rootOf2 = std::sqrt( 2.0 );
355 for(uint32_t row = 0; row < orientations.size(); row++ )
357 const double theta = orientations[row].getTheta();
358 const double phi = orientations[row].getPhi();
359 for(
int k = 0; k <= order; k+=2 )
362 for(
int m = 1; m <= k; m++ )
364 j = ( k * k + k + 2 ) / 2 + m;
365 B( row, j-1 ) = rootOf2 *
static_cast<double>( std::pow( static_cast<double>( -1 ), m + 1 ) )
366 * boost::math::spherical_harmonic_i( k, m, theta, phi );
369 B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
371 for(
int m = -k; m < 0; m++ )
373 j = ( k * k + k + 2 ) / 2 + m;
374 B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
386 for( std::size_t row = 0; row < orientations.size(); row++ )
388 const double theta = orientations[ row ].getTheta();
389 const double phi = orientations[ row ].getPhi();
392 for(
int k = 0; k <= order; k += 2 )
394 for(
int m = -k; m < 0; m++ )
396 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
399 B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
401 for(
int m = 1; m <= k; m++ )
403 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
413 size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
416 for(
size_t k = 0; k <= order; k += 2 )
418 for(
int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
420 result[ i ] = -
static_cast< double > ( k * ( k + 1 ) );
431 for( std::size_t i = 0; i < eigenvalues.
size(); ++i )
433 L( i, i ) = eigenvalues[ i ];
442 for( std::size_t i = 0; i < eigenvalues.
size(); ++i )
444 L( i, i ) = std::pow( eigenvalues[ i ], 2 );
451 size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
454 for(
size_t k = 0; k <= order; k += 2 )
456 double h = 2.0 * piDouble *
static_cast< double >( std::pow( static_cast< double >( -1 ), static_cast< double >( k / 2 ) ) ) *
457 static_cast< double >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast<double>( evenFactorial( k ) );
458 for(
int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
468 const std::vector< WUnitSphereCoordinates >& orientations )
470 std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
471 WPrecondEquals( order % 2, 0u );
472 WPrecondLess( numElements, orientations.size() + 1 );
475 std::vector< WVector3d > directions( numElements );
476 for( std::size_t i = 0; i < numElements; ++i )
478 directions[ i ] = orientations[ i ].getEuclidean();
482 std::vector< std::size_t > indices( order, 0 );
486 for( std::size_t j = 0; j < numElements; ++j )
489 std::size_t amount[ 3 ] = { 0, 0, 0 };
490 for( std::size_t k = 0; k < order; ++k )
492 ++amount[ indices[ k ] ];
496 std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
497 for( std::size_t i = 0; i < numElements; ++i )
499 TEMat( i, j ) = multiplicity;
500 for( std::size_t k = 0; k < order; ++k )
502 TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
507 positionIterateSortedOneStep( order, 3, indices );
512 std::vector< WUnitSphereCoordinates > ori2( orientations.begin(), orientations.begin() + numElements );