OpenWalnut  1.3.1
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 <vector>
29 
30 
31 #include "linearAlgebra/WLinearAlgebra.h"
32 #include "WMath.h"
33 #include "WUnitSphereCoordinates.h"
34 #include "WMatrix.h"
35 #include "WValue.h"
36 
37 /**
38  * Class for symmetric spherical harmonics
39  * The index scheme of the coefficients/basis values is like in the Descoteaux paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging".
40  */
42 {
43 // friend class WSymmetricSphericalHarmonicTest;
44 public:
45  /**
46  * Default constructor.
47  */
49 
50  /**
51  * Constructor.
52  * \param SHCoefficients the initial coefficients (stored like in the mentioned Descoteaux paper).
53  */
54  explicit WSymmetricSphericalHarmonic( const WValue<double>& SHCoefficients );
55 
56  /**
57  * Destructor.
58  */
60 
61  /**
62  * Return the value on the sphere.
63  * \param theta angle for the position on the unit sphere
64  * \param phi angle for the position on the unit sphere
65  *
66  * \return value on sphere
67  */
68  double getValue( double theta, double phi ) const;
69 
70  /**
71  * Return the value on the sphere.
72  * \param coordinates for the position on the unit sphere
73  *
74  * \return value on sphere
75  */
76  double getValue( const WUnitSphereCoordinates& coordinates ) const;
77 
78  /**
79  * Returns the used coefficients (stored like in the mentioned 2007 Descoteaux paper).
80  *
81  * \return coefficient list
82  */
83  const WValue<double>& getCoefficients() const;
84 
85  /**
86  * Returns the coefficients for Schultz' SH base.
87  *
88  * \return coefficient list
89  */
91 
92  /**
93  * Returns the coefficients for the complex base.
94  *
95  * \return coefficiend list
96  */
98 
99  /**
100  * Applies the Funk-Radon-Transformation. This is faster than matrix multiplication.
101  * ( O(n) instead of O(n²) )
102  *
103  * \param frtMat the frt matrix as calculated by calcFRTMatrix()
104  */
105  void applyFunkRadonTransformation( WMatrix< double > const& frtMat );
106 
107  /**
108  * Return the order of the spherical harmonic.
109  *
110  * \return order of SH
111  */
112  size_t getOrder() const;
113 
114  /**
115  * Calculate the generalized fractional anisotropy for this ODF.
116  *
117  * See: David S. Tuch, "Q-Ball Imaging", Magn. Reson. Med. 52, 2004, 1358-1372
118  *
119  * \note this only makes sense if this is an ODF, meaning funk-radon-transform was applied etc.
120  *
121  * \param orientations A vector of unit sphere coordinates.
122  *
123  * \return The generalized fractional anisotropy.
124  */
125  double calcGFA( std::vector< WUnitSphereCoordinates > const& orientations ) const;
126 
127  /**
128  * Calculate the generalized fractional anisotropy for this ODF. This version of
129  * the function uses precomputed base functions (because calculating the base function values
130  * is rather expensive). Use this version if you want to compute the GFA for multiple ODFs
131  * with the same base functions. The base function Matrix can be computed using \see calcBMatrix().
132  *
133  * See: David S. Tuch, "Q-Ball Imaging", Magn. Reson. Med. 52, 2004, 1358-1372
134  *
135  * \note this only makes sense if this is an ODF, meaning funk-radon-transform was applied etc.
136  *
137  * \param B The matrix of SH base functions.
138  *
139  * \return The generalized fractional anisotropy.
140  */
141  double calcGFA( WMatrix< double > const& B ) const;
142 
143  /**
144  * This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper. The orientations are given as WVector3d.
145  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
146  * \param order The order of the spherical harmonics intended to create
147  * \param lambda Regularization parameter for smoothing matrix
148  * \param withFRT include the Funk-Radon-Transformation?
149  * \return Transformation matrix
150  */
151  static WMatrix<double> getSHFittingMatrix( const std::vector< WVector3d >& orientations,
152  int order,
153  double lambda,
154  bool withFRT );
155 
156  /**
157  * This calculates the transformation/fitting matrix T like in the 2007 Descoteaux paper. The orientations are given as WUnitSphereCoordinates .
158  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
159  * \param order The order of the spherical harmonics intended to create
160  * \param lambda Regularization parameter for smoothing matrix
161  * \param withFRT include the Funk-Radon-Transformation?
162  * \return Transformation matrix
163  */
164  static WMatrix<double> getSHFittingMatrix( const std::vector< WUnitSphereCoordinates >& orientations,
165  int order,
166  double lambda,
167  bool withFRT );
168 
169  /**
170  * This calculates the transformation/fitting matrix T like in the 2010 Aganj paper. The orientations are given as WUnitSphereCoordinates .
171  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
172  * \param order The order of the spherical harmonics intended to create
173  * \param lambda Regularization parameter for smoothing matrix
174  * \return Transformation matrix
175  */
176  static WMatrix<double> getSHFittingMatrixForConstantSolidAngle( const std::vector< WVector3d >& orientations,
177  int order,
178  double lambda );
179 
180  /**
181  * This calculates the transformation/fitting matrix T like in the 2010 Aganj paper. The orientations are given as WUnitSphereCoordinates .
182  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
183  * \param order The order of the spherical harmonics intended to create
184  * \param lambda Regularization parameter for smoothing matrix
185  * \return Transformation matrix
186  */
187  static WMatrix<double> getSHFittingMatrixForConstantSolidAngle( const std::vector< WUnitSphereCoordinates >& orientations,
188  int order,
189  double lambda );
190 
191  /**
192  * Calculates the base matrix B like in the dissertation of Descoteaux.
193  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
194  * \param order The order of the spherical harmonics intended to create
195  * \return The base Matrix B
196  */
197  static WMatrix<double> calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations, int order );
198 
199  /**
200  * Calculates the base matrix B for the complex spherical harmonics.
201  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
202  * \param order The order of the spherical harmonics intended to create
203  * \return The base Matrix B
204  */
205  static WMatrix< std::complex< double > > calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations,
206  int order );
207  /**
208  * Calc eigenvalues for SH elements.
209  * \param order The order of the spherical harmonic
210  * \return The eigenvalues of the coefficients
211  */
212  static WValue<double> calcEigenvalues( size_t order );
213 
214  /**
215  * Calc matrix with the eigenvalues of the SH elements on its diagonal.
216  * \param order The order of the spherical harmonic
217  * \return The matrix with the eigenvalues of the coefficients
218  */
219  static WMatrix<double> calcMatrixWithEigenvalues( size_t order );
220 
221  /**
222  * This calcs the smoothing matrix L from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging"
223  * \param order The order of the spherical harmonic
224  * \return The smoothing matrix L
225  */
226  static WMatrix<double> calcSmoothingMatrix( size_t order );
227 
228  /**
229  * Calculates the Funk-Radon-Transformation-Matrix P from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging"
230  * \param order The order of the spherical harmonic
231  * \return The Funk-Radon-Matrix P
232  */
233  static WMatrix<double> calcFRTMatrix( size_t order );
234 
235  /**
236  * Calculates a matrix that converts spherical harmonics to symmetric tensors of equal or lower order.
237  *
238  * \param order The order of the symmetric tensor.
239  * \param orientations A vector of at least (orderTensor+1) * (orderTensor+2) / 2 orientations.
240  *
241  * \return the conversion matrix
242  */
243  static WMatrix< double > calcSHToTensorSymMatrix( std::size_t order, const std::vector< WUnitSphereCoordinates >& orientations );
244 
245  /**
246  * Normalize this SH in place.
247  */
248  void normalize();
249 
250 protected:
251 private:
252  /** order of the spherical harmonic */
253  size_t m_order;
254 
255  /** coefficients of the spherical harmonic */
257 };
258 
259 #endif // WSYMMETRICSPHERICALHARMONIC_H