Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkLocalVariationImageFilter.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 #ifndef _itkLocalVariationImageFilter_txx
23 #define _itkLocalVariationImageFilter_txx
24 #include "itkLocalVariationImageFilter.h"
25 
26 #include "itkConstShapedNeighborhoodIterator.h"
27 #include "itkImageRegionIterator.h"
28 #include "itkNeighborhoodAlgorithm.h"
29 #include "itkNeighborhoodInnerProduct.h"
30 #include "itkOffset.h"
31 #include "itkProgressReporter.h"
32 #include "itkVectorImage.h"
33 #include "itkZeroFluxNeumannBoundaryCondition.h"
34 
35 #include <algorithm>
36 #include <vector>
37 
38 namespace itk
39 {
40  template <class TInputImage, class TOutputImage>
41  LocalVariationImageFilter<TInputImage, TOutputImage>::LocalVariationImageFilter()
42  {
43  }
44 
45  template <class TInputImage, class TOutputImage>
46  void LocalVariationImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion() throw(
47  InvalidRequestedRegionError)
48  {
49  // call the superclass' implementation of this method
50  Superclass::GenerateInputRequestedRegion();
51 
52  // get pointers to the input and output
53  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
54  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
55 
56  if (!inputPtr || !outputPtr)
57  {
58  return;
59  }
60 
61  // get a copy of the input requested region (should equal the output
62  // requested region)
63  typename TInputImage::RegionType inputRequestedRegion;
64  inputRequestedRegion = inputPtr->GetRequestedRegion();
65 
66  // pad the input requested region by 1
67  inputRequestedRegion.PadByRadius(1);
68 
69  // crop the input requested region at the input's largest possible region
70  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
71  {
72  inputPtr->SetRequestedRegion(inputRequestedRegion);
73  return;
74  }
75  else
76  {
77  // Couldn't crop the region (requested region is outside the largest
78  // possible region). Throw an exception.
79 
80  // store what we tried to request (prior to trying to crop)
81  inputPtr->SetRequestedRegion(inputRequestedRegion);
82 
83  // build an exception
84  InvalidRequestedRegionError e(__FILE__, __LINE__);
85  e.SetLocation(ITK_LOCATION);
86  e.SetDescription("Requested region outside possible region.");
87  e.SetDataObject(inputPtr);
88  throw e;
89  }
90  }
91 
92  template <>
93  double SquaredEuclideanMetric<itk::VariableLengthVector<float>>::Calc(itk::VariableLengthVector<float> p)
94  {
95  return p.GetSquaredNorm();
96  }
97 
98  template <>
99  double SquaredEuclideanMetric<itk::VariableLengthVector<double>>::Calc(itk::VariableLengthVector<double> p)
100  {
101  return p.GetSquaredNorm();
102  }
103 
104  template <class TPixelType>
105  double SquaredEuclideanMetric<TPixelType>::Calc(TPixelType p)
106  {
107  return p * p;
108  }
109 
110  template <class TInputImage, class TOutputImage>
111  void LocalVariationImageFilter<TInputImage, TOutputImage>::ThreadedGenerateData(
112  const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId)
113  {
114  // Allocate output
115  typename OutputImageType::Pointer output = this->GetOutput();
116  typename InputImageType::ConstPointer input = this->GetInput();
117 
118  itk::Size<InputImageDimension> size;
119  for (int i = 0; i < InputImageDimension; i++)
120  size[i] = 1;
121 
122  // Find the data-set boundary "faces"
123  NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
124  typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType faceList =
125  bC(input, outputRegionForThread, size);
126 
127  // support progress methods/callbacks
128  ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
129 
130  ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
131  std::vector<InputPixelType> pixels;
132 
133  // Process each of the boundary faces. These are N-d regions which border
134  // the edge of the buffer.
135  for (auto fit = faceList.begin(); fit != faceList.end(); ++fit)
136  {
137  // iterators over output and input
138  ImageRegionIterator<OutputImageType> output_image_it(output, *fit);
139  ImageRegionConstIterator<InputImageType> input_image_it(input.GetPointer(), *fit);
140 
141  // neighborhood iterator for input image
142  ConstShapedNeighborhoodIterator<InputImageType> input_image_neighbors_it(size, input, *fit);
143  typename ConstShapedNeighborhoodIterator<InputImageType>::OffsetType offset;
144  input_image_neighbors_it.OverrideBoundaryCondition(&nbc);
145  input_image_neighbors_it.ClearActiveList();
146  for (int i = 0; i < InputImageDimension; i++)
147  {
148  offset.Fill(0);
149  offset[i] = -1;
150  input_image_neighbors_it.ActivateOffset(offset);
151  offset[i] = 1;
152  input_image_neighbors_it.ActivateOffset(offset);
153  }
154  input_image_neighbors_it.GoToBegin();
155  // const unsigned int neighborhoodSize = InputImageDimension*2;
156 
157  while (!input_image_neighbors_it.IsAtEnd())
158  {
159  // collect all the pixels in the neighborhood, note that we use
160  // GetPixel on the NeighborhoodIterator to honor the boundary conditions
161  typename OutputImageType::PixelType locVariation = 0;
162  typename ConstShapedNeighborhoodIterator<InputImageType>::ConstIterator input_neighbors_it;
163  for (input_neighbors_it = input_image_neighbors_it.Begin(); !input_neighbors_it.IsAtEnd(); input_neighbors_it++)
164  {
165  typename TInputImage::PixelType diffVec = input_neighbors_it.Get() - input_image_it.Get();
166  locVariation += SquaredEuclideanMetric<typename TInputImage::PixelType>::Calc(diffVec);
167  }
168  locVariation = sqrt(locVariation + 0.0001);
169  output_image_it.Set(locVariation);
170 
171  // update iterators
172  ++input_image_neighbors_it;
173  ++output_image_it;
174  ++input_image_it;
175 
176  // report progress
177  progress.CompletedPixel();
178  }
179  }
180  }
181 
182  /**
183  * Standard "PrintSelf" method
184  */
185  template <class TInputImage, class TOutput>
186  void LocalVariationImageFilter<TInputImage, TOutput>::PrintSelf(std::ostream &os, Indent indent) const
187  {
188  Superclass::PrintSelf(os, indent);
189  }
190 
191 } // end namespace itk
192 
193 #endif //_itkLocalVariationImageFilter_txx