1 // -------------------------------------------------------------------------
2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
3 // -------------------------------------------------------------------------
5 #ifndef __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__
6 #define __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__
9 #include <itkConstNeighborhoodIterator.h>
10 #include <itkImageRegionIterator.h>
11 #include <itkNeighborhoodAlgorithm.h>
12 #include <itkZeroFluxNeumannBoundaryCondition.h>
14 // -------------------------------------------------------------------------
15 template< class _TInputImage >
16 cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
17 DiscontinuityMapImageFilter( )
19 m_UseSquareRoot( false )
23 // -------------------------------------------------------------------------
24 template< class _TInputImage >
25 cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
26 ~DiscontinuityMapImageFilter( )
30 // -------------------------------------------------------------------------
31 template< class _TInputImage >
32 void cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
33 ThreadedGenerateData( const TRegion& region, itk::ThreadIdType threadId )
37 itk::ZeroFluxNeumannBoundaryCondition< TInputImage > nbc;
38 itk::ConstNeighborhoodIterator< TInputImage > bit;
39 itk::ImageRegionIterator< TOutputImage > it;
42 typename TOutputImage::Pointer output = this->GetOutput( );
43 typename TInputImage::ConstPointer input = this->GetInput( );
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( ) );
50 typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage >::FaceListType::iterator fit;
52 TOutputPixel sum1, sum2;
53 TOutputPixel neighborhoodSize;
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 )
59 bit = itk::ConstNeighborhoodIterator< TInputImage >( this->GetRadius( ), input, *fit );
60 neighborhoodSize = TOutputPixel( bit.Size( ) );
61 it = itk::ImageRegionIterator< TOutputImage >( output, *fit );
62 bit.OverrideBoundaryCondition( &nbc );
64 while( !bit.IsAtEnd( ) )
66 sum1 = itk::NumericTraits< TOutputPixel >::ZeroValue( );
67 sum2 = itk::NumericTraits< TOutputPixel >::ZeroValue( );
68 for ( i = 0; i < neighborhoodSize; ++i )
70 TOutputPixel v = TOutputPixel( bit.GetPixel( i ) );
75 ( sum2 - ( ( sum1 * sum1 ) / neighborhoodSize ) ) /
76 ( neighborhoodSize - TOutputPixel( 1 ) );
77 if( this->m_UseSquareRoot )
78 it.Set( std::sqrt( double( var ) ) );
89 #endif // __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__