Segmentation completely debugged (at least for 10 images)

diff --git a/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx b/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx
index ff821034776c889ffcd7cf53accbb754a63c90cb..deef6de805a8e711164b8313a4e4319905a30614 100644 (file)
--- a/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx
+++ b/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx
@@ -42,14 +42,14 @@ int main( int argc, char* argv[] )
   {
     std::cerr
       << "Usage: " << argv[ 0 ]
-      << " input_image thr_0 thr_1 n_samples" << std::endl;
+      << " input_image thr_0 thr_1 step" << std::endl;
     return( 1 );
 
   } // fi
   std::string input_image_fn = argv[ 1 ];
   TPixel thr_0 = TPixel( std::atof( argv[ 2 ] ) );
   TPixel thr_1 = TPixel( std::atof( argv[ 3 ] ) );
-  unsigned int n_samples = std::atoi( argv[ 4 ] );
+  unsigned int step = std::atoi( argv[ 4 ] );
   
   // Read image
   TImageReader::Pointer input_image_reader = TImageReader::New( );
@@ -99,13 +99,13 @@ int main( int argc, char* argv[] )
 
   // Configure algorithm
   TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( );
-  algorithm->AddThresholds( thr_0, thr_1, n_samples );
+  algorithm->AddThresholds( thr_0, thr_1, step );
   algorithm->AddSeed( seed_idx, 0 );
   algorithm->AddObserver( itk::AnyEvent( ), obs );
   algorithm->ThrowEventsOn( );
   algorithm->SetInput( input_image );
   algorithm->SetNeighborhoodOrder( 1 );
-  algorithm->SetDerivativeThreshold( double( 3 ) );
+  algorithm->SetDifferenceThreshold( double( 3 ) );
   algorithm->Update( );
 
   // Let some interaction and close program

diff --git a/lib/fpa/Image/RegionGrowWithMultipleThresholds.h b/lib/fpa/Image/RegionGrowWithMultipleThresholds.h
index 5156b33dbd9bce73614bd7c28734fd4be82aa37b..409b18bb90895d8e65c4a70b303f604f6be3891e 100644 (file)
--- a/lib/fpa/Image/RegionGrowWithMultipleThresholds.h
+++ b/lib/fpa/Image/RegionGrowWithMultipleThresholds.h
@@ -40,8 +40,8 @@ namespace fpa
       itkNewMacro( Self );
       itkTypeMacro( RegionGrowWithMultipleThresholds, RegionGrow );
 
-      itkGetConstMacro( DerivativeThreshold, double );
-      itkSetMacro( DerivativeThreshold, double );
+      itkGetConstMacro( DifferenceThreshold, double );
+      itkSetMacro( DifferenceThreshold, double );
 
     public:
       void AddThreshold( const TPixel& v );
@@ -59,6 +59,8 @@ namespace fpa
       virtual void _AfterMainLoop( );
       virtual void _AfterLoop( );
       virtual bool _UpdateResult( _TNode& n );
+      virtual void _Mark( const _TNode& n );
+      virtual bool _CheckStopCondition( );
 
     private:
       RegionGrowWithMultipleThresholds( const Self& ); // Not impl.
@@ -66,8 +68,12 @@ namespace fpa
 
     protected:
       TThresholds m_Thresholds;
-      double m_DerivativeThreshold;
+      double m_DifferenceThreshold;
       THistogram m_Histogram;
+      unsigned long m_TotalCount;
+      double m_LastDiff;
+      bool m_StopForced;
+      TPixel m_StopThreshold;
     };
 
   } // ecapseman

diff --git a/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx b/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx
index 0cb1b66598aa225f33799855c44a4435ed376a20..387a833bcab9bda89bba8d05e36eb9ec76fdeeae 100644 (file)
--- a/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx
+++ b/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx
@@ -2,6 +2,7 @@
 #define __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
 
 #include <limits>
+#include <itkConstNeighborhoodIterator.h>
 
 // -------------------------------------------------------------------------
 template< class I >
@@ -17,9 +18,8 @@ template< class I >
 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
 AddThresholds( const TPixel& t0, const TPixel& t1, const unsigned int& s )
 {
-  double o = ( double( t1 ) - double( t0 ) ) / ( double( s ) - double( 1 ) );
-  for( unsigned int i = 0; i < s; i++ )
-    this->AddThreshold( TPixel( ( double( i ) * o ) + double( t0 ) ) );
+  for( TPixel t = t0; t <= t1; t += TPixel( s ) )
+    this->AddThreshold( t );
 }
 
 // -------------------------------------------------------------------------
@@ -27,7 +27,11 @@ template< class I >
 fpa::Image::RegionGrowWithMultipleThresholds< I >::
 RegionGrowWithMultipleThresholds( )
   : Superclass( ),
-    m_DerivativeThreshold( double( 3 ) )
+    m_DifferenceThreshold( double( 3 ) ),
+    m_TotalCount( 0 ),
+    m_LastDiff( double( 0 ) ),
+    m_StopForced( false ),
+    m_StopThreshold( TPixel( 0 ) )
 {
 }
 
