...

[FrontAlgorithms.git] / lib / fpa / Common / OriginalRandomWalker.hxx

diff --git a/lib/fpa/Common/OriginalRandomWalker.hxx b/lib/fpa/Common/OriginalRandomWalker.hxx
index 46dd9354db02015ad64fdd718db5ea37dfd929b1..2d3dc73a58ea6592175395bf70beebfe678b0d4b 100644 (file)
--- a/lib/fpa/Common/OriginalRandomWalker.hxx
+++ b/lib/fpa/Common/OriginalRandomWalker.hxx
@@ -58,57 +58,80 @@ GenerateData( )
   // 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 );
 }