X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FImage%2FSkeletonFilter.hxx;h=aa02e5976163499fae67d1f13f05aab9f48d40f9;hb=6c0b77c2a8e3b821ccbe9c72c705fcd561bb90c2;hp=65a3fe2431187d8a8ed4a550fe9c6d87c1b3c53d;hpb=b9128432da0c510bb4ac9a165cc84c065b5e34e8;p=FrontAlgorithms.git diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx index 65a3fe2..aa02e59 100644 --- a/lib/fpa/Image/SkeletonFilter.hxx +++ b/lib/fpa/Image/SkeletonFilter.hxx @@ -1,244 +1,352 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + #ifndef __fpa__Image__SkeletonFilter__hxx__ #define __fpa__Image__SkeletonFilter__hxx__ +#include +#include +#include +#include + // ------------------------------------------------------------------------- -#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 ); +template< class _TInputImage, class _TDistanceMap > +fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra:: +_TDijkstra( ) + : Superclass( ) +{ + // 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 ); +} // ------------------------------------------------------------------------- -#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 \ - { \ - return( \ - dynamic_cast< const output_type* >( \ - this->Superclass::GetOutput( output_id ) \ - ) \ - ); \ - } - -__fpa__Image__SkeletonFilter__OutputMacro( Skeleton, TSkeleton, 0 ); -__fpa__Image__SkeletonFilter__OutputMacro( EndPoints, TIndices, 1 ); -__fpa__Image__SkeletonFilter__OutputMacro( Bifurcations, TIndices, 2 ); +template< class _TInputImage, class _TDistanceMap > +fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra:: +~_TDijkstra( ) +{ +} // ------------------------------------------------------------------------- -template< class _TDistanceMap, class _TCostMap > -fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: -SkeletonFilter( ) - : Superclass( ) +template< class _TInputImage, class _TDistanceMap > +void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra:: +_BeforeGenerateData( ) { - this->SetNumbeOfRequiredInputs( 3 ); - this->SetNumbeOfRequiredOutputs( 3 ); - - 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 ); + this->Superclass::_BeforeGenerateData( ); + this->m_SkeletonQueue.clear( ); } // ------------------------------------------------------------------------- -template< class _TDistanceMap, class _TCostMap > -fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: -~SkeletonFilter( ) +template< class _TInputImage, class _TDistanceMap > +void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra:: +_UpdateOutputValue( TNode& n ) { + typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue; + + 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.Value ) / ( d * d ); + this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) ); + + } // fi } // ------------------------------------------------------------------------- -template< class _TDistanceMap, class _TCostMap > -void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: -GenerateData( ) +template< class _TInputImage, class _TDistanceMap > +itk::ModifiedTimeType +fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: +GetMTime( ) const { - // 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 ); + itk::ModifiedTimeType t = this->Superclass::GetMTime( ); + itk::ModifiedTimeType q = this->m_DistanceMap->GetMTime( ); + return( ( q < t )? q: t ); } // ------------------------------------------------------------------------- -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 _TInputImage, class _TDistanceMap > +void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: +SetInput( TInputImage* input ) { - 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; + this->Superclass::SetNthInput( 0, input ); +} - // 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 ); - - // 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 ) +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TDistanceMap > +typename fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: +TInputImage* fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: +GetInput( ) +{ + return( dynamic_cast< TInputImage* >( this->Superclass::GetInput( 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 ) ) + ); +} + +// ------------------------------------------------------------------------- +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 ) ) ); +} + +// ------------------------------------------------------------------------- +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 ) ) + ); +} + +// ------------------------------------------------------------------------- +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( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TDistanceMap > +fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: +~SkeletonFilter( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TDistanceMap > +void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: +GenerateData( ) +{ + // Update distance map + this->m_DistanceMap->SetInput( this->GetInput( ) ); + this->m_DistanceMap->Update( ); + + // Correct seed + if( this->m_SeedFromMaximumDistance ) { - double d = double( dIt.Get( ) ); - if( d > 0 ) - { - double v = double( cIt.Get( ) ) / ( d * d ); - queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) ); + typedef itk::MinimumMaximumImageCalculator< TOutputImage > _TMinMax; + typename _TMinMax::Pointer minmax = _TMinMax::New( ); + minmax->SetImage( this->m_DistanceMap->GetOutput( ) ); + minmax->Compute( ); + this->m_Seed = minmax->GetIndexOfMaximum( ); - } // fi + } // fi - } // rof + // Compute MST + typename _TDijkstra::Pointer dijkstra = _TDijkstra::New( ); + dijkstra->SetInput( this->m_DistanceMap->GetOutput( ) ); + dijkstra->AddSeed( this->m_Seed ); + dijkstra->Update( ); - // BFS from maximum queue - while( queue.size( ) > 0 ) + // Compute end-points + const _TMST* mst = dijkstra->GetMinimumSpanningTree( ); + std::vector< TIndex > end_points; + this->_EndPoints( + end_points, this->m_DistanceMap->GetOutput( ), + mst, dijkstra->GetSkeletonQueue( ) + ); + + // Compute symbolic branches + typedef std::map< TIndex, TIndex, typename TIndex::LexicographicCompare > _TTags; + _TTags tags, branches; + typename std::vector< TIndex >::const_iterator eIt = end_points.begin( ); + for( ; eIt != end_points.end( ); ++eIt ) { - // Get node - auto nIt = queue.begin( ); - auto n = *nIt; - queue.erase( nIt ); + // Tag path + TIndex it = *eIt; + TIndex p = mst->GetParent( it ); + typename _TTags::iterator bIt = tags.end( ); + while( it != p && bIt == tags.end( ) ) + { + typename _TTags::iterator tIt = tags.find( it ); + if( tIt != tags.end( ) ) + { + // Ok, a bifurcation point has been found + // branch1: tIt->second <-> it (ok) + branches[ tIt->second ] = it; - unsigned char m = marks->GetPixel( n.second ); - if( ( m & 0x01 ) == 0x01 ) - continue; + // branch2: *eit <-> it (ok) + branches[ *eIt ] = it; - // 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 ) + // branch3: it <-> until next bifurcation + bIt = tIt; + } + else + tags[ it ] = *eIt; + it = p; + p = mst->GetParent( it ); + + } // elihw + if( bIt != tags.end( ) ) { - 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 pTag = bIt->second; + TIndex nTag = bIt->first; + it = bIt->first; + p = it; + while( tags[ it ] == pTag ) { - TIndex v = q.front( ); - q.pop( ); - unsigned char m = marks->GetPixel( v ); - if( ( m & 0x02 ) == 0x02 ) - continue; - m |= 0x03; - marks->SetPixel( v, m ); - - 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 ); - - } // fi - - } // rof - - } // rof + tags[ it ] = nTag; + p = it; + it = mst->GetParent( it ); } // elihw + tags[ it ] = nTag; + branches[ bIt->first ] = p; + } + else + { + tags[ it ] = *eIt; + branches[ *eIt ] = it; + + } // fi + + } // rof + + // Fill full branches + typedef typename _TMST::TVertices _TVertices; + typedef typename TSkeleton::TPath _TPath; + + TSkeleton* sk = this->GetOutput( ); + typename _TTags::const_iterator bIt = branches.begin( ); + for( ; bIt != branches.end( ); ++bIt ) + { + _TVertices v = mst->GetPath( bIt->first, bIt->second ); + typename _TPath::Pointer path = _TPath::New( ); + path->SetReferenceImage( this->GetInput( ) ); + typename _TVertices::const_reverse_iterator vIt = v.rbegin( ); + for( ; vIt != v.rend( ); ++vIt ) + path->AddVertex( *vIt ); + sk->AddBranch( path ); + + } // rof +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TDistanceMap > +template< class _TMarksPointer > +void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: +_MarkSphere( + _TMarksPointer& marks, const TOutputImage* dmap, const TIndex& center + ) +{ + typedef typename _TMarksPointer::ObjectType _TMarks; + typedef itk::ImageRegionIteratorWithIndex< _TMarks > _TMarksIt; + + static const double _eps = std::sqrt( double( Self::Dimension + 1 ) ); + typename _TMarks::SpacingType spac = dmap->GetSpacing( ); + typename _TMarks::RegionType region = dmap->GetRequestedRegion( ); - } // rof + typename _TMarks::PointType cnt; + dmap->TransformIndexToPhysicalPoint( center, cnt ); + double r = double( dmap->GetPixel( center ) ) * _eps; - } // elihw + 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 ] = center[ d ] - off; + i1[ d ] = center[ 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 < Self::Dimension; ++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 ) + mIt.Set( true ); } // ------------------------------------------------------------------------- -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 _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 ) + { + // 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 } #endif // __fpa__Image__SkeletonFilter__hxx__