]> Creatis software - FrontAlgorithms.git/commitdiff
...
authorLeonardo Flórez-Valencia <florez-l@javeriana.edu.co>
Tue, 14 Feb 2017 22:38:56 +0000 (17:38 -0500)
committerLeonardo Flórez-Valencia <florez-l@javeriana.edu.co>
Tue, 14 Feb 2017 22:38:56 +0000 (17:38 -0500)
17 files changed:
data/test_skeleton_2D.png [new file with mode: 0644]
examples/CMakeLists.txt
examples/Skeleton_00.cxx [new file with mode: 0644]
lib/fpa/Base/Algorithm.h
lib/fpa/Base/Dijkstra.h
lib/fpa/Base/Functors/GaussianModel.h
lib/fpa/Base/Functors/Inverse.h
lib/fpa/Base/RegionGrow.h
lib/fpa/Config.h
lib/fpa/Image/Algorithm.h
lib/fpa/Image/Dijkstra.h
lib/fpa/Image/Functors/SimpleDijkstraCost.h
lib/fpa/Image/Functors/SimpleNeighborhood.h
lib/fpa/Image/MinimumSpanningTree.h
lib/fpa/Image/RegionGrow.h
lib/fpa/Image/SkeletonFilter.h
lib/fpa/Image/SkeletonFilter.hxx

diff --git a/data/test_skeleton_2D.png b/data/test_skeleton_2D.png
new file mode 100644 (file)
index 0000000..17feca4
Binary files /dev/null and b/data/test_skeleton_2D.png differ
index 4a8ae35dc7fdca38df8861a05c558dffbc5df8d3..acb1bc5fc9d3be3a63d6e0a7d4c91b3708d34595 100644 (file)
@@ -5,6 +5,7 @@ IF(BUILD_Examples)
     _examples
     RegionGrow_00
     MoriRegionGrow_00
+    Skeleton_00
     )
   INCLUDE_DIRECTORIES(
     ${PROJECT_SOURCE_DIR}/lib
@@ -12,7 +13,7 @@ IF(BUILD_Examples)
     )
   FOREACH(_example ${_examples})
     ADD_EXECUTABLE(fpa_example_${_example} ${_example}.cxx)
-    TARGET_LINK_LIBRARIES(fpa_example_${_example} ${ITK_LIBRARIES})
+    TARGET_LINK_LIBRARIES(fpa_example_${_example} ${ITK_LIBRARIES} ${VTK_LIBRARIES})
   ENDFOREACH(_example)
 ENDIF(BUILD_Examples)
 
diff --git a/examples/Skeleton_00.cxx b/examples/Skeleton_00.cxx
new file mode 100644 (file)
index 0000000..5ddb5e4
--- /dev/null
@@ -0,0 +1,49 @@
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkSignedMaurerDistanceMapImageFilter.h>
+#include <itkMinimumMaximumImageCalculator.h>
+#include <fpa/Image/SkeletonFilter.h>
+
+// -------------------------------------------------------------------------
+typedef itk::Image< double, 2 >                                   TImage;
+typedef itk::ImageFileReader< TImage >                            TReader;
+typedef itk::SignedMaurerDistanceMapImageFilter< TImage, TImage > TDMap;
+typedef itk::MinimumMaximumImageCalculator< TImage >              TMinMax;
+typedef fpa::Image::SkeletonFilter< TImage >                      TFilter;
+
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+  if( argc < 2 )
+  {
+    std::cerr
+      << "Usage: " << argv[ 0 ]
+      << " input_filename" << std::endl;
+    return( 1 );
+
+  } // fi
+  std::string in_fname = argv[ 1 ];
+
+  TReader::Pointer reader = TReader::New( );
+  reader->SetFileName( in_fname );
+
+  TDMap::Pointer dmap = TDMap::New( );
+  dmap->SetInput( reader->GetOutput( ) );
+  dmap->InsideIsPositiveOn( );
+  dmap->SquaredDistanceOff( );
+  dmap->UseImageSpacingOn( );
+  dmap->Update( );
+
+  TMinMax::Pointer minmax = TMinMax::New( );
+  minmax->SetImage( dmap->GetOutput( ) );
+  minmax->Compute( );
+
+  TFilter::Pointer filter = TFilter::New( );
+  filter->SetInput( dmap->GetOutput( ) );
+  filter->AddSeed( minmax->GetIndexOfMaximum( ), 0 );
+  filter->Update( );
+
+  return( 0 );
+}
+
+// eof - $RCSfile$
index 297727ac22ef0a5019d4b71f60e793b49c846ae4..da7bee89127774915bcf49311806763e35fe2ba0 100644 (file)
@@ -75,7 +75,7 @@ namespace fpa
       Algorithm( );
       virtual ~Algorithm( );
 
