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