Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkDiffusionMultiShellQballReconstructionImageFilter.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 #ifndef __itkDiffusionMultiShellQballReconstructionImageFilter_cpp
17 #define __itkDiffusionMultiShellQballReconstructionImageFilter_cpp
18 
19 #include <itkDiffusionMultiShellQballReconstructionImageFilter.h>
20 
21 #include <itkTimeProbe.h>
22 #include <itkPointShell.h>
23 #include <mitkDiffusionFunctionCollection.h>
24 
25 namespace itk {
26 
27 template< class T, class TG, class TO, int L, int NODF>
28 DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
29 ::DiffusionMultiShellQballReconstructionImageFilter() :
30  m_ReconstructionType(Mode_Standard1Shell),
31  m_Interpolation_Flag(false),
32  m_Interpolation_SHT1_inv(NULL),
33  m_Interpolation_SHT2_inv(NULL),
34  m_Interpolation_SHT3_inv(NULL),
35  m_TARGET_SH_shell1(NULL),
36  m_TARGET_SH_shell2(NULL),
37  m_TARGET_SH_shell3(NULL),
38  m_MaxDirections(0),
39  m_CoeffReconstructionMatrix(NULL),
40  m_ODFSphericalHarmonicBasisMatrix(NULL),
41  m_GradientDirectionContainer(NULL),
42  m_NumberOfGradientDirections(0),
43  m_NumberOfBaselineImages(0),
44  m_Threshold(0),
45  m_BZeroImage(NULL),
46  m_CoefficientImage(NULL),
47  m_BValue(1.0),
48  m_Lambda(0.0),
49  m_IsHemisphericalArrangementOfGradientDirections(false),
50  m_IsArithmeticProgession(false)
51 {
52  // At least 1 inputs is necessary for a vector image.
53  // For images added one at a time we need at least six
54  this->SetNumberOfRequiredInputs( 1 );
55 }
56 
57 
58 template< class T, class TG, class TO, int L, int NODF>
59 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
60 ::SetGradientImage(const GradientDirectionContainerType *gradientDirection
61  , const GradientImagesType *gradientImage
62  , float bvalue)
63 {
64  m_BValue = bvalue;
65  m_NumberOfBaselineImages = 0;
66 
67  this->m_GradientDirectionContainer = GradientDirectionContainerType::New();
68  for(GradientDirectionContainerType::ConstIterator it = gradientDirection->Begin();
69  it != gradientDirection->End(); it++)
70  {
71  this->m_GradientDirectionContainer->push_back(it.Value());
72  }
73 
74 
75  if(m_BValueMap.size() == 0){
76  itkWarningMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : no GradientIndexMapAvalible");
77 
78  GradientDirectionContainerType::ConstIterator gdcit;
79  for( gdcit = m_GradientDirectionContainer->Begin(); gdcit != m_GradientDirectionContainer->End(); ++gdcit)
80  {
81  double bValueKey = int(((m_BValue * gdcit.Value().two_norm() * gdcit.Value().two_norm())+7.5)/10)*10;
82  m_BValueMap[bValueKey].push_back(gdcit.Index());
83  }
84 
85  }
86  if(m_BValueMap.find(0) == m_BValueMap.end())
87  {
88  itkExceptionMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : GradientIndxMap with no b-Zero indecies found: check input BValueMap");
89  }
90 
91 
92  m_NumberOfBaselineImages = m_BValueMap[0].size();
93  m_NumberOfGradientDirections = gradientDirection->Size() - m_NumberOfBaselineImages;
94  // ensure that the gradient image we received has as many components as
95  // the number of gradient directions
96  if( gradientImage->GetVectorLength() != m_NumberOfBaselineImages + m_NumberOfGradientDirections )
97  {
98  itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << m_NumberOfBaselineImages
99  << "baselines = " << m_NumberOfGradientDirections + m_NumberOfBaselineImages
100  << " directions specified but image has " << gradientImage->GetVectorLength()
101  << " components.");
102  }
103 
104 
105  ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) );
106 
107  std::string gradientImageClassName(ProcessObject::GetInput(0)->GetNameOfClass());
108  if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 )
109  itkExceptionMacro( << "There is only one Gradient image. I expect that to be a VectorImage. But its of type: " << gradientImageClassName );
110 
111 
112  m_BZeroImage = BZeroImageType::New();
113  typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >( ProcessObject::GetInput(0) );
114  m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing
115  m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin
116  m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction
117  m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
118  m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
119  m_BZeroImage->Allocate();
120 
121  m_CoefficientImage = CoefficientImageType::New();
122  m_CoefficientImage->SetSpacing( img->GetSpacing() ); // Set the image spacing
123  m_CoefficientImage->SetOrigin( img->GetOrigin() ); // Set the image origin
124  m_CoefficientImage->SetDirection( img->GetDirection() ); // Set the image direction
125  m_CoefficientImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
126  m_CoefficientImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
127  m_CoefficientImage->Allocate();
128 
129 }
130 
131 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
132 void DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
133 ::Normalize( OdfPixelType & out)
134 {
135  for(int i=0; i<NrOdfDirections; i++)
136  {
137  out[i] = out[i] < 0 ? 0 : out[i];
138  out[i] *= M_PI * 4 / NrOdfDirections;
139  }
140 }
141 
142 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
143 void DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
144 ::Projection1(vnl_vector<double> & vec, double delta)
145 {
146  if (delta==0){ //Clip attenuation values. If att<0 => att=0, if att>1 => att=1
147  for (unsigned int i=0; i<vec.size(); i++)
148  vec[i]=(vec[i]>=0 && vec[i]<=1)*vec[i]+(vec[i]>1);
149  }
150  else{ //Use function from Aganj et al, MRM, 2010
151  for (unsigned int i=0; i< vec.size(); i++)
152  vec[i]=CalculateThreashold(vec[i], delta);
153  }
154 }
155 
156 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
157 double DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
158 ::CalculateThreashold(const double value, const double delta)
159 {
160  return (value<0)*(0.5*delta) + (value>=0 && value<delta)*(0.5*delta+0.5*(value*value)/delta)
161  + (value>=delta && value<1-delta)*value+(value>=1-delta && value<1)*(1-0.5*delta-0.5*((1-value)*(1-value))/delta)
162  + (value>=1)*(1-0.5*delta);
163 }
164 
165 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
166 void DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
167 ::Projection2( vnl_vector<double> & E1,vnl_vector<double> & E2, vnl_vector<double> & E3, double delta )
168 {
169 
170  const double sF = sqrt(5.0);
171 
172  vnl_vector<double> vOnes(m_MaxDirections);
173  vOnes.fill(1.0);
174 
175  vnl_matrix<double> T0(m_MaxDirections, 3);
176  vnl_matrix<unsigned char> C(m_MaxDirections, 7);
177  vnl_matrix<double> A(m_MaxDirections, 7);
178  vnl_matrix<double> B(m_MaxDirections, 7);
179 
180  vnl_vector<double> s0(m_MaxDirections);
181  vnl_vector<double> a0(m_MaxDirections);
182  vnl_vector<double> b0(m_MaxDirections);
183  vnl_vector<double> ta(m_MaxDirections);
184  vnl_vector<double> tb(m_MaxDirections);
185  vnl_vector<double> e(m_MaxDirections);
186  vnl_vector<double> m(m_MaxDirections);
187  vnl_vector<double> a(m_MaxDirections);
188  vnl_vector<double> b(m_MaxDirections);
189 
190 
191  // logarithmierung aller werte in E
192  for(unsigned int i = 0 ; i < m_MaxDirections; i++)
193  {
194  T0(i,0) = -log(E1(i));
195  T0(i,1) = -log(E2(i));
196  T0(i,2) = -log(E3(i));
197  }
198 
199  //T0 = -T0.apply(std::log);
200 
201 
202  // Summeiere Zeilenweise über alle Shells sum = E1+E2+E3
203  for(unsigned int i = 0 ; i < m_MaxDirections; i++)
204  {
205  s0[i] = T0(i,0) + T0(i,1) + T0(i,2);
206  }
207 
208  for(unsigned int i = 0; i < m_MaxDirections; i ++)
209  {
210  // Alle Signal-Werte auf der Ersten shell E(N,0) normiert auf s0
211  a0[i] = T0(i,0) / s0[i];
212  // Alle Signal-Werte auf der Zweiten shell E(N,1) normiert auf s0
213  b0[i] = T0(i,1) / s0[i];
214  }
215 
216  ta = a0 * 3.0;
217  tb = b0 * 3.0;
218  e = tb - (ta * 2.0);
219  m = (tb * 2.0 ) + ta;
220 
221  for(unsigned int i = 0; i <m_MaxDirections; i++)
222  {
223  C(i,0) = (tb[i] < 1+3*delta && 0.5+1.5*(sF+1)*delta < ta[i] && ta[i] < 1-3* (sF+2) *delta);
224  C(i,1) = (e[i] <= -1+3*(2*sF+5)*delta) && (ta[i]>=1-3*(sF+2)*delta);
225  C(i,2) = (m[i] > 3-3*sF*delta) && (-1+3*(2*sF+5)*delta<e[i]) && (e[i]<-3*sF*delta);
226  C(i,3) = (m[i] >= 3-3*sF*delta && e[i] >= -3 *sF * delta);
227  C(i,4) = (2.5 + 1.5*(5+sF)*delta < m[i] && m[i] < 3-3*sF*delta && e[i] > -3*sF*delta);
228  C(i,5) = (ta[i] <= 0.5+1.5 *(sF+1)*delta && m[i] <= 2.5 + 1.5 *(5+sF) * delta);
229  C(i,6) = !((bool) C(i,0) ||(bool) C(i,1) ||(bool) C(i,2) ||(bool) C(i,3) ||(bool) C(i,4) ||(bool) C(i,5) ); // ~ANY(C(i,[0-5] ),2)
230 
231  A(i,0)=(bool)C(i,0) * a0(i);
232  A(i,1)=(bool)C(i,1) * (1.0/3.0-(sF+2)*delta);
233  A(i,2)=(bool)C(i,2) * (0.2+0.8*a0(i)-0.4*b0(i)-delta/sF);
234  A(i,3)=(bool)C(i,3) * (0.2+delta/sF);
235  A(i,4)=(bool)C(i,4) * (0.2*a0(i)+0.4*b0(i)+2*delta/sF);
236  A(i,5)=(bool)C(i,5) * (1.0/6.0+0.5*(sF+1)*delta);
237  A(i,6)=(bool)C(i,6) * a0(i);
238 
239  B(i,0)=(bool)C(i,0) * (1.0/3.0+delta);
240  B(i,1)=(bool)C(i,1) * (1.0/3.0+delta);
241  B(i,2)=(bool)C(i,2) * (0.4-0.4*a0(i)+0.2*b0(i)-2*delta/sF);
242  B(i,3)=(bool)C(i,3) * (0.4-3*delta/sF);
243  B(i,4)=(bool)C(i,4) * (0.4*a0(i)+0.8*b0(i)-delta/sF);
244  B(i,5)=(bool)C(i,5) * (1.0/3.0+delta);
245  B(i,6)=(bool)C(i,6) * b0(i);
246  }
247 
248  for(unsigned int i = 0 ; i < m_MaxDirections; i++)
249  {
250  double sumA = 0;
251  double sumB = 0;
252  for(int j = 0 ; j < 7; j++)
253  {
254  sumA += A(i,j);
255  sumB += B(i,j);
256  }
257  a[i] = sumA;
258  b[i] = sumB;
259  }
260 
261  for(unsigned int i = 0; i < m_MaxDirections; i++)
262  {
263  E1(i) = exp(-(a[i]*s0[i]));
264  E2(i) = exp(-(b[i]*s0[i]));
265  E3(i) = exp(-((1-a[i]-b[i])*s0[i]));
266  }
267 
268 }
269 
270 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
271 void DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
272 ::Projection3( vnl_vector<double> & A, vnl_vector<double> & a, vnl_vector<double> & b, double delta0)
273 {
274 
275  const double s6 = sqrt(6.0);
276  const double s15 = s6/2.0;
277 
278  vnl_vector<double> delta(a.size());
279  delta.fill(delta0);
280 
281  vnl_matrix<double> AM(a.size(), 15);
282  vnl_matrix<double> aM(a.size(), 15);
283  vnl_matrix<double> bM(a.size(), 15);
284  vnl_matrix<unsigned char> B(a.size(), 15);
285 
286  AM.set_column(0, A);
287  AM.set_column(1, A);
288  AM.set_column(2, A);
289  AM.set_column(3, delta);
290  AM.set_column(4, (A+a-b - (delta*s6))/3.0);
291  AM.set_column(5, delta);
292  AM.set_column(6, delta);
293  AM.set_column(7, delta);
294  AM.set_column(8, A);
295  AM.set_column(9, 0.2*(a*2+A-2*(s6+1)*delta));
296  AM.set_column(10,0.2*(b*(-2)+A+2-2*(s6+1)*delta));
297  AM.set_column(11, delta);
298  AM.set_column(12, delta);
299  AM.set_column(13, delta);
300  AM.set_column(14, 0.5-(1+s15)*delta);
301 
302 
303  aM.set_column(0, a);
304  aM.set_column(1, a);
305  aM.set_column(2, -delta + 1);
306  aM.set_column(3, a);
307  aM.set_column(4, (A*2+a*5+b+s6*delta)/6.0);
308  aM.set_column(5, a);
309  aM.set_column(6, -delta + 1);
310  aM.set_column(7, 0.5*(a+b)+(1+s15)*delta);
311  aM.set_column(8, -delta + 1);
312  aM.set_column(9, 0.2*(a*4+A*2+(s6+1)*delta));
313  aM.set_column(10, -delta + 1);
314  aM.set_column(11, (s6+3)*delta);
315  aM.set_column(12, -delta + 1);
316  aM.set_column(13, -delta + 1);
317  aM.set_column(14, -delta + 1);
318 
319  bM.set_column(0, b);
320  bM.set_column(1, delta);
321  bM.set_column(2, b);
322  bM.set_column(3, b);
323  bM.set_column(4, (A*(-2)+a+b*5-s6*delta)/6.0);
324  bM.set_column(5, delta);
325  bM.set_column(6, b);
326  bM.set_column(7, 0.5*(a+b)-(1+s15)*delta);
327  bM.set_column(8, delta);
328  bM.set_column(9, delta);
329  bM.set_column(10, 0.2*(b*4-A*2+1-(s6+1)*delta));
330  bM.set_column(11, delta);
331  bM.set_column(12, delta);
332  bM.set_column(13, -delta*(s6+3) + 1);
333  bM.set_column(14, delta);
334 
335  delta0 *= 0.99;
336 
337  vnl_matrix<double> R2(a.size(), 15);
338  std::vector<unsigned int> I(a.size());
339 
340  for (unsigned int i=0; i<AM.rows(); i++){
341  for (unsigned int j=0; j<AM.cols(); j++){
342  if (delta0 < AM(i,j) && 2*(AM(i,j)+delta0*s15)<aM(i,j)-bM(i,j) && bM(i,j)>delta0 && aM(i,j)<1-delta0)
343  R2(i,j) = (AM(i,j)-A(i))*(AM(i,j)-A(i))+ (aM(i,j)-a(i))*(aM(i,j)-a(i))+(bM(i,j)-b(i))*(bM(i,j)-b(i));
344  else
345  R2(i,j) = 1e20;
346  }
347  unsigned int index = 0;
348  double minvalue = 999;
349  for(int j = 0 ; j < 15 ; j++)
350  {
351  if(R2(i,j) < minvalue){
352  minvalue = R2(i,j);
353  index = j;
354  }
355  }
356  I[i] = index;
357  }
358 
359  for (unsigned int i=0; i < A.size(); i++){
360  A(i) = AM(i,(int)I[i]);
361  a(i) = aM(i,(int)I[i]);
362  b(i) = bM(i,(int)I[i]);
363  }
364 
365 }
366 
367 template<class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
368 void DiffusionMultiShellQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType,NOrderL, NrOdfDirections>
369 ::S_S0Normalization( vnl_vector<double> & vec, double S0 )
370 {
371  for(unsigned int i = 0; i < vec.size(); i++)
372  {
373  if (S0==0)
374  S0 = 0.01;
375  vec[i] /= S0;
376  }
377 
378 }
379 
380 template< class T, class TG, class TO, int L, int NODF>
381 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
382 ::DoubleLogarithm(vnl_vector<double> & vec)
383 {
384  for(unsigned int i = 0; i < vec.size(); i++)
385  {
386  vec[i] = log(-log(vec[i]));
387  }
388 }
389 
390 
391 
392 template< class T, class TG, class TO, int L, int NODF>
393 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
394 ::BeforeThreadedGenerateData()
395 {
396  m_ReconstructionType = Mode_Standard1Shell;
397 
398  if(m_BValueMap.size() == 4 ){
399 
400  BValueMapIteraotr it = m_BValueMap.begin();
401  it++; // skip b0 entry
402  const unsigned int bValue_shell1 = it->first;
403  const unsigned int size_shell1 = it->second.size();
404  IndiciesVector shell1 = it->second;
405  it++;
406  const unsigned int bValue_shell2 = it->first;
407  const unsigned int size_shell2 = it->second.size();
408  IndiciesVector shell2 = it->second;
409  it++;
410  const unsigned int bValue_shell3 = it->first;
411  const unsigned int size_shell3 = it->second.size();
412  IndiciesVector shell3 = it->second;
413 
414  // arithmetic progrssion
415  if(bValue_shell2 - bValue_shell1 == bValue_shell1 && bValue_shell3 - bValue_shell2 == bValue_shell1 )
416  {
417  // check if Interpolation is needed
418  // if shells with different numbers of directions exist
419  m_Interpolation_Flag = false;
420  if(size_shell1 != size_shell2 || size_shell2 != size_shell3 || size_shell1 != size_shell3)
421  {
422  m_Interpolation_Flag = true;
423  MITK_INFO << "Shell interpolation: shells with different numbers of directions";
424  }
425  else
426  {
427  // else if each shell holds same numbers of directions, but the gradient direction differ more than one 1 degree
428  m_Interpolation_Flag = CheckForDifferingShellDirections();
429  if(m_Interpolation_Flag) MITK_INFO << "Shell interpolation: gradient direction differ more than one 1 degree";
430  }
431 
432  m_ReconstructionType = Mode_Analytical3Shells;
433 
434  if(m_Interpolation_Flag)
435  {
436  unsigned int interp_SHOrder_shell1 = 12;
437  while( ((interp_SHOrder_shell1+1)*(interp_SHOrder_shell1+2)/2) > size_shell1 && interp_SHOrder_shell1 > L )
438  interp_SHOrder_shell1 -= 2 ;
439 
440  const int number_coeffs_shell1 = (int)(interp_SHOrder_shell1*interp_SHOrder_shell1 + interp_SHOrder_shell1 + 2.0)/2.0 + interp_SHOrder_shell1;
441 
442  unsigned int interp_SHOrder_shell2 = 12;
443  while( ((interp_SHOrder_shell2+1)*(interp_SHOrder_shell2+2)/2) > size_shell2 && interp_SHOrder_shell2 > L )
444  interp_SHOrder_shell2 -= 2 ;
445 
446  const int number_coeffs_shell2 = (int)(interp_SHOrder_shell2*interp_SHOrder_shell2 + interp_SHOrder_shell2 + 2.0)/2.0 + interp_SHOrder_shell2;
447 
448  unsigned int interp_SHOrder_shell3 = 12;
449  while( ((interp_SHOrder_shell3+1)*(interp_SHOrder_shell3+2)/2) > size_shell3 && interp_SHOrder_shell3 > L )
450  interp_SHOrder_shell3 -= 2 ;
451 
452  const int number_coeffs_shell3 = (int)(interp_SHOrder_shell3*interp_SHOrder_shell3 + interp_SHOrder_shell3 + 2.0)/2.0 + interp_SHOrder_shell3;
453 
454  // Create direction container for all directions (no duplicates, different directions from all shells)
455  IndiciesVector all_directions_container = GetAllDirections();
456  m_MaxDirections = all_directions_container.size();
457 
458  // create target SH-Basis
459  // initialize empty target matrix and set the wanted directions
460  vnl_matrix<double> * Q = new vnl_matrix<double>(3, m_MaxDirections);
461  ComputeSphericalFromCartesian(Q, all_directions_container);
462  // initialize a SH Basis, neede to interpolate from oldDirs -> newDirs
463  m_TARGET_SH_shell1 = new vnl_matrix<double>(m_MaxDirections, number_coeffs_shell1);
464  ComputeSphericalHarmonicsBasis(Q, m_TARGET_SH_shell1, interp_SHOrder_shell1);
465  delete Q;
466 
467  Q = new vnl_matrix<double>(3, m_MaxDirections);
468  ComputeSphericalFromCartesian(Q, all_directions_container);
469  m_TARGET_SH_shell2 = new vnl_matrix<double>(m_MaxDirections, number_coeffs_shell2);
470  ComputeSphericalHarmonicsBasis(Q, m_TARGET_SH_shell2, interp_SHOrder_shell2);
471  delete Q;
472 
473  Q = new vnl_matrix<double>(3, m_MaxDirections);
474  ComputeSphericalFromCartesian(Q, all_directions_container);
475  m_TARGET_SH_shell3 = new vnl_matrix<double>(m_MaxDirections, number_coeffs_shell3);
476  ComputeSphericalHarmonicsBasis(Q, m_TARGET_SH_shell3, interp_SHOrder_shell3);
477  delete Q;
478  // end creat target SH-Basis
479 
480 
481  // create measured-SHBasis
482  // Shell 1
483  vnl_matrix<double> * tempSHBasis;
484  vnl_matrix_inverse<double> * temp;
485 
486  // initialize empty matrix and set the measured directions of shell1
487  Q = new vnl_matrix<double>(3, shell1.size());
488  ComputeSphericalFromCartesian(Q, shell1);
489  // initialize a SH Basis, need to get the coeffs from measuredShell
490  tempSHBasis = new vnl_matrix<double>(shell1.size(), number_coeffs_shell1);
491  ComputeSphericalHarmonicsBasis(Q, tempSHBasis, interp_SHOrder_shell1);
492  // inversion of the SH=Basis
493  // c_s1 = B^-1 * shell1
494  // interp_Values = targetSHBasis * c_s1
495 
496  // Values = m_TARGET_SH_shell1 * Interpolation_SHT1_inv * DataShell1;
497  temp = new vnl_matrix_inverse<double>((*tempSHBasis));
498  m_Interpolation_SHT1_inv = new vnl_matrix<double>(temp->inverse());
499 
500  delete Q;
501  delete temp;
502  delete tempSHBasis;
503 
504  // Shell 2
505  Q = new vnl_matrix<double>(3, shell2.size());
506  ComputeSphericalFromCartesian(Q, shell2);
507 
508  tempSHBasis = new vnl_matrix<double>(shell2.size(), number_coeffs_shell2);
509  ComputeSphericalHarmonicsBasis(Q, tempSHBasis, interp_SHOrder_shell2);
510 
511  temp = new vnl_matrix_inverse<double>((*tempSHBasis));
512 
513  m_Interpolation_SHT2_inv = new vnl_matrix<double>(temp->inverse());
514 
515  delete Q;
516  delete temp;
517  delete tempSHBasis;
518 
519  // Shell 3
520  Q = new vnl_matrix<double>(3, shell3.size());
521  ComputeSphericalFromCartesian(Q, shell3);
522 
523  tempSHBasis = new vnl_matrix<double>(shell3.size(), number_coeffs_shell3);
524  ComputeSphericalHarmonicsBasis(Q, tempSHBasis, interp_SHOrder_shell3);
525 
526  temp = new vnl_matrix_inverse<double>((*tempSHBasis));
527 
528  m_Interpolation_SHT3_inv = new vnl_matrix<double>(temp->inverse());
529 
530  delete Q;
531  delete temp;
532  delete tempSHBasis;
533 
534  ComputeReconstructionMatrix(all_directions_container);
535 
536  MITK_INFO << "Reconstruction information: Multishell Reconstruction filter - Interpolation";
537  MITK_INFO << "Shell 1";
538  MITK_INFO << " SHOrder: " << interp_SHOrder_shell1;
539  MITK_INFO << " Number of Coeffs: " << number_coeffs_shell1;
540  MITK_INFO << " Number of Gradientdirections: " << size_shell1;
541 
542 
543  MITK_INFO << "Shell 2";
544  MITK_INFO << " SHOrder: " << interp_SHOrder_shell2;
545  MITK_INFO << " Number of Coeffs: " << number_coeffs_shell2;
546  MITK_INFO << " Number of Gradientdirections: " << size_shell2;
547 
548 
549  MITK_INFO << "Shell 3";
550  MITK_INFO << " SHOrder: " << interp_SHOrder_shell3;
551  MITK_INFO << " Number of Coeffs: " << number_coeffs_shell3;
552  MITK_INFO << " Number of Gradientdirections: " << size_shell3;
553 
554 
555  MITK_INFO << "Overall";
556  MITK_INFO << " SHOrder: " << L;
557  MITK_INFO << " Number of Coeffs: " << (L+1)*(L+2)*0.5;
558  MITK_INFO << " Number of Gradientdirections: " << m_MaxDirections;
559 
560  return;
561  }else
562  {
563  ComputeReconstructionMatrix(shell1);
564  }
565  }
566  }
567 
568  if(m_BValueMap.size() > 2 && m_ReconstructionType != Mode_Analytical3Shells)
569  {
570  m_ReconstructionType = Mode_NumericalNShells;
571  }
572  if(m_BValueMap.size() == 2){
573  BValueMapIteraotr it = m_BValueMap.begin();
574  it++; // skip b0 entry
575  IndiciesVector shell = it->second;
576  ComputeReconstructionMatrix(shell);
577  }
578 
579 }
580 
581 
582 template< class T, class TG, class TO, int L, int NODF>
583 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
584 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType /*NumberOfThreads*/)
585 {
586 
587  itk::TimeProbe clock;
588  clock.Start();
589  switch(m_ReconstructionType)
590  {
591  case Mode_Standard1Shell:
592  StandardOneShellReconstruction(outputRegionForThread);
593  break;
594  case Mode_Analytical3Shells:
595  AnalyticalThreeShellReconstruction(outputRegionForThread);
596  break;
597  case Mode_NumericalNShells:
598  break;
599  }
600  clock.Stop();
601  MITK_INFO << "Reconstruction in : " << clock.GetTotal() << " s";
602 }
603 
604 
605 template< class T, class TG, class TO, int L, int NODF>
606 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
607 ::StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread)
608 {
609  // Get output image pointer
610  typename OdfImageType::Pointer outputImage = static_cast< OdfImageType * >(ProcessObject::GetPrimaryOutput());
611  // Get input gradient image pointer
612  typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( ProcessObject::GetInput(0) );
613 
614  // ImageRegionIterator for the output image
615  ImageRegionIterator< OdfImageType > oit(outputImage, outputRegionForThread);
616  oit.GoToBegin();
617 
618  // ImageRegionIterator for the BZero (output) image
619  ImageRegionIterator< BZeroImageType > bzeroIterator(m_BZeroImage, outputRegionForThread);
620  bzeroIterator.GoToBegin();
621 
622  // Const ImageRegionIterator for input gradient image
623  typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
624  GradientIteratorType git(gradientImagePointer, outputRegionForThread );
625  git.GoToBegin();
626 
627  BValueMapIteraotr it = m_BValueMap.begin();
628  it++; // skip b0 entry
629  IndiciesVector SignalIndicies = it->second;
630  IndiciesVector BZeroIndicies = m_BValueMap[0];
631 
632  unsigned int NumbersOfGradientIndicies = SignalIndicies.size();
633 
634  typedef typename GradientImagesType::PixelType GradientVectorType;
635 
636  // iterate overall voxels of the gradient image region
637  while( ! git.IsAtEnd() )
638  {
639  GradientVectorType b = git.Get();
640  // ODF Vector
641  OdfPixelType odf(0.0);
642 
643  double b0average = 0;
644  const unsigned int b0size = BZeroIndicies.size();
645  for(unsigned int i = 0; i < b0size ; ++i)
646  {
647  b0average += b[BZeroIndicies[i]];
648  }
649  b0average /= b0size;
650  bzeroIterator.Set(b0average);
651  ++bzeroIterator;
652 
653  // Create the Signal Vector
654  vnl_vector<double> SignalVector(NumbersOfGradientIndicies);
655  if( (b0average != 0) && (b0average >= m_Threshold) )
656  {
657 
658  for( unsigned int i = 0; i< SignalIndicies.size(); i++ )
659  {
660  SignalVector[i] = static_cast<double>(b[SignalIndicies[i]]);
661  }
662 
663  // apply threashold an generate ln(-ln(E)) signal
664  // Replace SignalVector with PreNormalized SignalVector
665  S_S0Normalization(SignalVector, b0average);
666  Projection1(SignalVector);
667 
668  DoubleLogarithm(SignalVector);
669 
670  // approximate ODF coeffs
671  vnl_vector<double> coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector );
672  coeffs[0] = 1.0/(2.0*sqrt(M_PI));
673 
674  odf = element_cast<double, TO>(( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs )).data_block();
675  odf *= (M_PI*4/NODF);
676  }
677  // set ODF to ODF-Image
678  oit.Set( odf );
679  ++oit;
680  ++git;
681 
682  }
683 
684  MITK_INFO << "One Thread finished reconstruction";
685 }
686 
687 //#include "itkLevenbergMarquardtOptimizer.h"
688 
689 template< class T, class TG, class TO, int L, int NODF>
690 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
691 ::NumericalNShellReconstruction(const OutputImageRegionType& outputRegionForThread)
692 {
693 
694  /* itk::LevenbergMarquardtOptimizer::Pointer optimizer = itk::LevenbergMarquardtOptimizer::New();
695  optimizer->SetUseCostFunctionGradient(false);
696 
697  // Scale the translation components of the Transform in the Optimizer
698  itk::LevenbergMarquardtOptimizer::ScalesType scales(transform->GetNumberOfParameters());
699  scales.Fill(0.01);
700  unsigned long numberOfIterations = 80000;
701  double gradientTolerance = 1e-10; // convergence criterion
702  double valueTolerance = 1e-10; // convergence criterion
703  double epsilonFunction = 1e-10; // convergence criterion
704  optimizer->SetScales( scales );
705  optimizer->SetNumberOfIterations( numberOfIterations );
706  optimizer->SetValueTolerance( valueTolerance );
707  optimizer->SetGradientTolerance( gradientTolerance );
708  optimizer->SetEpsilonFunction( epsilonFunction );*/
709 
710 
711 }
712 
713 
714 template< class T, class TG, class TO, int L, int NODF>
715 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
716 ::AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread)
717 {
718  // Input Gradient Image and Output ODF Image
719  typedef typename GradientImagesType::PixelType GradientVectorType;
720  typename OdfImageType::Pointer outputImage = static_cast< OdfImageType * >(ProcessObject::GetPrimaryOutput());
721  typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( ProcessObject::GetInput(0) );
722 
723  // Define Image iterators
724  ImageRegionIterator< OdfImageType > odfOutputImageIterator(outputImage, outputRegionForThread);
725  ImageRegionConstIterator< GradientImagesType > gradientInputImageIterator(gradientImagePointer, outputRegionForThread );
726  ImageRegionIterator< BZeroImageType > bzeroIterator(m_BZeroImage, outputRegionForThread);
727  ImageRegionIterator< CoefficientImageType > coefficientImageIterator(m_CoefficientImage, outputRegionForThread);
728 
729  // All iterators seht to Begin of the specific OutputRegion
730  coefficientImageIterator.GoToBegin();
731  bzeroIterator.GoToBegin();
732  odfOutputImageIterator.GoToBegin();
733  gradientInputImageIterator.GoToBegin();
734 
735  // Get Shell Indicies for all non-BZero Gradients
736  // it MUST be a arithmetic progression eg.: 1000, 2000, 3000
737  BValueMapIteraotr it = m_BValueMap.begin();
738  it++;
739  // it = b-value = 1000
740  IndiciesVector Shell1Indiecies = it->second;
741  it++;
742  // it = b-value = 2000
743  IndiciesVector Shell2Indiecies = it->second;
744  it++;
745  // it = b-value = 3000
746  IndiciesVector Shell3Indiecies = it->second;
747  IndiciesVector BZeroIndicies = m_BValueMap[0];
748 
749  if(!m_Interpolation_Flag)
750  {
751  m_MaxDirections = Shell1Indiecies.size();
752  }// else: m_MaxDirection is set in BeforeThreadedGenerateData
753 
754  // Nx3 Signal Matrix with E(0) = Shell 1, E(1) = Shell 2, E(2) = Shell 3
755  vnl_vector< double > E1(m_MaxDirections);
756  vnl_vector< double > E2(m_MaxDirections);
757  vnl_vector< double > E3(m_MaxDirections);
758 
759  vnl_vector<double> AlphaValues(m_MaxDirections);
760  vnl_vector<double> BetaValues(m_MaxDirections);
761  vnl_vector<double> LAValues(m_MaxDirections);
762  vnl_vector<double> PValues(m_MaxDirections);
763 
764  vnl_vector<double> DataShell1(Shell1Indiecies.size());
765  vnl_vector<double> DataShell2(Shell2Indiecies.size());
766  vnl_vector<double> DataShell3(Shell3Indiecies.size());
767 
768  vnl_matrix<double> tempInterpolationMatrixShell1,tempInterpolationMatrixShell2,tempInterpolationMatrixShell3;
769 
770  if(m_Interpolation_Flag)
771  {
772  tempInterpolationMatrixShell1 = (*m_TARGET_SH_shell1) * (*m_Interpolation_SHT1_inv);
773  tempInterpolationMatrixShell2 = (*m_TARGET_SH_shell2) * (*m_Interpolation_SHT2_inv);
774  tempInterpolationMatrixShell3 = (*m_TARGET_SH_shell3) * (*m_Interpolation_SHT3_inv);
775  }
776 
777  OdfPixelType odf(0.0);
778  typename CoefficientImageType::PixelType coeffPixel(0.0);
779 
780  double P2,A,B2,B,P,alpha,beta,lambda, ER1, ER2;
781 
782 
783 
784  // iterate overall voxels of the gradient image region
785  while( ! gradientInputImageIterator.IsAtEnd() )
786  {
787 
788  odf = 0.0;
789  coeffPixel = 0.0;
790 
791  GradientVectorType b = gradientInputImageIterator.Get();
792 
793  // calculate for each shell the corresponding b0-averages
794  double shell1b0Norm =0;
795  double shell2b0Norm =0;
796  double shell3b0Norm =0;
797  double b0average = 0;
798  const unsigned int b0size = BZeroIndicies.size();
799 
800  if(b0size == 1)
801  {
802  shell1b0Norm = b[BZeroIndicies[0]];
803  shell2b0Norm = b[BZeroIndicies[0]];
804  shell3b0Norm = b[BZeroIndicies[0]];
805  b0average = b[BZeroIndicies[0]];
806  }else if(b0size % 3 ==0)
807  {
808  for(unsigned int i = 0; i < b0size ; ++i)
809  {
810  if(i < b0size / 3) shell1b0Norm += b[BZeroIndicies[i]];
811  if(i >= b0size / 3 && i < (b0size / 3)*2) shell2b0Norm += b[BZeroIndicies[i]];
812  if(i >= (b0size / 3) * 2) shell3b0Norm += b[BZeroIndicies[i]];
813  }
814  shell1b0Norm /= (b0size/3);
815  shell2b0Norm /= (b0size/3);
816  shell3b0Norm /= (b0size/3);
817  b0average = (shell1b0Norm + shell2b0Norm+ shell3b0Norm)/3;
818  }else
819  {
820  for(unsigned int i = 0; i <b0size ; ++i)
821  {
822  shell1b0Norm += b[BZeroIndicies[i]];
823  }
824  shell1b0Norm /= b0size;
825  shell2b0Norm = shell1b0Norm;
826  shell3b0Norm = shell1b0Norm;
827  b0average = shell1b0Norm;
828  }
829 
830  bzeroIterator.Set(b0average);
831  ++bzeroIterator;
832 
833  if( (b0average != 0) && ( b0average >= m_Threshold) )
834  {
835  // Get the Signal-Value for each Shell at each direction (specified in the ShellIndicies Vector .. this direction corresponse to this shell...)
836 
837  /*//fsl fix ---------------------------------------------------
838  for(int i = 0 ; i < Shell1Indiecies.size(); i++)
839  DataShell1[i] = static_cast<double>(b[Shell1Indiecies[i]]);
840  for(int i = 0 ; i < Shell2Indiecies.size(); i++)
841  DataShell2[i] = static_cast<double>(b[Shell2Indiecies[i]]);
842  for(int i = 0 ; i < Shell3Indiecies.size(); i++)
843  DataShell3[i] = static_cast<double>(b[Shell2Indiecies[i]]);
844 
845  // Normalize the Signal: Si/S0
846  S_S0Normalization(DataShell1, shell1b0Norm);
847  S_S0Normalization(DataShell2, shell2b0Norm);
848  S_S0Normalization(DataShell3, shell2b0Norm);
849  *///fsl fix -------------------------------------------ende--
850 
852  for(unsigned int i = 0 ; i < Shell1Indiecies.size(); i++)
853  DataShell1[i] = static_cast<double>(b[Shell1Indiecies[i]]);
854  for(unsigned int i = 0 ; i < Shell2Indiecies.size(); i++)
855  DataShell2[i] = static_cast<double>(b[Shell2Indiecies[i]]);
856  for(unsigned int i = 0 ; i < Shell3Indiecies.size(); i++)
857  DataShell3[i] = static_cast<double>(b[Shell3Indiecies[i]]);
858 
859 
860 
861  // Normalize the Signal: Si/S0
862  S_S0Normalization(DataShell1, shell1b0Norm);
863  S_S0Normalization(DataShell2, shell2b0Norm);
864  S_S0Normalization(DataShell3, shell3b0Norm);
865 
866 
867  if(m_Interpolation_Flag)
868  {
869  E1 = tempInterpolationMatrixShell1 * DataShell1;
870  E2 = tempInterpolationMatrixShell2 * DataShell2;
871  E3 = tempInterpolationMatrixShell3 * DataShell3;
872  }else{
873  E1 = (DataShell1);
874  E2 = (DataShell2);
875  E3 = (DataShell3);
876  }
877 
878  //Implements Eq. [19] and Fig. 4.
879  Projection1(E1);
880  Projection1(E2);
881  Projection1(E3);
882  //inqualities [31]. Taking the lograithm of th first tree inqualities
883  //convert the quadratic inqualities to linear ones.
884  Projection2(E1,E2,E3);
885 
886  for( unsigned int i = 0; i< m_MaxDirections; i++ )
887  {
888  double e1 = E1.get(i);
889  double e2 = E2.get(i);
890  double e3 = E3.get(i);
891 
892  P2 = e2-e1*e1;
893  A = (e3 -e1*e2) / ( 2* P2);
894  B2 = A * A -(e1 * e3 - e2 * e2) /P2;
895  B = 0;
896  if(B2 > 0) B = sqrt(B2);
897  P = 0;
898  if(P2 > 0) P = sqrt(P2);
899 
900  alpha = A + B;
901  beta = A - B;
902 
903  PValues.put(i, P);
904  AlphaValues.put(i, alpha);
905  BetaValues.put(i, beta);
906 
907  }
908 
909  Projection3(PValues, AlphaValues, BetaValues);
910 
911  for(unsigned int i = 0 ; i < m_MaxDirections; i++)
912  {
913  const double fac = (PValues[i] * 2 ) / (AlphaValues[i] - BetaValues[i]);
914  lambda = 0.5 + 0.5 * std::sqrt(1 - fac * fac);;
915  ER1 = std::fabs(lambda * (AlphaValues[i] - BetaValues[i]) + (BetaValues[i] - E1.get(i) ))
916  + std::fabs(lambda * (AlphaValues[i] * AlphaValues[i] - BetaValues[i] * BetaValues[i]) + (BetaValues[i] * BetaValues[i] - E2.get(i) ))
917  + std::fabs(lambda * (AlphaValues[i] * AlphaValues[i] * AlphaValues[i] - BetaValues[i] * BetaValues[i] * BetaValues[i]) + (BetaValues[i] * BetaValues[i] * BetaValues[i] - E3.get(i) ));
918  ER2 = std::fabs((1-lambda) * (AlphaValues[i] - BetaValues[i]) + (BetaValues[i] - E1.get(i) ))
919  + std::fabs((1-lambda) * (AlphaValues[i] * AlphaValues[i] - BetaValues[i] * BetaValues[i]) + (BetaValues[i] * BetaValues[i] - E2.get(i) ))
920  + std::fabs((1-lambda) * (AlphaValues[i] * AlphaValues[i] * AlphaValues[i] - BetaValues[i] * BetaValues[i] * BetaValues[i]) + (BetaValues[i] * BetaValues[i] * BetaValues[i] - E3.get(i)));
921  if(ER1 < ER2)
922  LAValues.put(i, lambda);
923  else
924  LAValues.put(i, 1-lambda);
925 
926  }
927 
928  DoubleLogarithm(AlphaValues);
929  DoubleLogarithm(BetaValues);
930 
931  vnl_vector<double> SignalVector(element_product((LAValues) , (AlphaValues)-(BetaValues)) + (BetaValues));
932 
933  vnl_vector<double> coeffs((*m_CoeffReconstructionMatrix) *SignalVector );
934 
935  // the first coeff is a fix value
936  coeffs[0] = 1.0/(2.0*sqrt(M_PI));
937  coeffPixel = element_cast<double, TO>(coeffs).data_block();
938 
939  // Cast the Signal-Type from double to float for the ODF-Image
940  odf = element_cast<double, TO>( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs ).data_block();
941  odf *= ((M_PI*4)/NODF);
942  }
943 
944  // set ODF to ODF-Image
945  coefficientImageIterator.Set(coeffPixel);
946  odfOutputImageIterator.Set( odf );
947  ++odfOutputImageIterator;
948  ++coefficientImageIterator;
949  ++gradientInputImageIterator;
950  }
951 
952 }
953 
954 
955 
956 template< class T, class TG, class TO, int L, int NODF>
957 void DiffusionMultiShellQballReconstructionImageFilter<T, TG, TO, L, NODF>::
958 ComputeSphericalHarmonicsBasis(vnl_matrix<double> * QBallReference, vnl_matrix<double> *SHBasisOutput, int LOrder , vnl_matrix<double>* LaplaciaBaltramiOutput, vnl_vector<int>* SHOrderAssociation, vnl_matrix<double>* SHEigenvalues)
959 {
960 
961  // MITK_INFO << *QBallReference;
962 
963  for(unsigned int i=0; i< (*SHBasisOutput).rows(); i++)
964  {
965  for(int k = 0; k <= LOrder; k += 2)
966  {
967  for(int m =- k; m <= k; m++)
968  {
969  int j = ( k * k + k + 2 ) / 2 + m - 1;
970 
971  // Compute SHBasisFunctions
972  if(QBallReference){
973  double phi = (*QBallReference)(0,i);
974  double th = (*QBallReference)(1,i);
975  (*SHBasisOutput)(i,j) = mitk::sh::Yj(m,k,th,phi);
976  }
977 
978  // Laplacian Baltrami Order Association
979  if(LaplaciaBaltramiOutput)
980  (*LaplaciaBaltramiOutput)(j,j) = k*k*(k + 1)*(k+1);
981 
982  // SHEigenvalues with order Accosiation kj
983  if(SHEigenvalues)
984  (*SHEigenvalues)(j,j) = -k* (k+1);
985 
986  // Order Association
987  if(SHOrderAssociation)
988  (*SHOrderAssociation)[j] = k;
989 
990  }
991  }
992  }
993 }
994 
995 template< class T, class TG, class TO, int L, int NOdfDirections>
996 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NOdfDirections>
997 ::ComputeReconstructionMatrix(IndiciesVector const & refVector)
998 {
999 
1000  typedef std::unique_ptr< vnl_matrix< double> > MatrixDoublePtr;
1001  typedef std::unique_ptr< vnl_vector< int > > VectorIntPtr;
1002  typedef std::unique_ptr< vnl_matrix_inverse< double > > InverseMatrixDoublePtr;
1003 
1004  int numberOfGradientDirections = refVector.size();
1005 
1006  if( numberOfGradientDirections <= (((L+1)*(L+2))/2) || numberOfGradientDirections < 6 )
1007  {
1008  itkExceptionMacro( << "At least (L+1)(L+2)/2 gradient directions for each shell are required; current : " << numberOfGradientDirections );
1009  }
1010 
1011  CheckDuplicateDiffusionGradients();
1012 
1013  const int LOrder = L;
1014  int NumberOfCoeffs = (int)(LOrder*LOrder + LOrder + 2.0)/2.0 + LOrder;
1015  MatrixDoublePtr SHBasisMatrix(new vnl_matrix<double>(numberOfGradientDirections,NumberOfCoeffs));
1016  SHBasisMatrix->fill(0.0);
1017  VectorIntPtr SHOrderAssociation(new vnl_vector<int>(NumberOfCoeffs));
1018  SHOrderAssociation->fill(0.0);
1019  MatrixDoublePtr LaplacianBaltrami(new vnl_matrix<double>(NumberOfCoeffs,NumberOfCoeffs));
1020  LaplacianBaltrami->fill(0.0);
1021  MatrixDoublePtr FRTMatrix(new vnl_matrix<double>(NumberOfCoeffs,NumberOfCoeffs));
1022  FRTMatrix->fill(0.0);
1023  MatrixDoublePtr SHEigenvalues(new vnl_matrix<double>(NumberOfCoeffs,NumberOfCoeffs));
1024  SHEigenvalues->fill(0.0);
1025 
1026  MatrixDoublePtr Q(new vnl_matrix<double>(3, numberOfGradientDirections));
1027 
1028  // Convert Cartesian to Spherical Coordinates refVector -> Q
1029  ComputeSphericalFromCartesian(Q.get(), refVector);
1030 
1031  // SHBasis-Matrix + LaplacianBaltrami-Matrix + SHOrderAssociationVector
1032  ComputeSphericalHarmonicsBasis(Q.get() ,SHBasisMatrix.get() , LOrder , LaplacianBaltrami.get(), SHOrderAssociation.get(), SHEigenvalues.get());
1033 
1034  // Compute FunkRadon Transformation Matrix Associated to SHBasis Order lj
1035  for(int i=0; i<NumberOfCoeffs; i++)
1036  {
1037  (*FRTMatrix)(i,i) = 2.0 * M_PI * mitk::sh::legendre0((*SHOrderAssociation)[i]);
1038  }
1039 
1040  MatrixDoublePtr temp(new vnl_matrix<double>(((SHBasisMatrix->transpose()) * (*SHBasisMatrix)) + (m_Lambda * (*LaplacianBaltrami))));
1041 
1042  InverseMatrixDoublePtr pseudo_inv(new vnl_matrix_inverse<double>((*temp)));
1043  MatrixDoublePtr inverse(new vnl_matrix<double>(NumberOfCoeffs,NumberOfCoeffs));
1044  (*inverse) = pseudo_inv->inverse();
1045 
1046  const double factor = (1.0/(16.0*M_PI*M_PI));
1047  MatrixDoublePtr SignalReonstructionMatrix (new vnl_matrix<double>((*inverse) * (SHBasisMatrix->transpose())));
1048  m_CoeffReconstructionMatrix = new vnl_matrix<double>(( factor * ((*FRTMatrix) * ((*SHEigenvalues) * (*SignalReonstructionMatrix))) ));
1049 
1050 
1051  // SH Basis for ODF-reconstruction
1052  vnl_matrix_fixed<double, 3, NOdfDirections>* U = PointShell<NOdfDirections, vnl_matrix_fixed<double, 3, NOdfDirections> >::DistributePointShell();
1053  for(int i=0; i<NOdfDirections; i++)
1054  {
1055  double x = (*U)(0,i);
1056  double y = (*U)(1,i);
1057  double z = (*U)(2,i);
1058  double cart[3];
1059  mitk::sh::Cart2Sph(x,y,z,cart);
1060  (*U)(0,i) = cart[0];
1061  (*U)(1,i) = cart[1];
1062  (*U)(2,i) = cart[2];
1063  }
1064 
1065  MatrixDoublePtr tempPtr (new vnl_matrix<double>( U->as_matrix() ));
1066  m_ODFSphericalHarmonicBasisMatrix = new vnl_matrix<double>(NOdfDirections,NumberOfCoeffs);
1067  ComputeSphericalHarmonicsBasis(tempPtr.get(), m_ODFSphericalHarmonicBasisMatrix, LOrder);
1068 }
1069 
1070 template< class T, class TG, class TO, int L, int NOdfDirections>
1071 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NOdfDirections>
1072 ::ComputeSphericalFromCartesian(vnl_matrix<double> * Q, IndiciesVector const & refShell)
1073 {
1074  for(unsigned int i = 0; i < refShell.size(); i++)
1075  {
1076  double x = m_GradientDirectionContainer->ElementAt(refShell[i]).normalize().get(0);
1077  double y = m_GradientDirectionContainer->ElementAt(refShell[i]).normalize().get(1);
1078  double z = m_GradientDirectionContainer->ElementAt(refShell[i]).normalize().get(2);
1079  double cart[3];
1080  mitk::sh::Cart2Sph(x,y,z,cart);
1081  (*Q)(0,i) = cart[0];
1082  (*Q)(1,i) = cart[1];
1083  (*Q)(2,i) = cart[2];
1084  }
1085 }
1086 
1087 
1088 template< class T, class TG, class TO, int L, int NODF>
1089 bool DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
1090 ::CheckDuplicateDiffusionGradients()
1091 {
1092  bool value = false;
1093 
1094  BValueMapIteraotr mapIterator = m_BValueMap.begin();
1095  mapIterator++;
1096  while(mapIterator != m_BValueMap.end())
1097  {
1098  std::vector<unsigned int>::const_iterator it1 = mapIterator->second.begin();
1099  std::vector<unsigned int>::const_iterator it2 = mapIterator->second.begin();
1100 
1101  for(; it1 != mapIterator->second.end(); ++it1)
1102  {
1103  for(; it2 != mapIterator->second.end(); ++it2)
1104  {
1105  if(m_GradientDirectionContainer->ElementAt(*it1) == m_GradientDirectionContainer->ElementAt(*it2) && it1 != it2)
1106  {
1107  itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." );
1108  value = true;
1109  }
1110  }
1111  }
1112  ++mapIterator;
1113  }
1114  return value;
1115 }
1116 
1117 template< class T, class TG, class TO, int L, int NODF>
1118 std::vector<unsigned int> DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
1119 ::GetAllDirections()
1120 {
1121  IndiciesVector directioncontainer;
1122 
1123  BValueMapIteraotr mapIterator = m_BValueMap.begin();
1124  mapIterator++;
1125  IndiciesVector shell1 = mapIterator->second;
1126  mapIterator++;
1127  IndiciesVector shell2 = mapIterator->second;
1128  mapIterator++;
1129  IndiciesVector shell3 = mapIterator->second;
1130 
1131  while(shell1.size()>0)
1132  {
1133  unsigned int wntIndex = shell1.back();
1134  shell1.pop_back();
1135 
1136  IndiciesVector::iterator containerIt = directioncontainer.begin();
1137  bool directionExist = false;
1138  while(containerIt != directioncontainer.end())
1139  {
1140  if (fabs(dot(m_GradientDirectionContainer->ElementAt(*containerIt), m_GradientDirectionContainer->ElementAt(wntIndex))) > 0.9998)
1141  {
1142  directionExist = true;
1143  break;
1144  }
1145  containerIt++;
1146  }
1147  if(!directionExist)
1148  {
1149  directioncontainer.push_back(wntIndex);
1150  }
1151  }
1152 
1153  while(shell2.size()>0)
1154  {
1155  unsigned int wntIndex = shell2.back();
1156  shell2.pop_back();
1157 
1158  IndiciesVector::iterator containerIt = directioncontainer.begin();
1159  bool directionExist = false;
1160  while(containerIt != directioncontainer.end())
1161  {
1162  if (fabs(dot(m_GradientDirectionContainer->ElementAt(*containerIt), m_GradientDirectionContainer->ElementAt(wntIndex))) > 0.9998)
1163  {
1164  directionExist = true;
1165  break;
1166  }
1167  containerIt++;
1168  }
1169  if(!directionExist)
1170  {
1171  directioncontainer.push_back(wntIndex);
1172  }
1173  }
1174 
1175  while(shell3.size()>0)
1176  {
1177  unsigned int wntIndex = shell3.back();
1178  shell3.pop_back();
1179 
1180  IndiciesVector::iterator containerIt = directioncontainer.begin();
1181  bool directionExist = false;
1182  while(containerIt != directioncontainer.end())
1183  {
1184  if (fabs(dot(m_GradientDirectionContainer->ElementAt(*containerIt), m_GradientDirectionContainer->ElementAt(wntIndex))) > 0.9998)
1185  {
1186  directionExist = true;
1187  break;
1188  }
1189  containerIt++;
1190  }
1191  if(!directionExist)
1192  {
1193  directioncontainer.push_back(wntIndex);
1194  }
1195  }
1196 
1197  return directioncontainer;
1198 }
1199 
1200 // corresponding directions between shells (e.g. dir1_shell1 vs dir1_shell2) differ more than 1 degree.
1201 template< class T, class TG, class TO, int L, int NODF>
1202 bool DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
1203 ::CheckForDifferingShellDirections()
1204 {
1205  bool interp_flag = false;
1206 
1207  BValueMapIteraotr mapIterator = m_BValueMap.begin();
1208  mapIterator++;
1209  IndiciesVector shell1 = mapIterator->second;
1210  mapIterator++;
1211  IndiciesVector shell2 = mapIterator->second;
1212  mapIterator++;
1213  IndiciesVector shell3 = mapIterator->second;
1214 
1215  for (unsigned int i=0; i< shell1.size(); i++)
1216  if (fabs(dot(m_GradientDirectionContainer->ElementAt(shell1[i]), m_GradientDirectionContainer->ElementAt(shell2[i]))) <= 0.9998) {interp_flag=true; break;}
1217  for (unsigned int i=0; i< shell1.size(); i++)
1218  if (fabs(dot(m_GradientDirectionContainer->ElementAt(shell1[i]), m_GradientDirectionContainer->ElementAt(shell3[i]))) <= 0.9998) {interp_flag=true; break;}
1219  for (unsigned int i=0; i< shell1.size(); i++)
1220  if (fabs(dot(m_GradientDirectionContainer->ElementAt(shell2[i]), m_GradientDirectionContainer->ElementAt(shell3[i]))) <= 0.9998) {interp_flag=true; break;}
1221  return interp_flag;
1222 }
1223 
1224 
1225 template< class T, class TG, class TO, int L, int NODF>
1226 void DiffusionMultiShellQballReconstructionImageFilter<T,TG,TO,L,NODF>
1227 ::PrintSelf(std::ostream& os, Indent indent) const
1228 {
1229  std::locale C("C");
1230  std::locale originalLocale = os.getloc();
1231  os.imbue(C);
1232 
1233  Superclass::PrintSelf(os,indent);
1234 
1235  //os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl;
1236  if ( m_GradientDirectionContainer )
1237  {
1238  os << indent << "GradientDirectionContainer: "
1239  << m_GradientDirectionContainer << std::endl;
1240  }
1241  else
1242  {
1243  os << indent <<
1244  "GradientDirectionContainer: (Gradient directions not set)" << std::endl;
1245  }
1246  os << indent << "NumberOfGradientDirections: " <<
1247  m_NumberOfGradientDirections << std::endl;
1248  os << indent << "NumberOfBaselineImages: " <<
1249  m_NumberOfBaselineImages << std::endl;
1250  os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl;
1251  os << indent << "BValue: " << m_BValue << std::endl;
1252 
1253  os.imbue( originalLocale );
1254 }
1255 
1256 
1257 
1258 
1259 }
1260 
1261 #endif // __itkDiffusionMultiShellQballReconstructionImageFilter_cpp
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
MITK_INFO
#define MITK_INFO
Definition: mitkLogMacros.h:22
mitk::sh::legendre0
static double legendre0(int l)
Definition: mitkDiffusionFunctionCollection.cpp:65
mitk::sh::Yj
static double Yj(int m, int k, double theta, double phi)
Definition: mitkDiffusionFunctionCollection.cpp:82
mitk::sh::Cart2Sph
static void Cart2Sph(double x, double y, double z, double *cart)
Definition: mitkDiffusionFunctionCollection.cpp:46
itk::DiffusionMultiShellQballReconstructionImageFilter::SetGradientImage
void SetGradientImage(const GradientDirectionContainerType *gradientDirectionContainer, const GradientImagesType *gradientImage, float bvalue)
Definition: itkDiffusionMultiShellQballReconstructionImageFilter.cpp:60
PixelType
unsigned short PixelType
Definition: TumorProgressionMapping.cpp:55
itk::DiffusionMultiShellQballReconstructionImageFilter::GradientDirectionContainerType
VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > GradientDirectionContainerType
Definition: itkDiffusionMultiShellQballReconstructionImageFilter.h:59
itk::DiffusionMultiShellQballReconstructionImageFilter::ThreadedGenerateData
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, ThreadIdType NumberOfThreads)
Definition: itkDiffusionMultiShellQballReconstructionImageFilter.cpp:584
itk::DiffusionMultiShellQballReconstructionImageFilter::BValueMapIteraotr
std::map< unsigned int, std::vector< unsigned int > >::iterator BValueMapIteraotr
Definition: itkDiffusionMultiShellQballReconstructionImageFilter.h:64
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.