00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include <stdint.h>
00026
00027 #include <vector>
00028
00029 #include <boost/math/special_functions/spherical_harmonic.hpp>
00030
00031 #include "../exceptions/WPreconditionNotMet.h"
00032 #include "linearAlgebra/WLinearAlgebra.h"
00033 #include "WLinearAlgebraFunctions.h"
00034 #include "WMatrix.h"
00035 #include "WTensorSym.h"
00036 #include "WUnitSphereCoordinates.h"
00037
00038 #include "WSymmetricSphericalHarmonic.h"
00039
00040 WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic() :
00041 m_order( 0 ),
00042 m_SHCoefficients( 0 )
00043 {
00044 }
00045
00046 WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic( const WValue<double>& SHCoefficients ) :
00047 m_SHCoefficients( SHCoefficients )
00048 {
00049
00050
00051 double q = std::sqrt( 0.25 + 2.0 * static_cast<double>( m_SHCoefficients.size() ) ) - 1.5;
00052 m_order = static_cast<size_t>( q );
00053 }
00054
00055 WSymmetricSphericalHarmonic::~WSymmetricSphericalHarmonic()
00056 {
00057 }
00058
00059 double WSymmetricSphericalHarmonic::getValue( double theta, double phi ) const
00060 {
00061 double result = 0.0;
00062 int j = 0;
00063 const double rootOf2 = std::sqrt( 2.0 );
00064 for( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
00065 {
00066
00067 for( int m = 1; m <= k; m++ )
00068 {
00069 j = ( k*k+k+2 ) / 2 + m;
00070 result += m_SHCoefficients[ j-1 ] * rootOf2 *
00071 static_cast<double>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
00072 }
00073
00074 result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
00075
00076 for( int m = -k; m < 0; m++ )
00077 {
00078 j = ( k*k+k+2 ) / 2 + m;
00079 result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
00080 }
00081 }
00082 return result;
00083 }
00084
00085 double WSymmetricSphericalHarmonic::getValue( const WUnitSphereCoordinates& coordinates ) const
00086 {
00087 return getValue( coordinates.getTheta(), coordinates.getPhi() );
00088 }
00089
00090 const WValue<double>& WSymmetricSphericalHarmonic::getCoefficients() const
00091 {
00092 return m_SHCoefficients;
00093 }
00094
00095 WValue< double > WSymmetricSphericalHarmonic::getCoefficientsSchultz() const
00096 {
00097 WValue< double > res( m_SHCoefficients.size() );
00098 size_t k = 0;
00099 for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
00100 {
00101 for( int m = -l; m <= l; ++m )
00102 {
00103 res[ k ] = m_SHCoefficients[ k ];
00104 if( m < 0 && ( ( -m ) % 2 == 1 ) )
00105 {
00106 res[ k ] *= -1.0;
00107 }
00108 else if( m > 0 && ( m % 2 == 0 ) )
00109 {
00110 res[ k ] *= -1.0;
00111 }
00112 ++k;
00113 }
00114 }
00115 return res;
00116 }
00117
00118 WValue< std::complex< double > > WSymmetricSphericalHarmonic::getCoefficientsComplex() const
00119 {
00120 WValue< std::complex< double > > res( m_SHCoefficients.size() );
00121 size_t k = 0;
00122 double r = 1.0 / sqrt( 2.0 );
00123 std::complex< double > i( 0.0, -1.0 );
00124 for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
00125 {
00126 for( int m = -l; m <= l; ++m )
00127 {
00128 if( m == 0 )
00129 {
00130 res[ k ] = m_SHCoefficients[ k ];
00131 }
00132 else if( m < 0 )
00133 {
00134 res[ k ] += i * r * m_SHCoefficients[ k - 2 * m ];
00135 res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
00136 }
00137 else if( m > 0 )
00138 {
00139 res[ k ] += i * ( -m % 2 == 0 ? -r : r ) * m_SHCoefficients[ k ];
00140 res[ k ] += r * m_SHCoefficients[ k - 2 * m ];
00141 }
00142 ++k;
00143 }
00144 }
00145 return res;
00146 }
00147
00148 double WSymmetricSphericalHarmonic::calcGFA( std::vector< WUnitSphereCoordinates > const& orientations ) const
00149 {
00150 double n = static_cast< double >( orientations.size() );
00151 double d = 0.0;
00152 double gfa = 0.0;
00153 double mean = 0.0;
00154 double v[ 15 ];
00155
00156 for( std::size_t i = 0; i < orientations.size(); ++i )
00157 {
00158 v[ i ] = getValue( orientations[ i ] );
00159 mean += v[ i ];
00160 }
00161 mean /= n;
00162
00163 for( std::size_t i = 0; i < orientations.size(); ++i )
00164 {
00165 double f = v[ i ];
00166
00167 gfa += f * f;
00168 f -= mean;
00169 d += f * f;
00170 }
00171
00172 if( gfa == 0.0 || n <= 1.0 )
00173 {
00174 return 0.0;
00175 }
00176
00177 gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
00178
00179 WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
00180 if( gfa < 0.0 )
00181 {
00182 return 0.0;
00183 }
00184 else if( gfa > 1.0 )
00185 {
00186 return 1.0;
00187 }
00188 return gfa;
00189 }
00190
00191 double WSymmetricSphericalHarmonic::calcGFA( WMatrix< double > const& B ) const
00192 {
00193
00194 WMatrix< double > s( B.getNbCols(), 1 );
00195 WAssert( B.getNbCols() == m_SHCoefficients.size(), "" );
00196 for( std::size_t k = 0; k < B.getNbCols(); ++k )
00197 {
00198 s( k, 0 ) = m_SHCoefficients[ k ];
00199 }
00200 s = B * s;
00201 WAssert( s.getNbRows() == B.getNbRows(), "" );
00202 WAssert( s.getNbCols() == 1u, "" );
00203
00204 double n = static_cast< double >( s.getNbRows() );
00205 double d = 0.0;
00206 double gfa = 0.0;
00207 double mean = 0.0;
00208 for( std::size_t i = 0; i < s.getNbRows(); ++i )
00209 {
00210 mean += s( i, 0 );
00211 }
00212 mean /= n;
00213
00214 for( std::size_t i = 0; i < s.getNbRows(); ++i )
00215 {
00216 double f = s( i, 0 );
00217
00218 gfa += f * f;
00219 f -= mean;
00220 d += f * f;
00221 }
00222
00223 if( gfa == 0.0 || n <= 1.0 )
00224 {
00225 return 0.0;
00226 }
00227
00228 gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
00229
00230 WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
00231 if( gfa < 0.0 )
00232 {
00233 return 0.0;
00234 }
00235 else if( gfa > 1.0 )
00236 {
00237 return 1.0;
00238 }
00239 return gfa;
00240 }
00241
00242 void WSymmetricSphericalHarmonic::applyFunkRadonTransformation( WMatrix< double > const& frtMat )
00243 {
00244 WAssert( frtMat.getNbCols() == m_SHCoefficients.size(), "" );
00245 WAssert( frtMat.getNbRows() == m_SHCoefficients.size(), "" );
00246
00247 for( size_t j = 0; j < m_SHCoefficients.size(); j++ )
00248 {
00249 m_SHCoefficients[ j ] *= frtMat( j, j );
00250 }
00251 }
00252
00253 size_t WSymmetricSphericalHarmonic::getOrder() const
00254 {
00255 return m_order;
00256 }
00257
00258 WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WVector3d >& orientations,
00259 int order,
00260 double lambda,
00261 bool withFRT )
00262 {
00263
00264 std::vector< WUnitSphereCoordinates > sphereCoordinates;
00265 for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
00266 {
00267 sphereCoordinates.push_back( WUnitSphereCoordinates( *it ) );
00268 }
00269 return WSymmetricSphericalHarmonic::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
00270 }
00271
00272 WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
00273 int order,
00274 double lambda,
00275 bool withFRT )
00276 {
00277 WMatrix<double> B( WSymmetricSphericalHarmonic::calcBaseMatrix( orientations, order ) );
00278 WMatrix<double> Bt( B.transposed() );
00279 WMatrix<double> result( Bt*B );
00280 if( lambda != 0.0 )
00281 {
00282 WMatrix<double> L( WSymmetricSphericalHarmonic::calcSmoothingMatrix( order ) );
00283 result += lambda*L;
00284 }
00285
00286 result = pseudoInverse( result )*Bt;
00287
00288 if( withFRT )
00289 {
00290 WMatrix<double> P( WSymmetricSphericalHarmonic::calcFRTMatrix( order ) );
00291 return ( P * result );
00292 }
00293 return result;
00294 }
00295
00296 WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
00297 int order )
00298 {
00299 WMatrix<double> B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
00300
00301
00302 int j = 0;
00303 const double rootOf2 = std::sqrt( 2.0 );
00304 for(uint32_t row = 0; row < orientations.size(); row++ )
00305 {
00306 const double theta = orientations[row].getTheta();
00307 const double phi = orientations[row].getPhi();
00308 for( int k = 0; k <= order; k+=2 )
00309 {
00310
00311 for( int m = 1; m <= k; m++ )
00312 {
00313 j = ( k * k + k + 2 ) / 2 + m;
00314 B( row, j-1 ) = rootOf2 * static_cast<double>( std::pow( static_cast<double>( -1 ), m + 1 ) )
00315 * boost::math::spherical_harmonic_i( k, m, theta, phi );
00316 }
00317
00318 B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
00319
00320 for( int m = -k; m < 0; m++ )
00321 {
00322 j = ( k * k + k + 2 ) / 2 + m;
00323 B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
00324 }
00325 }
00326 }
00327 return B;
00328 }
00329
00330 WMatrix< std::complex< double > >
00331 WSymmetricSphericalHarmonic::calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations, int order )
00332 {
00333 WMatrix< std::complex< double > > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
00334
00335 for( std::size_t row = 0; row < orientations.size(); row++ )
00336 {
00337 const double theta = orientations[ row ].getTheta();
00338 const double phi = orientations[ row ].getPhi();
00339
00340 int j = 0;
00341 for( int k = 0; k <= order; k += 2 )
00342 {
00343 for( int m = -k; m < 0; m++ )
00344 {
00345 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
00346 ++j;
00347 }
00348 B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
00349 ++j;
00350 for( int m = 1; m <= k; m++ )
00351 {
00352 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
00353 ++j;
00354 }
00355 }
00356 }
00357 return B;
00358 }
00359
00360 WMatrix<double> WSymmetricSphericalHarmonic::calcSmoothingMatrix( size_t order )
00361 {
00362 size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
00363 std::size_t i = 0;
00364 WMatrix<double> L( R, R );
00365 for( size_t k = 0; k <= order; k += 2 )
00366 {
00367 for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
00368 {
00369 L( i, i ) = static_cast< double > ( k * k * ( k + 1 ) * ( k + 1 ) );
00370 ++i;
00371 }
00372 }
00373 return L;
00374 }
00375
00376 WMatrix<double> WSymmetricSphericalHarmonic::calcFRTMatrix( size_t order )
00377 {
00378 size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
00379 std::size_t i = 0;
00380 WMatrix< double > result( R, R );
00381 for( size_t k = 0; k <= order; k += 2 )
00382 {
00383 double h = 2.0 * piDouble * static_cast< double >( std::pow( static_cast< double >( -1 ), static_cast< double >( k / 2 ) ) ) *
00384 static_cast< double >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast<double>( evenFactorial( k ) );
00385 for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
00386 {
00387 result( i, i ) = h;
00388 ++i;
00389 }
00390 }
00391 return result;
00392 }
00393
00394 WMatrix< double > WSymmetricSphericalHarmonic::calcSHToTensorSymMatrix( std::size_t order,
00395 const std::vector< WUnitSphereCoordinates >& orientations )
00396 {
00397 std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
00398 WPrecondEquals( order % 2, 0u );
00399 WPrecondLess( numElements, orientations.size() + 1 );
00400
00401
00402 std::vector< WVector3d > directions( numElements );
00403 for( std::size_t i = 0; i < numElements; ++i )
00404 {
00405 directions[ i ] = orientations[ i ].getEuclidean();
00406 }
00407
00408
00409 std::vector< std::size_t > indices( order, 0 );
00410
00411
00412 WMatrix< double > TEMat( numElements, numElements );
00413 for( std::size_t j = 0; j < numElements; ++j )
00414 {
00415
00416 std::size_t amount[ 3 ] = { 0, 0, 0 };
00417 for( std::size_t k = 0; k < order; ++k )
00418 {
00419 ++amount[ indices[ k ] ];
00420 }
00421
00422
00423 std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
00424 for( std::size_t i = 0; i < numElements; ++i )
00425 {
00426 TEMat( i, j ) = multiplicity;
00427 for( std::size_t k = 0; k < order; ++k )
00428 {
00429 TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
00430 }
00431 }
00432
00433
00434 positionIterateSortedOneStep( order, 3, indices );
00435 }
00436 directions.clear();
00437
00438
00439 std::vector< WUnitSphereCoordinates > ori2( orientations.begin(), orientations.begin() + numElements );
00440
00441 WMatrix< double > p = pseudoInverse( TEMat );
00442
00443 return p * calcBaseMatrix( ori2, order );
00444 }
00445
00446 void WSymmetricSphericalHarmonic::normalize()
00447 {
00448 double scale = 0.0;
00449 if( m_SHCoefficients.size() > 0 )
00450 scale = std::sqrt( 4.0 * piDouble ) * m_SHCoefficients[0];
00451 for( size_t i = 0; i < m_SHCoefficients.size(); i++ )
00452 {
00453 m_SHCoefficients[ i ] /= scale;
00454 }
00455 }
00456