ANIMA  4.0
animaLegendreDerivatives.cxx
Go to the documentation of this file.
2 
3 #include <itkMacro.h>
4 #include <boost/math/special_functions/legendre.hpp>
5 
6 namespace anima
7 {
8 
9 double legendre_first_derivative(int L, int M, double value)
10 {
11  if (L < 1)
12  {
13  throw itk::ExceptionObject(__FILE__, __LINE__,
14  "Legendre polynomials first derivative only implemented for L > 1",ITK_LOCATION);
15  }
16 
17  // Handle negative M, not sure if derivative is valid for negative M
18  double factor = 1;
19  if (M < 0)
20  {
21  M = abs(M);
22  factor = std::tgamma(L - M + 1) / std::tgamma (L + M + 1);
23  if (M%2 != 0)
24  factor *= -1;
25  }
26 
27  double sqValue = value * value;
28  if (fabs (value) == 1)
29  {
30  if (value > 0)
31  value -= 1.0e-16;
32  else
33  value += 1.0e-16;
34 
35  sqValue = value * value;
36  }
37 
38  double resVal = L * value * boost::math::legendre_p(L,M,value) - (L+M) * boost::math::legendre_p(L-1,M,value);
39 
40  resVal *= factor / (sqValue - 1.0);
41 
42  return resVal;
43 }
44 
45 double legendre_second_derivative(int L, int M, double value)
46 {
47  if (L < 2)
48  {
49  throw itk::ExceptionObject(__FILE__, __LINE__,
50  "Legendre polynomials second derivative only implemented for L > 2",ITK_LOCATION);
51  }
52 
53  // Handle negative M, not sure if derivative is valid for negative M
54  double factor = 1;
55  if (M < 0)
56  {
57  M = abs(M);
58  factor = std::tgamma (L - M + 1) / std::tgamma (L + M + 1);
59  if (M%2 != 0)
60  factor *= -1;
61  }
62 
63  double sqValue = value * value;
64  if (fabs (value) == 1)
65  {
66  if (value > 0)
67  value -= 1.0e-16;
68  else
69  value += 1.0e-16;
70 
71  sqValue = value * value;
72  }
73 
74  double resVal = L * ((L - 1) * sqValue - 1) * boost::math::legendre_p(L,M,value);
75  resVal += (L + M) * (3 - 2*L) * value * boost::math::legendre_p(L-1,M,value);
76  resVal += (L + M) * (L + M - 1) * boost::math::legendre_p(L-2,M,value);
77 
78  resVal *= factor / ((sqValue - 1.0) * (sqValue - 1.0));
79 
80  return resVal;
81 }
82 
83 } // end of namespace anima
double legendre_second_derivative(int L, int M, double value)
double legendre_first_derivative(int L, int M, double value)