+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Image__SkeletonFilter__hxx__
+#define __fpa__Image__SkeletonFilter__hxx__
+
+#include <itkMinimumMaximumImageCalculator.h>
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+SetInput( const TScalarImage* image )
+{
+ // Do nothing
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+SetInput( const TImage* image )
+{
+ this->m_DistanceMap->SetInput( image );
+ this->Superclass::SetInput( this->m_DistanceMap->GetOutput( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+_TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+GetInput( )
+{
+ return( this->m_DistanceMap->GetInput( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+const _TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+GetInput( ) const
+{
+ return( this->m_DistanceMap->GetInput( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+typename fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+TSkeleton* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+GetSkeleton( )
+{
+ return(
+ dynamic_cast< TSkeleton* >(
+ this->itk::ProcessObject::GetOutput( this->m_SkeletonIdx )
+ )
+ );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+typename fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+TMarks* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+GetMarks( )
+{
+ return(
+ dynamic_cast< TMarks* >(
+ this->itk::ProcessObject::GetOutput( this->m_MarksIdx )
+ )
+ );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+SkeletonFilter( )
+ : Superclass( ),
+ m_SeedFromMaximumDistance( false )
+{
+ this->m_DistanceMap = TDistanceMap::New( );
+ this->m_DistanceMap->InsideIsPositiveOn( );
+ this->m_DistanceMap->SquaredDistanceOff( );
+ this->m_DistanceMap->UseImageSpacingOn( );
+
+ unsigned int nOutputs = this->GetNumberOfRequiredOutputs( );
+ this->SetNumberOfRequiredOutputs( nOutputs + 2 );
+ this->m_SkeletonIdx = nOutputs;
+ this->m_MarksIdx = nOutputs + 1;
+
+ typename TSkeleton::Pointer skeleton = TSkeleton::New( );
+ typename TMarks::Pointer marks = TMarks::New( );
+ this->SetNthOutput( this->m_SkeletonIdx, skeleton.GetPointer( ) );
+ this->SetNthOutput( this->m_MarksIdx, marks.GetPointer( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+~SkeletonFilter( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+GenerateData( )
+{
+ // Prepare distance map
+ this->m_DistanceMap->Update( );
+
+ // Update start seed
+ TVertex seed;
+ if( this->m_SeedFromMaximumDistance )
+ {
+ typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
+ typename _TMinMax::Pointer minmax = _TMinMax::New( );
+ minmax->SetImage( this->m_DistanceMap->GetOutput( ) );
+ minmax->Compute( );
+ seed = minmax->GetIndexOfMaximum( );
+ }
+ else
+ seed = *( this->BeginSeeds( ) );
+ this->m_Seeds.clear( );
+ this->m_Seeds.insert( seed );
+
+ // Prepare skeleton candidates queue
+ this->m_SkeletonQueue.clear( );
+
+ // Go!
+ this->Superclass::GenerateData( );
+
+ // Backtracking
+ _TAdjacencies A;
+ std::vector< TVertex > end_points;
+ this->_EndPoints( end_points, A );
+ this->_Skeleton( end_points, A );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+_SetOutputValue( const TVertex& vertex, const TOutputValue& value )
+{
+ /* TODO
+ typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
+
+ this->Superclass::_UpdateResult( n );
+ double d = double( this->GetInput( )->GetPixel( n.Vertex ) );
+ if( d >= double( 0 ) )
+ {
+ // Update skeleton candidates
+ d += double( 1e-5 );
+ double v = double( n.Result ) / ( d * d );
+ this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) );
+
+ } // fi
+ */
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+_EndPoints( std::vector< TVertex >& end_points, _TAdjacencies& A )
+{
+ /* TODO
+ auto marks = this->GetMarks( );
+ auto mst = this->GetMinimumSpanningTree( );
+ auto spac = marks->GetSpacing( );
+
+ // Some values
+ marks->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
+ marks->SetRequestedRegion( mst->GetRequestedRegion( ) );
+ marks->SetBufferedRegion( mst->GetBufferedRegion( ) );
+ marks->SetSpacing( mst->GetSpacing( ) );
+ marks->SetOrigin( mst->GetOrigin( ) );
+ marks->SetDirection( mst->GetDirection( ) );
+ marks->Allocate( );
+ marks->FillBuffer( 0 );
+
+ // BFS from maximum queue
+ while( this->m_SkeletonQueue.size( ) > 0 )
+ {
+ // Get node
+ auto nIt = this->m_SkeletonQueue.begin( );
+ auto n = *nIt;
+ this->m_SkeletonQueue.erase( nIt );
+
+ // Mark it and update end-points
+ unsigned char m = marks->GetPixel( n.second );
+ if( m != 0 )
+ continue;
+ marks->SetPixel( n.second, 1 );
+ end_points.push_back( n.second );
+
+ // Mark path
+ TVertex it = n.second;
+ TVertex p = mst->GetParent( it );
+ while( it != p )
+ {
+ this->_MarkSphere( it );
+ it = p;
+ p = mst->GetParent( it );
+
+ } // elihw
+ this->_MarkSphere( it );
+ A[ n.second ] = it;
+
+ } // elihw
+ */
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+_Skeleton( const std::vector< TVertex >& end_points, _TAdjacencies& A )
+{
+ /* TODO
+ typedef itk::Image< unsigned long, Self::Dimension > _TTagsImage;
+ typedef typename TMST::TPath _TPath;
+
+ auto mst = this->GetMinimumSpanningTree( );
+ auto skeleton = this->GetSkeleton( );
+
+ // Tag branches
+ typename _TTagsImage::Pointer tags = _TTagsImage::New( );
+ tags->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
+ tags->SetRequestedRegion( mst->GetRequestedRegion( ) );
+ tags->SetBufferedRegion( mst->GetBufferedRegion( ) );
+ tags->Allocate( );
+ tags->FillBuffer( 0 );
+ for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
+ {
+ TVertex it = *eIt;
+ TVertex p = mst->GetParent( it );
+ while( it != p )
+ {
+ tags->SetPixel( it, tags->GetPixel( it ) + 1 );
+ it = p;
+ p = mst->GetParent( it );
+
+ } // elihw
+ tags->SetPixel( it, tags->GetPixel( it ) + 1 );
+
+ } // rof
+
+ // Build paths (branches)
+ for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
+ {
+ TVertex it = *eIt;
+ TVertex p = mst->GetParent( it );
+ TVertex sIdx = *eIt;
+ typename _TPath::Pointer path = _TPath::New( );
+ path->SetReferenceImage( mst );
+ while( it != p )
+ {
+ if( tags->GetPixel( sIdx ) != tags->GetPixel( it ) )
+ {
+ // Ok, a new branch has been added
+ path->AddVertex( it );
+ skeleton->AddBranch( path );
+
+ // Update a new starting index
+ path = _TPath::New( );
+ path->SetReferenceImage( mst );
+ sIdx = it;
+ }
+ else
+ path->AddVertex( it );
+ it = p;
+ p = mst->GetParent( it );
+
+ } // elihw
+
+ // Finally, add last branch
+ path->AddVertex( it );
+ skeleton->AddBranch( path );
+
+ } // rof
+ */
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TScalar >
+void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+_MarkSphere( const TVertex& idx )
+{
+ /* TODO
+ typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
+
+ static const double _eps = std::sqrt( double( Self::Dimension + 1 ) );
+ auto input = this->GetInput( );
+ auto marks = this->GetMarks( );
+ auto spac = input->GetSpacing( );
+ auto region = input->GetRequestedRegion( );
+
+ typename _TImage::PointType cnt;
+ input->TransformIndexToPhysicalPoint( idx, cnt );
+ double r = double( input->GetPixel( idx ) ) * _eps;
+
+ TVertex i0, i1;
+ for( unsigned int d = 0; d < Self::Dimension; ++d )
+ {
+ long off = long( std::ceil( r / double( spac[ d ] ) ) );
+ if( off < 3 )
+ off = 3;
+ i0[ d ] = idx[ d ] - off;
+ i1[ d ] = idx[ d ] + off;
+
+ if( i0[ d ] < region.GetIndex( )[ d ] )
+ i0[ d ] = region.GetIndex( )[ d ];
+
+ if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
+ i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
+
+ } // rof
+
+ typename _TImage::SizeType size;
+ for( unsigned int d = 0; d < Self::Dimension; ++d )
+ size[ d ] = i1[ d ] - i0[ d ] + 1;
+
+ typename _TImage::RegionType neighRegion;
+ neighRegion.SetIndex( i0 );
+ neighRegion.SetSize( size );
+
+ _TMarksIt mIt( marks, neighRegion );
+ for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
+ {
+ TVertex w = mIt.GetIndex( );
+ typename _TImage::PointType p;
+ marks->TransformIndexToPhysicalPoint( w, p );
+ mIt.Set( 1 );
+
+ } // rof
+ */
+}
+
+#endif // __fpa__Image__SkeletonFilter__hxx__
+
+// eof - $RCSfile$