--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+#ifndef __fpa__Common__OriginalRandomWalker__h__
+#define __fpa__Common__OriginalRandomWalker__h__
+
+#include <fpa/Config.h>
+#include <vector>
+#include <itkImage.h>
+#include <itkImageToImageFilter.h>
+
+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 <fpa/Common/OriginalRandomWalker.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+#endif // __fpa__Common__OriginalRandomWalker__h__
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+#ifndef __fpa__Common__OriginalRandomWalker__hxx__
+#define __fpa__Common__OriginalRandomWalker__hxx__
+
+#include <cmath>
+#include <map>
+#include <itkImageRegionConstIteratorWithIndex.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( );
+}
+
+// -------------------------------------------------------------------------
+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$
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: \
// @email florez-l@javeriana.edu.co
// =========================================================================
-#include <algorithm>
-#include <map>
-#include <vector>
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
-#include <itkImageRegionIteratorWithIndex.h>
-#include <itkVariableSizeMatrix.h>
-
-#include <Eigen/Sparse>
-#include <Eigen/SparseQR>
-#include <Eigen/OrderingMethods>
+#include <fpa/Common/OriginalRandomWalker.h>
// -------------------------------------------------------------------------
-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 );
}
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff