]> Creatis software - FrontAlgorithms.git/commitdiff
...
authorLeonardo Flórez-Valencia <florez-l@javeriana.edu.co>
Thu, 23 Mar 2017 21:38:44 +0000 (16:38 -0500)
committerLeonardo Flórez-Valencia <florez-l@javeriana.edu.co>
Thu, 23 Mar 2017 21:38:44 +0000 (16:38 -0500)
13 files changed:
data/binary_test_2D_00.png [new file with mode: 0644]
examples/CMakeLists.txt
examples/Dijkstra_Maurer.cxx [new file with mode: 0644]
libs/fpa/Base/Dijkstra.h
libs/fpa/Base/Dijkstra.hxx
libs/fpa/Base/Functors/InvertValue.h [new file with mode: 0644]
libs/fpa/Base/Functors/VertexParentBase.h [new file with mode: 0644]
libs/fpa/Base/MinimumSpanningTree.h [new file with mode: 0644]
libs/fpa/Base/MinimumSpanningTree.hxx [new file with mode: 0644]
libs/fpa/Image/Dijkstra.h [new file with mode: 0644]
libs/fpa/Image/Functors/VertexIdentity.h [new file with mode: 0644]
libs/fpa/Image/Functors/VertexParentBase.h [new file with mode: 0644]
libs/fpa/Image/MinimumSpanningTree.h [new file with mode: 0644]

diff --git a/data/binary_test_2D_00.png b/data/binary_test_2D_00.png
new file mode 100644 (file)
index 0000000..fbe9d39
Binary files /dev/null and b/data/binary_test_2D_00.png differ
index f0c88265d1000a4e2bf26d87941092e5643cb49e..297746a359eea3a2f597ff25fe16789583739238 100644 (file)
@@ -5,6 +5,7 @@ IF(BUILD_EXAMPLES)
     RegionGrow_Tautology
     RegionGrow_BinaryThreshold
     RegionGrow_Mori
+    Dijkstra_Maurer
     #CreateMoriInputImage
     #BronchiiInitialSegmentationWithMori
     #BronchiiInitialSegmentationWithBinaryThresholdRegionGrow
diff --git a/examples/Dijkstra_Maurer.cxx b/examples/Dijkstra_Maurer.cxx
new file mode 100644 (file)
index 0000000..8875959
--- /dev/null
@@ -0,0 +1,120 @@
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <itkMinimumMaximumImageCalculator.h>
+#include <itkSignedMaurerDistanceMapImageFilter.h>
+
+#include <fpa/Image/Dijkstra.h>
+#include <fpa/Image/Functors/VertexIdentity.h>
+#include <fpa/Base/Functors/InvertValue.h>
+
+// -------------------------------------------------------------------------
+static const unsigned int VDim = 2;
+typedef short                                              TPixel;
+typedef double                                             TScalar;
+typedef itk::Image< TPixel, VDim >                         TImage;
+typedef itk::Image< TScalar, VDim >                        TScalarImage;
+typedef itk::ImageFileReader< TImage >                     TReader;
+typedef itk::ImageFileWriter< TScalarImage >               TWriter;
+typedef fpa::Image::Dijkstra< TScalarImage, TScalarImage > TFilter;
+typedef itk::MinimumMaximumImageCalculator< TImage >       TMinMax;
+typedef itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage > TDMap;
+
+typedef fpa::Image::Functors::VertexIdentity< TScalarImage, TScalar > TVertexFunc;
+typedef fpa::Base::Functors::InvertValue< TScalar, TScalar > TValueFunc;
+
+typedef TFilter::TMST                TMST;
+typedef itk::ImageFileWriter< TMST > TMSTWriter;
+
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+  // Get arguments
+  if( argc < 5 + VDim )
+  {
+    std::cerr
+      << "Usage: " << argv[ 0 ]
+      << " input_image output_image output_mst stop_at_one_front";
+    for( unsigned int i = 0; i < VDim; ++i )
+      std::cerr << " s_" << i;
+    std::cerr << " ..." << std::endl;
+    return( 1 );
+
+  } // fi
+  std::string input_image_filename = argv[ 1 ];
+  std::string output_image_filename = argv[ 2 ];
+  std::string output_mst_filename = argv[ 3 ];
+  bool stop_at_one_front = ( std::atoi( argv[ 4 ] ) == 1 );
+
+  TReader::Pointer reader = TReader::New( );
+  reader->SetFileName( input_image_filename );
+  try
+  {
+    reader->Update( );
+  }
+  catch( std::exception& err )
+  {
+    std::cerr << "Error caught: " << err.what( ) << std::endl;
+    return( 1 );
+
+  } // yrt
+
+  TMinMax::Pointer minmax = TMinMax::New( );
+  minmax->SetImage( reader->GetOutput( ) );
+  minmax->Compute( );
+
+  TDMap::Pointer dmap = TDMap::New( );
+  dmap->SetInput( reader->GetOutput( ) );
+  dmap->SetBackgroundValue( minmax->GetMinimum( ) );
+  dmap->InsideIsPositiveOn( );
+  dmap->UseImageSpacingOn( );
+  dmap->SquaredDistanceOff( );
+
+  TFilter::Pointer filter = TFilter::New( );
+  filter->SetInput( dmap->GetOutput( ) );
+  filter->SetStopAtOneFront( stop_at_one_front );
+
+  for( int i = 5; i < argc; i += VDim )
+  {
+    if( i + VDim <= argc )
+    {
+      TImage::IndexType seed;
+      for( int j = 0; j < VDim; ++j )
+        seed[ j ] = std::atoi( argv[ i + j ] );
+      filter->AddSeed( seed );
+
+    } // fi
+
+  } // rof
+
+  TVertexFunc::Pointer vertex_func = TVertexFunc::New( );
+  filter->SetFunctor( vertex_func );
+
+  TValueFunc::Pointer value_func = TValueFunc::New( );
+  value_func->SetAlpha( 1 );
+  value_func->SetBeta( 1 );
+  filter->SetFunctor( value_func );
+
+  TWriter::Pointer writer = TWriter::New( );
+  writer->SetInput( filter->GetOutput( ) );
+  writer->SetFileName( output_image_filename );
+
+  TMSTWriter::Pointer mst_writer = TMSTWriter::New( );
+  mst_writer->SetInput( filter->GetMinimumSpanningTree( ) );
+  mst_writer->SetFileName( output_mst_filename );
+
+  try
+  {
+    writer->Update( );
+    mst_writer->Update( );
+  }
+  catch( std::exception& err )
+  {
+    std::cerr << "ERROR: " << err.what( ) << std::endl;
+    return( 1 );
+
+  } // yrt
+  return( 0 );
+}
+
+// eof - $RCSfile$
index 2c49c06b3ff9b79201c0c55dd83568fd03a0c042..d231c7c41da897886b69c63d97f0fca117908646 100644 (file)
@@ -7,6 +7,8 @@
 #define __fpa__Base__Dijkstra__h__
 
 #include <itkFunctionBase.h>