-      virtual void GenerateData( ) fpa_OVERRIDE;
+      virtual void GenerateData( ) override;
 
       // Particular methods
       virtual bool _ContinueGenerateData( );
index 5633f221f5b3bddedda29ab40575c6bcba12f397..aa62f90bbe2dcd07eae54fd684c306d277f95c57 100644 (file)
@@ -56,16 +56,16 @@ namespace fpa
       Dijkstra( );
       virtual ~Dijkstra( );
 
-      virtual void _AfterGenerateData( ) fpa_OVERRIDE;
+      virtual void _AfterGenerateData( ) override;
 
-      virtual void _UpdateResult( const _TQueueNode& n ) fpa_OVERRIDE;
+      virtual void _UpdateResult( const _TQueueNode& n ) override;
       virtual bool _UpdateValue(
         _TQueueNode& v, const _TQueueNode& p
-        ) fpa_OVERRIDE;
-      virtual unsigned long _QueueSize( ) const fpa_OVERRIDE;
-      virtual void _QueueClear( ) fpa_OVERRIDE;
-      virtual void _QueuePush( const _TQueueNode& node ) fpa_OVERRIDE;
-      virtual _TQueueNode _QueuePop( ) fpa_OVERRIDE;
+        ) override;
+      virtual unsigned long _QueueSize( ) const override;
+      virtual void _QueueClear( ) override;
+      virtual void _QueuePush( const _TQueueNode& node ) override;
+      virtual _TQueueNode _QueuePop( ) override;
 
     private:
       // Purposely not defined
index a74e5176bde7d1a0699d762cff0c4b1f90e08792..8d32f51de050d0fc4ee26f4a32dbc2f4aedd9e7b 100644 (file)
@@ -45,7 +45,7 @@ namespace fpa
         */
 
       public:
-        virtual TOutput Evaluate( const TInput& x ) const fpa_OVERRIDE;
+        virtual TOutput Evaluate( const TInput& x ) const override;
 
       protected:
         GaussianModel( );
index dc8c7a24fd43845586395c9fb20e2a4784e2d2ad..af68022246189edf414cfa37b26c906a09671370 100644 (file)
@@ -33,7 +33,7 @@ namespace fpa
         itkSetMacro( NegativeValue, _TOutput );
 
       public:
-        virtual TOutput Evaluate( const TInput& x ) const fpa_OVERRIDE;
+        virtual TOutput Evaluate( const TInput& x ) const override;
 
       protected:
         Inverse( );
index 2ef070399aa6416c5ed310b92657c4d26ea90ba9..a8fcc932aa5dbba92883bb983b872698ad52f4c6 100644 (file)
@@ -47,11 +47,11 @@ namespace fpa
 
       virtual bool _UpdateValue(
         _TQueueNode& v, const _TQueueNode& p
-        ) fpa_OVERRIDE;
-      virtual unsigned long _QueueSize( ) const fpa_OVERRIDE;
-      virtual void _QueueClear( ) fpa_OVERRIDE;
-      virtual void _QueuePush( const _TQueueNode& node ) fpa_OVERRIDE;
-      virtual _TQueueNode _QueuePop( ) fpa_OVERRIDE;
+        ) override;
+      virtual unsigned long _QueueSize( ) const override;
+      virtual void _QueueClear( ) override;
+      virtual void _QueuePush( const _TQueueNode& node ) override;
+      virtual _TQueueNode _QueuePop( ) override;
 
     private:
       // Purposely not defined
