X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FImage%2FDijkstraWithEndPointDetection.hxx;h=cb417f027092de357d1353ff18b4cd367f860729;hb=eb4acd3dde87a3e33593c3ce87d0d351dec23f69;hp=d8e344a26da44ca93c3c2f64d5bcb1516cdabde7;hpb=8480ad8f9b3e3b556aa8f3f25ff07daf1db6dc78;p=FrontAlgorithms.git

diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx
index d8e344a..cb417f0 100644
--- a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx
+++ b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx
@@ -4,6 +4,7 @@
 #include <vector>
 #include <set>
 #include <itkImageRegionConstIteratorWithIndex.h>
+#include <itkImageRegionIteratorWithIndex.h>
 
 // -------------------------------------------------------------------------
 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 );
-
-     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
-
-     } // rof
-     this->m_Seeds[ s ] = seed;
-
-     } // rof
-  */
+  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( );
+
+        } // fi
+
+      } // rof
+      this->m_Seeds.erase( seed );
+      n.Parent = this->m_SeedVertices[ s ];
+      this->m_Seeds[ this->m_SeedVertices[ s ] ] = n;
+
+    } // rof
+
+  } // fi
 
   // End initialization
   this->Superclass::_BeforeGenerateData( );
@@ -105,276 +254,246 @@ template< class I, class O >
 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 _AfterGenerateData( )
 {
+  // Finish base algorithm
   this->Superclass::_AfterGenerateData( );
 
-  // Finish base algorithm
-  /* TODO
-     this->m_FullTree.clear( );
-     this->m_ReducedTree.clear( );
-  */
-  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( ) );
 
+  // Detect endpoints and bifurcations
+  this->_EndPointsAndBifurcations( );
+
+  // Find branches
+  this->_FindBranches( );
+
+  // Label pixels and branches
+  this->_LabelAll( );
+}
+
+// -------------------------------------------------------------------------
+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 >
+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( );
-  typename I::SpacingType spac = input->GetSpacing( );
-  double max_spac = spac[ 0 ];
-  for( unsigned int d = 1; d < I::ImageDimension; ++d )
-    max_spac =
-      ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
-  max_spac *= double( 3 );
-  TLabelImage* label = this->GetLabelImage( );
-  label->FillBuffer( 0 );
+  const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
+
+  // Prepare a pixel marks structure
+  typedef std::map< TVertex, bool, TVertexCompare > _TPixelMarks;
+  _TPixelMarks pixel_marks;
 
-  // Prepare an object to hold marks
-  std::set< TVertex, TVertexCompare > tree_marks;
-  /* TODO
-     typename TMarkImage::Pointer marks = this->GetOutputMarkImage( );
-     marks->FillBuffer( 0 );
-  */
+  // 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( );
-  this->m_NumberOfBranches = 1;
-  std::map< TLabel, std::pair< TVertex, TVertex > > branches;
   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
-    _TRegion region = this->_Region( v, max_spac );
-    itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
-    iIt.GoToBegin( );
-    TVertex max_vertex = iIt.GetIndex( );
-    _TPixel max_value = iIt.Get( );
-    for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
-    {
-      _TPixel value = iIt.Get( );
-      if( value > max_value )
-      {
-        max_value = value;
-        max_vertex = iIt.GetIndex( );
+    // 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 ) );
 
-      } // fi
+    // Now, check for real marking conditions
+    //   1. Has it been visited by dijkstra?
+    if( this->_Node( v ).Label == Self::FarLabel )
+      continue;
+    //   2. Is it already marked?
+    if( pixel_marks[ v ] )
+      continue;
+    //   3. Is it already an endpoint?
+    if( endpoints->Find( v ) )
+      continue;
+    //   4. Ok, it is completely new!
+    endpoints->Insert( v );
 
-    } // rof
+    // Get the path all the way to global seed
+    TVertices path;
+    mst->GetPath( path, v, seed );
 
-    // Re-check labelling
-    if( label->GetPixel( max_vertex ) != 0 )
-      continue;
-    this->m_EndPoints.push_back( max_vertex );
-
-    // Get the path all the way to seed
-    std::vector< TVertex > path;
-    this->GetMinimumSpanningTree( )->
-      GetPath( path, max_vertex, this->GetSeed( 0 ) );
-
-    // Mark branches
-    bool start = true;
-    bool change = false;
-    TVertex last_start = max_vertex;
-    typename std::vector< TVertex >::const_iterator pIt = path.begin( );
-    for( ; pIt != path.end( ); ++pIt )
+    // 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 )
         {
-          tree_marks.insert( *pIt );
-
-          // Mark a sphere around current point as visited
-          double dist = std::sqrt( double( input->GetPixel( *pIt ) ) );
-          region = this->_Region( max_vertex, dist * double( 1.5 ) );
-          itk::ImageRegionIteratorWithIndex< TLabelImage >
-            lIt( label, region );
-          for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
-            lIt.Set( this->m_NumberOfBranches );
-
-          // Next vertex in current path
-          // TODO: this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
-          /*
-            this->m_FullTree[ max_vertex ] =
-            TTreeNode( this->_Parent( max_vertex ), 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!
-          branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt );
-          last_start = *pIt;
-          this->m_BifurcationPoints.push_back( *pIt );
-          this->m_NumberOfBranches++;
-          start = false;
-
-        } // 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
-          double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
-          region = this->_Region( max_vertex, dist * double( 1.5 ) );
-          itk::ImageRegionIteratorWithIndex< TLabelImage >
-            lIt( label, region );
-          for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
-            lIt.Set( this->m_NumberOfBranches );
-
-          // Next vertex in current path
-          /* TODO
-             this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
-             this->m_FullTree[ max_vertex ] =
-             TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
-          */
-          change = true;
-        }
-        else
+        // A bifurcation point has been reached!
+        if( *pIt != seed )
         {
-          // A bifurcation point has been reached!
-          // TODO: this->m_BifurcationPoints.push_back( max_vertex );
-          branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt );
-
-          last_start = *pIt;
-          this->m_NumberOfBranches++;
-          /* TODO
-             this->m_FullTree[ max_vertex ] =
-             TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
-          */
+          bifurcations->Insert( *pIt );
+          adding_new_points = false;
 
         } // fi
 
       } // fi
 
     } // rof
