-// =========================================================================
-// @author Leonardo Florez Valencia
-// @email florez-l@javeriana.edu.co
-// =========================================================================
-
-#ifndef __fpa__Image__MoriRegionGrow__hxx__
-#define __fpa__Image__MoriRegionGrow__hxx__
-
-#include <queue>
-#include <itkProgressReporter.h>
-
-// -------------------------------------------------------------------------
-template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
-typename
-fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-TAuxiliaryImage*
-fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-GetAuxiliaryImage( )
-{
- return(
- dynamic_cast< TAuxiliaryImage* >(
- this->itk::ProcessObject::GetOutput( 1 )
- )
- );
-}
-
-// -------------------------------------------------------------------------
-template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
-const typename
-fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-TAuxiliaryImage*
-fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-GetAuxiliaryImage( ) const
-{
- return(
- dynamic_cast< const TAuxiliaryImage* >(
- this->itk::ProcessObject::GetOutput( 1 )
- )
- );
-}
-
-// -------------------------------------------------------------------------
-template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
-void
-fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-SetThresholdRange(
- const TInputPixel& lower, const TInputPixel& upper,
- const TInputPixel& delta
- )
-{
- this->SetLowerThreshold( lower );
- this->SetUpperThreshold( upper );
- this->SetDeltaThreshold( delta );
-}
-
-// -------------------------------------------------------------------------
-template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
-fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-MoriRegionGrow( )
- : Superclass( )
-{
- this->m_InsideValue = TOutputPixel( 1 );
- this->m_OutsideValue = TOutputPixel( 0 );
- this->m_UpperThreshold = std::numeric_limits< TInputPixel >::max( );
- if( std::numeric_limits< TInputPixel >::is_integer )
- this->m_LowerThreshold = std::numeric_limits< TInputPixel >::min( );
- else
- this->m_LowerThreshold = -this->m_UpperThreshold;
- this->m_DeltaThreshold = TInputPixel( 1 );
-
- this->SetNumberOfRequiredOutputs( 2 );
- this->SetNthOutput( 0, TOutputImage::New( ) );
- this->SetNthOutput( 1, TAuxiliaryImage::New( ) );
-
- this->m_ThresholdFilter = TThresholdFilter::New( );
-}
-
-// -------------------------------------------------------------------------
-template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
-fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-~MoriRegionGrow( )
-{
-}
-
-// -------------------------------------------------------------------------
-template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
-void
-fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-GenerateInputRequestedRegion( )
-{
- this->Superclass::GenerateInputRequestedRegion( );
- if( this->GetInput( ) )
- {
- TInputImage* input =
- const_cast< TInputImage* >( this->GetInput( ) );
- input->SetRequestedRegionToLargestPossibleRegion( );
-
- } // fi
-}
-
-// -------------------------------------------------------------------------
-template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
-void
-fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-EnlargeOutputRequestedRegion( itk::DataObject* output )
-{
- this->Superclass::EnlargeOutputRequestedRegion( output );
- output->SetRequestedRegionToLargestPossibleRegion( );
-}
-
-// -------------------------------------------------------------------------
-template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
-void fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
-GenerateData( )
-{
- 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( ) );
-}
-
-#endif // __fpa__Image__MoriRegionGrow__hxx__
-
-// eof - $RCSfile$