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