CBICA Toolkit  1.0.0
itkDiffusionTensor3DReconstructionImageFilter.h
Go to the documentation of this file.
1 /**
2 \file itkDiffusionTensor3DReconstructionImageFilter.h
3 
4 \brief Declaration of DiffusionTensor3DReconstructionImageFilter
5 
6 https://www.cbica.upenn.edu/sbia/software/ <br>
7 software@cbica.upenn.edu
8 
9 Copyright (c) 2015 University of Pennsylvania. All rights reserved. <br>
10 See COPYING file or https://www.cbica.upenn.edu/sbia/software/license.html
11 
12 */
13 
14 #pragma once
15 
16 #include "itkImageToImageFilter.h"
17 #include "itkDiffusionTensor3D.h"
18 #include "vnl/vnl_matrix.h"
19 #include "vnl/vnl_vector_fixed.h"
20 #include "vnl/vnl_matrix_fixed.h"
21 #include "vnl/algo/vnl_svd.h"
22 #include "itkVectorContainer.h"
23 #include "itkVectorImage.h"
24 #include "itkSpatialObject.h"
25 #include "itkNumericTraits.h"
26 
27 #if WIN32
28 __declspec(dllexport) inline void getRidOfLNK4221(){};
29 #endif
30 
31 namespace itk
32 {
33 /**
34 \class DiffusionTensor3DReconstructionImageFilter
35 \brief This class takes as input one or more reference image (acquired in the
36 absence of diffusion sensitizing gradients) and 'n' diffusion
37 weighted images and their gradient directions and computes an image of
38 tensors. (with DiffusionTensor3D as the pixel type). Once that is done, you
39 can apply filters on this tensor image to compute FA, ADC, RGB weighted
40 maps etc.
41 
42  * sbia additions to itk DiffusionTensor3DRecon Filter
43  added a spatial Object as a mask.
44  added cludge to fix nonSPD tensors. at the moment only this just
45  takes abs of any negative eigenvalues. Posible future direction is to
46  include spd constrain when this occurs.
47 
48 \par Inputs and Usage
49 There are two ways to use this class. When you have one reference image and \c n
50 gradient images, you would use the class as
51 \code
52  filter->SetReferenceImage( image0 );
53  filter->AddGradientImage( direction1, image1 );
54  filter->AddGradientImage( direction2, image2 );
55  ...
56 \endcode
57 
58 \par
59 When you have the 'n' gradient and one or more reference images in a single
60 multi-component image (VectorImage), you can specify the images simply as
61 \code
62  filter->SetGradientImage( directionsContainer, vectorImage );
63 \endcode
64 Note that this method is used to specify both the reference and gradient images.
65 This is convenient when the DWI images are read in using the
66 <a href="http://wiki.na-mic.org/Wiki/index.php/NAMIC_Wiki:DTI:Nrrd_format">NRRD</a>
67 format. Like the Nrrd format, the reference images are those components of the
68 vectorImage whose gradient direction is (0,0,0). If more than one reference image
69 is present, they are averaged prior to applying the Stejskal-Tanner equations.
70 
71 \par Outputs
72 The output image is an image of Tensors:
73 \code
74  Image< DiffusionTensor3D< TTensorPixelType >, 3 >
75 \endcode
76 
77 \par Parameters
78 \li Threshold - Threshold on the reference image data. The output tensor will
79 be a null tensor for pixels in the reference image that have a value less
80 than this.
81 \li BValue - See the documentation of SetBValue().
82 \li At least 6 gradient images must be specified for the filter to be able
83 to run.
84 
85 
86 \par Template parameters
87 The class is templated over the pixel type of the reference and gradient
88 images (expected to be scalar data types) and the internal representation
89 of the DiffusionTensor3D pixel (double, float etc).
90 
91 \par References:
92 \li<a href="http://lmi.bwh.harvard.edu/papers/pdfs/2002/westinMEDIA02.pdf">[1]</a>
93 <em>C.F.Westin, S.E.Maier, H.Mamata, A.Nabavi, F.A.Jolesz, R.Kikinis,
94 "Processing and visualization for Diffusion tensor MRI", Medical image
95 Analysis, 2002, pp 93-108.</em>
96 \li<a href="splweb.bwh.harvard.edu:8000/pages/papers/westin/ISMRM2002.pdf">[2]</a>
97 <em>A Dual Tensor Basis Solution to the Stejskal-Tanner Equations for DT-MRI</em>
98 
99 \par WARNING:
100 Although this filter has been written to support multiple threads, please
101 set the number of threads to 1.
102 \code
103  filter->SetNumberOfThreads(1);
104 \endcode
105 This is due to buggy code in netlib/dsvdc, that is called by vnl_svd.
106 (used to compute the psudo-inverse to find the dual tensor basis).
107 
108 \author Thanks to Xiaodong Tao, GE, for contributing parts of this class. Also
109 thanks to Casey Goodlet, UNC for patches to support multiple baseline images
110 and other improvements.
111 
112 \note
113 This work is part of the National Alliance for Medical image Computing
114 (NAMIC), funded by the National Institutes of Health through the NIH Roadmap
115 for Medical Research, Grant U54 EB005149.
116 
117 \par Examples and Datasets
118 See Examples/Filtering/DiffusionTensor3DReconstructionImageFilter.cxx
119 Sample DTI datasets may be obtained from
120 \begin verbatim
121  ftp://public.kitware.com/pub/namic/DTI/Data/dwi.nhdr
122  ftp://public.kitware.com/pub/namic/DTI/Data/dwi.img.gz ( gunzip this )
123 \end verbatim
124 
125 
126 \sa DiffusionTensor3D SymmetricSecondRankTensor
127 \ingroup Multithreaded TensorObjects
128  */
129 
130 template< class TReferenceImagePixelType,
131  class TGradientImagePixelType=TReferenceImagePixelType,
132  class TTensorPixelType=double >
134  public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >,
135  Image< DiffusionTensor3D< TTensorPixelType >, 3 > >
136 {
137 
138 public:
139 
141  typedef SmartPointer<Self> Pointer;
142  typedef SmartPointer<const Self> ConstPointer;
143  typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>,
144  Image< DiffusionTensor3D< TTensorPixelType >, 3 > >
145  Superclass;
146 
147  /** Method for creation through the object factory. */
148  itkNewMacro(Self);
149 
150  /** Runtime information support. */
152  ImageToImageFilter);
153 
154  typedef TReferenceImagePixelType ReferencePixelType;
155 
156  typedef TGradientImagePixelType GradientPixelType;
157 
158  typedef DiffusionTensor3D< TTensorPixelType > TensorPixelType;
159 
160  typedef typename TensorPixelType::EigenValuesArrayType
161  TensorEigenValuesType;
162 
163  typedef typename TensorPixelType::EigenVectorsMatrixType
164  TensorEigenVectorsType;
165 
166  typedef typename TensorPixelType::MatrixType
167  TensorMatrixType;
168 
169 
170  /** Reference image data, This image is acquired in the absence
171  * of a diffusion sensitizing field gradient */
172  typedef typename Superclass::InputImageType ReferenceImageType;
173 
174  typedef Image< TensorPixelType, 3 > TensorImageType;
175 
176  typedef TensorImageType OutputImageType;
177 
178  typedef typename Superclass::OutputImageRegionType
179  OutputImageRegionType;
180 
181  /** Typedef defining one (of the many) gradient images. */
182  typedef Image< GradientPixelType, 3 > GradientImageType;
183 
184  /** An alternative typedef defining one (of the many) gradient images.
185  * It will be assumed that the vectorImage has the same dimension as the
186  * Reference image and a vector length parameter of \c n (number of
187  * gradient directions) */
188  typedef VectorImage< GradientPixelType, 3 > GradientImagesType;
189 
190  /** Typedef defining the residual image.*/
191 // typedef VectorImage< typename NumericTraits<GradientPixelType>::RealType, 3 >
192  typedef VectorImage< GradientPixelType, 3 >
194  typedef typename ResidualImageType::PixelType ResidualPixelType;
195 
196  /** Type for the mask of the fixed image. Only pixels that are "inside"
197  this mask will be considered for the computation of the metric */
198  typedef SpatialObject< 3 > ImageMaskType;
199  typedef typename ImageMaskType::Pointer ImageMaskPointer;
200 
201 
202  /** Holds the tensor basis coefficients G_k */
203  typedef vnl_matrix_fixed< double, 6, 6 > TensorBasisMatrixType;
204  typedef vnl_matrix_fixed< double, 6, 6 > TensorBasisInverseMatrixType;
205  typedef vnl_matrix< double > CoefficientMatrixType;
206 
207  /** Holds each magnetic field gradient used to acquire one DWImage */
208  typedef vnl_vector_fixed< double, 3 > GradientDirectionType;
209 
210  /** Container to hold gradient directions of the 'n' DW measurements */
211  typedef VectorContainer< unsigned int,
213 
214 
215  /** Set method to add a gradient direction and its corresponding image. */
216  void AddGradientImage( const GradientDirectionType &, const GradientImageType *image);
217 
218  /** Another set method to add a gradient directions and its corresponding
219  * image. The image here is a VectorImage. The user is expected to pass the
220  * gradient directions in a container. The ith element of the container
221  * corresponds to the gradient direction of the ith component image the
222  * VectorImage. For the baseline image, a vector of all zeros
223  * should be set. */
224  void SetGradientImage( GradientDirectionContainerType *,
225  const GradientImagesType *image);
226 
227  /** Set method to set the reference image. */
228  void SetReferenceImage( ReferenceImageType *referenceImage )
229  {
230  if( m_GradientImageTypeEnumeration == GradientIsInASingleImage)
231  {
232  itkExceptionMacro( << "Cannot call both methods:"
233  << "AddGradientImage and SetGradientImage. Please call only one of them.");
234  }
235 
236  this->ProcessObject::SetNthInput( 0, referenceImage );
237 
238  m_GradientImageTypeEnumeration = GradientIsInManyImages;
239  }
240 
241  /** Get reference image */
243  { return ( static_cast< ReferenceImageType *>(this->ProcessObject::GetInput(0)) ); }
244 
245  /** Return the gradient direction. idx is 0 based */
246  virtual GradientDirectionType GetGradientDirection( unsigned int idx) const
247  {
248  if( idx >= m_NumberOfGradientDirections )
249  {
250  itkExceptionMacro( << "Gradient direction " << idx << "does not exist" );
251  }
252  return m_GradientDirectionContainer->ElementAt( idx+1 );
253  }
254 
255  /** Threshold on the reference image data. The output tensor will be a null
256  * tensor for pixels in the reference image that have a value less than this
257  * threshold. */
258  itkSetMacro( Threshold, ReferencePixelType );
259  itkGetMacro( Threshold, ReferencePixelType );
260 
261  /** Set/Get the image mask. */
262  itkSetObjectMacro( ImageMask, ImageMaskType );
263  itkGetConstObjectMacro( ImageMask, ImageMaskType );
264 
265  /** Get/set compute residuals flag. */
266  itkSetMacro( CalculateResidualImage, bool );
267  itkGetMacro( CalculateResidualImage, bool );
268 
269  /** Get the Residual image . */
270  itkGetConstObjectMacro( ResidualImage, ResidualImageType );
271 
272  /**
273  * The BValue \f$ (s/mm^2) \f$ value used in normalizing the tensors to
274  * physically meaningful units. See equation (24) of the first reference for
275  * a description of how this is applied to the tensor estimation.
276  * Equation (1) of the same reference describes the physical significance.
277  */
278  itkSetMacro( BValue, TTensorPixelType);
279 #ifdef GetBValue
280 #undef GetBValue
281 #endif
282  itkGetConstReferenceMacro( BValue, TTensorPixelType);
283 
284 #ifdef ITK_USE_CONCEPT_CHECKING
285  /** Begin concept checking */
286  itkConceptMacro(ReferenceEqualityComparableCheck,
287  (Concept::EqualityComparable<ReferencePixelType>));
288  itkConceptMacro(TensorEqualityComparableCheck,
289  (Concept::EqualityComparable<TensorPixelType>));
290  itkConceptMacro(GradientConvertibleToDoubleCheck,
291  (Concept::Convertible<GradientPixelType, double>));
292  itkConceptMacro(DoubleConvertibleToTensorCheck,
293  (Concept::Convertible<double, TensorPixelType>));
294  itkConceptMacro(GradientReferenceAdditiveOperatorsCheck,
295  (Concept::AdditiveOperators<GradientPixelType, GradientPixelType,
296  ReferencePixelType>));
297  itkConceptMacro(ReferenceOStreamWritableCheck,
298  (Concept::OStreamWritable<ReferencePixelType>));
299  itkConceptMacro(TensorOStreamWritableCheck,
300  (Concept::OStreamWritable<TensorPixelType>));
301  /** End concept checking */
302 #endif
303 
304 protected:
307  void PrintSelf(std::ostream& os, Indent indent) const;
308 
309  void ComputeTensorBasis();
310 
311  void BeforeThreadedGenerateData();
312  void ThreadedGenerateData( const
313  OutputImageRegionType &outputRegionForThread, int);
314 
315  /** enum to indicate if the gradient image is specified as a single multi-
316  * component image or as several separate images */
317  typedef enum
318  {
319  GradientIsInASingleImage = 1,
320  GradientIsInManyImages,
321  Else
322  } GradientImageTypeEnumeration;
323 
324  /* method to compute the residual */
325  ResidualPixelType ComputeResidual( vnl_vector<double>, TensorPixelType, double);
326 
327 private:
328 
329  /* Tensor basis coefficients */
330  TensorBasisMatrixType m_TensorBasis;
331 
332  /* Tensor basis coefficients */
333  TensorBasisInverseMatrixType m_TensorBasisInverse;
334 
335  CoefficientMatrixType m_BMatrix;
336 
337  /** container to hold gradient directions */
338  GradientDirectionContainerType::Pointer m_GradientDirectionContainer;
339 
340  /** Number of gradient measurements */
341  unsigned int m_NumberOfGradientDirections;
342 
343  /** Number of baseline images */
344  unsigned int m_NumberOfBaselineImages;
345 
346  /** Threshold on the reference image data */
347  ReferencePixelType m_Threshold;
348 
349  /** LeBihan's b-value for normalizing tensors */
350  TTensorPixelType m_BValue;
351 
352  /** Image Mask */
353  mutable ImageMaskPointer m_ImageMask;
354 
355  /** Gradient image was specified in a single image or in multiple images */
356  GradientImageTypeEnumeration m_GradientImageTypeEnumeration;
357 
358  typename ResidualImageType::Pointer m_ResidualImage;
359  bool m_CalculateResidualImage;
360 
361 };
362 
363 }
364 
365 #ifndef ITK_MANUAL_INSTANTIATION
366 #include "itkDiffusionTensor3DReconstructionImageFilter.cpp"
367 #endif
VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType
Container to hold gradient directions of the 'n' DW measurements.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:212
Image< GradientPixelType, 3 > GradientImageType
Typedef defining one (of the many) gradient images.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:182
VectorImage< GradientPixelType, 3 > GradientImagesType
An alternative typedef defining one (of the many) gradient images.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:188
vnl_vector_fixed< double, 3 > GradientDirectionType
Holds each magnetic field gradient used to acquire one DWImage.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:208
SpatialObject< 3 > ImageMaskType
Type for the mask of the fixed image.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:198
virtual GradientDirectionType GetGradientDirection(unsigned int idx) const
Return the gradient direction.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:246
This class takes as input one or more reference image (acquired in the absence of diffusion sensitizi...
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:133
void SetReferenceImage(ReferenceImageType *referenceImage)
Set method to set the reference image.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:228
vnl_matrix_fixed< double, 6, 6 > TensorBasisMatrixType
Holds the tensor basis coefficients G_k.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:203
Superclass::InputImageType ReferenceImageType
Reference image data, This image is acquired in the absence of a diffusion sensitizing field gradient...
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:172
VectorImage< GradientPixelType, 3 > ResidualImageType
Typedef defining the residual image.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:193
virtual ReferenceImageType * GetReferenceImage()
Get reference image.
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:242