#ifndef __fpa__Common__OriginalRandomWalker__hxx__
#define __fpa__Common__OriginalRandomWalker__hxx__
-#include <cmath>
#include <itkImageRegionConstIteratorWithIndex.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <Eigen/Sparse>
// 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 = input->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 )
itkExceptionMacro( << "Error decomposing matrix." );
- _TMatrix x = solver.solve( R * B );
+ itkDebugMacro( << "Solving problem..." );
+ _TMatrix x = solver.solve( C );
if( solver.info( ) != Eigen::Success )
itkExceptionMacro( << "Error solving system." );
// Fill outputs
+ itkDebugMacro( << "Filling output..." );
this->_Output( x, St, invLabels );
}
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(