#ifndef __FPA__BASE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
#define __FPA__BASE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__

#include <cmath>
#include <map>
#include <vector>

// -------------------------------------------------------------------------
template< class _TMST >
const typename
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
TMinimumSpanningTree*
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
GetMinimumSpanningTree( )
{
  return(
    dynamic_cast< const _TMST* >( this->itk::ProcessObject::GetInput( 0 ) )
    );
}

// -------------------------------------------------------------------------
template< class _TMST >
void
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
SetMinimumSpanningTree( TMinimumSpanningTree* mst )
{
  this->itk::ProcessObject::SetNthInput( 0, const_cast< _TMST* >( mst ) );
}

// -------------------------------------------------------------------------
template< class _TMST >
typename
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
TVertices*
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
GetEndPoints( )
{
  return(
    dynamic_cast< TVertices* >( this->itk::ProcessObject::GetOutput( 0 ) )
    );
}

// -------------------------------------------------------------------------
template< class _TMST >
typename
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
TVertices*
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
GetBifurcations( )
{
  return(
    dynamic_cast< TVertices* >( this->itk::ProcessObject::GetOutput( 1 ) )
    );
}

// -------------------------------------------------------------------------
template< class _TMST >
typename
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
TVertices*
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
GetCollisions( )
{
  return(
    dynamic_cast< TVertices* >( this->itk::ProcessObject::GetOutput( 2 ) )
    );
}

// -------------------------------------------------------------------------
template< class _TMST >
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
  : Superclass( )
{
  this->SetNumberOfRequiredInputs( 1 );
  this->SetNumberOfRequiredOutputs( 3 );
  typename TVertices::Pointer ep = TVertices::New( );
  typename TVertices::Pointer bf = TVertices::New( );
  typename TVertices::Pointer co = TVertices::New( );
  this->itk::ProcessObject::SetNthOutput( 0, ep.GetPointer( ) );
  this->itk::ProcessObject::SetNthOutput( 1, bf.GetPointer( ) );
  this->itk::ProcessObject::SetNthOutput( 2, co.GetPointer( ) );
}

// -------------------------------------------------------------------------
template< class _TMST >
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
~ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
{
}

// -------------------------------------------------------------------------
template< class _TMST >
void
fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
GenerateData( )
{
  typedef std::multimap< double, TVertex > _TQueue;
  typedef std::pair< TVertex, TVertex >    _TBranch;
  typedef std::vector< _TBranch >          _TBranches;
  typedef typename _TQueue::value_type     _TQueueValue;

  // 0. Useful objects and values
  auto mst = this->GetMinimumSpanningTree( );
  auto endpoints = this->GetEndPoints( );
  auto bifurcations = this->GetBifurcations( );
  auto collisions = this->GetCollisions( );
  _TBranches branches;
  endpoints->Get( ).clear( );
  bifurcations->Get( ).clear( );
  collisions->Get( ).clear( );

  // 1. Get priority queue
  auto& q = mst->GetNodeQueue( );

  // 2. Main loop: iterate in inverse order (max-priority)
  unsigned long label = 0;
  for( auto qIt = q.rbegin( ); qIt != q.rend( ); ++qIt )
  {
    auto gCost = qIt->first;
    auto vertex = qIt->second;

    // 2.1. Check if the vertex has already been visited
    if( this->_Mark( vertex ) > 0 )
      continue;

    // 2.2. Get path from front seed
    auto path = mst->GetPath( vertex );

    // 2.3. Prepare new branch data and prepare new end-point
    label++;
    branches.push_back( _TBranch( vertex, vertex ) );
    endpoints->Get( ).push_back( vertex );

    // 2.4. Backtracking
    auto pIt = path.begin( );
    TVertex last_collision;
    bool inCollision = false;
    while( pIt != path.end( ) )
    {
      auto mark = this->_SkeletonMark( *pIt );
      if( mark > 0 && mark != label )
      {
        // 2.4.1. A bifurcation has been reached
        bifurcations->Get( ).push_back( *pIt );

        // Reorder labels
        auto coll_branch = branches[ mark ];
        branches[ mark  ] = _TBranch( coll_branch.first, *pIt );
        branches[ label - 1 ] = _TBranch( qIt->second, *pIt );
        branches.push_back( _TBranch( *pIt, coll_branch.second ) );

        // Mark skeleton (b,coll_branch_second) with the new label
        label++;
        while( *pIt != coll_branch.second && pIt != path.end( ) )
        {
          this->_MarkSkeleton( *pIt, label );
          pIt++;

        } // elihw

        // Force inner loop termination
        pIt = path.end( );
      }
      else
      {
        if( !inCollision )
        {
          mark = this->_Mark( *pIt );
          if( mark > 0 && mark != label )
          {
            // 2.4.2. A collision has been reached
            last_collision = *pIt;
            collisions->Get( ).push_back( *pIt );
            inCollision = true;
          }
          else
          {
            // 2.4.3. Just mark sphere and skeleton
            double r = this->_Radius( *pIt );
            this->_MarkSkeleton( *pIt, label );
            this->_MarkSphere( *pIt, r, label );
            branches[ label - 1 ].second = *pIt;

            // 2.4.4. Is this a seed? -> add it to endpoints
            auto sIt = pIt;
            sIt++;
            if( sIt == path.end( ) )
              endpoints->Get( ).push_back( *pIt );

          } // fi

        } // fi
        pIt++;

      } // fi

    } // elihw

  } // rof

  endpoints->SetReferenceImage( mst );
  bifurcations->SetReferenceImage( mst );
  collisions->SetReferenceImage( mst );
}

#endif // __FPA__BASE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__

// eof - $RCSfile$