From 1f269742b9c1575ee931b76ebe1c001b7d6596a3 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Tue, 3 Oct 2017 12:08:07 -0500 Subject: [PATCH] ... --- .../algorithms/RandomWalkLabelling.h | 13 +- .../algorithms/RandomWalkLabelling.hxx | 57 +- .../algorithms/RandomWalkSegmentation.h | 51 +- .../algorithms/RandomWalkSegmentation.hxx | 558 ++++++++++-------- lib/fpa/Config.h.in | 49 +- lib/fpa/Filters/Algorithm.h | 4 + lib/fpa/Filters/Algorithm.hxx | 4 +- lib/fpa/Filters/Image/Interface.hxx | 2 +- lib/fpa/Filters/Mori.hxx | 1 + lib/fpa/Filters/RandomWalker.hxx | 3 + lib/fpa/Filters/RegionGrow.hxx | 1 + 11 files changed, 435 insertions(+), 308 deletions(-) diff --git a/appli/CTArteries/algorithms/RandomWalkLabelling.h b/appli/CTArteries/algorithms/RandomWalkLabelling.h index 06667a4..1f3602a 100644 --- a/appli/CTArteries/algorithms/RandomWalkLabelling.h +++ b/appli/CTArteries/algorithms/RandomWalkLabelling.h @@ -9,13 +9,15 @@ /** */ -template< class _TRawImage, class _TLabelsImage > +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > class RandomWalkLabelling : public fpa::Filters::Image::RegionGrow< _TRawImage, _TLabelsImage, typename _TLabelsImage::PixelType > { public: - typedef _TRawImage TRawImage; + typedef _TRawImage TRawImage; + typedef _TCostsImage TCostsImage; typedef _TLabelsImage TLabelsImage; + typedef typename TCostsImage::PixelType TScalar; typedef typename TLabelsImage::PixelType TLabel; typedef fpa::Filters::Image::RegionGrow< TRawImage, TLabelsImage, TLabel > Superclass; @@ -50,7 +52,10 @@ public: itkGetConstMacro( UpperThreshold, double ); itkSetMacro( UpperThreshold, double ); - fpaFilterInputMacro( InputMarks, TLabelsImage ); + itkGetConstMacro( MaxCost, TScalar ); + itkSetMacro( MaxCost, TScalar ); + + fpaFilterInputMacro( InputCosts, TCostsImage ); fpaFilterInputMacro( InputPath, TPath ); public: @@ -84,6 +89,8 @@ protected: double m_LowerThreshold; double m_UpperThreshold; + TScalar m_MaxCost; + unsigned long m_CurrIdx; }; diff --git a/appli/CTArteries/algorithms/RandomWalkLabelling.hxx b/appli/CTArteries/algorithms/RandomWalkLabelling.hxx index 573e88e..3cd2c6a 100644 --- a/appli/CTArteries/algorithms/RandomWalkLabelling.hxx +++ b/appli/CTArteries/algorithms/RandomWalkLabelling.hxx @@ -4,52 +4,54 @@ #ifndef __RandomWalkerLabelling__hxx__ #define __RandomWalkerLabelling__hxx__ +#include + // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -void RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +void RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: SetInputImage( const TRawImage* i ) { this->SetInput( i ); } // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -typename RandomWalkLabelling< _TRawImage, _TLabelsImage >:: -TLabelsImage* RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +typename RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: +TLabelsImage* RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: GetOutputLabels( ) { return( this->GetOutput( ) ); } // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -const typename RandomWalkLabelling< _TRawImage, _TLabelsImage >:: -TLabelsImage* RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +const typename RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: +TLabelsImage* RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: GetOutputLabels( ) const { return( this->GetOutput( ) ); } // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -typename RandomWalkLabelling< _TRawImage, _TLabelsImage >:: -TLabel RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +typename RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: +TLabel RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: GetOutsideLabel( ) const { return( this->GetInitValue( ) ); } // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -void RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +void RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: SetOutsideLabel( const TLabel& v ) { this->SetInitValue( v ); } // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: RandomWalkLabelling( ) : Superclass( ), m_InsideLabel( TLabel( 1 ) ), @@ -57,23 +59,24 @@ RandomWalkLabelling( ) m_UpperLabel( TLabel( 4 ) ), m_Radius( double( 1 ) ), m_LowerThreshold( double( 0 ) ), - m_UpperThreshold( double( 0 ) ) + m_UpperThreshold( double( 0 ) ), + m_MaxCost( std::numeric_limits< TScalar >::max( ) ) { - fpaFilterInputConfigureMacro( InputMarks, TLabelsImage ); + fpaFilterInputConfigureMacro( InputCosts, TCostsImage ); fpaFilterInputConfigureMacro( InputPath, TPath ); this->SetOutsideLabel( TLabel( 2 ) ); } // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: ~RandomWalkLabelling( ) { } // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -void RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +void RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: _PrepareSeeds( const itk::DataObject* reference ) { const TPath* path = this->GetInputPath( ); @@ -90,12 +93,12 @@ _PrepareSeeds( const itk::DataObject* reference ) } // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -void RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +void RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: _PostComputeOutputValue( TNode& n ) { const TRawImage* raw = this->GetInput( ); - const TLabelsImage* labels = this->GetInputMarks( ); + const TCostsImage* costs = this->GetInputCosts( ); const TPath* path = this->GetInputPath( ); TPoint c, p; @@ -105,7 +108,7 @@ _PostComputeOutputValue( TNode& n ) if( d <= this->m_Radius ) { n.FrontId = 1; - if( labels->GetPixel( n.Vertex ) == 0 ) + if( costs->GetPixel( n.Vertex ) == this->m_MaxCost ) { double v = double( raw->GetPixel( n.Vertex ) ); if( v <= this->m_LowerThreshold ) @@ -127,8 +130,8 @@ _PostComputeOutputValue( TNode& n ) } // ------------------------------------------------------------------------- -template< class _TRawImage, class _TLabelsImage > -void RandomWalkLabelling< _TRawImage, _TLabelsImage >:: +template< class _TRawImage, class _TCostsImage, class _TLabelsImage > +void RandomWalkLabelling< _TRawImage, _TCostsImage, _TLabelsImage >:: _Reinitialize( ) { const TPath* path = this->GetInputPath( ); diff --git a/appli/CTArteries/algorithms/RandomWalkSegmentation.h b/appli/CTArteries/algorithms/RandomWalkSegmentation.h index 8d5570d..6644dba 100644 --- a/appli/CTArteries/algorithms/RandomWalkSegmentation.h +++ b/appli/CTArteries/algorithms/RandomWalkSegmentation.h @@ -71,8 +71,55 @@ protected: virtual void GenerateData( ) override; private: - template< class _TImage > - void _save( _TImage* image, const std::string& fname ); + template< class _TIn, class _TSeed > + void _SynchSeed( const _TIn* in, _TSeed& seed ); + + template< class _TIn, class _TOutPtr > + void _Smooth( const _TIn* in, _TOutPtr& out ); + + template< class _TIn, class _TOutPtr, class _TAxisPtr, class _TSeeds > + typename _TIn::RegionType _RawSegmentation( + const _TIn* in, const _TSeeds& seeds, + const typename _TIn::IndexType& s0, + const typename _TIn::IndexType& s1, + const typename _TOutPtr::ObjectType::PixelType& beta, + _TOutPtr& out, _TAxisPtr& out_axis, + typename _TOutPtr::ObjectType::PixelType& oMean, + typename _TOutPtr::ObjectType::PixelType& oSTD + ); + + template< class _TIn, class _TOutPtr > + typename _TIn::RegionType _ROI( + const _TIn* in, const typename _TIn::RegionType& roi, unsigned int pad, + _TOutPtr& out + ); + + template< class _TAxisPtr, class _TRegion > + void _AxisROI( + const typename _TAxisPtr::ObjectType* in, const _TRegion& roi, + _TAxisPtr& out + ); + + template< class _TInRaw, class _TInCosts, class _TAxis, class _TScalar, class _TOutPtr > + void _Label( + const _TInRaw* raw, const _TInCosts* costs, const _TAxis* axis, + const _TScalar& mean, const _TScalar& dev, const _TScalar& radius, + _TOutPtr& out + ); + + template< class _TIn, class _TOutPtr > + void _RandomWalker( + const _TIn* in, const typename _TOutPtr::ObjectType* labels, + const typename _TIn::PixelType& beta, + _TOutPtr& out + ); + + template< class _TIn, class _TSeeds, class _TSeed, class _TOutPtr, class _TAxisPtr > + void _DistanceAndAxis( + const _TIn* in, const _TSeeds& seeds, + const _TSeed& p0, const _TSeed& p1, + _TOutPtr& out_dist, _TAxisPtr& out_axis + ); private: // Purposely not implemented diff --git a/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx b/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx index 66e1018..15fb9af 100644 --- a/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx +++ b/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx @@ -5,20 +5,22 @@ #define __RandomWalkSegmentation__hxx__ #include -#include -#include #include - #include - -#include #include -#include - +#include #include "DijkstraWithMeanAndVariance.h" #include "RandomWalkLabelling.h" -#include +/* TODO + #include + #include + + + #include + + #include +*/ // ------------------------------------------------------------------------- template< class _TInputImage, class _TOutputImage > @@ -125,274 +127,328 @@ template< class _TInputImage, class _TOutputImage > void RandomWalkSegmentation< _TInputImage, _TOutputImage >:: GenerateData( ) { + typedef unsigned char _TLabel; typedef typename TOutputImage::PixelType _TScalar; - typedef DijkstraWithMeanAndVariance< TOutputImage, TOutputImage > _TInit; - typedef typename _TInit::TMarksImage _TMarksImage; - typedef fpa::Functors::Dijkstra::Image::Gaussian< TOutputImage, _TScalar > _TGaussFun; - typedef RandomWalkLabelling< TOutputImage, _TMarksImage > _TLabelling; - typedef typename _TLabelling::TLabelsImage _TLabelsImage; - typedef fpa::Filters::Image::RandomWalker< TOutputImage, _TLabelsImage, _TScalar > _TRandomWalker; - typedef itk::BinaryThresholdImageFilter< _TLabelsImage, _TLabelsImage > _TLabelExtract; - typedef fpa::Filters::Image::ExtractAxis< _TLabelsImage, _TScalar > _TAxisExtract; - typedef ivq::ITK::RegionOfInterestWithPaddingImageFilter< TOutputImage, TOutputImage > _TInputROI; - typedef ivq::ITK::RegionOfInterestWithPaddingImageFilter< _TMarksImage, _TMarksImage > _TMarksROI; - typedef itk::SmoothingRecursiveGaussianImageFilter< TInputImage, TOutputImage > _TSmooth; + typedef itk::Image< _TLabel, TInputImage::ImageDimension > _TLabels; // Prepare initial seeds const TInputImage* input = this->GetInput( ); - if( this->m_StartSeed.IsPoint ) - input->TransformPhysicalPointToIndex( - this->m_StartSeed.Point, - this->m_StartSeed.Index - ); - else - input->TransformIndexToPhysicalPoint( - this->m_StartSeed.Index, - this->m_StartSeed.Point - ); - this->m_StartSeed.IsPoint = true; - if( this->m_EndSeed.IsPoint ) - input->TransformPhysicalPointToIndex( - this->m_EndSeed.Point, - this->m_EndSeed.Index - ); - else - input->TransformIndexToPhysicalPoint( - this->m_EndSeed.Index, - this->m_EndSeed.Point - ); - this->m_EndSeed.IsPoint = true; + this->_SynchSeed( input, this->m_StartSeed ); + this->_SynchSeed( input, this->m_EndSeed ); for( TSeed& seed: this->m_Seeds ) - { - if( seed.IsPoint ) - input->TransformPhysicalPointToIndex( seed.Point, seed.Index ); - else - input->TransformIndexToPhysicalPoint( seed.Index, seed.Point ); - seed.IsPoint = true; - - } // rof - - // Intermediary objects - typename TOutputImage::Pointer smooth_input; - typename TOutputImage::Pointer input_roi; - typename TPath::Pointer init_axis; - typename TPath::Pointer init_axis_roi; - typename TPath::Pointer output_axis_roi; - typename _TMarksImage::Pointer init_marks; - typename _TMarksImage::Pointer init_marks_roi; - typename _TLabelsImage::Pointer init_labels; - typename _TLabelsImage::Pointer final_labels; - typename _TLabelsImage::Pointer inside_labels; - typename TInputImage::RegionType roi; - typename TOutputImage::Pointer output_roi; - double init_mean, init_std; + this->_SynchSeed( input, seed ); // Smooth input - std::cout << "0" << std::endl; - { // begin - typename _TSmooth::Pointer smooth = _TSmooth::New( ); - smooth->SetInput( input ); - smooth->SetNormalizeAcrossScale( true ); - smooth->SetSigmaArray( input->GetSpacing( ) * double( 2 ) ); - smooth->Update( ); - smooth_input = smooth->GetOutput( ); - smooth_input->DisconnectPipeline( ); - - } // end + typename TOutputImage::Pointer smooth_in; + this->_Smooth( this->GetInput( ), smooth_in ); // Initial segmentation - std::cout << "1" << std::endl; - { // begin - typename _TGaussFun::Pointer init_fun = _TGaussFun::New( ); - init_fun->SetBeta( this->m_Beta ); - - typename _TInit::Pointer init = _TInit::New( ); - init->SetInput( smooth_input ); - init->SetWeightFunction( init_fun ); - init->StopAtOneFrontOn( ); - for( TSeed seed: this->m_Seeds ) - init->AddSeed( seed.Point ); - init->Update( ); - - // Get initial values - init_mean = init->GetMean( ); - init_std = init->GetDeviation( ) * this->m_Sigma; - - // Get initial objects - init->GetMinimumSpanningTree( )-> - GetPath( init_axis, this->m_StartSeed.Index, this->m_EndSeed.Index ); - init_marks = init->GetMarks( ); - init_marks->DisconnectPipeline( ); - - typename TInputImage::IndexType min_idx = init->GetMinVertex( ); - typename TInputImage::IndexType max_idx = init->GetMaxVertex( ); - typename TInputImage::SizeType roi_size; - for( unsigned int i = 0; i < TInputImage::ImageDimension; ++i ) - roi_size[ i ] = max_idx[ i ] - min_idx[ i ] + 1; - roi.SetIndex( min_idx ); - roi.SetSize( roi_size ); - - } // end + typename TOutputImage::Pointer init_seg; + typename TPath::Pointer init_axis; + _TScalar init_mean, init_std; + typename TInputImage::RegionType roi = + this->_RawSegmentation( + smooth_in.GetPointer( ), this->m_Seeds, + this->m_StartSeed.Index, this->m_EndSeed.Index, + this->m_Beta, + init_seg, init_axis, init_mean, init_std + ); + init_std *= _TScalar( this->m_Sigma ); // Extract input ROIs - unsigned int pad = 10; - { // begin - typename _TInputROI::Pointer input_roi_filter = _TInputROI::New( ); - input_roi_filter->SetInput( smooth_input ); - input_roi_filter->SetRegionOfInterest( roi ); - input_roi_filter->SetPadding( pad ); - input_roi_filter->Update( ); - input_roi = input_roi_filter->GetOutput( ); - input_roi->DisconnectPipeline( ); - roi = input_roi_filter->GetRegionOfInterest( ); - - } // end - - { // begin - typename _TMarksROI::Pointer init_marks_roi_filter = _TMarksROI::New( ); - init_marks_roi_filter->SetInput( init_marks ); - init_marks_roi_filter->SetRegionOfInterest( roi ); - init_marks_roi_filter->SetPadding( 0 ); - init_marks_roi_filter->Update( ); - init_marks_roi = init_marks_roi_filter->GetOutput( ); - init_marks_roi->DisconnectPipeline( ); - - } // end - - // Convert initial axis - { // begin - init_axis_roi = TPath::New( ); - init_axis_roi->SetReferenceImage( input_roi.GetPointer( ) ); - for( unsigned long i = 0; i < init_axis->GetSize( ); ++i ) - { - TIndex v = init_axis->GetVertex( i ); - for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d ) - v[ d ] -= roi.GetIndex( )[ d ]; - init_axis_roi->AddVertex( v ); - - } // rof - - } // end + typename TOutputImage::Pointer smooth_in_roi, init_seg_roi; + roi = this->_ROI( smooth_in.GetPointer( ), roi, 10, smooth_in_roi ); + this->_ROI( init_seg.GetPointer( ), roi, 0, init_seg_roi ); + typename TPath::Pointer init_axis_roi = TPath::New( ); + init_axis_roi->SetReferenceImage( smooth_in_roi.GetPointer( ) ); + this->_AxisROI( init_axis.GetPointer( ), roi, init_axis_roi ); // Labelling - std::cout << "2" << std::endl; - { // begin - typename _TLabelling::Pointer labelling = _TLabelling::New( ); - labelling->SetInputImage( input_roi ); - labelling->SetInputMarks( init_marks_roi ); - labelling->SetInputPath( init_axis_roi ); - labelling->SetInsideLabel( 1 ); - labelling->SetOutsideLabel( 2 ); - labelling->SetLowerLabel( 3 ); - labelling->SetUpperLabel( 4 ); - labelling->SetRadius( this->m_Radius ); - labelling->SetLowerThreshold( init_mean - init_std ); - labelling->SetUpperThreshold( init_mean + init_std ); - labelling->Update( ); - init_labels = labelling->GetOutputLabels( ); - init_labels->DisconnectPipeline( ); - - } // end + typename _TLabels::Pointer init_labels; + _TScalar radius = _TScalar( this->m_Radius ); + this->_Label( + smooth_in_roi.GetPointer( ), + init_seg_roi.GetPointer( ), + init_axis_roi.GetPointer( ), + init_mean, init_std, radius, + init_labels + ); // Random walker - std::cout << "3" << std::endl; - { // begin - typename _TGaussFun::Pointer rw_fun = _TGaussFun::New( ); - rw_fun->SetBeta( init_std / double( 2 ) ); - - typename _TRandomWalker::Pointer rw = _TRandomWalker::New( ); - rw->SetInputImage( input_roi ); - rw->SetInputLabels( init_labels ); - rw->SetWeightFunction( rw_fun ); - rw->Update( ); - final_labels = rw->GetOutputLabels( ); - final_labels->DisconnectPipeline( ); - - } // end - - // Extract inside label - std::cout << "4" << std::endl; - { // begin - typename _TLabelExtract::Pointer label_extract = _TLabelExtract::New( ); - label_extract->SetInput( final_labels ); - label_extract->SetInsideValue( 1 ); - label_extract->SetOutsideValue( 0 ); - label_extract->SetLowerThreshold( 1 ); - label_extract->SetUpperThreshold( 1 ); - label_extract->Update( ); - inside_labels = label_extract->GetOutput( ); - inside_labels->DisconnectPipeline( ); - - } // end - - // Prepare output values - std::cout << "5" << std::endl; - { // begin - TIndex start_seed, end_seed; - inside_labels->TransformPhysicalPointToIndex( this->m_StartSeed.Point, start_seed ); - inside_labels->TransformPhysicalPointToIndex( this->m_EndSeed.Point, end_seed ); - - typename _TAxisExtract::Pointer axis_extract = _TAxisExtract::New( ); - axis_extract->SetInput( inside_labels ); - axis_extract->AddSeed( start_seed ); - axis_extract->AddSeed( end_seed ); - for( TSeed seed: this->m_Seeds ) - axis_extract->AddSeed( seed.Point ); - axis_extract->SetStartIndex( start_seed ); - axis_extract->SetEndIndex( end_seed ); - axis_extract->Update( ); - output_axis_roi = TPath::New( ); - output_axis_roi->Graft( axis_extract->GetOutput( ) ); - output_roi = axis_extract->GetCenterness( )->GetOutput( ); - output_roi->DisconnectPipeline( ); - - } // end + typename _TLabels::Pointer rw_seg; + this->_RandomWalker( + smooth_in_roi.GetPointer( ), + init_labels.GetPointer( ), + init_std / _TScalar( 2 ), + rw_seg + ); + + // ROI outputs + typename TOutputImage::Pointer out_dist; + typename TPath::Pointer out_axis; + this->_DistanceAndAxis( + rw_seg.GetPointer( ), + this->m_Seeds, this->m_StartSeed, this->m_EndSeed, + out_dist, out_axis + ); // Put everything back to requested region - std::cout << "6" << std::endl; - { // begin - TOutputImage* output = this->GetOutput( ); - output->SetBufferedRegion( output->GetRequestedRegion( ) ); - output->Allocate( ); - output->FillBuffer( -std::numeric_limits< _TScalar >::max( ) ); - - itk::ImageRegionConstIterator< TOutputImage > rIt( - output_roi, output_roi->GetRequestedRegion( ) - ); - itk::ImageRegionIterator< TOutputImage > oIt( output, roi ); - rIt.GoToBegin( ); - oIt.GoToBegin( ); - for( ; !rIt.IsAtEnd( ); ++rIt, ++oIt ) - oIt.Set( rIt.Get( ) ); - - TPath* output_axis = this->GetOutputAxis( ); - output_axis->SetReferenceImage( output ); - for( unsigned long i = 0; i < output_axis_roi->GetSize( ); ++i ) - { - TIndex v = output_axis_roi->GetVertex( i ); - for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d ) - v[ d ] += roi.GetIndex( )[ d ]; - output_axis->AddVertex( v ); - - } // rof - - } // end - std::cout << "7" << std::endl; + /* TODO + std::cout << "6" << std::endl; + { // begin + TOutputImage* output = this->GetOutput( ); + output->SetBufferedRegion( output->GetRequestedRegion( ) ); + output->Allocate( ); + output->FillBuffer( -std::numeric_limits< _TScalar >::max( ) ); + + itk::ImageRegionConstIterator< TOutputImage > rIt( + output_roi, output_roi->GetRequestedRegion( ) + ); + itk::ImageRegionIterator< TOutputImage > oIt( output, roi ); + rIt.GoToBegin( ); + oIt.GoToBegin( ); + for( ; !rIt.IsAtEnd( ); ++rIt, ++oIt ) + oIt.Set( rIt.Get( ) ); + + TPath* output_axis = this->GetOutputAxis( ); + output_axis->SetReferenceImage( output ); + for( unsigned long i = 0; i < output_axis_roi->GetSize( ); ++i ) + { + TIndex v = output_axis_roi->GetVertex( i ); + for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d ) + v[ d ] += roi.GetIndex( )[ d ]; + output_axis->AddVertex( v ); + + } // rof + + } // end + std::cout << "7" << std::endl; + */ } // ------------------------------------------------------------------------- template< class _TInputImage, class _TOutputImage > -template< class _TImage > +template< class _TIn, class _TSeed > void RandomWalkSegmentation< _TInputImage, _TOutputImage >:: -_save( _TImage* image, const std::string& fname ) +_SynchSeed( const _TIn* in, _TSeed& seed ) { - typename itk::ImageFileWriter< _TImage >::Pointer w = - itk::ImageFileWriter< _TImage >::New( ); - w->SetInput( image ); - w->SetFileName( fname ); - w->Update( ); + if( seed.IsPoint ) + in->TransformPhysicalPointToIndex( seed.Point, seed.Index ); + else + in->TransformIndexToPhysicalPoint( seed.Index, seed.Point ); + seed.IsPoint = true; } +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +template< class _TIn, class _TOutPtr > +void RandomWalkSegmentation< _TInputImage, _TOutputImage >:: +_Smooth( const _TIn* in, _TOutPtr& out ) +{ + typedef typename _TOutPtr::ObjectType _TOut; + typedef itk::SmoothingRecursiveGaussianImageFilter< _TIn, _TOut > _TSmooth; + + typename _TSmooth::Pointer smooth = _TSmooth::New( ); + smooth->SetInput( in ); + smooth->SetNormalizeAcrossScale( true ); + smooth->SetSigmaArray( in->GetSpacing( ) * double( 2 ) ); + smooth->Update( ); + out = smooth->GetOutput( ); + out->DisconnectPipeline( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +template< class _TIn, class _TOutPtr, class _TAxisPtr, class _TSeeds > +typename _TIn::RegionType +RandomWalkSegmentation< _TInputImage, _TOutputImage >:: +_RawSegmentation( + const _TIn* in, const _TSeeds& seeds, + const typename _TIn::IndexType& s0, + const typename _TIn::IndexType& s1, + const typename _TOutPtr::ObjectType::PixelType& beta, + _TOutPtr& out, _TAxisPtr& out_axis, + typename _TOutPtr::ObjectType::PixelType& oMean, + typename _TOutPtr::ObjectType::PixelType& oSTD + ) +{ + typedef typename _TOutPtr::ObjectType _TOut; + typedef typename _TOut::PixelType _TScalar; + typedef DijkstraWithMeanAndVariance< _TIn, _TOut > _TInit; + typedef fpa::Functors::Dijkstra::Image::Gaussian< _TIn, _TScalar > _TFun; + + typename _TFun::Pointer fun = _TFun::New( ); + fun->SetBeta( beta ); + + typename _TInit::Pointer init = _TInit::New( ); + init->SetInput( in ); + init->SetWeightFunction( fun ); + init->StopAtOneFrontOn( ); + for( typename _TSeeds::value_type seed: seeds ) + init->AddSeed( seed.Point ); + init->Update( ); + + // Get initial values + oMean = _TScalar( init->GetMean( ) ); + oSTD = _TScalar( init->GetDeviation( ) ); + + // Get initial objects + init->GetMinimumSpanningTree( )->GetPath( out_axis, s0, s1 ); + out = init->GetOutput( ); + out->DisconnectPipeline( ); + + // Get ROI + typename _TIn::IndexType min_idx = init->GetMinVertex( ); + typename _TIn::IndexType max_idx = init->GetMaxVertex( ); + typename _TIn::SizeType roi_size; + for( unsigned int i = 0; i < TInputImage::ImageDimension; ++i ) + roi_size[ i ] = max_idx[ i ] - min_idx[ i ] + 1; + typename _TIn::RegionType roi; + roi.SetIndex( min_idx ); + roi.SetSize( roi_size ); + return( roi ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +template< class _TIn, class _TOutPtr > +typename _TIn::RegionType +RandomWalkSegmentation< _TInputImage, _TOutputImage >:: +_ROI( + const _TIn* in, const typename _TIn::RegionType& roi, unsigned int pad, + _TOutPtr& out + ) +{ + typedef typename _TOutPtr::ObjectType _TOut; + typedef ivq::ITK::RegionOfInterestWithPaddingImageFilter< _TIn, _TOut > _TROI; + + typename _TROI::Pointer filter = _TROI::New( ); + filter->SetInput( in ); + filter->SetRegionOfInterest( roi ); + filter->SetPadding( pad ); + filter->Update( ); + out = filter->GetOutput( ); + out->DisconnectPipeline( ); + return( filter->GetRegionOfInterest( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +template< class _TAxisPtr, class _TRegion > +void RandomWalkSegmentation< _TInputImage, _TOutputImage >:: +_AxisROI( + const typename _TAxisPtr::ObjectType* in, const _TRegion& roi, + _TAxisPtr& out + ) +{ + typedef typename _TAxisPtr::ObjectType _TAxis; + + for( unsigned long i = 0; i < in->GetSize( ); ++i ) + { + typename _TAxis::TIndex v = in->GetVertex( i ); + for( unsigned int d = 0; d < _TAxis::Dimension; ++d ) + v[ d ] -= roi.GetIndex( )[ d ]; + out->AddVertex( v ); + + } // rof +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +template< class _TInRaw, class _TInCosts, class _TAxis, class _TScalar, class _TOutPtr > +void RandomWalkSegmentation< _TInputImage, _TOutputImage >:: +_Label( + const _TInRaw* raw, const _TInCosts* costs, const _TAxis* axis, + const _TScalar& mean, const _TScalar& dev, const _TScalar& radius, + _TOutPtr& out + ) +{ + typedef typename _TOutPtr::ObjectType _TOut; + typedef RandomWalkLabelling< _TInRaw, _TInCosts, _TOut > _TLabel; + + typename _TLabel::Pointer label = _TLabel::New( ); + label->SetInputImage( raw ); + label->SetInputCosts( costs ); + label->SetInputPath( axis ); + label->SetInsideLabel( 1 ); + label->SetOutsideLabel( 2 ); + label->SetLowerLabel( 3 ); + label->SetUpperLabel( 4 ); + label->SetRadius( radius ); + label->SetLowerThreshold( mean - dev ); + label->SetUpperThreshold( mean + dev ); + label->Update( ); + out = label->GetOutputLabels( ); + out->DisconnectPipeline( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +template< class _TIn, class _TOutPtr > +void RandomWalkSegmentation< _TInputImage, _TOutputImage >:: +_RandomWalker( + const _TIn* in, const typename _TOutPtr::ObjectType* labels, + const typename _TIn::PixelType& beta, + _TOutPtr& out + ) +{ + typedef typename _TIn::PixelType _TScalar; + typedef typename _TOutPtr::ObjectType _TOut; + typedef fpa::Functors::Dijkstra::Image::Gaussian< _TIn, _TScalar > _TFun; + typedef fpa::Filters::Image::RandomWalker< _TIn, _TOut, _TScalar > _TRandomWalker; + typedef itk::BinaryThresholdImageFilter< _TOut, _TOut > _TExtract; + + typename _TFun::Pointer fun = _TFun::New( ); + fun->SetBeta( beta ); + + typename _TRandomWalker::Pointer rw = _TRandomWalker::New( ); + rw->SetInputImage( in ); + rw->SetInputLabels( labels ); + rw->SetWeightFunction( fun ); + + typename _TExtract::Pointer extract = _TExtract::New( ); + extract->SetInput( rw->GetOutputLabels( ) ); + extract->SetInsideValue( 1 ); + extract->SetOutsideValue( 0 ); + extract->SetLowerThreshold( 1 ); + extract->SetUpperThreshold( 1 ); + extract->Update( ); + out = extract->GetOutput( ); + out->DisconnectPipeline( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +template< class _TIn, class _TSeeds, class _TSeed, class _TOutPtr, class _TAxisPtr > +void RandomWalkSegmentation< _TInputImage, _TOutputImage >:: +_DistanceAndAxis( + const _TIn* in, const _TSeeds& seeds, + const _TSeed& p0, const _TSeed& p1, + _TOutPtr& out_dist, _TAxisPtr& out_axis + ) +{ + typedef typename _TOutPtr::ObjectType _TOut; + typedef typename _TOut::PixelType _TScalar; + typedef fpa::Filters::Image::ExtractAxis< _TIn, _TScalar > _TExtract; + + // Prepare output values + typename _TIn::IndexType s0, s1; + in->TransformPhysicalPointToIndex( p0.Point, s0 ); + in->TransformPhysicalPointToIndex( p1.Point, s1 ); + + typename _TExtract::Pointer extract = _TExtract::New( ); + extract->SetInput( const_cast< _TIn* >( in ) ); + extract->AddSeed( s0 ); + extract->AddSeed( s1 ); + for( typename _TSeeds::value_type seed: seeds ) + extract->AddSeed( seed.Point ); + extract->SetStartIndex( s0 ); + extract->SetEndIndex( s1 ); + extract->Update( ); + out_axis = TPath::New( ); + out_axis->Graft( extract->GetOutput( ) ); + this->_Smooth( extract->GetCenterness( )->GetOutput( ), out_dist ); +} #endif // __RandomWalkSegmentation__hxx__ diff --git a/lib/fpa/Config.h.in b/lib/fpa/Config.h.in index d00ed7f..f93fb23 100644 --- a/lib/fpa/Config.h.in +++ b/lib/fpa/Config.h.in @@ -25,29 +25,32 @@ typedef __t__::TInternalTraits::TVertex TVertex // ------------------------------------------------------------------------- -#define fpaFilterInputMacro( __n__, __t__ ) \ - private: \ - unsigned int m_##__n__##Idx; \ - public: \ - __t__* Get##__n__( ) \ - { \ - return( \ - dynamic_cast< __t__* >( \ - this->itk::ProcessObject::GetInput( \ - this->m_##__n__##Idx \ - ) ) ); \ - } \ - const __t__* Get##__n__( ) const \ - { \ - return( \ - dynamic_cast< const __t__* >( \ - this->itk::ProcessObject::GetInput( \ - this->m_##__n__##Idx \ - ) ) ); \ - } \ - void Set##__n__( __t__* i ) \ - { \ - this->itk::ProcessObject::SetNthInput( this->m_##__n__##Idx, i ); \ +#define fpaFilterInputMacro( __n__, __t__ ) \ + private: \ + unsigned int m_##__n__##Idx; \ + public: \ + __t__* Get##__n__( ) \ + { \ + return( \ + dynamic_cast< __t__* >( \ + this->itk::ProcessObject::GetInput( \ + this->m_##__n__##Idx \ + ) ) ); \ + } \ + const __t__* Get##__n__( ) const \ + { \ + return( \ + dynamic_cast< const __t__* >( \ + this->itk::ProcessObject::GetInput( \ + this->m_##__n__##Idx \ + ) ) ); \ + } \ + void Set##__n__( const __t__* i ) \ + { \ + this->itk::ProcessObject::SetNthInput( \ + this->m_##__n__##Idx, \ + const_cast< __t__* >( i ) \ + ); \ } // ------------------------------------------------------------------------- diff --git a/lib/fpa/Filters/Algorithm.h b/lib/fpa/Filters/Algorithm.h index bf0a44f..0d9eb90 100644 --- a/lib/fpa/Filters/Algorithm.h +++ b/lib/fpa/Filters/Algorithm.h @@ -43,6 +43,9 @@ namespace fpa itkGetConstMacro( InitValue, TOutputValue ); itkSetMacro( InitValue, TOutputValue ); + itkGetConstMacro( FillValue, TOutputValue ); + itkSetMacro( FillValue, TOutputValue ); + public: virtual itk::ModifiedTimeType GetMTime( ) const override; virtual void InvokeEvent( const itk::EventObject& e ); @@ -100,6 +103,7 @@ namespace fpa protected: bool m_VisualDebug; TOutputValue m_InitValue; + TOutputValue m_FillValue; std::set< itk::Object* > m_AssociatedObjects; }; diff --git a/lib/fpa/Filters/Algorithm.hxx b/lib/fpa/Filters/Algorithm.hxx index be4eacc..b591581 100644 --- a/lib/fpa/Filters/Algorithm.hxx +++ b/lib/fpa/Filters/Algorithm.hxx @@ -5,6 +5,7 @@ #ifndef __fpa__Filters__Algorithm__hxx__ #define __fpa__Filters__Algorithm__hxx__ +#include #include // ------------------------------------------------------------------------- @@ -65,7 +66,8 @@ Algorithm( ) TTraits::TMarksInterface( this ), TTraits::TSeedsInterface( this ), m_VisualDebug( false ), - m_InitValue( TOutputValue( 0 ) ) + m_InitValue( TOutputValue( 0 ) ), + m_FillValue( std::numeric_limits< TOutputValue >::max( ) ) { } diff --git a/lib/fpa/Filters/Image/Interface.hxx b/lib/fpa/Filters/Image/Interface.hxx index 2ea6e03..aaafad2 100644 --- a/lib/fpa/Filters/Image/Interface.hxx +++ b/lib/fpa/Filters/Image/Interface.hxx @@ -45,7 +45,7 @@ _ConfigureOutputs( ) out->SetOrigin( in->GetOrigin( ) ); out->SetDirection( in->GetDirection( ) ); out->Allocate( ); - out->FillBuffer( this->GetInitValue( ) ); + out->FillBuffer( this->GetFillValue( ) ); TMarksImage* marks = this->GetMarks( ); marks->SetLargestPossibleRegion( in->GetLargestPossibleRegion( ) ); diff --git a/lib/fpa/Filters/Mori.hxx b/lib/fpa/Filters/Mori.hxx index 68048f1..e888172 100644 --- a/lib/fpa/Filters/Mori.hxx +++ b/lib/fpa/Filters/Mori.hxx @@ -157,6 +157,7 @@ Mori( ) : Superclass( true ) { this->SetOutsideValue( TOutputValue( 0 ) ); + this->SetFillValue( TOutputValue( 0 ) ); this->m_InsideValue = std::numeric_limits< TOutputValue >::max( ); this->m_Predicate = TPredicate::New( ); this->m_Predicate->StrictOff( ); diff --git a/lib/fpa/Filters/RandomWalker.hxx b/lib/fpa/Filters/RandomWalker.hxx index 18b8d2b..f5c18c7 100644 --- a/lib/fpa/Filters/RandomWalker.hxx +++ b/lib/fpa/Filters/RandomWalker.hxx @@ -5,6 +5,8 @@ #ifndef __fpa__Filters__RandomWalker__hxx__ #define __fpa__Filters__RandomWalker__hxx__ +#include + // ------------------------------------------------------------------------- template< class _TDataInterface > void fpa::Filters::RandomWalker< _TDataInterface >:: @@ -42,6 +44,7 @@ RandomWalker( ) : Superclass( false ) { this->SetInitValue( TCost( 0 ) ); + this->SetFillValue( std::numeric_limits< TCost >::max( ) ); /* TODO this->SetWeight( fpa::Functors::RegionGrow::Tautology< TInputValue >::New( ) diff --git a/lib/fpa/Filters/RegionGrow.hxx b/lib/fpa/Filters/RegionGrow.hxx index 50ae38d..8028376 100644 --- a/lib/fpa/Filters/RegionGrow.hxx +++ b/lib/fpa/Filters/RegionGrow.hxx @@ -62,6 +62,7 @@ RegionGrow( ) : Superclass( false ) { this->SetOutsideValue( TOutputValue( 0 ) ); + this->SetFillValue( TOutputValue( 0 ) ); this->m_InsideValue = std::numeric_limits< TOutputValue >::max( ); this->SetPredicate( fpa::Functors::RegionGrow::Tautology< TInputValue >::New( ) -- 2.47.1