X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FImage%2FSkeletonFilter.hxx;fp=lib%2Ffpa%2FImage%2FSkeletonFilter.hxx;h=0000000000000000000000000000000000000000;hb=3c639e5da479c7216a0a302ffa156ac6762caeed;hp=aa02e5976163499fae67d1f13f05aab9f48d40f9;hpb=5bf766068f54d061d3816f4950a076c3cf3a4d8b;p=FrontAlgorithms.git diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx deleted file mode 100644 index aa02e59..0000000 --- a/lib/fpa/Image/SkeletonFilter.hxx +++ /dev/null @@ -1,354 +0,0 @@ -// ========================================================================= -// @author Leonardo Florez Valencia -// @email florez-l@javeriana.edu.co -// ========================================================================= - -#ifndef __fpa__Image__SkeletonFilter__hxx__ -#define __fpa__Image__SkeletonFilter__hxx__ - -#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 _TInputImage, class _TDistanceMap > -fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra:: -~_TDijkstra( ) -{ -} - -// ------------------------------------------------------------------------- -template< class _TInputImage, class _TDistanceMap > -void fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >::_TDijkstra:: -_BeforeGenerateData( ) -{ - this->Superclass::_BeforeGenerateData( ); - this->m_SkeletonQueue.clear( ); -} - -// ------------------------------------------------------------------------- -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 _TInputImage, class _TDistanceMap > -itk::ModifiedTimeType -fpa::Image::SkeletonFilter< _TInputImage, _TDistanceMap >:: -GetMTime( ) const -{ - itk::ModifiedTimeType t = this->Superclass::GetMTime( ); - itk::ModifiedTimeType q = this->m_DistanceMap->GetMTime( ); - return( ( q < t )? q: t ); -} - -// ------------------------------------------------------------------------- -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< 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 ) - { - 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 - - // 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( ) - ); - - // 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 ) - { - // 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; - - // branch2: *eit <-> it (ok) - branches[ *eIt ] = it; - - // branch3: it <-> until next bifurcation - bIt = tIt; - } - else - tags[ it ] = *eIt; - it = p; - p = mst->GetParent( it ); - - } // elihw - 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 - - // 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( ); - - 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 ] = 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 _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__ - -// eof - $RCSfile$