#include <itkImage.h>
#include <itkImageToImageFilter.h>
#include <itkSimpleFastMutexLock.h>
+#include <fpa/Functors/VertexFunction.h>
namespace fpa
{
typedef typename TImage::RegionType TRegion;
typedef typename TLabels::PixelType TLabel;
+ typedef fpa::Functors::VertexFunction< TIndex, TScalar > TEdgeFunction;
+
protected:
struct _TBoundaryThreadStruct
{
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( );
virtual void GenerateData( ) override;
- virtual TScalar _W( const TIndex& i, const TIndex& j );
virtual TScalar _L( const TIndex& i, const TIndex& j );
// Boundary construction
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;
};
#include <Eigen/Sparse>
// -------------------------------------------------------------------------
-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 );
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( );
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( );
_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 >::
TIndex k = i;
k[ d ] += n;
if( r.IsInside( k ) )
- s += this->_W( i, k );
+ s += this->m_EdgeFunction->Evaluate( i, k );
} // rof
return( s );
}
else
- return( this->_W( i, j ) * TScalar( -1 ) );
+ return( -( this->m_EdgeFunction->Evaluate( i, j ) ) );
}
// -------------------------------------------------------------------------
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
{
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 )
{
this->m_Mutex.Unlock( );
} // fi
-
+
} // fi
} // fi
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 );
{
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;
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* >(
);
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 ) );
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( )
Self& operator=( const Self& other );
protected:
- double m_Alpha;
- double m_Beta;
+ TValue m_Beta;
+ TValue m_Epsilon;
+ bool m_TreatAsWeight;
};
} // ecapseman
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <fpa/Common/OriginalRandomWalker.h>
+#include <fpa/Functors/Dijkstra/Image/Gaussian.h>
// -------------------------------------------------------------------------
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;
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;