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