]> Creatis software - FrontAlgorithms.git/blobdiff - lib/fpa/Base/MoriRegionGrow.hxx
...
[FrontAlgorithms.git] / lib / fpa / Base / MoriRegionGrow.hxx
index 9e1c300201888b4bf81e6bbf3e373f16e32f8fbb..3b7d83949c9cf3262c5f5c1590c433730aecdd2b 100644 (file)
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
 #ifndef __fpa__Base__MoriRegionGrow__hxx__
 #define __fpa__Base__MoriRegionGrow__hxx__
 
+#include <queue>
+#include <itkProgressReporter.h>
+
 // -------------------------------------------------------------------------
-template< class _TSuperclass >
-fpa::Base::MoriRegionGrow< _TSuperclass >::
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
+void
+fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
+SetThresholdRange(
+  const TInputValue& lower, const TInputValue& upper,
+  const TInputValue& delta
+  )
+{
+  this->SetLowerThreshold( lower );
+  this->SetUpperThreshold( upper );
+  this->SetDeltaThreshold( delta );
+}
+
+// -------------------------------------------------------------------------
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
+fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
 MoriRegionGrow( )
-  : Superclass( )
+  : Superclass( ),
+    _TMarksInterface( this ),
+    _TSeedsInterface( this )
 {
+  this->m_UpperThreshold = std::numeric_limits< TInputValue >::max( );
+  if( std::numeric_limits< TInputValue >::is_integer )
+    this->m_LowerThreshold = std::numeric_limits< TInputValue >::min( );
+  else
+    this->m_LowerThreshold = -this->m_UpperThreshold;
+  this->m_DeltaThreshold = TInputValue( 1 );
 }
 
 // -------------------------------------------------------------------------
-template< class _TSuperclass >
-fpa::Base::MoriRegionGrow< _TSuperclass >::
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
+fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
 ~MoriRegionGrow( )
 {
 }
 
 // -------------------------------------------------------------------------
-template< class _TSuperclass >
-bool fpa::Base::MoriRegionGrow< _TSuperclass >::
-_UpdateValue( _TQueueNode& v, const _TQueueNode& p )
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
+void
+fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
+GenerateData( )
 {
-  bool ret = this->Superclass::_UpdateValue( v, p );
-  if( !ret )
-    this->m_AuxilaryQueue.push( v );
-  return( ret );
-}
+  // Prepare progress counter
+  double prog_size = double( this->m_UpperThreshold - this->m_LowerThreshold );
+  prog_size /= double( this->m_DeltaThreshold );
+  itk::ProgressReporter progress( this, 0, prog_size, 100, 0, 1 );
+
+  // Init objects
+  this->_ConfigureOutputs( std::numeric_limits< TOutputValue >::max( ) );
+  this->_InitMarks( this->GetNumberOfSeeds( ) );
+
+  // Init queues
+  typedef std::pair< TVertex, unsigned long > _TNode;
+  std::queue< _TNode > queues[ 2 ];
+  unsigned long frontId = 1;
+  typename TSeedsInterface::TSeeds::const_iterator sIt = this->BeginSeeds( );
+  for( ; sIt != this->EndSeeds( ); ++sIt )
+    queues[ 0 ].push( _TNode( *sIt, frontId++ ) );
+  unsigned int cur_queue = 0;
+  unsigned int aux_queue = 1;
+
+  // Incremental strategy
+  TInputValue upper = this->m_LowerThreshold + this->m_DeltaThreshold;
+  this->m_Curve.clear( );
+  unsigned long count = 0;
+  while( upper <= this->m_UpperThreshold )
+  {
+    // Main growing strategy
+    while( queues[ cur_queue ].size( ) > 0 )
+    {
+      // Get next candidate
+      _TNode node = queues[ cur_queue ].front( );
+      queues[ cur_queue ].pop( );
+      if( this->_IsMarked( node.first ) )
+        continue;
+      this->_Mark( node.first, node.second );
+
+      // Apply inclusion predicate
+      TInputValue value = this->_GetInputValue( node.first );
+      bool in = ( ( this->m_LowerThreshold < value ) && ( value < upper ) );
+      if( !in )
+      {
+        if( value < this->m_UpperThreshold )
+          queues[ aux_queue ].push( node );
+        this->_Mark( node.first, 0 );
+        continue;
+
+      } // fi
+
+      // Ok, value lays inside region
+      this->_SetOutputValue( node.first, this->m_Curve.size( ) + 1 );
+      count++;
 
+      // Add neighborhood
+      TVertices neighbors = this->_GetNeighbors( node.first );
+      typename TVertices::const_iterator neighIt = neighbors.begin( );
+      for( ; neighIt != neighbors.end( ); ++neighIt )
+      {
+        TVertex neigh = *neighIt;
+        if( this->_IsMarked( neigh ) )
+        {
+          // Invoke stop at collisions
+          unsigned long nColl = this->_Collisions( node.first, neigh );
+          if(
+            this->StopAtOneFront( ) &&
+            this->GetNumberOfSeeds( ) > 1 &&
+            nColl == 1
+            )
+          {
+            while( queues[ 0 ].size( ) > 0 )
+              queues[ 0 ].pop( );
+            while( queues[ 1 ].size( ) > 0 )
+              queues[ 1 ].pop( );
+
+          } // fi
+        }
+        else
+          queues[ cur_queue ].push( _TNode( neigh, node.second ) );
+
+      } // rof
+
+    } // elihw
+
+    // Update curve
+    if( this->m_Curve.size( ) > 0 )
+    {
+      if( this->m_Curve.back( ).YValue < count )
+        this->m_Curve.push_back( TCurveData( upper, count ) );
+      if( this->m_Curve.size( ) > 2 )
+      {
+        long j = this->m_Curve.size( ) - 2;
+        double yp = double( this->m_Curve[ j + 1 ].YValue );
+        double yn = double( this->m_Curve[ j - 1 ].YValue );
+        double xp = double( this->m_Curve[ j + 1 ].XValue );
+        double xn = double( this->m_Curve[ j - 1 ].XValue );
+        this->m_Curve[ j ].Diff1 = ( yp - yn ) / ( xp - xn );
+
+      } // fi
+    }
+    else
+      this->m_Curve.push_back( TCurveData( upper, count ) );
+
+    // Update queue
+    cur_queue = aux_queue;
+    aux_queue = ( cur_queue + 1 ) % 2;
+
+    // Update threshold
+    upper += this->m_DeltaThreshold;
+    progress.CompletedPixel( );
+
+  } // elihw
+
+  // Compute optimum threshold
+  double dmax = -std::numeric_limits< double >::max( );
+  long jmax = 0;
+  for( long j = 1; j < this->m_Curve.size( ) - 1; ++j )
+  {
+    double dp = this->m_Curve[ j + 1 ].Diff1;
+    double dn = this->m_Curve[ j - 1 ].Diff1;
+    double xp = double( this->m_Curve[ j + 1 ].XValue );
+    double xn = double( this->m_Curve[ j - 1 ].XValue );
+    double d2 = ( dp - dn ) / ( xp - xn );
+    if( d2 > dmax )
+    {
+      dmax = d2;
+      jmax = j;
+
+    } // fi
+
+  } // rof
+  this->m_OptimumThreshold = TOutputValue( jmax );
+  this->_FreeMarks( );
+}
 
 #endif // __fpa__Base__MoriRegionGrow__hxx__