#ifndef __CPM__DATASTRUCTURES__QUADEDGEMESH__HXX__
#define __CPM__DATASTRUCTURES__QUADEDGEMESH__HXX__

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
TPrimalEdge* cpm::DataStructures::QuadEdgeMesh< P, D, T >::
FindEdge( const PointIdentifier& a ) const
{
  if( a < this->m_OnextRings.size( ) )
    return( this->m_OnextRings[ a ] );
  else
    return( NULL );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
TPrimalEdge* cpm::DataStructures::QuadEdgeMesh< P, D, T >::
FindEdge( const PointIdentifier& a, const PointIdentifier& b ) const
{
  TPrimalEdge* found = NULL;
  TPrimalEdge* e = this->FindEdge( a );
  if( e != NULL )
  {
    typename TPrimalEdge::Iterator eIt = e->BeginOnext( );
    for( ; eIt != e->EndOnext( ) && found == NULL; eIt++ )
      if( ( *eIt )->GetDestination( ) == b )
        found = *eIt;

  } // fi
  return( found );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
TPrimalEdge* cpm::DataStructures::QuadEdgeMesh< P, D, T >::
FindEntryEdge( const CellIdentifier& i ) const
{
  if( i < this->GetNumberOfCells( ) )
  {
    CellAutoPointer c;
    this->GetCell( i, c );

    TQuadEdgeCell* ec = dynamic_cast< TQuadEdgeCell* >( c.GetPointer( ) );
    if( ec != NULL )
      return( ec->GetEntryPrimalEdge( ) );
    else
      return( NULL );
  }
  else
    return( NULL );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
const typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
CntNormals& cpm::DataStructures::QuadEdgeMesh< P, D, T >::
GetPointNormalsContainer( ) const
{
  return( this->m_PointNormals );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
NormalsIterator cpm::DataStructures::QuadEdgeMesh< P, D, T >::
BeginPointNormals( ) const
{
  return( this->m_PointNormals.begin( ) );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
NormalsIterator cpm::DataStructures::QuadEdgeMesh< P, D, T >::
EndPointNormals( ) const
{
  return( this->m_PointNormals.end( ) );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
const typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
VectorType& cpm::DataStructures::QuadEdgeMesh< P, D, T >::
GetPointNormal( const PointIdentifier& id ) const
{
  static const VectorType zero( TScalar( 0 ) );
  if( id < this->m_PointNormals.size( ) )
    return( this->m_PointNormals[ id ] );
  else
    return( zero );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
bool cpm::DataStructures::QuadEdgeMesh< P, D, T >::
RequestedRegionIsOutsideOfTheBufferedRegion( )
{
  // Mesh regions don't have meaning with QuadEdges, this is important
  // in order to guarantee correct pipeline execution
  return( false );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
SetPoints( PointsContainer* points )
{
  this->_ReleaseQuadEdgeObjects( );
  this->m_OnextRings.resize( points->Size( ), NULL );
  this->m_PointNormals.resize(
    points->Size( ), VectorType( TScalar( 0 ) )
    );
  this->Superclass::SetPoints( points );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
SetPoint( PointIdentifier id , PointType point )
{
  // If the point is brand new, add some non-existent edges and normals
  unsigned long N = this->GetNumberOfPoints( );
  N = ( N < id + 1 )? id + 1: N;
  if( this->m_OnextRings.size( ) < N )
  {
    this->m_OnextRings.resize( N, NULL );
    this->m_PointNormals.resize( N, VectorType( TScalar( 0 ) ) );

  } // fi

  // Ok pass it to itk
  this->Superclass::SetPoint( id, point );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
Initialize( )
{
  this->Superclass::Initialize( );
  this->_ReleaseQuadEdgeObjects( );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
Graft( const itk::DataObject* data )
{
  this->Superclass::Graft( data );
  const Self* mesh = dynamic_cast< const Self* >( data );
  if( mesh != NULL )
  {
    this->_ReleaseQuadEdgeObjects( );
    this->m_PrimalEdges    = mesh->m_PrimalEdges;
    this->m_DualEdges      = mesh->m_DualEdges;
    this->m_OnextRings     = mesh->m_OnextRings;
    this->m_PointNormals   = mesh->m_PointNormals;

  } // fi
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
SetCells( CellsContainer* cells )
{
  // Clear cells
  if( this->GetCells( ) != NULL )
    this->SetCells( CellsContainer::New( ) );
  this->_ReleaseQuadEdgeObjects( );
  this->m_OnextRings.resize( this->GetNumberOfPoints( ), NULL );

  // Copy cells
  if( cells != NULL )
  {
    for( unsigned long cId = 0; cId < cells->Size( ); ++cId )
    {
      CellAutoPointer cell( cells->GetElement( cId ), false );
      this->SetCell( cId, cell );

    } // rof

  } // fi
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
SetCellData( CellDataContainer* data )
{
  /* TODO
     No need for CellData in initial application:
     how should this be overridden?
  */
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
SetCell( CellIdentifier id, CellAutoPointer& ptr )
{
  // Overwrite an existing cell?
  if( id < this->GetNumberOfCells( ) )
  {
    itkExceptionMacro(
      << "TODO: What do I do here? Trying to overwrite a cell!"
      );

  } // fi

  // Do not allow to add vertices or edges as cells for the sake of
  // correct mesh construction
  if( ptr->GetNumberOfPoints( ) < 3 )
  {
    itkExceptionMacro(
      << "Face have less than 3 points ("
      << ptr->GetNumberOfPoints( )
      << ")"
      );

  } // fi

  // Fill up geometrical prerrequisites
  if( !( this->_CheckPoints( ptr ) ) )
  {
    itkExceptionMacro( << "Points missing in the mesh" );

  } // fi

  // Get (or create) individual edges
  _TEdges edges;
  this->_ConstructEdges( edges, ptr );

  // Configure new face's geometry
  edges[ 0 ]->SetLnextRingGeometry( id );

  // Update normals
  typename _TEdges::const_iterator eIt = edges.begin( );
  for( ; eIt != edges.end( ); ++eIt )
    this->m_PointNormals[ ( *eIt )->GetOrigin( ) ] =
      this->_ComputePointNormal( *eIt );

  // Save cell
  TQuadEdgeCell* qcell = new TQuadEdgeCell( edges[ 0 ] );
  CellAutoPointer cell;
  cell.TakeOwnership( qcell );
  this->Superclass::SetCell( id, cell );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
SetCellData( CellIdentifier id, CellPixelType data )
{
  /* TODO
     No need for CellData in initial application:
     how should this be overridden?
  */
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
BuildCellLinks( ) const
{
  // This is no longer required when using QuadEdges
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
cpm::DataStructures::QuadEdgeMesh< P, D, T >::
QuadEdgeMesh( )
  : Superclass( )
{
  this->_ReleaseQuadEdgeObjects( );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
cpm::DataStructures::QuadEdgeMesh< P, D, T >::
~QuadEdgeMesh( )
{
  this->_ReleaseQuadEdgeObjects( );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
bool cpm::DataStructures::QuadEdgeMesh< P, D, T >::
_CheckPoints( const CellAutoPointer& ptr ) const
{
  const typename Superclass::PointsContainer* pnts = this->GetPoints( );

  bool all_exist = true;
  typename CellType::PointIdConstIterator pIt = ptr->PointIdsBegin( );
  for( ; pIt != ptr->PointIdsEnd( ) && all_exist; ++pIt )
    all_exist &= pnts->IndexExists( *pIt );
  return( all_exist );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
_DeletePoint( const PointIdentifier& pId )
{
  PointIdentifier lastId = this->GetNumberOfPoints( ) - 1;
  if( lastId < pId )
    return;

  // Modify point container
  PointType to_delete = this->GetPoint( pId );
  PointType last = this->GetPoint( lastId );
  this->Superclass::SetPoint( pId, last );

  // Modify Onext ring and normals containers
  this->m_OnextRings[ pId ] = this->m_OnextRings[ lastId ];
  this->m_PointNormals[ pId ] = this->m_PointNormals[ lastId ];
  this->m_OnextRings.pop_back( );
  this->m_PointNormals.pop_back( );

  // Update origins
  typename TPrimalEdge::Iterator oIt =
    this->m_OnextRings[ pId ]->BeginOnext( );
  for( ; oIt != this->m_OnextRings[ pId ]->EndOnext( ); ++oIt )
    ( *oIt )->SetOrigin( pId );

  // Erase the point
  this->GetPoints( )->CastToSTLContainer( ).pop_back( );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
TPrimalEdge* cpm::DataStructures::QuadEdgeMesh< P, D, T >::
_CreateQuadEdge( const PointIdentifier& a, const PointIdentifier& b )
{
  // Create brand new object
  typename TPrimalEdge::Pointer edge = TPrimalEdge::New( );
  edge->CreateRings( );

  // Configure geometry
  edge->SetOrigin( a );
  edge->SetDestination( b );

  // Keep trace of all reserved memory
  this->m_PrimalEdges.insert( edge );
  this->m_PrimalEdges.insert( edge->GetSym( ) );
  this->m_DualEdges.insert( edge->GetRot( ) );
  this->m_DualEdges.insert( edge->GetInvRot( ) );

  return( edge );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
_DeleteEdge( TPrimalEdge* edge )
{
  this->m_DualEdges.erase( edge->GetInvRot( ) );
  this->m_PrimalEdges.erase( edge->GetSym( ) );
  this->m_DualEdges.erase( edge->GetRot( ) );
  this->m_PrimalEdges.erase( edge );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
_DeleteFace( const CellIdentifier& f )
{
  unsigned long cellId = this->GetNumberOfCells( ) - 1;
  if( cellId < f )
    return;

  // Exchange desired face with last one
  CellAutoPointer to_delete, last;
  this->GetCell( f, to_delete );
  this->GetCell( cellId, last );
  this->Superclass::SetCell( f, last );
  this->Superclass::SetCell( cellId, to_delete );

  // This will delete cell's memory when to_delete is no longer used
  to_delete.TakeOwnership( );

  // Squeeze container
  this->GetCells( )->CastToSTLContainer( ).pop_back( );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
_ConstructEdges( _TEdges& edges, const CellAutoPointer& ptr )
{
  edges.clear( );
  typename CellType::PointIdConstIterator pIt = ptr->PointIdsBegin( );
  for( pIt = ptr->PointIdsBegin( ); pIt != ptr->PointIdsEnd( ); ++pIt )
  {
    typename CellType::PointIdConstIterator nIt = pIt;
    ++nIt;
    if( nIt == ptr->PointIdsEnd( ) )
      nIt = ptr->PointIdsBegin( );

    TPrimalEdge* edge = this->FindEdge( *pIt, *nIt );
    if( edge == NULL )
    {
      edge = this->_CreateQuadEdge( *pIt, *nIt );
      if( this->m_OnextRings[ *pIt ] != NULL )
      {
        if( this->m_OnextRings[ *pIt ]->IsOriginInternal( ) )
        {
          itkExceptionMacro(
            << "Edge ["
            << *pIt << ", " << *nIt << "] is already internal."
            );

        } // fi

        // Put edge in Onext ring
        typename TPrimalEdge::Iterator eIt =
          this->m_OnextRings[ *pIt ]->BeginOnext( );
        while( eIt != this->m_OnextRings[ *pIt ]->EndOnext( ) )
        {
          if( !( ( *eIt )->IsLeftSet( ) ) )
          {
            TPrimalEdge::Splice( *eIt, edge );
            eIt = this->m_OnextRings[ *pIt ]->EndOnext( );
          }
          else
            ++eIt;

        } // elihw
      }
      else
        this->m_OnextRings[ *pIt ] = edge;
    }
    else
    {
      if( edge->IsLeftSet( ) )
      {
        itkExceptionMacro(
          << "Edge ["
          << edge->GetOrigin( ) << ", "
          << edge->GetDestination( ) << "] already have a left face ("
          << edge->GetLeft( )
          << "). Face already exits (at least in part)."
          );

      } // fi

    } // fi
    edges.push_back( edge );

  } // rof

  // Reorder Onext rings
  typename _TEdges::iterator ecIt;
  for( ecIt = edges.begin( ); ecIt != edges.end( ); ++ecIt )
  {
    TPrimalEdge* e = *ecIt;
    TPrimalEdge* eOprev;
    if( ecIt != edges.begin( ) )
    {
      typename _TEdges::iterator eOprevIt = ecIt;
      --eOprevIt;
      eOprev = ( *eOprevIt )->GetSym( );
    }
    else
      eOprev = edges.back( )->GetSym( );
    TPrimalEdge::ReorderOnextRing( e, eOprev );

  } // rof
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
typename cpm::DataStructures::QuadEdgeMesh< P, D, T >::
VectorType cpm::DataStructures::QuadEdgeMesh< P, D, T >::
_ComputePointNormal( const TPrimalEdge* e ) const
{
  static const TScalar zero = TScalar( 0 );

  VectorType n( zero );
  if( Superclass::PointDimension == 3 )
  {
    PointType p0 = this->GetPoint( e->GetOrigin( ) );
    typename TPrimalEdge::ConstIterator eIt = e->BeginOnext( );
    unsigned int count = 0;
    for( ; eIt != e->EndOnext( ); eIt++ )
    {
      if( ( *eIt )->IsLeftSet( ) )
      {
        typename TPrimalEdge::ConstIterator nIt = eIt;
        nIt++;
        if( nIt == e->EndOnext( ) )
          nIt = e->BeginOnext( );

        n += itk::CrossProduct(
          this->GetPoint( ( *eIt )->GetDestination( ) ) - p0,
          this->GetPoint( ( *nIt )->GetDestination( ) ) - p0
          );
        count++;

      } // fi

    } // rof
    TScalar nn = n.GetNorm( );
    if( nn > zero && count > 0 )
      n /= nn * TScalar( count );

  } // fi
  return( n );
}

// -------------------------------------------------------------------------
template< typename P, unsigned int D, typename T >
void cpm::DataStructures::QuadEdgeMesh< P, D, T >::
_ReleaseQuadEdgeObjects( )
{
  this->m_PrimalEdges.clear( );
  this->m_OnextRings.clear( );
  this->m_DualEdges.clear( );
  this->m_PointNormals.clear( );
}

#endif // __CPM__DATASTRUCTURES__QUADEDGEMESH__HXX__

// eof - $RCSfile$