X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FImage%2FDijkstraWithEndPointDetection.hxx;h=cb417f027092de357d1353ff18b4cd367f860729;hb=6468f1fda5fed2788fbaef1a7925c91ecff83d13;hp=76695c267e2c888faf0d97b10af17cdf3e941fd3;hpb=b4ed6ddfa7e90e892f07cad9a760961bd4e84e6b;p=FrontAlgorithms.git

diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx
index 76695c2..cb417f0 100644
--- a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx
+++ b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx
@@ -170,7 +170,10 @@ GraftBranches( itk::DataObject* obj )
 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->m_EndPointsIndex = this->m_LabelImageIndex + 1;
@@ -203,39 +206,43 @@ template< class I, class O >
 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 _BeforeGenerateData( )
 {
-  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 )
+  if( this->m_CorrectSeeds )
   {
-    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 )
+    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 )
     {
-      if( iIt.Get( ) > max_value )
+      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 )
       {
-        this->m_SeedVertices[ s ] = iIt.GetIndex( );
-        max_value = iIt.Get( );
+        if( iIt.Get( ) > max_value )
+        {
+          this->m_SeedVertices[ s ] = iIt.GetIndex( );
+          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.erase( seed );
-    n.Parent = this->m_SeedVertices[ s ];
-    this->m_Seeds[ this->m_SeedVertices[ s ] ] = n;
 
-  } // rof
+  } // fi
 
   // End initialization
   this->Superclass::_BeforeGenerateData( );
@@ -312,7 +319,8 @@ _EndPointsAndBifurcations( )
 
     // Correct it to nearest start candidate (high distance value)
     double vr = std::sqrt( double( input->GetPixel( v ) ) );
-    v = this->_MaxInRegion( input, v, vr * double( 1.5 ) );
+    if( this->m_CorrectEndPoints )
+      v = this->_MaxInRegion( input, v, vr * double( 1.5 ) );
 
     // Now, check for real marking conditions
     //   1. Has it been visited by dijkstra?
@@ -420,6 +428,7 @@ _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 );
@@ -465,28 +474,26 @@ _LabelAll( )
   this->m_NumberOfBranches = actual_label;
 
   // 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 );
-
-     } // rof
-
-     } // rof
-  */
+  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
 }
 
 // -------------------------------------------------------------------------
@@ -507,7 +514,9 @@ _Region( const TVertex& c, const double& r )
   for( unsigned int d = 0; d < I::ImageDimension; ++d )
   {
     // NOTE: 3 is a minimum neighborhood size
-    long s = long( std::ceil( r / double( spac[ d ] ) ) ) + 3;
+    long s =
+      long( std::ceil( r / double( spac[ d ] ) ) ) +
+      long( this->m_SafetyNeighborhoodSize );
     i0[ d ] = c[ d ] - s;
     i1[ d ] = c[ d ] + s;