X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FImage%2FSkeletonFilter.hxx;h=aa02e5976163499fae67d1f13f05aab9f48d40f9;hb=6c0b77c2a8e3b821ccbe9c72c705fcd561bb90c2;hp=f1e315fe6ca3506f83956194614b9bf75b7f2e27;hpb=d5fe8fd4bac61fafea412c358323ad90fe2b034b;p=FrontAlgorithms.git diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx index f1e315f..aa02e59 100644 --- a/lib/fpa/Image/SkeletonFilter.hxx +++ b/lib/fpa/Image/SkeletonFilter.hxx @@ -6,314 +6,276 @@ #ifndef __fpa__Image__SkeletonFilter__hxx__ #define __fpa__Image__SkeletonFilter__hxx__ +#include #include #include +#include -#include -#include +// ------------------------------------------------------------------------- +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 ); +} // ------------------------------------------------------------------------- -template< class _TImage, class _TScalar > -void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -SetInput( const TScalarImage* image ) +template< class _TInputImage, class _TDistanceMap > +fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra:: +~_TDijkstra( ) { - // Do nothing } // ------------------------------------------------------------------------- -template< class _TImage, class _TScalar > -void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -SetInput( const TImage* image ) +template< class _TInputImage, class _TDistanceMap > +void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra:: +_BeforeGenerateData( ) { - this->m_DistanceMap->SetInput( image ); - this->Superclass::SetInput( this->m_DistanceMap->GetOutput( ) ); + this->Superclass::_BeforeGenerateData( ); + this->m_SkeletonQueue.clear( ); } // ------------------------------------------------------------------------- -template< class _TImage, class _TScalar > -_TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -GetInput( ) +template< class _TInputImage, class _TDistanceMap > +void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra:: +_UpdateOutputValue( TNode& n ) { - return( this->m_DistanceMap->GetInput( ) ); + 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 _TImage, class _TScalar > -const _TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -GetInput( ) const +template< class _TInputImage, class _TDistanceMap > +itk::ModifiedTimeType +fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: +GetMTime( ) const { - return( this->m_DistanceMap->GetInput( ) ); + itk::ModifiedTimeType t = this->Superclass::GetMTime( ); + itk::ModifiedTimeType q = this->m_DistanceMap->GetMTime( ); + return( ( q < t )? q: t ); } // ------------------------------------------------------------------------- -template< class _TImage, class _TScalar > -typename fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -TSkeleton* fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -GetSkeleton( ) +template< class _TInputImage, class _TDistanceMap > +void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: +SetInput( TInputImage* input ) +{ + this->Superclass::SetNthInput( 0, input ); +} + +// ------------------------------------------------------------------------- +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< TSkeleton* >( - this->itk::ProcessObject::GetOutput( this->m_SkeletonIdx ) - ) + dynamic_cast< const TInputImage* >( this->Superclass::GetInput( 0 ) ) ); } // ------------------------------------------------------------------------- -template< class _TImage, class _TScalar > -typename fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -TMarks* fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -GetMarks( ) +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< TMarks* >( - this->itk::ProcessObject::GetOutput( this->m_MarksIdx ) - ) + dynamic_cast< const TSkeleton* >( this->Superclass::GetOutput( 0 ) ) ); } // ------------------------------------------------------------------------- -template< class _TImage, class _TScalar > -fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +template< class _TInputImage, class _TDistanceMap > +fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: SkeletonFilter( ) : Superclass( ), m_SeedFromMaximumDistance( false ) { - typedef fpa::Image::Functors::VertexIdentity< TScalarImage, _TScalar > _TVertexFunc; - typedef fpa::Base::Functors::InvertValue< _TScalar, _TScalar > _TValueFunc; - - typename _TVertexFunc::Pointer vertex_func = _TVertexFunc::New( ); - typename _TValueFunc::Pointer value_func = _TValueFunc::New( ); - value_func->SetAlpha( 1 ); - value_func->SetBeta( 1 ); - - this->SetFunctor( vertex_func ); - this->SetFunctor( value_func ); + this->SetNumberOfRequiredInputs( 1 ); + this->SetNumberOfRequiredOutputs( 1 ); + this->SetNthOutput( 0, TSkeleton::New( ) ); 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 >:: +template< class _TInputImage, class _TDistanceMap > +fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: ~SkeletonFilter( ) { } // ------------------------------------------------------------------------- -template< class _TImage, class _TScalar > -void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +template< class _TInputImage, class _TDistanceMap > +void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: GenerateData( ) { - // Prepare distance map + // Update distance map + this->m_DistanceMap->SetInput( this->GetInput( ) ); this->m_DistanceMap->Update( ); - // Update start seed - TVertex seed; + // Correct seed if( this->m_SeedFromMaximumDistance ) { - typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax; + typedef itk::MinimumMaximumImageCalculator< TOutputImage > _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 ) -{ - typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue; - - this->Superclass::_SetOutputValue( vertex, value ); - double d = double( this->m_DistanceMap->GetOutput( )->GetPixel( vertex ) ); - if( d >= double( 0 ) ) - { - // Update skeleton candidates - d += double( 1e-5 ); - double v = double( value ) / ( d * d ); - this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, vertex ) ); + this->m_Seed = minmax->GetIndexOfMaximum( ); } // fi -} - -// ------------------------------------------------------------------------- -template< class _TImage, class _TScalar > -void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -_EndPoints( std::vector< TVertex >& end_points, _TAdjacencies& A ) -{ - typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue; - - TMarks* marks = this->GetMarks( ); - TMST* mst = this->GetMinimumSpanningTree( ); - typename TImage::SpacingType 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 ); + // Compute MST + typename _TDijkstra::Pointer dijkstra = _TDijkstra::New( ); + dijkstra->SetInput( this->m_DistanceMap->GetOutput( ) ); + dijkstra->AddSeed( this->m_Seed ); + dijkstra->Update( ); + + // Compute end-points + const _TMST* mst = dijkstra->GetMinimumSpanningTree( ); + std::vector< TIndex > end_points; + this->_EndPoints( + end_points, this->m_DistanceMap->GetOutput( ), + mst, dijkstra->GetSkeletonQueue( ) + ); - // BFS from maximum queue - while( this->m_SkeletonQueue.size( ) > 0 ) + // 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 - typename _TSkeletonQueue::iterator nIt = this->m_SkeletonQueue.begin( ); - _TSkeletonQueueValue 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 ) + // Tag path + TIndex it = *eIt; + TIndex p = mst->GetParent( it ); + typename _TTags::iterator bIt = tags.end( ); + while( it != p && bIt == tags.end( ) ) { - this->_MarkSphere( it ); - it = p; - p = mst->GetParent( it ); - - } // elihw - this->_MarkSphere( it ); - A[ n.second ] = it; + 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; - } // elihw -} + // branch2: *eit <-> it (ok) + branches[ *eIt ] = it; -// ------------------------------------------------------------------------- -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 ); + // branch3: it <-> until next bifurcation + bIt = tIt; + } + else + tags[ it ] = *eIt; it = p; p = mst->GetParent( it ); } // elihw - tags->SetPixel( it, tags->GetPixel( it ) + 1 ); + if( bIt != tags.end( ) ) + { + TIndex pTag = bIt->second; + TIndex nTag = bIt->first; + it = bIt->first; + p = it; + while( tags[ it ] == pTag ) + { + 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 - // Build paths (branches) - eIt = end_points.begin( ); - for( ; eIt != end_points.end( ); ++eIt ) + // 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 ) { - TVertex it = *eIt; - TVertex p = mst->GetParent( it ); - TVertex sIdx = it; + _TVertices v = mst->GetPath( bIt->first, bIt->second ); 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 ); + 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 _TImage, class _TScalar > -void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: -_MarkSphere( const TVertex& 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( TImage::ImageDimension + 1 ) ); - const TScalarImage* dmap = this->m_DistanceMap->GetOutput( ); - TMarks* marks = this->GetMarks( ); - typename TImage::SpacingType spac = dmap->GetSpacing( ); - typename TImage::RegionType region = dmap->GetRequestedRegion( ); + static const double _eps = std::sqrt( double( Self::Dimension + 1 ) ); + typename _TMarks::SpacingType spac = dmap->GetSpacing( ); + typename _TMarks::RegionType region = dmap->GetRequestedRegion( ); - typename TImage::PointType cnt; - dmap->TransformIndexToPhysicalPoint( idx, cnt ); - double r = double( dmap->GetPixel( idx ) ) * _eps; + typename _TMarks::PointType cnt; + dmap->TransformIndexToPhysicalPoint( center, cnt ); + double r = double( dmap->GetPixel( center ) ) * _eps; - TVertex i0, i1; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + 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 ]; @@ -323,21 +285,66 @@ _MarkSphere( const TVertex& idx ) } // rof - typename TImage::SizeType size; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + 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 ) { - TVertex 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 }