OpenWalnut 1.3.1
WSymmetricSphericalHarmonic.cpp
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 #include <stdint.h>
00026 
00027 #include <vector>
00028 
00029 #include <boost/math/special_functions/spherical_harmonic.hpp>
00030 
00031 #include "core/common/WLogger.h"
00032 #include "../exceptions/WPreconditionNotMet.h"
00033 #include "linearAlgebra/WLinearAlgebra.h"
00034 #include "WLinearAlgebraFunctions.h"
00035 #include "WMatrix.h"
00036 #include "WTensorSym.h"
00037 #include "WUnitSphereCoordinates.h"
00038 
00039 #include "WSymmetricSphericalHarmonic.h"
00040 
00041 WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic() :
00042     m_order( 0 ),
00043     m_SHCoefficients( 0 )
00044 {
00045 }
00046 
00047 WSymmetricSphericalHarmonic::WSymmetricSphericalHarmonic( const WValue<double>& SHCoefficients ) :
00048   m_SHCoefficients( SHCoefficients )
00049 {
00050   // calculate order from SHCoefficients.size:
00051   // - this is done by reversing the R=(l+1)*(l+2)/2 formula from the Descoteaux paper
00052   double q = std::sqrt( 0.25 + 2.0 * static_cast<double>( m_SHCoefficients.size() ) ) - 1.5;
00053   m_order = static_cast<size_t>( q );
00054 }
00055 
00056 WSymmetricSphericalHarmonic::~WSymmetricSphericalHarmonic()
00057 {
00058 }
00059 
00060 double WSymmetricSphericalHarmonic::getValue( double theta, double phi ) const
00061 {
00062   double result = 0.0;
00063   int j = 0;
00064   const double rootOf2 = std::sqrt( 2.0 );
00065   for( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
00066   {
00067     // m = 1 .. k
00068     for( int m = 1; m <= k; m++ )
00069     {
00070       j = ( k*k+k+2 ) / 2 + m;
00071       result += m_SHCoefficients[ j-1 ] * rootOf2 *
00072                   static_cast<double>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
00073     }
00074     // m = 0
00075     result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
00076     // m = -k .. -1
00077     for( int m = -k; m < 0; m++ )
00078     {
00079       j = ( k*k+k+2 ) / 2 + m;
00080       result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
00081     }
00082   }
00083   return result;
00084 }
00085 
00086 double WSymmetricSphericalHarmonic::getValue( const WUnitSphereCoordinates& coordinates ) const
00087 {
00088   return getValue( coordinates.getTheta(), coordinates.getPhi() );
00089 }
00090 
00091 const WValue<double>& WSymmetricSphericalHarmonic::getCoefficients() const
00092 {
00093   return m_SHCoefficients;
00094 }
00095 
00096 WValue< double > WSymmetricSphericalHarmonic::getCoefficientsSchultz() const
00097 {
00098     WValue< double > res( m_SHCoefficients.size() );
00099     size_t k = 0;
00100     for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
00101     {
00102         for( int m = -l; m <= l; ++m )
00103         {
00104             res[ k ] = m_SHCoefficients[ k ];
00105             if( m < 0 && ( ( -m ) % 2 == 1 ) )
00106             {
00107                 res[ k ] *= -1.0;
00108             }
00109             else if( m > 0 && ( m % 2 == 0 ) )
00110             {
00111                 res[ k ] *= -1.0;
00112             }
00113             ++k;
00114         }
00115     }
00116     return res;
00117 }
00118 
00119 WValue< std::complex< double > > WSymmetricSphericalHarmonic::getCoefficientsComplex() const
00120 {
00121     WValue< std::complex< double > > res( m_SHCoefficients.size() );
00122     size_t k = 0;
00123     double r = 1.0 / sqrt( 2.0 );
00124     std::complex< double > i( 0.0, -1.0 );
00125     for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
00126     {
00127         for( int m = -l; m <= l; ++m )
00128         {
00129             if( m == 0 )
00130             {
00131                 res[ k ] = m_SHCoefficients[ k ];
00132             }
00133             else if( m < 0 )
00134             {
00135                 res[ k ] += i * r * m_SHCoefficients[ k - 2 * m ];
00136                 res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
00137             }
00138             else if( m > 0 )
00139             {
00140                 res[ k ] += i * ( -m % 2 == 0 ? -r : r ) * m_SHCoefficients[ k ];
00141                 res[ k ] += r * m_SHCoefficients[ k - 2 * m ];
00142             }
00143             ++k;
00144         }
00145     }
00146     return res;
00147 }
00148 
00149 double WSymmetricSphericalHarmonic::calcGFA( std::vector< WUnitSphereCoordinates > const& orientations ) const
00150 {
00151     double n = static_cast< double >( orientations.size() );
00152     double d = 0.0;
00153     double gfa = 0.0;
00154     double mean = 0.0;
00155     double v[ 15 ];
00156 
00157     for( std::size_t i = 0; i < orientations.size(); ++i )
00158     {
00159         v[ i ] = getValue( orientations[ i ] );
00160         mean += v[ i ];
00161     }
00162     mean /= n;
00163 
00164     for( std::size_t i = 0; i < orientations.size(); ++i )
00165     {
00166         double f = v[ i ];
00167         // we use gfa as a temporary here
00168         gfa += f * f;
00169         f -= mean;
00170         d += f * f;
00171     }
00172 
00173     if( gfa == 0.0 || n <= 1.0 )
00174     {
00175         return 0.0;
00176     }
00177     // this is the real gfa
00178     gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
00179 
00180     WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
00181     if( gfa < 0.0 )
00182     {
00183         return 0.0;
00184     }
00185     else if( gfa > 1.0 )
00186     {
00187         return 1.0;
00188     }
00189     return gfa;
00190 }
00191 
00192 double WSymmetricSphericalHarmonic::calcGFA( WMatrix< double > const& B ) const
00193 {
00194     // sh coeffs
00195     WMatrix< double > s( B.getNbCols(), 1 );
00196     WAssert( B.getNbCols() == m_SHCoefficients.size(), "" );
00197     for( std::size_t k = 0; k < B.getNbCols(); ++k )
00198     {
00199         s( k, 0 ) = m_SHCoefficients[ k ];
00200     }
00201     s = B * s;
00202     WAssert( s.getNbRows() == B.getNbRows(), "" );
00203     WAssert( s.getNbCols() == 1u, "" );
00204 
00205     double n = static_cast< double >( s.getNbRows() );
00206     double d = 0.0;
00207     double gfa = 0.0;
00208     double mean = 0.0;
00209     for( std::size_t i = 0; i < s.getNbRows(); ++i )
00210     {
00211         mean += s( i, 0 );
00212     }
00213     mean /= n;
00214 
00215     for( std::size_t i = 0; i < s.getNbRows(); ++i )
00216     {
00217         double f = s( i, 0 );
00218         // we use gfa as a temporary here
00219         gfa += f * f;
00220         f -= mean;
00221         d += f * f;
00222     }
00223 
00224     if( gfa == 0.0 || n <= 1.0 )
00225     {
00226         return 0.0;
00227     }
00228     // this is the real gfa
00229     gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
00230 
00231     WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
00232     if( gfa < 0.0 )
00233     {
00234         return 0.0;
00235     }
00236     else if( gfa > 1.0 )
00237     {
00238         return 1.0;
00239     }
00240     return gfa;
00241 }
00242 
00243 void WSymmetricSphericalHarmonic::applyFunkRadonTransformation( WMatrix< double > const& frtMat )
00244 {
00245     WAssert( frtMat.getNbCols() == m_SHCoefficients.size(), "" );
00246     WAssert( frtMat.getNbRows() == m_SHCoefficients.size(), "" );
00247     // Funk-Radon-Transformation as in Descoteaux's thesis
00248     for( size_t j = 0; j < m_SHCoefficients.size(); j++ )
00249     {
00250         m_SHCoefficients[ j ] *= frtMat( j, j );
00251     }
00252 }
00253 
00254 size_t WSymmetricSphericalHarmonic::getOrder() const
00255 {
00256   return m_order;
00257 }
00258 
00259 WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WVector3d >& orientations,
00260                                                                         int order,
00261                                                                         double lambda,
00262                                                                         bool withFRT )
00263 {
00264   // convert euclid 3d orientations/gradients to sphere coordinates
00265   std::vector< WUnitSphereCoordinates > sphereCoordinates;
00266   for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
00267   {
00268     sphereCoordinates.push_back( WUnitSphereCoordinates( *it ) );
00269   }
00270   return WSymmetricSphericalHarmonic::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
00271 }
00272 
00273 WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
00274                                                                  int order,
00275                                                                  double lambda,
00276                                                                  bool withFRT )
00277 {
00278   WMatrix<double> B( WSymmetricSphericalHarmonic::calcBaseMatrix( orientations, order ) );
00279   WMatrix<double> Bt( B.transposed() );
00280   WMatrix<double> result( Bt*B );
00281   if( lambda != 0.0 )
00282   {
00283     WMatrix<double> L( WSymmetricSphericalHarmonic::calcSmoothingMatrix( order ) );
00284     result += lambda*L;
00285   }
00286 
00287   result = pseudoInverse( result )*Bt;
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::getSHFittingMatrixForConstantSolidAngle( const std::vector< WVector3d >& orientations,
00297                                                                                       int order,
00298                                                                                       double lambda )
00299 {
00300   // convert euclid 3d orientations/gradients to sphere coordinates
00301   std::vector< WUnitSphereCoordinates > sphereCoordinates;
00302   for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
00303   {
00304     sphereCoordinates.push_back( WUnitSphereCoordinates( *it ) );
00305   }
00306   return WSymmetricSphericalHarmonic::getSHFittingMatrixForConstantSolidAngle( sphereCoordinates, order, lambda );
00307 }
00308 
00309 WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrixForConstantSolidAngle(
00310     const std::vector< WUnitSphereCoordinates >& orientations,
00311     int order,
00312     double lambda )
00313 {
00314   WMatrix<double> B( WSymmetricSphericalHarmonic::calcBaseMatrix( orientations, order ) );
00315   WMatrix<double> Bt( B.transposed() );
00316   WMatrix<double> result( Bt*B );
00317 
00318   // regularisation
00319   if( lambda != 0.0 )
00320   {
00321     WMatrix<double> L( WSymmetricSphericalHarmonic::calcSmoothingMatrix( order ) );
00322     result += lambda*L;
00323   }
00324 
00325   result = pseudoInverse( result )*Bt;
00326 
00327   // multiply with eigenvalues of coefficients / Laplace-Beltrami operator
00328   WMatrix<double> LB( WSymmetricSphericalHarmonic::calcMatrixWithEigenvalues( order ) );
00329   wlog::debug( "" ) << "LB: " << LB;
00330   result = LB*result;
00331   wlog::debug( "" ) << "LB*result: " << result;
00332 
00333   // apply FRT
00334   WMatrix<double> P( WSymmetricSphericalHarmonic::calcFRTMatrix( order ) );
00335   result = P * result;
00336   wlog::debug( "" ) << "P: " << P;
00337   wlog::debug( "" ) << "P*result: " << result;
00338 
00339   // correction factor
00340   double correctionFactor = 1.0 / ( 16.0 * std::pow( piDouble, 2 ) );
00341   result *= correctionFactor;
00342 
00343   return result;
00344 }
00345 
00346 
00347 WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
00348                                                                     int order )
00349 {
00350   WMatrix<double> B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
00351 
00352   // calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
00353   int j = 0;
00354   const double rootOf2 = std::sqrt( 2.0 );
00355   for(uint32_t row = 0; row < orientations.size(); row++ )
00356   {
00357     const double theta = orientations[row].getTheta();
00358     const double phi = orientations[row].getPhi();
00359     for( int k = 0; k <= order; k+=2 )
00360     {
00361       // m = 1 .. k
00362       for( int m = 1; m <= k; m++ )
00363       {
00364         j = ( k * k + k + 2 ) / 2 + m;
00365         B( row, j-1 ) = rootOf2 * static_cast<double>( std::pow( static_cast<double>( -1 ), m + 1 ) )
00366                         * boost::math::spherical_harmonic_i( k, m, theta, phi );
00367       }
00368       // m = 0
00369       B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
00370       // m = -k .. -1
00371       for( int m = -k; m < 0; m++ )
00372       {
00373         j = ( k * k + k + 2 ) / 2 + m;
00374         B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
00375       }
00376     }
00377   }
00378   return B;
00379 }
00380 
00381 WMatrix< std::complex< double > >
00382 WSymmetricSphericalHarmonic::calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations, int order )
00383 {
00384     WMatrix< std::complex< double > > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
00385 
00386     for( std::size_t row = 0; row < orientations.size(); row++ )
00387     {
00388         const double theta = orientations[ row ].getTheta();
00389         const double phi = orientations[ row ].getPhi();
00390 
00391         int j = 0;
00392         for( int k = 0; k <= order; k += 2 )
00393         {
00394             for( int m = -k; m < 0; m++ )
00395             {
00396                 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
00397                 ++j;
00398             }
00399             B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
00400             ++j;
00401             for( int m = 1; m <= k; m++ )
00402             {
00403                 B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
00404                 ++j;
00405             }
00406         }
00407     }
00408     return B;
00409 }
00410 
00411 WValue<double> WSymmetricSphericalHarmonic::calcEigenvalues( size_t order )
00412 {
00413     size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
00414     std::size_t i = 0;
00415     WValue<double> result( R );
00416     for( size_t k = 0; k <= order; k += 2 )
00417     {
00418         for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
00419         {
00420             result[ i ] = -static_cast< double > ( k * ( k + 1 ) );
00421             ++i;
00422         }
00423     }
00424     return result;
00425 }
00426 
00427 WMatrix<double> WSymmetricSphericalHarmonic::calcMatrixWithEigenvalues( size_t order )
00428 {
00429     WValue<double> eigenvalues( WSymmetricSphericalHarmonic::calcEigenvalues( order ) );
00430     WMatrix<double> L( eigenvalues.size(), eigenvalues.size() );
00431     for( std::size_t i = 0; i < eigenvalues.size(); ++i )
00432     {
00433         L( i, i ) = eigenvalues[ i ];
00434     }
00435     return L;
00436 }
00437 
00438 WMatrix<double> WSymmetricSphericalHarmonic::calcSmoothingMatrix( size_t order )
00439 {
00440     WValue<double> eigenvalues( WSymmetricSphericalHarmonic::calcEigenvalues( order ) );
00441     WMatrix<double> L( eigenvalues.size(), eigenvalues.size() );
00442     for( std::size_t i = 0; i < eigenvalues.size(); ++i )
00443     {
00444         L( i, i ) = std::pow( eigenvalues[ i ], 2 );
00445     }
00446     return L;
00447 }
00448 
00449 WMatrix<double> WSymmetricSphericalHarmonic::calcFRTMatrix( size_t order )
00450 {
00451     size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
00452     std::size_t i = 0;
00453     WMatrix< double > result( R, R );
00454     for( size_t k = 0; k <= order; k += 2 )
00455     {
00456         double h = 2.0 * piDouble * static_cast< double >( std::pow( static_cast< double >( -1 ), static_cast< double >( k / 2 ) ) ) *
00457                     static_cast< double >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast<double>( evenFactorial( k ) );
00458         for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
00459         {
00460             result( i, i ) = h;
00461             ++i;
00462         }
00463     }
00464     return result;
00465 }
00466 
00467 WMatrix< double > WSymmetricSphericalHarmonic::calcSHToTensorSymMatrix( std::size_t order,
00468                                                                                const std::vector< WUnitSphereCoordinates >& orientations )
00469 {
00470     std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
00471     WPrecondEquals( order % 2, 0u );
00472     WPrecondLess( numElements, orientations.size() + 1 );
00473 
00474     // store first numElements orientations as 3d-vectors
00475     std::vector< WVector3d > directions( numElements );
00476     for( std::size_t i = 0; i < numElements; ++i )
00477     {
00478         directions[ i ] = orientations[ i ].getEuclidean();
00479     }
00480 
00481     // initialize an index array
00482     std::vector< std::size_t > indices( order, 0 );
00483 
00484     // calc tensor evaluation matrix
00485     WMatrix< double > TEMat( numElements, numElements );
00486     for( std::size_t j = 0; j < numElements; ++j ) // j is the 'permutation index'
00487     {
00488         // stores how often each value is represented in the index array
00489         std::size_t amount[ 3 ] = { 0, 0, 0 };
00490         for( std::size_t k = 0; k < order; ++k )
00491         {
00492             ++amount[ indices[ k ] ];
00493         }
00494 
00495         // from WTensorSym.h
00496         std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
00497         for( std::size_t i = 0; i < numElements; ++i ) // i is the 'direction index'
00498         {
00499             TEMat( i, j ) = multiplicity;
00500             for( std::size_t k = 0; k < order; ++k )
00501             {
00502                 TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
00503             }
00504         }
00505 
00506         // from TensorBase.h
00507         positionIterateSortedOneStep( order, 3, indices );
00508     }
00509     directions.clear();
00510 
00511     // we do not want more orientations than nessessary
00512     std::vector< WUnitSphereCoordinates > ori2( orientations.begin(), orientations.begin() + numElements );
00513 
00514     WMatrix< double > p = pseudoInverse( TEMat );
00515 
00516     return p * calcBaseMatrix( ori2, order );
00517 }
00518 
00519 void WSymmetricSphericalHarmonic::normalize()
00520 {
00521   double scale = 0.0;
00522   if( m_SHCoefficients.size() > 0 )
00523     scale = std::sqrt( 4.0 * piDouble ) * m_SHCoefficients[0];
00524   for( size_t i = 0; i < m_SHCoefficients.size(); i++ )
00525   {
00526     m_SHCoefficients[ i ] /= scale;
00527   }
00528 }
00529 // NOLINT