+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
#ifndef __fpa__Image__SkeletonFilter__hxx__
#define __fpa__Image__SkeletonFilter__hxx__
-#include <fpa/Base/Functors/Inverse.h>
-#include <fpa/Image/Functors/SimpleNeighborhood.h>
+#include <itkImage.h>
#include <itkImageRegionIteratorWithIndex.h>
+#include <itkMinimumMaximumImageCalculator.h>
+#include <fpa/Image/Functors/Dijkstra/Invert.h>
// -------------------------------------------------------------------------
-#define fpa_Image_SkeletonFilter_OutputMacro( o_n, o_t ) \
- template< class _TImage > \
- typename fpa::Image::SkeletonFilter< _TImage >:: \
- o_t* fpa::Image::SkeletonFilter< _TImage >:: \
- Get##o_n( ) \
- { \
- return( \
- dynamic_cast< o_t* >( \
- this->itk::ProcessObject::GetOutput( this->m_##o_n##Idx ) \
- ) \
- ); \
- }
-
-fpa_Image_SkeletonFilter_OutputMacro( Skeleton, TSkeleton );
-fpa_Image_SkeletonFilter_OutputMacro( Marks, TMarks );
-
-// -------------------------------------------------------------------------
-template< class _TImage >
-fpa::Image::SkeletonFilter< _TImage >::
-SkeletonFilter( )
+template< class _TInputImage, class _TDistanceMap >
+fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra::
+_TDijkstra( )
: Superclass( )
{
- typedef typename _TImage::PixelType _TPixel;
- typedef fpa::Image::Functors::SimpleNeighborhood< _TImage::ImageDimension > _TNeighFunc;
- typedef fpa::Base::Functors::Inverse< _TPixel, _TPixel > _TInvFunc;
-
- 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( ) );
-
- typename _TNeighFunc::Pointer nfunc = _TNeighFunc::New( );
- nfunc->SetOrder( 2 );
- this->SetNeighborhoodFunction( nfunc );
-
- typename _TInvFunc::Pointer ifunc = _TInvFunc::New( );
- this->SetConversionFunction( ifunc );
+ // Prepare weight function
+ typedef fpa::Image::Functors::Dijkstra::Invert< TOutputImage, TScalar > _TWeight;
+ typename _TWeight::Pointer weight = _TWeight::New( );
+ weight->SetAlpha( 1 );
+ weight->SetBeta( 1 );
+ this->SetWeightFunction( weight );
}
// -------------------------------------------------------------------------
-template< class _TImage >
-fpa::Image::SkeletonFilter< _TImage >::
-~SkeletonFilter( )
+template< class _TInputImage, class _TDistanceMap >
+fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra::
+~_TDijkstra( )
{
}
// -------------------------------------------------------------------------
-template< class _TImage >
-void fpa::Image::SkeletonFilter< _TImage >::
+template< class _TInputImage, class _TDistanceMap >
+void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra::
_BeforeGenerateData( )
{
this->Superclass::_BeforeGenerateData( );
}
// -------------------------------------------------------------------------
-template< class _TImage >
-void fpa::Image::SkeletonFilter< _TImage >::
-_UpdateResult( const _TQueueNode& n )
+template< class _TInputImage, class _TDistanceMap >
+void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra::
+_UpdateOutputValue( const TNode& n )
{
typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
- this->Superclass::_UpdateResult( n );
+ this->Superclass::_UpdateOutputValue( 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 );
+ double v = double( n.Value ) / ( d * d );
this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) );
} // fi
}
// -------------------------------------------------------------------------
-template< class _TImage >
-void fpa::Image::SkeletonFilter< _TImage >::
-_AfterGenerateData( )
+template< class _TInputImage, class _TDistanceMap >
+itk::ModifiedTimeType
+fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+GetMTime( ) const
{
- this->Superclass::_AfterGenerateData( );
+ itk::ModifiedTimeType t = this->Superclass::GetMTime( );
+ itk::ModifiedTimeType q = this->m_DistanceMap->GetMTime( );
+ return( ( q < t )? q: t );
+}
- _TAdjacencies A;
- std::vector< TIndex > end_points;
- this->_EndPoints( end_points, A );
- this->_Skeleton( end_points, A );
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TDistanceMap >
+void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+SetInput( TInputImage* input )
+{
+ this->Superclass::SetNthInput( 0, input );
}
// -------------------------------------------------------------------------
-template< class _TImage >
-void fpa::Image::SkeletonFilter< _TImage >::
-_EndPoints( std::vector< TIndex >& end_points, _TAdjacencies& A )
+template< class _TInputImage, class _TDistanceMap >
+typename fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+TInputImage* fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+GetInput( )
{
- auto marks = this->GetMarks( );
- auto mst = this->GetMinimumSpanningTree( );
- auto spac = marks->GetSpacing( );
+ return( dynamic_cast< TInputImage* >( this->Superclass::GetInput( 0 ) ) );
+}
- // 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 );
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TDistanceMap >
+const typename fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+TInputImage* fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+GetInput( ) const
+{
+ return(
+ dynamic_cast< const TInputImage* >( this->Superclass::GetInput( 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 );
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TDistanceMap >
+typename fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+TSkeleton* fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+GetOutput( )
+{
+ return( dynamic_cast< TSkeleton* >( this->Superclass::GetOutput( 0 ) ) );
+}
- // 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
- TIndex it = n.second;
- TIndex p = mst->GetParent( it );
- while( it != p )
- {
- this->_MarkSphere( it );
- it = p;
- p = mst->GetParent( it );
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TDistanceMap >
+const typename fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+TSkeleton* fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+GetOutput( ) const
+{
+ return(
+ dynamic_cast< const TSkeleton* >( this->Superclass::GetOutput( 0 ) )
+ );
+}
- } // elihw
- this->_MarkSphere( it );
- A[ n.second ] = it;
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TDistanceMap >
+fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+SkeletonFilter( )
+ : Superclass( ),
+ m_SeedFromMaximumDistance( false )
+{
+ this->SetNumberOfRequiredInputs( 1 );
+ this->SetNumberOfRequiredOutputs( 1 );
+ this->SetNthOutput( 0, TSkeleton::New( ) );
+
+ this->m_DistanceMap = TDistanceMap::New( );
+}
- } // elihw
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TDistanceMap >
+fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+~SkeletonFilter( )
+{
}
// -------------------------------------------------------------------------
-template< class _TImage >
-void fpa::Image::SkeletonFilter< _TImage >::
-_Skeleton( const std::vector< TIndex >& end_points, _TAdjacencies& A )
+template< class _TInputImage, class _TDistanceMap >
+void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+GenerateData( )
{
- 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 )
+ // Update distance map
+ this->m_DistanceMap->SetInput( this->GetInput( ) );
+ this->m_DistanceMap->Update( );
+
+ // Correct seed
+ if( this->m_SeedFromMaximumDistance )
{
- TIndex it = *eIt;
- TIndex p = mst->GetParent( it );
- while( it != p )
- {
- tags->SetPixel( it, tags->GetPixel( it ) + 1 );
- it = p;
- p = mst->GetParent( it );
+ typedef itk::MinimumMaximumImageCalculator< TOutputImage > _TMinMax;
+ typename _TMinMax::Pointer minmax = _TMinMax::New( );
+ minmax->SetImage( this->m_DistanceMap->GetOutput( ) );
+ minmax->Compute( );
+ this->m_Seed = minmax->GetIndexOfMaximum( );
- } // elihw
- tags->SetPixel( it, tags->GetPixel( it ) + 1 );
+ } // fi
- } // rof
+ // Compute MST
+ typename _TDijkstra::Pointer dijkstra = _TDijkstra::New( );
+ dijkstra->SetInput( this->m_DistanceMap->GetOutput( ) );
+ dijkstra->AddSeed( this->m_Seed );
+ dijkstra->Update( );
- // Build paths (branches)
- for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
+ // Compute end-points
+ const _TMST* mst = dijkstra->GetMinimumSpanningTree( );
+ std::vector< TIndex > end_points;
+ this->_EndPoints(
+ end_points, this->m_DistanceMap->GetOutput( ),
+ mst, dijkstra->GetSkeletonQueue( )
+ );
+
+ // Create branches
+ TSkeleton* sk = this->GetOutput( );
+ typename std::vector< TIndex >::const_iterator eIt = end_points.begin( );
+ std::map< TIndex, unsigned long, typename TIndex::LexicographicCompare > tags;
+ unsigned long t = 1;
+ for( ; eIt != end_points.end( ); ++eIt )
{
+ // Tag path
TIndex it = *eIt;
TIndex p = mst->GetParent( it );
- TIndex 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 );
+ tags[ it ] = t;
it = p;
p = mst->GetParent( it );
} // elihw
-
- // Finally, add last branch
- path->AddVertex( it );
- skeleton->AddBranch( path );
+ // TODO: More on it
} // rof
}
// -------------------------------------------------------------------------
-template< class _TImage >
-void fpa::Image::SkeletonFilter< _TImage >::
-_MarkSphere( const TIndex& idx )
+template< class _TInputImage, class _TDistanceMap >
+template< class _TMarksPointer >
+void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+_MarkSphere(
+ _TMarksPointer& marks, const TOutputImage* dmap, const TIndex& center
+ )
{
- typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
+ typedef typename _TMarksPointer::ObjectType _TMarks;
+ 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 _TMarks::SpacingType spac = dmap->GetSpacing( );
+ typename _TMarks::RegionType region = dmap->GetRequestedRegion( );
- typename _TImage::PointType cnt;
- input->TransformIndexToPhysicalPoint( idx, cnt );
- double r = double( input->GetPixel( idx ) ) * _eps;
+ typename _TMarks::PointType cnt;
+ dmap->TransformIndexToPhysicalPoint( center, cnt );
+ double r = double( dmap->GetPixel( center ) ) * _eps;
TIndex 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;
+ i0[ d ] = center[ d ] - off;
+ i1[ d ] = center[ d ] + off;
if( i0[ d ] < region.GetIndex( )[ d ] )
i0[ d ] = region.GetIndex( )[ d ];
} // rof
- typename _TImage::SizeType size;
+ typename _TMarks::SizeType size;
for( unsigned int d = 0; d < Self::Dimension; ++d )
size[ d ] = i1[ d ] - i0[ d ] + 1;
- typename _TImage::RegionType neighRegion;
+ typename _TMarks::RegionType neighRegion;
neighRegion.SetIndex( i0 );
neighRegion.SetSize( size );
_TMarksIt mIt( marks, neighRegion );
for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
+ mIt.Set( true );
+}
+
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TDistanceMap >
+void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::
+_EndPoints(
+ std::vector< TIndex >& end_points,
+ const TOutputImage* dmap,
+ const _TMST* mst,
+ const _TSkeletonQueue& queue
+ )
+{
+ typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
+
+ // Some values
+ typedef itk::Image< bool, Self::Dimension > _TMarks;
+ typename _TMarks::Pointer marks = _TMarks::New( );
+ 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( false );
+
+ // Get candidates in maximum to minimum iteration
+ typename _TSkeletonQueue::const_reverse_iterator nIt = queue.rbegin( );
+ for( ; nIt != queue.rend( ); ++nIt )
{
- TIndex w = mIt.GetIndex( );
- typename _TImage::PointType p;
- marks->TransformIndexToPhysicalPoint( w, p );
- mIt.Set( 1 );
+ // Mark it and update end-points
+ if( !( marks->GetPixel( nIt->second ) ) )
+ {
+ marks->SetPixel( nIt->second, true );
+ end_points.push_back( nIt->second );
+
+ // Mark path
+ TIndex it = nIt->second;
+ TIndex p = mst->GetParent( it );
+ while( it != p )
+ {
+ this->_MarkSphere( marks, dmap, it );
+ it = p;
+ p = mst->GetParent( it );
+
+ } // elihw
+ this->_MarkSphere( marks, dmap, it );
+
+ } // fi
} // rof
}
+
+
+/* TODO
+ this->m_DistanceMap = TDistanceMap::New( );
+ this->m_DistanceMap->InsideIsPositiveOn( );
+ this->m_DistanceMap->SquaredDistanceOff( );
+ this->m_DistanceMap->UseImageSpacingOn( );
+
+ template< class _TImage, class _TScalar >
+ void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
+ _Skeleton( const std::vector< TVertex >& end_points, _TAdjacencies& A )
+ {
+ typedef typename TSkeleton::TPath _TPath;
+ typedef itk::Image< unsigned long, TImage::ImageDimension > _TTagsImage;
+
+ TMST* mst = this->GetMinimumSpanningTree( );
+ TSkeleton* 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 );
+ typename std::vector< TVertex >::const_iterator eIt = end_points.begin( );
+ for( ; 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)
+ eIt = end_points.begin( );
+ for( ; eIt != end_points.end( ); ++eIt )
+ {
+ TVertex it = *eIt;
+ TVertex p = mst->GetParent( it );
+ TVertex sIdx = it;
+ 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
+ }
+*/
+
#endif // __fpa__Image__SkeletonFilter__hxx__
// eof - $RCSfile$
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / blobdiff