+#include <fpa/Base/MinimumSpanningTree.h>
+#include <fpa/Base/Functors/VertexParentBase.h>
 
 namespace fpa
 {
@@ -14,7 +16,7 @@ namespace fpa
   {
     /**
      */
-    template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
+    template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
     class Dijkstra
       : public _TFilter,
         public _TMarksInterface,
@@ -25,6 +27,7 @@ namespace fpa
       typedef _TFilter                        Superclass;
       typedef _TMarksInterface                TMarksInterface;
       typedef _TSeedsInterface                TSeedsInterface;
+      typedef _TMST                           TMST;
       typedef itk::SmartPointer< Self >       Pointer;
       typedef itk::SmartPointer< const Self > ConstPointer;
 
@@ -34,7 +37,7 @@ namespace fpa
       typedef typename Superclass::TVertices    TVertices;
 
       typedef itk::FunctionBase< TInputValue, TOutputValue > TIntensityFunctor;
-      typedef itk::FunctionBase< TVertex, TOutputValue >     TVertexFunctor;
+      typedef fpa::Base::Functors::VertexParentBase< TVertex, TOutputValue > TVertexFunctor;
 
     protected:
       struct _TNode
@@ -42,15 +45,17 @@ namespace fpa
         TVertex Vertex;
         TVertex Parent;
         TOutputValue Cost;
-        _TNode( const TVertex& v, const TVertex& p )
+        unsigned long FrontId;
+        _TNode( const TVertex& v, const TVertex& p, const unsigned long& fId )
           {
             this->Vertex = v;
             this->Parent = p;
+            this->FrontId = fId;
             this->Cost = TOutputValue( 0 );
           }
         bool operator<( const _TNode& b ) const
           {
-            return( this->Cost < b.Cost );
+            return( b.Cost < this->Cost );
           }
       };
 
@@ -58,6 +63,9 @@ namespace fpa
       itkTypeMacro( Dijkstra, TFilter );
 
     public:
+      TMST* GetMinimumSpanningTree( );
+      const TMST* GetMinimumSpanningTree( ) const;
+
       const TIntensityFunctor* GetIntensityFunctor( ) const;
       const TVertexFunctor* GetVertexFunctor( ) const;
 
@@ -77,6 +85,7 @@ namespace fpa
     protected:
       typename TIntensityFunctor::Pointer m_IntensityFunctor;
       typename TVertexFunctor::Pointer m_VertexFunctor;
+      unsigned long m_MSTIndex;
     };
 
   } // ecapseman
index 9e8e5290566bcfff0a2e6a8ea9acabef1bc8193c..c5563c92d29702a7738a7e9ebd46483b1da95156 100644 (file)
@@ -7,30 +7,58 @@
 #define __fpa__Base__Dijkstra__hxx__
 
 // -------------------------------------------------------------------------
