OpenWalnut
1.4.0
|
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