X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FImage%2FDijkstraWithEndPointDetection.hxx;h=dc51a6442c0a66f7ed77663c25e8bb3c161923a7;hb=f28d460c6d9ce93e55c073bd04dba8de16d55d3a;hp=8ddd9baa107d81d8cb7263efc6e656c88b6440d5;hpb=5e326afe442245572b6c3ec98ebeec8b45f9012f;p=FrontAlgorithms.git diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx index 8ddd9ba..dc51a64 100644 --- a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx +++ b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx @@ -4,6 +4,7 @@ #include #include #include +#include // ------------------------------------------------------------------------- template< class I, class O > @@ -44,17 +45,153 @@ GraftLabelImage( itk::DataObject* obj ) this->GraftNthOutput( this->m_LabelImageIndex, lbl ); } + + + +// ------------------------------------------------------------------------- +template< class I, class O > +typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetEndPoints( ) +{ + return( + dynamic_cast< TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +const typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetEndPoints( ) const +{ + return( + dynamic_cast< const TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GraftEndPoints( itk::DataObject* obj ) +{ + TUniqueVertices* lbl = + dynamic_cast< TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex ) + ); + if( lbl != NULL ) + this->GraftNthOutput( this->m_EndPointsIndex, lbl ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetBifurcations( ) +{ + return( + dynamic_cast< TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +const typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetBifurcations( ) const +{ + return( + dynamic_cast< const TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GraftBifurcations( itk::DataObject* obj ) +{ + TUniqueVertices* lbl = + dynamic_cast< TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex ) + ); + if( lbl != NULL ) + this->GraftNthOutput( this->m_BifurcationsIndex, lbl ); +} + + +// ------------------------------------------------------------------------- +template< class I, class O > +typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetBranches( ) +{ + return( + dynamic_cast< TBranches* >( + this->itk::ProcessObject::GetOutput( this->m_BranchesIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +const typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetBranches( ) const +{ + return( + dynamic_cast< const TBranches* >( + this->itk::ProcessObject::GetOutput( this->m_BranchesIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GraftBranches( itk::DataObject* obj ) +{ + TBranches* lbl = + dynamic_cast< TBranches* >( + this->itk::ProcessObject::GetOutput( this->m_BranchesIndex ) + ); + if( lbl != NULL ) + this->GraftNthOutput( this->m_BranchesIndex, lbl ); +} + // ------------------------------------------------------------------------- template< class I, class O > fpa::Image::DijkstraWithEndPointDetection< I, O >:: DijkstraWithEndPointDetection( ) - : Superclass( ) + : Superclass( ), + m_CorrectSeeds( true ), + m_CorrectEndPoints( true ), + m_SafetyNeighborhoodSize( 3 ) { this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( ); - this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 1 ); + this->m_EndPointsIndex = this->m_LabelImageIndex + 1; + this->m_BifurcationsIndex = this->m_LabelImageIndex + 2; + this->m_BranchesIndex = this->m_LabelImageIndex + 3; + this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 4 ); this->itk::ProcessObject::SetNthOutput( this->m_LabelImageIndex, TLabelImage::New( ) ); + this->itk::ProcessObject::SetNthOutput( + this->m_EndPointsIndex, TUniqueVertices::New( ) + ); + this->itk::ProcessObject::SetNthOutput( + this->m_BifurcationsIndex, TUniqueVertices::New( ) + ); + this->itk::ProcessObject::SetNthOutput( + this->m_BranchesIndex, TBranches::New( ) + ); } // ------------------------------------------------------------------------- @@ -69,31 +206,43 @@ template< class I, class O > void fpa::Image::DijkstraWithEndPointDetection< I, O >:: _BeforeGenerateData( ) { - // Correct seeds - /* TODO - for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s ) - { - _TNode seed = this->m_Seeds[ s ]; - _TRegion region = this->_Region( seed.Vertex, max_spac ); - itk::ImageRegionConstIteratorWithIndex< I > iIt( img, region ); + if( this->m_CorrectSeeds ) + { + const I* input = this->GetInput( ); + typename I::SpacingType spac = input->GetSpacing( ); + double max_spac = double( spac[ 0 ] ); + for( unsigned int d = 1; d < I::ImageDimension; ++d ) + max_spac = + ( double( spac[ d ] ) > max_spac )? double( spac[ d ] ): max_spac; + max_spac *= double( 3 ); + + // Correct seeds + for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s ) + { + TVertex seed = this->m_SeedVertices[ s ]; + _TNode n = this->m_Seeds[ seed ]; + _TRegion region = this->_Region( seed, max_spac ); + itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); + + iIt.GoToBegin( ); + _TPixel max_value = iIt.Get( ); + for( ++iIt; !iIt.IsAtEnd( ); ++iIt ) + { + if( iIt.Get( ) > max_value ) + { + this->m_SeedVertices[ s ] = iIt.GetIndex( ); + max_value = iIt.Get( ); - iIt.GoToBegin( ); - _TPixel max_value = iIt.Get( ); - for( ++iIt; !iIt.IsAtEnd( ); ++iIt ) - { - if( iIt.Get( ) > max_value ) - { - seed.Vertex = iIt.GetIndex( ); - seed.Parent = seed.Vertex; - max_value = iIt.Get( ); + } // fi - } // fi + } // rof + this->m_Seeds.erase( seed ); + n.Parent = this->m_SeedVertices[ s ]; + this->m_Seeds[ this->m_SeedVertices[ s ] ] = n; - } // rof - this->m_Seeds[ s ] = seed; + } // rof - } // rof - */ + } // fi // End initialization this->Superclass::_BeforeGenerateData( ); @@ -108,240 +257,246 @@ _AfterGenerateData( ) // Finish base algorithm this->Superclass::_AfterGenerateData( ); - // Prepare backtracking objects - this->m_EndPoints.clear( ); - this->m_BifurcationPoints.clear( ); + // Check if backtracking has any sense if( this->m_Candidates.size( ) == 0 ) return; this->InvokeEvent( TEndEvent( ) ); this->InvokeEvent( TStartBacktrackingEvent( ) ); - // Get some input values - const I* input = this->GetInput( ); - typename I::SpacingType spac = input->GetSpacing( ); - double ms = double( spac[ 0 ] ); - for( unsigned int d = 1; d < I::ImageDimension; ++d ) - ms =( ms < double( spac[ d ] ) )? double( spac[ d ] ): ms; + // Detect endpoints and bifurcations + this->_EndPointsAndBifurcations( ); - // Prepare labels - TLabelImage* label = this->GetLabelImage( ); - label->FillBuffer( 0 ); + // Find branches + this->_FindBranches( ); - // Object to hold marks - std::set< TVertex, TVertexCompare > tree_marks; + // Label pixels and branches + this->_LabelAll( ); +} - // Object to hold all branches - this->m_Branches.clear( ); +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_SetResult( const TVertex& v, const _TNode& n ) +{ + this->Superclass::_SetResult( v, n ); + this->m_Candidates.insert( _TCandidate( n.Result, v ) ); +} - // First label - this->m_NumberOfBranches = 1; +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_EndPointsAndBifurcations( ) +{ + // Prepare backtracking objects + TUniqueVertices* endpoints = this->GetEndPoints( ); + TUniqueVertices* bifurcations = this->GetBifurcations( ); + endpoints->Clear( ); + bifurcations->Clear( ); + + // Get some input values + TVertex seed = this->GetSeed( 0 ); + const I* input = this->GetInput( ); + const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( ); + + // Prepare a pixel marks structure + typedef std::map< TVertex, bool, TVertexCompare > _TPixelMarks; + _TPixelMarks pixel_marks; + + // Object to hold tree-nodes marks + typedef std::set< TVertex, TVertexCompare > _TTreeMarks; + _TTreeMarks tree_marks; // Iterate over the candidates, starting from the candidates that - // are near thin branches + // are near thin branches (high accumulated cost) typename _TCandidates::const_reverse_iterator cIt = this->m_Candidates.rbegin( ); for( ; cIt != this->m_Candidates.rend( ); ++cIt ) { - // If pixel has been already labelled, pass + // If pixel has already been visited, pass TVertex v = cIt->second; - if( label->GetPixel( v ) != 0 ) + if( pixel_marks[ v ] ) continue; - // Compute nearest start candidate - TVertex endpoint = - this->_MaxInRegion( - input, v, - double( std::sqrt( input->GetPixel( v ) ) ) * double( 1.5 ) - ); + // Correct it to nearest start candidate (high distance value) + double vr = std::sqrt( double( input->GetPixel( v ) ) ); + if( this->m_CorrectEndPoints ) + v = this->_MaxInRegion( input, v, vr * double( 1.5 ) ); - // Re-check labelling - if( this->_Node( endpoint ).Label == Self::FarLabel ) + // Now, check for real marking conditions + // 1. Has it been visited by dijkstra? + if( this->_Node( v ).Label == Self::FarLabel ) continue; - if( label->GetPixel( endpoint ) != 0 ) + // 2. Is it already marked? + if( pixel_marks[ v ] ) continue; - if( this->m_EndPoints.find( endpoint ) != this->m_EndPoints.end( ) ) + // 3. Is it already an endpoint? + if( endpoints->Find( v ) ) continue; - this->m_EndPoints.insert( endpoint ); - std::cout << "endpoint " << endpoint << " inserted " << this->m_EndPoints.size( ) << std::endl; - - // Get the path all the way to seed - std::vector< TVertex > path; - this->GetMinimumSpanningTree( )-> - GetPath( path, endpoint, this->GetSeed( 0 ) ); - - // Mark branches - bool start = true; - bool change = false; - TVertex last_start = endpoint; - typename std::vector< TVertex >::const_iterator pIt = path.begin( ); - for( ; pIt != path.end( ); ++pIt ) + // 4. Ok, it is completely new! + endpoints->Insert( v ); + + // Get the path all the way to global seed + TVertices path; + mst->GetPath( path, v, seed ); + + // Backtracking to find endpoints and bifurcations + bool adding_new_points = true; + typename TVertices::const_iterator pIt = path.begin( ); + for( ; pIt != path.end( ) && adding_new_points; ++pIt ) { - this->InvokeEvent( - TBacktrackingEvent( *pIt, this->m_NumberOfBranches ) - ); - if( start ) + this->InvokeEvent( TBacktrackingEvent( *pIt, 1 ) ); + if( tree_marks.find( *pIt ) == tree_marks.end( ) ) { - if( tree_marks.find( *pIt ) == tree_marks.end( ) ) + // Mark current point as a tree point + tree_marks.insert( *pIt ); + + // Mark a region around current point as visited + vr = std::sqrt( double( input->GetPixel( *pIt ) ) ); + _TRegion region = this->_Region( *pIt, vr ); + if( region.GetNumberOfPixels( ) > 0 ) { - // Mark a region around current point as visited - tree_marks.insert( *pIt ); - this->_Label( *pIt, TLabel( this->m_NumberOfBranches ) ); + itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); + for( iIt.GoToBegin( ); !iIt.IsAtEnd( ); ++iIt ) + pixel_marks[ iIt.GetIndex( ) ] = true; } else - { - // A bifurcation point has been reached! - if( *pIt != this->GetSeed( 0 ) ) - { - this->m_Branches[ *pIt ][ last_start ] = this->m_NumberOfBranches; - last_start = *pIt; - this->m_BifurcationPoints.insert( *pIt ); - std::cout << "bifurcation = " << *pIt << " " << this->GetSeed( 0 ) << std::endl; - this->m_NumberOfBranches++; - start = false; - - } // fi - - } // fi + pixel_marks[ *pIt ] = true; } else { - if( - std::find( - this->m_BifurcationPoints.begin( ), - this->m_BifurcationPoints.end( ), - *pIt - ) == this->m_BifurcationPoints.end( ) - ) - { - // Mark a sphere around current point as visited - this->_Label( *pIt, TLabel( this->m_NumberOfBranches ) ); - change = true; - } - else + // A bifurcation point has been reached! + if( *pIt != seed ) { - // A bifurcation point has been reached! - if( *pIt != this->GetSeed( 0 ) ) - { - this->m_Branches[ *pIt ][ last_start ] = this->m_NumberOfBranches; - last_start = *pIt; - this->m_NumberOfBranches++; - - } // fi + bifurcations->Insert( *pIt ); + adding_new_points = false; } // fi } // fi } // rof - if( start || change ) - this->m_NumberOfBranches++; this->InvokeEvent( TEndBacktrackingEvent( ) ); } // rof +} - typename TBranches::const_iterator bIt = this->m_Branches.begin( ); - unsigned int leo = 0; - for( ; bIt != this->m_Branches.end( ); ++bIt ) +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_FindBranches( ) +{ + // Configure and obtain inputs + TVertex seed = this->GetSeed( 0 ); + const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( ); + const TUniqueVertices* endpoints = this->GetEndPoints( ); + const TUniqueVertices* bifurcations = this->GetBifurcations( ); + TBranches* branches = this->GetBranches( ); + branches->Clear( ); + + // Reconstruct pixels + typename TUniqueVertices::ConstIterator eIt = endpoints->Begin( ); + for( ; eIt != endpoints->End( ); ++eIt ) { - typename TBranch::const_iterator brIt = bIt->second.begin( ); - for( ; brIt != bIt->second.end( ); ++brIt ) + // Get the path all the way to global seed + TVertices path; + mst->GetPath( path, *eIt, seed ); + + TVertex start_vertex = *eIt; + typename TVertices::const_iterator pIt = path.begin( ); + for( ; pIt != path.end( ); ++pIt ) { - std::cout << bIt->first << " " << brIt->first << std::endl; - leo++; - } - } - - // Re-enumerate labels - std::map< TLabel, unsigned long > histo; - for( - typename std::set< TVertex, TVertexCompare >::iterator treeIt = tree_marks.begin( ); - treeIt != tree_marks.end( ); - ++treeIt - ) - histo[ label->GetPixel( *treeIt ) ]++; - - std::map< TLabel, TLabel > changes; - TLabel last_change = 1; - for( TLabel i = 1; i <= this->m_NumberOfBranches; ++i ) - { - if( histo[ i ] != 0 ) - changes[ i ] = last_change++; + if( bifurcations->Find( *pIt ) ) + { + branches->SetValue( start_vertex, *pIt, 1 ); + start_vertex = *pIt; - } // rof - this->m_NumberOfBranches = changes.size( ); + } // fi - std::cout << leo << " - " << this->m_EndPoints.size( ) << " - " << this->m_BifurcationPoints.size( ) << " - " << this->m_NumberOfBranches << std::endl; + } // rof - /* - for( - typename TTree::iterator treeIt = this->m_FullTree.begin( ); - treeIt != this->m_FullTree.end( ); - ++treeIt - ) - { - TMark old = treeIt->second.second; - treeIt->second.second = changes[ old ]; + // Finish with branches to global seed + branches->SetValue( start_vertex, seed, 1 ); - } // fi + } // rof +} - // Construct reduced tree - for( - typename TVertices::const_iterator eIt = this->m_EndPoints.begin( ); - eIt != this->m_EndPoints.end( ); - ++eIt - ) - { - typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt ); - if( tIt != this->m_FullTree.end( ) ) - { - TMark id = tIt->second.second; - do +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_LabelAll( ) +{ + // Configure and obtain inputs + const I* input = this->GetInput( ); + const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( ); + TBranches* branches = this->GetBranches( ); + TLabelImage* label = this->GetLabelImage( ); + label->FillBuffer( 0 ); + + // Label branches + typename TBranches::Iterator bIt = branches->Begin( ); + TLabel actual_label = 0; + for( ; bIt != branches->End( ); ++bIt ) + { + typename TBranches::RowIterator brIt = branches->Begin( bIt ); + for( ; brIt != branches->End( bIt ); ++brIt ) { - tIt = this->m_FullTree.find( tIt->second.first ); - if( tIt == this->m_FullTree.end( ) ) - break; + actual_label++; + brIt->second = actual_label; - } while( tIt->second.second == id ); - this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id ); + TVertices path; + mst->GetPath( path, bIt->first, brIt->first ); + typename TVertices::const_iterator pIt = path.begin( ); + for( ; pIt != path.end( ); ++pIt ) + { + double d = std::sqrt( double( input->GetPixel( *pIt ) ) ); + _TRegion region = this->_Region( *pIt, d ); + itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); + itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region ); + iIt.GoToBegin( ); + lIt.GoToBegin( ); + for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt ) + { + // Mask real input image + if( + !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) && + lIt.Get( ) == TLabel( 0 ) + ) + lIt.Set( actual_label ); + + } // rof - } // fi + } // rof } // rof - for( - typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( ); - bIt != this->m_BifurcationPoints.end( ); - ++bIt - ) - { - typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt ); - if( tIt != this->m_FullTree.end( ) ) - { - TMark id = tIt->second.second; - do - { - tIt = this->m_FullTree.find( tIt->second.first ); - if( tIt == this->m_FullTree.end( ) ) - break; + } // rof + this->m_NumberOfBranches = actual_label; - } while( tIt->second.second == id ); - this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id ); + // Label bifurcations + /* TODO: Is this necessary? + typename TUniqueVertices::const_iterator bifIt = + this->m_Bifurcations.begin( ); + for( ; bifIt != this->m_Bifurcations.end( ); ++bifIt ) + { + actual_label++; + double d = std::sqrt( double( input->GetPixel( *bifIt ) ) ); + _TRegion region = this->_Region( *bifIt, d ); + itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); + itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region ); + iIt.GoToBegin( ); + lIt.GoToBegin( ); + for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt ) + { + // Mask real input image + if( !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) ) + lIt.Set( actual_label ); - } // fi + } // rof - } // rof + } // rof */ } -// ------------------------------------------------------------------------- -template< class I, class O > -void fpa::Image::DijkstraWithEndPointDetection< I, O >:: -_SetResult( const TVertex& v, const _TNode& n ) -{ - this->Superclass::_SetResult( v, n ); - this->m_Candidates.insert( _TCandidate( n.Result, v ) ); -} - // ------------------------------------------------------------------------- template< class I, class O > typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: @@ -359,7 +514,10 @@ _Region( const TVertex& c, const double& r ) _TSize size; for( unsigned int d = 0; d < I::ImageDimension; ++d ) { - long s = long( std::ceil( r / double( spac[ d ] ) ) ) + 3; + // NOTE: 3 is a minimum neighborhood size + long s = + long( std::ceil( r / double( spac[ d ] ) ) ) + + long( this->m_SafetyNeighborhoodSize ); i0[ d ] = c[ d ] - s; i1[ d ] = c[ d ] + s; @@ -405,25 +563,6 @@ _MaxInRegion( const _T* image, const TVertex& v, const double& r ) return( max_vertex ); } -// ------------------------------------------------------------------------- -template< class I, class O > -void fpa::Image::DijkstraWithEndPointDetection< I, O >:: -_Label( const TVertex& v, const TLabel& l ) -{ - typedef itk::ImageRegionIteratorWithIndex< TLabelImage > _TIt; - - double d = std::sqrt( double( this->GetInput( )->GetPixel( v ) ) ); - _TRegion region = this->_Region( v, d ); - if( region.GetNumberOfPixels( ) > 0 ) - { - _TIt lIt( this->GetLabelImage( ), region ); - for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt ) - lIt.Set( l ); - } - else - this->GetLabelImage( )->SetPixel( v, l ); -} - #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__ // eof - $RCSfile$