#ifndef __CPM__ALGORITHMS__QUADEDGE__MESHPLANECUTTERFILTER__HXX__
#define __CPM__ALGORITHMS__QUADEDGE__MESHPLANECUTTERFILTER__HXX__

#include <cmath>
#include <queue>
#include <set>

// -------------------------------------------------------------------------
template< class M >
void cpm::Algorithms::QuadEdge::MeshPlaneCutterFilter< M >::
SetPlaneNormal( const TVector& n )
{
  TScalar nn = n.GetNorm( );
  if( TScalar( 0 ) < nn )
  {
    this->m_PlaneNormal = n / nn;
    this->Modified( );

  } // fi
}

// -------------------------------------------------------------------------
template< class M >
cpm::Algorithms::QuadEdge::MeshPlaneCutterFilter< M >::
MeshPlaneCutterFilter( )
  : Superclass( )
{
  this->m_PlanePoint.Fill( TScalar( 0 ) );
  this->m_PlaneNormal.Fill( TScalar( 0 ) );
  this->m_PlaneNormal[ 0 ] = TScalar( 1 );
}

// -------------------------------------------------------------------------
template< class M >
cpm::Algorithms::QuadEdge::MeshPlaneCutterFilter< M >::
~MeshPlaneCutterFilter( )
{
}

// -------------------------------------------------------------------------
template< class M >
void cpm::Algorithms::QuadEdge::MeshPlaneCutterFilter< M >::
GenerateData( )
{
  typedef typename M::TPrimalEdge   _TEdge;
  typedef typename M::TQuadEdgeCell _TCell;
  typedef typename _TEdge::Iterator _TEdgeIt;

  const M* input = this->GetInput( );
  M* output = this->GetOutput( );

  // Clear output data
  output->Initialize( );

  // Marked edges
  std::set< _TEdge* > marks;

  unsigned long nCells = input->GetNumberOfCells( );
  unsigned long actualCell = 0;

  TVector pp = this->m_PlanePoint;
  TVector np = this->m_PlaneNormal;
  TVector pl, pm, pml;
  TScalar num, den;
  while( actualCell < nCells )
  {
    std::queue< _TEdge* > q;

    // Is the next cell a QuadEdgeCell?
    typename M::CellAutoPointer cell_ptr;
    input->GetCell( actualCell++, cell_ptr );
    _TCell* cell = dynamic_cast< _TCell* >( cell_ptr.GetPointer( ) );
    if( cell == NULL )
      continue;

    // Ok, it is... now, does its face needs to be analyzed?
    _TEdge* edge = cell->GetEntryPrimalEdge( );
    _TEdgeIt fIt = edge->BeginLnext( );
    for( ; fIt != edge->EndLnext( ); ++fIt )
    {
      // Has the edge already been visited?
      if( marks.find( *fIt ) != marks.end( ) )
        continue;

      // No, ok -> intersection candidate
      q.push( *fIt );

    } // rof

    // Main loop
    while( !( q.empty( ) ) )
    {
      // Pick next candidate
      edge = q.front( );
      q.pop( );

      // Check it
      if( marks.find( edge ) != marks.end( ) )
        continue;

      // Mark it (with its symmetric)
      marks.insert( edge );
      marks.insert( edge->GetSym( ) );

      // Now, edge is a true candidate: check intersection
      pl = input->GetPoint( edge->GetOrigin( ) ).GetVectorFromOrigin( );
      pm = input->GetPoint( edge->GetDestination( ) ).GetVectorFromOrigin( );
      pml = pm - pl;
      num = ( pp - pl ) * pml;
      den = pml * np;
      // std::cout << num << " " << den << " " << pp << " " << np << " " << pl << " " << pm << " {" << edge->GetOrigin( ) << "," << edge->GetDestination( ) << "}" << std::endl;
      if( !( double( 0 ) < std::fabs( double( den ) ) ) )
      {
        if( !( double( 0 ) < std::fabs( double( num ) ) ) )
        {
          // Line lies on plane: both points intersect
          // std::cout << pl << " " << pm << " ";

        } // fi
        // else: Line parallel to plane -> no intersection
      }
      else
      {
        num /= den;
        if( TScalar( 0 ) <= num && num <= TScalar( 1 ) )
        {
          // One point intersection
          // TPoint pnt( TScalar( 0 )
          output->SetPoint(
            output->GetNumberOfPoints( ),
            ( ( pm - pl ) * num ) + pl
            );

          // Check over face
          fIt = edge->BeginLnext( );
          for( ; fIt != edge->EndLnext( ); ++fIt )
          {
            // Has the edge already been visited?
            if( marks.find( *fIt ) != marks.end( ) )
              continue;

            // No, ok -> intersection candidate
            q.push( *fIt );

          } // rof
        }
        // else: no intersection between pm and pl

      } // fi

    } // elihw

  } // elihw
}

#endif // __CPM__ALGORITHMS__QUADEDGE__MESHPLANECUTTERFILTER__HXX__

// eof - $RCSfile$