ANIMA  4.0
animaWatsonDistribution.hxx
Go to the documentation of this file.
1 #pragma once
2 #include <cmath>
3 
4 #include "animaWatsonDistribution.h"
5 #include <animaVectorOperations.h>
6 #include <animaErrorFunctions.h>
7 
8 #include <itkMacro.h>
9 #include <itkObjectFactory.h>
10 
11 #include <boost/math/special_functions/bessel.hpp>
12 
13 namespace anima
14 {
15 
16 template <class VectorType, class ScalarType>
17 double EvaluateWatsonPDF(const VectorType &v, const VectorType &meanAxis, const ScalarType &kappa)
18 {
19  /************************************************************************************************
20  * \fn template <class VectorType, class ScalarType>
21  * double
22  * EvaluateWatsonPDF(const VectorType &v,
23  * const VectorType &meanAxis,
24  * const ScalarType &kappa)
25  *
26  * \brief Evaluate the Watson probability density function using the definition of
27  * Fisher et al., Statistical Analysis of Spherical Data, 1993, p.89.
28  *
29  * \author Aymeric Stamm
30  * \date July 2014
31  *
32  * \param v Sample on which evaluating the PDF.
33  * \param meanAxis Mean axis of the Watson distribution.
34  * \param kappa Concentration parameter of the Watson distribution.
35  **************************************************************************************************/
36 
37  if (std::abs(kappa) < 1.0e-6)
38  return 1.0 / (4.0 * M_PI);
39  else if (kappa > 0)
40  {
41  double kappaSqrt = std::sqrt(kappa);
42  double c = anima::ComputeScalarProduct(v, meanAxis);
43  double inExp = kappa * (c * c - 1.0);
44  return kappaSqrt * std::exp(inExp) / (4.0 * M_PI * anima::EvaluateDawsonFunctionNR(kappaSqrt));
45  }
46  else
47  {
48  double Ck = std::sqrt(-kappa / M_PI) / (2.0 * M_PI * std::erf(std::sqrt(-kappa)));
49  double c = anima::ComputeScalarProduct(v, meanAxis);
50  double inExp = kappa * c * c;
51  return Ck * std::exp(inExp);
52  }
53 }
54 
55 template <class ScalarType>
56 double EvaluateWatsonPDF(const vnl_vector_fixed <ScalarType,3> &v, const vnl_vector_fixed <ScalarType,3> &meanAxis, const ScalarType &kappa)
57 {
58  if (std::abs(v.squared_magnitude() - 1.0) > 1.0e-6 || std::abs(meanAxis.squared_magnitude() - 1.0) > 1.0e-6)
59  throw itk::ExceptionObject(__FILE__, __LINE__,"The Watson distribution is on the 2-sphere.",ITK_LOCATION);
60 
61  return EvaluateWatsonPDF<vnl_vector_fixed <ScalarType,3>, ScalarType>(v,meanAxis,kappa);
62 }
63 
64 template <class ScalarType>
65 double EvaluateWatsonPDF(const itk::Point <ScalarType,3> &v, const itk::Point <ScalarType,3> &meanAxis, const ScalarType &kappa)
66 {
67  if (std::abs(v.GetVnlVector().squared_magnitude() - 1.0) > 1.0e-6 || std::abs(meanAxis.GetVnlVector().squared_magnitude() - 1.0) > 1.0e-6)
68  throw itk::ExceptionObject(__FILE__, __LINE__,"The Watson distribution is on the 2-sphere.",ITK_LOCATION);
69 
70  return EvaluateWatsonPDF<itk::Point <ScalarType,3>, ScalarType>(v,meanAxis,kappa);
71 }
72 
73 template <class ScalarType>
74 double EvaluateWatsonPDF(const itk::Vector <ScalarType,3> &v, const itk::Vector <ScalarType,3> &meanAxis, const ScalarType &kappa)
75 {
76  if (std::abs(v.GetSquaredNorm() - 1.0) > 1.0e-6 || std::abs(meanAxis.GetSquaredNorm() - 1.0) > 1.0e-6)
77  throw itk::ExceptionObject(__FILE__, __LINE__,"The Watson distribution is on the 2-sphere.",ITK_LOCATION);
78 
79  return EvaluateWatsonPDF<itk::Vector <ScalarType,3>, ScalarType>(v,meanAxis,kappa);
80 }
81 
82 template <class ScalarType>
83  void GetStandardWatsonSHCoefficients(const ScalarType k, std::vector<ScalarType> &coefficients, std::vector<ScalarType> &derivatives)
84 {
85  // Computes the first 7 non-zero SH coefficients of the standard Watson PDF (multiplied by 4 M_PI).
86  const unsigned int nbCoefs = 7;
87  coefficients.resize(nbCoefs);
88  derivatives.resize(nbCoefs);
89 
90  double sqrtPi = std::sqrt(M_PI);
91  double dawsonValue = anima::EvaluateDawsonIntegral(std::sqrt(k), true);
92  double k2 = k * k;
93  double k3 = k2 * k;
94  double k4 = k3 * k;
95  double k5 = k4 * k;
96  double k6 = k5 * k;
97  double k7 = k6 * k;
98 
99  coefficients[0] = 2.0 * sqrtPi;
100  derivatives[0] = 0.0;
101 
102  if (k < 1.0e-4) // Handles small k values by Taylor series expansion
103  {
104  coefficients[1] = k / 3.0;
105  coefficients[2] = k2 / 35.0;
106  coefficients[3] = k3 / 693.0;
107  coefficients[4] = k4 / 19305.0;
108  coefficients[5] = k5 / 692835.0;
109  coefficients[6] = k6 / 30421755.0;
110  derivatives[1] = 1.0 / 3.0;
111  derivatives[2] = 2.0 * k / 35.0;
112  derivatives[3] = k2 / 231.0;
113  derivatives[4] = 4.0 * k3 / 19305.0;
114  derivatives[5] = k4 / 138567.0;
115  derivatives[6] = 6.0 * k5 / 30421755.0;
116 
117  for (unsigned int i = 1;i < nbCoefs;++i)
118  {
119  double tmpVal = std::pow(2.0, i + 1.0) * sqrtPi / std::sqrt(4.0 * i + 1.0);
120  coefficients[i] *= tmpVal;
121  derivatives[i] *= tmpVal;
122  }
123 
124  return;
125  }
126 
127  coefficients[1] = (3.0 - (3.0 + 2.0 * k) * dawsonValue) / (2.0 * k);
128  coefficients[2] = (5.0 * (-21.0 + 2.0 * k) + 3.0 * (35.0 + 4.0 * k * (5.0 + k)) * dawsonValue) / (16.0 * k2);
129  coefficients[3] = (21.0 * (165.0 + 4.0 * (-5.0 + k) * k) - 5.0 * (693.0 + 378.0 * k + 84.0 * k2 + 8.0 * k3) * dawsonValue) / (64.0 * k3);
130  coefficients[4] = (3.0 * (-225225.0 + 2.0 * k * (15015.0 + 2.0 * k * (-1925.0 + 62.0 * k))) + 35.0 * (19305.0 + 8.0 * k * (1287.0 + k * (297.0 + 2.0 * k * (18.0 + k)))) * dawsonValue) / (1024.0 * k4);
131  coefficients[5] = (11.0 * (3968055.0 + 8.0 * k * (-69615.0 + 2.0 * k * (9828.0 + k * (-468.0 + 29.0 * k)))) - 63.0 * (692835.0 + 2.0 * k * (182325.0 + 4.0 * k * (10725.0 + 2.0 * k * (715.0 + k * (55.0 + 2.0 * k))))) * dawsonValue) / (4096.0 * k5);
132  coefficients[6] = (13.0 * (-540571185.0 + 2.0 * k * (39171825.0 + 4.0 * k * (-2909907.0 + 2.0 * k * (82467.0 + k * (-7469.0 + 122.0 * k))))) + 231.0 * (30421755.0 + 4.0 * k * (3968055.0 + k * (944775.0 + 4.0 * k * (33150.0 + k * (2925.0 + 4.0 * k * (39.0 + k)))))) * dawsonValue) / (32768.0 * k6);
133 
134  derivatives[1] = 3.0 * (-1.0 + (-1.0 + 2.0 * k) * dawsonValue + 2.0 * dawsonValue * dawsonValue) / (4.0 * k2);
135  derivatives[2] = 5.0 * ((21.0 - 2.0 * k) + (63.0 + 4.0 * (-11.0 + k) * k) * dawsonValue - 12.0 * (7.0 + 2.0 * k) * dawsonValue * dawsonValue) / (32.0 * k3);
136  derivatives[3] = 21.0 * ((-165.0 - 4.0 * (-5.0 + k) * k) + (-5.0 + 2.0 * k) * (165.0 + 4.0 * (-3.0 + k) * k) * dawsonValue + 10.0 * (99.0 + 4.0 * k * (9.0 + k)) * dawsonValue * dawsonValue) / (128.0 * k4);
137  derivatives[4] = 3.0 * ((225225.0 - 2.0 * k * (15015.0 + 2.0 * k * (-1925.0 + 62.0 * k))) + (1576575.0 + 8.0 * k * (-75075.0 + k * (10395.0 + 2.0 * k * (-978.0 + 31.0 * k)))) * dawsonValue - 840.0 * (2145.0 + 2.0 * k * (429.0 + 66.0 * k + 4.0 * k2)) * dawsonValue * dawsonValue) / (2048.0 * k5);
138  derivatives[5] = 11.0 * ((-3968055.0 - 8.0 * k * (-69615.0 + 2.0 * k * (9828.0 + k * (-468.0 + 29.0 * k)))) + (-35712495.0 + 2.0 * k * (5917275.0 + 8.0 * k * (-118755.0 + k * (21060.0 + k * (-965.0 + 58.0 * k))))) * dawsonValue + 630.0 * (62985.0 + 8.0 * k * (3315.0 + k * (585.0 + 2.0 * k * (26.0 + k)))) * dawsonValue * dawsonValue) / (8192.0 * k6);
139  derivatives[6] = 13.0 * ((540571185.0 - 2.0 * k * (39171825.0 + 4.0 * k * (-2909907.0 + 2.0 * k * (82467.0 + k * (-7469.0 + 122.0 * k))))) + (5946283035.0 + 4.0 * k * (-446558805.0 + k * (79910523.0 + 4.0 * k * (-3322242.0 + k * (187341.0 + 4.0 * k * (-3765.0 + 61.0 * k)))))) * dawsonValue - 2772.0 * (2340135.0 + 2.0 * k * (508725.0 + 4.0 * k * (24225.0 + 2.0 * k * (1275.0 + k * (75.0 + 2.0 * k))))) * dawsonValue * dawsonValue) / (65536.0 * k7);
140 
141  for (unsigned int i = 1;i < nbCoefs;++i)
142  {
143  double sqrtVal = std::sqrt(1.0 + 4.0 * i);
144  coefficients[i] *= (sqrtPi * sqrtVal / dawsonValue);
145  derivatives[i] *= (sqrtPi * sqrtVal / (dawsonValue * dawsonValue));
146  }
147 }
148 
149 } // end namespace anima
anima::EvaluateWatsonPDF
double EvaluateWatsonPDF(const VectorType &v, const VectorType &meanAxis, const ScalarType &kappa)
Definition: animaWatsonDistribution.hxx:17
anima::GetStandardWatsonSHCoefficients
void GetStandardWatsonSHCoefficients(const ScalarType k, std::vector< ScalarType > &coefficients, std::vector< ScalarType > &derivatives)
Definition: animaWatsonDistribution.hxx:83
anima::ComputeScalarProduct
double ComputeScalarProduct(const VectorType &v1, const VectorType &v2, const unsigned int NDimension)
Definition: animaVectorOperations.hxx:196
anima::EvaluateDawsonIntegral
double EvaluateDawsonIntegral(const double x, const bool scaled)
Definition: animaErrorFunctions.cxx:11
anima::EvaluateDawsonFunctionNR
double EvaluateDawsonFunctionNR(double x)
Numerical recipes implementation of Dawson integral.
Definition: animaErrorFunctions.cxx:20
Generated by   doxygen 1.8.13