+ const TInputImage* input = this->GetInput( );
+ TAuxiliaryImage* auxiliary = this->GetAuxiliaryImage( );
+ TRegion region = input->GetRequestedRegion( );
+
+ // 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, 0.7 );
+
+ // Configure outputs
+ auxiliary->SetBufferedRegion( region );
+ auxiliary->Allocate( );
+ auxiliary->FillBuffer( 0 );
+
+ // Init marks
+ typedef itk::Image< bool, TInputImage::ImageDimension > _TMarks;
+ typename _TMarks::Pointer marks = _TMarks::New( );
+ marks->CopyInformation( input );
+ marks->SetRequestedRegion( region );
+ marks->SetBufferedRegion( input->GetBufferedRegion( ) );
+ marks->Allocate( );
+ marks->FillBuffer( false );
+
+ // Prepare queues
+ std::queue< TIndex > queues[ 2 ];
+ queues[ 0 ].push( this->m_Seed );
+ unsigned int cur_queue = 0;
+ unsigned int aux_queue = 1;
+
+ // Incremental strategy
+ TInputPixel 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
+ TIndex node = queues[ cur_queue ].front( );
+ queues[ cur_queue ].pop( );
+ if( marks->GetPixel( node ) )
+ continue;
+ marks->SetPixel( node, true );
+
+ // Apply inclusion predicate
+ TInputPixel value = input->GetPixel( node );
+ bool in = ( ( this->m_LowerThreshold < value ) && ( value < upper ) );
+ if( !in )
+ {
+ if( value < this->m_UpperThreshold )
+ queues[ aux_queue ].push( node );
+ marks->SetPixel( node, false );
+ continue;
+
+ } // fi
+
+ // Ok, pixel lays inside region
+ auxiliary->SetPixel( node, this->m_Curve.size( ) + 1 );
+ count++;
+
+ // Add neighborhood
+ for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d )
+ {
+ TIndex neigh = node;
+ for( int i = -1; i <= 1; i += 2 )
+ {
+ neigh[ d ] = node[ d ] + i;
+ if( region.IsInside( neigh ) )
+ queues[ cur_queue ].push( neigh );
+
+ } // rof
+
+ } // 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( );
+ int jmax = 0;
+ for( int 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 = this->m_Curve[ jmax ].XValue;
+
+ // Final threshold
+ this->m_ThresholdFilter->SetInput( auxiliary );
+ this->m_ThresholdFilter->SetInsideValue( this->m_InsideValue );
+ this->m_ThresholdFilter->SetOutsideValue( this->m_OutsideValue );
+ this->m_ThresholdFilter->SetLowerThreshold( 1 );
+ this->m_ThresholdFilter->SetUpperThreshold( jmax );
+ this->m_ThresholdFilter->Update( );
+ this->GetOutput( )->Graft( this->m_ThresholdFilter->GetOutput( ) );