Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkDiffusionFunctionCollection.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
18 #include <math.h>
19 #include "mitkNumericTypes.h"
20 
21 // for Windows
22 #ifndef M_PI
23 #define M_PI 3.14159265358979323846
24 #endif
25 
26 // Namespace ::SH
27 #include <boost/math/special_functions/legendre.hpp>
28 #include <boost/math/special_functions/spherical_harmonic.hpp>
29 #include <boost/version.hpp>
30 
31 
32 // Namespace ::Gradients
33 #include "itkVectorContainer.h"
34 #include "vnl/vnl_vector.h"
35 
36 //------------------------- SH-function ------------------------------------
37 
38 double mitk::sh::factorial(int number) {
39  if(number <= 1) return 1;
40  double result = 1.0;
41  for(int i=1; i<=number; i++)
42  result *= i;
43  return result;
44 }
45 
46 void mitk::sh::Cart2Sph(double x, double y, double z, double *cart)
47 {
48  double phi, th, rad;
49  rad = sqrt(x*x+y*y+z*z);
50  if( rad < mitk::eps )
51  {
52  th = M_PI/2;
53  phi = M_PI/2;
54  }
55  else
56  {
57  th = acos(z/rad);
58  phi = atan2(y, x);
59  }
60  cart[0] = phi;
61  cart[1] = th;
62  cart[2] = rad;
63 }
64 
65 double mitk::sh::legendre0(int l)
66 {
67  if( l%2 != 0 )
68  {
69  return 0;
70  }
71  else
72  {
73  double prod1 = 1.0;
74  for(int i=1;i<l;i+=2) prod1 *= i;
75  double prod2 = 1.0;
76  for(int i=2;i<=l;i+=2) prod2 *= i;
77  return pow(-1.0,l/2.0)*(prod1/prod2);
78  }
79 }
80 
81 
82 double mitk::sh::Yj(int m, int l, double theta, double phi)
83 {
84  if (m<0)
85  return sqrt(2.0)*::boost::math::spherical_harmonic_r(l, -m, theta, phi);
86  else if (m==0)
87  return ::boost::math::spherical_harmonic_r(l, m, theta, phi);
88  else
89  return pow(-1.0,m)*sqrt(2.0)*::boost::math::spherical_harmonic_i(l, m, theta, phi);
90 
91  return 0;
92 }
93 
94 
95 
96 //------------------------- gradients-function ------------------------------------
97 
98 std::vector<unsigned int> mitk::gradients::GetAllUniqueDirections(const BValueMap & refBValueMap, GradientDirectionContainerType *refGradientsContainer )
99 {
100 
101  IndiciesVector directioncontainer;
102  auto mapIterator = refBValueMap.begin();
103 
104  if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1)
105  mapIterator++; //skip bzero Values
106 
107  for( ; mapIterator != refBValueMap.end(); mapIterator++){
108 
109  IndiciesVector currentShell = mapIterator->second;
110 
111  while(currentShell.size()>0)
112  {
113  unsigned int wntIndex = currentShell.back();
114  currentShell.pop_back();
115 
116  auto containerIt = directioncontainer.begin();
117  bool directionExist = false;
118  while(containerIt != directioncontainer.end())
119  {
120  if (fabs(dot(refGradientsContainer->ElementAt(*containerIt), refGradientsContainer->ElementAt(wntIndex))) > 0.9998)
121  {
122  directionExist = true;
123  break;
124  }
125  containerIt++;
126  }
127  if(!directionExist)
128  {
129  directioncontainer.push_back(wntIndex);
130  }
131  }
132  }
133 
134  return directioncontainer;
135 }
136 
137 
139 {
140  auto mapIterator = refBValueMap.begin();
141 
142  if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1)
143  mapIterator++; //skip bzero Values
144 
145  for( ; mapIterator != refBValueMap.end(); mapIterator++){
146 
147  auto mapIterator_2 = refBValueMap.begin();
148  if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1)
149  mapIterator_2++; //skip bzero Values
150 
151  for( ; mapIterator_2 != refBValueMap.end(); mapIterator_2++){
152 
153  if(mapIterator_2 == mapIterator) continue;
154 
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; }
159 
160  }
161  }
162  return false;
163 }
164 
165 
166 template<typename type>
167 double mitk::gradients::dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 )
168 {
169  double result = (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]) / (v1.two_norm() * v2.two_norm());
170  return result ;
171 }
172 
173 vnl_matrix<double> mitk::gradients::ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer)
174 {
175 
176  vnl_matrix<double> Q(3, refShell.size());
177  Q.fill(0.0);
178 
179  for(unsigned int i = 0; i < refShell.size(); i++)
180  {
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);
185  double cart[3];
186  mitk::sh::Cart2Sph(x,y,z,cart);
187  Q(0,i) = cart[0];
188  Q(1,i) = cart[1];
189  Q(2,i) = cart[2];
190  }
191  return Q;
192 }
193 
194 vnl_matrix<double> mitk::gradients::ComputeSphericalHarmonicsBasis(const vnl_matrix<double> & QBallReference, const unsigned int & LOrder)
195 {
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++)
201  {
202  int j = ( k * k + k + 2 ) / 2.0 + m - 1;
203  double phi = QBallReference(0,i);
204  double th = QBallReference(1,i);
205  double val = mitk::sh::Yj(m,k,th,phi);
206  SHBasisOutput(i,j) = val;
207  }
208  return SHBasisOutput;
209 }
210 
212  const GradientDirectionContainerType *origninalGradentcontainer)
213 {
215  auto mapIterator = bValueMap.begin();
216 
217  if(bValueMap.find(0) != bValueMap.end() && bValueMap.size() > 1){
218  mapIterator++; //skip bzero Values
219  vnl_vector_fixed<double, 3> vec;
220  vec.fill(0.0);
221  directioncontainer->push_back(vec);
222  }
223 
224  for( ; mapIterator != bValueMap.end(); mapIterator++){
225 
226  IndiciesVector currentShell = mapIterator->second;
227 
228  while(currentShell.size()>0)
229  {
230  unsigned int wntIndex = currentShell.back();
231  currentShell.pop_back();
232 
233  mitk::gradients::GradientDirectionContainerType::Iterator containerIt = directioncontainer->Begin();
234  bool directionExist = false;
235  while(containerIt != directioncontainer->End())
236  {
237  if (fabs(dot(containerIt.Value(), origninalGradentcontainer->ElementAt(wntIndex))) > 0.9998)
238  {
239  directionExist = true;
240  break;
241  }
242  containerIt++;
243  }
244  if(!directionExist)
245  {
246  GradientDirectionType dir(origninalGradentcontainer->ElementAt(wntIndex));
247  directioncontainer->push_back(dir.normalize());
248  }
249  }
250  }
251 
252  return directioncontainer;
253 }
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.