]> Creatis software - FrontAlgorithms.git/commitdiff
...
authorLeonardo Flórez-Valencia <leonardo.florez@gmail.com>
Thu, 21 Sep 2017 09:15:33 +0000 (11:15 +0200)
committerLeonardo Flórez-Valencia <leonardo.florez@gmail.com>
Thu, 21 Sep 2017 09:15:33 +0000 (11:15 +0200)
lib/fpa/Common/OriginalRandomWalker.h
lib/fpa/Common/OriginalRandomWalker.hxx
lib/fpa/Functors/Dijkstra/Image/Gaussian.h
tests/image/RandomWalker/Original.cxx

index d4af732998028d0b9714391972d5432b098669c1..0c3a9f82f08a8219a10480b1062415a9d526dcf7 100644 (file)
@@ -11,6 +11,7 @@
 #include <itkImage.h>
 #include <itkImageToImageFilter.h>
 #include <itkSimpleFastMutexLock.h>
+#include <fpa/Functors/VertexFunction.h>
 
 namespace fpa
 {
@@ -38,6 +39,8 @@ namespace fpa
       typedef typename TImage::RegionType TRegion;
       typedef typename TLabels::PixelType TLabel;
 
+      typedef fpa::Functors::VertexFunction< TIndex, TScalar > TEdgeFunction;
+
     protected:
       struct _TBoundaryThreadStruct
       {
@@ -68,21 +71,17 @@ namespace fpa
         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 );
+      itkGetConstObjectMacro( EdgeFunction, TEdgeFunction );
+      itkGetObjectMacro( EdgeFunction, TEdgeFunction );
+      itkSetObjectMacro( EdgeFunction, TEdgeFunction );
 
       fpaFilterInputMacro( InputLabels, TLabels );
       fpaFilterOutputMacro( OutputProbabilities, TScalarImage );
 
-    public:
-      void AddSeed( const TIndex& seed, const TLabel& label );
+      /* TODO
+         public:
+         void AddSeed( const TIndex& seed, const TLabel& label );
+      */
 
     protected:
       OriginalRandomWalker( );
@@ -90,7 +89,6 @@ namespace fpa
 
       virtual void GenerateData( ) override;
 
-      virtual TScalar _W( const TIndex& i, const TIndex& j );
       virtual TScalar _L( const TIndex& i, const TIndex& j );
 
       // Boundary construction
@@ -161,12 +159,12 @@ namespace fpa
       Self& operator=( const Self& other );
 
     protected:
-      std::vector< TIndex > m_Seeds;
-      std::vector< TLabel > m_Labels;
+      /* TODO
+         std::vector< TIndex > m_Seeds;
+         std::vector< TLabel > m_Labels;
+      */
 
-      TScalar m_Beta;
-      TScalar m_Epsilon;
-      bool m_NormalizeWeights;
+      typename TEdgeFunction::Pointer m_EdgeFunction;
 
       itk::SimpleFastMutexLock m_Mutex;
     };
index 9bc08d4c5c461d802deb2de6200f4101d47518e7..4b60c1bb65302529b9ccf0d304f71b2c126210aa 100644 (file)
 #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( );