-template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
+typename
+fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
+TMST* fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
+GetMinimumSpanningTree( )
+{
+  return(
+    dynamic_cast< TMST* >(
+      this->itk::ProcessObject::GetOutput( this->m_MSTIndex )
+      )
+    );
+}
+
+// -------------------------------------------------------------------------
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
+const typename
+fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
+TMST* fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
+GetMinimumSpanningTree( ) const
+{
+  return(
+    dynamic_cast< const TMST* >(
+      this->itk::ProcessObject::GetOutput( this->m_MSTIndex )
+      )
+    );
+}
+
+// -------------------------------------------------------------------------
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
 const typename
-fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >::
+fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
 TIntensityFunctor*
-fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >::
+fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
 GetIntensityFunctor( ) const
 {
   return( this->m_IntensityFunctor );
 }
 
 // -------------------------------------------------------------------------
-template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
 const typename
-fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >::
+fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
 TVertexFunctor*
-fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >::
+fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
 GetVertexFunctor( ) const
 {
   return( this->m_VertexFunctor );
 }
 
 // -------------------------------------------------------------------------
-template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
-void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >::
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
+void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
 SetFunctor( TIntensityFunctor* functor )
 {
   if( this->m_IntensityFunctor.GetPointer( ) != functor )
@@ -42,8 +70,8 @@ SetFunctor( TIntensityFunctor* functor )
 }
 
 // -------------------------------------------------------------------------
-template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
-void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >::
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
+void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
 SetFunctor( TVertexFunctor* functor )
 {
   if( this->m_VertexFunctor.GetPointer( ) != functor )
@@ -55,27 +83,96 @@ SetFunctor( TVertexFunctor* functor )
 }
 
 // -------------------------------------------------------------------------
-template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
-fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >::
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
+fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
 Dijkstra( )
   : Superclass( ),
     _TMarksInterface( this ),
     _TSeedsInterface( this )
 {
+  this->m_MSTIndex = this->GetNumberOfRequiredOutputs( );
+  this->SetNumberOfRequiredOutputs( this->m_MSTIndex + 1 );
+  this->itk::ProcessObject::SetNthOutput( this->m_MSTIndex, TMST::New( ) );
 }
 
 // -------------------------------------------------------------------------
-template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
-fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >::
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
+fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
 ~Dijkstra( )
 {
 }
 
 // -------------------------------------------------------------------------
-template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
-void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >::
+template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST >
+void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >::
 GenerateData( )
 {
+  // Init objects
+  this->_ConfigureOutputs( TOutputValue( 0 ) );
+  this->_InitMarks( this->GetNumberOfSeeds( ) );
+  TMST* mst = this->GetMinimumSpanningTree( );
+
+  // Init queue
+  std::vector< _TNode > q;
+  unsigned long frontId = 1;
+  for( TVertex seed: this->GetSeeds( ) )
+    q.push_back( _TNode( seed, seed, frontId++ ) );
+
+  // Main loop
+  while( q.size( ) > 0 )
+  {
+    // Get next candidate
+    std::pop_heap( q.begin( ), q.end( ) );
+    _TNode node = q.back( );
+    q.pop_back( );
+    if( this->_IsMarked( node.Vertex ) )
+      continue;
+    this->_Mark( node.Vertex, node.FrontId );
+
+    // Ok, pixel lays inside region
+    this->_SetOutputValue( node.Vertex, node.Cost );
+    mst->SetParent( node.Vertex, node.Parent );
+
+    // Add neighborhood
+    TVertices neighbors = this->_GetNeighbors( node.Vertex );
+    for( TVertex neigh: neighbors )
+    {
+      if( this->_IsMarked( neigh ) )
+      {
+        // Invoke stop at collisions
+        unsigned long nColl = this->_Collisions( node.Vertex, neigh );
+        if(
+          this->StopAtOneFront( ) &&
+          this->GetNumberOfSeeds( ) > 1 &&
+          nColl == 1
+          )
+          q.clear( );
+      }
+      else
+      {
+        // Compute new cost
+        TOutputValue ncost =
+          this->m_VertexFunctor->Evaluate( neigh, node.Vertex );
+        if( this->m_IntensityFunctor.IsNotNull( ) )
+          ncost = this->m_IntensityFunctor->Evaluate( ncost );
+
+        // This algorithm only supports positive values
+        if( ncost >= TOutputValue( 0 ) )
+        {
+          // Insert new node
+          _TNode nn( neigh, node.Vertex, node.FrontId );
+          nn.Cost = node.Cost + ncost;
+          q.push_back( nn );
+          std::push_heap( q.begin( ), q.end( ) );
+
+        } // fi
+
+      } // fi
+
+    } // rof
+
+  } // elihw
+  this->_FreeMarks( );
 }
 
 #endif // __fpa__Base__Dijkstra__hxx__
