OpenWalnut  1.2.5
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 #include "../WExportCommon.h"
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  */
41 class OWCOMMON_EXPORT WSymmetricSphericalHarmonic // NOLINT
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  */
59  virtual ~WSymmetricSphericalHarmonic();
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  */
90  WValue< double > getCoefficientsSchultz() const;
91 
92  /**
93  * Returns the coefficients for the complex base.
94  *
95  * \return coefficiend list
96  */
97  WValue< std::complex< double > > getCoefficientsComplex() const;
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  * Calculates the base matrix B like in the dissertation of Descoteaux.
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  * \return The base Matrix B
174  */
175  static WMatrix<double> calcBaseMatrix( const std::vector< WUnitSphereCoordinates >& orientations, int order );
176 
177  /**
178  * Calculates the base matrix B for the complex spherical harmonics.
179  * \param orientations The vector with the used orientation on the unit sphere (usually the gradients of the HARDI)
180  * \param order The order of the spherical harmonics intended to create
181  * \return The base Matrix B
182  */
183  static WMatrix< std::complex< double > > calcComplexBaseMatrix( std::vector< WUnitSphereCoordinates > const& orientations,
184  int order );
185 
186  /**
187  * This calcs the smoothing matrix L from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging"
188  * \param order The order of the spherical harmonic
189  * \return The smoothing matrix L
190  */
191  static WMatrix<double> calcSmoothingMatrix( size_t order );
192 
193  /**
194  * Calculates the Funk-Radon-Transformation-Matrix P from the 2007 Descoteaux Paper "Regularized, Fast, and Robust Analytical Q-Ball Imaging"
195  * \param order The order of the spherical harmonic
196  * \return The Funk-Radon-Matrix P
197  */
198  static WMatrix<double> calcFRTMatrix( size_t order );
199 
200  /**
201  * Calculates a matrix that converts spherical harmonics to symmetric tensors of equal or lower order.
202  *
203  * \param order The order of the symmetric tensor.
204  * \param orientations A vector of at least (orderTensor+1) * (orderTensor+2) / 2 orientations.
205  *
206  * \return the conversion matrix
207  */
208  static WMatrix< double > calcSHToTensorSymMatrix( std::size_t order, const std::vector< WUnitSphereCoordinates >& orientations );
209 
210  /**
211  * Normalize this SH in place.
212  */
213  void normalize();
214 
215 protected:
216 
217 private:
218  /** order of the spherical harmonic */
219  size_t m_order;
220 
221  /** coefficients of the spherical harmonic */
223 };
224 
225 #endif // WSYMMETRICSPHERICALHARMONIC_H