X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FImage%2FSkeletonFilter.hxx;h=069eea071af885b3ae432a8edebfa66d327e9ec1;hb=7955ba9683b8032e8cd3ca7aa361568fdbb218d2;hp=010c6945391c518eedd08a3fade1dca64723cd14;hpb=026b2fe203089e1917ab78ebafb3131f147223f5;p=FrontAlgorithms.git diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx index 010c694..069eea0 100644 --- a/lib/fpa/Image/SkeletonFilter.hxx +++ b/lib/fpa/Image/SkeletonFilter.hxx @@ -1,257 +1,268 @@ #ifndef __fpa__Image__SkeletonFilter__hxx__ #define __fpa__Image__SkeletonFilter__hxx__ -#include -#include +#include +#include #include // ------------------------------------------------------------------------- -#define fpa_Image_SkeletonFilter_InputMacro( i_n, i_t, i_i ) \ - template< class _TDistanceMap, class _TCostMap > \ - typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - i_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - Get##i_n( ) \ - { \ - return( \ - dynamic_cast< i_t* >( this->Superclass::GetInput( i_i ) ) \ - ); \ - } \ - template< class _TDistanceMap, class _TCostMap > \ - const \ - typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - i_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - Get##i_n( ) const \ - { \ - return( \ - dynamic_cast< const i_t* >( this->Superclass::GetInput( i_i ) ) \ - ); \ - } \ - template< class _TDistanceMap, class _TCostMap > \ - void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - Set##i_n( i_t* input ) \ - { \ - this->Superclass::SetNthInput( i_i, 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( o_n, o_t, o_i ) \ - template< class _TDistanceMap, class _TCostMap > \ - typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - o_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ +#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->Superclass::GetOutput( o_i ) ) \ - ); \ - } \ - template< class _TDistanceMap, class _TCostMap > \ - const typename \ - fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - o_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - Get##o_n( ) const \ - { \ - return( \ - dynamic_cast< const o_t* >( this->Superclass::GetOutput( o_i ) ) \ + 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( Marks, TMarks, 3 ); +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->SetNumberOfRequiredInputs( 3 ); - this->SetNumberOfRequiredOutputs( 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( ); typename TMarks::Pointer marks = TMarks::New( ); - this->SetNthOutput( 0, skeleton.GetPointer( ) ); - this->SetNthOutput( 1, end_points.GetPointer( ) ); - this->SetNthOutput( 2, bifurcations.GetPointer( ) ); - this->SetNthOutput( 3, marks.GetPointer( ) ); + 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( ) +{ + this->Superclass::_BeforeGenerateData( ); + this->m_SkeletonQueue.clear( ); +} + +// ------------------------------------------------------------------------- +template< class _TImage > +void fpa::Image::SkeletonFilter< _TImage >:: +_UpdateResult( const _TQueueNode& n ) { - // 0. I/O objects - const TDistanceMap* dmap = this->GetDistanceMap( ); - const TCostMap* cmap = this->GetCostMap( ); - const TMST* mst = this->GetMinimumSpanningTree( ); - TIndices* ep = this->GetEndPoints( ); - TIndices* bi = this->GetBifurcations( ); - TSkeleton* sk = this->GetSkeleton( ); - - // 1. Check input correspondance - // TODO - - // 2. Detect end-points - this->_EndPoints( dmap, cmap, mst, ep->Get( ) ); - std::cout << "endpoints" << std::endl; - - // 3. Build skeleton and keep track of bifurcations - this->_Skeleton( dmap, cmap, mst, ep->Get( ), bi->Get( ), sk ); + 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 _TDistanceMap, class _TCostMap > -void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: -_EndPoints( - const TDistanceMap* dmap, - const TCostMap* cmap, - const TMST* mst, - TIndicesData& end_points - ) +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::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::ImageRegionIteratorWithIndex< TMarks > _TMarksIt; - static const double _eps = std::sqrt( double( TMarks::ImageDimension + 1 ) ); + static const double _eps = std::sqrt( double( Self::Dimension + 1 ) ); + auto input = this->GetInput( ); + auto mst = this->GetMinimumSpanningTree( ); + auto marks = this->GetMarks( ); // Some values - auto marks = this->GetMarks( ); - 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->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 ); - - // 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 > double( 0 ) ) - { - double v = double( cIt.Get( ) ) / ( d * d ); - queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) ); - - } // fi - - } // rof + auto spac = marks->GetSpacing( ); // BFS from maximum queue - auto region = dmap->GetRequestedRegion( ); - double init_v = queue.begin( )->first; - 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 != 0 ) continue; - - std::cout << queue.size( ) << std::endl; - - // Mark it and update end-points - m |= 0x01; - marks->SetPixel( n.second, m ); - end_points.insert( n.second ); + marks->SetPixel( n.second, 1 ); + end_points.push_back( n.second ); // Mark path - auto spac = marks->GetSpacing( ); - typename TMST::TPath::Pointer path; - mst->GetPath( path, n.second ); - for( unsigned long i = 0; i < path->GetSize( ); ++i ) + TIndex it = n.second; + TIndex p = mst->GetParent( it ); + while( it != p ) { - TIndex idx = path->GetVertex( i ); - typename _TCostMap::PointType cnt; - cmap->TransformIndexToPhysicalPoint( idx, cnt ); - double r = double( dmap->GetPixel( idx ) ) * _eps; - - TIndex i0, i1; - for( unsigned int d = 0; d < TMarks::ImageDimension; ++d ) - { - long off = long( std::ceil( r / double( spac[ d ] ) ) ); - if( off == 0 ) - off = 1; - 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 TMarks::SizeType size; - for( unsigned int d = 0; d < TMarks::ImageDimension; ++d ) - size[ d ] = i1[ d ] - i0[ d ] + 1; - - typename TMarks::RegionType neighRegion; - neighRegion.SetIndex( i0 ); - neighRegion.SetSize( size ); - - _TMarksIt mIt( marks, neighRegion ); - for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) - { - TIndex w = mIt.GetIndex( ); - typename _TCostMap::PointType p; - marks->TransformIndexToPhysicalPoint( w, p ); - if( cnt.EuclideanDistanceTo( p ) <= r ) - mIt.Set( mIt.Get( ) | 0x02 ); - - } // rof - - } // rof + // 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; + + 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; + + 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 ) + { + 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, - const TIndicesData& end_points, - TIndicesData& 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 ) { - typename _TPath::Pointer path; - mst->GetPath( path, *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