From 524d4f16dc12dc2f07df32a55192c76d810ac8c0 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Mon, 18 Sep 2017 17:28:39 +0200 Subject: [PATCH] ... --- lib/fpa/Common/OriginalRandomWalker.h | 92 +++++++ lib/fpa/Common/OriginalRandomWalker.hxx | 328 ++++++++++++++++++++++++ lib/fpa/Config.h.in | 7 + tests/image/RandomWalker/Original.cxx | 317 +++++------------------ 4 files changed, 487 insertions(+), 257 deletions(-) create mode 100644 lib/fpa/Common/OriginalRandomWalker.h create mode 100644 lib/fpa/Common/OriginalRandomWalker.hxx diff --git a/lib/fpa/Common/OriginalRandomWalker.h b/lib/fpa/Common/OriginalRandomWalker.h new file mode 100644 index 0000000..915cb19 --- /dev/null +++ b/lib/fpa/Common/OriginalRandomWalker.h @@ -0,0 +1,92 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= +#ifndef __fpa__Common__OriginalRandomWalker__h__ +#define __fpa__Common__OriginalRandomWalker__h__ + +#include +#include +#include +#include + +namespace fpa +{ + namespace Common + { + /** + */ + template< class _TImage, class _TLabels, class _TScalar = float > + class OriginalRandomWalker + : public itk::ImageToImageFilter< _TImage, _TLabels > + { + public: + typedef _TImage TImage; + typedef _TLabels TLabels; + typedef _TScalar TScalar; + + typedef OriginalRandomWalker Self; + typedef itk::ImageToImageFilter< TImage, TLabels > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef itk::Image< TScalar, TImage::ImageDimension > TScalarImage; + typedef typename TImage::IndexType TIndex; + typedef typename TImage::PixelType TPixel; + typedef typename TImage::RegionType TRegion; + typedef typename TLabels::PixelType TLabel; + + public: + itkNewMacro( Self ); + itkTypeMacro( + 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 ); + + fpaFilterInputMacro( InputLabels, TLabels ); + fpaFilterOutputMacro( OutputProbabilities, TScalarImage ); + + public: + void AddSeed( const TIndex& seed, const TLabel& label ); + + protected: + OriginalRandomWalker( ); + virtual ~OriginalRandomWalker( ); + + virtual void GenerateData( ) override; + + static unsigned long _1D( const TIndex& idx, const TRegion& region ); + static TIndex _ND( const unsigned long& i, const TRegion& region ); + + private: + // Purposely not implemented + OriginalRandomWalker( const Self& other ); + Self& operator=( const Self& other ); + + protected: + std::vector< TIndex > m_Seeds; + std::vector< TLabel > m_Labels; + + TScalar m_Beta; + TScalar m_Epsilon; + bool m_NormalizeWeights; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION +#endif // __fpa__Common__OriginalRandomWalker__h__ +// eof - $RCSfile$ diff --git a/lib/fpa/Common/OriginalRandomWalker.hxx b/lib/fpa/Common/OriginalRandomWalker.hxx new file mode 100644 index 0000000..7c29595 --- /dev/null +++ b/lib/fpa/Common/OriginalRandomWalker.hxx @@ -0,0 +1,328 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= +#ifndef __fpa__Common__OriginalRandomWalker__hxx__ +#define __fpa__Common__OriginalRandomWalker__hxx__ + +#include +#include +#include +#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( ); +} + +// ------------------------------------------------------------------------- +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 ) +{ + fpaFilterOptionalInputConfigureMacro( InputLabels, TLabels ); + fpaFilterOutputConfigureMacro( OutputProbabilities, TScalarImage ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TLabels, class _TScalar > +fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: +~OriginalRandomWalker( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TLabels, class _TScalar > +void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: +GenerateData( ) +{ + typedef Eigen::SparseMatrix< TScalar > _TSparseMatrix; + typedef Eigen::SimplicialLDLT< _TSparseMatrix > _TSparseSolver; + typedef Eigen::Triplet< TScalar > _TTriplet; + typedef std::vector< _TTriplet > _TTriplets; + typedef std::map< unsigned long, unsigned long > _TMap; + typedef _TMap::value_type _TMapValue; + + // Some input values + const TImage* in = this->GetInput( ); + TRegion region = in->GetRequestedRegion( ); + unsigned long N = region.GetNumberOfPixels( ); + + // Prepare seeds for linear algebra + _TMap seeds; + for( unsigned int i = 0; i < this->m_Seeds.size( ); ++i ) + seeds.insert( + _TMapValue( + Self::_1D( this->m_Seeds[ i ], region ), + this->m_Labels[ i ] + ) + ); + + // Use input labels + /* TODO + const TLabels* in_labels = this->GetInputLabels( ); + if( in_labels != NULL ) + { + itk::ImageRegionConstIteratorWithIndex< TLabels > lIt( in_labels, region ); + for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt ) + if( lIt.Get( ) != 0 ) + seeds.insert( + _TMapValue( Self::_1D( lIt.GetIndex( ), region ), lIt.Get( ) ) + ); + + } // fi + */ + + // Prepare label tables + _TMap labels, inv_labels, seeds_indexes; + _TMap::const_iterator hIt = seeds.begin( ); + for( ; hIt != seeds.end( ); ++hIt ) + { + seeds_indexes[ hIt->first ] = seeds_indexes.size( ) - 1; + if( labels.find( hIt->second ) == labels.end( ) ) + { + labels[ hIt->second ] = labels.size( ) - 1; + inv_labels[ labels[ hIt->second ] ] = hIt->second; + + } // fi + + } // rof + + // Prepare matrix/image index conversion + _TMap indexes; + unsigned long o = 0; + for( unsigned long n = 0; n < N; ++n ) + { + if( seeds.find( n ) == seeds.end( ) ) + indexes.insert( _TMap::value_type( n, n - o ) ); + else + o++; + + } // rof + unsigned long nLabels = labels.size( ); + + // Boundary matrix + _TTriplets Bt; + _TMap::const_iterator mIt = seeds.begin( ); + for( unsigned long i = 0; mIt != seeds.end( ); ++mIt, ++i ) + Bt.push_back( _TTriplet( i, labels[ mIt->second ], TScalar( 1 ) ) ); + _TSparseMatrix B( seeds.size( ), nLabels ); + B.setFromTriplets( Bt.begin( ), Bt.end( ) ); + B.makeCompressed( ); + + // Laplacian matrix + _TTriplets Lt; + itk::ImageRegionConstIteratorWithIndex< TImage > it( in, region ); + TScalar maxV = TScalar( -1 ); + for( it.GoToBegin( ); !it.IsAtEnd( ); ++it ) + { + TIndex idx = it.GetIndex( ); + TScalar vidx = TScalar( it.Get( ) ); + unsigned long iidx = Self::_1D( idx, region ); + + // Neighbors + TScalar s = TScalar( 0 ); + for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + { + TIndex jdx; + for( int l = -1; l <= 1; l += 2 ) + { + jdx = idx; + jdx[ d ] += l; + if( region.IsInside( jdx ) ) + { + TScalar vjdx = TScalar( in->GetPixel( jdx ) ); + unsigned long ijdx = Self::_1D( jdx, region ); + TScalar v = std::fabs( vidx - vjdx ); + Lt.push_back( _TTriplet( iidx, ijdx, v ) ); + if( maxV < v ) + maxV = v; + + } // fi + + } // rof + + } // rof + + } // rof + std::vector< TScalar > diag( N, TScalar( 0 ) ); + typename _TTriplets::iterator tIt; + TScalar betaV = -this->m_Beta; + if( this->m_NormalizeWeights ) + betaV /= maxV; + for( tIt = Lt.begin( ); tIt != Lt.end( ); ++tIt ) + { + TScalar v = std::exp( betaV * tIt->value( ) ); + if( v < this->m_Epsilon ) + v = this->m_Epsilon; + *tIt = _TTriplet( tIt->row( ), tIt->col( ), -v ); + diag[ tIt->col( ) ] += v; + + } // rof + for( unsigned long i = 0; i < diag.size( ); ++i ) + Lt.push_back( _TTriplet( i, i, diag[ i ] ) ); + + // Compute R and A + _TTriplets Rt, At; + for( tIt = Lt.begin( ); tIt != Lt.end( ); ++tIt ) + { + _TMap::const_iterator cIt, rIt; + cIt = seeds.find( tIt->col( ) ); + rIt = seeds.find( tIt->row( ) ); + if( cIt != seeds.end( ) ) + { + if( rIt == seeds.end( ) ) + { + _TMap::const_iterator iIt = indexes.find( tIt->row( ) ); + _TMap::const_iterator jIt = seeds_indexes.find( cIt->first ); + Rt.push_back( _TTriplet( iIt->second, jIt->second, -tIt->value( ) ) ); + + } // fi + } + else + { + if( rIt == seeds.end( ) ) + { + _TMap::const_iterator iIt = indexes.find( tIt->row( ) ); + _TMap::const_iterator jIt = indexes.find( tIt->col( ) ); + At.push_back( _TTriplet( iIt->second, jIt->second, tIt->value( ) ) ); + + } // fi + + } // fi + + } // rof + + _TSparseMatrix R( N - seeds.size( ), seeds.size( ) ); + R.setFromTriplets( Rt.begin( ), Rt.end( ) ); + R.makeCompressed( ); + + _TSparseMatrix A( N - seeds.size( ), N - seeds.size( ) ); + A.setFromTriplets( At.begin( ), At.end( ) ); + A.makeCompressed( ); + + // Solve dirichlet problem + _TSparseSolver solver; + solver.compute( A ); + if( solver.info( ) != Eigen::Success ) + { + std::cerr << "Error computing." << std::endl; + } // fi + _TSparseMatrix x = solver.solve( R * B ); + if( solver.info( ) != Eigen::Success ) + { + std::cerr << "Error solving." << std::endl; + } // fi + + // Fill outputs + TLabels* out_labels = this->GetOutput( ); + out_labels->SetLargestPossibleRegion( in->GetLargestPossibleRegion( ) ); + out_labels->SetRequestedRegion( in->GetRequestedRegion( ) ); + out_labels->SetBufferedRegion( in->GetBufferedRegion( ) ); + out_labels->SetSpacing( in->GetSpacing( ) ); + out_labels->SetOrigin( in->GetOrigin( ) ); + out_labels->SetDirection( in->GetDirection( ) ); + out_labels->Allocate( ); + + TScalarImage* out_probs = this->GetOutputProbabilities( ); + out_probs->SetLargestPossibleRegion( in->GetLargestPossibleRegion( ) ); + out_probs->SetRequestedRegion( in->GetRequestedRegion( ) ); + out_probs->SetBufferedRegion( in->GetBufferedRegion( ) ); + out_probs->SetSpacing( in->GetSpacing( ) ); + out_probs->SetOrigin( in->GetOrigin( ) ); + out_probs->SetDirection( in->GetDirection( ) ); + out_probs->Allocate( ); + + mIt = seeds.begin( ); + o = 0; + unsigned long j = 0; + for( unsigned long j = 0; j < N; ++j ) + { + TIndex idx; + TLabel lbl; + TScalar p; + if( mIt->first != j ) + { + idx = Self::_ND( j, region ); + TScalar maxP = x.coeffRef( j - o, 0 ); + unsigned long maxL = 0; + for( unsigned long l = 1; l < nLabels; ++l ) + { + TScalar vp = x.coeffRef( j - o, l ); + if( maxP < vp ) + { + maxP = vp; + maxL = l; + + } // fi + + } // rof + lbl = inv_labels[ maxL ]; + p = maxP; + } + else + { + idx = Self::_ND( mIt->first, region ); + lbl = mIt->second; + p = TScalar( 1 ); + mIt++; + o++; + + } // fi + out_labels->SetPixel( idx, lbl ); + out_probs->SetPixel( idx, p ); + + } // rof +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TLabels, class _TScalar > +unsigned long +fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: +_1D( const TIndex& idx, const TRegion& region ) +{ + unsigned long i = idx[ 0 ]; + unsigned long off = 1; + typename TRegion::SizeType size = region.GetSize( ); + for( unsigned int d = 1; d < TIndex::Dimension; ++d ) + { + off *= size[ d - 1 ]; + i += idx[ d ] * off; + + } // rof + return( i ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TLabels, class _TScalar > +typename fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: +TIndex fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >:: +_ND( const unsigned long& i, const TRegion& region ) +{ + typename TRegion::SizeType size = region.GetSize( ); + + unsigned long j = i; + TIndex idx; + if( TIndex::Dimension == 3 ) + { + unsigned long z = size[ 0 ] * size[ 1 ]; + idx[ 2 ] = j / z; + j -= idx[ 2 ] * z; + + } // fi + idx[ 1 ] = j / size[ 0 ]; + idx[ 0 ] = j % size[ 0 ]; + return( idx ); +} + +#endif // __fpa__Common__OriginalRandomWalker__hxx__ +// eof - $RCSfile$ diff --git a/lib/fpa/Config.h.in b/lib/fpa/Config.h.in index 831fb72..e850355 100644 --- a/lib/fpa/Config.h.in +++ b/lib/fpa/Config.h.in @@ -53,6 +53,13 @@ this->m_##__n__##Idx + 1 \ ) +// ------------------------------------------------------------------------- +#define fpaFilterOptionalInputConfigureMacro( __n__, __t__ ) \ + this->m_##__n__##Idx = this->GetNumberOfIndexedInputs( ); \ + this->itk::ProcessObject::SetNumberOfIndexedInputs( \ + this->m_##__n__##Idx + 1 \ + ) + // ------------------------------------------------------------------------- #define fpaFilterOutputMacro( __n__, __t__ ) \ private: \ diff --git a/tests/image/RandomWalker/Original.cxx b/tests/image/RandomWalker/Original.cxx index 8b7bca1..5a6c6df 100644 --- a/tests/image/RandomWalker/Original.cxx +++ b/tests/image/RandomWalker/Original.cxx @@ -3,283 +3,86 @@ // @email florez-l@javeriana.edu.co // ========================================================================= -#include -#include -#include #include #include #include -#include -#include - -#include -#include -#include +#include // ------------------------------------------------------------------------- -const unsigned int Dim = 2; +const unsigned int Dim = 3; typedef unsigned char TPixel; typedef unsigned char TLabel; -typedef double TScalar; -typedef itk::Image< TPixel, Dim > TImage; -typedef itk::Image< TLabel, Dim > TLabelImage; - -typedef Eigen::SparseMatrix< TScalar > TSparseMatrix; - -/* TODO - class Eigen::AMDOrdering< StorageIndex > - class Eigen::COLAMDOrdering< StorageIndex > - class Eigen::NaturalOrdering< StorageIndex > -*/ -typedef Eigen::AMDOrdering< int > TSolverStorage; - -typedef Eigen::SparseQR< TSparseMatrix, TSolverStorage > TSparseSolver; -typedef Eigen::Triplet< TScalar > TTriplet; -typedef std::vector< TTriplet > TTriplets; - -// ------------------------------------------------------------------------- -template< class _TIndex, class _TRegion > -unsigned long GetIndex( const _TIndex& idx, const _TRegion& region ) -{ - unsigned long i = 0; - unsigned long off = 1; - for( unsigned int d = 0; d < _TIndex::Dimension; ++d ) - { - i += idx[ d ] * off; - off *= region.GetSize( )[ d ]; - } // rof - return( i ); -} +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::Image< TLabel, Dim > TLabels; // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) { // Get arguments - /* TODO - if( argc < 2 ) - { - std::cerr - << "Usage: " << argv[ 0 ] << " input_image" - << std::endl; - return( 1 ); - - } // fi - std::string input_image_filename = argv[ 1 ]; - */ - TScalar beta = 90; - TScalar eps = 1e-5; + if( argc < 4 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image output_labels output_probabilities" + << std::endl; + return( 1 ); - TPixel img[ 5 ][ 5 ] = - { - { 1, 2 , 100, 2 , 1 }, - { 1, 100, 100, 100, 1 }, - { 1, 100, 100, 100, 1 }, - { 1, 100, 100, 100, 1 }, - { 1, 2 , 100, 2 , 1 }, - }; + } // fi + std::string input_image_filename = argv[ 1 ]; + std::string output_labels_filename = argv[ 2 ]; + std::string output_probabilities_filename = argv[ 3 ]; + double beta = 90; + double eps = 1e-5; + + TImage::IndexType s1, s2; + s1.Fill( 0 ); + s2.Fill( 0 ); + s1[ 0 ] = 131; + s1[ 1 ] = 150; + s2[ 0 ] = 200; + s2[ 1 ] = 200; // Read image TImage::Pointer input; - { // begin - TImage::SizeType size; - size.Fill( 5 ); - - input = TImage::New( ); - input->SetRegions( size ); - input->Allocate( ); - itk::ImageRegionIteratorWithIndex< TImage > it( input, input->GetLargestPossibleRegion( ) ); - for( it.GoToBegin( ); !it.IsAtEnd( ); ++it ) - it.Set( img[ it.GetIndex( )[ 1 ] ][ it.GetIndex( )[ 0 ] ] ); - - /* TODO - typedef itk::ImageFileReader< TImage > TImageReader; - TImageReader::Pointer input_image_reader = TImageReader::New( ); - input_image_reader->SetFileName( input_image_filename ); - input_image_reader->Update( ); - input = input_image_reader->GetOutput( ); - input->DisconnectPipeline( ); - */ - } // end - TImage::RegionType region = input->GetLargestPossibleRegion( ); - unsigned long N = region.GetNumberOfPixels( ); - - // Seeds - TImage::IndexType s1, s2, s3; - s1[ 0 ] = 0; - s1[ 1 ] = 1; - s2[ 0 ] = 2; - s2[ 1 ] = 3; - s3[ 0 ] = 4; - s3[ 1 ] = 4; - - /* - std::vector< unsigned long > seeds; - std::vector< unsigned long > labels; - */ - std::map< unsigned long, unsigned long > seeds; - seeds[ GetIndex( s1, region ) ] = 1; - seeds[ GetIndex( s2, region ) ] = 2; - seeds[ GetIndex( s3, region ) ] = 2; - unsigned long nLabels = 2; - - // Construct L - TTriplets Lt; - itk::ImageRegionIteratorWithIndex< TImage > it( input, region ); - TScalar maxV = -1; - for( it.GoToBegin( ); !it.IsAtEnd( ); ++it ) + typedef itk::ImageFileReader< TImage > TImageReader; + TImageReader::Pointer input_image_reader = TImageReader::New( ); + input_image_reader->SetFileName( input_image_filename ); + + // Random walker + typedef fpa::Common::OriginalRandomWalker< TImage, TLabels > TFilter; + TFilter::Pointer rw = TFilter::New( ); + rw->SetInput( input_image_reader->GetOutput( ) ); + rw->SetBeta( beta ); + rw->SetEpsilon( eps ); + rw->AddSeed( s1, 100 ); + rw->AddSeed( s2, 150 ); + rw->NormalizeWeightsOn( ); + + // Save labels + typedef itk::ImageFileWriter< TFilter::TLabels > TLabelsWriter; + TLabelsWriter::Pointer output_labels_writer = TLabelsWriter::New( ); + output_labels_writer->SetFileName( output_labels_filename ); + output_labels_writer->SetInput( rw->GetOutput( ) ); + + // Save probabilities + typedef itk::ImageFileWriter< TFilter::TScalarImage > TProbabilitesWriter; + TProbabilitesWriter::Pointer output_probabilities_writer = + TProbabilitesWriter::New( ); + output_probabilities_writer->SetFileName( output_probabilities_filename ); + output_probabilities_writer->SetInput( rw->GetOutputProbabilities( ) ); + + try { - TImage::IndexType idx = it.GetIndex( ); - TScalar vidx = TScalar( it.Get( ) ); - unsigned long iidx = GetIndex( idx, region ); - - // Neighbors - TScalar s = TScalar( 0 ); - for( unsigned int d = 0; d < Dim; ++d ) - { - TImage::IndexType jdx; - for( int l = -1; l <= 1; l += 2 ) - { - jdx = idx; - jdx[ d ] += l; - if( region.IsInside( jdx ) ) - { - TScalar vjdx = TScalar( input->GetPixel( jdx ) ); - unsigned long ijdx = GetIndex( jdx, region ); - TScalar v = std::fabs( vidx - vjdx ); - Lt.push_back( TTriplet( iidx, ijdx, v ) ); - if( maxV < v ) - maxV = v; - - } // fi - - } // rof - - } // rof - - } // rof - - std::vector< TScalar > diag( N, TScalar( 0 ) ); - for( TTriplets::iterator tIt = Lt.begin( ); tIt != Lt.end( ); ++tIt ) - { - TScalar v = std::exp( -beta * tIt->value( ) / maxV ); - if( v < eps ) - v = eps; - *tIt = TTriplet( tIt->row( ), tIt->col( ), -v ); - diag[ tIt->col( ) ] += v; - - } // rof - for( unsigned long i = 0; i < diag.size( ); ++i ) - Lt.push_back( TTriplet( i, i, diag[ i ] ) ); - TSparseMatrix L( N, N ); - L.setFromTriplets( Lt.begin( ), Lt.end( ) ); - L.makeCompressed( ); - - // Boundary - TTriplets boundaryt; - // TODO: for( unsigned int i = 0; i < seeds.size( ); ++i ) - std::map< unsigned long, unsigned long >::const_iterator mIt = seeds.begin( ); - for( unsigned long i = 0; mIt != seeds.end( ); ++mIt, ++i ) - boundaryt.push_back( TTriplet( i, mIt->second - 1, TScalar( 1 ) ) ); - TSparseMatrix boundary( seeds.size( ), nLabels ); - boundary.setFromTriplets( boundaryt.begin( ), boundaryt.end( ) ); - boundary.makeCompressed( ); - - // Compute RHS - TTriplets RHSt; - unsigned long x = 0; - std::map< unsigned long, unsigned long >::const_iterator sIt = seeds.begin( ); - for( ; sIt != seeds.end( ); ++sIt ) - { - unsigned long y = 0; - for( unsigned long n = 0; n < N; ++n ) - { - std::map< unsigned long, unsigned long >::const_iterator nIt = - seeds.find( n ); - if( nIt == seeds.end( ) ) - RHSt.push_back( TTriplet( y++, x, L.coeff( sIt->first, n ) ) ); - - } // rof - x++; - - } // rof - TSparseMatrix RHS( N - seeds.size( ), seeds.size( ) ); - RHS.setFromTriplets( RHSt.begin( ), RHSt.end( ) ); - RHS.makeCompressed( ); - - // Compute A - TTriplets At; - x = 0; - for( unsigned long m = 0; m < N; ++m ) - { - if( seeds.find( m ) == seeds.end( ) ) - { - unsigned long y = 0; - for( unsigned long n = 0; n < N; ++n ) - if( seeds.find( n ) == seeds.end( ) ) - At.push_back( TTriplet( y++, x, L.coeff( m, n ) ) ); - x++; - - } // fi - - } // rof - TSparseMatrix A( N - seeds.size( ), N - seeds.size( ) ); - A.setFromTriplets( At.begin( ), At.end( ) ); - A.makeCompressed( ); - - // Solve dirichlet problem - TSparseSolver solver; - solver.compute( A ); - if( solver.info( ) != Eigen::Success ) + output_labels_writer->Update( ); + output_probabilities_writer->Update( ); + } + catch( std::exception& err ) { - std::cerr << "Error computing." << std::endl; - } // fi - TSparseMatrix sol = solver.solve( ( RHS * TScalar( -1 ) ) * boundary ); - if( solver.info( ) != Eigen::Success ) - { - std::cerr << "Error solving." << std::endl; - } // fi - - TLabelImage::Pointer output = TLabelImage::New( ); - output->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) ); - output->SetRequestedRegion( input->GetRequestedRegion( ) ); - output->SetBufferedRegion( input->GetBufferedRegion( ) ); - output->SetSpacing( input->GetSpacing( ) ); - output->SetOrigin( input->GetOrigin( ) ); - output->SetDirection( input->GetDirection( ) ); - output->Allocate( ); - - /* TODO - itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( output, output->GetRequestedRegion( ) ); - for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt ) - { - unsigned long i = GetIndex( lIt.GetIndex( ), region ); - std::map< unsigned long, unsigned long >::const_iterator sIt = - seeds.find( i ); - if( sIt == seeds.end( ) ) - { - std::cout << i << std::endl; - TScalar maxProb = TScalar( -1 ); - for( unsigned long k = 0; k < nLabels; ++k ) - { - std::cout << "\t\t" << i << " " << k << std::endl; - TScalar p = sol.coeff( i, k ); - if( maxProb < p ) - maxProb = p; - - } // rof - std::cout << "\t" << maxProb << std::endl; - } - else - { - std::cout << "---> " << i << std::endl; - lIt.Set( sIt->second ); - - } // fi - - } // rof - */ + std::cerr << "Error caught: " << err.what( ) << std::endl; + return( 1 ); + } // yrt return( 0 ); } -- 2.45.1