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