#ifndef __fpa__Image__SkeletonFilter__hxx__
#define __fpa__Image__SkeletonFilter__hxx__

// -------------------------------------------------------------------------
#define __fpa__Image__SkeletonFilter__InputMacro( input_name, input_type, input_id ) \
  template< class _TDistanceMap, class _TCostMap >                      \
  typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::     \
  input_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::  \
  Get##input_name( )                                                    \
  {                                                                     \
    return(                                                             \
      dynamic_cast< input_type* >(                                      \
        this->Superclass::GetInput( input_id )                          \
        )                                                               \
      );                                                                \
  }                                                                     \
  template< class _TDistanceMap, class _TCostMap >                      \
  const                                                                 \
  typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::     \
  input_type*                                                           \
  fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::              \
  Get##input_name( ) const                                              \
  {                                                                     \
    return(                                                             \
      dynamic_cast< const input_type* >(                                \
        this->Superclass::GetInput( input_id )                          \
        )                                                               \
      );                                                                \
  }                                                                     \
  template< class _TDistanceMap, class _TCostMap >                      \
  void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::         \
  Set##input_name( input_type* input )                                  \
  {                                                                     \
    this->Superclass::SetInput( input_id, input );                      \
  }

__fpa__Image__SkeletonFilter__InputMacro( DistanceMap, TDistanceMap, 0 );
__fpa__Image__SkeletonFilter__InputMacro( CostMap, TCostMap, 1 );
__fpa__Image__SkeletonFilter__InputMacro( MinimumSpanningTree, TMST, 2 );

// -------------------------------------------------------------------------
#define __fpa__Image__SkeletonFilter__OutputMacro( output_name, output_type, output_id ) \
  template< class _TDistanceMap, class _TCostMap >                      \
  typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::     \
  output_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \
  Get##output_name( )                                                   \
  {                                                                     \
    return(                                                             \
      dynamic_cast< output_type* >(                                     \
        this->Superclass::GetOutput( output_id )                        \
        )                                                               \
      );                                                                \
  }                                                                     \
  template< class _TDistanceMap, class _TCostMap >                      \
  const typename                                                        \
  fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::              \
  output_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \
  Get##output_name( ) const                                             \
  {                                                                     \
    return(                                                             \
      dynamic_cast< const output_type* >(                               \
        this->Superclass::GetOutput( output_id )                        \
        )                                                               \
      );                                                                \
  }

__fpa__Image__SkeletonFilter__OutputMacro( Skeleton, TSkeleton, 0 );
__fpa__Image__SkeletonFilter__OutputMacro( EndPoints, TIndices, 1 );
__fpa__Image__SkeletonFilter__OutputMacro( Bifurcations, TIndices, 2 );

// -------------------------------------------------------------------------
template< class _TDistanceMap, class _TCostMap >
fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
SkeletonFilter( )
  : Superclass( )
{
  this->SetNumbeOfRequiredInputs( 3 );
  this->SetNumbeOfRequiredOutputs( 3 );

  typename TIndices::Pointer end_points = TIndices::New( );
  typename TIndices::Pointer bifurcations = TIndices::New( );
  typename TSkeleton::Pointer skeleton = TSkeleton::New( );
  this->SetNthOutput( 0, skeleton );
  this->SetNthOutput( 1, end_points );
  this->SetNthOutput( 2, bifurcations );
}

// -------------------------------------------------------------------------
template< class _TDistanceMap, class _TCostMap >
fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
~SkeletonFilter( )
{
}

// -------------------------------------------------------------------------
template< class _TDistanceMap, class _TCostMap >
void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
GenerateData( )
{
  // 0. I/O objects
  const TDistanceMap* dmap = this->GetDistanceMap( );
  const TCostMap* cmap = this->GetCostMap( );
  const TMST* mst = this->GetMinimumSpanningTree( );
  TIndices* end_points = this->GetEndPoints( );
  TIndices* bifurcations = this->GetBifurcations( );
  TSkeleton* skeleton = this->GetSkeleton( );

  // 1. Check input correspondance
  // TODO

  // 2. Detect end-points
  this->_EndPoints( dmap, cmap, mst, end_points );

  // 3. Build skeleton and keep track of bifurcations
  this->_Skeleton( dmap, cmap, mst, end_points, bifurcations, skeleton );
}

