From 51f387425fa52c4a92ebf01c9dde3156580f8d40 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Fri, 15 Sep 2017 17:11:45 +0200 Subject: [PATCH] ... --- CMakeLists.txt | 3 + tests/image/RandomWalker/CMakeLists.txt | 1 + tests/image/RandomWalker/Original.cxx | 286 ++++++++++++++++++++++++ 3 files changed, 290 insertions(+) create mode 100644 tests/image/RandomWalker/Original.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index 3e9783b..b11c77c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -103,6 +103,9 @@ if(ITKVtkGlue_LOADED) include(${VTK_USE_FILE}) endif(ITKVtkGlue_LOADED) +find_package(Eigen3 CONFIG REQUIRED) +include(${EIGEN3_USE_FILE}) + ## ========================= ## == Installation values == ## ========================= diff --git a/tests/image/RandomWalker/CMakeLists.txt b/tests/image/RandomWalker/CMakeLists.txt index e4bd49a..859b3a6 100644 --- a/tests/image/RandomWalker/CMakeLists.txt +++ b/tests/image/RandomWalker/CMakeLists.txt @@ -7,6 +7,7 @@ set(_pfx test_fpa_Image_RandomWalker_) set( _tests Gaussian + Original ) include_directories(${PROJECT_SOURCE_DIR}/lib ${PROJECT_BINARY_DIR}/lib) diff --git a/tests/image/RandomWalker/Original.cxx b/tests/image/RandomWalker/Original.cxx new file mode 100644 index 0000000..8b7bca1 --- /dev/null +++ b/tests/image/RandomWalker/Original.cxx @@ -0,0 +1,286 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +// ------------------------------------------------------------------------- +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 > 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 ); +} + +// ------------------------------------------------------------------------- +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; + + 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 }, + }; + + // 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 ) + { + 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 ) + { + 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 + */ + + return( 0 ); +} + +// eof - $RCSfile$ -- 2.47.1