Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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