index c59b1beb4d82b890c452b5cc880b1cd1448746fc..9728f0bb4fa1723e788e4000206daf7340180fde 100644 (file)
@@ -7,20 +7,6 @@
  * =========================================================================
  */
 
-#if __cplusplus >= 201103L
-#  define fpa_OVERRIDE override
-#  define fpa_DELETE_FUNCTION =delete
-#  define fpa_NULLPTR  nullptr
-#  define fpa_NOEXCEPT noexcept
-#  define fpa_HAS_CXX11_STATIC_ASSERT
-#  define fpa_HAS_CXX11_RVREF
-#else
-#  define fpa_OVERRIDE
-#  define fpa_DELETE_FUNCTION
-#  define fpa_NULLPTR  NULL
-#  define fpa_NOEXCEPT throw()
-#endif
-
 #endif // __fpa__Config__h__
 
 // eof - $RCSfile$
index 70a87a611543fe78b7d3076a44df96eba271355a..d613ce3d8afce33b090f781c9abc04df658cc9d5 100644 (file)
@@ -40,13 +40,13 @@ namespace fpa
       Algorithm( );
       virtual ~Algorithm( );
 
-      virtual void _BeforeGenerateData( ) fpa_OVERRIDE;
-      virtual void _InitMarks( ) fpa_OVERRIDE;
-      virtual void _InitResults( const TOutput& init_value ) fpa_OVERRIDE;
-      virtual bool _IsMarked( const TVertex& v ) const fpa_OVERRIDE;
-      virtual void _Mark( const _TQueueNode& n ) fpa_OVERRIDE;
-      virtual TFrontId _GetMark( const TVertex& v ) const fpa_OVERRIDE;
-      virtual void _UpdateResult( const _TQueueNode& n ) fpa_OVERRIDE;
+      virtual void _BeforeGenerateData( ) override;
+      virtual void _InitMarks( ) override;
+      virtual void _InitResults( const TOutput& init_value ) override;
+      virtual bool _IsMarked( const TVertex& v ) const override;
+      virtual void _Mark( const _TQueueNode& n ) override;
+      virtual TFrontId _GetMark( const TVertex& v ) const override;
+      virtual void _UpdateResult( const _TQueueNode& n ) override;
 
     private:
       // Purposely not defined
index 253f97380b312e4e252ef5944fcc2f2d348084ba..96c642b273be7fce59a3d6f51bde113b3ed2f827 100644 (file)
@@ -37,7 +37,7 @@ namespace fpa
       Dijkstra( );
       virtual ~Dijkstra( );
 
-      virtual void _BeforeGenerateData( ) fpa_OVERRIDE;
+      virtual void _BeforeGenerateData( ) override;
 
     private:
       // Purposely not defined
index d9a0405b5d2d6ff30747bb16114a7d9c51185ff9..bad95001c3a80f4e1b62395d1508ca24b6cb888a 100644 (file)
@@ -38,7 +38,7 @@ namespace fpa
       public:
         virtual TOutput Evaluate(
           const TIndex& a, const TIndex& b
-          ) const fpa_OVERRIDE;
+          ) const override;
 
       protected:
         SimpleDijkstraCost( );
index 1fd854757003c21efe03231c5d230429c2785b1d..9f777b01f19326584a64de06f81bb32aa0028f2c 100644 (file)
@@ -36,7 +36,7 @@ namespace fpa
         itkSetMacro( Order, unsigned int );
 
       public:
-        virtual TOutput Evaluate( const TIndex& center ) const fpa_OVERRIDE;
+        virtual TOutput Evaluate( const TIndex& center ) const override;
 
       protected:
         SimpleNeighborhood( );
