Discontinuity map (Wan Yu) filter added.

[cpPlugins.git] / lib / cpExtensions / Algorithms / DiscontinuityMapImageFilter.hxx

diff --git a/lib/cpExtensions/Algorithms/DiscontinuityMapImageFilter.hxx b/lib/cpExtensions/Algorithms/DiscontinuityMapImageFilter.hxx
new file mode 100644 (file)
index 0000000..877edbc
--- /dev/null
+++ b/lib/cpExtensions/Algorithms/DiscontinuityMapImageFilter.hxx
@@ -0,0 +1,92 @@
+// -------------------------------------------------------------------------
+// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+// -------------------------------------------------------------------------
+
+#ifndef __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__
+#define __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__
+
+#include <cmath>
+#include <itkConstNeighborhoodIterator.h>
+#include <itkImageRegionIterator.h>
+#include <itkNeighborhoodAlgorithm.h>
+#include <itkZeroFluxNeumannBoundaryCondition.h>
+
+// -------------------------------------------------------------------------
+template< class _TInputImage >
+cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
+DiscontinuityMapImageFilter( )
+  : Superclass( ),
+    m_UseSquareRoot( false )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class _TInputImage >
+cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
+~DiscontinuityMapImageFilter( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class _TInputImage >
+void cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
+ThreadedGenerateData( const TRegion& region, itk::ThreadIdType threadId )
+{
+  unsigned int i;
+
+  itk::ZeroFluxNeumannBoundaryCondition< TInputImage > nbc;
+  itk::ConstNeighborhoodIterator< TInputImage > bit;
+  itk::ImageRegionIterator< TOutputImage > it;
+
+  // Allocate output
+  typename TOutputImage::Pointer output = this->GetOutput( );
+  typename TInputImage::ConstPointer input = this->GetInput( );
+
+  // Find the data-set boundary "faces"
+  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage >::FaceListType faceList;
+  itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage > bC;
+  faceList = bC( input, region, this->GetRadius( ) );
+
+  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage >::FaceListType::iterator fit;
+
+  TOutputPixel sum1, sum2;
+  TOutputPixel neighborhoodSize;
+
+  // Process each of the boundary faces.  These are N-d regions which border
+  // the edge of the buffer.
+  for( fit = faceList.begin( ); fit != faceList.end( ); ++fit )
+  {
+    bit = itk::ConstNeighborhoodIterator< TInputImage >( this->GetRadius( ), input, *fit );
+    neighborhoodSize = bit.Size( );
+    it = itk::ImageRegionIterator< TOutputImage >( output, *fit );
+    bit.OverrideBoundaryCondition( &nbc );
+    bit.GoToBegin( );
+    while( !bit.IsAtEnd( ) )
+    {
+      sum1 = itk::NumericTraits< TOutputPixel >::ZeroValue( );
+      sum2 = itk::NumericTraits< TOutputPixel >::ZeroValue( );
+      for ( i = 0; i < neighborhoodSize; ++i )
+      {
+        TOutputPixel v = TOutputPixel( bit.GetPixel( i ) );
+        sum1 += v;
+        sum2 += v * v;
+      }
+      TOutputPixel mean = ( sum1 / TOutputPixel( neighborhoodSize ) );
+      TOutputPixel var =
+        ( sum2 - ( ( sum1 * sum1 ) / TOutputPixel( neighborhoodSize ) ) ) /
+        ( TOutputPixel( neighborhoodSize ) - TOutputPixel( 1 ) );
+      if( this->m_UseSquareRoot )
+        it.Set( std::sqrt( double( var ) ) );
+      else
+        it.Set( var );
+      ++bit;
+      ++it;
+
+    } // elihw
+
+  } // rof
+}
+
+#endif // __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__
+
+// eof - $RCSfile$