// -------------------------------------------------------------------------
// @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 )
{
  this->m_Polygon->GetPoints( )->InsertNextPoint( x, y, 0 );
  this->Modified( );
}

// -------------------------------------------------------------------------
template< class _TImage >
void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
AddPoint( double p[ 2 ] )
{
  this->m_Polygon->GetPoints( )->InsertNextPoint( p[ 0 ], p[ 1 ], 0 );
  this->Modified( );
}

// -------------------------------------------------------------------------
template< class _TImage >
void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
AddPoints( vtkPoints* points )
{
  double p[ 3 ];
  for( unsigned long i = 0; i < points->GetNumberOfPoints( ); ++i )
  {
    points->GetPoint( i, p );
    this->m_Polygon->GetPoints( )->InsertNextPoint( p[ 0 ], p[ 1 ], 0 );

  } // rof
  this->Modified( );
}

// -------------------------------------------------------------------------
template< class _TImage >
void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
ClearPoints( )
{
  this->m_Polygon = vtkSmartPointer< vtkPolygon >::New( );
}

// -------------------------------------------------------------------------
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( )
{
}

// -------------------------------------------------------------------------
template< class _TImage >
void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
AfterThreadedGenerateData( )
{
}

// -------------------------------------------------------------------------
template< class _TImage >
void cpExtensions::Algorithms::RasterContourFilter< _TImage >::
ThreadedGenerateData( const TRegion& region, itk::ThreadIdType id )
{
  vtkPoints* points = this->m_Polygon->GetPoints( );
  unsigned long nPoints = points->GetNumberOfPoints( );
  double* arr =
    static_cast< double* >( points->GetData()->GetVoidPointer( 0 ) );
  double p[ 3 ], n[ 3 ], b[ 6 ];
  this->m_Polygon->GetPoints( )->GetBounds( b );
  n[ 0 ] = n[ 1 ] = double( 0 );
  n[ 2 ] = double( 1 );

  _TImage* out = this->GetOutput( );
  itk::ImageRegionIteratorWithIndex< _TImage > iIt( out, region );
  for( iIt.GoToBegin( ); !iIt.IsAtEnd( ); ++iIt )
  {
    TIndex idx = iIt.GetIndex( );
    TPoint pnt;
    out->TransformIndexToPhysicalPoint( idx, pnt );
    p[ 0 ] = pnt[ 0 ];
    p[ 1 ] = pnt[ 1 ];
    p[ 2 ] = double( 0 );

    int i, j, c = 0;
    int nvert = nPoints;
    for( i = 0, j = nvert - 1; i < nvert; j = i++ )
    {
      double pi[ 3 ], pj[ 3 ];
      points->GetPoint( i, pi );
      points->GetPoint( j, pj );

      if(
        ( ( pi[ 1 ] > p[ 1 ] ) != ( pj[ 1 ] > p[ 1 ] ) ) &&
        ( p[ 0 ] < ( pj[ 0 ] - pi[ 0 ] ) * ( p[ 1 ] - pi[ 1 ] ) / ( pj[ 1 ] - pi[ 1 ] ) + pi[ 0 ] )
        )
        c = !c;

    } // rof

    if( c != 0 )
      iIt.Set( this->m_InsideValue );
    else
      iIt.Set( this->m_OutsideValue );

  } // rof
}

#endif // __cpExtensions__Algorithms__RasterContourFilter__hxx__

// eof - $RCSfile$