29 #include <boost/math/special_functions/spherical_harmonic.hpp>
31 #include "../exceptions/WPreconditionNotMet.h"
32 #include "linearAlgebra/WLinearAlgebra.h"
33 #include "WLinearAlgebraFunctions.h"
35 #include "WTensorSym.h"
36 #include "WUnitSphereCoordinates.h"
38 #include "WSymmetricSphericalHarmonic.h"
47 m_SHCoefficients( SHCoefficients )
52 m_order =
static_cast<size_t>( q );
63 const double rootOf2 = std::sqrt( 2.0 );
64 for(
int k = 0; k <= static_cast<int>(
m_order ); k+=2 )
67 for(
int m = 1; m <= k; m++ )
69 j = ( k*k+k+2 ) / 2 + m;
71 static_cast<double>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
74 result +=
m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
76 for(
int m = -k; m < 0; m++ )
78 j = ( k*k+k+2 ) / 2 + m;
79 result +=
m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
99 for(
int l = 0; l <= static_cast< int >(
m_order ); l += 2 )
101 for(
int m = -l; m <= l; ++m )
104 if( m < 0 && ( ( -m ) % 2 == 1 ) )
108 else if( m > 0 && ( m % 2 == 0 ) )
122 double r = 1.0 / sqrt( 2.0 );
123 std::complex< double > i( 0.0, -1.0 );
124 for(
int l = 0; l <= static_cast< int >(
m_order ); l += 2 )
126 for(
int m = -l; m <= l; ++m )
135 res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
150 double n =
static_cast< double >( orientations.size() );
156 for( std::size_t i = 0; i < orientations.size(); ++i )
158 v[ i ] =
getValue( orientations[ i ] );
163 for( std::size_t i = 0; i < orientations.size(); ++i )
172 if( gfa == 0.0 || n <= 1.0 )
177 gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
179 WAssert( gfa >= -0.01 && gfa <= 1.01,
"" );
196 for( std::size_t k = 0; k < B.
getNbCols(); ++k )
201 WAssert( s.getNbRows() == B.
getNbRows(),
"" );
202 WAssert( s.getNbCols() == 1u,
"" );
204 double n =
static_cast< double >( s.getNbRows() );
208 for( std::size_t i = 0; i < s.getNbRows(); ++i )
214 for( std::size_t i = 0; i < s.getNbRows(); ++i )
216 double f = s( i, 0 );
223 if( gfa == 0.0 || n <= 1.0 )
228 gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
230 WAssert( gfa >= -0.01 && gfa <= 1.01,
"" );
264 std::vector< WUnitSphereCoordinates > sphereCoordinates;
265 for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
286 result = pseudoInverse( result )*Bt;
291 return ( P * result );
299 WMatrix<double> B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
303 const double rootOf2 = std::sqrt( 2.0 );
304 for(uint32_t row = 0; row < orientations.size(); row++ )
306 const double theta = orientations[row].getTheta();
307 const double phi = orientations[row].getPhi();
308 for(
int k = 0; k <= order; k+=2 )
311 for(
int m = 1; m <= k; m++ )
313 j = ( k * k + k + 2 ) / 2 + m;
314 B( row, j-1 ) = rootOf2 *
static_cast<double>( std::pow( static_cast<double>( -1 ), m + 1 ) )
315 * boost::math::spherical_harmonic_i( k, m, theta, phi );
318 B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
320 for(
int m = -k; m < 0; m++ )
322 j = ( k * k + k + 2 ) / 2 + m;
323 B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
335 for( std::size_t row = 0; row < orientations.size(); row++ )
337 const double theta = orientations[ row ].getTheta();
338 const double phi = orientations[ row ].getPhi();
341 for(
int k = 0; k <= order; k += 2 )
343 for(
int m = -k; m < 0; m++ )
345 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
348 B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
350 for(
int m = 1; m <= k; m++ )
352 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
362 size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
365 for(
size_t k = 0; k <= order; k += 2 )
367 for(
int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
369 L( i, i ) =
static_cast< double > ( k * k * ( k + 1 ) * ( k + 1 ) );
378 size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
381 for(
size_t k = 0; k <= order; k += 2 )
383 double h = 2.0 * piDouble *
static_cast< double >( std::pow( static_cast< double >( -1 ), static_cast< double >( k / 2 ) ) ) *
384 static_cast< double >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast<double>( evenFactorial( k ) );
385 for(
int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
395 const std::vector< WUnitSphereCoordinates >& orientations )
397 std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
398 WPrecondEquals( order % 2, 0u );
399 WPrecondLess( numElements, orientations.size() + 1 );
402 std::vector< WVector3d > directions( numElements );
403 for( std::size_t i = 0; i < numElements; ++i )
405 directions[ i ] = orientations[ i ].getEuclidean();
409 std::vector< std::size_t > indices( order, 0 );
413 for( std::size_t j = 0; j < numElements; ++j )
416 std::size_t amount[ 3 ] = { 0, 0, 0 };
417 for( std::size_t k = 0; k < order; ++k )
419 ++amount[ indices[ k ] ];
423 std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
424 for( std::size_t i = 0; i < numElements; ++i )
426 TEMat( i, j ) = multiplicity;
427 for( std::size_t k = 0; k < order; ++k )
429 TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
434 positionIterateSortedOneStep( order, 3, indices );
439 std::vector< WUnitSphereCoordinates > ori2( orientations.begin(), orientations.begin() + numElements );