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