OpenWalnut  1.4.0
WSymmetricSphericalHarmonic.h
00001 //---------------------------------------------------------------------------
00002 //
00003 // Project: OpenWalnut ( http://www.openwalnut.org )
00004 //
00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
00006 // For more information see http://www.openwalnut.org/copying
00007 //
00008 // This file is part of OpenWalnut.
00009 //
00010 // OpenWalnut is free software: you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as published by
00012 // the Free Software Foundation, either version 3 of the License, or
00013 // (at your option) any later version.
00014 //
00015 // OpenWalnut is distributed in the hope that it will be useful,
00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 // GNU Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public License
00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
00022 //
00023 //---------------------------------------------------------------------------
00024 
00025 #ifndef WSYMMETRICSPHERICALHARMONIC_H
00026 #define WSYMMETRICSPHERICALHARMONIC_H
00027 
00028 #include <stdint.h>
00029 
00030 #include <cmath>
00031 #include <vector>
00032 
00033 #include <boost/math/special_functions/spherical_harmonic.hpp>
00034 
00035 #include "core/common/WLogger.h"
00036 #include "core/common/math/WGeometryFunctions.h"
00037 #include "../exceptions/WPreconditionNotMet.h"
00038 #include "linearAlgebra/WVectorFixed.h"
00039 #include "WLinearAlgebraFunctions.h"
00040 #include "WMath.h"
00041 #include "WMatrix.h"
00042 #include "WTensorSym.h"
00043 #include "WUnitSphereCoordinates.h"
00044 #include "WValue.h"
00045 
00046 
00047 /**
00048  * Class for symmetric spherical harmonics
00049  * The index scheme of the coefficients/basis values is like in the Descoteaux paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".
00050  */
00051 template< typename T > class WSymmetricSphericalHarmonic  // NOLINT
00052 {
00053 // friend class WSymmetricSphericalHarmonicTest;
00054 public:
00055     /**
00056      * Default constructor.
00057      */
00058     WSymmetricSphericalHarmonic();
00059 
00060     /**
00061      * Constructor.
00062      * \param SHCoefficients the initial coefficients (stored like in the mentioned Descoteaux paper).
00063      */
00064     explicit WSymmetricSphericalHarmonic( const WValue< T >& SHCoefficients );
00065 
00066     /**
00067      * Destructor.
00068      */
00069     virtual ~WSymmetricSphericalHarmonic();
00070 
00071     /**
00072      * Return the value on the sphere.
00073      * \param theta angle for the position on the unit sphere
00074      * \param phi angle for the position on the unit sphere
00075      *
00076      * \return value on sphere
00077      */
00078     T getValue( T theta, T phi ) const;
00079 
00080     /**
00081      * Return the value on the sphere.
00082      * \param coordinates for the position on the unit sphere
00083      *
00084      * \return value on sphere
00085      */
00086     T getValue( const WUnitSphereCoordinates< T >& coordinates ) const;
00087 
00088     /**
00089      * Returns the used coefficients (stored like in the mentioned 2007 Descoteaux paper).
00090      *
00091      * \return coefficient list
00092      */
00093     const WValue< T >& getCoefficients() const;
00094 
00095     /**
00096      * Returns the coefficients for Schultz' SH base.
00097      *
00098      * \return coefficient list
00099      */
00100     WValue< T > getCoefficientsSchultz() const;
00101 
00102     /**
00103      * Returns the coefficients for the complex base.
00104      *
00105      * \return coefficiend list
00106      */
00107     WValue< std::complex< T > > getCoefficientsComplex() const;
00108 
00109     /**
00110      * Applies the Funk-Radon-Transformation. This is faster than matrix multiplication.
00111      * ( O(n) instead of O(n²) )
00112      *
00113      * \param frtMat the frt matrix as calculated by calcFRTMatrix()
00114      */
00115     void applyFunkRadonTransformation( WMatrix< T > const& frtMat );
00116 
00117     /**
00118      * Return the order of the spherical harmonic.
00119      *
00120      * \return order of SH
00121      */
00122     size_t getOrder() const;
00123 
00124     /**
00125      * Calculate the generalized fractional anisotropy for this ODF.
00126      *
00127      * See: David S. Tuch, "Q-Ball Imaging", Magn. Reson. Med. 52, 2004, 1358-1372
00128      *
00129      * \note this only makes sense if this is an ODF, meaning funk-radon-transform was applied etc.
00130      *
00131      * \param orientations A vector of unit sphere coordinates.
00132      *
00133      * \return The generalized fractional anisotropy.
00134      */
00135     T calcGFA( std::vector< WUnitSphereCoordinates< T > > const& orientations ) const;
00136 
00137     /**
00138      * Calculate the generalized fractional anisotropy for this ODF. This version of
00139      * the function uses precomputed base functions (because calculating the base function values
00140      * is rather expensive). Use this version if you want to compute the GFA for multiple ODFs
00141      * with the same base functions. The base function Matrix can be computed using \see calcBMatrix().
00142      *
00143      * See: David S. Tuch, "Q-Ball Imaging", Magn. Reson. Med. 52, 2004, 1358-1372
00144      *
00145      * \note this only makes sense if this is an ODF, meaning funk-radon-transform was applied etc.
00146      *
00147      * \param B The matrix of SH base functions.
00148      *
00149      * \return The generalized fractional anisotropy.
00150      */
00151     T calcGFA( WMatrix< T > const& B ) const;
00152 
00153     /**
00154     * This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper. The orientations are given as WMatrixFixed< T, 3, 1 >.
00155     * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
00156     * \param order The order of the spherical harmonics intended to create
00157     * \param lambda Regularization parameter for smoothing matrix
00158     * \param withFRT include the Funk-Radon-Transformation?
00159     * \return Transformation matrix
00160     */
00161     static WMatrix< T > getSHFittingMatrix( const std::vector< WMatrixFixed< T, 3, 1 > >& orientations,
00162                                                int order,
00163                                                T lambda,
00164                                                bool withFRT );
00165 
00166     /**
00167     * This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper. The orientations are given as WUnitSphereCoordinates< T >.
00168     * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
00169     * \param order The order of the spherical harmonics intended to create
00170     * \param lambda Regularization parameter for smoothing matrix
00171     * \param withFRT include the Funk-Radon-Transformation?
00172     * \return Transformation matrix
00173     */
00174     static WMatrix< T > getSHFittingMatrix( const std::vector< WUnitSphereCoordinates< T > >& orientations,
00175                                                int order,
00176                                                T lambda,
00177                                                bool withFRT );
00178 
00179     /**
00180     * This calculates the transformation/fitting matrix T like in the 2010 Aganj paper. The orientations are given as WUnitSphereCoordinates< T >.
00181     * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
00182     * \param order The order of the spherical harmonics intended to create
00183     * \param lambda Regularization parameter for smoothing matrix
00184     * \return Transformation matrix
00185     */
00186     static WMatrix< T > getSHFittingMatrixForConstantSolidAngle( const std::vector< WMatrixFixed< T, 3, 1 > >& orientations,
00187                                                                     int order,
00188                                                                     T lambda );
00189 
00190     /**
00191     * This calculates the transformation/fitting matrix T like in the 2010 Aganj paper. The orientations are given as WUnitSphereCoordinates< T >.
00192     * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
00193     * \param order The order of the spherical harmonics intended to create
00194     * \param lambda Regularization parameter for smoothing matrix
00195     * \return Transformation matrix
00196     */
00197     static WMatrix< T > getSHFittingMatrixForConstantSolidAngle( const std::vector< WUnitSphereCoordinates< T > >& orientations,
00198                                                                  int order,
00199                                                                  T lambda );
00200 
00201     /**
00202     * Calculates the base matrix B like in the dissertation of Descoteaux.
00203     * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
00204     * \param order The order of the spherical harmonics intended to create
00205     * \return The base Matrix B
00206     */
00207     static WMatrix< T > calcBaseMatrix( const std::vector< WUnitSphereCoordinates< T > >& orientations, int order );
00208 
00209     /**
00210     * Calculates the base matrix B for the complex spherical harmonics.
00211     * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
00212     * \param order The order of the spherical harmonics intended to create
00213     * \return The base Matrix B
00214     */
00215     static WMatrix< std::complex< T > > calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates< T > > const& orientations,
00216                                                                            int order );
00217     /**
00218     * Calc eigenvalues for SH elements.
00219     * \param order The order of the spherical harmonic
00220     * \return The eigenvalues of the coefficients
00221     */
00222     static WValue< T > calcEigenvalues( size_t order );
00223 
00224     /**
00225     * Calc matrix with the eigenvalues of the SH elements on its diagonal.
00226     * \param order The order of the spherical harmonic
00227     * \return The matrix with the eigenvalues of the coefficients
00228     */
00229     static WMatrix< T > calcMatrixWithEigenvalues( size_t order );
00230 
00231     /**
00232     * This calcs the smoothing matrix L from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging"
00233     * \param order The order of the spherical harmonic
00234     * \return The smoothing matrix L
00235     */
00236     static WMatrix< T > calcSmoothingMatrix( size_t order );
00237 
00238     /**
00239     * Calculates the Funk-Radon-Transformation-Matrix P from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging"
00240     * \param order The order of the spherical harmonic
00241     * \return The Funk-Radon-Matrix P
00242     */
00243     static WMatrix< T > calcFRTMatrix( size_t order );
00244 
00245    /**
00246      * Calculates a matrix that converts spherical harmonics to symmetric tensors of equal or lower order.
00247      *
00248      * \param order The order of the symmetric tensor.
00249      *
00250      * \return the conversion matrix
00251      */
00252     static WMatrix< T > calcSHToTensorSymMatrix( std::size_t order );
00253 
00254    /**
00255      * Calculates a matrix that converts spherical harmonics to symmetric tensors of equal or lower order.
00256      *
00257      * \param order The order of the symmetric tensor.
00258      * \param orientations A vector of at least (orderTensor+1) * (orderTensor+2) / 2 orientations.
00259      *
00260      * \return the conversion matrix
00261      */
00262     static WMatrix< T > calcSHToTensorSymMatrix( std::size_t order, const std::vector< WUnitSphereCoordinates< T > >& orientations );
00263 
00264     /**
00265      * Normalize this SH in place.
00266      */
00267     void normalize();
00268 
00269 protected:
00270 private:
00271     /** order of the spherical harmonic */
00272     size_t m_order;
00273 
00274     /** coefficients of the spherical harmonic */
00275     WValue< T > m_SHCoefficients;
00276 };
00277 
00278 template< typename T >
00279 WSymmetricSphericalHarmonic< T >::WSymmetricSphericalHarmonic() :
00280     m_order( 0 ),
00281     m_SHCoefficients( 0 )
00282 {
00283 }
00284 
00285 template< typename T >
00286 WSymmetricSphericalHarmonic< T >::WSymmetricSphericalHarmonic( const WValue< T >& SHCoefficients ) :
00287   m_SHCoefficients( SHCoefficients )
00288 {
00289   // calculate order from SHCoefficients.size:
00290   // - this is done by reversing the R=(l+1)*(l+2)/2 formula from the Descoteaux paper
00291   T q = std::sqrt( 0.25 + 2.0 * static_cast< T >( m_SHCoefficients.size() ) ) - 1.5;
00292   m_order = static_cast<size_t>( q );
00293 }
00294 
00295 template< typename T >
00296 WSymmetricSphericalHarmonic< T >::~WSymmetricSphericalHarmonic()
00297 {
00298 }
00299 
00300 template< typename T >
00301 T WSymmetricSphericalHarmonic< T >::getValue( T theta, T phi ) const
00302 {
00303   T result = 0.0;
00304   int j = 0;
00305   const T rootOf2 = std::sqrt( 2.0 );
00306   for( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
00307   {
00308     // m = 1 .. k
00309     for( int m = 1; m <= k; m++ )
00310     {
00311       j = ( k*k+k+2 ) / 2 + m;
00312       result += m_SHCoefficients[ j-1 ] * rootOf2 *
00313                   static_cast< T >( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
00314     }
00315     // m = 0
00316     result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
00317     // m = -k .. -1
00318     for( int m = -k; m < 0; m++ )
00319     {
00320       j = ( k*k+k+2 ) / 2 + m;
00321       result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
00322     }
00323   }
00324   return result;
00325 }
00326 
00327 template< typename T >
00328 T WSymmetricSphericalHarmonic< T >::getValue( const WUnitSphereCoordinates< T >& coordinates ) const
00329 {
00330   return getValue( coordinates.getTheta(), coordinates.getPhi() );
00331 }
00332 
00333 template< typename T >
00334 const WValue< T >& WSymmetricSphericalHarmonic< T >::getCoefficients() const
00335 {
00336   return m_SHCoefficients;
00337 }
00338 
00339 template< typename T >
00340 WValue< T > WSymmetricSphericalHarmonic< T >::getCoefficientsSchultz() const
00341 {
00342     WValue< T > res( m_SHCoefficients.size() );
00343     size_t k = 0;
00344     for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
00345     {
00346         for( int m = -l; m <= l; ++m )
00347         {
00348             res[ k ] = m_SHCoefficients[ k ];
00349             if( m < 0 && ( ( -m ) % 2 == 1 ) )
00350             {
00351                 res[ k ] *= -1.0;
00352             }
00353             else if( m > 0 && ( m % 2 == 0 ) )
00354             {
00355                 res[ k ] *= -1.0;
00356             }
00357             ++k;
00358         }
00359     }
00360     return res;
00361 }
00362 
00363 template< typename T >
00364 WValue< std::complex< T > > WSymmetricSphericalHarmonic< T >::getCoefficientsComplex() const
00365 {
00366     WValue< std::complex< T > > res( m_SHCoefficients.size() );
00367     size_t k = 0;
00368     T r = 1.0 / sqrt( 2.0 );
00369     std::complex< T > i( 0.0, -1.0 );
00370     for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
00371     {
00372         for( int m = -l; m <= l; ++m )
00373         {
00374             if( m == 0 )
00375             {
00376                 res[ k ] = m_SHCoefficients[ k ];
00377             }
00378             else if( m < 0 )
00379             {
00380                 res[ k ] += i * r * m_SHCoefficients[ k - 2 * m ];
00381                 res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
00382             }
00383             else if( m > 0 )
00384             {
00385                 res[ k ] += i * ( -m % 2 == 0 ? -r : r ) * m_SHCoefficients[ k ];
00386                 res[ k ] += r * m_SHCoefficients[ k - 2 * m ];
00387             }
00388             ++k;
00389         }
00390     }
00391     return res;
00392 }
00393 
00394 template< typename T >
00395 T WSymmetricSphericalHarmonic< T >::calcGFA( std::vector< WUnitSphereCoordinates< T > > const& orientations ) const
00396 {
00397     T n = static_cast< T >( orientations.size() );
00398     T d = 0.0;
00399     T gfa = 0.0;
00400     T mean = 0.0;
00401     std::vector< T > v( orientations.size() );
00402 
00403     for( std::size_t i = 0; i < orientations.size(); ++i )
00404     {
00405         v[ i ] = getValue( orientations[ i ] );
00406         mean += v[ i ];
00407     }
00408     mean /= n;
00409 
00410     for( std::size_t i = 0; i < orientations.size(); ++i )
00411     {
00412         T f = v[ i ];
00413         // we use gfa as a temporary here
00414         gfa += f * f;
00415         f -= mean;
00416         d += f * f;
00417     }
00418 
00419     if( gfa == 0.0 || n <= 1.0 )
00420     {
00421         return 0.0;
00422     }
00423     // this is the real gfa
00424     gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
00425 
00426     WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
00427     if( gfa < 0.0 )
00428     {
00429         return 0.0;
00430     }
00431     else if( gfa > 1.0 )
00432     {
00433         return 1.0;
00434     }
00435     return gfa;
00436 }
00437 
00438 template< typename T >
00439 T WSymmetricSphericalHarmonic< T >::calcGFA( WMatrix< T > const& B ) const
00440 {
00441     // sh coeffs
00442     WMatrix< T > s( B.getNbCols(), 1 );
00443     WAssert( B.getNbCols() == m_SHCoefficients.size(), "" );
00444     for( std::size_t k = 0; k < B.getNbCols(); ++k )
00445     {
00446         s( k, 0 ) = m_SHCoefficients[ k ];
00447     }
00448     s = B * s;
00449     WAssert( s.getNbRows() == B.getNbRows(), "" );
00450     WAssert( s.getNbCols() == 1u, "" );
00451 
00452     T n = static_cast< T >( s.getNbRows() );
00453     T d = 0.0;
00454     T gfa = 0.0;
00455     T mean = 0.0;
00456     for( std::size_t i = 0; i < s.getNbRows(); ++i )
00457     {
00458         mean += s( i, 0 );
00459     }
00460     mean /= n;
00461 
00462     for( std::size_t i = 0; i < s.getNbRows(); ++i )
00463     {
00464         T f = s( i, 0 );
00465         // we use gfa as a temporary here
00466         gfa += f * f;
00467         f -= mean;
00468         d += f * f;
00469     }
00470 
00471     if( gfa == 0.0 || n <= 1.0 )
00472     {
00473         return 0.0;
00474     }
00475     // this is the real gfa
00476     gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
00477 
00478     WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
00479     if( gfa < 0.0 )
00480     {
00481         return 0.0;
00482     }
00483     else if( gfa > 1.0 )
00484     {
00485         return 1.0;
00486     }
00487     return gfa;
00488 }
00489 
00490 template< typename T >
00491 void WSymmetricSphericalHarmonic< T >::applyFunkRadonTransformation( WMatrix< T > const& frtMat )
00492 {
00493     WAssert( frtMat.getNbCols() == m_SHCoefficients.size(), "" );
00494     WAssert( frtMat.getNbRows() == m_SHCoefficients.size(), "" );
00495     // Funk-Radon-Transformation as in Descoteaux's thesis
00496     for( size_t j = 0; j < m_SHCoefficients.size(); j++ )
00497     {
00498         m_SHCoefficients[ j ] *= frtMat( j, j );
00499     }
00500 }
00501 
00502 template< typename T >
00503 size_t WSymmetricSphericalHarmonic< T >::getOrder() const
00504 {
00505   return m_order;
00506 }
00507 
00508 template< typename T >
00509 WMatrix< T > WSymmetricSphericalHarmonic< T >::getSHFittingMatrix( const std::vector< WMatrixFixed< T, 3, 1 > >& orientations,
00510                                                                         int order,
00511                                                                         T lambda,
00512                                                                         bool withFRT )
00513 {
00514     // convert euclid 3d orientations/gradients to sphere coordinates
00515     std::vector< WUnitSphereCoordinates< T > > sphereCoordinates;
00516     for( typename std::vector< WMatrixFixed< T, 3, 1 > >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
00517     {
00518         sphereCoordinates.push_back( WUnitSphereCoordinates< T >( *it ) );
00519     }
00520     return WSymmetricSphericalHarmonic< T >::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
00521 }
00522 
00523 template< typename T >
00524 WMatrix< T > WSymmetricSphericalHarmonic< T >::getSHFittingMatrix( const std::vector< WUnitSphereCoordinates< T > >& orientations,
00525                                                                  int order,
00526                                                                  T lambda,
00527                                                                  bool withFRT )
00528 {
00529   WMatrix< T > B( WSymmetricSphericalHarmonic< T >::calcBaseMatrix( orientations, order ) );
00530   WMatrix< T > Bt( B.transposed() );
00531   WMatrix< T > result( Bt*B );
00532   if( lambda != 0.0 )
00533   {
00534     WMatrix< T > L( WSymmetricSphericalHarmonic< T >::calcSmoothingMatrix( order ) );
00535     result += lambda*L;
00536   }
00537 
00538   result = pseudoInverse( result )*Bt;
00539   if( withFRT )
00540   {
00541     WMatrix< T > P( WSymmetricSphericalHarmonic< T >::calcFRTMatrix( order ) );
00542     return ( P * result );
00543   }
00544   return result;
00545 }
00546 
00547 template< typename T >
00548 WMatrix< T > WSymmetricSphericalHarmonic< T >::getSHFittingMatrixForConstantSolidAngle( const std::vector< WMatrixFixed< T, 3, 1 > >& orientations,
00549                                                                                       int order,
00550                                                                                       T lambda )
00551 {
00552     // convert euclid 3d orientations/gradients to sphere coordinates
00553     std::vector< WUnitSphereCoordinates< T > > sphereCoordinates;
00554     for( typename std::vector< WMatrixFixed< T, 3, 1 > >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
00555     {
00556         sphereCoordinates.push_back( WUnitSphereCoordinates< T >( *it ) );
00557     }
00558     return WSymmetricSphericalHarmonic< T >::getSHFittingMatrixForConstantSolidAngle( sphereCoordinates, order, lambda );
00559 }
00560 
00561 template< typename T >
00562 WMatrix< T > WSymmetricSphericalHarmonic< T >::getSHFittingMatrixForConstantSolidAngle(
00563     const std::vector< WUnitSphereCoordinates< T > >& orientations,
00564     int order,
00565     T lambda )
00566 {
00567   WMatrix< T > B( WSymmetricSphericalHarmonic< T >::calcBaseMatrix( orientations, order ) );
00568   WMatrix< T > Bt( B.transposed() );
00569   WMatrix< T > result( Bt*B );
00570 
00571   // regularisation
00572   if( lambda != 0.0 )
00573   {
00574     WMatrix< T > L( WSymmetricSphericalHarmonic< T >::calcSmoothingMatrix( order ) );
00575     result += lambda*L;
00576   }
00577 
00578   result = pseudoInverse( result )*Bt;
00579 
00580   // multiply with eigenvalues of coefficients / Laplace-Beltrami operator
00581   WMatrix< T > LB( WSymmetricSphericalHarmonic< T >::calcMatrixWithEigenvalues( order ) );
00582   wlog::debug( "" ) << "LB: " << LB;
00583   result = LB*result;
00584   wlog::debug( "" ) << "LB*result: " << result;
00585 
00586   // apply FRT
00587   WMatrix< T > P( WSymmetricSphericalHarmonic< T >::calcFRTMatrix( order ) );
00588   result = P * result;
00589   wlog::debug( "" ) << "P: " << P;
00590   wlog::debug( "" ) << "P*result: " << result;
00591 
00592   // correction factor
00593   T correctionFactor = 1.0 / ( 16.0 * std::pow( static_cast< T >( piDouble ), 2 ) );
00594   result *= correctionFactor;
00595 
00596   return result;
00597 }
00598 
00599 
00600 template< typename T >
00601 WMatrix< T > WSymmetricSphericalHarmonic< T >::calcBaseMatrix( const std::vector< WUnitSphereCoordinates< T > >& orientations,
00602                                                                     int order )
00603 {
00604   WMatrix< T > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
00605 
00606   // calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
00607   int j = 0;
00608   const T rootOf2 = std::sqrt( 2.0 );
00609   for(uint32_t row = 0; row < orientations.size(); row++ )
00610   {
00611     const T theta = orientations[row].getTheta();
00612     const T phi = orientations[row].getPhi();
00613     for( int k = 0; k <= order; k+=2 )
00614     {
00615       // m = 1 .. k
00616       for( int m = 1; m <= k; m++ )
00617       {
00618         j = ( k * k + k + 2 ) / 2 + m;
00619         B( row, j-1 ) = rootOf2 * static_cast< T >( std::pow( static_cast< T >( -1 ), m + 1 ) )
00620                         * boost::math::spherical_harmonic_i( k, m, theta, phi );
00621       }
00622       // m = 0
00623       B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
00624       // m = -k .. -1
00625       for( int m = -k; m < 0; m++ )
00626       {
00627         j = ( k * k + k + 2 ) / 2 + m;
00628         B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
00629       }
00630     }
00631   }
00632   return B;
00633 }
00634 
00635 template< typename T >
00636 WMatrix< std::complex< T > >
00637 WSymmetricSphericalHarmonic< T >::calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates< T > > const& orientations, int order )
00638 {
00639     WMatrix< std::complex< T > > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
00640 
00641     for( std::size_t row = 0; row < orientations.size(); row++ )
00642     {
00643         const T theta = orientations[ row ].getTheta();
00644         const T phi = orientations[ row ].getPhi();
00645 
00646         int j = 0;
00647         for( int k = 0; k <= order; k += 2 )
00648         {
00649             for( int m = -k; m < 0; m++ )
00650             {
00651                 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
00652                 ++j;
00653             }
00654             B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
00655             ++j;
00656             for( int m = 1; m <= k; m++ )
00657             {
00658                 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
00659                 ++j;
00660             }
00661         }
00662     }
00663     return B;
00664 }
00665 
00666 template< typename T >
00667 WValue< T > WSymmetricSphericalHarmonic< T >::calcEigenvalues( size_t order )
00668 {
00669     size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
00670     std::size_t i = 0;
00671     WValue< T > result( R );
00672     for( size_t k = 0; k <= order; k += 2 )
00673     {
00674         for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
00675         {
00676             result[ i ] = -static_cast< T > ( k * ( k + 1 ) );
00677             ++i;
00678         }
00679     }
00680     return result;
00681 }
00682 
00683 template< typename T >
00684 WMatrix< T > WSymmetricSphericalHarmonic< T >::calcMatrixWithEigenvalues( size_t order )
00685 {
00686     WValue< T > eigenvalues( WSymmetricSphericalHarmonic< T >::calcEigenvalues( order ) );
00687     WMatrix< T > L( eigenvalues.size(), eigenvalues.size() );
00688     for( std::size_t i = 0; i < eigenvalues.size(); ++i )
00689     {
00690         L( i, i ) = eigenvalues[ i ];
00691     }
00692     return L;
00693 }
00694 
00695 template< typename T >
00696 WMatrix< T > WSymmetricSphericalHarmonic< T >::calcSmoothingMatrix( size_t order )
00697 {
00698     WValue< T > eigenvalues( WSymmetricSphericalHarmonic< T >::calcEigenvalues( order ) );
00699     WMatrix< T > L( eigenvalues.size(), eigenvalues.size() );
00700     for( std::size_t i = 0; i < eigenvalues.size(); ++i )
00701     {
00702         L( i, i ) = std::pow( eigenvalues[ i ], 2 );
00703     }
00704     return L;
00705 }
00706 
00707 template< typename T >
00708 WMatrix< T > WSymmetricSphericalHarmonic< T >::calcFRTMatrix( size_t order )
00709 {
00710     size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
00711     std::size_t i = 0;
00712     WMatrix< T > result( R, R );
00713     for( size_t k = 0; k <= order; k += 2 )
00714     {
00715         T h = 2.0 * static_cast< T >( piDouble ) * static_cast< T >( std::pow( static_cast< T >( -1 ), static_cast< T >( k / 2 ) ) ) *
00716                     static_cast< T >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast< T >( evenFactorial( k ) );
00717         for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
00718         {
00719             result( i, i ) = h;
00720             ++i;
00721         }
00722     }
00723     return result;
00724 }
00725 
00726 template< typename T >
00727 WMatrix< T > WSymmetricSphericalHarmonic< T >::calcSHToTensorSymMatrix( std::size_t order )
00728 {
00729     std::vector< WVector3d > vertices;
00730     std::vector< unsigned int > triangles;
00731     // calc necessary tesselation level
00732     int level = static_cast< int >( log( static_cast< float >( ( order + 1 ) * ( order + 2 ) ) ) / 2.0f + 1.0f );
00733     tesselateIcosahedron( &vertices, &triangles, level );
00734     std::vector< WUnitSphereCoordinates< T > > orientations;
00735     for( typename std::vector< WVector3d >::const_iterator cit = vertices.begin(); cit != vertices.end(); cit++ )
00736     {
00737         // avoid linear dependent orientations
00738         if( ( *cit )[ 0 ] >= 0.0001 )
00739         {
00740             orientations.push_back( WUnitSphereCoordinates< T >( WMatrixFixed< T, 3, 1 >( *cit ) ) );
00741         }
00742     }
00743     return WSymmetricSphericalHarmonic< T >::calcSHToTensorSymMatrix( order, orientations );
00744 }
00745 
00746 template< typename T >
00747 WMatrix< T > WSymmetricSphericalHarmonic< T >::calcSHToTensorSymMatrix( std::size_t order,
00748                                                                         const std::vector< WUnitSphereCoordinates< T > >& orientations )
00749 {
00750     std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
00751     WPrecondEquals( order % 2, 0u );
00752     WPrecondLess( numElements, orientations.size() + 1 );
00753 
00754     // store first numElements orientations as 3d-vectors
00755     std::vector< WMatrixFixed< T, 3, 1 > > directions( numElements );
00756     for( std::size_t i = 0; i < numElements; ++i )
00757     {
00758         directions[ i ] = orientations[ i ].getEuclidean();
00759     }
00760 
00761     // initialize an index array
00762     std::vector< std::size_t > indices( order, 0 );
00763 
00764     // calc tensor evaluation matrix
00765     WMatrix< T > TEMat( numElements, numElements );
00766     for( std::size_t j = 0; j < numElements; ++j ) // j is the 'permutation index'
00767     {
00768         // stores how often each value is represented in the index array
00769         std::size_t amount[ 3 ] = { 0, 0, 0 };
00770         for( std::size_t k = 0; k < order; ++k )
00771         {
00772             ++amount[ indices[ k ] ];
00773         }
00774 
00775         // from WTensorSym.h
00776         std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
00777         for( std::size_t i = 0; i < numElements; ++i ) // i is the 'direction index'
00778         {
00779             TEMat( i, j ) = multiplicity;
00780             for( std::size_t k = 0; k < order; ++k )
00781             {
00782                 TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
00783             }
00784         }
00785 
00786         // from TensorBase.h
00787         positionIterateSortedOneStep( order, 3, indices );
00788     }
00789     directions.clear();
00790 
00791     // we do not want more orientations than nessessary
00792     std::vector< WUnitSphereCoordinates< T > > ori2( orientations.begin(), orientations.begin() + numElements );
00793 
00794     WMatrix< T > p = pseudoInverse( TEMat );
00795 
00796     WAssert( ( TEMat*p ).isIdentity( 1e-3 ), "Test of inverse matrix failed. Are the given orientations linear independent?" );
00797 
00798     return p * calcBaseMatrix( ori2, order );
00799 }
00800 
00801 template< typename T >
00802 void WSymmetricSphericalHarmonic< T >::normalize()
00803 {
00804   T scale = 0.0;
00805   if( m_SHCoefficients.size() > 0 )
00806   {
00807       scale = std::sqrt( 4.0 * static_cast< T >( piDouble ) ) * m_SHCoefficients[0];
00808   }
00809   for( size_t i = 0; i < m_SHCoefficients.size(); i++ )
00810   {
00811       m_SHCoefficients[ i ] /= scale;
00812   }
00813 }
00814 
00815 #endif  // WSYMMETRICSPHERICALHARMONIC_H