#ifndef __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
#define __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__

#include <cmath>
#include <itkImageRegionIteratorWithIndex.h>

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
const typename
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
TImage*
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
GetCostsImage( )
{
  return(
    dynamic_cast< const _TImage* >( this->itk::ProcessObject::GetInput( 1 ) )
    );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
const typename
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
TImage*
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
GetDistanceMap( )
{
  return(
    dynamic_cast< const _TImage* >( this->itk::ProcessObject::GetInput( 2 ) )
    );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
void
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
SetCostsImage( TImage* image )
{
  this->itk::ProcessObject::SetNthInput( 1, const_cast< _TImage* >( image ) );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
void
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
SetDistanceMap( TImage* image )
{
  this->itk::ProcessObject::SetNthInput( 2, const_cast< _TImage* >( image ) );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
  : Superclass( ),
    m_SquaredDistanceMap( false )
{
  this->SetNumberOfRequiredInputs( 3 );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
~ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
{
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
void
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
GenerateData( )
{
  auto image = this->GetCostsImage( );

  this->m_MarkImage = TMarkImage::New( );
  this->m_MarkImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
  this->m_MarkImage->SetRequestedRegion( image->GetRequestedRegion( ) );
  this->m_MarkImage->SetBufferedRegion( image->GetBufferedRegion( ) );
  this->m_MarkImage->SetDirection( image->GetDirection( ) );
  this->m_MarkImage->SetOrigin( image->GetOrigin( ) );
  this->m_MarkImage->SetSpacing( image->GetSpacing( ) );
  this->m_MarkImage->Allocate( );
  this->m_MarkImage->FillBuffer( 0 );

  this->m_SkeletonImage = TMarkImage::New( );
  this->m_SkeletonImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
  this->m_SkeletonImage->SetRequestedRegion( image->GetRequestedRegion( ) );
  this->m_SkeletonImage->SetBufferedRegion( image->GetBufferedRegion( ) );
  this->m_SkeletonImage->SetDirection( image->GetDirection( ) );
  this->m_SkeletonImage->SetOrigin( image->GetOrigin( ) );
  this->m_SkeletonImage->SetSpacing( image->GetSpacing( ) );
  this->m_SkeletonImage->Allocate( );
  this->m_SkeletonImage->FillBuffer( 0 );

  this->Superclass::GenerateData( );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
void
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
_MarkSkeleton( const TVertex& v, const unsigned long& l )
{
  if( this->m_SkeletonImage->GetRequestedRegion( ).IsInside( v ) )
    this->m_SkeletonImage->SetPixel( v, l );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
void
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
_MarkSphere(
  const TVertex& v,
  const double& r,
  const unsigned long& l
  )
{
  double rr = r * double( 1.5 );

  // Get marking region
  auto rreg = this->m_MarkImage->GetRequestedRegion( );
  auto spac = this->m_MarkImage->GetSpacing( );
  TVertex ireg = rreg.GetIndex( );
  TVertex jreg = ireg + rreg.GetSize( );
  TVertex idx, jdx;
  typename _TImage::SizeType size;
  bool all = false;
  for( unsigned int d = 0; d < _TImage::ImageDimension; ++d )
  {
    unsigned long s = std::ceil( rr / double( spac[ d ] ) );
    if( s < 3 )
    {
      s = 3;
      all = true;

    } // fi
    s += 3;
    idx[ d ] = v[ d ] - s;
    jdx[ d ] = v[ d ] + s;

    jreg[ d ]--;
    if( idx[ d ] < ireg[ d ] ) idx[ d ] = ireg[ d ];
    if( idx[ d ] > jreg[ d ] ) idx[ d ] = jreg[ d ];
    if( jdx[ d ] < ireg[ d ] ) jdx[ d ] = ireg[ d ];
    if( jdx[ d ] > jreg[ d ] ) jdx[ d ] = jreg[ d ];
    size[ d ] = jdx[ d ] - idx[ d ] + 1;

  } // rof
  typename _TImage::RegionType region;
  region.SetIndex( idx );
  region.SetSize( size );

  // Mark region
  typename _TImage::PointType c;
  this->m_MarkImage->TransformIndexToPhysicalPoint( v, c );
  itk::ImageRegionIteratorWithIndex< TMarkImage >
    spIt( this->m_MarkImage, region );
  for( spIt.GoToBegin( ); !spIt.IsAtEnd( ); ++spIt )
  {
    typename _TImage::PointType pnt;
    this->m_MarkImage->TransformIndexToPhysicalPoint( spIt.GetIndex( ), pnt );
    if( double( pnt.EuclideanDistanceTo( c ) ) <= rr || all )
      if( spIt.Get( ) == 0 )
        spIt.Set( l );

  } // rof
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
unsigned long
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
_Mark( const TVertex& v )
{
  return( this->m_MarkImage->GetPixel( v ) );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
unsigned long fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
_SkeletonMark( const TVertex& v )
{
  return( this->m_SkeletonImage->GetPixel( v ) );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TMST >
double
fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
_Radius( const TVertex& v )
{
  double d = double( this->GetDistanceMap( )->GetPixel( v ) );
  double s = std::fabs( d );
  if( s > double( 0 ) )
    s /= d;
  if( this->m_SquaredDistanceMap )
    return( s * std::sqrt( std::fabs( d ) ) );
  else
    return( d );
}

#endif // __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__

// eof - $RCSfile$