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