@@ -45,10 +49,17 @@ _BeforeMainLoop( )
 {
   const I* img = this->GetInput( );
 
+  // Clear all states
   this->ClearMembershipFunctions( );
+  this->m_Histogram.clear( );
+  this->m_TotalCount = img->GetLargestPossibleRegion( ).GetNumberOfPixels( );
+  this->m_LastDiff = double( -1 );
+  this->m_StopForced = false;
+
+  // Initialize thresholding functors
   typename TThresholds::const_iterator tIt = this->m_Thresholds.begin( );
-  TPixel min_thr = *tIt;
-  for( ++tIt; tIt != this->m_Thresholds.end( ); ++tIt )
+  TPixel min_thr = std::numeric_limits< TPixel >::min( );
+  for( ; tIt != this->m_Thresholds.end( ); ++tIt )
   {
     typename TFunction::Pointer function = TFunction::New( );
     function->SetInputImage( img );
@@ -58,9 +69,48 @@ _BeforeMainLoop( )
 
   } // rof
 
+  // Correct seeds
+  typename I::SizeType radius;
+  radius.Fill( 3 );
+  itk::ConstNeighborhoodIterator< I > it(
+    radius, img, img->GetRequestedRegion( )
+    );
+  for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
+  {
+    _TNode seed = this->m_Seeds[ s ];
+    it.SetLocation( seed.Vertex );
+
+    typename I::SizeType size = it.GetSize( );
+    unsigned int nN = 1;
+    for( unsigned int d = 0; d < I::ImageDimension; ++d )
+      nN *= size[ d ];
+    TPixel min_value = img->GetPixel( seed.Vertex );
+    for( unsigned int i = 0; i < nN; ++i )
+    {
+      if( it.GetPixel( i ) < min_value )
+      {
+        seed.Vertex = it.GetIndex( i );
+        min_value = it.GetPixel( i );
+
+      } // fi
+
+    } // rof
+    this->m_Seeds[ s ] = seed;
+
+  } // rof
+
+  // Continue all initializations
   this->Superclass::_BeforeMainLoop( );
 }
 
+// -------------------------------------------------------------------------
+template< class I >
+void fpa::Image::RegionGrowWithMultipleThresholds< I >::
+_AfterMainLoop( )
+{
+  this->Superclass::_AfterMainLoop( );
+}
+
 // -------------------------------------------------------------------------
 template< class I >
 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
@@ -70,36 +120,111 @@ _AfterLoop( )
   {
     TFunction* f =
       dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
-    // std::cout << f->GetUpperThreshold( ) << " " << this->m_Marks.size( ) << std::endl;
+    this->m_Histogram[ f->GetUpperThreshold( ) ] = this->m_Marks.size( );
+
+    /* TODO: code to print the signal for demo purposes
+       std::cout
+       << f->GetUpperThreshold( ) << " "
+       << this->m_Marks.size( )
+       << std::endl;
+    */
+
+    // Get the two last complete count, if any, and compute last
+    // finite difference
+    if( this->m_Histogram.size( ) > 1 )
+    {
+      typename THistogram::const_reverse_iterator hIt =
+        this->m_Histogram.rbegin( );
+      typename THistogram::const_reverse_iterator gIt = hIt;
+      gIt++;
 
+      this->m_LastDiff  = double( hIt->second ) - double( gIt->second );
+      this->m_LastDiff /= double( hIt->first )  - double( gIt->first );
+      this->m_LastDiff *= m_DifferenceThreshold;
+    }
+    else
+      this->m_LastDiff = double( -1 );
+    
   } // fi
   this->Superclass::_AfterLoop( );
 }
 
-// -------------------------------------------------------------------------
-template< class I >
-void fpa::Image::RegionGrowWithMultipleThresholds< I >::
-_AfterMainLoop( )
-{
-  this->Superclass::_AfterMainLoop( );
-}
-
 // -------------------------------------------------------------------------
 template< class I >
 bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
 _UpdateResult( _TNode& n )
 {
-  #error ACA VOY -> explicar esto!!!! cambiar salida por una especie de curva de nivel
+  // Here, the output is changed to the threshold used to compute the
+  // growing condition of the actual vertex n
   bool ret = this->Superclass::_UpdateResult( n );
   if( this->m_ActualFunction != this->m_Functions.end( ) )
   {
     TFunction* f =
       dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
-    this->GetOutput( )->SetPixel( n.Vertex, f->GetUpperThreshold( ) );
+    if( f != NULL )
+      this->GetOutput( )->SetPixel( n.Vertex, f->GetUpperThreshold( ) );
 
   } // fi
   return( ret );
 }
+
+// -------------------------------------------------------------------------
+template< class I >
+void fpa::Image::RegionGrowWithMultipleThresholds< I >::
+_Mark( const _TNode& n )
+{
+  this->Superclass::_Mark( n );
+
+  // Check if the histogram's support is enough
+  if( this->m_Histogram.size( ) < 2 )
+    return;
+  typename THistogram::const_reverse_iterator hIt =
+    this->m_Histogram.rbegin( );
+
+  // Get the actual function
+  TFunction* f =
+    dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
+  if( f == NULL )
+    return;
+
+  // Previous pixel count
+  unsigned long prev_aCount = hIt->second;
+  double prev_p = double( prev_aCount ) / double( this->m_TotalCount );
+
+  // Actual pixel count
+  unsigned long aCount = this->m_Marks.size( );
+  double p = double( aCount ) / double( this->m_TotalCount );
+
+  // Loop over, at least, 0.1% of the total number of pixels before
+  // performing stop analysis 
+  if( double( 1e-3 ) < p && double( 1e-3 ) < prev_p )
+  {
+    // Does the difference is worthy to be analyzed?
+    if( aCount > hIt->second )
+    {
+      // Compute finite difference and compare it to the previous one
+      double diff = double( aCount ) - double( hIt->second );
+      diff /= double( f->GetUpperThreshold( ) ) - double( hIt->first );
+      if( diff > this->m_LastDiff )
+      {
+        this->m_StopForced = true;
+        this->m_StopThreshold = hIt->first;
+
+      } // fi
+
+    } // fi
+    
+  } // fi
+}
+
+// -------------------------------------------------------------------------
+template< class I >
+bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
+_CheckStopCondition( )
+{
+  return( this->m_StopForced );
+}
+
 #endif // __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
 
 // eof - $RCSfile$