From 1f5281211c3972084b3a90f6f48b7aae454ecae7 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Thu, 21 Sep 2017 11:15:33 +0200 Subject: [PATCH] ... --- lib/fpa/Common/OriginalRandomWalker.h | 32 +++++----- lib/fpa/Common/OriginalRandomWalker.hxx | 71 +++++++++------------- lib/fpa/Functors/Dijkstra/Image/Gaussian.h | 38 ++++++++---- tests/image/RandomWalker/Original.cxx | 15 +++-- 4 files changed, 80 insertions(+), 76 deletions(-) diff --git a/lib/fpa/Common/OriginalRandomWalker.h b/lib/fpa/Common/OriginalRandomWalker.h index d4af732..0c3a9f8 100644 --- a/lib/fpa/Common/OriginalRandomWalker.h +++ b/lib/fpa/Common/OriginalRandomWalker.h @@ -11,6 +11,7 @@ #include #include #include +#include namespace fpa { @@ -38,6 +39,8 @@ namespace fpa typedef typename TImage::RegionType TRegion; typedef typename TLabels::PixelType TLabel; + typedef fpa::Functors::VertexFunction< TIndex, TScalar > TEdgeFunction; + protected: struct _TBoundaryThreadStruct { @@ -68,21 +71,17 @@ namespace fpa fpa::Common::OriginalRandomWalker, itk::ImageToImageFilter ); - itkGetConstMacro( Beta, TScalar ); - itkSetMacro( Beta, TScalar ); - - itkGetConstMacro( Epsilon, TScalar ); - itkSetMacro( Epsilon, TScalar ); - - itkBooleanMacro( NormalizeWeights ); - itkGetConstMacro( NormalizeWeights, bool ); - itkSetMacro( NormalizeWeights, bool ); + itkGetConstObjectMacro( EdgeFunction, TEdgeFunction ); + itkGetObjectMacro( EdgeFunction, TEdgeFunction ); + itkSetObjectMacro( EdgeFunction, TEdgeFunction ); fpaFilterInputMacro( InputLabels, TLabels ); fpaFilterOutputMacro( OutputProbabilities, TScalarImage ); - public: - void AddSeed( const TIndex& seed, const TLabel& label ); + /* TODO + public: + void AddSeed( const TIndex& seed, const TLabel& label ); + */ protected: OriginalRandomWalker( ); @@ -90,7 +89,6 @@ namespace fpa virtual void GenerateData( ) override; - virtual TScalar _W( const TIndex& i, const TIndex& j ); virtual TScalar _L( const TIndex& i, const TIndex& j ); // Boundary construction @@ -161,12 +159,12 @@ namespace fpa Self& operator=( const Self& other ); protected: - std::vector< TIndex > m_Seeds; - std::vector< TLabel > m_Labels; + /* TODO + std::vector< TIndex > m_Seeds; + std::vector< TLabel > m_Labels; + */ - TScalar m_Beta; - TScalar m_Epsilon; - bool m_NormalizeWeights; + typename TEdgeFunction::Pointer m_EdgeFunction; itk::SimpleFastMutexLock m_Mutex; }; diff --git a/lib/fpa/Common/OriginalRandomWalker.hxx b/lib/fpa/Common/OriginalRandomWalker.hxx index 9bc08d4..4b60c1b 100644 --- a/lib/fpa/Common/OriginalRandomWalker.hxx +++ b/lib/fpa/Common/OriginalRandomWalker.hxx @@ -11,23 +11,22 @@ #include // ------------------------------------------------------------------------- -template< class _TImage, class _TLabels, class _TScalar > -void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: -AddSeed( const TIndex& seed, const TLabel& label ) -{ - this->m_Seeds.push_back( seed ); - this->m_Labels.push_back( label ); - this->Modified( ); -} +/* TODO + template< class _TImage, class _TLabels, class _TScalar > + void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: + AddSeed( const TIndex& seed, const TLabel& label ) + { + this->m_Seeds.push_back( seed ); + this->m_Labels.push_back( label ); + this->Modified( ); + } +*/ // ------------------------------------------------------------------------- template< class _TImage, class _TLabels, class _TScalar > fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: OriginalRandomWalker( ) - : Superclass( ), - m_Beta( TScalar( 90 ) ), - m_Epsilon( TScalar( 1e-5 ) ), - m_NormalizeWeights( true ) + : Superclass( ) { fpaFilterInputConfigureMacro( InputLabels, TLabels ); fpaFilterOutputConfigureMacro( OutputProbabilities, TScalarImage ); @@ -51,6 +50,12 @@ GenerateData( ) typedef Eigen::SparseMatrix< TScalar > _TMatrix; typedef Eigen::SimplicialLDLT< _TMatrix > _TSolver; + // Configure edge function + if( this->m_EdgeFunction.IsNull( ) ) + itkExceptionMacro( << "Undefined edge function." ); + const TImage* input = this->GetInput( ); + this->m_EdgeFunction->SetDataObject( input ); + // Allocate outputs this->AllocateOutputs( ); @@ -74,7 +79,7 @@ GenerateData( ) this->_Laplacian( At, Rt, St ); // Matrices - TRegion region = this->GetInput( )->GetRequestedRegion( ); + TRegion region = input->GetRequestedRegion( ); unsigned long nSeeds = St.size( ); unsigned long nLabels = labels.size( ); unsigned long N = region.GetNumberOfPixels( ); @@ -99,34 +104,15 @@ GenerateData( ) _TSolver solver; solver.compute( A ); if( solver.info( ) != Eigen::Success ) - { - std::cerr << "Error computing." << std::endl; - } // fi + itkExceptionMacro( << "Error decomposing matrix." ); _TMatrix x = solver.solve( R * B ); if( solver.info( ) != Eigen::Success ) - { - std::cerr << "Error solving." << std::endl; - } // fi + itkExceptionMacro( << "Error solving system." ); // Fill outputs this->_Output( x, St, invLabels ); } -// ------------------------------------------------------------------------- -template< class _TImage, class _TLabels, class _TScalar > -_TScalar fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: -_W( const TIndex& i, const TIndex& j ) -{ - const TImage* in = this->GetInput( ); - TScalar a = TScalar( in->GetPixel( i ) ); - TScalar b = TScalar( in->GetPixel( j ) ); - TScalar v = std::exp( -this->m_Beta * std::fabs( a - b ) ); - if( v < this->m_Epsilon ) - return( this->m_Epsilon ); - else - return( v ); -} - // ------------------------------------------------------------------------- template< class _TImage, class _TLabels, class _TScalar > _TScalar fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: @@ -143,7 +129,7 @@ _L( const TIndex& i, const TIndex& j ) TIndex k = i; k[ d ] += n; if( r.IsInside( k ) ) - s += this->_W( i, k ); + s += this->m_EdgeFunction->Evaluate( i, k ); } // rof @@ -151,7 +137,7 @@ _L( const TIndex& i, const TIndex& j ) return( s ); } else - return( this->_W( i, j ) * TScalar( -1 ) ); + return( -( this->m_EdgeFunction->Evaluate( i, j ) ) ); } // ------------------------------------------------------------------------- @@ -317,13 +303,13 @@ _ThreadedLaplacian( const TImage* in = this->GetInput( ); const TLabels* in_labels = this->GetInputLabels( ); - TRegion rqRegion = in->GetRequestedRegion( ); + TRegion reqRegion = in->GetRequestedRegion( ); _TIt it( in, region ); for( it.GoToBegin( ); !it.IsAtEnd( ); ++it ) { TIndex idx = it.GetIndex( ); bool iSeed = ( in_labels->GetPixel( idx ) != 0 ); - unsigned long i = Self::_1D( idx, rqRegion ); + unsigned long i = Self::_1D( idx, reqRegion ); unsigned long si; // A's diagonal values @@ -344,10 +330,10 @@ _ThreadedLaplacian( { TIndex jdx = idx; jdx[ d ] += s; - if( rqRegion.IsInside( jdx ) ) + if( reqRegion.IsInside( jdx ) ) { TScalar L = this->_L( idx, jdx ); - unsigned long j = Self::_1D( jdx, rqRegion ); + unsigned long j = Self::_1D( jdx, reqRegion ); bool jSeed = ( in_labels->GetPixel( jdx ) != 0 ); if( !jSeed ) { @@ -365,7 +351,7 @@ _ThreadedLaplacian( this->m_Mutex.Unlock( ); } // fi - + } // fi } // fi @@ -451,6 +437,7 @@ _ThreadedOutput( const TLabels* in_labels = this->GetInputLabels( ); TLabels* out_labels = this->GetOutput( ); TScalarImage* out_probs = this->GetOutputProbabilities( ); + TRegion reqRegion = out_labels->GetRequestedRegion( ); itk::ImageRegionConstIteratorWithIndex< TLabels > iIt( in_labels, region ); itk::ImageRegionIteratorWithIndex< TLabels > oIt( out_labels, region ); itk::ImageRegionIteratorWithIndex< TScalarImage > pIt( out_probs, region ); @@ -461,7 +448,7 @@ _ThreadedOutput( { if( iIt.Get( ) == 0 ) { - unsigned long i = Self::_1D( iIt.GetIndex( ), region ); + unsigned long i = Self::_1D( iIt.GetIndex( ), reqRegion ); unsigned long j = Self::_NearSeedIndex( i, *S ); TScalar maxP = X->coeff( j, 0 ); unsigned long maxL = 0; diff --git a/lib/fpa/Functors/Dijkstra/Image/Gaussian.h b/lib/fpa/Functors/Dijkstra/Image/Gaussian.h index e917e95..2effb51 100644 --- a/lib/fpa/Functors/Dijkstra/Image/Gaussian.h +++ b/lib/fpa/Functors/Dijkstra/Image/Gaussian.h @@ -38,14 +38,20 @@ namespace fpa fpa::Functors::VertexFunction ); - itkGetConstMacro( Alpha, double ); - itkSetMacro( Alpha, double ); + itkGetConstMacro( Beta, TValue ); + itkSetMacro( Beta, TValue ); - itkGetConstMacro( Beta, double ); - itkSetMacro( Beta, double ); + itkGetConstMacro( Epsilon, TValue ); + itkSetMacro( Epsilon, TValue ); + + itkBooleanMacro( TreatAsWeight ); + itkGetConstMacro( TreatAsWeight, bool ); + itkSetMacro( TreatAsWeight, bool ); public: - virtual TValue Evaluate( const TVertex& v, const TVertex& p ) const override + virtual TValue Evaluate( + const TVertex& v, const TVertex& p + ) const override { const TImage* image = dynamic_cast< const TImage* >( @@ -53,11 +59,15 @@ namespace fpa ); if( image != NULL ) { - double d = double( image->GetPixel( v ) ); - d -= double( image->GetPixel( p ) ); + TValue d = TValue( image->GetPixel( v ) ); + d -= TValue( image->GetPixel( p ) ); d /= this->m_Beta; - d = std::exp( d * d ) - double( 1 ); - return( TValue( std::pow( d, this->m_Alpha ) ) ); + d *= d; + if( this->m_TreatAsWeight ) d = std::exp( d ); + else d = std::exp( -d ); + + if( d < this->m_Epsilon ) return( this->m_Epsilon ); + else return( d ); } else return( TValue( -1 ) ); @@ -66,8 +76,9 @@ namespace fpa protected: Gaussian( ) : Superclass( ), - m_Alpha( double( 1 ) ), - m_Beta( double( 1 ) ) + m_Beta( TValue( 1 ) ), + m_Epsilon( TValue( 1e-5 ) ), + m_TreatAsWeight( true ) { } virtual ~Gaussian( ) @@ -80,8 +91,9 @@ namespace fpa Self& operator=( const Self& other ); protected: - double m_Alpha; - double m_Beta; + TValue m_Beta; + TValue m_Epsilon; + bool m_TreatAsWeight; }; } // ecapseman diff --git a/tests/image/RandomWalker/Original.cxx b/tests/image/RandomWalker/Original.cxx index 46ba198..0ea2bd0 100644 --- a/tests/image/RandomWalker/Original.cxx +++ b/tests/image/RandomWalker/Original.cxx @@ -7,11 +7,13 @@ #include #include #include +#include // ------------------------------------------------------------------------- const unsigned int Dim = 2; typedef unsigned char TPixel; typedef unsigned char TLabel; +typedef double TScalar; typedef itk::Image< TPixel, Dim > TImage; typedef itk::Image< TLabel, Dim > TLabels; @@ -52,14 +54,19 @@ int main( int argc, char* argv[] ) TLabelsReader::Pointer input_labels_reader = TLabelsReader::New( ); input_labels_reader->SetFileName( input_labels_filename ); + // Random walker function + typedef fpa::Functors::Dijkstra::Image::Gaussian< TImage, TScalar > TFunction; + TFunction::Pointer edge = TFunction::New( ); + edge->SetBeta( beta ); + edge->SetEpsilon( eps ); + edge->TreatAsWeightOff( ); + // Random walker - typedef fpa::Common::OriginalRandomWalker< TImage, TLabels > TFilter; + typedef fpa::Common::OriginalRandomWalker< TImage, TLabels, TScalar > TFilter; TFilter::Pointer rw = TFilter::New( ); rw->SetInput( input_image_reader->GetOutput( ) ); rw->SetInputLabels( input_labels_reader->GetOutput( ) ); - rw->SetBeta( beta ); - rw->SetEpsilon( eps ); - rw->SetNormalizeWeights( normalize ); + rw->SetEdgeFunction( edge ); // Save labels typedef itk::ImageFileWriter< TFilter::TLabels > TLabelsWriter; -- 2.47.1