23 #define M_PI 3.14159265358979323846
27 #include <boost/math/special_functions/legendre.hpp>
28 #include <boost/math/special_functions/spherical_harmonic.hpp>
29 #include <boost/version.hpp>
33 #include "itkVectorContainer.h"
34 #include "vnl/vnl_vector.h"
39 if(number <= 1)
return 1;
41 for(
int i=1; i<=number; i++)
49 rad = sqrt(x*x+y*y+z*z);
74 for(
int i=1;i<l;i+=2) prod1 *= i;
76 for(
int i=2;i<=l;i+=2) prod2 *= i;
77 return pow(-1.0,l/2.0)*(prod1/prod2);
85 return sqrt(2.0)*::boost::math::spherical_harmonic_r(l, -m, theta, phi);
87 return ::boost::math::spherical_harmonic_r(l, m, theta, phi);
89 return pow(-1.0,m)*sqrt(2.0)*::boost::math::spherical_harmonic_i(l, m, theta, phi);
101 IndiciesVector directioncontainer;
102 auto mapIterator = refBValueMap.begin();
104 if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1)
107 for( ; mapIterator != refBValueMap.end(); mapIterator++){
109 IndiciesVector currentShell = mapIterator->second;
111 while(currentShell.size()>0)
113 unsigned int wntIndex = currentShell.back();
114 currentShell.pop_back();
116 auto containerIt = directioncontainer.begin();
117 bool directionExist =
false;
118 while(containerIt != directioncontainer.end())
120 if (fabs(dot(refGradientsContainer->ElementAt(*containerIt), refGradientsContainer->ElementAt(wntIndex))) > 0.9998)
122 directionExist =
true;
129 directioncontainer.push_back(wntIndex);
134 return directioncontainer;
140 auto mapIterator = refBValueMap.begin();
142 if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1)
145 for( ; mapIterator != refBValueMap.end(); mapIterator++){
147 auto mapIterator_2 = refBValueMap.begin();
148 if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1)
151 for( ; mapIterator_2 != refBValueMap.end(); mapIterator_2++){
153 if(mapIterator_2 == mapIterator)
continue;
155 IndiciesVector currentShell = mapIterator->second;
156 IndiciesVector testShell = mapIterator_2->second;
157 for (
unsigned int i = 0; i< currentShell.size(); i++)
158 if (fabs(dot(refGradientsContainer->ElementAt(currentShell[i]), refGradientsContainer->ElementAt(testShell[i]))) <= 0.9998) {
return true; }
166 template<
typename type>
169 double result = (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]) / (v1.two_norm() * v2.two_norm());
176 vnl_matrix<double> Q(3, refShell.size());
179 for(
unsigned int i = 0; i < refShell.size(); i++)
181 GradientDirectionType dir = refGradientsContainer->ElementAt(refShell[i]);
182 double x = dir.normalize().get(0);
183 double y = dir.normalize().get(1);
184 double z = dir.normalize().get(2);
196 vnl_matrix<double> SHBasisOutput(QBallReference.cols(), (LOrder+1)*(LOrder+2)*0.5);
197 SHBasisOutput.fill(0.0);
198 for(
int i=0; i< (int)SHBasisOutput.rows(); i++)
199 for(
int k = 0; k <= (int)LOrder; k += 2)
200 for(
int m =- k; m <= k; m++)
202 int j = ( k * k + k + 2 ) / 2.0 + m - 1;
203 double phi = QBallReference(0,i);
204 double th = QBallReference(1,i);
206 SHBasisOutput(i,j) = val;
208 return SHBasisOutput;
212 const GradientDirectionContainerType *origninalGradentcontainer)
215 auto mapIterator = bValueMap.begin();
217 if(bValueMap.find(0) != bValueMap.end() && bValueMap.size() > 1){
219 vnl_vector_fixed<double, 3> vec;
221 directioncontainer->push_back(vec);
224 for( ; mapIterator != bValueMap.end(); mapIterator++){
226 IndiciesVector currentShell = mapIterator->second;
228 while(currentShell.size()>0)
230 unsigned int wntIndex = currentShell.back();
231 currentShell.pop_back();
233 mitk::gradients::GradientDirectionContainerType::Iterator containerIt = directioncontainer->Begin();
234 bool directionExist =
false;
235 while(containerIt != directioncontainer->End())
237 if (fabs(dot(containerIt.Value(), origninalGradentcontainer->ElementAt(wntIndex))) > 0.9998)
239 directionExist =
true;
246 GradientDirectionType dir(origninalGradentcontainer->ElementAt(wntIndex));
247 directioncontainer->push_back(dir.normalize());
252 return directioncontainer;
itk::SmartPointer< Self > Pointer
static double legendre0(int l)
static double Yj(int m, int k, double theta, double phi)
static double dot(vnl_vector_fixed< type, 3 > const &v1, vnl_vector_fixed< type, 3 > const &v2)
static void Cart2Sph(double x, double y, double z, double *cart)
static double factorial(int number)
itk::SmartPointer< const Self > ConstPointer
static vnl_matrix< double > ComputeSphericalHarmonicsBasis(const vnl_matrix< double > &QBallReference, const unsigned int &LOrder)
static bool CheckForDifferingShellDirections(const BValueMap &bValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer)
static mitk::gradients::GradientDirectionContainerType::Pointer CreateNormalizedUniqueGradientDirectionContainer(const BValueMap &bValueMap, const GradientDirectionContainerType *origninalGradentcontainer)
MITKCORE_EXPORT const ScalarType eps
static std::vector< unsigned int > GetAllUniqueDirections(const BValueMap &bValueMap, GradientDirectionContainerType *refGradientsContainer)
static vnl_matrix< double > ComputeSphericalFromCartesian(const IndiciesVector &refShell, const GradientDirectionContainerType *refGradientsContainer)
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.