diff --git a/libs/fpa/Base/Functors/InvertValue.h b/libs/fpa/Base/Functors/InvertValue.h
new file mode 100644 (file)
index 0000000..9da3221
--- /dev/null
@@ -0,0 +1,80 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__Functors__InvertValue__h__
+#define __fpa__Base__Functors__InvertValue__h__
+
+#include <cmath>
+#include <itkFunctionBase.h>
+
+namespace fpa
+{
+  namespace Base
+  {
+    namespace Functors
+    {
+      /**
+       */
+      template< class _TInputValue, class _TOutputValue >
+      class InvertValue
+        : public itk::FunctionBase< _TInputValue, _TOutputValue >
+      {
+      public:
+        typedef InvertValue                                      Self;
+        typedef itk::FunctionBase< _TInputValue, _TOutputValue > Superclass;
+        typedef itk::SmartPointer< Self >                        Pointer;
+        typedef itk::SmartPointer< const Self >                  ConstPointer;
+
+        typedef _TInputValue  TInputValue;
+        typedef _TOutputValue TOutputValue;
+
+      public:
+        itkNewMacro( Self );
+        itkTypeMacro( InvertValue, itk::FunctionBase );
+
+        itkGetConstMacro( Alpha, double );
+        itkGetConstMacro( Beta, double );
+        itkSetMacro( Alpha, double );
+        itkSetMacro( Beta, double );
+
+      public:
+        virtual TOutputValue Evaluate( const TInputValue& a ) const override
+          {
+            double d = this->m_Alpha;
+            d += std::pow( double( a ), this->m_Beta );
+            double x = -1;
+            if( std::fabs( d ) > double( 0 ) )
+              x =
+                double( 1 ) /
+                ( this->m_Alpha + std::pow( double( a ), this->m_Beta ) );
+            return( TOutputValue( x ) );
+          }
+
+      protected:
+        InvertValue( )
+          : Superclass( ),
+            m_Alpha( double( 1 ) ),
+            m_Beta( double( 1 ) )
+          { }
+        virtual ~InvertValue( ) { }
+
+      private:
+        InvertValue( const Self& other );
+        Self& operator=( const Self& other );
+
+      protected:
+        double m_Alpha;
+        double m_Beta;
+      };
+
+    } // ecapseman
+
+  } // ecapseman
+
+} // ecapseman
+
+#endif // __fpa__Base__Functors__InvertValue__h__
+
+// eof - $RCSfile$
diff --git a/libs/fpa/Base/Functors/VertexParentBase.h b/libs/fpa/Base/Functors/VertexParentBase.h
new file mode 100644 (file)
index 0000000..2725063
--- /dev/null
@@ -0,0 +1,58 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__Functors__VertexParentBase__h__
+#define __fpa__Base__Functors__VertexParentBase__h__
+
+#include <itkObject.h>
+#include <itkObjectFactory.h>
+
+namespace fpa
+{
+  namespace Base
+  {
+    namespace Functors
+    {
+      /**
+       */
+      template< class _TVertex, class _TOutputValue >
+      class VertexParentBase
+        : public itk::Object
+      {
+      public:
+        typedef VertexParentBase                Self;
+        typedef itk::Object                     Superclass;
+        typedef itk::SmartPointer< Self >       Pointer;
+        typedef itk::SmartPointer< const Self > ConstPointer;
+
+        typedef _TVertex      TVertex;
+        typedef _TOutputValue TOutputValue;
+
+      public:
+        itkTypeMacro( VertexParentBase, TFilter );
+
+      public:
+        virtual TOutputValue Evaluate( const TVertex& v, const TVertex& p ) const = 0;
+
+      protected:
+        VertexParentBase( )
+          : Superclass( )
+          { }
+        virtual ~VertexParentBase( ) { }
+
+      private:
+        VertexParentBase( const Self& other );
+        Self& operator=( const Self& other );
+      };
+
+    } // ecapseman
+
+  } // ecapseman
+
+} // ecapseman
+
+#endif // __fpa__Base__Functors__VertexParentBase__h__
+
+// eof - $RCSfile$
diff --git a/libs/fpa/Base/MinimumSpanningTree.h b/libs/fpa/Base/MinimumSpanningTree.h
new file mode 100644 (file)
index 0000000..35cb6e4
--- /dev/null
@@ -0,0 +1,77 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__MinimumSpanningTree__h__
+#define __fpa__Base__MinimumSpanningTree__h__
+
+#include <vector>
+
+namespace fpa
+{
+  namespace Base
+  {
+    /**
+     */
+    template< class _TVertex, class _Superclass >
+    class MinimumSpanningTree
+      : public _Superclass
+    {
+    public:
+      typedef MinimumSpanningTree             Self;
+      typedef _Superclass                     Superclass;
+      typedef itk::SmartPointer< Self >       Pointer;
+      typedef itk::SmartPointer< const Self > ConstPointer;
+
+      typedef _TVertex                      TVertex;
+      typedef std::pair< TVertex, bool >    TCollision;
+      typedef std::vector< TCollision >     TCollisionsRow;
+      typedef std::vector< TCollisionsRow > TCollisions;
+      typedef std::vector< TVertex >        TVertices;
+
+    protected:
+      typedef std::vector< unsigned long > _TRow;
+      typedef std::vector< _TRow >         _TMatrix;
+
+    public:
+      itkTypeMacro( fpa::Base::MinimumSpanningTree, _Superclass );
+
+    public:
+      const TCollisions& GetCollisions( ) const;
+      void SetCollisions( const TCollisions& collisions );
+
+      void ClearSeeds( );
+      void AddSeed( const TVertex& seed );
+
+      virtual TVertex GetParent( const TVertex& v ) const = 0;
+      virtual void SetParent( const TVertex& v, const TVertex& p ) = 0;
+
+      virtual TVertices GetPath( const TVertex& a ) const;
+      virtual TVertices GetPath( const TVertex& a, const TVertex& b ) const;
+
+    protected:
+      MinimumSpanningTree( );
+      virtual ~MinimumSpanningTree( );
+
+    private:
+      MinimumSpanningTree( const Self& other );
+      Self& operator=( const Self& other );
+
+    protected:
+      TCollisions m_Collisions;
+      _TMatrix    m_FrontPaths;
+      std::vector< TVertex > m_Seeds;
+    };
+
+  } // ecapseman
+
+} // ecapseman
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#  include <fpa/Base/MinimumSpanningTree.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __fpa__Base__MinimumSpanningTree__h__
+
+// eof - $RCSfile$
diff --git a/libs/fpa/Base/MinimumSpanningTree.hxx b/libs/fpa/Base/MinimumSpanningTree.hxx
new file mode 100644 (file)
index 0000000..79ddebe
--- /dev/null
@@ -0,0 +1,239 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__MinimumSpanningTree__hxx__
+#define __fpa__Base__MinimumSpanningTree__hxx__
+
+// -------------------------------------------------------------------------
+template< class _TVertex, class _Superclass >
+const typename fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+TCollisions& fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+GetCollisions( ) const
+{
+  return( this->m_Collisions );
+}
+
+// -------------------------------------------------------------------------
+template< class _TVertex, class _Superclass >
+void fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+SetCollisions( const TCollisions& collisions )
+{
+  static const unsigned long _inf =
+    std::numeric_limits< unsigned long >::max( );
+  if( this->m_Collisions == collisions )
+    return;
+
+  this->m_Collisions = collisions;
+
+  // Prepare a front graph
+  unsigned long N = this->m_Collisions.size( );
+  _TMatrix dist( N, _TRow( N, _inf ) );
+  this->m_FrontPaths = dist;
+  for( unsigned long i = 0; i < N; ++i )
+  {
+    for( unsigned long j = 0; j < N; ++j )
+    {
+      if( this->m_Collisions[ i ][ j ].second )
+      {
+        dist[ i ][ j ] = 1;
+        dist[ j ][ i ] = 1;
+        this->m_FrontPaths[ i ][ j ] = j;
+        this->m_FrontPaths[ j ][ i ] = i;
+
+      } // fi
+
+    } // rof
+    dist[ i ][ i ] = 0;
+    this->m_FrontPaths[ i ][ i ] = i;
+
+  } // rof
+
+  // Use Floyd-Warshall to compute all possible paths between fronts
+  for( unsigned long k = 0; k < N; ++k )
+  {
+    for( unsigned long i = 0; i < N; ++i )
+    {
+      for( unsigned long j = 0; j < N; ++j )
+      {
+        // WARNING: you don't want a numeric overflow!!!
+        unsigned long dik = dist[ i ][ k ];
+        unsigned long dkj = dist[ k ][ j ];
+        unsigned long sum = _inf;
+        if( dik < _inf && dkj < _inf )
+          sum = dik + dkj;
+
+        // Ok, continue Floyd-Warshall
+        if( sum < dist[ i ][ j ] )
+        {
+          dist[ i ][ j ] = sum;
+          this->m_FrontPaths[ i ][ j ] = this->m_FrontPaths[ i ][ k ];
+
+        } // fi
+
+      } // rof
+
+    } // rof
+
+  } // rof
+  this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TVertex, class _Superclass >
+void fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+ClearSeeds( )
+{
+  this->m_Seeds.clear( );
+  if( this->m_DataObject != NULL )
+    this->m_DataObject->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TVertex, class _Superclass >
+void fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+AddSeed( const _TVertex& seed )
+{
+  this->m_Seeds.push_back( seed );
+  if( this->m_DataObject != NULL )
+    this->m_DataObject->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TVertex, class _Superclass >
+typename fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+TVertices fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+GetPath( const _TVertex& a ) const
+{
+  TVertices vertices;
+  _TVertex it = a;
+  _TVertex p = this->GetParent( it );
+  while( it != p )
+  {
+    vertices.push_back( it );
+    it = p;
+    p = this->GetParent( it );
+
+  } // elihw
+  vertices.push_back( it );
+  return( vertices );
+}
+
+// -------------------------------------------------------------------------
+template< class _TVertex, class _Superclass >
+typename fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+TVertices fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+GetPath( const _TVertex& a, const _TVertex& b ) const
+{
+  static const unsigned long _inf =
+    std::numeric_limits< unsigned long >::max( );
+
+  TVertices vertices;
+  TVertices pa = this->GetPath( a );
+  TVertices pb = this->GetPath( b );
+  if( pa.size( ) > 0 && pb.size( ) > 0 )
+  {
+    // Find front identifiers
+    unsigned long ia = _inf, ib = _inf;
+    unsigned long N = this->m_Seeds.size( );
+    for( unsigned long i = 0; i < N; ++i )
+    {
+      if( this->m_Seeds[ i ] == pa[ pa.size( ) - 1 ] )
+        ia = i;
+      if( this->m_Seeds[ i ] == pb[ pb.size( ) - 1 ] )
+        ib = i;
+
+    } // rof
+
+    // Check if there is a front-jump between given seeds
+    if( ia != ib )
+    {
+      // Compute front path
+      std::vector< long > fpath;
+      fpath.push_back( ia );
+      while( ia != ib )
+      {
+        ia = this->m_FrontPaths[ ia ][ ib ];
+        fpath.push_back( ia );
+
+      } // elihw
+
+      // Continue only if both fronts are connected
+      unsigned int N = fpath.size( );
+      if( N > 0 )
+      {
+        // First path: from start vertex to first collision
+        vertices = this->GetPath(
+          a, this->m_Collisions[ fpath[ 0 ] ][ fpath[ 1 ] ].first
+          );
+
+        // Intermediary paths
+        for( unsigned int i = 1; i < N - 1; ++i )
+        {
+          TVertices ipath =
+            this->GetPath(
+              this->m_Collisions[ fpath[ i ] ][ fpath[ i - 1 ] ].first,
+              this->m_Collisions[ fpath[ i ] ][ fpath[ i + 1 ] ].first
+              );
+          for( long id = 0; id < ipath.size( ); ++id )
+            vertices.push_back( ipath[ id ] );
+
+        } // rof
+
+        // Final path: from last collision to end point
+        TVertices lpath =
+          this->GetPath(
+            this->m_Collisions[ fpath[ N - 1 ] ][ fpath[ N - 2 ] ].first, b
+            );
+        for( long id = 0; id < lpath.size( ); ++id )
+          vertices.push_back( lpath[ id ] );
+
+      } // fi
+    }
+    else
+    {
+      // Ignore common part: find common ancestor
+      long aIt = pa.size( ) - 1;
+      long bIt = pb.size( ) - 1;
+      bool cont = true;
+      while( aIt >= 0 && bIt >= 0 && cont )
+      {
+        cont = ( pa[ aIt ] == pb[ bIt ] );
+        aIt--;
+        bIt--;
+
+      } // elihw
+      aIt++;
+      bIt++;
+
+      // Glue both parts
+      for( long cIt = 0; cIt <= aIt; ++cIt )
+        vertices.push_back( pa[ cIt ] );
+      for( ; bIt >= 0; --bIt )
+        vertices.push_back( pb[ bIt ] );
+
+    } // fi
+
+  } // fi
+  return( vertices );
+}
+
+// -------------------------------------------------------------------------
+template< class _TVertex, class _Superclass >
+fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+MinimumSpanningTree( )
+  : Superclass( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class _TVertex, class _Superclass >
+fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >::
+~MinimumSpanningTree( )
+{
+}
+
+#endif // __fpa__Base__MinimumSpanningTree__hxx__
+
+// eof - $RCSfile$
diff --git a/libs/fpa/Image/Dijkstra.h b/libs/fpa/Image/Dijkstra.h
new file mode 100644 (file)
index 0000000..8516d7b
--- /dev/null
@@ -0,0 +1,91 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Image__Dijkstra__h__
+#define __fpa__Image__Dijkstra__h__
+
+#include <fpa/Base/Dijkstra.h>
+#include <fpa/Base/SeedsInterface.h>
+#include <fpa/Image/MarksInterface.h>
+#include <fpa/Image/Filter.h>
+#include <fpa/Image/MinimumSpanningTree.h>
+#include <fpa/Image/Functors/VertexParentBase.h>
+
+namespace fpa
+{
+  namespace Image
+  {
+    /**
+     */
+    template< class _TInputImage, class _TOutputImage >
+    class Dijkstra
+      : public fpa::Base::Dijkstra< fpa::Image::Filter< _TInputImage, _TOutputImage >, fpa::Image::MarksInterface< _TInputImage::ImageDimension >, fpa::Base::SeedsInterface< typename _TInputImage::IndexType, typename _TInputImage::IndexType::LexicographicCompare >, fpa::Image::MinimumSpanningTree< _TInputImage::ImageDimension > >
+    {
+    public:
+      // Interfaces
+      typedef fpa::Image::Filter< _TInputImage, _TOutputImage > TFilter;
+      typedef fpa::Image::MarksInterface< _TInputImage::ImageDimension > TMarksInterface;
+      typedef fpa::Base::SeedsInterface< typename _TInputImage::IndexType, typename _TInputImage::IndexType::LexicographicCompare > TSeedsInterface;
+      typedef fpa::Image::MinimumSpanningTree< _TInputImage::ImageDimension > TMST;
+
+      // Smart pointers
+      typedef Dijkstra Self;
+      typedef fpa::Base::Dijkstra< TFilter, TMarksInterface, TSeedsInterface, TMST > Superclass;
+      typedef itk::SmartPointer< Self > Pointer;
+      typedef itk::SmartPointer< const Self > ConstPointer;
+
+      typedef typename TFilter::TInputImage TInputImage;
+      typedef typename TFilter::TOutputValue TOutputValue;
+      typedef typename TFilter::TVertex TVertex;
+
+    public:
+      itkNewMacro( Self );
+      itkTypeMacro( fpa::Image::Dijkstra, fpa::Base::Dijkstra );
+
+    protected:
+      Dijkstra( ) : Superclass( ) { }
+      virtual ~Dijkstra( )        { }
+
+      virtual void _ConfigureOutputs( const TOutputValue& init_value ) override
+        {
+          this->Superclass::_ConfigureOutputs( init_value );
+
+          typename TVertex::OffsetType o;
+          o.Fill( 0 );
+          const TInputImage* input = this->GetInput( );
+          TMST* mst = this->GetMinimumSpanningTree( );
+          mst->CopyInformation( input );
+          mst->SetBufferedRegion( input->GetRequestedRegion( ) );
+          mst->Allocate( );
+          mst->FillBuffer( o );
+        }
+
+      virtual void GenerateData( ) override
+        {
+          // Configure functors with input image
+          typedef typename TFilter::TInputImage  _TInputIage;
+          typedef typename TFilter::TOutputValue _TOutputValue;
+          typedef fpa::Image::Functors::VertexParentBase< _TInputImage, _TOutputValue > _TVFunc;
+          _TVFunc* vfunc =
+            dynamic_cast< _TVFunc* >( this->m_VertexFunctor.GetPointer( ) );
+          if( vfunc != NULL )
+            vfunc->SetImage( this->GetInput( ) );
+
+          // Ok, continue
+          this->Superclass::GenerateData( );
+        }
+
+    private:
+      Dijkstra( const Self& other );
+      Self& operator=( const Self& other );
+    };
+
+  } // ecapseman
+
+} // ecapseman
+
+#endif // __fpa__Image__Dijkstra__h__
+
+// eof - $RCSfile$
diff --git a/libs/fpa/Image/Functors/VertexIdentity.h b/libs/fpa/Image/Functors/VertexIdentity.h
new file mode 100644 (file)
index 0000000..716c2e7
--- /dev/null
@@ -0,0 +1,67 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Image__Functors__VertexIdentity__h__
+#define __fpa__Image__Functors__VertexIdentity__h__
+
+#include <fpa/Image/Functors/VertexParentBase.h>
+
+namespace fpa
+{
+  namespace Image
+  {
+    namespace Functors
+    {
+      /**
+       */
+      template< class _TInputImage, class _TOutputValue >
+      class VertexIdentity
+        : public fpa::Image::Functors::VertexParentBase< _TInputImage, _TOutputValue >
+      {
+      public:
+        typedef _TInputImage                    TInputImage;
+        typedef _TOutputValue                   TOutputValue;
+        typedef VertexIdentity                  Self;
+        typedef itk::SmartPointer< Self >       Pointer;
+        typedef itk::SmartPointer< const Self > ConstPointer;
+        typedef fpa::Image::Functors::VertexParentBase< TInputImage, TOutputValue > Superclass;
+
+        typedef typename Superclass::TVertex TVertex;
+
+      public:
+        itkNewMacro( Self );
+        itkTypeMacro(
+          fpa::Image::Functors::VertexIdentity,
+          fpa::Image::Functors::VertexParentBase
+          );
+
+      public:
+        virtual TOutputValue Evaluate(
+          const TVertex& a, const TVertex& p
+          ) const override
+          {
+            return( TOutputValue( this->m_Image->GetPixel( a ) ) );
+          }
+
+      protected:
+        VertexIdentity( )
+          : Superclass( )
+          { }
+        virtual ~VertexIdentity( ) { }
+
+      private:
+        VertexIdentity( const Self& other );
+        Self& operator=( const Self& other );
+      };
+
+    } // ecapseman
+
+  } // ecapseman
+
+} // ecapseman
+
+#endif // __fpa__Image__Functors__VertexIdentity__h__
+
+// eof - $RCSfile$
diff --git a/libs/fpa/Image/Functors/VertexParentBase.h b/libs/fpa/Image/Functors/VertexParentBase.h
new file mode 100644 (file)
index 0000000..20b76e0
--- /dev/null
@@ -0,0 +1,63 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Image__Functors__VertexParentBase__h__
+#define __fpa__Image__Functors__VertexParentBase__h__
+
+#include <fpa/Base/Functors/VertexParentBase.h>
+
+namespace fpa
+{
+  namespace Image
+  {
+    namespace Functors
+    {
+      /**
+       */
+      template< class _TInputImage, class _TOutputValue >
+      class VertexParentBase
+        : public fpa::Base::Functors::VertexParentBase< typename _TInputImage::IndexType, _TOutputValue >
+      {
+      public:
+        typedef _TInputImage                    TInputImage;
+        typedef _TOutputValue                   TOutputValue;
+        typedef typename TInputImage::IndexType TVertex;
+        typedef VertexParentBase                Self;
+        typedef itk::SmartPointer< Self >       Pointer;
+        typedef itk::SmartPointer< const Self > ConstPointer;
+        typedef fpa::Base::Functors::VertexParentBase< TVertex, TOutputValue > Superclass;
+
+      public:
+        itkTypeMacro(
+          fpa::Image::Functors::VertexParentBase,
+          fpa::Base::Functors::VertexParentBase
+          );
+
+        itkGetConstObjectMacro( Image, TInputImage );
+        itkSetConstObjectMacro( Image, TInputImage );
+
+      protected:
+        VertexParentBase( )
+          : Superclass( )
+          { }
+        virtual ~VertexParentBase( ) { }
+
+      private:
+        VertexParentBase( const Self& other );
+        Self& operator=( const Self& other );
+
+      protected:
+        typename TInputImage::ConstPointer m_Image;
+      };
+
+    } // ecapseman
+
+  } // ecapseman
+
+} // ecapseman
+
+#endif // __fpa__Image__Functors__VertexParentBase__h__
+
+// eof - $RCSfile$
diff --git a/libs/fpa/Image/MinimumSpanningTree.h b/libs/fpa/Image/MinimumSpanningTree.h
new file mode 100644 (file)
index 0000000..a93157c
--- /dev/null
@@ -0,0 +1,71 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Image__MinimumSpanningTree__h__
+#define __fpa__Image__MinimumSpanningTree__h__
+
+#include <fpa/Base/MinimumSpanningTree.h>
+#include <itkImage.h>
+
+namespace fpa
+{
+  namespace Image
+  {
+    /**
+     */
+    template< unsigned int _VDim >
+    class MinimumSpanningTree
+      : public fpa::Base::MinimumSpanningTree< itk::Index< _VDim >, itk::Image< itk::Offset< _VDim >, _VDim > >
+    {
+    public:
+      typedef itk::Index< _VDim > TVertex;
+      typedef itk::Image< itk::Offset< _VDim >, _VDim > TBaseImage;
+
+      typedef MinimumSpanningTree             Self;
+      typedef itk::SmartPointer< Self >       Pointer;
+      typedef itk::SmartPointer< const Self > ConstPointer;
+      typedef fpa::Base::MinimumSpanningTree< TVertex, TBaseImage > Superclass;
+
+      typedef typename Superclass::TCollision     TCollision;
+      typedef typename Superclass::TCollisionsRow TCollisionsRow;
+      typedef typename Superclass::TCollisions    TCollisions;
+      typedef typename Superclass::TVertices      TVertices;
+
+    public:
+      itkNewMacro( Self );
+      itkTypeMacro(
+        fpa::Image::MinimumSpanningTree,
+        fpa::Base::MinimumSpanningTree
+        );
+
+    public:
+      virtual TVertex GetParent( const TVertex& v ) const override
+        {
+          return( v + this->GetPixel( v ) );
+        }
+      virtual void SetParent( const TVertex& v, const TVertex& p ) override
+        {
+          this->SetPixel( v, p - v );
+        }
+
+    protected:
+      MinimumSpanningTree( )
+        : Superclass( )
+        { }
+      virtual ~MinimumSpanningTree( )
+        { }
+
+    private:
+      MinimumSpanningTree( const Self& other );
+      Self& operator=( const Self& other );
+    };
+
+  } // ecapseman
+
+} // ecapseman
+
+#endif // __fpa__Image__MinimumSpanningTree__h__
+
+// eof - $RCSfile$