Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkDiffusionQballReconstructionImageFilter.txx
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 __itkDiffusionQballReconstructionImageFilter_txx
17 #define __itkDiffusionQballReconstructionImageFilter_txx
18 
19 #include "itkDiffusionQballReconstructionImageFilter.h"
20 #include "itkImageRegionConstIterator.h"
21 #include "itkImageRegionConstIteratorWithIndex.h"
22 #include "itkImageRegionIterator.h"
23 #include "itkArray.h"
24 #include "vnl/vnl_vector.h"
25 #include "itkPointShell.h"
26 
27 #define _USE_MATH_DEFINES
28 #include <math.h>
29 
30 namespace itk {
31 
32 #define QBALL_RECON_PI M_PI
33 
34  template< class TReferenceImagePixelType,
35  class TGradientImagePixelType,
36  class TOdfPixelType,
37  int NrOdfDirections,
38  int NrBasisFunctionCenters>
39  DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
40  TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
41  NrBasisFunctionCenters>
42  ::DiffusionQballReconstructionImageFilter() :
43  m_GradientDirectionContainer(NULL),
44  m_NumberOfGradientDirections(0),
45  m_NumberOfEquatorSamplingPoints(0),
46  m_NumberOfBaselineImages(1),
47  m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()),
48  m_BValue(1.0),
49  m_GradientImageTypeEnumeration(Else),
50  m_DirectionsDuplicated(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  template< class TReferenceImagePixelType,
58  class TGradientImagePixelType,
59  class TOdfPixelType,
60  int NrOdfDirections,
61  int NrBasisFunctionCenters>
62  void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
63  TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
64  NrBasisFunctionCenters>
65  ::BeforeThreadedGenerateData()
66  {
67  // If we have more than 2 inputs, then each input, except the first is a
68  // gradient image. The number of gradient images must match the number of
69  // gradient directions.
70  const unsigned int numberOfInputs = this->GetNumberOfInputs();
71 
72  // There need to be at least 6 gradient directions to be able to compute the
73  // tensor basis
74  if( m_NumberOfGradientDirections < 1 )
75  {
76  itkExceptionMacro( << "Your image contains no diffusion gradients!" );
77  }
78 
79  // If there is only 1 gradient image, it must be an itk::VectorImage. Otherwise
80  // we must have a container of (numberOfInputs-1) itk::Image. Check to make sure
81  if ( numberOfInputs == 1
82  && m_GradientImageTypeEnumeration != GradientIsInASingleImage )
83  {
84  std::string gradientImageClassName(
85  this->ProcessObject::GetInput(0)->GetNameOfClass());
86  if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 )
87  {
88  itkExceptionMacro( <<
89  "There is only one Gradient image. I expect that to be a VectorImage. "
90  << "But its of type: " << gradientImageClassName );
91  }
92  }
93 
94  // Compute reconstruction matrix that is multiplied to the data-vector
95  // each voxel in order to reconstruct the ODFs
96  this->ComputeReconstructionMatrix();
97 
98  // Allocate the b-zero image
99  m_BZeroImage = BZeroImageType::New();
100  if( m_GradientImageTypeEnumeration == GradientIsInManyImages )
101  {
102  typename ReferenceImageType::Pointer img = static_cast< ReferenceImageType * >(this->ProcessObject::GetInput(0));
103  m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing
104  m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin
105  m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction
106  m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
107  m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
108  }
109  // The gradients are specified in a single multi-component image
110  else if( m_GradientImageTypeEnumeration == GradientIsInASingleImage )
111  {
112  typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >(
113  this->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  }
120 
121  m_BZeroImage->Allocate();
122 
123  }
124 
125  template< class TReferenceImagePixelType,
126  class TGradientImagePixelType,
127  class TOdfPixelType,
128  int NrOdfDirections,
129  int NrBasisFunctionCenters>
130  typename itk::DiffusionQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters>::OdfPixelType
131  itk::DiffusionQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters>::Normalize( OdfPixelType odf, typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
132  {
133  switch( m_NormalizationMethod )
134  {
135 
136  // divide by sum to retreive a PDF
137  case QBR_STANDARD:
138  {
139  odf.Normalize();
140  return odf;
141  break;
142  }
143  // ADC style
144  case QBR_B_ZERO_B_VALUE:
145  {
146  for(int i=0; i<NrOdfDirections; i++)
147  {
148  odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
149  }
150  return odf;
151  break;
152  }
153  // divide by b-zero value of voxel
154  case QBR_B_ZERO:
155  {
156  odf *= 1.0/b0;
157  return odf;
158  break;
159  }
160  // no normalization of ODF
161  case QBR_NONE:
162  {
163  return odf;
164  break;
165  }
166  }
167 
168  return odf;
169  }
170 
171 
172  template< class TReferenceImagePixelType,
173  class TGradientImagePixelType,
174  class TOdfPixelType,
175  int NrOdfDirections,
176  int NrBasisFunctionCenters >
177  vnl_vector<TOdfPixelType> itk::DiffusionQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters>::PreNormalize( vnl_vector<TOdfPixelType> vec )
178  {
179  switch( m_NormalizationMethod )
180  {
181  // standard: no normalization before reconstruction
182  case QBR_STANDARD:
183  {
184  return vec;
185  break;
186  }
187  // log of signal
188  case QBR_B_ZERO_B_VALUE:
189  {
190  int n = vec.size();
191  for(int i=0; i<n; i++)
192  {
193  vec[i] = log(vec[i]);
194  }
195  return vec;
196  break;
197  }
198  // no normalization before reconstruction here
199  case QBR_B_ZERO:
200  {
201  return vec;
202  break;
203  }
204  // no normalization before reconstruction here
205  case QBR_NONE:
206  {
207  return vec;
208  break;
209  }
210  }
211 
212  return vec;
213  }
214 
215  template< class TReferenceImagePixelType,
216  class TGradientImagePixelType,
217  class TOdfPixelType,
218  int NrOdfDirections,
219  int NrBasisFunctionCenters>
220  void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
221  TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
222  NrBasisFunctionCenters>
223  ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
224  ThreadIdType )
225  {
226  // init output and b-zero iterators
227  typename OutputImageType::Pointer outputImage =
228  static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
229  ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
230  oit.GoToBegin();
231  ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
232  oit2.GoToBegin();
233  vnl_vector<TOdfPixelType> B(m_NumberOfGradientDirections);
234 
235  // Two cases here:
236  // 1. Gradients specified in multiple images
237  // 'n' iterators for each of the gradient images
238  // 2. Gradients specified in a single multi-component image
239  // one iterator for all gradient directions
240 
241  if( m_GradientImageTypeEnumeration == GradientIsInManyImages )
242  {
243  // b-zero reference image iterator
244  ImageRegionConstIterator< ReferenceImageType >
245  it(static_cast< ReferenceImageType * >(this->ProcessObject::GetInput(0)),
246  outputRegionForThread);
247  it.GoToBegin();
248 
249  // fill vector with gradient iterators
250  typedef ImageRegionConstIterator< GradientImageType > GradientIteratorType;
251  std::vector< GradientIteratorType * > gradientItContainer;
252  for( unsigned int i = 1; i<= m_NumberOfGradientDirections; i++ )
253  {
254  // adapt index in case we have a duplicated shell
255  int index = i;
256  if(m_DirectionsDuplicated)
257  index = i % (m_NumberOfGradientDirections/2);
258 
259  // init and pushback current input image iterator
260  typename GradientImageType::Pointer gradientImagePointer = NULL;
261  // dynamic_cast would be nice, static because of SGI
262  gradientImagePointer = static_cast< GradientImageType * >(
263  this->ProcessObject::GetInput(index) );
264  GradientIteratorType *git = new GradientIteratorType(
265  gradientImagePointer, outputRegionForThread );
266  git->GoToBegin();
267  gradientItContainer.push_back(git);
268  }
269 
270  // Following loop does the actual reconstruction work in each voxel
271  // (Tuch, Q-Ball Reconstruction [1])
272  while( !it.IsAtEnd() )
273  {
274 
275  // b-zero reference value
276  ReferencePixelType b0 = it.Get();
277 
278  // init ODF
279  OdfPixelType odf(0.0);
280 
281  // threshold on reference value to suppress noisy regions
282  if( (b0 != 0) && (b0 >= m_Threshold) )
283  {
284 
285  // fill array of diffusion measurements
286  for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
287  {
288  GradientPixelType b = gradientItContainer[i]->Get();
289  B[i] = static_cast<TOdfPixelType>(b);
290  ++(*gradientItContainer[i]);
291  }
292 
293  // pre-normalization according to m_NormalizationMethod
294  B = PreNormalize(B);
295 
296  // actual reconstruction
297  odf = ( (*m_ReconstructionMatrix) * B ).data_block();
298 
299  // post-normalization according to m_NormalizationMethod
300  odf.Normalize();
301  }
302  else
303  {
304  // in case we fall below threshold, we just increment to next voxel
305  for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
306  {
307  ++(*gradientItContainer[i]);
308  }
309  }
310 
311  for (unsigned int i=0; i<odf.Size(); i++)
312  if (odf.GetElement(i)!=odf.GetElement(i))
313  odf.Fill(0.0);
314 
315  // set and increment output iterators
316  oit.Set( odf );
317  ++oit;
318  oit2.Set( b0 );
319  ++oit2;
320  ++it;
321  }
322 
323  // clean up
324  for( unsigned int i = 0; i< gradientItContainer.size(); i++ )
325  {
326  delete gradientItContainer[i];
327  }
328  }
329  // The gradients are specified in a single multi-component image
330  else if( m_GradientImageTypeEnumeration == GradientIsInASingleImage )
331  {
332 
333  // init input iterator
334  typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
335  typedef typename GradientImagesType::PixelType GradientVectorType;
336  typename GradientImagesType::Pointer gradientImagePointer = NULL;
337  // dynamic_cast would be nice, static because of SGI
338  gradientImagePointer = static_cast< GradientImagesType * >(
339  this->ProcessObject::GetInput(0) );
340  GradientIteratorType git(gradientImagePointer, outputRegionForThread );
341  git.GoToBegin();
342 
343  // set of indicies each for the baseline images and gradient images
344  std::vector<unsigned int> baselineind; // contains baseline indicies
345  std::vector<unsigned int> gradientind; // contains gradient indicies
346  for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
347  gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
348  {
349  if(gdcit.Value().one_norm() <= 0.0)
350  {
351  baselineind.push_back(gdcit.Index());
352  }
353  else
354  {
355  gradientind.push_back(gdcit.Index());
356  }
357  }
358 
359  // in case we have a duplicated shell, we also duplicate or indices
360  if( m_DirectionsDuplicated )
361  {
362  int gradIndSize = gradientind.size();
363  for(int i=0; i<gradIndSize; i++)
364  gradientind.push_back(gradientind[i]);
365  }
366 
367  // Following loop does the actual reconstruction work in each voxel
368  // (Tuch, Q-Ball Reconstruction [1])
369  while( !git.IsAtEnd() )
370  {
371  // current vector of diffusion measurements
372  GradientVectorType b = git.Get();
373 
374  // average of current b-zero reference values
375  typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
376  for(unsigned int i = 0; i < baselineind.size(); ++i)
377  {
378  b0 += b[baselineind[i]];
379  }
380  b0 /= this->m_NumberOfBaselineImages;
381 
382  // init resulting ODF
383  OdfPixelType odf(0.0);
384 
385  // threshold on reference value to suppress noisy regions
386  if( (b0 != 0) && (b0 >= m_Threshold) )
387  {
388  for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
389  {
390  B[i] = static_cast<TOdfPixelType>(b[gradientind[i]]);
391  }
392 
393  // pre-normalization according to m_NormalizationMethod
394  B = PreNormalize(B);
395 
396  // actual reconstruction
397  odf = ( (*m_ReconstructionMatrix) * B ).data_block();
398 
399  // post-normalization according to m_NormalizationMethod
400  odf = Normalize(odf, b0);
401 
402  }
403 
404  for (unsigned int i=0; i<odf.Size(); i++)
405  if (odf.GetElement(i)!=odf.GetElement(i))
406  odf.Fill(0.0);
407 
408  // set and increment output iterators
409  oit.Set( odf );
410  ++oit;
411  oit2.Set( b0 );
412  ++oit2;
413  ++git; // Gradient image iterator
414  }
415  }
416 
417  std::cout << "One Thread finished reconstruction" << std::endl;
418  }
419 
420  template< class TReferenceImagePixelType,
421  class TGradientImagePixelType,
422  class TOdfPixelType,
423  int NrOdfDirections,
424  int NrBasisFunctionCenters>
425  void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
426  TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
427  NrBasisFunctionCenters>
428  ::ComputeReconstructionMatrix()
429  {
430 
431  if( m_NumberOfGradientDirections < 1 )
432  {
433  itkExceptionMacro( << "Your image contains no diffusion gradients!" );
434  }
435 
436  {
437  // check for duplicate diffusion gradients
438  bool warning = false;
439  for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
440  gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
441  {
442  for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin();
443  gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2)
444  {
445  if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index())
446  {
447  itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." );
448  warning = true;
449  break;
450  }
451  }
452  if (warning) break;
453  }
454 
455  // handle acquisition schemes where only half of the spherical
456  // shell is sampled by the gradient directions. In this case,
457  // each gradient direction is duplicated in negative direction.
458  vnl_vector<double> centerMass(3);
459  centerMass.fill(0.0);
460  int count = 0;
461  for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
462  gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
463  {
464  if(gdcit1.Value().one_norm() > 0.0)
465  {
466  centerMass += gdcit1.Value();
467  count ++;
468  }
469  }
470  centerMass /= count;
471  if(centerMass.two_norm() > 0.1)
472  {
473  m_DirectionsDuplicated = true;
474  m_NumberOfGradientDirections *= 2;
475  }
476  }
477 
478  // set default number of equator sampling points if needed
479  if(!this->m_NumberOfEquatorSamplingPoints)
480  this->m_NumberOfEquatorSamplingPoints
481  = (int) ceil((double)sqrt(8*QBALL_RECON_PI*this->m_NumberOfGradientDirections));
482 
483  vnl_matrix<double>* Q =
484  new vnl_matrix<double>(3, m_NumberOfGradientDirections);
485 
486  {
487  // Fill matrix Q with gradient directions
488  int i = 0;
489  for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
490  gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
491  {
492  if(gdcit.Value().one_norm() > 0.0)
493  {
494  (*Q)(0,i) = gdcit.Value().get(0);
495  (*Q)(1,i) = gdcit.Value().get(1);
496  (*Q)(2,i++) = gdcit.Value().get(2);
497  }
498  }
499  if(m_DirectionsDuplicated)
500  {
501  for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
502  gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
503  {
504  if(gdcit.Value().one_norm() > 0.0)
505  {
506  (*Q)(0,i) = -gdcit.Value().get(0);
507  (*Q)(1,i) = -gdcit.Value().get(1);
508  (*Q)(2,i++) = -gdcit.Value().get(2);
509  }
510  }
511  }
512 
513  }
514 
515  vnl_matrix_fixed<double, 3, NOdfDirections>* U =
516  itk::PointShell<NOdfDirections, vnl_matrix_fixed<double, 3, NOdfDirections> >::DistributePointShell();
517 
518  vnl_matrix_fixed<double, 3, NBasisFunctionCenters>* V =
519  itk::PointShell<NBasisFunctionCenters, vnl_matrix_fixed<double, 3, NBasisFunctionCenters> >::DistributePointShell();
520 
521  // calculate sampling points on the equator perpendicular to z-axis
522  vnl_matrix<double> *C
523  = new vnl_matrix<double>(3, m_NumberOfEquatorSamplingPoints);
524  for(unsigned int i=0; i<m_NumberOfEquatorSamplingPoints; i++)
525  {
526  double theta = i * (2*QBALL_RECON_PI / m_NumberOfEquatorSamplingPoints);
527  (*C)(0,i) = cos(theta);
528  (*C)(1,i) = sin(theta);
529  (*C)(2,i) = NumericTraits<double>::Zero;
530  }
531 
532  // rotate the sampling points to each directions of interest
533  vnl_matrix<double> *S
534  = new vnl_matrix<double>(3,m_NumberOfEquatorSamplingPoints*NOdfDirections);
535 
536  {
537  vnl_vector_fixed<double,3> z(NumericTraits<double>::Zero);
538  z.put(2,NumericTraits<double>::One);
539  vnl_matrix_fixed<double,3,3> eye(NumericTraits<double>::Zero);
540  eye.fill_diagonal(NumericTraits<double>::One);
541  for(int i=0; i<NOdfDirections; i++)
542  {
543  vnl_vector_fixed<double,3> ui = (*U).get_column(i);
544  vnl_matrix<double> *RC
545  = new vnl_matrix<double>(3,m_NumberOfEquatorSamplingPoints);
546 
547  if( (z(0)*ui(0)+z(1)*ui(1)+z(2)*ui(2)+1) != 0 )
548  {
549  vnl_matrix_fixed<double,3,3> R;
550  R.set_column(0, (z+ui)*(z(0)+ui(0)));
551  R.set_column(1, (z+ui)*(z(1)+ui(1)));
552  R.set_column(2, (z+ui)*(z(2)+ui(2)));
553  R /= (z(0)*ui(0)+z(1)*ui(1)+z(2)*ui(2)+1);
554  R -= eye;
555  (*RC) = R*(*C);
556  }
557  else
558  {
559  RC = C;
560  }
561  (*S).set_columns(i*m_NumberOfEquatorSamplingPoints, *RC);
562  }
563  }
564 
565  // determine interpolation kernel width first
566  // use to determine diffusion measurement contribution to each of the kernels
567  vnl_matrix<double> *H_plus
568  = new vnl_matrix<double>(NBasisFunctionCenters,m_NumberOfGradientDirections);
569 
570  double maxSigma = QBALL_RECON_PI/6;
571  double bestSigma = maxSigma;
572 
573  {
574  double stepsize = 0.01;
575  double start = 0.01;
576  double minCondition = NumericTraits<double>::max();
577 
578  vnl_matrix<double> *H
579  = new vnl_matrix<double>(m_NumberOfGradientDirections,NBasisFunctionCenters,(double)0);
580 
581  {
582 
583  int increasing = 0;
584  for( double sigma=start; sigma<maxSigma; sigma+=stepsize)
585  {
586  vnl_matrix<double> *tmpH
587  = new vnl_matrix<double>(m_NumberOfGradientDirections,NBasisFunctionCenters);
588 
589  for(unsigned int r=0; r<m_NumberOfGradientDirections; r++)
590  {
591  for(int c=0; c<NBasisFunctionCenters; c++)
592  {
593  double qtv = (*Q)(0,r)*(*V)(0,c)
594  + (*Q)(1,r)*(*V)(1,c)
595  + (*Q)(2,r)*(*V)(2,c);
596  qtv = (qtv<-1.0) ? -1.0 : ( (qtv>1.0) ? 1.0 : qtv);
597  double x = acos(qtv);
598  (*tmpH)(r,c) = (1.0/(sigma*sqrt(2.0*QBALL_RECON_PI)))
599  *exp((-x*x)/(2*sigma*sigma));
600  }
601  }
602 
603  vnl_svd<double> *solver;
604  if(m_NumberOfGradientDirections>NBasisFunctionCenters)
605  {
606  solver = new vnl_svd<double>(*tmpH);
607  }
608  else
609  {
610  solver = new vnl_svd<double>(tmpH->transpose());
611  }
612  double condition = solver->sigma_max() / solver->sigma_min();
613 
614  std::cout << sigma << ": " << condition << std::endl;
615 
616  if( condition < minCondition )
617  {
618  minCondition = condition;
619  bestSigma = sigma;
620  H->update(*tmpH);
621  }
622  else
623  {
624  // optimum assumed to be hit after condition increased 3 times
625  if (++increasing>3) break;
626  }
627  }
628  }
629 
630 
631  vnl_matrix_inverse<double> *pseudoInverse
632  = new vnl_matrix_inverse<double>(*H);
633  (*H_plus) = pseudoInverse->pinverse();
634 
635  std::cout << "choosing sigma = " << bestSigma << std::endl;
636 
637  }
638 
639  // this is the contribution of each kernel to each sampling point on the
640  // equator
641  vnl_matrix<double> *G
642  = new vnl_matrix<double>(m_NumberOfEquatorSamplingPoints*NOdfDirections,NBasisFunctionCenters);
643 
644  {
645  for(unsigned int r=0; r<m_NumberOfEquatorSamplingPoints*NOdfDirections; r++)
646  {
647  for(int c=0; c<NBasisFunctionCenters; c++)
648  {
649  double stv = (*S)(0,r)*(*V)(0,c)
650  + (*S)(1,r)*(*V)(1,c)
651  + (*S)(2,r)*(*V)(2,c);
652  stv = (stv<-1.0) ? -1.0 : ( (stv>1.0) ? 1.0 : stv);
653  double x = acos(stv);
654  (*G)(r,c) = (1.0/(bestSigma*sqrt(2.0*QBALL_RECON_PI)))
655  *exp((-x*x)/(2*bestSigma*bestSigma));
656  }
657  }
658  }
659 
660  vnl_matrix<double> *GH_plus =
661  new vnl_matrix<double>(m_NumberOfEquatorSamplingPoints*NOdfDirections,m_NumberOfGradientDirections);
662 
663  // simple matrix multiplication, manual cause of stack overflow using operator
664  for (unsigned i = 0; i < m_NumberOfEquatorSamplingPoints*NOdfDirections; ++i)
665  for (unsigned j = 0; j < m_NumberOfGradientDirections; ++j)
666  {
667  double accum = (*G)(i,0) * (*H_plus)(0,j);
668  for (unsigned k = 1; k < NOdfDirections; ++k)
669  accum += (*G)(i,k) * (*H_plus)(k,j);
670  (*GH_plus)(i,j) = accum;
671  }
672 
673  typename vnl_matrix<double>::iterator it3;
674  for( it3 = (*GH_plus).begin(); it3 != (*GH_plus).end(); it3++)
675  {
676  if(*it3<0.0)
677  *it3 = 0;
678  }
679 
680  // this is an addition to the original tuch algorithm
681  for(unsigned int i=0; i<NOdfDirections*m_NumberOfEquatorSamplingPoints; i++)
682  {
683  vnl_vector< double > r = GH_plus->get_row(i);
684  r /= r.sum();
685  GH_plus->set_row(i,r);
686  }
687 
688  m_ReconstructionMatrix
689  = new vnl_matrix<TOdfPixelType>(NOdfDirections,m_NumberOfGradientDirections,0.0);
690  for(int i=0; i<NOdfDirections; i++)
691  {
692  for(unsigned int j=0; j<m_NumberOfGradientDirections; j++)
693  {
694  for(unsigned int k=0; k<m_NumberOfEquatorSamplingPoints; k++)
695  {
696  (*m_ReconstructionMatrix)(i,j) += (TOdfPixelType)(*GH_plus)(m_NumberOfEquatorSamplingPoints*i+k,j);
697  }
698  }
699  }
700 
701  // this is also an addition to the original tuch algorithm
702  for(int i=0; i<NOdfDirections; i++)
703  {
704  vnl_vector< TOdfPixelType > r = m_ReconstructionMatrix->get_row(i);
705  r /= r.sum();
706  m_ReconstructionMatrix->set_row(i,r);
707  }
708  std::cout << "Reconstruction Matrix computed." << std::endl;
709 
710  }
711 
712  template< class TReferenceImagePixelType,
713  class TGradientImagePixelType,
714  class TOdfPixelType,
715  int NrOdfDirections,
716  int NrBasisFunctionCenters>
717  void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
718  TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
719  NrBasisFunctionCenters>
720  ::AddGradientImage( const GradientDirectionType &gradientDirection,
721  const GradientImageType *gradientImage )
722  {
723  // Make sure crazy users did not call both AddGradientImage and
724  // SetGradientImage
725  if( m_GradientImageTypeEnumeration == GradientIsInASingleImage)
726  {
727  itkExceptionMacro( << "Cannot call both methods:"
728  << "AddGradientImage and SetGradientImage. Please call only one of them.");
729  }
730 
731  // If the container to hold the gradient directions hasn't been allocated
732  // yet, allocate it.
733  if( !this->m_GradientDirectionContainer )
734  {
735  this->m_GradientDirectionContainer = GradientDirectionContainerType::New();
736  }
737 
738  this->m_NumberOfGradientDirections = m_GradientDirectionContainer->Size();
739 
740  m_GradientDirectionContainer->InsertElement( this->m_NumberOfGradientDirections,
741  gradientDirection / gradientDirection.two_norm() );
742 
743  this->ProcessObject::SetNthInput( this->m_NumberOfGradientDirections,
744  const_cast< GradientImageType* >(gradientImage) );
745 
746  m_GradientImageTypeEnumeration = GradientIsInManyImages;
747  }
748 
749  template< class TReferenceImagePixelType,
750  class TGradientImagePixelType,
751  class TOdfPixelType,
752  int NrOdfDirections,
753  int NrBasisFunctionCenters>
754  void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
755  TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
756  NrBasisFunctionCenters>
757  ::SetGradientImage(const GradientDirectionContainerType *gradientDirection,
758  const GradientImagesType *gradientImage )
759  {
760  // Make sure crazy users did not call both AddGradientImage and
761  // SetGradientImage
762  if( m_GradientImageTypeEnumeration == GradientIsInManyImages )
763  {
764  itkExceptionMacro( << "Cannot call both methods:"
765  << "AddGradientImage and SetGradientImage. Please call only one of them.");
766  }
767 
768  this->m_GradientDirectionContainer = GradientDirectionContainerType::New();
769  for(GradientDirectionContainerType::ConstIterator it = gradientDirection->Begin();
770  it != gradientDirection->End(); it++)
771  {
772  this->m_GradientDirectionContainer->push_back(it.Value());
773  }
774 
775  unsigned int numImages = gradientDirection->Size();
776  this->m_NumberOfBaselineImages = 0;
777  for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin();
778  it != this->m_GradientDirectionContainer->End(); it++)
779  {
780  if(it.Value().one_norm() <= 0.0)
781  {
782  this->m_NumberOfBaselineImages++;
783  }
784  else // Normalize non-zero gradient directions
785  {
786  it.Value() = it.Value() / it.Value().two_norm();
787  }
788  }
789 
790  this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages;
791 
792  // ensure that the gradient image we received has as many components as
793  // the number of gradient directions
794  if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections )
795  {
796  itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages
797  << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages
798  << " directions specified but image has " << gradientImage->GetVectorLength()
799  << " components.");
800  }
801 
802  this->ProcessObject::SetNthInput( 0,
803  const_cast< GradientImagesType* >(gradientImage) );
804  m_GradientImageTypeEnumeration = GradientIsInASingleImage;
805  }
806 
807 
808  template< class TReferenceImagePixelType,
809  class TGradientImagePixelType,
810  class TOdfPixelType,
811  int NrOdfDirections,
812  int NrBasisFunctionCenters>
813  void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
814  TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
815  NrBasisFunctionCenters>
816  ::PrintSelf(std::ostream& os, Indent indent) const
817  {
818  std::locale C("C");
819  std::locale originalLocale = os.getloc();
820  os.imbue(C);
821 
822  Superclass::PrintSelf(os,indent);
823 
824  os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl;
825  if ( m_GradientDirectionContainer )
826  {
827  os << indent << "GradientDirectionContainer: "
828  << m_GradientDirectionContainer << std::endl;
829  }
830  else
831  {
832  os << indent <<
833  "GradientDirectionContainer: (Gradient directions not set)" << std::endl;
834  }
835  os << indent << "NumberOfGradientDirections: " <<
836  m_NumberOfGradientDirections << std::endl;
837  os << indent << "NumberOfBaselineImages: " <<
838  m_NumberOfBaselineImages << std::endl;
839  os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl;
840  os << indent << "BValue: " << m_BValue << std::endl;
841  if ( this->m_GradientImageTypeEnumeration == GradientIsInManyImages )
842  {
843  os << indent << "Gradient images haven been supplied " << std::endl;
844  }
845  else if ( this->m_GradientImageTypeEnumeration == GradientIsInManyImages )
846  {
847  os << indent << "A multicomponent gradient image has been supplied" << std::endl;
848  }
849 
850  os.imbue( originalLocale );
851  }
852 
853 }
854 
855 #endif // __itkDiffusionQballReconstructionImageFilter_txx