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