libs/fpa/Image/MoriRegionGrow.hxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #ifndef __fpa__Image__MoriRegionGrow__hxx__
   7 #define __fpa__Image__MoriRegionGrow__hxx__
   8
   9 #include <queue>
  10 #include <itkProgressReporter.h>
  11
  12 // -------------------------------------------------------------------------
  13 template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
  14 typename
  15 fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
  16 TAuxiliaryImage*
  17 fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
  18 GetAuxiliaryImage( )
  19 {
  20   return(
  21     dynamic_cast< TAuxiliaryImage* >(
  22       this->itk::ProcessObject::GetOutput( 1 )
  23       )
  24     );
  25 }
  26
  27 // -------------------------------------------------------------------------
  28 template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
  29 const typename
  30 fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
  31 TAuxiliaryImage*
  32 fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
  33 GetAuxiliaryImage( ) const
  34 {
  35   return(
  36     dynamic_cast< const TAuxiliaryImage* >(
  37       this->itk::ProcessObject::GetOutput( 1 )
  38       )
  39     );
  40 }
  41
  42 // -------------------------------------------------------------------------
  43 template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
  44 void
  45 fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
  46 SetThresholdRange(
  47   const TInputPixel& lower, const TInputPixel& upper,
  48   const TInputPixel& delta
  49   )
  50 {
  51   this->SetLowerThreshold( lower );
  52   this->SetUpperThreshold( upper );
  53   this->SetDeltaThreshold( delta );
  54 }
  55
  56 // -------------------------------------------------------------------------
  57 template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
  58 fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
  59 MoriRegionGrow( )
  60   : Superclass( )
  61 {
  62   this->m_InsideValue = TOutputPixel( 1 );
  63   this->m_OutsideValue = TOutputPixel( 0 );
  64   this->m_UpperThreshold = std::numeric_limits< TInputPixel >::max( );
  65   if( std::numeric_limits< TInputPixel >::is_integer )
  66     this->m_LowerThreshold = std::numeric_limits< TInputPixel >::min( );
  67   else
  68     this->m_LowerThreshold = -this->m_UpperThreshold;
  69   this->m_DeltaThreshold = TInputPixel( 1 );
  70
  71   this->SetNumberOfRequiredOutputs( 2 );
  72   this->SetNthOutput( 0, TOutputImage::New( ) );
  73   this->SetNthOutput( 1, TAuxiliaryImage::New( ) );
  74
  75   this->m_ThresholdFilter = TThresholdFilter::New( );
  76 }
  77
  78 // -------------------------------------------------------------------------
  79 template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
  80 fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
  81 ~MoriRegionGrow( )
  82 {
  83 }
  84
  85 // -------------------------------------------------------------------------
  86 template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
  87 void
  88 fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
  89 GenerateInputRequestedRegion( )
  90 {
  91   this->Superclass::GenerateInputRequestedRegion( );
  92   if( this->GetInput( ) )
  93   {
  94     TInputImage* input =
  95       const_cast< TInputImage* >( this->GetInput( ) );
  96     input->SetRequestedRegionToLargestPossibleRegion( );
  97
  98   } // fi
  99 }
 100
 101 // -------------------------------------------------------------------------
 102 template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
 103 void
 104 fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
 105 EnlargeOutputRequestedRegion( itk::DataObject* output )
 106 {
 107   this->Superclass::EnlargeOutputRequestedRegion( output );
 108   output->SetRequestedRegionToLargestPossibleRegion( );
 109 }
 110
 111 // -------------------------------------------------------------------------
 112 template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
 113 void fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
 114 GenerateData( )
 115 {
 116   const TInputImage* input = this->GetInput( );
 117   TAuxiliaryImage* auxiliary = this->GetAuxiliaryImage( );
 118   TRegion region = input->GetRequestedRegion( );
 119
 120   // Prepare progress counter
 121   double prog_size = double( this->m_UpperThreshold - this->m_LowerThreshold );
 122   prog_size /= double( this->m_DeltaThreshold );
 123   itk::ProgressReporter progress( this, 0, prog_size, 100, 0, 0.7 );
 124
 125   // Configure outputs
 126   auxiliary->SetBufferedRegion( region );
 127   auxiliary->Allocate( );
 128   auxiliary->FillBuffer( 0 );
 129
 130   // Init marks
 131   typedef itk::Image< bool, TInputImage::ImageDimension > _TMarks;
 132   typename _TMarks::Pointer marks = _TMarks::New( );
 133   marks->CopyInformation( input );
 134   marks->SetRequestedRegion( region );
 135   marks->SetBufferedRegion( input->GetBufferedRegion( ) );
 136   marks->Allocate( );
 137   marks->FillBuffer( false );
 138
 139   // Prepare queues
 140   std::queue< TIndex > queues[ 2 ];
 141   queues[ 0 ].push( this->m_Seed );
 142   unsigned int cur_queue = 0;
 143   unsigned int aux_queue = 1;
 144
 145   // Incremental strategy
 146   TInputPixel upper = this->m_LowerThreshold + this->m_DeltaThreshold;
 147   this->m_Curve.clear( );
 148   unsigned long count = 0;
 149   while( upper <= this->m_UpperThreshold )
 150   {
 151     // Main growing strategy
 152     while( queues[ cur_queue ].size( ) > 0 )
 153     {
 154       // Get next candidate
 155       TIndex node = queues[ cur_queue ].front( );
 156       queues[ cur_queue ].pop( );
 157       if( marks->GetPixel( node ) )
 158         continue;
 159       marks->SetPixel( node, true );
 160
 161       // Apply inclusion predicate
 162       TInputPixel value = input->GetPixel( node );
 163       bool in = ( ( this->m_LowerThreshold < value ) && ( value < upper ) );
 164       if( !in )
 165       {
 166         if( value < this->m_UpperThreshold )
 167           queues[ aux_queue ].push( node );
 168         marks->SetPixel( node, false );
 169         continue;
 170
 171       } // fi
 172
 173       // Ok, pixel lays inside region
 174       auxiliary->SetPixel( node, this->m_Curve.size( ) + 1 );
 175       count++;
 176
 177       // Add neighborhood
 178       for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d )
 179       {
 180         TIndex neigh = node;
 181         for( int i = -1; i <= 1; i += 2 )
 182         {
 183           neigh[ d ] = node[ d ] + i;
 184           if( region.IsInside( neigh ) )
 185             queues[ cur_queue ].push( neigh );
 186
 187         } // rof
 188
 189       } // rof
 190
 191     } // elihw
 192
 193     // Update curve
 194     if( this->m_Curve.size( ) > 0 )
 195     {
 196       if( this->m_Curve.back( ).YValue < count )
 197         this->m_Curve.push_back( TCurveData( upper, count ) );
 198       if( this->m_Curve.size( ) > 2 )
 199       {
 200         long j = this->m_Curve.size( ) - 2;
 201         double yp = double( this->m_Curve[ j + 1 ].YValue );
 202         double yn = double( this->m_Curve[ j - 1 ].YValue );
 203         double xp = double( this->m_Curve[ j + 1 ].XValue );
 204         double xn = double( this->m_Curve[ j - 1 ].XValue );
 205         this->m_Curve[ j ].Diff1 = ( yp - yn ) / ( xp - xn );
 206
 207       } // fi
 208     }
 209     else
 210       this->m_Curve.push_back( TCurveData( upper, count ) );
 211
 212     // Update queue
 213     cur_queue = aux_queue;
 214     aux_queue = ( cur_queue + 1 ) % 2;
 215
 216     // Update threshold
 217     upper += this->m_DeltaThreshold;
 218     progress.CompletedPixel( );
 219
 220   } // elihw
 221
 222   // Compute optimum threshold
 223   double dmax = -std::numeric_limits< double >::max( );
 224   int jmax = 0;
 225   for( int j = 1; j < this->m_Curve.size( ) - 1; ++j )
 226   {
 227     double dp = this->m_Curve[ j + 1 ].Diff1;
 228     double dn = this->m_Curve[ j - 1 ].Diff1;
 229     double xp = double( this->m_Curve[ j + 1 ].XValue );
 230     double xn = double( this->m_Curve[ j - 1 ].XValue );
 231     double d2 = ( dp - dn ) / ( xp - xn );
 232     if( d2 > dmax )
 233     {
 234       dmax = d2;
 235       jmax = j;
 236
 237     } // fi
 238
 239   } // rof
 240   this->m_OptimumThreshold = this->m_Curve[ jmax ].XValue;
 241
 242   // Final threshold
 243   this->m_ThresholdFilter->SetInput( auxiliary );
 244   this->m_ThresholdFilter->SetInsideValue( this->m_InsideValue );
 245   this->m_ThresholdFilter->SetOutsideValue( this->m_OutsideValue );
 246   this->m_ThresholdFilter->SetLowerThreshold( 1 );
 247   this->m_ThresholdFilter->SetUpperThreshold( jmax );
 248   this->m_ThresholdFilter->Update( );
 249   this->GetOutput( )->Graft( this->m_ThresholdFilter->GetOutput( ) );
 250 }
 251
 252 #endif // __fpa__Image__MoriRegionGrow__hxx__
 253
 254 // eof - $RCSfile$