-    if( start || change )
-      this->m_NumberOfBranches++;
     this->InvokeEvent( TEndBacktrackingEvent( ) );
 
-
-
-    /*
-      this->InvokeEvent( TEndBacktrackingEvent( ) );
-    */
-
   } // rof
+}
 
-  // 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 ) ]++;
-
-  for(
-    typename std::map< TLabel, unsigned long >::iterator hIt = histo.begin( );
-    hIt != histo.end( );
-    ++hIt
-    )
-    std::cout << hIt->first << " " << hIt->second << std::endl;
-
-  /*
-  std::map< TMark, TMark > changes;
-  TMark last_change = 1;
-  for( TMark i = 1; i <= this->m_NumberOfBranches; ++i )
+// -------------------------------------------------------------------------
+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 )
   {
-    if( histo[ i ] != 0 )
-      changes[ i ] = last_change++;
+    // Get the path all the way to global seed
+    TVertices path;
+    mst->GetPath( path, *eIt, seed );
 
-  } // rof
-  this->m_NumberOfBranches = changes.size( );
-  */
-
-  /*
-    for(
-    typename TTree::iterator treeIt = this->m_FullTree.begin( );
-    treeIt != this->m_FullTree.end( );
-    ++treeIt
-    )
+    TVertex start_vertex = *eIt;
+    typename TVertices::const_iterator pIt = path.begin( );
+    for( ; pIt != path.end( ); ++pIt )
     {
-    TMark old = treeIt->second.second;
-    treeIt->second.second = changes[ old ];
+      if( bifurcations->Find( *pIt ) )
+      {
+        branches->SetValue( start_vertex, *pIt, 1 );
+        start_vertex = *pIt;
 
-    } // fi
+      } // fi
 
-    // 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
-    {
-    tIt = this->m_FullTree.find( tIt->second.first );
-    if( tIt == this->m_FullTree.end( ) )
-    break;
+    } // rof
 
-    } while( tIt->second.second == id );
-    this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id );
+    // Finish with branches to global seed
+    branches->SetValue( start_vertex, seed, 1 );
 
-    } // fi
+  } // rof
+}
 
-    } // rof
+// -------------------------------------------------------------------------
+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( );
+  const TUniqueVertices* bifurcations = this->GetBifurcations( );
+  TBranches* branches = this->GetBranches( );
+  TLabelImage* label = this->GetLabelImage( );
+  label->FillBuffer( 0 );
 
-    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
+  // 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[ *bIt ] = 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
-  */
-}
 
-// -------------------------------------------------------------------------
-template< class I, class O >
-void fpa::Image::DijkstraWithEndPointDetection< I, O >::
-_SetResult( const TVertex& v, const _TNode& n )
-{
-  this->Superclass::_SetResult( v, n );
+  } // rof
+  this->m_NumberOfBranches = actual_label;
 
-  TResult vv = TResult( this->_VertexValue( v ) );
-  if( TResult( 0 ) < vv )
-    this->m_Candidates.insert( _TCandidate( n.Result / vv, v ) );
+  // Label bifurcations
+  typename TUniqueVertices::ConstIterator bifIt =
+    bifurcations->Begin( );
+  for( ; bifIt != bifurcations->End( ); ++bifIt )
+  {
+    actual_label++;
+    double d = std::sqrt( double( input->GetPixel( *bifIt ) ) );
+    _TRegion region = this->_Region( *bifIt, d * 1.5 );
+    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 );
+
+    } // rof
+
+  } // rof
 }
 
 // -------------------------------------------------------------------------
@@ -394,7 +513,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 ] ) ) );
+    // 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;
 
@@ -409,9 +531,37 @@ _Region( const TVertex& c, const double& r )
   // Prepare region and return it
   region.SetIndex( i0 );
   region.SetSize( size );
+
   return( region );
 }
 
+// -------------------------------------------------------------------------
+template< class I, class O >
+template< class _T > 
+typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
+TVertex fpa::Image::DijkstraWithEndPointDetection< I, O >::
+_MaxInRegion( const _T* image, const TVertex& v, const double& r )
+{
+  typedef itk::ImageRegionConstIteratorWithIndex< _T > _TIt;
+
+  _TIt iIt( image, this->_Region( v, r ) );
+  iIt.GoToBegin( );
+  TVertex max_vertex = iIt.GetIndex( );
+  typename _T::PixelType max_value = iIt.Get( );
+  for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
+  {
+    typename _T::PixelType value = iIt.Get( );
+    if( value > max_value )
+    {
+      max_value = value;
+      max_vertex = iIt.GetIndex( );
+
+    } // fi
+
+  } // rof
+  return( max_vertex );
+}
+
 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
 
 // eof - $RCSfile$