-}
+/* TODO
+   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 )
+  : Superclass( )
 {
   fpaFilterInputConfigureMacro( InputLabels, TLabels );
   fpaFilterOutputConfigureMacro( OutputProbabilities, TScalarImage );
@@ -51,6 +50,12 @@ GenerateData( )
   typedef Eigen::SparseMatrix< TScalar >    _TMatrix;
   typedef Eigen::SimplicialLDLT< _TMatrix > _TSolver;
 
+  // Configure edge function
+  if( this->m_EdgeFunction.IsNull( ) )
+    itkExceptionMacro( << "Undefined edge function." );
+  const TImage* input = this->GetInput( );
+  this->m_EdgeFunction->SetDataObject( input );
+
   // Allocate outputs
   this->AllocateOutputs( );
 
@@ -74,7 +79,7 @@ GenerateData( )
   this->_Laplacian( At, Rt, St );
 
   // Matrices
-  TRegion region = this->GetInput( )->GetRequestedRegion( );
+  TRegion region = input->GetRequestedRegion( );
   unsigned long nSeeds = St.size( );
   unsigned long nLabels = labels.size( );
   unsigned long N = region.GetNumberOfPixels( );
@@ -99,34 +104,15 @@ GenerateData( )
   _TSolver solver;
   solver.compute( A );
   if( solver.info( ) != Eigen::Success )
-  {
-    std::cerr << "Error computing." << std::endl;
-  } // fi
+    itkExceptionMacro( << "Error decomposing matrix." );
   _TMatrix x = solver.solve( R * B );
   if( solver.info( ) != Eigen::Success )
-  {
-    std::cerr << "Error solving." << std::endl;
-  } // fi
+    itkExceptionMacro( << "Error solving system." );
 
   // Fill outputs
   this->_Output( x, St, invLabels );
 }
 
-// -------------------------------------------------------------------------
-template< class _TImage, class _TLabels, class _TScalar >
-_TScalar fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
-_W( const TIndex& i, const TIndex& j )
-{
-  const TImage* in = this->GetInput( );
-  TScalar a = TScalar( in->GetPixel( i ) );
-  TScalar b = TScalar( in->GetPixel( j ) );
-  TScalar v = std::exp( -this->m_Beta * std::fabs( a - b ) );
-  if( v < this->m_Epsilon )
-    return( this->m_Epsilon );
-  else
-    return( v );
-}
-
 // -------------------------------------------------------------------------
 template< class _TImage, class _TLabels, class _TScalar >
 _TScalar fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
@@ -143,7 +129,7 @@ _L( const TIndex& i, const TIndex& j )
         TIndex k = i;
         k[ d ] += n;
         if( r.IsInside( k ) )
-          s += this->_W( i, k );
+          s += this->m_EdgeFunction->Evaluate( i, k );
 
       } // rof
 
@@ -151,7 +137,7 @@ _L( const TIndex& i, const TIndex& j )
     return( s );
   }
   else
-    return( this->_W( i, j ) * TScalar( -1 ) );
+    return( -( this->m_EdgeFunction->Evaluate( i, j ) ) );
 }
 
 // -------------------------------------------------------------------------
@@ -317,13 +303,13 @@ _ThreadedLaplacian(
 
   const TImage* in = this->GetInput( );
   const TLabels* in_labels = this->GetInputLabels( );
-  TRegion rqRegion = in->GetRequestedRegion( );
+  TRegion reqRegion = in->GetRequestedRegion( );
   _TIt it( in, region );
   for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
   {
     TIndex idx = it.GetIndex( );
     bool iSeed = ( in_labels->GetPixel( idx ) != 0 );
-    unsigned long i = Self::_1D( idx, rqRegion );
+    unsigned long i = Self::_1D( idx, reqRegion );
     unsigned long si;
 
     // A's diagonal values
@@ -344,10 +330,10 @@ _ThreadedLaplacian(
       {
         TIndex jdx = idx;
         jdx[ d ] += s;
-        if( rqRegion.IsInside( jdx ) )
+        if( reqRegion.IsInside( jdx ) )
         {
           TScalar L = this->_L( idx, jdx );
-          unsigned long j = Self::_1D( jdx, rqRegion );
+          unsigned long j = Self::_1D( jdx, reqRegion );
           bool jSeed = ( in_labels->GetPixel( jdx ) != 0 );
           if( !jSeed )
           {
@@ -365,7 +351,7 @@ _ThreadedLaplacian(
               this->m_Mutex.Unlock( );
 
             } // fi
-            
+
           } // fi
 
         } // fi
@@ -451,6 +437,7 @@ _ThreadedOutput(
   const TLabels* in_labels = this->GetInputLabels( );
   TLabels* out_labels = this->GetOutput( );
   TScalarImage* out_probs = this->GetOutputProbabilities( );
+  TRegion reqRegion = out_labels->GetRequestedRegion( );
   itk::ImageRegionConstIteratorWithIndex< TLabels > iIt( in_labels, region );
   itk::ImageRegionIteratorWithIndex< TLabels > oIt( out_labels, region );
   itk::ImageRegionIteratorWithIndex< TScalarImage > pIt( out_probs, region );
@@ -461,7 +448,7 @@ _ThreadedOutput(
   {
     if( iIt.Get( ) == 0 )
     {
-      unsigned long i = Self::_1D( iIt.GetIndex( ), region );
+      unsigned long i = Self::_1D( iIt.GetIndex( ), reqRegion );
       unsigned long j = Self::_NearSeedIndex( i, *S );
       TScalar maxP = X->coeff( j, 0 );
       unsigned long maxL = 0;
index e917e95bff25897eacb45e78019ca1b15a6628ae..2effb513f4f707236c3203be83514168b9700420 100644 (file)
@@ -38,14 +38,20 @@ namespace fpa
             fpa::Functors::VertexFunction
             );
 
-        itkGetConstMacro( Alpha, double );
-        itkSetMacro( Alpha, double );
+          itkGetConstMacro( Beta, TValue );
+          itkSetMacro( Beta, TValue );
 
-        itkGetConstMacro( Beta, double );
-        itkSetMacro( Beta, double );
+          itkGetConstMacro( Epsilon, TValue );
+          itkSetMacro( Epsilon, TValue );
+
+          itkBooleanMacro( TreatAsWeight );
+          itkGetConstMacro( TreatAsWeight, bool );
+          itkSetMacro( TreatAsWeight, bool );
 
         public:
-          virtual TValue Evaluate( const TVertex& v, const TVertex& p ) const override
+          virtual TValue Evaluate(
+            const TVertex& v, const TVertex& p
+            ) const override
             {
               const TImage* image =
                 dynamic_cast< const TImage* >(
@@ -53,11 +59,15 @@ namespace fpa
                   );
               if( image != NULL )
               {
-                double d = double( image->GetPixel( v ) );
-                d       -= double( image->GetPixel( p ) );
+                TValue d = TValue( image->GetPixel( v ) );
+                d       -= TValue( image->GetPixel( p ) );
                 d       /= this->m_Beta;
-                d = std::exp( d * d ) - double( 1 );
-                return( TValue( std::pow( d, this->m_Alpha ) ) );
+                d       *= d;
+                if( this->m_TreatAsWeight ) d = std::exp(  d );
+                else                        d = std::exp( -d );
+
+                if( d < this->m_Epsilon ) return( this->m_Epsilon );
+                else                      return( d );
               }
               else
                 return( TValue( -1 ) );
@@ -66,8 +76,9 @@ namespace fpa
         protected:
           Gaussian( )
             : Superclass( ),
-              m_Alpha( double( 1 ) ),
-              m_Beta( double( 1 ) )
+              m_Beta( TValue( 1 ) ),
+              m_Epsilon( TValue( 1e-5 ) ),
+              m_TreatAsWeight( true )
             {
             }
           virtual ~Gaussian( )
@@ -80,8 +91,9 @@ namespace fpa
           Self& operator=( const Self& other );
 
         protected:
-          double m_Alpha;
-          double m_Beta;
+          TValue m_Beta;
+          TValue m_Epsilon;
+          bool   m_TreatAsWeight;
         };
 
       } // ecapseman
index 46ba198914cf8fb822a820984c6561ee865aa52c..0ea2bd0c652ae7043d56e76f01e75f8cf0eb8092 100644 (file)
@@ -7,11 +7,13 @@
 #include <itkImageFileReader.h>
 #include <itkImageFileWriter.h>
 #include <fpa/Common/OriginalRandomWalker.h>
+#include <fpa/Functors/Dijkstra/Image/Gaussian.h>
 
 // -------------------------------------------------------------------------
 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 > TLabels;
@@ -52,14 +54,19 @@ int main( int argc, char* argv[] )
   TLabelsReader::Pointer input_labels_reader = TLabelsReader::New( );
   input_labels_reader->SetFileName( input_labels_filename );
 
+  // Random walker function
+  typedef fpa::Functors::Dijkstra::Image::Gaussian< TImage, TScalar > TFunction;
+  TFunction::Pointer edge = TFunction::New( );
+  edge->SetBeta( beta );
+  edge->SetEpsilon( eps );
+  edge->TreatAsWeightOff( );
+
   // Random walker
-  typedef fpa::Common::OriginalRandomWalker< TImage, TLabels > TFilter;
+  typedef fpa::Common::OriginalRandomWalker< TImage, TLabels, TScalar > TFilter;
   TFilter::Pointer rw = TFilter::New( );
   rw->SetInput( input_image_reader->GetOutput( ) );
   rw->SetInputLabels( input_labels_reader->GetOutput( ) );
-  rw->SetBeta( beta );
-  rw->SetEpsilon( eps );
-  rw->SetNormalizeWeights( normalize );
+  rw->SetEdgeFunction( edge );
 
   // Save labels
   typedef itk::ImageFileWriter< TFilter::TLabels > TLabelsWriter;