index 4a4e0e62bbbc53ebbc4cd5830ca6ec1d005c144e..7d3941d23bbe57d904389af0161cffb74002f20a 100644 (file)
@@ -32,15 +32,15 @@ namespace fpa
         );
 
     public:
-      virtual TVertex GetParent( const TVertex& v ) const fpa_OVERRIDE;
-      virtual void SetParent( const TVertex& v, const TVertex& p ) fpa_OVERRIDE;
+      virtual TVertex GetParent( const TVertex& v ) const override;
+      virtual void SetParent( const TVertex& v, const TVertex& p ) override;
 
       virtual void GetPath(
         typename TPath::Pointer& path, const TVertex& a
-        ) const fpa_OVERRIDE;
+        ) const override;
       virtual void GetPath(
         typename TPath::Pointer& path, const TVertex& a, const TVertex& b
-        ) const fpa_OVERRIDE;
+        ) const override;
 
     protected:
       MinimumSpanningTree( );
index ad41e80a2dbea689573b407e36c3ea2c23b03b22..dbc57f00759be5d017c1c49bb0c42f0eba583c71 100644 (file)
@@ -34,7 +34,7 @@ namespace fpa
       RegionGrow( );
       virtual ~RegionGrow( );
 
-      virtual void _BeforeGenerateData( ) fpa_OVERRIDE;
+      virtual void _BeforeGenerateData( ) override;
 
     private:
       // Purposely not defined
index 5a2563444877ccdbf09cddf5aedef50b3a509db2..0fbb3d5ade737d233e150e986aa67bc792f16a73 100644 (file)
@@ -1,6 +1,7 @@
 #ifndef __fpa__Image__SkeletonFilter__h__
 #define __fpa__Image__SkeletonFilter__h__
 
+#include <map>
 #include <vector>
 #include <fpa/Image/Dijkstra.h>
 #include <cpExtensions/DataStructures/Skeleton.h>
@@ -35,6 +36,7 @@ namespace fpa
     protected:
       typedef typename Superclass::_TQueueNode _TQueueNode;
       typedef std::multimap< double, TIndex, std::greater< double > > _TSkeletonQueue;
+      typedef std::map< TIndex, TIndex, typename TIndex::LexicographicCompare > _TAdjacencies;
 
     public:
       itkNewMacro( Self );
@@ -48,12 +50,13 @@ namespace fpa
       SkeletonFilter( );
       virtual ~SkeletonFilter( );
 
-      virtual void _BeforeGenerateData( ) fpa_OVERRIDE;
-      virtual void _UpdateResult( const _TQueueNode& n ) fpa_OVERRIDE;
-      virtual void _AfterGenerateData( ) fpa_OVERRIDE;
+      virtual void _BeforeGenerateData( ) override;
+      virtual void _UpdateResult( const _TQueueNode& n ) override;
+      virtual void _AfterGenerateData( ) override;
 
-      void _EndPoints( std::vector< TIndex >& end_points );
-      void _Skeleton( const std::vector< TIndex >& end_points );
+      void _EndPoints( std::vector< TIndex >& end_points, _TAdjacencies& A );
+      void _Skeleton( const std::vector< TIndex >& end_points, _TAdjacencies& A );
+      void _MarkSphere( const TIndex& idx );
 
     private:
       // Purposely not defined
index 069eea071af885b3ae432a8edebfa66d327e9ec1..d8dc66a1f0841216b02c7ddeeeba8ee941a5ddac 100644 (file)
@@ -93,36 +93,32 @@ _AfterGenerateData( )
 {
   this->Superclass::_AfterGenerateData( );
 
+  _TAdjacencies A;
   std::vector< TIndex > end_points;
-  this->_EndPoints( end_points );
-  this->_Skeleton( end_points );
+  this->_EndPoints( end_points, A );
+  this->_Skeleton( end_points, A );
 }
 
 // -------------------------------------------------------------------------
 template< class _TImage >
 void fpa::Image::SkeletonFilter< _TImage >::
