Skeleton reconstruction added

[FrontAlgorithms.git] / lib / fpa / Image / DijkstraWithEndPointDetection.hxx

diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx
index 76695c267e2c888faf0d97b10af17cdf3e941fd3..dc51a6442c0a66f7ed77663c25e8bb3c161923a7 100644 (file)
--- 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?
@@ -507,7 +515,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;