#ifndef __fpa__Common__OriginalRandomWalker__hxx__
#define __fpa__Common__OriginalRandomWalker__hxx__
-#include <cmath>
#include <itkImageRegionConstIteratorWithIndex.h>
#include <itkImageRegionIteratorWithIndex.h>
#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( );
- // Build boundary triplets and count labels
- _TTriplets St, Bt;
- std::map< TLabel, unsigned long > labels;
- this->_Boundary( St, labels );
- struct _TTripletsOrd
- {
- bool operator()( const _TTriplet& a, const _TTriplet& b )
- {
- return( a.row( ) < b.row( ) );
- }
- };
- std::sort( St.begin( ), St.end( ), _TTripletsOrd( ) );
- for( unsigned long i = 0; i < St.size( ); ++i )
- Bt.push_back( _TTriplet( i, labels[ St[ i ].col( ) ], St[ i ].value( ) ) );
-
- // Laplacian triplets
- _TTriplets At, Rt;
- this->_Laplacian( At, Rt, St );
-
- // Matrices
- TRegion region = this->GetInput( )->GetRequestedRegion( );
- unsigned long nSeeds = St.size( );
- unsigned long nLabels = labels.size( );
- unsigned long N = region.GetNumberOfPixels( );
-
- std::vector< TLabel > invLabels( nLabels );
- for( typename std::map< TLabel, unsigned long >::value_type s: labels )
- invLabels[ s.second ] = s.first;
-
- _TMatrix B( nSeeds, nLabels );
- B.setFromTriplets( Bt.begin( ), Bt.end( ) );
- B.makeCompressed( );
-
- _TMatrix R( N - nSeeds, nSeeds );
- R.setFromTriplets( Rt.begin( ), Rt.end( ) );
- R.makeCompressed( );
-
- _TMatrix A( N - nSeeds, N - nSeeds );
- A.setFromTriplets( At.begin( ), At.end( ) );
- A.makeCompressed( );
+ // Persisting objects
+ _TMatrix A( 1, 1 ), C( 1, 1 );
+ _TTriplets St;
+ std::vector< TLabel > invLabels;
+
+ { // begin
+ // Build boundary triplets and count labels
+ _TTriplets Bt;
+ std::map< TLabel, unsigned long > labels;
+ itkDebugMacro( << "Building boundary matrix..." );
+ this->_Boundary( St, labels );
+ struct _TTripletsOrd
+ {
+ bool operator()( const _TTriplet& a, const _TTriplet& b )
+ {
+ return( a.row( ) < b.row( ) );
+ }
+ };
+ itkDebugMacro( << "Sorting boundary pixels..." );
+ std::sort( St.begin( ), St.end( ), _TTripletsOrd( ) );
+ itkDebugMacro( << "Assigning boundary pixels..." );
+ for( unsigned long i = 0; i < St.size( ); ++i )
+ Bt.push_back(
+ _TTriplet( i, labels[ St[ i ].col( ) ], St[ i ].value( ) )
+ );
+
+ // Laplacian triplets
+ itkDebugMacro( << "Building laplacian matrix..." );
+ _TTriplets At, Rt;
+ this->_Laplacian( At, Rt, St );
+
+ // Matrices
+ TRegion region = input->GetRequestedRegion( );
+ unsigned long nSeeds = St.size( );
+ unsigned long nLabels = labels.size( );
+ unsigned long N = region.GetNumberOfPixels( );
+
+ itkDebugMacro( << "Creating inverse labels..." );
+ invLabels.resize( nLabels );
+ for( typename std::map< TLabel, unsigned long >::value_type s: labels )
+ invLabels[ s.second ] = s.first;
+
+ itkDebugMacro( << "Creating B matrix..." );
+ _TMatrix B( nSeeds, nLabels );
+ B.setFromTriplets( Bt.begin( ), Bt.end( ) );
+ B.makeCompressed( );
+
+ itkDebugMacro( << "Creating R matrix..." );
+ _TMatrix R( N - nSeeds, nSeeds );
+ R.setFromTriplets( Rt.begin( ), Rt.end( ) );
+ R.makeCompressed( );
+
+ itkDebugMacro( << "Creating C matrix..." );
+ C = R * B;
+
+ itkDebugMacro( << "Creating A matrix..." );
+ A.resize( N - nSeeds, N - nSeeds );
+ A.setFromTriplets( At.begin( ), At.end( ) );
+ A.makeCompressed( );
+ } // end
// Solve dirichlet problem
_TSolver solver;
+ itkDebugMacro( << "Factorizing problem..." );
solver.compute( A );
if( solver.info( ) != Eigen::Success )
- {
- std::cerr << "Error computing." << std::endl;
- } // fi
- _TMatrix x = solver.solve( R * B );
+ itkExceptionMacro( << "Error decomposing matrix." );
+ itkDebugMacro( << "Solving problem..." );
+ _TMatrix x = solver.solve( C );
if( solver.info( ) != Eigen::Success )
- {
- std::cerr << "Error solving." << std::endl;
- } // fi
+ itkExceptionMacro( << "Error solving system." );
// Fill outputs
+ itkDebugMacro( << "Filling output..." );
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 ) ) );
}
// -------------------------------------------------------------------------
thrStr.B = reinterpret_cast< const void* >( &B );
// Configure threader
- const TImage* in = this->GetInputLabels( );
+ const TImage* in = this->GetInput( );
const itk::ImageRegionSplitterBase* split = this->GetImageRegionSplitter( );
const unsigned int nThreads =
split->GetNumberOfSplits(
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;
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / blobdiff