-_EndPoints( std::vector< TIndex >& end_points )
+_EndPoints( std::vector< TIndex >& end_points, _TAdjacencies& A )
 {
-  typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
-
-  static const double _eps = std::sqrt( double( Self::Dimension + 1 ) );
-  auto input = this->GetInput( );
-  auto mst = this->GetMinimumSpanningTree( );
   auto marks = this->GetMarks( );
+  auto mst = this->GetMinimumSpanningTree( );
+  auto spac = marks->GetSpacing( );
 
   // Some values
-  marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
-  marks->SetRequestedRegion( input->GetRequestedRegion( ) );
-  marks->SetBufferedRegion( input->GetBufferedRegion( ) );
-  marks->SetSpacing( input->GetSpacing( ) );
-  marks->SetOrigin( input->GetOrigin( ) );
-  marks->SetDirection( input->GetDirection( ) );
+  marks->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
+  marks->SetRequestedRegion( mst->GetRequestedRegion( ) );
+  marks->SetBufferedRegion( mst->GetBufferedRegion( ) );
+  marks->SetSpacing( mst->GetSpacing( ) );
+  marks->SetOrigin( mst->GetOrigin( ) );
+  marks->SetDirection( mst->GetDirection( ) );
   marks->Allocate( );
   marks->FillBuffer( 0 );
-  auto spac = marks->GetSpacing( );
 
   // BFS from maximum queue
-  auto region = input->GetRequestedRegion( );
   while( this->m_SkeletonQueue.size( ) > 0 )
   {
     // Get node
@@ -136,66 +132,20 @@ _EndPoints( std::vector< TIndex >& end_points )
       continue;
     marks->SetPixel( n.second, 1 );
     end_points.push_back( n.second );
+    std::cout << this->m_SkeletonQueue.size( ) << std::endl;
 
     // Mark path
     TIndex it = n.second;
     TIndex p = mst->GetParent( it );
     while( it != p )
     {
-      // TODO: tags->SetPixel( it, tags->GetPixel( it ) + 1 );
+      this->_MarkSphere( it );
       it = p;
       p = mst->GetParent( it );
 
     } // elihw
-    // TODO: tags->SetPixel( it, tags->GetPixel( it ) + 1 );
-
-    /* TODO: use mst directly rather than computing paths
-       typename TMST::TPath::Pointer path;
-       mst->GetPath( path, n.second );
-       for( unsigned long i = 0; i < path->GetSize( ); ++i )
-       {
-       TIndex idx = path->GetVertex( i );
-       typename _TImage::PointType cnt;
-       input->TransformIndexToPhysicalPoint( idx, cnt );
-       double r = double( input->GetPixel( idx ) ) * _eps;
-
-       TIndex i0, i1;
-       for( unsigned int d = 0; d < Self::Dimension; ++d )
-       {
-       long off = long( std::ceil( r / double( spac[ d ] ) ) );
-       if( off < 3 )
-       off = 3;
-       i0[ d ] = idx[ d ] - off;
-       i1[ d ] = idx[ d ] + off;
-
-       if( i0[ d ] < region.GetIndex( )[ d ] )
-       i0[ d ] = region.GetIndex( )[ d ];
-
-       if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
-       i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
-
-       } // rof
-
-       typename _TImage::SizeType size;
-       for( unsigned int d = 0; d < Self::Dimension; ++d )
-       size[ d ] = i1[ d ] - i0[ d ] + 1;
-
-       typename _TImage::RegionType neighRegion;
-       neighRegion.SetIndex( i0 );
-       neighRegion.SetSize( size );
-
-       _TMarksIt mIt( marks, neighRegion );
-       for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
-       {
-       TIndex w = mIt.GetIndex( );
-       typename _TImage::PointType p;
-       marks->TransformIndexToPhysicalPoint( w, p );
-       mIt.Set( 1 );
-
-       } // rof
-
-       } // rof
-    */
+    this->_MarkSphere( it );
+    A[ n.second ] = it;
 
   } // elihw
 }
