OpenWalnut  1.3.1
WSymmetricSphericalHarmonic.cpp
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #include <stdint.h>
26 
27 #include <vector>
28 
29 #include <boost/math/special_functions/spherical_harmonic.hpp>
30 
31 #include "core/common/WLogger.h"
32 #include "../exceptions/WPreconditionNotMet.h"
33 #include "linearAlgebra/WLinearAlgebra.h"
34 #include "WLinearAlgebraFunctions.h"
35 #include "WMatrix.h"
36 #include "WTensorSym.h"
37 #include "WUnitSphereCoordinates.h"
38 
39 #include "WSymmetricSphericalHarmonic.h"
40 
42  m_order( 0 ),
43  m_SHCoefficients( 0 )
44 {
45 }
46 
48  m_SHCoefficients( SHCoefficients )
49 {
50  // calculate order from SHCoefficients.size:
51  // - this is done by reversing the R=(l+1)*(l+2)/2 formula from the Descoteaux paper
52  double q = std::sqrt( 0.25 + 2.0 * static_cast<double>( m_SHCoefficients.size() ) ) - 1.5;
53  m_order = static_cast<size_t>( q );
54 }
55 
57 {
58 }
59 
60 double WSymmetricSphericalHarmonic::getValue( double theta, double phi ) const
61 {
62  double result = 0.0;
63  int j = 0;
64  const double rootOf2 = std::sqrt( 2.0 );
65  for( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
66  {
67  // m = 1 .. k
68  for( int m = 1; m <= k; m++ )
69  {
70  j = ( k*k+k+2 ) / 2 + m;
71  result += m_SHCoefficients[ j-1 ] * rootOf2 *
72  static_cast<double>( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
73  }
74  // m = 0
75  result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
76  // m = -k .. -1
77  for( int m = -k; m < 0; m++ )
78  {
79  j = ( k*k+k+2 ) / 2 + m;
80  result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
81  }
82  }
83  return result;
84 }
85 
87 {
88  return getValue( coordinates.getTheta(), coordinates.getPhi() );
89 }
90 
92 {
93  return m_SHCoefficients;
94 }
95 
97 {
99  size_t k = 0;
100  for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
101  {
102  for( int m = -l; m <= l; ++m )
103  {
104  res[ k ] = m_SHCoefficients[ k ];
105  if( m < 0 && ( ( -m ) % 2 == 1 ) )
106  {
107  res[ k ] *= -1.0;
108  }
109  else if( m > 0 && ( m % 2 == 0 ) )
110  {
111  res[ k ] *= -1.0;
112  }
113  ++k;
114  }
115  }
116  return res;
117 }
118 
120 {
122  size_t k = 0;
123  double r = 1.0 / sqrt( 2.0 );
124  std::complex< double > i( 0.0, -1.0 );
125  for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
126  {
127  for( int m = -l; m <= l; ++m )
128  {
129  if( m == 0 )
130  {
131  res[ k ] = m_SHCoefficients[ k ];
132  }
133  else if( m < 0 )
134  {
135  res[ k ] += i * r * m_SHCoefficients[ k - 2 * m ];
136  res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
137  }
138  else if( m > 0 )
139  {
140  res[ k ] += i * ( -m % 2 == 0 ? -r : r ) * m_SHCoefficients[ k ];
141  res[ k ] += r * m_SHCoefficients[ k - 2 * m ];
142  }
143  ++k;
144  }
145  }
146  return res;
147 }
148 
149 double WSymmetricSphericalHarmonic::calcGFA( std::vector< WUnitSphereCoordinates > const& orientations ) const
150 {
151  double n = static_cast< double >( orientations.size() );
152  double d = 0.0;
153  double gfa = 0.0;
154  double mean = 0.0;
155  double v[ 15 ];
156 
157  for( std::size_t i = 0; i < orientations.size(); ++i )
158  {
159  v[ i ] = getValue( orientations[ i ] );
160  mean += v[ i ];
161  }
162  mean /= n;
163 
164  for( std::size_t i = 0; i < orientations.size(); ++i )
165  {
166  double f = v[ i ];
167  // we use gfa as a temporary here
168  gfa += f * f;
169  f -= mean;
170  d += f * f;
171  }
172 
173  if( gfa == 0.0 || n <= 1.0 )
174  {
175  return 0.0;
176  }
177  // this is the real gfa
178  gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
179 
180  WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
181  if( gfa < 0.0 )
182  {
183  return 0.0;
184  }
185  else if( gfa > 1.0 )
186  {
187  return 1.0;
188  }
189  return gfa;
190 }
191 
193 {
194  // sh coeffs
195  WMatrix< double > s( B.getNbCols(), 1 );
196  WAssert( B.getNbCols() == m_SHCoefficients.size(), "" );
197  for( std::size_t k = 0; k < B.getNbCols(); ++k )
198  {
199  s( k, 0 ) = m_SHCoefficients[ k ];
200  }
201  s = B * s;
202  WAssert( s.getNbRows() == B.getNbRows(), "" );
203  WAssert( s.getNbCols() == 1u, "" );
204 
205  double n = static_cast< double >( s.getNbRows() );
206  double d = 0.0;
207  double gfa = 0.0;
208  double mean = 0.0;
209  for( std::size_t i = 0; i < s.getNbRows(); ++i )
210  {
211  mean += s( i, 0 );
212  }
213  mean /= n;
214 
215  for( std::size_t i = 0; i < s.getNbRows(); ++i )
216  {
217  double f = s( i, 0 );
218  // we use gfa as a temporary here
219  gfa += f * f;
220  f -= mean;
221  d += f * f;
222  }
223 
224  if( gfa == 0.0 || n <= 1.0 )
225  {
226  return 0.0;
227  }
228  // this is the real gfa
229  gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
230 
231  WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
232  if( gfa < 0.0 )
233  {
234  return 0.0;
235  }
236  else if( gfa > 1.0 )
237  {
238  return 1.0;
239  }
240  return gfa;
241 }
242 
244 {
245  WAssert( frtMat.getNbCols() == m_SHCoefficients.size(), "" );
246  WAssert( frtMat.getNbRows() == m_SHCoefficients.size(), "" );
247  // Funk-Radon-Transformation as in Descoteaux's thesis
248  for( size_t j = 0; j < m_SHCoefficients.size(); j++ )
249  {
250  m_SHCoefficients[ j ] *= frtMat( j, j );
251  }
252 }
253 
255 {
256  return m_order;
257 }
258 
259 WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WVector3d >& orientations,
260  int order,
261  double lambda,
262  bool withFRT )
263 {
264  // convert euclid 3d orientations/gradients to sphere coordinates
265  std::vector< WUnitSphereCoordinates > sphereCoordinates;
266  for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
267  {
268  sphereCoordinates.push_back( WUnitSphereCoordinates( *it ) );
269  }
270  return WSymmetricSphericalHarmonic::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
271 }
272 
273 WMatrix<double> WSymmetricSphericalHarmonic::getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
274  int order,
275  double lambda,
276  bool withFRT )
277 {
279  WMatrix<double> Bt( B.transposed() );
280  WMatrix<double> result( Bt*B );
281  if( lambda != 0.0 )
282  {
284  result += lambda*L;
285  }
286 
287  result = pseudoInverse( result )*Bt;
288  if( withFRT )
289  {
291  return ( P * result );
292  }
293  return result;
294 }
295 
297  int order,
298  double lambda )
299 {
300  // convert euclid 3d orientations/gradients to sphere coordinates
301  std::vector< WUnitSphereCoordinates > sphereCoordinates;
302  for( std::vector< WVector3d >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
303  {
304  sphereCoordinates.push_back( WUnitSphereCoordinates( *it ) );
305  }
306  return WSymmetricSphericalHarmonic::getSHFittingMatrixForConstantSolidAngle( sphereCoordinates, order, lambda );
307 }
308 
310  const std::vector< WUnitSphereCoordinates >& orientations,
311  int order,
312  double lambda )
313 {
315  WMatrix<double> Bt( B.transposed() );
316  WMatrix<double> result( Bt*B );
317 
318  // regularisation
319  if( lambda != 0.0 )
320  {
322  result += lambda*L;
323  }
324 
325  result = pseudoInverse( result )*Bt;
326 
327  // multiply with eigenvalues of coefficients / Laplace-Beltrami operator
329  wlog::debug( "" ) << "LB: " << LB;
330  result = LB*result;
331  wlog::debug( "" ) << "LB*result: " << result;
332 
333  // apply FRT
335  result = P * result;
336  wlog::debug( "" ) << "P: " << P;
337  wlog::debug( "" ) << "P*result: " << result;
338 
339  // correction factor
340  double correctionFactor = 1.0 / ( 16.0 * std::pow( piDouble, 2 ) );
341  result *= correctionFactor;
342 
343  return result;
344 }
345 
346 
347 WMatrix<double> WSymmetricSphericalHarmonic::calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
348  int order )
349 {
350  WMatrix<double> B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
351 
352  // calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
353  int j = 0;
354  const double rootOf2 = std::sqrt( 2.0 );
355  for(uint32_t row = 0; row < orientations.size(); row++ )
356  {
357  const double theta = orientations[row].getTheta();
358  const double phi = orientations[row].getPhi();
359  for( int k = 0; k <= order; k+=2 )
360  {
361  // m = 1 .. k
362  for( int m = 1; m <= k; m++ )
363  {
364  j = ( k * k + k + 2 ) / 2 + m;
365  B( row, j-1 ) = rootOf2 * static_cast<double>( std::pow( static_cast<double>( -1 ), m + 1 ) )
366  * boost::math::spherical_harmonic_i( k, m, theta, phi );
367  }
368  // m = 0
369  B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
370  // m = -k .. -1
371  for( int m = -k; m < 0; m++ )
372  {
373  j = ( k * k + k + 2 ) / 2 + m;
374  B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
375  }
376  }
377  }
378  return B;
379 }
380 
382 WSymmetricSphericalHarmonic::calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations, int order )
383 {
384  WMatrix< std::complex< double > > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
385 
386  for( std::size_t row = 0; row < orientations.size(); row++ )
387  {
388  const double theta = orientations[ row ].getTheta();
389  const double phi = orientations[ row ].getPhi();
390 
391  int j = 0;
392  for( int k = 0; k <= order; k += 2 )
393  {
394  for( int m = -k; m < 0; m++ )
395  {
396  B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
397  ++j;
398  }
399  B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
400  ++j;
401  for( int m = 1; m <= k; m++ )
402  {
403  B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
404  ++j;
405  }
406  }
407  }
408  return B;
409 }
410 
412 {
413  size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
414  std::size_t i = 0;
415  WValue<double> result( R );
416  for( size_t k = 0; k <= order; k += 2 )
417  {
418  for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
419  {
420  result[ i ] = -static_cast< double > ( k * ( k + 1 ) );
421  ++i;
422  }
423  }
424  return result;
425 }
426 
428 {
430  WMatrix<double> L( eigenvalues.size(), eigenvalues.size() );
431  for( std::size_t i = 0; i < eigenvalues.size(); ++i )
432  {
433  L( i, i ) = eigenvalues[ i ];
434  }
435  return L;
436 }
437 
439 {
441  WMatrix<double> L( eigenvalues.size(), eigenvalues.size() );
442  for( std::size_t i = 0; i < eigenvalues.size(); ++i )
443  {
444  L( i, i ) = std::pow( eigenvalues[ i ], 2 );
445  }
446  return L;
447 }
448 
450 {
451  size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
452  std::size_t i = 0;
453  WMatrix< double > result( R, R );
454  for( size_t k = 0; k <= order; k += 2 )
455  {
456  double h = 2.0 * piDouble * static_cast< double >( std::pow( static_cast< double >( -1 ), static_cast< double >( k / 2 ) ) ) *
457  static_cast< double >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast<double>( evenFactorial( k ) );
458  for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
459  {
460  result( i, i ) = h;
461  ++i;
462  }
463  }
464  return result;
465 }
466 
468  const std::vector< WUnitSphereCoordinates >& orientations )
469 {
470  std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
471  WPrecondEquals( order % 2, 0u );
472  WPrecondLess( numElements, orientations.size() + 1 );
473 
474  // store first numElements orientations as 3d-vectors
475  std::vector< WVector3d > directions( numElements );
476  for( std::size_t i = 0; i < numElements; ++i )
477  {
478  directions[ i ] = orientations[ i ].getEuclidean();
479  }
480 
481  // initialize an index array
482  std::vector< std::size_t > indices( order, 0 );
483 
484  // calc tensor evaluation matrix
485  WMatrix< double > TEMat( numElements, numElements );
486  for( std::size_t j = 0; j < numElements; ++j ) // j is the 'permutation index'
487  {
488  // stores how often each value is represented in the index array
489  std::size_t amount[ 3 ] = { 0, 0, 0 };
490  for( std::size_t k = 0; k < order; ++k )
491  {
492  ++amount[ indices[ k ] ];
493  }
494 
495  // from WTensorSym.h
496  std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
497  for( std::size_t i = 0; i < numElements; ++i ) // i is the 'direction index'
498  {
499  TEMat( i, j ) = multiplicity;
500  for( std::size_t k = 0; k < order; ++k )
501  {
502  TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
503  }
504  }
505 
506  // from TensorBase.h
507  positionIterateSortedOneStep( order, 3, indices );
508  }
509  directions.clear();
510 
511  // we do not want more orientations than nessessary
512  std::vector< WUnitSphereCoordinates > ori2( orientations.begin(), orientations.begin() + numElements );
513 
514  WMatrix< double > p = pseudoInverse( TEMat );
515 
516  return p * calcBaseMatrix( ori2, order );
517 }
518 
520 {
521  double scale = 0.0;
522  if( m_SHCoefficients.size() > 0 )
523  scale = std::sqrt( 4.0 * piDouble ) * m_SHCoefficients[0];
524  for( size_t i = 0; i < m_SHCoefficients.size(); i++ )
525  {
526  m_SHCoefficients[ i ] /= scale;
527  }
528 }
529 // NOLINT