From 4b880a6cf8104f7176a1701c01151bd861ba78ba Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Tue, 21 Mar 2017 15:37:35 -0500 Subject: [PATCH] ... --- examples/CMakeLists.txt | 1 + examples/RegionGrow_Mori.cxx | 82 ++++++++++ libs/fpa/Image/MoriRegionGrow.h | 120 ++++++++++++++ libs/fpa/Image/MoriRegionGrow.hxx | 254 ++++++++++++++++++++++++++++++ 4 files changed, 457 insertions(+) create mode 100644 examples/RegionGrow_Mori.cxx create mode 100644 libs/fpa/Image/MoriRegionGrow.h create mode 100644 libs/fpa/Image/MoriRegionGrow.hxx diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index cbbfcf9..f0c8826 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -4,6 +4,7 @@ IF(BUILD_EXAMPLES) _examples RegionGrow_Tautology RegionGrow_BinaryThreshold + RegionGrow_Mori #CreateMoriInputImage #BronchiiInitialSegmentationWithMori #BronchiiInitialSegmentationWithBinaryThresholdRegionGrow diff --git a/examples/RegionGrow_Mori.cxx b/examples/RegionGrow_Mori.cxx new file mode 100644 index 0000000..7fd51b4 --- /dev/null +++ b/examples/RegionGrow_Mori.cxx @@ -0,0 +1,82 @@ +#include +#include +#include +#include + +#include + +// ------------------------------------------------------------------------- +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$ diff --git a/libs/fpa/Image/MoriRegionGrow.h b/libs/fpa/Image/MoriRegionGrow.h new file mode 100644 index 0000000..606d4f7 --- /dev/null +++ b/libs/fpa/Image/MoriRegionGrow.h @@ -0,0 +1,120 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Image__MoriRegionGrow__h__ +#define __fpa__Image__MoriRegionGrow__h__ + +#include +#include + +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 +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__MoriRegionGrow__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/MoriRegionGrow.hxx b/libs/fpa/Image/MoriRegionGrow.hxx new file mode 100644 index 0000000..5efdf65 --- /dev/null +++ b/libs/fpa/Image/MoriRegionGrow.hxx @@ -0,0 +1,254 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Image__MoriRegionGrow__hxx__ +#define __fpa__Image__MoriRegionGrow__hxx__ + +#include +#include + +// ------------------------------------------------------------------------- +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$ -- 2.45.0