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   typename TSeedsInterface::TSeeds::const_iterator sIt = this->BeginSeeds( );
  69   for( ; sIt != this->EndSeeds( ); ++sIt )
  70     queues[ 0 ].push( _TNode( *sIt, frontId++ ) );
  71   unsigned int cur_queue = 0;
  72   unsigned int aux_queue = 1;
  73
  74   // Incremental strategy
  75   TInputValue upper = this->m_LowerThreshold + this->m_DeltaThreshold;
  76   this->m_Curve.clear( );
  77   unsigned long count = 0;
  78   while( upper <= this->m_UpperThreshold )
  79   {
  80     // Main growing strategy
  81     while( queues[ cur_queue ].size( ) > 0 )
  82     {
  83       // Get next candidate
  84       _TNode node = queues[ cur_queue ].front( );
  85       queues[ cur_queue ].pop( );
  86       if( this->_IsMarked( node.first ) )
  87         continue;
  88       this->_Mark( node.first, node.second );
  89
  90       // Apply inclusion predicate
  91       TInputValue value = this->_GetInputValue( node.first );
  92       bool in = ( ( this->m_LowerThreshold < value ) && ( value < upper ) );
  93       if( !in )
  94       {
  95         if( value < this->m_UpperThreshold )
  96           queues[ aux_queue ].push( node );
  97         this->_Mark( node.first, 0 );
  98         continue;
  99
 100       } // fi
 101
 102       // Ok, value lays inside region
 103       this->_SetOutputValue( node.first, this->m_Curve.size( ) + 1 );
 104       count++;
 105
 106       // Add neighborhood
 107       TVertices neighbors = this->_GetNeighbors( node.first );
 108       typename TVertices::const_iterator neighIt = neighbors.begin( );
 109       for( ; neighIt != neighbors.end( ); ++neighIt )
 110       {
 111         TVertex neigh = *neighIt;
 112         if( this->_IsMarked( neigh ) )
 113         {
 114           // Invoke stop at collisions
 115           unsigned long nColl = this->_Collisions( node.first, neigh );
 116           if(
 117             this->StopAtOneFront( ) &&
 118             this->GetNumberOfSeeds( ) > 1 &&
 119             nColl == 1
 120             )
 121           {
 122             while( queues[ 0 ].size( ) > 0 )
 123               queues[ 0 ].pop( );
 124             while( queues[ 1 ].size( ) > 0 )
 125               queues[ 1 ].pop( );
 126
 127           } // fi
 128         }
 129         else
 130           queues[ cur_queue ].push( _TNode( neigh, node.second ) );
 131
 132       } // rof
 133
 134     } // elihw
 135
 136     // Update curve
 137     if( this->m_Curve.size( ) > 0 )
 138     {
 139       if( this->m_Curve.back( ).YValue < count )
 140         this->m_Curve.push_back( TCurveData( upper, count ) );
 141       if( this->m_Curve.size( ) > 2 )
 142       {
 143         long j = this->m_Curve.size( ) - 2;
 144         double yp = double( this->m_Curve[ j + 1 ].YValue );
 145         double yn = double( this->m_Curve[ j - 1 ].YValue );
 146         double xp = double( this->m_Curve[ j + 1 ].XValue );
 147         double xn = double( this->m_Curve[ j - 1 ].XValue );
 148         this->m_Curve[ j ].Diff1 = ( yp - yn ) / ( xp - xn );
 149
 150       } // fi
 151     }
 152     else
 153       this->m_Curve.push_back( TCurveData( upper, count ) );
 154
 155     // Update queue
 156     cur_queue = aux_queue;
 157     aux_queue = ( cur_queue + 1 ) % 2;
 158
 159     // Update threshold
 160     upper += this->m_DeltaThreshold;
 161     progress.CompletedPixel( );
 162
 163   } // elihw
 164
 165   // Compute optimum threshold
 166   double dmax = -std::numeric_limits< double >::max( );
 167   long jmax = 0;
 168   for( long j = 1; j < this->m_Curve.size( ) - 1; ++j )
 169   {
 170     double dp = this->m_Curve[ j + 1 ].Diff1;
 171     double dn = this->m_Curve[ j - 1 ].Diff1;
 172     double xp = double( this->m_Curve[ j + 1 ].XValue );
 173     double xn = double( this->m_Curve[ j - 1 ].XValue );
 174     double d2 = ( dp - dn ) / ( xp - xn );
 175     if( d2 > dmax )
 176     {
 177       dmax = d2;
 178       jmax = j;
 179
 180     } // fi
 181
 182   } // rof
 183   this->m_OptimumThreshold = TOutputValue( jmax );
 184   this->_FreeMarks( );
 185 }
 186
 187 #endif // __fpa__Base__MoriRegionGrow__hxx__
 188
 189 // eof - $RCSfile$