DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkSphericalHarmonicsGenerator.hxx
Go to the documentation of this file.
1 
17 #ifndef __itkSphericalHarmonicsGenerator_hxx
18 #define __itkSphericalHarmonicsGenerator_hxx
19 
20 #include "itkSphericalHarmonicsGenerator.h"
21 #include <complex>
22 #include <gsl/gsl_sf_legendre.h>
23 #include "itkImage.h"
24 #include <gsl/gsl_sf_coupling.h>
25 
26 // #include "utl.h"
27 #include "DMRITOOLConfigure.h"
28 #include "utlCore.h"
29 #include "utlDMRI.h"
30 #include "utlITK.h"
31 #include "utlDMRIStoredTables.h"
32 
33 namespace utl
34 {
35 double Gamma(const double);
36 }
37 
38 namespace itk
39 {
40 
41 template<class PreciseType>
42 SphericalHarmonicsGenerator<PreciseType>
43 ::SphericalHarmonicsGenerator()
44 {
45 }
46 
47 template <class PreciseType>
48 SphericalHarmonicsGenerator<PreciseType>
49 ::~SphericalHarmonicsGenerator()
50 {
51 }
52 
53 template < class PreciseType >
54 void
55 SphericalHarmonicsGenerator<PreciseType>
56 ::PrintSelf ( std::ostream& os, Indent indent ) const
57 {
58  Superclass::PrintSelf(os, indent);
59 }
60 
61 template<class PreciseType>
62 std::complex<PreciseType>
63 SphericalHarmonicsGenerator<PreciseType>
64 ::ComplexSH(const int l, const int m, const PreciseType theta, const PreciseType phi)
65 {
66  int absm = (m<0 ? -m : m);
67  PreciseType sign = utl::IsEven(absm)? 1.0 : -1.0;
68 
69  std::complex<PreciseType> retval(0.0,(PreciseType)(absm*phi));
70  retval = std::exp(retval) * gsl_sf_legendre_sphPlm(l,absm,cos(theta));
71 
72  if (m<0)
73  retval = sign * std::conj(retval);
74 
75  return retval;
76 }
77 
78 
79 
80 template < class PreciseType >
81 PreciseType
82 SphericalHarmonicsGenerator<PreciseType>
83 ::RealSH (const int l, const int m, const PreciseType theta, const PreciseType phi )
84 {
85  std::complex<PreciseType> cplx;
86  if(m > 0)
87  {
88  cplx = ComplexSH(l, m, theta, phi);
89  return utl::SQRT2*std::imag(cplx);
90  }
91  else if( m == 0 )
92  {
93  cplx = ComplexSH(l, m, theta, phi);
94  return std::real(cplx);
95  }
96  else
97  {
98  // NOTE: use m%2==1 is wrong, because it can be -1.
99  PreciseType sign = utl::IsEven(m)? 1.0 : -1.0;
100  // utlPrintVar5(true, l, m, utl::IsEven(m), m%2, -m%2);
101  cplx = sign * ComplexSH(l, m, theta, phi);
102  return utl::SQRT2*std::real(cplx);
103  }
104 }
105 
106 template < class PreciseType >
107 PreciseType
108 SphericalHarmonicsGenerator<PreciseType>
109 ::ComplexTripleIntegration (const int l1, const int m1, const int l2, const int m2, const int l3, const int m3)
110 {
111  utlException(l1<std::abs(m1) || l2<std::abs(m2) || l3<std::abs(m3), "wrong input! m should be less than l");
112  PreciseType result = 0.0;
113  if (m1+m2+m3==0 && l1+l2>=l3 && l1+l3>=l2 && l2+l3>=l1)
114  {
115  result = std::sqrt( (2*l1+1)*(2*l2+1)*(2*l3+1)/(4*M_PI) ) * gsl_sf_coupling_3j(2*l1,2*l2,2*l3,2*m1,2*m2,2*m3) * gsl_sf_coupling_3j(2*l1,2*l2,2*l3,0,0,0);
116  }
117  else
118  result = 0;
119  return result;
120 }
121 
122 template < class PreciseType >
123 PreciseType
124 SphericalHarmonicsGenerator<PreciseType>
125 ::RealTripleIntegration (const int l1, const int m1, const int l2, const int m2, const int l3, const int m3, const bool is_precalculated)
126 {
127  utlException(l1<std::abs(m1) || l2<std::abs(m2) || l3<std::abs(m3), "wrong input! m should be less than l");
128  utlException((l1%2!=0) || (l2%2!=0) || (l3%2!=0), "l1, l2, l3 need to be even integers" );
129  PreciseType result = 0;
130  // Sqrt[2]*Re[SphericalHarmonicY[l, Abs[m], \[Theta], \[Phi]]] = 1/Sqrt[2]*(SphericalHarmonicY[l,-m,Theta,Phi]+SphericalHarmonicY[l,m,Theta,Phi])
131  if (is_precalculated)
132  {
133  utlException(!utl::IsFileExist(utl::SH3Itegralhdr), "no SH3Itegralhdr.");
134  utl::InitializeSHTripleIntegrationTable();
135  int rank = utl::DimToRankSH(utl::SH3IntegralTable->GetShape()[0]);
136  if (l1<=rank && l2<=rank && l3<=rank)
137  {
138  unsigned index[3];
139  index[0] = utl::GetIndexSHj(l1,m1);
140  index[1] = utl::GetIndexSHj(l2,m2);
141  index[2] = utl::GetIndexSHj(l3,m3);
142  result = (*utl::SH3IntegralTable)(index);
143  }
144  else
145  result = RealTripleIntegration(l1,m1,l2,m2,l3,m3,false);
146  }
147  else
148  {
149  // NOTE: order 0 -2 -1 0 1 2 -4 -3 -2 -1 0 1 2 3 4 ...
150  // if m<0, Y[l,m] = Sqrt[2] * Re[SHY[l,Abs[m]]]
151  // if m=0, Y[l,m] = Re[SHY[l,m]], (imaginary part==0)
152  // if m>0, Y[l,m] = Sqrt[2] * Im[SHY[l,m]]
153  if (m1<0 && m2<0 && m3<0) // result = Sqrt[2]/4* Integrate[ ((-1)^(-m1)*Y[l1,m1]+Y[l1,-m1]) * ((-1)^(-m2)*Y[l2,m2]+Y[l2,-m2]) * ((-1)^(-m3)*Y[l3,m3]+Y[l3,-m3]) ]
154  {
155  if (m1+m2==m3)
156  result = ( m3%2==0?1:-1 ) * ComplexTripleIntegration(l1,m1,l2,m2,l3,-m3) / utl::SQRT2; // (-1)^(m1)*(-1)^(m2)*Integrate[Y[l1,m1]*Y[l2,m2]*Y[l3,-m3]]==(-1)^m3*Integrate[Y[l1,-m1]*Y[l2,-m2]*Y[l3,m3]]
157  else if (m1+m3==m2)
158  result = ( m2%2==0?1:-1 ) * ComplexTripleIntegration(l1,m1,l2,-m2,l3,m3) / utl::SQRT2;
159  else if (m2+m3==m1)
160  {
161  // utlPrintVar3(true,l1,l2,l3);
162  // utlPrintVar3(true,m1,m2,m3);
163  // utlPrintVar3(true,m1%2, ( m1%2==1?-1:1 ),ComplexTripleIntegration(l1,-m1,l2,m2,l3,m3));
164  // utlPrintVar1(true, ComplexTripleIntegration(l1,-m1,l2,m2,l3,m3));
165  result = ( m1%2==0?1:-1 ) * ComplexTripleIntegration(l1,-m1,l2,m2,l3,m3) / utl::SQRT2;
166  }
167  }
168  else if ( (m1>0 && m2<0 && m3<0) || (m1<0 && m2>0 && m3<0) || (m1<0 && m2<0 && m3>0) ) // Integrate[Sin[m1*Phi]*Cos[m2*Phi]*Cos[m3*Phi]]=0
169  {
170  result = 0;
171  }
172  else if (m1>0 && m2>0 && m3<0) // result = -Sqrt[2]/4* Integrate[ (Y[l1,m1]-(-1)^(m1)*Y[l1,-m1]) * (Y[l2,m2]-(-1)^(-m2)*Y[l2,-m2]) * ((-1)^(-m3)*Y[l3,m3]+Y[l3,-m3]) ]
173  {
174  if (m1+m2+m3==0)
175  result = (-1.0) * ( m3%2==0?1:-1 ) * ComplexTripleIntegration(l1,m1,l2,m2,l3,m3) / utl::SQRT2; // (-1)^(m1)*(-1)^(m2)*Integrate[Y[l1,m1]*Y[l2,m2]*Y[l3,-m3]]==(-1)^m3*Integrate[Y[l1,-m1]*Y[l2,-m2]*Y[l3,m3]]
176  else if (m1==m2+m3)
177  result = ( m2%2==0?1:-1 ) * ComplexTripleIntegration(l1,m1,l2,-m2,l3,-m3) / utl::SQRT2;
178  else if (m2==m1+m3)
179  result = ( m1%2==0?1:-1 ) * ComplexTripleIntegration(l1,-m1,l2,m2,l3,-m3) / utl::SQRT2;
180  }
181  else if (m1>0 && m2<0 && m3>0)
182  {
183  if (m1+m2+m3==0)
184  result = (-1.0) * ( m2%2==0?1:-1 ) * ComplexTripleIntegration(l1,m1,l2,m2,l3,m3) / utl::SQRT2; // (-1)^(m1)*(-1)^(m2)*Integrate[Y[l1,m1]*Y[l2,m2]*Y[l3,-m3]]==(-1)^m3*Integrate[Y[l1,-m1]*Y[l2,-m2]*Y[l3,m3]]
185  else if (m1==m2+m3)
186  result = ( m3%2==0?1:-1 ) * ComplexTripleIntegration(l1,m1,l2,-m2,l3,-m3) / utl::SQRT2;
187  else if (m3==m1+m2)
188  result = ( m1%2==0?1:-1 ) * ComplexTripleIntegration(l1,-m1,l2,-m2,l3,m3) / utl::SQRT2;
189  }
190  else if (m1<0 && m2>0 && m3>0)
191  {
192  if (m1+m2+m3==0)
193  result = (-1.0) * ( m1%2==0?1:-1 ) * ComplexTripleIntegration(l1,m1,l2,m2,l3,m3) / utl::SQRT2; // (-1)^(m1)*(-1)^(m2)*Integrate[Y[l1,m1]*Y[l2,m2]*Y[l3,-m3]]==(-1)^m3*Integrate[Y[l1,-m1]*Y[l2,-m2]*Y[l3,m3]]
194  else if (m3==m2+m1)
195  result = ( m2%2==0?1:-1 ) * ComplexTripleIntegration(l1,-m1,l2,-m2,l3,m3) / utl::SQRT2;
196  else if (m2==m1+m3)
197  result = ( m3%2==0?1:-1 ) * ComplexTripleIntegration(l1,-m1,l2,m2,l3,-m3) / utl::SQRT2;
198  }
199  else if (m1>0 && m2>0 && m3>0) // Integrate[Sin[m1*Phi]*Sin[m2*Phi]*Sin[m3*Phi]]=0
200  result = 0;
201  else if (m1==0 && m2==0 && m3==0)
202  result = ComplexTripleIntegration(l1,0,l2,0,l3,0);
203  else if (m1==0 && m2==m3)
204  result = ( m3%2==0?1:-1 ) * ComplexTripleIntegration(l1,0,l2,m2,l3,-m3);
205  else if (m2==0 && m3==m1)
206  result = ( m1%2==0?1:-1 ) * ComplexTripleIntegration(l1,m1,l2,0,l3,-m3);
207  else if (m3==0 && m1==m2)
208  result = ( m2%2==0?1:-1 ) * ComplexTripleIntegration(l1,m1,l2,-m2,l3,0);
209  }
210  return result;
211 }
212 
213 template < class PreciseType >
214 std::complex<double>
215 SphericalHarmonicsGenerator<PreciseType>
216 ::ComplexDerivativeOfTheta ( const int l, const int m, const double theta, const double phi)
217 {
218  std::complex<double> result(0,0), ee(0,-phi);
219  std::complex<double> sh_value = ComplexSH(l, m, theta, phi);
220  double theta_used = std::abs(theta)>1e-6? theta : (theta>=0?1e-6:-1e-6);
221  if (m==l)
222  return l * 1.0/std::tan(theta_used) * sh_value;
223  std::complex<double> sh_value_2 = ComplexSH(l, m+1, theta_used, phi);
224  result = m * 1.0/std::tan(theta_used) * sh_value + std::sqrt( utl::Gamma(1+l-m) * utl::Gamma(2+l+m) / (utl::Gamma(l-m)*utl::Gamma(l+m+1)) ) * sh_value_2 * ee;
225  return result;
226 }
227 
228 template < class PreciseType >
229 std::complex<double>
230 SphericalHarmonicsGenerator<PreciseType>
231 ::ComplexDerivativeOfPhi ( const int l, const int m, const double theta, const double phi)
232 {
233  std::complex<double> sh_value = ComplexSH(l, m, theta, phi);
234  std::complex<double> r1 (0,m);
235  return r1*sh_value;
236 }
237 
238 template < class PreciseType >
239 double
240 SphericalHarmonicsGenerator<PreciseType>
241 ::RealDerivativeOfTheta ( const int l, const int m, const double theta, const double phi)
242 {
243  if (m==0)
244  return std::real(ComplexDerivativeOfTheta(l,m,theta,phi));
245  else if (m>0)
246  return utl::SQRT2 * std::imag(ComplexDerivativeOfTheta(l,m,theta,phi));
247  else
248  return utl::SQRT2 * std::real(ComplexDerivativeOfTheta(l,-m,theta,phi));
249 }
250 
251 template < class PreciseType >
252 double
253 SphericalHarmonicsGenerator<PreciseType>
254 ::RealDerivativeOfPhi ( const int l, const int m, const double theta, const double phi)
255 {
256  if (m==0)
257  return std::real(ComplexDerivativeOfPhi(l,m,theta,phi));
258  else if (m>0)
259  return utl::SQRT2 * std::imag(ComplexDerivativeOfPhi(l,m,theta,phi));
260  else
261  return utl::SQRT2 * std::real(ComplexDerivativeOfPhi(l,-m,theta,phi));
262 }
263 
264 }
265 
266 #endif
utlDMRIStoredTables.h
Created "11-12-2015.
utl::InitializeSHTripleIntegrationTable
void InitializeSHTripleIntegrationTable(const int rank0=-1, const int rank1=-1, const int rank2=-1, const bool useExactSize=false)
Definition: itkSphericalHarmonicsGenerator.h:111
utl::SH3IntegralTable
static std::shared_ptr< NDArray< double, 3 > > SH3IntegralTable(new NDArray< double, 3 >)
utl::GetIndexSHj
int GetIndexSHj(const int l, const int m)
Definition: utlDMRI.h:204
utl::IsEven
bool IsEven(const int value)
Definition: utlCore.h:819
utlException
#define utlException(cond, expout)
Definition: utlCoreMacro.h:548
utl::IsFileExist
bool IsFileExist(const std::string &file)
Definition: utlCore.h:529
utl::DimToRankSH
int DimToRankSH(const int dimm)
Definition: utlDMRI.h:182
spams::abs
T abs(const T x)
template version of the fabs function
M_PI
#define M_PI
Definition: utlCoreMacro.h:57
utl
Definition: utl.h:90
itk::SphericalHarmonicsGenerator
Generate complex and real Spherical Harmonic Basis.
Definition: itkSphericalHarmonicsGenerator.h:54
utl::SH3Itegralhdr
static const std::string SH3Itegralhdr
Definition: DMRITOOLConfigure.h:60
utl::sign
int sign(const T &x)
Return the sign of x.
Definition: utlCore.h:285
utl::SQRT2
static constexpr double SQRT2
Definition: utlConstants.h:31
utl::Gamma
double Gamma(const double x)
gamma function.
Definition: itkSpecialFunctionGenerator.hxx:65
Generated by   doxygen 1.8.11