@@ -203,67 +153,124 @@ _EndPoints( std::vector< TIndex >& end_points )
 // -------------------------------------------------------------------------
 template< class _TImage >
 void fpa::Image::SkeletonFilter< _TImage >::
-_Skeleton( const std::vector< TIndex >& end_points )
+_Skeleton( const std::vector< TIndex >& end_points, _TAdjacencies& A )
 {
-  typedef itk::Image< unsigned long, Self::Dimension > _TTagsImage;
-  typedef typename TMST::TPath _TPath;
+  std::cout << end_points.size( ) << " " << A.size( ) << std::endl;
+
+  /* TODO
+     typedef itk::Image< unsigned long, Self::Dimension > _TTagsImage;
+     typedef typename TMST::TPath _TPath;
+
+     auto mst = this->GetMinimumSpanningTree( );
+     auto 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 );
+     for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
+     {
+     TIndex it = *eIt;
+     TIndex 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)
+     for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
+     {
+     TIndex it = *eIt;
+     TIndex p = mst->GetParent( it );
+     TIndex sIdx = *eIt;
+     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
+  */
+}
 
-  auto mst = this->GetMinimumSpanningTree( );
-  auto 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 );
-  for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
+// -------------------------------------------------------------------------
+template< class _TImage >
+void fpa::Image::SkeletonFilter< _TImage >::
+_MarkSphere( const TIndex& idx )
+{
+  typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
+
+  static const double _eps = std::sqrt( double( Self::Dimension + 1 ) );
+  auto input = this->GetInput( );
+  auto marks = this->GetMarks( );
+  auto spac = input->GetSpacing( );
+  auto region = input->GetRequestedRegion( );
+
+  typename _TImage::PointType cnt;
+  input->TransformIndexToPhysicalPoint( idx, cnt );
+  double r = double( input->GetPixel( idx ) ) * _eps;
+
+  TIndex i0, i1;
+  for( unsigned int d = 0; d < Self::Dimension; ++d )
   {
-    TIndex it = *eIt;
-    TIndex p = mst->GetParent( it );
-    while( it != p )
-    {
-      tags->SetPixel( it, tags->GetPixel( it ) + 1 );
-      it = p;
-      p = mst->GetParent( it );
+    long off = long( std::ceil( r / double( spac[ d ] ) ) );
+    if( off < 3 )
+      off = 3;
+    i0[ d ] = idx[ d ] - off;
+    i1[ d ] = idx[ d ] + off;
 
-    } // elihw
-    tags->SetPixel( it, tags->GetPixel( it ) + 1 );
+    if( i0[ d ] < region.GetIndex( )[ d ] )
+      i0[ d ] = region.GetIndex( )[ d ];
+
+    if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
+      i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
 
   } // rof
 
-  // Build paths (branches)
-  for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
-  {
-    TIndex it = *eIt;
-    TIndex p = mst->GetParent( it );
-    TIndex sIdx = *eIt;
-    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 );
+  typename _TImage::SizeType size;
+  for( unsigned int d = 0; d < Self::Dimension; ++d )
+    size[ d ] = i1[ d ] - i0[ d ] + 1;
 
-    } // elihw
+  typename _TImage::RegionType neighRegion;
+  neighRegion.SetIndex( i0 );
+  neighRegion.SetSize( size );
 
-    // Finally, add last branch
-    path->AddVertex( it );
-    skeleton->AddBranch( path );
+  _TMarksIt mIt( marks, neighRegion );
+  for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
+  {
+    TIndex w = mIt.GetIndex( );
+    typename _TImage::PointType p;
+    marks->TransformIndexToPhysicalPoint( w, p );
+    mIt.Set( 1 );
 
   } // rof
 }