--- /dev/null
+// -------------------------------------------------------------------------
+// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+// ------------------------------------------------------------------------
+
+// Inclusion test taken from:
+// https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
+
+#ifndef __cpExtensions__Algorithms__RasterContourFilter__hxx__
+#define __cpExtensions__Algorithms__RasterContourFilter__hxx__
+
+#include <itkImageRegionIteratorWithIndex.h>
+#include <vtkPolygon.h>
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
+AddPoint( double x, double y )
+{
+ TPoint pnt;
+ pnt[ 0 ] = x;
+ pnt[ 1 ] = y;
+ this->m_Contour.push_back( pnt );
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
+AddPoint( double p[ 2 ] )
+{
+ TPoint pnt;
+ pnt[ 0 ] = p[ 0 ];
+ pnt[ 1 ] = p[ 2 ];
+ this->m_Contour.push_back( pnt );
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
+ClearPoints( )
+{
+ this->m_Contour.clear( );
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+cpExtensions::Algorithms::RasterContourFilter< _TImage >::
+RasterContourFilter( )
+ : Superclass( ),
+ m_InsideValue( TPixel( 1 ) ),
+ m_OutsideValue( TPixel( 0 ) )
+{
+ this->ClearPoints( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+cpExtensions::Algorithms::RasterContourFilter< _TImage >::
+~RasterContourFilter( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
+AllocateOutputs( )
+{
+ _TImage* out = this->GetOutput( 0 );
+ out->SetSpacing( this->m_Template->GetSpacing( ) );
+ out->SetRegions( this->m_Template->GetRequestedRegion( ) );
+ out->SetOrigin( this->m_Template->GetOrigin( ) );
+ out->SetDirection( this->m_Template->GetDirection( ) );
+ out->Allocate( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
+BeforeThreadedGenerateData( )
+{
+ // Keep just indices, not points
+ this->m_Polygon.clear( );
+ _TImage* out = this->GetOutput( 0 );
+ TIndex minIdx, maxIdx;
+ for( auto c = this->m_Contour.begin( ); c != this->m_Contour.end( ); ++c )
+ {
+ TIndex idx;
+ out->TransformPhysicalPointToIndex( *c, idx );
+ bool added = true;
+ if( this->m_Polygon.size( ) > 0 )
+ {
+ if( this->m_Polygon.back( ) != idx )
+ {
+ this->m_Polygon.push_back( idx );
+ minIdx[ 0 ] = ( idx[ 0 ] < minIdx[ 0 ] )? idx[ 0 ]: minIdx[ 0 ];
+ minIdx[ 1 ] = ( idx[ 1 ] < minIdx[ 1 ] )? idx[ 1 ]: minIdx[ 1 ];
+ maxIdx[ 0 ] = ( idx[ 0 ] > maxIdx[ 0 ] )? idx[ 0 ]: maxIdx[ 0 ];
+ maxIdx[ 1 ] = ( idx[ 1 ] > maxIdx[ 1 ] )? idx[ 1 ]: maxIdx[ 1 ];
+
+ } // fi
+ }
+ else
+ {
+ this->m_Polygon.push_back( idx );
+ minIdx = maxIdx = idx;
+
+ } // fi
+
+ } // rof
+
+ // Set ROI
+ typename _TImage::SizeType size;
+ size[ 0 ] = maxIdx[ 0 ] - minIdx[ 0 ] + 1;
+ size[ 1 ] = maxIdx[ 1 ] - minIdx[ 1 ] + 1;
+ this->m_ROI.SetIndex( minIdx );
+ this->m_ROI.SetSize( size );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
+AfterThreadedGenerateData( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
+ThreadedGenerateData( const TRegion& region, itk::ThreadIdType id )
+{
+ long nVerts = this->m_Polygon.size( );
+ _TImage* out = this->GetOutput( );
+ itk::ImageRegionIteratorWithIndex< _TImage > iIt( out, region );
+ for( iIt.GoToBegin( ); !iIt.IsAtEnd( ); ++iIt )
+ {
+ TIndex p = iIt.GetIndex( );
+ bool inside = false;
+ if( this->m_ROI.IsInside( p ) )
+ {
+ long i, j;
+ for( i = 0, j = nVerts - 1; i < nVerts; j = i++ )
+ {
+ TIndex pi = this->m_Polygon[ i ];
+ TIndex pj = this->m_Polygon[ j ];
+ double pi0 = double( pi[ 0 ] );
+ double pi1 = double( pi[ 1 ] );
+ double pj0 = double( pj[ 0 ] );
+ double pj1 = double( pj[ 1 ] );
+ double p0 = double( p[ 0 ] );
+ double p1 = double( p[ 1 ] );
+ double ji0 = pj0 - pi0;
+ double ji1 = pj1 - pi1;
+ double i1 = p1 - pi1;
+ if(
+ ( ( pi1 > p1 ) != ( pj1 > p1 ) ) &&
+ ( p0 < ( ( ji0 * i1 ) / ji1 ) + pi0 )
+ )
+ inside = !inside;
+
+ } // rof
+
+ } // fi
+ iIt.Set( ( inside )? this->m_InsideValue: this->m_OutsideValue );
+
+ } // rof
+}
+
+#endif // __cpExtensions__Algorithms__RasterContourFilter__hxx__
+
+// eof - $RCSfile$