OpenWalnut  1.4.0
WSymmetricSphericalHarmonic.h
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 #ifndef WSYMMETRICSPHERICALHARMONIC_H
26 #define WSYMMETRICSPHERICALHARMONIC_H
27 
28 #include <stdint.h>
29 
30 #include <cmath>
31 #include <vector>
32 
33 #ifndef Q_MOC_RUN
34 #include <boost/math/special_functions/spherical_harmonic.hpp>
35 #endif
36 
37 #include "core/common/WLogger.h"
38 #include "core/common/math/WGeometryFunctions.h"
39 #include "../exceptions/WPreconditionNotMet.h"
40 #include "linearAlgebra/WVectorFixed.h"
41 #include "WLinearAlgebraFunctions.h"
42 #include "WMath.h"
43 #include "WMatrix.h"
44 #include "WTensorSym.h"
45 #include "WUnitSphereCoordinates.h"
46 #include "WValue.h"
47 
48 
49 /**
50  * Class for symmetric spherical harmonics
51  * The index scheme of the coefficients/basis values is like in the Descoteaux paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".
52  */
53 template< typename T > class WSymmetricSphericalHarmonic // NOLINT
54 {
55 // friend class WSymmetricSphericalHarmonicTest;
56 public:
57  /**
58  * Default constructor.
59  */
61 
62  /**
63  * Constructor.
64  * \param SHCoefficients the initial coefficients (stored like in the mentioned Descoteaux paper).
65  */
66  explicit WSymmetricSphericalHarmonic( const WValue< T >& SHCoefficients );
67 
68  /**
69  * Destructor.
70  */
72 
73  /**
74  * Return the value on the sphere.
75  * \param theta angle for the position on the unit sphere
76  * \param phi angle for the position on the unit sphere
77  *
78  * \return value on sphere
79  */
80  T getValue( T theta, T phi ) const;
81 
82  /**
83  * Return the value on the sphere.
84  * \param coordinates for the position on the unit sphere
85  *
86  * \return value on sphere
87  */
88  T getValue( const WUnitSphereCoordinates< T >& coordinates ) const;
89 
90  /**
91  * Returns the used coefficients (stored like in the mentioned 2007 Descoteaux paper).
92  *
93  * \return coefficient list
94  */
95  const WValue< T >& getCoefficients() const;
96 
97  /**
98  * Returns the coefficients for Schultz' SH base.
99  *
100  * \return coefficient list
101  */
103 
104  /**
105  * Returns the coefficients for the complex base.
106  *
107  * \return coefficiend list
108  */
110 
111  /**
112  * Applies the Funk-Radon-Transformation. This is faster than matrix multiplication.
113  * ( O(n) instead of O(n²) )
114  *
115  * \param frtMat the frt matrix as calculated by calcFRTMatrix()
116  */
117  void applyFunkRadonTransformation( WMatrix< T > const& frtMat );
118 
119  /**
120  * Return the order of the spherical harmonic.
121  *
122  * \return order of SH
123  */
124  size_t getOrder() const;
125 
126  /**
127  * Calculate the generalized fractional anisotropy for this ODF.
128  *
129  * See: David S. Tuch, "Q-Ball Imaging", Magn. Reson. Med. 52, 2004, 1358-1372
130  *
131  * \note this only makes sense if this is an ODF, meaning funk-radon-transform was applied etc.
132  *
133  * \param orientations A vector of unit sphere coordinates.
134  *
135  * \return The generalized fractional anisotropy.
136  */
137  T calcGFA( std::vector< WUnitSphereCoordinates< T > > const& orientations ) const;
138 
139  /**
140  * Calculate the generalized fractional anisotropy for this ODF. This version of
141  * the function uses precomputed base functions (because calculating the base function values
142  * is rather expensive). Use this version if you want to compute the GFA for multiple ODFs
143  * with the same base functions. The base function Matrix can be computed using \see calcBMatrix().
144  *
145  * See: David S. Tuch, "Q-Ball Imaging", Magn. Reson. Med. 52, 2004, 1358-1372
146  *
147  * \note this only makes sense if this is an ODF, meaning funk-radon-transform was applied etc.
148  *
149  * \param B The matrix of SH base functions.
150  *
151  * \return The generalized fractional anisotropy.
152  */
153  T calcGFA( WMatrix< T > const& B ) const;
154 
155  /**
156  * This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper. The orientations are given as WMatrixFixed< T, 3, 1 >.
157  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
158  * \param order The order of the spherical harmonics intended to create
159  * \param lambda Regularization parameter for smoothing matrix
160  * \param withFRT include the Funk-Radon-Transformation?
161  * \return Transformation matrix
162  */
163  static WMatrix< T > getSHFittingMatrix( const std::vector< WMatrixFixed< T, 3, 1 > >& orientations,
164  int order,
165  T lambda,
166  bool withFRT );
167 
168  /**
169  * This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper. The orientations are given as WUnitSphereCoordinates< T >.
170  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
171  * \param order The order of the spherical harmonics intended to create
172  * \param lambda Regularization parameter for smoothing matrix
173  * \param withFRT include the Funk-Radon-Transformation?
174  * \return Transformation matrix
175  */
176  static WMatrix< T > getSHFittingMatrix( const std::vector< WUnitSphereCoordinates< T > >& orientations,
177  int order,
178  T lambda,
179  bool withFRT );
180 
181  /**
182  * This calculates the transformation/fitting matrix T like in the 2010 Aganj paper. The orientations are given as WUnitSphereCoordinates< T >.
183  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
184  * \param order The order of the spherical harmonics intended to create
185  * \param lambda Regularization parameter for smoothing matrix
186  * \return Transformation matrix
187  */
188  static WMatrix< T > getSHFittingMatrixForConstantSolidAngle( const std::vector< WMatrixFixed< T, 3, 1 > >& orientations,
189  int order,
190  T lambda );
191 
192  /**
193  * This calculates the transformation/fitting matrix T like in the 2010 Aganj paper. The orientations are given as WUnitSphereCoordinates< T >.
194  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
195  * \param order The order of the spherical harmonics intended to create
196  * \param lambda Regularization parameter for smoothing matrix
197  * \return Transformation matrix
198  */
199  static WMatrix< T > getSHFittingMatrixForConstantSolidAngle( const std::vector< WUnitSphereCoordinates< T > >& orientations,
200  int order,
201  T lambda );
202 
203  /**
204  * Calculates the base matrix B like in the dissertation of Descoteaux.
205  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
206  * \param order The order of the spherical harmonics intended to create
207  * \return The base Matrix B
208  */
209  static WMatrix< T > calcBaseMatrix( const std::vector< WUnitSphereCoordinates< T > >& orientations, int order );
210 
211  /**
212  * Calculates the base matrix B for the complex spherical harmonics.
213  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
214  * \param order The order of the spherical harmonics intended to create
215  * \return The base Matrix B
216  */
217  static WMatrix< std::complex< T > > calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates< T > > const& orientations,
218  int order );
219  /**
220  * Calc eigenvalues for SH elements.
221  * \param order The order of the spherical harmonic
222  * \return The eigenvalues of the coefficients
223  */
224  static WValue< T > calcEigenvalues( size_t order );
225 
226  /**
227  * Calc matrix with the eigenvalues of the SH elements on its diagonal.
228  * \param order The order of the spherical harmonic
229  * \return The matrix with the eigenvalues of the coefficients
230  */
231  static WMatrix< T > calcMatrixWithEigenvalues( size_t order );
232 
233  /**
234  * This calcs the smoothing matrix L from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging"
235  * \param order The order of the spherical harmonic
236  * \return The smoothing matrix L
237  */
238  static WMatrix< T > calcSmoothingMatrix( size_t order );
239 
240  /**
241  * Calculates the Funk-Radon-Transformation-Matrix P from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging"
242  * \param order The order of the spherical harmonic
243  * \return The Funk-Radon-Matrix P
244  */
245  static WMatrix< T > calcFRTMatrix( size_t order );
246 
247  /**
248  * Calculates a matrix that converts spherical harmonics to symmetric tensors of equal or lower order.
249  *
250  * \param order The order of the symmetric tensor.
251  *
252  * \return the conversion matrix
253  */
254  static WMatrix< T > calcSHToTensorSymMatrix( std::size_t order );
255 
256  /**
257  * Calculates a matrix that converts spherical harmonics to symmetric tensors of equal or lower order.
258  *
259  * \param order The order of the symmetric tensor.
260  * \param orientations A vector of at least (orderTensor+1) * (orderTensor+2) / 2 orientations.
261  *
262  * \return the conversion matrix
263  */
264  static WMatrix< T > calcSHToTensorSymMatrix( std::size_t order, const std::vector< WUnitSphereCoordinates< T > >& orientations );
265 
266  /**
267  * Normalize this SH in place.
268  */
269  void normalize();
270 
271 protected:
272 private:
273  /** order of the spherical harmonic */
274  size_t m_order;
275 
276  /** coefficients of the spherical harmonic */
278 };
279 
280 template< typename T >
282  m_order( 0 ),
283  m_SHCoefficients( 0 )
284 {
285 }
286 
287 template< typename T >
289  m_SHCoefficients( SHCoefficients )
290 {
291  // calculate order from SHCoefficients.size:
292  // - this is done by reversing the R=(l+1)*(l+2)/2 formula from the Descoteaux paper
293  T q = std::sqrt( 0.25 + 2.0 * static_cast< T >( m_SHCoefficients.size() ) ) - 1.5;
294  m_order = static_cast<size_t>( q );
295 }
296 
297 template< typename T >
299 {
300 }
301 
302 template< typename T >
304 {
305  T result = 0.0;
306  int j = 0;
307  const T rootOf2 = std::sqrt( 2.0 );
308  for( int k = 0; k <= static_cast<int>( m_order ); k+=2 )
309  {
310  // m = 1 .. k
311  for( int m = 1; m <= k; m++ )
312  {
313  j = ( k*k+k+2 ) / 2 + m;
314  result += m_SHCoefficients[ j-1 ] * rootOf2 *
315  static_cast< T >( std::pow( -1.0, m+1 ) ) * boost::math::spherical_harmonic_i( k, m, theta, phi );
316  }
317  // m = 0
318  result += m_SHCoefficients[ ( k*k+k+2 ) / 2 - 1 ] * boost::math::spherical_harmonic_r( k, 0, theta, phi );
319  // m = -k .. -1
320  for( int m = -k; m < 0; m++ )
321  {
322  j = ( k*k+k+2 ) / 2 + m;
323  result += m_SHCoefficients[ j-1 ] * rootOf2 * boost::math::spherical_harmonic_r( k, -m, theta, phi );
324  }
325  }
326  return result;
327 }
328 
329 template< typename T >
331 {
332  return getValue( coordinates.getTheta(), coordinates.getPhi() );
333 }
334 
335 template< typename T >
337 {
338  return m_SHCoefficients;
339 }
340 
341 template< typename T >
343 {
344  WValue< T > res( m_SHCoefficients.size() );
345  size_t k = 0;
346  for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
347  {
348  for( int m = -l; m <= l; ++m )
349  {
350  res[ k ] = m_SHCoefficients[ k ];
351  if( m < 0 && ( ( -m ) % 2 == 1 ) )
352  {
353  res[ k ] *= -1.0;
354  }
355  else if( m > 0 && ( m % 2 == 0 ) )
356  {
357  res[ k ] *= -1.0;
358  }
359  ++k;
360  }
361  }
362  return res;
363 }
364 
365 template< typename T >
367 {
368  WValue< std::complex< T > > res( m_SHCoefficients.size() );
369  size_t k = 0;
370  T r = 1.0 / sqrt( 2.0 );
371  std::complex< T > i( 0.0, -1.0 );
372  for( int l = 0; l <= static_cast< int >( m_order ); l += 2 )
373  {
374  for( int m = -l; m <= l; ++m )
375  {
376  if( m == 0 )
377  {
378  res[ k ] = m_SHCoefficients[ k ];
379  }
380  else if( m < 0 )
381  {
382  res[ k ] += i * r * m_SHCoefficients[ k - 2 * m ];
383  res[ k ] += ( -m % 2 == 1 ? -r : r ) * m_SHCoefficients[ k ];
384  }
385  else if( m > 0 )
386  {
387  res[ k ] += i * ( -m % 2 == 0 ? -r : r ) * m_SHCoefficients[ k ];
388  res[ k ] += r * m_SHCoefficients[ k - 2 * m ];
389  }
390  ++k;
391  }
392  }
393  return res;
394 }
395 
396 template< typename T >
397 T WSymmetricSphericalHarmonic< T >::calcGFA( std::vector< WUnitSphereCoordinates< T > > const& orientations ) const
398 {
399  T n = static_cast< T >( orientations.size() );
400  T d = 0.0;
401  T gfa = 0.0;
402  T mean = 0.0;
403  std::vector< T > v( orientations.size() );
404 
405  for( std::size_t i = 0; i < orientations.size(); ++i )
406  {
407  v[ i ] = getValue( orientations[ i ] );
408  mean += v[ i ];
409  }
410  mean /= n;
411 
412  for( std::size_t i = 0; i < orientations.size(); ++i )
413  {
414  T f = v[ i ];
415  // we use gfa as a temporary here
416  gfa += f * f;
417  f -= mean;
418  d += f * f;
419  }
420 
421  if( gfa == 0.0 || n <= 1.0 )
422  {
423  return 0.0;
424  }
425  // this is the real gfa
426  gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
427 
428  WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
429  if( gfa < 0.0 )
430  {
431  return 0.0;
432  }
433  else if( gfa > 1.0 )
434  {
435  return 1.0;
436  }
437  return gfa;
438 }
439 
440 template< typename T >
442 {
443  // sh coeffs
444  WMatrix< T > s( B.getNbCols(), 1 );
445  WAssert( B.getNbCols() == m_SHCoefficients.size(), "" );
446  for( std::size_t k = 0; k < B.getNbCols(); ++k )
447  {
448  s( k, 0 ) = m_SHCoefficients[ k ];
449  }
450  s = B * s;
451  WAssert( s.getNbRows() == B.getNbRows(), "" );
452  WAssert( s.getNbCols() == 1u, "" );
453 
454  T n = static_cast< T >( s.getNbRows() );
455  T d = 0.0;
456  T gfa = 0.0;
457  T mean = 0.0;
458  for( std::size_t i = 0; i < s.getNbRows(); ++i )
459  {
460  mean += s( i, 0 );
461  }
462  mean /= n;
463 
464  for( std::size_t i = 0; i < s.getNbRows(); ++i )
465  {
466  T f = s( i, 0 );
467  // we use gfa as a temporary here
468  gfa += f * f;
469  f -= mean;
470  d += f * f;
471  }
472 
473  if( gfa == 0.0 || n <= 1.0 )
474  {
475  return 0.0;
476  }
477  // this is the real gfa
478  gfa = sqrt( ( n * d ) / ( ( n - 1.0 ) * gfa ) );
479 
480  WAssert( gfa >= -0.01 && gfa <= 1.01, "" );
481  if( gfa < 0.0 )
482  {
483  return 0.0;
484  }
485  else if( gfa > 1.0 )
486  {
487  return 1.0;
488  }
489  return gfa;
490 }
491 
492 template< typename T >
494 {
495  WAssert( frtMat.getNbCols() == m_SHCoefficients.size(), "" );
496  WAssert( frtMat.getNbRows() == m_SHCoefficients.size(), "" );
497  // Funk-Radon-Transformation as in Descoteaux's thesis
498  for( size_t j = 0; j < m_SHCoefficients.size(); j++ )
499  {
500  m_SHCoefficients[ j ] *= frtMat( j, j );
501  }
502 }
503 
504 template< typename T >
506 {
507  return m_order;
508 }
509 
510 template< typename T >
512  int order,
513  T lambda,
514  bool withFRT )
515 {
516  // convert euclid 3d orientations/gradients to sphere coordinates
517  std::vector< WUnitSphereCoordinates< T > > sphereCoordinates;
518  for( typename std::vector< WMatrixFixed< T, 3, 1 > >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
519  {
520  sphereCoordinates.push_back( WUnitSphereCoordinates< T >( *it ) );
521  }
522  return WSymmetricSphericalHarmonic< T >::getSHFittingMatrix( sphereCoordinates, order, lambda, withFRT );
523 }
524 
525 template< typename T >
527  int order,
528  T lambda,
529  bool withFRT )
530 {
532  WMatrix< T > Bt( B.transposed() );
533  WMatrix< T > result( Bt*B );
534  if( lambda != 0.0 )
535  {
537  result += lambda*L;
538  }
539 
540  result = pseudoInverse( result )*Bt;
541  if( withFRT )
542  {
544  return ( P * result );
545  }
546  return result;
547 }
548 
549 template< typename T >
551  int order,
552  T lambda )
553 {
554  // convert euclid 3d orientations/gradients to sphere coordinates
555  std::vector< WUnitSphereCoordinates< T > > sphereCoordinates;
556  for( typename std::vector< WMatrixFixed< T, 3, 1 > >::const_iterator it = orientations.begin(); it != orientations.end(); it++ )
557  {
558  sphereCoordinates.push_back( WUnitSphereCoordinates< T >( *it ) );
559  }
560  return WSymmetricSphericalHarmonic< T >::getSHFittingMatrixForConstantSolidAngle( sphereCoordinates, order, lambda );
561 }
562 
563 template< typename T >
565  const std::vector< WUnitSphereCoordinates< T > >& orientations,
566  int order,
567  T lambda )
568 {
570  WMatrix< T > Bt( B.transposed() );
571  WMatrix< T > result( Bt*B );
572 
573  // regularisation
574  if( lambda != 0.0 )
575  {
577  result += lambda*L;
578  }
579 
580  result = pseudoInverse( result )*Bt;
581 
582  // multiply with eigenvalues of coefficients / Laplace-Beltrami operator
584  wlog::debug( "" ) << "LB: " << LB;
585  result = LB*result;
586  wlog::debug( "" ) << "LB*result: " << result;
587 
588  // apply FRT
590  result = P * result;
591  wlog::debug( "" ) << "P: " << P;
592  wlog::debug( "" ) << "P*result: " << result;
593 
594  // correction factor
595  T correctionFactor = 1.0 / ( 16.0 * std::pow( static_cast< T >( piDouble ), 2 ) );
596  result *= correctionFactor;
597 
598  return result;
599 }
600 
601 
602 template< typename T >
604  int order )
605 {
606  WMatrix< T > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
607 
608  // calc B Matrix like in the 2007 Descoteaux paper ("Regularized, Fast, and Robust Analytical Q-Ball Imaging")
609  int j = 0;
610  const T rootOf2 = std::sqrt( 2.0 );
611  for(uint32_t row = 0; row < orientations.size(); row++ )
612  {
613  const T theta = orientations[row].getTheta();
614  const T phi = orientations[row].getPhi();
615  for( int k = 0; k <= order; k+=2 )
616  {
617  // m = 1 .. k
618  for( int m = 1; m <= k; m++ )
619  {
620  j = ( k * k + k + 2 ) / 2 + m;
621  B( row, j-1 ) = rootOf2 * static_cast< T >( std::pow( static_cast< T >( -1 ), m + 1 ) )
622  * boost::math::spherical_harmonic_i( k, m, theta, phi );
623  }
624  // m = 0
625  B( row, ( k * k + k + 2 ) / 2 - 1 ) = boost::math::spherical_harmonic_r( k, 0, theta, phi );
626  // m = -k .. -1
627  for( int m = -k; m < 0; m++ )
628  {
629  j = ( k * k + k + 2 ) / 2 + m;
630  B( row, j-1 ) = rootOf2*boost::math::spherical_harmonic_r( k, -m, theta, phi );
631  }
632  }
633  }
634  return B;
635 }
636 
637 template< typename T >
640 {
641  WMatrix< std::complex< T > > B( orientations.size(), ( ( order + 1 ) * ( order + 2 ) ) / 2 );
642 
643  for( std::size_t row = 0; row < orientations.size(); row++ )
644  {
645  const T theta = orientations[ row ].getTheta();
646  const T phi = orientations[ row ].getPhi();
647 
648  int j = 0;
649  for( int k = 0; k <= order; k += 2 )
650  {
651  for( int m = -k; m < 0; m++ )
652  {
653  B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
654  ++j;
655  }
656  B( row, j ) = boost::math::spherical_harmonic( k, 0, theta, phi );
657  ++j;
658  for( int m = 1; m <= k; m++ )
659  {
660  B( row, j ) = boost::math::spherical_harmonic( k, m, theta, phi );
661  ++j;
662  }
663  }
664  }
665  return B;
666 }
667 
668 template< typename T >
670 {
671  size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
672  std::size_t i = 0;
673  WValue< T > result( R );
674  for( size_t k = 0; k <= order; k += 2 )
675  {
676  for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
677  {
678  result[ i ] = -static_cast< T > ( k * ( k + 1 ) );
679  ++i;
680  }
681  }
682  return result;
683 }
684 
685 template< typename T >
687 {
689  WMatrix< T > L( eigenvalues.size(), eigenvalues.size() );
690  for( std::size_t i = 0; i < eigenvalues.size(); ++i )
691  {
692  L( i, i ) = eigenvalues[ i ];
693  }
694  return L;
695 }
696 
697 template< typename T >
699 {
701  WMatrix< T > L( eigenvalues.size(), eigenvalues.size() );
702  for( std::size_t i = 0; i < eigenvalues.size(); ++i )
703  {
704  L( i, i ) = std::pow( eigenvalues[ i ], 2 );
705  }
706  return L;
707 }
708 
709 template< typename T >
711 {
712  size_t R = ( ( order + 1 ) * ( order + 2 ) ) / 2;
713  std::size_t i = 0;
714  WMatrix< T > result( R, R );
715  for( size_t k = 0; k <= order; k += 2 )
716  {
717  T h = 2.0 * static_cast< T >( piDouble ) * static_cast< T >( std::pow( static_cast< T >( -1 ), static_cast< T >( k / 2 ) ) ) *
718  static_cast< T >( oddFactorial( ( k <= 1 ) ? 1 : k - 1 ) ) / static_cast< T >( evenFactorial( k ) );
719  for( int m = -static_cast< int >( k ); m <= static_cast< int >( k ); ++m )
720  {
721  result( i, i ) = h;
722  ++i;
723  }
724  }
725  return result;
726 }
727 
728 template< typename T >
730 {
731  std::vector< WVector3d > vertices;
732  std::vector< unsigned int > triangles;
733  // calc necessary tesselation level
734  int level = static_cast< int >( log( static_cast< float >( ( order + 1 ) * ( order + 2 ) ) ) / 2.0f + 1.0f );
735  tesselateIcosahedron( &vertices, &triangles, level );
736  std::vector< WUnitSphereCoordinates< T > > orientations;
737  for( typename std::vector< WVector3d >::const_iterator cit = vertices.begin(); cit != vertices.end(); cit++ )
738  {
739  // avoid linear dependent orientations
740  if( ( *cit )[ 0 ] >= 0.0001 )
741  {
742  orientations.push_back( WUnitSphereCoordinates< T >( WMatrixFixed< T, 3, 1 >( *cit ) ) );
743  }
744  }
745  return WSymmetricSphericalHarmonic< T >::calcSHToTensorSymMatrix( order, orientations );
746 }
747 
748 template< typename T >
750  const std::vector< WUnitSphereCoordinates< T > >& orientations )
751 {
752  std::size_t numElements = ( order + 1 ) * ( order + 2 ) / 2;
753  WPrecondEquals( order % 2, 0u );
754  WPrecondLess( numElements, orientations.size() + 1 );
755 
756  // store first numElements orientations as 3d-vectors
757  std::vector< WMatrixFixed< T, 3, 1 > > directions( numElements );
758  for( std::size_t i = 0; i < numElements; ++i )
759  {
760  directions[ i ] = orientations[ i ].getEuclidean();
761  }
762 
763  // initialize an index array
764  std::vector< std::size_t > indices( order, 0 );
765 
766  // calc tensor evaluation matrix
767  WMatrix< T > TEMat( numElements, numElements );
768  for( std::size_t j = 0; j < numElements; ++j ) // j is the 'permutation index'
769  {
770  // stores how often each value is represented in the index array
771  std::size_t amount[ 3 ] = { 0, 0, 0 };
772  for( std::size_t k = 0; k < order; ++k )
773  {
774  ++amount[ indices[ k ] ];
775  }
776 
777  // from WTensorSym.h
778  std::size_t multiplicity = calcSupersymmetricTensorMultiplicity( order, amount[ 0 ], amount[ 1 ], amount[ 2 ] );
779  for( std::size_t i = 0; i < numElements; ++i ) // i is the 'direction index'
780  {
781  TEMat( i, j ) = multiplicity;
782  for( std::size_t k = 0; k < order; ++k )
783  {
784  TEMat( i, j ) *= directions[ i ][ indices[ k ] ];
785  }
786  }
787 
788  // from TensorBase.h
789  positionIterateSortedOneStep( order, 3, indices );
790  }
791  directions.clear();
792 
793  // we do not want more orientations than nessessary
794  std::vector< WUnitSphereCoordinates< T > > ori2( orientations.begin(), orientations.begin() + numElements );
795 
796  WMatrix< T > p = pseudoInverse( TEMat );
797 
798  WAssert( ( TEMat*p ).isIdentity( 1e-3 ), "Test of inverse matrix failed. Are the given orientations linear independent?" );
799 
800  return p * calcBaseMatrix( ori2, order );
801 }
802 
803 template< typename T >
805 {
806  T scale = 0.0;
807  if( m_SHCoefficients.size() > 0 )
808  {
809  scale = std::sqrt( 4.0 * static_cast< T >( piDouble ) ) * m_SHCoefficients[0];
810  }
811  for( size_t i = 0; i < m_SHCoefficients.size(); i++ )
812  {
813  m_SHCoefficients[ i ] /= scale;
814  }
815 }
816 
817 #endif // WSYMMETRICSPHERICALHARMONIC_H
WSymmetricSphericalHarmonic()
Default constructor.
static WMatrix< T > calcFRTMatrix(size_t order)
Calculates the Funk-Radon-Transformation-Matrix P from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".
WValue< T > m_SHCoefficients
coefficients of the spherical harmonic
static WMatrix< T > getSHFittingMatrixForConstantSolidAngle(const std::vector< WMatrixFixed< T, 3, 1 > > &orientations, int order, T lambda)
This calculates the transformation/fitting matrix T like in the 2010 Aganj paper. ...
static WMatrix< T > calcSmoothingMatrix(size_t order)
This calcs the smoothing matrix L from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".
T getValue(T theta, T phi) const
Return the value on the sphere.
void normalize()
Normalize this SH in place.
Base class for all higher level values like tensors, vectors, matrices and so on. ...
Definition: WValue.h:40
size_t getNbCols() const
Get number of columns.
Definition: WMatrix.h:383
virtual ~WSymmetricSphericalHarmonic()
Destructor.
This class stores coordinates on the unit sphere.
T getPhi() const
Return the phi angle.
Matrix template class with variable number of rows and columns.
const WValue< T > & getCoefficients() const
Returns the used coefficients (stored like in the mentioned 2007 Descoteaux paper).
WValue< std::complex< T > > getCoefficientsComplex() const
Returns the coefficients for the complex base.
size_t getOrder() const
Return the order of the spherical harmonic.
WMatrix transposed() const
Returns the transposed matrix.
Definition: WMatrix.h:444
static WValue< T > calcEigenvalues(size_t order)
Calc eigenvalues for SH elements.
Class for symmetric spherical harmonics The index scheme of the coefficients/basis values is like in ...
T calcGFA(std::vector< WUnitSphereCoordinates< T > > const &orientations) const
Calculate the generalized fractional anisotropy for this ODF.
static WMatrix< T > calcBaseMatrix(const std::vector< WUnitSphereCoordinates< T > > &orientations, int order)
Calculates the base matrix B like in the dissertation of Descoteaux.
static WMatrix< T > calcSHToTensorSymMatrix(std::size_t order)
Calculates a matrix that converts spherical harmonics to symmetric tensors of equal or lower order...
size_t m_order
order of the spherical harmonic
size_t size() const
Get number of components the value consists of.
Definition: WValue.h:116
static WMatrix< T > calcMatrixWithEigenvalues(size_t order)
Calc matrix with the eigenvalues of the SH elements on its diagonal.
size_t getNbRows() const
Get number of rows.
Definition: WMatrix.h:375
T getTheta() const
Return the theta angle.
static WMatrix< std::complex< T > > calcComplexBaseMatrix(std::vector< WUnitSphereCoordinates< T > > const &orientations, int order)
Calculates the base matrix B for the complex spherical harmonics.
WValue< T > getCoefficientsSchultz() const
Returns the coefficients for Schultz' SH base.
WStreamedLogger debug(const std::string &source)
Logging a debug message.
Definition: WLogger.h:335
void applyFunkRadonTransformation(WMatrix< T > const &frtMat)
Applies the Funk-Radon-Transformation.
static WMatrix< T > getSHFittingMatrix(const std::vector< WMatrixFixed< T, 3, 1 > > &orientations, int order, T lambda, bool withFRT)
This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper.