lib/fpa/Base/MoriRegionGrow.hxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #ifndef __fpa__Base__MoriRegionGrow__hxx__
   7 #define __fpa__Base__MoriRegionGrow__hxx__
   8
   9 #include <queue>
  10 #include <itkProgressReporter.h>
  11
  12 // -------------------------------------------------------------------------
  13 template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
  14 void
  15 fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
  16 SetThresholdRange(
  17   const TInputValue& lower, const TInputValue& upper,
  18   const TInputValue& delta
  19   )
  20 {
  21   this->SetLowerThreshold( lower );
  22   this->SetUpperThreshold( upper );
  23   this->SetDeltaThreshold( delta );
  24 }
  25
  26 // -------------------------------------------------------------------------
  27 template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
  28 fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
  29 MoriRegionGrow( )
  30   : Superclass( ),
  31     _TMarksInterface( this ),
  32     _TSeedsInterface( this )
  33 {
  34   this->m_UpperThreshold = std::numeric_limits< TInputValue >::max( );
  35   if( std::numeric_limits< TInputValue >::is_integer )
  36     this->m_LowerThreshold = std::numeric_limits< TInputValue >::min( );
  37   else
  38     this->m_LowerThreshold = -this->m_UpperThreshold;
  39   this->m_DeltaThreshold = TInputValue( 1 );
  40 }
  41
  42 // -------------------------------------------------------------------------
  43 template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
  44 fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
  45 ~MoriRegionGrow( )
  46 {
  47 }
  48
  49 // -------------------------------------------------------------------------
  50 template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
  51 void
  52 fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
  53 GenerateData( )
  54 {
  55   // Prepare progress counter
  56   double prog_size = double( this->m_UpperThreshold - this->m_LowerThreshold );
  57   prog_size /= double( this->m_DeltaThreshold );
  58   itk::ProgressReporter progress( this, 0, prog_size, 100, 0, 1 );
  59
  60   // Init objects
  61   this->_ConfigureOutputs( std::numeric_limits< TOutputValue >::max( ) );
  62   this->_InitMarks( this->GetNumberOfSeeds( ) );
  63
  64   // Init queues
  65   typedef std::pair< TVertex, unsigned long > _TNode;
  66   std::queue< _TNode > queues[ 2 ];
  67   unsigned long frontId = 1;
  68   for( TVertex seed: this->GetSeeds( ) )
  69     queues[ 0 ].push( _TNode( seed, frontId++ ) );
  70   unsigned int cur_queue = 0;
  71   unsigned int aux_queue = 1;
  72
  73   // Incremental strategy
  74   TInputValue upper = this->m_LowerThreshold + this->m_DeltaThreshold;
  75   this->m_Curve.clear( );
  76   unsigned long count = 0;
  77   while( upper <= this->m_UpperThreshold )
  78   {
  79     // Main growing strategy
  80     while( queues[ cur_queue ].size( ) > 0 )
  81     {
  82       // Get next candidate
  83       _TNode node = queues[ cur_queue ].front( );
  84       queues[ cur_queue ].pop( );
  85       if( this->_IsMarked( node.first ) )
  86         continue;
  87       this->_Mark( node.first, node.second );
  88
  89       // Apply inclusion predicate
  90       TInputValue value = this->_GetInputValue( node.first );
  91       bool in = ( ( this->m_LowerThreshold < value ) && ( value < upper ) );
  92       if( !in )
  93       {
  94         if( value < this->m_UpperThreshold )
  95           queues[ aux_queue ].push( node );
  96         this->_Mark( node.first, 0 );
  97         continue;
  98
  99       } // fi
 100
 101       // Ok, value lays inside region
 102       this->_SetOutputValue( node.first, this->m_Curve.size( ) + 1 );
 103       count++;
 104
 105       // Add neighborhood
 106       TVertices neighbors = this->_GetNeighbors( node.first );
 107       for( TVertex neigh: neighbors )
 108       {
 109         if( this->_IsMarked( neigh ) )
 110         {
 111           // Invoke stop at collisions
 112           unsigned long nColl = this->_Collisions( node.first, neigh );
 113           if(
 114             this->StopAtOneFront( ) &&
 115             this->GetNumberOfSeeds( ) > 1 &&
 116             nColl == 1
 117             )
 118           {
 119             while( queues[ 0 ].size( ) > 0 )
 120               queues[ 0 ].pop( );
 121             while( queues[ 1 ].size( ) > 0 )
 122               queues[ 1 ].pop( );
 123
 124           } // fi
 125         }
 126         else
 127           queues[ cur_queue ].push( _TNode( neigh, node.second ) );
 128
 129       } // rof
 130
 131     } // elihw
 132
 133     // Update curve
 134     if( this->m_Curve.size( ) > 0 )
 135     {
 136       if( this->m_Curve.back( ).YValue < count )
 137         this->m_Curve.push_back( TCurveData( upper, count ) );
 138       if( this->m_Curve.size( ) > 2 )
 139       {
 140         long j = this->m_Curve.size( ) - 2;
 141         double yp = double( this->m_Curve[ j + 1 ].YValue );
 142         double yn = double( this->m_Curve[ j - 1 ].YValue );
 143         double xp = double( this->m_Curve[ j + 1 ].XValue );
 144         double xn = double( this->m_Curve[ j - 1 ].XValue );
 145         this->m_Curve[ j ].Diff1 = ( yp - yn ) / ( xp - xn );
 146
 147       } // fi
 148     }
 149     else
 150       this->m_Curve.push_back( TCurveData( upper, count ) );
 151
 152     // Update queue
 153     cur_queue = aux_queue;
 154     aux_queue = ( cur_queue + 1 ) % 2;
 155
 156     // Update threshold
 157     upper += this->m_DeltaThreshold;
 158     progress.CompletedPixel( );
 159
 160   } // elihw
 161
 162   // Compute optimum threshold
 163   double dmax = -std::numeric_limits< double >::max( );
 164   long jmax = 0;
 165   for( long j = 1; j < this->m_Curve.size( ) - 1; ++j )
 166   {
 167     double dp = this->m_Curve[ j + 1 ].Diff1;
 168     double dn = this->m_Curve[ j - 1 ].Diff1;
 169     double xp = double( this->m_Curve[ j + 1 ].XValue );
 170     double xn = double( this->m_Curve[ j - 1 ].XValue );
 171     double d2 = ( dp - dn ) / ( xp - xn );
 172     if( d2 > dmax )
 173     {
 174       dmax = d2;
 175       jmax = j;
 176
 177     } // fi
 178
 179   } // rof
 180   this->m_OptimumThreshold = TOutputValue( jmax );
 181   this->_FreeMarks( );
 182 }
 183
 184 #endif // __fpa__Base__MoriRegionGrow__hxx__
 185
 186 // eof - $RCSfile$