#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 );
+#include <fpa/Base/Functors/Inverse.h>
+#include <fpa/Image/Functors/SimpleNeighborhood.h>
+#include <itkImageRegionIteratorWithIndex.h>
// -------------------------------------------------------------------------
-#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 \
+#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< const output_type* >( \
- this->Superclass::GetOutput( output_id ) \
+ dynamic_cast< o_t* >( \
+ this->itk::ProcessObject::GetOutput( this->m_##o_n##Idx ) \
) \
); \
}
-__fpa__Image__SkeletonFilter__OutputMacro( Skeleton, TSkeleton, 0 );
-__fpa__Image__SkeletonFilter__OutputMacro( EndPoints, TIndices, 1 );
-__fpa__Image__SkeletonFilter__OutputMacro( Bifurcations, TIndices, 2 );
+fpa_Image_SkeletonFilter_OutputMacro( Skeleton, TSkeleton );
+fpa_Image_SkeletonFilter_OutputMacro( Marks, TMarks );
// -------------------------------------------------------------------------
-template< class _TDistanceMap, class _TCostMap >
-fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
+template< class _TImage >
+fpa::Image::SkeletonFilter< _TImage >::
SkeletonFilter( )
: Superclass( )
{
- this->SetNumbeOfRequiredInputs( 3 );
- this->SetNumbeOfRequiredOutputs( 3 );
+ typedef typename _TImage::PixelType _TPixel;
+ typedef typename _TImage::Superclass _TImageBase;
+ typedef fpa::Image::Functors::SimpleNeighborhood< _TImageBase > _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 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 );
+ 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->SetCostConversionFunction( ifunc );
}
// -------------------------------------------------------------------------
-template< class _TDistanceMap, class _TCostMap >
-fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
+template< class _TImage >
+fpa::Image::SkeletonFilter< _TImage >::
~SkeletonFilter( )
{
}
// -------------------------------------------------------------------------
-template< class _TDistanceMap, class _TCostMap >
-void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
-GenerateData( )
+template< class _TImage >
+void fpa::Image::SkeletonFilter< _TImage >::
+_BeforeGenerateData( )
{
- // 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 );
+ this->Superclass::_BeforeGenerateData( );
+ this->m_SkeletonQueue.clear( );
}
// -------------------------------------------------------------------------
-template< class _TDistanceMap, class _TCostMap >
-void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
-_EndPoints(
- const TDistanceMap* dmap,
- const TCostMap* cmap,
- const TMST* mst,
- TIndices* end_points
- )
+template< class _TImage >
+void fpa::Image::SkeletonFilter< _TImage >::
+_UpdateResult( const _TQueueNode& n )
{
- 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 );
+ typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
- // 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 )
+ this->Superclass::_UpdateResult( n );
+ double d = double( this->GetInput( )->GetPixel( n.Vertex ) );
+ if( d >= double( 0 ) )
{
- double d = double( dIt.Get( ) );
- if( d > 0 )
- {
- double v = double( cIt.Get( ) ) / ( d * d );
- queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
+ // Update skeleton candidates
+ d += double( 1e-5 );
+ double v = double( n.Result ) / ( d * d );
+ this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) );
- } // fi
+ } // fi
+}
- } // rof
+// -------------------------------------------------------------------------
+template< class _TImage >
+void fpa::Image::SkeletonFilter< _TImage >::
+_AfterGenerateData( )
+{
+ this->Superclass::_AfterGenerateData( );
+
+ std::vector< TIndex > end_points;
+ this->_EndPoints( end_points );
+ this->_Skeleton( end_points );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void fpa::Image::SkeletonFilter< _TImage >::
+_EndPoints( std::vector< TIndex >& end_points )
+{
+ typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
+
+ static const double _eps = std::sqrt( double( Self::Dimension + 1 ) );
+ auto input = this->GetInput( );
+ auto mst = this->GetMinimumSpanningTree( );
+ auto marks = this->GetMarks( );
+
+ // Some values
+ marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
+ marks->SetRequestedRegion( input->GetRequestedRegion( ) );
+ marks->SetBufferedRegion( input->GetBufferedRegion( ) );
+ marks->SetSpacing( input->GetSpacing( ) );
+ marks->SetOrigin( input->GetOrigin( ) );
+ marks->SetDirection( input->GetDirection( ) );
+ marks->Allocate( );
+ marks->FillBuffer( 0 );
+ auto spac = marks->GetSpacing( );
// BFS from maximum queue
- while( queue.size( ) > 0 )
+ auto region = input->GetRequestedRegion( );
+ while( this->m_SkeletonQueue.size( ) > 0 )
{
// Get node
- auto nIt = queue.begin( );
+ auto nIt = this->m_SkeletonQueue.begin( );
auto n = *nIt;
- queue.erase( nIt );
+ this->m_SkeletonQueue.erase( nIt );
+ // Mark it and update end-points
unsigned char m = marks->GetPixel( n.second );
- if( ( m & 0x01 ) == 0x01 )
+ if( m != 0 )
continue;
+ marks->SetPixel( n.second, 1 );
+ end_points.push_back( n.second );
- // 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 )
+ // Mark path
+ TIndex it = n.second;
+ TIndex p = mst->GetParent( it );
+ while( it != p )
{
- 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 );
+ // TODO: tags->SetPixel( it, tags->GetPixel( it ) + 1 );
+ it = p;
+ p = mst->GetParent( it );
+
+ } // elihw
+ // TODO: tags->SetPixel( it, tags->GetPixel( it ) + 1 );
+
+ /* TODO: use mst directly rather than computing paths
+ typename TMST::TPath::Pointer path;
+ mst->GetPath( path, n.second );
+ for( unsigned long i = 0; i < path->GetSize( ); ++i )
+ {
+ TIndex idx = path->GetVertex( i );
+ typename _TImage::PointType cnt;
+ input->TransformIndexToPhysicalPoint( idx, cnt );
+ double r = double( input->GetPixel( idx ) ) * _eps;
- 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 );
+ 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;
- } // fi
+ if( i0[ d ] < region.GetIndex( )[ d ] )
+ i0[ d ] = region.GetIndex( )[ d ];
- } // rof
+ if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
+ i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
- } // rof
+ } // rof
- } // elihw
+ typename _TImage::SizeType size;
+ for( unsigned int d = 0; d < Self::Dimension; ++d )
+ size[ d ] = i1[ d ] - i0[ d ] + 1;
- } // rof
+ typename _TImage::RegionType neighRegion;
+ neighRegion.SetIndex( i0 );
+ neighRegion.SetSize( size );
+
+ _TMarksIt mIt( marks, neighRegion );
+ for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
+ {
+ TIndex w = mIt.GetIndex( );
+ typename _TImage::PointType p;
+ marks->TransformIndexToPhysicalPoint( w, p );
+ mIt.Set( 1 );
+
+ } // rof
+
+ } // 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
- )
+template< class _TImage >
+void fpa::Image::SkeletonFilter< _TImage >::
+_Skeleton( const std::vector< TIndex >& end_points )
{
+ 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 )
+ {
+ TIndex it = *eIt;
+ TIndex 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 )
+ {
+ 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 );
+ it = p;
+ p = mst->GetParent( it );
+
+ } // elihw
+
+ // Finally, add last branch
+ path->AddVertex( it );
+ skeleton->AddBranch( path );
+
+ } // rof
}
#endif // __fpa__Image__SkeletonFilter__hxx__
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / blobdiff