OpenWalnut  1.2.5
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends
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 "../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   // calculate order from SHCoefficients.size:
00050   // - this is done by reversing the R=(l+1)*(l+2)/2 formula from the Descoteaux paper
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     // m = 1 .. k
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     // m = 0
00074     result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
00075     // m = -k .. -1
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         // we use gfa as a temporary here
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     // this is the real gfa
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     // sh coeffs
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         // we use gfa as a temporary here
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     // this is the real gfa
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     // Funk-Radon-Transformation as in Descoteaux's thesis
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   // convert euclid 3d orientations/gradients to sphere coordinates
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 //   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::calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
00297                                                                     int order )
00298 {
00299   WMatrix<double> B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
00300 
00301   // calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
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       // m = 1 .. k
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       // m = 0
00318       B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
00319       // m = -k .. -1
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     // store first numElements orientations as 3d-vectors
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     // initialize an index array
00409     std::vector< std::size_t > indices( order, 0 );
00410 
00411     // calc tensor evaluation matrix
00412     WMatrix< double > TEMat( numElements, numElements );
00413     for( std::size_t j = 0; j < numElements; ++j ) // j is the 'permutation index'
00414     {
00415         // stores how often each value is represented in the index array
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         // from WTensorSym.h
00423         std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
00424         for( std::size_t i = 0; i < numElements; ++i ) // i is the 'direction index'
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         // from TensorBase.h
00434         positionIterateSortedOneStep( order, 3, indices );
00435     }
00436     directions.clear();
00437 
00438     // we do not want more orientations than nessessary
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 // NOLINT
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends