#ifndef __CPM__ALGORITHMS__QUADEDGEMESHTO__DUALFILTER__HXX__
#define __CPM__ALGORITHMS__QUADEDGEMESHTO__DUALFILTER__HXX__

#include <queue>
#include <map>
#include <itkPolygonCell.h>

// -------------------------------------------------------------------------
template< class I, class O >
cpm::Algorithms::QuadEdge::MeshToDualFilter< I, O >::
MeshToDualFilter( )
  : Superclass( )
{
}

// -------------------------------------------------------------------------
template< class I, class O >
cpm::Algorithms::QuadEdge::MeshToDualFilter< I, O >::
~MeshToDualFilter( )
{
}

// -------------------------------------------------------------------------
template< class I, class O >
void cpm::Algorithms::QuadEdge::MeshToDualFilter< I, O >::
GenerateData( )
{
  typedef itk::PolygonCell< typename O::CellType > TPolygon;
  typedef typename I::TPrimalEdge                  TInPrimalEdge;
  typedef typename I::PixelType                    TInScalar;
  typedef typename O::PixelType                    TOutScalar;
  typedef typename I::PointType                    TInPoint;
  typedef typename O::PointType                    TOutPoint;
  typedef typename TInPoint::VectorType            TInVector;
  typedef std::map< TPointId, TPointId >          _TBorderColumn;
  typedef std::map< TPointId, _TBorderColumn >    _TBorderMatrix;

  const I* in = this->GetInput( );
  O* out = this->GetOutput( );
  // TODO: this->m_Function->SetInputMesh( in );

  // Convert faces to points
  unsigned long nCells = in->GetNumberOfCells( );
  out->GetPoints( )->Reserve( nCells );
  for( unsigned long cId = 0; cId < nCells; cId++ )
  {
    const typename I::TPrimalEdge* e = in->FindEntryEdge( cId );
    if( e == NULL )
      continue;
    if( !( e->IsLeftSet( ) ) )
      continue;

    // Compute new point
    TInVector bc( TInScalar( 0 ) );
    typename TInPrimalEdge::PointConstIterator pIt = e->BeginPointLnext( );
    unsigned long pCount = 0;
    for( ; pIt != e->EndPointLnext( ); pIt++ )
    {
      bc += in->GetPoint( *pIt ).GetVectorFromOrigin( );
      pCount++;

    } // rof
    if( pCount != 0 )
      bc /= TInScalar( pCount );

    TOutPoint pnt;
    for( unsigned int d = 0; d < I::PointDimension; d++ )
      pnt[ d ] = TOutScalar( bc[ d ] );
    out->SetPoint( e->GetLeft( ), pnt );

  } // rof

  _TBorderMatrix borderIds;
  unsigned long nPoints = in->GetNumberOfPoints( );
  for( unsigned long pId = 0; pId < nPoints; pId++ )
  {
    const TInPrimalEdge* e = in->FindEdge( pId );

    // Add dual faces: each point becomes a face
    typename O::CellAutoPointer cell;
    TPolygon* face = new TPolygon( );
    typename TInPrimalEdge::ConstIterator eIt = e->BeginOnext( );
    for( ; eIt != e->EndOnext( ); eIt++ )
    {
      if( ( *eIt )->IsAtBorder( ) )
      {
        // Keep track of created border points
        TPointId ori = ( *eIt )->GetOrigin( );
        TPointId des = ( *eIt )->GetDestination( );
        typename _TBorderMatrix::const_iterator bmIt =
          borderIds.find( ori );
        bool found = false;
        if( bmIt != borderIds.end( ) )
          found = ( bmIt->second.find( des ) != bmIt->second.end( ) );
        if( !found )
        {
          // Compute border point
          TInVector bc = in->GetPoint( ( *eIt )->GetOrigin( ) ).
            GetVectorFromOrigin( );
          bc += in->GetPoint( ( *eIt )->GetDestination( ) ).
            GetVectorFromOrigin( );
          bc /= TInScalar( 2 );

          TOutPoint pnt;
          for( unsigned int d = 0; d < I::PointDimension; d++ )
            pnt[ d ] = TOutScalar( bc[ d ] );

          unsigned long nOutPoints = out->GetNumberOfPoints( );
          out->SetPoint( nOutPoints, pnt );
          borderIds[ ori ][ des ] = nOutPoints;
          borderIds[ des ][ ori ] = nOutPoints;

        } // fi
        face->AddPointId( borderIds[ ori ][ des ] );

      } // fi
      if( ( *eIt )->IsLeftSet( ) )
        face->AddPointId( ( *eIt )->GetLeft( ) );

    } // rof
    cell.TakeOwnership( face );
    out->SetCell( out->GetNumberOfCells( ), cell );

  } // rof
}

#endif // __CPM__ALGORITHMS__QUADEDGEMESHTO__DUALFILTER__HXX__

// eof - $RCSfile$