...

[FrontAlgorithms.git] / lib / fpa / Image / SkeletonFilter.hxx

diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx
index 0f33388f80d04b8202fa3791a787beb66531e4ba..6a322371c3a5a4d0b6b4d447ce422417de12a554 100644 (file)
--- a/lib/fpa/Image/SkeletonFilter.hxx
+++ b/lib/fpa/Image/SkeletonFilter.hxx
@@ -167,13 +167,86 @@ GenerateData( )
   dijkstra->Update( );
 
   // Compute end-points
+  const _TMST* mst = dijkstra->GetMinimumSpanningTree( );
   std::vector< TIndex > end_points;
   this->_EndPoints(
-    end_points,
-    this->m_DistanceMap->GetOutput( ),
-    dijkstra->GetMinimumSpanningTree( ),
-    dijkstra->GetSkeletonQueue( )
+    end_points, this->m_DistanceMap->GetOutput( ),
+    mst, dijkstra->GetSkeletonQueue( )
     );
+
+  // Compute symbolic branches
+  typedef std::map< TIndex, TIndex, typename TIndex::LexicographicCompare > _TTags;
+  _TTags tags, branches;
+  typename std::vector< TIndex >::const_iterator eIt = end_points.begin( );
+  for( ; eIt != end_points.end( ); ++eIt )
+  {
+    // Tag path
+    TIndex it = *eIt;
+    TIndex p = mst->GetParent( it );
+    typename _TTags::iterator bIt = tags.end( );
+    while( it != p && bIt == tags.end( ) )
+    {
+      typename _TTags::iterator tIt = tags.find( it );
+      if( tIt != tags.end( ) )
+      {
+        // Ok, a bifurcation point has been found
+        // branch1: tIt->second <-> it (ok)
+        branches[ tIt->second ] = it;
+
+        // branch2: *eit <-> it (ok)
+        branches[ *eIt ] = it;
+
+        // branch3: it <-> until next bifurcation
+        bIt = tIt;
+      }
+      else
+        tags[ it ] = *eIt;
+      it = p;
+      p = mst->GetParent( it );
+
+    } // elihw
+    if( bIt != tags.end( ) )
+    {
+      TIndex pTag = bIt->second;
+      TIndex nTag = bIt->first;
+      it = bIt->first;
+      p = it;
+      while( tags[ it ] == pTag )
+      {
+        tags[ it ] = nTag;
+        p = it;
+        it = mst->GetParent( it );
+
+      } // elihw
+      tags[ it ] = nTag;
+      branches[ bIt->first ] = p;
+    }
+    else
+    {
+      tags[ it ] = *eIt;
+      branches[ *eIt ] = it;
+
+    } // fi
+
+  } // rof
+
+  // Fill full branches
+  typedef typename _TMST::TVertices _TVertices;
+  typedef typename TSkeleton::TPath _TPath;
+
+  TSkeleton* sk = this->GetOutput( );
+  typename _TTags::const_iterator bIt = branches.begin( );
+  for( ; bIt != branches.end( ); ++bIt )
+  {
+    _TVertices v = mst->GetPath( bIt->first, bIt->second );
+    typename _TPath::Pointer path = _TPath::New( );
+    path->SetReferenceImage( this->GetInput( ) );
+    typename _TVertices::const_reverse_iterator vIt = v.rbegin( );
+    for( ; vIt != v.rend( ); ++vIt )
+      path->AddVertex( *vIt );
+    sk->AddBranch( path );
+
+  } // rof
 }
 
 // -------------------------------------------------------------------------
@@ -276,84 +349,6 @@ _EndPoints(
   } // rof
 }
 
-
-
-/* TODO
-   this->m_DistanceMap = TDistanceMap::New( );
-   this->m_DistanceMap->InsideIsPositiveOn( );
-   this->m_DistanceMap->SquaredDistanceOff( );
-   this->m_DistanceMap->UseImageSpacingOn( );
-
-   template< class _TImage, class _TScalar >
-   void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
-   _Skeleton( const std::vector< TVertex >& end_points, _TAdjacencies& A )
-   {
-   typedef typename TSkeleton::TPath _TPath;
-   typedef itk::Image< unsigned long, TImage::ImageDimension > _TTagsImage;
-
-   TMST* mst = this->GetMinimumSpanningTree( );
-   TSkeleton* skeleton = this->GetSkeleton( );
-
-   // Tag branches
-   typename _TTagsImage::Pointer tags = _TTagsImage::New( );
-   tags->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
-   tags->SetRequestedRegion( mst->GetRequestedRegion( ) );
-   tags->SetBufferedRegion( mst->GetBufferedRegion( ) );
-   tags->Allocate( );
-   tags->FillBuffer( 0 );
-   typename std::vector< TVertex >::const_iterator eIt = end_points.begin( );
-   for( ; eIt != end_points.end( ); ++eIt )
-   {
-   TVertex it = *eIt;
-   TVertex p = mst->GetParent( it );
-   while( it != p )
-   {
-   tags->SetPixel( it, tags->GetPixel( it ) + 1 );
-   it = p;
-   p = mst->GetParent( it );
-   
-   } // elihw
-   tags->SetPixel( it, tags->GetPixel( it ) + 1 );
-
-   } // rof
-
-   // Build paths (branches)
-   eIt = end_points.begin( );
-   for( ; eIt != end_points.end( ); ++eIt )
-   {
-   TVertex it = *eIt;
-   TVertex p = mst->GetParent( it );
-   TVertex sIdx = it;
-   typename _TPath::Pointer path = _TPath::New( );
-   path->SetReferenceImage( mst );
-   while( it != p )
-   {
-   if( tags->GetPixel( sIdx ) != tags->GetPixel( it ) )
-   {
-   // Ok, a new branch has been added
-   path->AddVertex( it );
-   skeleton->AddBranch( path );
-   
-   // Update a new starting index
-   path = _TPath::New( );
-   path->SetReferenceImage( mst );
-   sIdx = it;
-   }
-   else
-   path->AddVertex( it );
-   it = p;
-   p = mst->GetParent( it );
-
-   } // elihw
-   
-   // Finally, add last branch
-   path->AddVertex( it );
-   skeleton->AddBranch( path );
-   
-   } // rof
-   }
-*/
-
 #endif // __fpa__Image__SkeletonFilter__hxx__
 
 // eof - $RCSfile$