--- /dev/null
+#include <itkCommand.h>
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+
+#include <fpa/Image/MoriRegionGrow.h>
+
+// -------------------------------------------------------------------------
+static const unsigned int VDim = 3;
+typedef short TPixel;
+typedef itk::Image< TPixel, VDim > TImage;
+typedef itk::ImageFileReader< TImage > TReader;
+typedef itk::ImageFileWriter< TImage > TWriter;
+typedef fpa::Image::MoriRegionGrow< TImage, TImage > TFilter;
+typedef itk::ImageFileWriter< TFilter::TAuxiliaryImage > TAuxWriter;
+
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+ // Get arguments
+ if( argc < 7 + VDim )
+ {
+ std::cerr
+ << "Usage: " << argv[ 0 ]
+ << " input_image output_image auxiliary_image lower upper delta";
+ for( unsigned int i = 0; i < VDim; ++i )
+ std::cerr << " s_" << i;
+ std::cerr << std::endl;
+ return( 1 );
+
+ } // fi
+ std::string input_image_filename = argv[ 1 ];
+ std::string output_image_filename = argv[ 2 ];
+ std::string auxiliary_image_filename = argv[ 3 ];
+ TPixel lower = std::atof( argv[ 4 ] );
+ TPixel upper = std::atof( argv[ 5 ] );
+ TPixel delta = std::atof( argv[ 6 ] );
+
+ TReader::Pointer reader = TReader::New( );
+ reader->SetFileName( input_image_filename );
+
+ TFilter::Pointer filter = TFilter::New( );
+ filter->SetInput( reader->GetOutput( ) );
+ filter->SetThresholdRange( lower, upper, delta );
+ TImage::IndexType seed;
+ for( int i = 0; i < VDim; ++i )
+ seed[ i ] = std::atoi( argv[ i + 7 ] );
+ filter->SetSeed( seed );
+ filter->SetInsideValue( 255 );
+ filter->SetOutsideValue( 0 );
+
+ TWriter::Pointer writer = TWriter::New( );
+ writer->SetInput( filter->GetOutput( ) );
+ writer->SetFileName( output_image_filename );
+ try
+ {
+ writer->Update( );
+ }
+ catch( std::exception& err )
+ {
+ std::cerr << "ERROR: " << err.what( ) << std::endl;
+ return( 1 );
+
+ } // yrt
+
+ TAuxWriter::Pointer aux_writer = TAuxWriter::New( );
+ aux_writer->SetInput( filter->GetAuxiliaryImage( ) );
+ aux_writer->SetFileName( auxiliary_image_filename );
+ try
+ {
+ aux_writer->Update( );
+ }
+ catch( std::exception& err )
+ {
+ std::cerr << "ERROR: " << err.what( ) << std::endl;
+ return( 1 );
+
+ } // yrt
+ return( 0 );
+}
+
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Image__MoriRegionGrow__h__
+#define __fpa__Image__MoriRegionGrow__h__
+
+#include <itkImageToImageFilter.h>
+#include <itkBinaryThresholdImageFilter.h>
+
+namespace fpa
+{
+ namespace Image
+ {
+ /**
+ */
+ template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel = unsigned short >
+ class MoriRegionGrow
+ : public itk::ImageToImageFilter< _TInputImage, _TOutputImage >
+ {
+ public:
+ typedef _TInputImage TInputImage;
+ typedef _TOutputImage TOutputImage;
+ typedef _TAuxiliaryPixel TAuxiliaryPixel;
+ typedef itk::Image< TAuxiliaryPixel, TInputImage::ImageDimension > TAuxiliaryImage;
+
+ typedef MoriRegionGrow Self;
+ typedef itk::ImageToImageFilter< TInputImage, TOutputImage > Superclass;
+ typedef itk::SmartPointer< Self > Pointer;
+ typedef itk::SmartPointer< const Self > ConstPointer;
+
+ typedef typename TInputImage::IndexType TIndex;
+ typedef typename TInputImage::RegionType TRegion;
+ typedef typename TInputImage::PixelType TInputPixel;
+ typedef typename TOutputImage::PixelType TOutputPixel;
+
+ struct TCurveData
+ {
+ TInputPixel XValue;
+ unsigned long YValue;
+ double Diff1;
+ TCurveData( TInputPixel v, unsigned long c )
+ {
+ this->XValue = v;
+ this->YValue = c;
+ this->Diff1 = double( 0 );
+ }
+ };
+ typedef std::vector< TCurveData > TCurve;
+ typedef itk::BinaryThresholdImageFilter< TAuxiliaryImage, TOutputImage > TThresholdFilter;
+
+ public:
+ itkNewMacro( Self );
+ itkTypeMacro( fpa::Image::MoriRegionGrow, itk::ImageToImageFilter );
+
+ itkGetConstMacro( Seed, TIndex );
+ itkGetConstMacro( InsideValue, TOutputPixel );
+ itkGetConstMacro( OutsideValue, TOutputPixel );
+ itkGetConstMacro( LowerThreshold, TInputPixel );
+ itkGetConstMacro( UpperThreshold, TInputPixel );
+ itkGetConstMacro( DeltaThreshold, TInputPixel );
+ itkGetConstMacro( OptimumThreshold, TInputPixel );
+ itkGetConstMacro( Curve, TCurve );
+
+ itkSetMacro( Seed, TIndex );
+ itkSetMacro( InsideValue, TOutputPixel );
+ itkSetMacro( OutsideValue, TOutputPixel );
+ itkSetMacro( LowerThreshold, TInputPixel );
+ itkSetMacro( UpperThreshold, TInputPixel );
+ itkSetMacro( DeltaThreshold, TInputPixel );
+
+ public:
+ TAuxiliaryImage* GetAuxiliaryImage( );
+ const TAuxiliaryImage* GetAuxiliaryImage( ) const;
+
+ void SetThresholdRange(
+ const TInputPixel& lower, const TInputPixel& upper,
+ const TInputPixel& delta = TInputPixel( 1 )
+ );
+
+ protected:
+ MoriRegionGrow( );
+ virtual ~MoriRegionGrow( );
+
+ virtual void GenerateInputRequestedRegion( ) override;
+ virtual void EnlargeOutputRequestedRegion(
+ itk::DataObject* output
+ ) override;
+ virtual void GenerateData( ) override;
+
+ private:
+ // Purposely not implemented
+ MoriRegionGrow( const Self& other );
+ Self& operator=( const Self& other );
+
+ protected:
+ TIndex m_Seed;
+ TOutputPixel m_InsideValue;
+ TOutputPixel m_OutsideValue;
+ TInputPixel m_LowerThreshold;
+ TInputPixel m_UpperThreshold;
+ TInputPixel m_DeltaThreshold;
+ TInputPixel m_OptimumThreshold;
+
+ TCurve m_Curve;
+ typename TThresholdFilter::Pointer m_ThresholdFilter;
+ };
+
+ } // ecapseman
+
+} // ecapseman
+
+#ifndef ITK_MANUAL_INSTANTIATION
+# include <fpa/Image/MoriRegionGrow.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __fpa__Image__MoriRegionGrow__h__
+
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @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$
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff