lib/cpExtensions/Algorithms/DiscontinuityMapImageFilter.hxx

   1 // -------------------------------------------------------------------------
   2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
   3 // -------------------------------------------------------------------------
   4
   5 #ifndef __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__
   6 #define __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__
   7
   8 #include <cmath>
   9 #include <itkConstNeighborhoodIterator.h>
  10 #include <itkImageRegionIterator.h>
  11 #include <itkNeighborhoodAlgorithm.h>
  12 #include <itkZeroFluxNeumannBoundaryCondition.h>
  13
  14 // -------------------------------------------------------------------------
  15 template< class _TInputImage >
  16 cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
  17 DiscontinuityMapImageFilter( )
  18   : Superclass( ),
  19     m_UseSquareRoot( false )
  20 {
  21 }
  22
  23 // -------------------------------------------------------------------------
  24 template< class _TInputImage >
  25 cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
  26 ~DiscontinuityMapImageFilter( )
  27 {
  28 }
  29
  30 // -------------------------------------------------------------------------
  31 template< class _TInputImage >
  32 void cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
  33 ThreadedGenerateData( const TRegion& region, itk::ThreadIdType threadId )
  34 {
  35   unsigned int i;
  36
  37   itk::ZeroFluxNeumannBoundaryCondition< TInputImage > nbc;
  38   itk::ConstNeighborhoodIterator< TInputImage > bit;
  39   itk::ImageRegionIterator< TOutputImage > it;
  40
  41   // Allocate output
  42   typename TOutputImage::Pointer output = this->GetOutput( );
  43   typename TInputImage::ConstPointer input = this->GetInput( );
  44
  45   // Find the data-set boundary "faces"
  46   typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage >::FaceListType faceList;
  47   itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage > bC;
  48   faceList = bC( input, region, this->GetRadius( ) );
  49
  50   typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage >::FaceListType::iterator fit;
  51
  52   TOutputPixel sum1, sum2;
  53   TOutputPixel neighborhoodSize;
  54
  55   // Process each of the boundary faces.  These are N-d regions which border
  56   // the edge of the buffer.
  57   for( fit = faceList.begin( ); fit != faceList.end( ); ++fit )
  58   {
  59     bit = itk::ConstNeighborhoodIterator< TInputImage >( this->GetRadius( ), input, *fit );
  60     neighborhoodSize = bit.Size( );
  61     it = itk::ImageRegionIterator< TOutputImage >( output, *fit );
  62     bit.OverrideBoundaryCondition( &nbc );
  63     bit.GoToBegin( );
  64     while( !bit.IsAtEnd( ) )
  65     {
  66       sum1 = itk::NumericTraits< TOutputPixel >::ZeroValue( );
  67       sum2 = itk::NumericTraits< TOutputPixel >::ZeroValue( );
  68       for ( i = 0; i < neighborhoodSize; ++i )
  69       {
  70         TOutputPixel v = TOutputPixel( bit.GetPixel( i ) );
  71         sum1 += v;
  72         sum2 += v * v;
  73       }
  74       TOutputPixel mean = ( sum1 / TOutputPixel( neighborhoodSize ) );
  75       TOutputPixel var =
  76         ( sum2 - ( ( sum1 * sum1 ) / TOutputPixel( neighborhoodSize ) ) ) /
  77         ( TOutputPixel( neighborhoodSize ) - TOutputPixel( 1 ) );
  78       if( this->m_UseSquareRoot )
  79         it.Set( std::sqrt( double( var ) ) );
  80       else
  81         it.Set( var );
  82       ++bit;
  83       ++it;
  84
  85     } // elihw
  86
  87   } // rof
  88 }
  89
  90 #endif // __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__
  91
  92 // eof - $RCSfile$