...

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3e9783b4945d0d66b65d5cc1db9657158424e1fa..b11c77cfe0293f6002d4837e1df9d860b351d492 100644 (file)
--- 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 e4bd49a86ad0abdf4284076d5f21b7f85f5674e8..859b3a6aaffd196fda90d81562c216a4dad528a5 100644 (file)
--- 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 (file)
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 <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>
+
+// -------------------------------------------------------------------------
+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$