Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkTensorImageToDiffusionImageFilter.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 /*=========================================================================
17  2
18  3 Program: Tensor ToolKit - TTK
19  4 Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.txx $
20  5 Language: C++
21  6 Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
22  7 Version: $Revision: 68 $
23  8
24  9 Copyright (c) INRIA 2010. All rights reserved.
25  10 See LICENSE.txt for details.
26  11
27  12 This software is distributed WITHOUT ANY WARRANTY; without even
28  13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
29  14 PURPOSE. See the above copyright notices for more information.
30  15
31  16 =========================================================================*/
32 #ifndef _itk_TensorImageToDiffusionImageFilter_txx_
33 #define _itk_TensorImageToDiffusionImageFilter_txx_
34 #endif
35 
36 #include "itkTensorImageToDiffusionImageFilter.h"
37 #include "itkTensorToL2NormImageFilter.h"
38 #include "itkRescaleIntensityImageFilter.h"
39 #include <itkImageRegionIterator.h>
40 #include <itkImageRegionConstIterator.h>
41 
42 namespace itk
43 {
44 
45 template <class TInputScalarType, class TOutputScalarType>
46 void
47 TensorImageToDiffusionImageFilter<TInputScalarType, TOutputScalarType>
48 ::BeforeThreadedGenerateData()
49 {
50 
51  if( m_GradientList->Size()==0 )
52  {
53  throw itk::ExceptionObject (__FILE__,__LINE__,"Error: gradient list is empty, cannot generate DWI.");
54  }
55 
56  if( m_BaselineImage.IsNull() )
57  {
58  // create a B0 image by taking the norm of the tensor field * scale:
59  typedef itk::TensorToL2NormImageFilter<InputImageType, itk::Image<InputScalarType,3> >
60  TensorToL2NormFilterType;
61 
62  typename TensorToL2NormFilterType::Pointer myFilter1 = TensorToL2NormFilterType::New();
63  myFilter1->SetInput (this->GetInput());
64 
65  try
66  {
67  myFilter1->Update();
68  }
69  catch (itk::ExceptionObject &e)
70  {
71  std::cerr << e;
72  return;
73  }
74 
75  typename itk::RescaleIntensityImageFilter< itk::Image<InputScalarType,3>, BaselineImageType>::Pointer rescaler=
76  itk::RescaleIntensityImageFilter<itk::Image<InputScalarType,3>, BaselineImageType>::New();
77 
78  rescaler->SetOutputMinimum ( m_Min );
79  rescaler->SetOutputMaximum ( m_Max );
80  rescaler->SetInput ( myFilter1->GetOutput() );
81  try
82  {
83  rescaler->Update();
84  }
85  catch (itk::ExceptionObject &e)
86  {
87  std::cerr << e;
88  return;
89  }
90 
91  m_BaselineImage = rescaler->GetOutput();
92  }
93 
94  typename OutputImageType::Pointer outImage = OutputImageType::New();
95  outImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing
96  outImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin
97  outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction
98  outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion());
99  outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() );
100  outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() );
101  outImage->SetVectorLength(m_GradientList->Size());
102  outImage->Allocate();
103 
104  this->SetNumberOfRequiredOutputs (1);
105  this->SetNthOutput (0, outImage);
106 
107 }
108 
109 template <class TInputScalarType, class TOutputScalarType>
110 void TensorImageToDiffusionImageFilter<TInputScalarType, TOutputScalarType>
111 ::ThreadedGenerateData (const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId )
112 {
113  typedef ImageRegionIterator<OutputImageType> IteratorOutputType;
114  typedef ImageRegionConstIterator<InputImageType> IteratorInputType;
115  typedef ImageRegionConstIterator<BaselineImageType> IteratorBaselineType;
116 
117  unsigned long numPixels = outputRegionForThread.GetNumberOfPixels();
118  unsigned long step = numPixels/100;
119  unsigned long progress = 0;
120 
121  IteratorOutputType itOut (this->GetOutput(), outputRegionForThread);
122  IteratorInputType itIn (this->GetInput(), outputRegionForThread);
123  IteratorBaselineType itB0 (m_BaselineImage, outputRegionForThread);
124 
125  typedef ImageRegionConstIterator< MaskImageType > IteratorMaskImageType;
126  IteratorMaskImageType itMask;
127 
128  if( m_MaskImage.IsNotNull() )
129  {
130  itMask = IteratorMaskImageType( m_MaskImage, outputRegionForThread);
131  itMask.GoToBegin();
132  }
133 
134  if( threadId==0 )
135  {
136  this->UpdateProgress (0.0);
137  }
138 
139  while(!itIn.IsAtEnd())
140  {
141 
142  if( this->GetAbortGenerateData() )
143  {
144  throw itk::ProcessAborted(__FILE__,__LINE__);
145  }
146 
147  InputPixelType T = itIn.Get();
148 
149  BaselinePixelType b0 = itB0.Get();
150 
151  OutputPixelType out;
152  out.SetSize(m_GradientList->Size());
153  out.Fill(0);
154 
155  short maskvalue = 1;
156  if( m_MaskImage.IsNotNull() )
157  {
158  maskvalue = itMask.Get();
159  ++itMask;
160  }
161 
162  std::vector<unsigned int> b0_indices;
163  if( b0 > 0)
164  {
165  for( unsigned int i=0; i<m_GradientList->Size(); i++)
166  {
167  GradientType g = m_GradientList->at(i);
168 
169  // normalize vector so the following computations work
170  const double twonorm = g.two_norm();
171  if( twonorm < vnl_math::eps )
172  {
173  b0_indices.push_back(i);
174  continue;
175  }
176 
177  GradientType gn = g.normalize();
178 
179  InputPixelType S;
180  S[0] = gn[0]*gn[0];
181  S[1] = gn[1]*gn[0];
182  S[2] = gn[2]*gn[0];
183  S[3] = gn[1]*gn[1];
184  S[4] = gn[2]*gn[1];
185  S[5] = gn[2]*gn[2];
186 
187  const double res =
188  T[0]*S[0] +
189  2 * T[1]*S[1] + T[3]*S[3] +
190  2 * T[2]*S[2] + 2 * T[4]*S[4] + T[5]*S[5];
191 
192  // check for corrupted tensor
193  if (res>=0)
194  {
195  // estimate the bvalue from the base value and the norm of the gradient
196  // - because of this estimation the vector have to be normalized beforehand
197  // otherwise the modelled signal is wrong ( i.e. not scaled properly )
198  const double bval = m_BValue * twonorm * twonorm;
199  out[i] = static_cast<OutputScalarType>( maskvalue * 1.0 * b0 * exp ( -1.0 * bval * res ) );
200  }
201  }
202  }
203 
204  for(unsigned int idx = 0; idx < b0_indices.size(); idx++ )
205  {
206  out[b0_indices.at(idx)] = b0;
207  }
208 
209 
210  itOut.Set(out);
211 
212  if( threadId==0 && step>0)
213  {
214  if( (progress%step)==0 )
215  {
216  this->UpdateProgress ( double(progress)/double(numPixels) );
217  }
218  }
219 
220  ++progress;
221  ++itB0;
222  ++itIn;
223  ++itOut;
224 
225  }
226 
227  if( threadId==0 )
228  {
229  this->UpdateProgress (1.0);
230  }
231 }
232 
233 template <class TInputScalarType, class TOutputScalarType>
234 void
235 TensorImageToDiffusionImageFilter<TInputScalarType, TOutputScalarType>
236 ::UpdateOutputInformation()
237 {
238  // Calls to superclass updateoutputinformation
239  Superclass::UpdateOutputInformation();
240 
241  this->GetOutput()->SetVectorLength( m_GradientList->size() );
242 }
243 
244 } // end of namespace