Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkRegularizedIVIMReconstructionSingleIteration.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 #ifndef _itkRegularizedIVIMReconstructionSingleIteration_txx
18 #define _itkRegularizedIVIMReconstructionSingleIteration_txx
19 
20 // itk includes
21 #include "itkConstShapedNeighborhoodIterator.h"
22 #include "itkNeighborhoodInnerProduct.h"
23 #include "itkImageRegionIterator.h"
24 #include "itkNeighborhoodAlgorithm.h"
25 #include "itkZeroFluxNeumannBoundaryCondition.h"
26 #include "itkOffset.h"
27 #include "itkProgressReporter.h"
28 #include "itkRegularizedIVIMLocalVariationImageFilter.h"
29 
30 // other includes
31 #include <vector>
32 #include <algorithm>
33 
34 #define IVIM_FOO -100000
35 
36 namespace itk
37 {
38 
39  /**
40  * constructor
41  */
42  template <class TInputPixel, class TOutputPixel, class TRefPixelType>
43  RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
44  ::RegularizedIVIMReconstructionSingleIteration()
45  {
46  m_Lambda = 1.0;
47  m_LocalVariation = LocalVariationImageType::New();
48  }
49 
50  /**
51  * generate requested region
52  */
53  template <class TInputPixel, class TOutputPixel, class TRefPixelType>
54  void
55  RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
56  ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
57  {
58  // call the superclass' implementation of this method
59  Superclass::GenerateInputRequestedRegion();
60 
61  // get pointers to the input and output
62  typename Superclass::InputImagePointer inputPtr =
63  const_cast< InputImageType * >( this->GetInput() );
64  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
65 
66  if ( !inputPtr || !outputPtr )
67  {
68  return;
69  }
70 
71  // get a copy of the input requested region (should equal the output
72  // requested region)
73  typename InputImageType::RegionType inputRequestedRegion;
74  inputRequestedRegion = inputPtr->GetRequestedRegion();
75 
76  // pad the input requested region by 1
77  inputRequestedRegion.PadByRadius( 1 );
78 
79  // crop the input requested region at the input's largest possible region
80  if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
81  {
82  inputPtr->SetRequestedRegion( inputRequestedRegion );
83  return;
84  }
85  else
86  {
87  // Couldn't crop the region (requested region is outside the largest
88  // possible region). Throw an exception.
89 
90  // store what we tried to request (prior to trying to crop)
91  inputPtr->SetRequestedRegion( inputRequestedRegion );
92 
93  // build an exception
94  InvalidRequestedRegionError e(__FILE__, __LINE__);
95  e.SetLocation(ITK_LOCATION);
96  e.SetDescription("Requested region outside possible region.");
97  e.SetDataObject(inputPtr);
98  throw e;
99  }
100  }
101 
102 
103  /**
104  * generate output
105  */
106  template <class TInputPixel, class TOutputPixel, class TRefPixelType>
107  void
108  RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
109  ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
110  ThreadIdType)
111  {
112 
113  typename OutputImageType::Pointer output = this->GetOutput();
114  typename InputImageType::ConstPointer input = this->GetInput();
115 
116  // Find the data-set boundary "faces"
117  itk::Size<InputImageDimension> size;
118  for( int i=0; i<InputImageDimension; i++)
119  size[i] = 1;
120 
121  NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
122  typename NeighborhoodAlgorithm::
123  ImageBoundaryFacesCalculator<InputImageType>::FaceListType
124  faceList = bC(input, outputRegionForThread, size);
125 
126  NeighborhoodAlgorithm::
127  ImageBoundaryFacesCalculator<LocalVariationImageType> lv_bC;
128  typename NeighborhoodAlgorithm::
129  ImageBoundaryFacesCalculator<LocalVariationImageType>::FaceListType
130  lv_faceList = lv_bC(m_LocalVariation, outputRegionForThread, size);
131 
132  ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
133  ZeroFluxNeumannBoundaryCondition<LocalVariationImageType> lv_nbc;
134  std::vector<double> ws;
135  std::vector<double> hs;
136 
137  typename NeighborhoodAlgorithm::
138  ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator
139  lv_fit=lv_faceList.begin();
140 
141  // Process each of the boundary faces. These are N-d regions which border
142  // the edge of the buffer.
143  for ( typename NeighborhoodAlgorithm::
144  ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator
145  fit=faceList.begin(); fit != faceList.end(); ++fit)
146  {
147 
148  // iterators over output, input, original and local variation image
149  ImageRegionIterator<OutputImageType> output_image_it =
150  ImageRegionIterator<OutputImageType>(output, *fit);
151  ImageRegionConstIterator<InputImageType> input_image_it =
152  ImageRegionConstIterator<InputImageType>(input, *fit);
153  ImageRegionConstIterator<RefImageType> orig_image_it =
154  ImageRegionConstIterator<RefImageType>(m_OriginalImage, *fit);
155  ImageRegionConstIterator<LocalVariationImageType> loc_var_image_it =
156  ImageRegionConstIterator<LocalVariationImageType>(
157  m_LocalVariation, *fit);
158 
159  // neighborhood in input image
160  ConstShapedNeighborhoodIterator<InputImageType>
161  input_image_neighbors_it(size, input, *fit);
162  typename ConstShapedNeighborhoodIterator<InputImageType>::
163  OffsetType offset;
164  input_image_neighbors_it.OverrideBoundaryCondition(&nbc);
165  input_image_neighbors_it.ClearActiveList();
166  for(int i=0; i<InputImageDimension; i++)
167  {
168  offset.Fill(0);
169  offset[i] = -1;
170  input_image_neighbors_it.ActivateOffset(offset);
171  offset[i] = 1;
172  input_image_neighbors_it.ActivateOffset(offset);
173  }
174  input_image_neighbors_it.GoToBegin();
175 
176  // neighborhood in local variation image
177  ConstShapedNeighborhoodIterator<LocalVariationImageType>
178  loc_var_image_neighbors_it(size, m_LocalVariation, *lv_fit);
179  loc_var_image_neighbors_it.OverrideBoundaryCondition(&lv_nbc);
180  loc_var_image_neighbors_it.ClearActiveList();
181  for(int i=0; i<InputImageDimension; i++)
182  {
183  offset.Fill(0);
184  offset[i] = -1;
185  loc_var_image_neighbors_it.ActivateOffset(offset);
186  offset[i] = 1;
187  loc_var_image_neighbors_it.ActivateOffset(offset);
188  }
189  loc_var_image_neighbors_it.GoToBegin();
190 
191  const unsigned int neighborhoodSize = InputImageDimension*2;
192  ws.resize(neighborhoodSize);
193 
194  while ( ! output_image_it.IsAtEnd() )
195  {
196 
197  //MITK_INFO << "looking at voxel " << output_image_it.GetIndex();
198 
199  // 1 / ||nabla_alpha(u)||_a
200  double locvar_alpha_inv = 1.0/loc_var_image_it.Get();
201  //MITK_INFO << "locvar: " << loc_var_image_it.Get();
202  // compute w_alphabetas
203  int count = 0;
204  double wsum = 0;
205  typename ConstShapedNeighborhoodIterator<LocalVariationImageType>::
206  ConstIterator loc_var_neighbors_it;
207  for (loc_var_neighbors_it = loc_var_image_neighbors_it.Begin();
208  ! loc_var_neighbors_it.IsAtEnd();
209  loc_var_neighbors_it++)
210  {
211  // w_alphabeta(u) =
212  // 1 / ||nabla_alpha(u)||_a + 1 / ||nabla_beta(u)||_a
213  ws[count] =
214  locvar_alpha_inv + (1.0/(double)loc_var_neighbors_it.Get());
215  wsum += ws[count++];
216 
217  //MITK_INFO << "nb var: " << count << ": " << loc_var_neighbors_it.Get();
218  }
219  //MITK_INFO << "wsum: " << wsum;
220 
221  // h_alphaalpha * u_alpha^zero
222  typename RefImageType::PixelType orig = orig_image_it.Get();
223 // vnl_vector<double> estim(orig.GetSize());
224 // vnl_matrix<double> diff(orig.GetSize(),1);
225 // vnl_matrix<double> estimdash(orig.GetSize(),2);
226  vnl_vector_fixed<double,2> step;
227  step[0] = 0; step[1]=0;
228  for(size_t ind=0; ind<m_BValues.size(); ind++)
229  {
230  //MITK_INFO << "refval: " << orig[ind];
231  double estim = (1-input_image_it.Get()[0])*exp(-m_BValues[ind]*input_image_it.Get()[1]);
232  double estimdash1 = exp(-m_BValues[ind]*input_image_it.Get()[1]);
233  double estimdash2 = (-1.0) * (1.0-input_image_it.Get()[0]) * m_BValues[ind] * exp(-m_BValues[ind]*input_image_it.Get()[1]);
234  //MITK_INFO << "estimdash1: " << estimdash1 << "; estimdash2: " << estimdash2 << "; diff: " << (double) orig[ind] - estim;
235  if(orig[ind] != IVIM_FOO)
236  {
237  step[0] += ((double) orig[ind] - estim) * estimdash2;
238  step[1] += ((double) orig[ind] - estim) * estimdash1;
239  }
240  }
241 
242  step[1] *= m_Lambda / (m_Lambda+wsum);
243 
244  // add the different h_alphabeta * u_beta
245  count = 0;
246  typename ConstShapedNeighborhoodIterator<InputImageType>::
247  ConstIterator input_neighbors_it;
248  for (input_neighbors_it = input_image_neighbors_it.Begin();
249  ! input_neighbors_it.IsAtEnd();
250  input_neighbors_it++)
251  {
252  step[1] += (input_neighbors_it.Get()[1] - input_image_it.Get()[1]) * (ws[count++] / (m_Lambda+wsum));
253  }
254 
255  //MITK_INFO << "stepfinal: " << step[0] << "; " << step[1];
256 
257  // set output result
258  OutputPixelType out;
259  out[0] = input_image_it.Get()[0] + .001*step[0];
260  out[1] = input_image_it.Get()[1] + .00001*step[1];
261  output_image_it.Set( out );
262 
263  //MITK_INFO << "(" << input_image_it.Get()[0] << " ; " << input_image_it.Get()[1] << ") => (" << out[0] << " ; " << out[1] << ")";
264  // increment iterators
265  ++output_image_it;
266  ++input_image_it;
267  ++orig_image_it;
268  ++loc_var_image_it;
269  ++input_image_neighbors_it;
270  ++loc_var_image_neighbors_it;
271 
272  }
273 
274  ++lv_fit;
275  }
276  }
277 
278  /**
279  * first calculate local variation in the image
280  */
281  template <class TInputPixel, class TOutputPixel, class TRefPixelType>
282  void
283  RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
284  ::BeforeThreadedGenerateData()
285  {
286  typedef typename itk::RegularizedIVIMLocalVariationImageFilter
287  <InputImageType,LocalVariationImageType> FilterType;
288  typename FilterType::Pointer filter = FilterType::New();
289  filter->SetInput(this->GetInput(0));
290  filter->SetNumberOfThreads(this->GetNumberOfThreads());
291  filter->Update();
292  this->m_LocalVariation = filter->GetOutput();
293  }
294 
295  /**
296  * Standard "PrintSelf" method
297  */
298  template <class TInputPixel, class TOutputPixel, class TRefPixelType>
299  void
300  RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
301  ::PrintSelf(
302  std::ostream& os,
303  Indent indent) const
304  {
305  Superclass::PrintSelf( os, indent );
306  }
307 
308 } // end namespace itk
309 
310 #endif