1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
16 #ifndef _itkRegularizedIVIMLocalVariationImageFilter_txx
17 #define _itkRegularizedIVIMLocalVariationImageFilter_txx
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"
34 template <class TInputImage, class TOutputImage>
35 RegularizedIVIMLocalVariationImageFilter<TInputImage, TOutputImage>
36 ::RegularizedIVIMLocalVariationImageFilter()
39 template <class TInputImage, class TOutputImage>
41 RegularizedIVIMLocalVariationImageFilter<TInputImage, TOutputImage>
42 ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
44 // call the superclass' implementation of this method
45 Superclass::GenerateInputRequestedRegion();
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();
52 if ( !inputPtr || !outputPtr )
57 // get a copy of the input requested region (should equal the output
59 typename TInputImage::RegionType inputRequestedRegion;
60 inputRequestedRegion = inputPtr->GetRequestedRegion();
62 // pad the input requested region by 1
63 inputRequestedRegion.PadByRadius( 1 );
65 // crop the input requested region at the input's largest possible region
66 if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
68 inputPtr->SetRequestedRegion( inputRequestedRegion );
73 // Couldn't crop the region (requested region is outside the largest
74 // possible region). Throw an exception.
76 // store what we tried to request (prior to trying to crop)
77 inputPtr->SetRequestedRegion( inputRequestedRegion );
80 InvalidRequestedRegionError e(__FILE__, __LINE__);
81 e.SetLocation(ITK_LOCATION);
82 e.SetDescription("Requested region outside possible region.");
83 e.SetDataObject(inputPtr);
88 template< class TInputImage, class TOutputImage>
90 RegularizedIVIMLocalVariationImageFilter< TInputImage, TOutputImage>
91 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
92 ThreadIdType threadId)
96 typename OutputImageType::Pointer output = this->GetOutput();
97 typename InputImageType::ConstPointer input = this->GetInput();
99 itk::Size<InputImageDimension> size;
100 for( int i=0; i<InputImageDimension; i++)
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);
109 // support progress methods/callbacks
110 ProgressReporter progress(
111 this, threadId, outputRegionForThread.GetNumberOfPixels());
113 ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
114 std::vector<InputPixelType> pixels;
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)
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);
129 // neighborhood iterator for input image
130 ConstShapedNeighborhoodIterator<InputImageType>
131 input_image_neighbors_it(size, input, *fit);
132 typename ConstShapedNeighborhoodIterator<InputImageType>::
134 input_image_neighbors_it.OverrideBoundaryCondition(&nbc);
135 input_image_neighbors_it.ClearActiveList();
136 for(int i=0; i<InputImageDimension; i++)
140 input_image_neighbors_it.ActivateOffset(offset);
142 input_image_neighbors_it.ActivateOffset(offset);
144 input_image_neighbors_it.GoToBegin();
145 //const unsigned int neighborhoodSize = InputImageDimension*2;
147 while ( ! input_image_neighbors_it.IsAtEnd() )
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++)
158 typename TInputImage::PixelType diffVec =
159 input_neighbors_it.Get()-input_image_it.Get();
160 locVariation += IVIMSquaredEuclideanMetric
161 <typename TInputImage::PixelType>::Calc(diffVec);
163 locVariation = sqrt(locVariation + 0.0001);
164 output_image_it.Set(locVariation);
167 ++input_image_neighbors_it;
172 progress.CompletedPixel();
178 * Standard "PrintSelf" method
180 template <class TInputImage, class TOutput>
182 RegularizedIVIMLocalVariationImageFilter<TInputImage, TOutput>
187 Superclass::PrintSelf( os, indent );
190 } // end namespace itk
192 #endif //_itkRegularizedIVIMLocalVariationImageFilter_txx