// -------------------------------------------------------------------------
template< class _TDistanceMap, class _TCostMap >
void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
_EndPoints(
  const TDistanceMap* dmap,
  const TCostMap* cmap,
  const TMST* mst,
  TIndices* end_points
  )
{
  typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt;
  typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt;
  typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
  typedef typename _TQueue::value_type _TQueueValue;
  typedef itk::Image< unsigned char, _TCostMap::ImageDimension > _TMarks;

  // Some values
  typename _TMarks::Pointer marks = _TMarks::New( );
  marks->SetLargestPossibleRegion( dmap->GetLargestPossibleRegion( ) );
  marks->SetRequestedRegion( dmap->GetRequestedRegion( ) );
  marks->SetBufferedRegion( dmap->GetBufferedRegion( ) );
  marks->SetSpacing( dmap->GetSpacing( ) );
  marks->SetOrigin( dmap->GetOrigin( ) );
  marks->SetDirection( dmap->GetDirection( ) );
  marks->Allocate( );
  marks->FillBuffer( 0 );

  // Create queue
  _TQueue queue;
  _TDistMapIt dIt( dmap, dmap->GetRequestedRegion( ) );
  _TCostMapIt cIt( cmap, cmap->GetRequestedRegion( ) );
  dIt.GoToBegin( );
  cIt.GoToBegin( );
  for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
  {
    double d = double( dIt.Get( ) );
    if( d > 0 )
    {
      double v = double( cIt.Get( ) ) / ( d * d );
      queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );

    } // fi

  } // rof

  // BFS from maximum queue
  while( queue.size( ) > 0 )
  {
    // Get node
    auto nIt = queue.begin( );
    auto n = *nIt;
    queue.erase( nIt );

    unsigned char m = marks->GetPixel( n.second );
    if( ( m & 0x01 ) == 0x01 )
      continue;

    // Mark it and update end-points
    m |= 0x01;
    marks->SetPixel( n.second, m );
    end_points->Get( ).insert( n.second );

    // Get path
    typename TMST::TPath::Pointer path;
    mst->GetPath( path, n.second );
    for( unsigned long i = 0; i < path->GetSize( ); ++i )
    {
      TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i );
      typename _TCostMap::PointType cnt;
      cmap->TransformIndexToPhysicalPoint( idx, cnt );
      double d = double( dmap->GetPixel( idx ) );
      d *= double( 2 );

      // Mark sphere
      std::queue< TIndex > q;
      q.push( idx );
      while( q.size( ) > 0 )
      {
        TIndex v = q.front( );
        q.pop( );
        unsigned char m = marks->GetPixel( v );
        if( ( m & 0x02 ) == 0x02 )
          continue;
        m |= 0x03;
        marks->SetPixel( v, m );

        for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x )
        {
          for( int y = -1; y <= 1; y += 2 )
          {
            TIndex w = v;
            w[ x ] += y;
            if( region.IsInside( w ) )
            {
              typename _TCostMap::PointType p;
              cmap->TransformIndexToPhysicalPoint( w, p );
              if( cnt.EuclideanDistanceTo( p ) <= d )
                q.push( w );

            } // fi

          } // rof

        } // rof

      } // elihw

    } // rof

  } // elihw
}

// -------------------------------------------------------------------------
template< class _TDistanceMap, class _TCostMap >
void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
_Skeleton(
  const TDistanceMap* dmap,
  const TCostMap* cmap,
  const TMST* mst,
  TIndices* end_points
  TIndices* bifurcations,
  TSkeleton* skeleton
  )
{
}

#endif // __fpa__Image__SkeletonFilter__hxx__

// eof - $RCSfile$