...

[FrontAlgorithms.git] / libs / fpa / Base / Dijkstra.hxx

diff --git a/libs/fpa/Base/Dijkstra.hxx b/libs/fpa/Base/Dijkstra.hxx
index 9e8e5290566bcfff0a2e6a8ea9acabef1bc8193c..c5563c92d29702a7738a7e9ebd46483b1da95156 100644 (file)
--- a/libs/fpa/Base/Dijkstra.hxx
+++ b/libs/fpa/Base/Dijkstra.hxx
@@ -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__