// -------------------------------------------------------------------------
// @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$