X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FImage%2FSkeletonFilter.hxx;h=9cb3019c85aed4ae34be9dd76d6b6cf80a080aef;hb=e8548d1529b25ab74bca646333612fbe3e16e73f;hp=2b1cb6591e06751e9a5224b3f292b95441ee6803;hpb=cd6b3f8433deaf2ba7c6bb0ece8e5912c760db17;p=FrontAlgorithms.git

diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx
index 2b1cb65..9cb3019 100644
--- a/lib/fpa/Image/SkeletonFilter.hxx
+++ b/lib/fpa/Image/SkeletonFilter.hxx
@@ -6,6 +6,7 @@
 #ifndef __fpa__Image__SkeletonFilter__hxx__
 #define __fpa__Image__SkeletonFilter__hxx__
 
+#include <itkImageRegionIteratorWithIndex.h>
 #include <itkMinimumMaximumImageCalculator.h>
 
 // -------------------------------------------------------------------------
@@ -138,20 +139,18 @@ template< class _TImage, class _TScalar >
 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 _SetOutputValue( const TVertex& vertex, const TOutputValue& value )
 {
-  /* TODO
-     typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
-
-     this->Superclass::_UpdateResult( n );
-     double d = double( this->GetInput( )->GetPixel( n.Vertex ) );
-     if( d >= double( 0 ) )
-     {
-     // Update skeleton candidates
-     d += double( 1e-5 );
-     double v = double( n.Result ) / ( d * d );
-     this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) );
-
-     } // fi
-  */
+  typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
+
+  this->Superclass::_SetOutputValue( vertex, value );
+  double d = double( this->m_DistanceMap->GetOutput( )->GetPixel( vertex ) );
+  if( d >= double( 0 ) )
+  {
+    // Update skeleton candidates
+    d += double( 1e-5 );
+    double v = double( value ) / ( d * d );
+    this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, vertex ) );
+
+  } // fi
 }
 
 // -------------------------------------------------------------------------
@@ -159,51 +158,51 @@ template< class _TImage, class _TScalar >
 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 _EndPoints( std::vector< TVertex >& end_points, _TAdjacencies& A )
 {
-  /* TODO
-     auto marks = this->GetMarks( );
-     auto mst = this->GetMinimumSpanningTree( );
-     auto spac = marks->GetSpacing( );
-
-     // Some values
-     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 );
-
-     // BFS from maximum queue
-     while( this->m_SkeletonQueue.size( ) > 0 )
-     {
-     // Get node
-     auto nIt = this->m_SkeletonQueue.begin( );
-     auto n = *nIt;
-     this->m_SkeletonQueue.erase( nIt );
-
-     // Mark it and update end-points
-     unsigned char m = marks->GetPixel( n.second );
-     if( m != 0 )
-     continue;
-     marks->SetPixel( n.second, 1 );
-     end_points.push_back( n.second );
-
-     // Mark path
-     TVertex it = n.second;
-     TVertex p = mst->GetParent( it );
-     while( it != p )
-     {
-     this->_MarkSphere( it );
-     it = p;
-     p = mst->GetParent( it );
-
-     } // elihw
-     this->_MarkSphere( it );
-     A[ n.second ] = it;
-
-     } // elihw
-  */
+  typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
+
+  TMarks* marks = this->GetMarks( );
+  TMST* mst = this->GetMinimumSpanningTree( );
+  typename TImage::SpacingType spac = marks->GetSpacing( );
+
+  // Some values
+  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 );
+
+  // BFS from maximum queue
+  while( this->m_SkeletonQueue.size( ) > 0 )
+  {
+    // Get node
+    typename _TSkeletonQueue::iterator nIt = this->m_SkeletonQueue.begin( );
+    _TSkeletonQueueValue n = *nIt;
+    this->m_SkeletonQueue.erase( nIt );
+
+    // Mark it and update end-points
+    unsigned char m = marks->GetPixel( n.second );
+    if( m != 0 )
+      continue;
+    marks->SetPixel( n.second, 1 );
+    end_points.push_back( n.second );
+
+    // Mark path
+    TVertex it = n.second;
+    TVertex p = mst->GetParent( it );
+    while( it != p )
+    {
+      this->_MarkSphere( it );
+      it = p;
+      p = mst->GetParent( it );
+
+    } // elihw
+    this->_MarkSphere( it );
+    A[ n.second ] = it;
+
+  } // elihw
 }
 
 // -------------------------------------------------------------------------
@@ -211,69 +210,65 @@ template< class _TImage, class _TScalar >
 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 _Skeleton( const std::vector< TVertex >& end_points, _TAdjacencies& A )
 {
-  /* 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 )
-     {
-     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)
-     for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
-     {
-     TVertex it = *eIt;
-     TVertex p = mst->GetParent( it );
-     TVertex 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
-  */
+  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 );
+  for( TVertex it: end_points )
+  {
+    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)
+  for( TVertex it: end_points )
+  {
+    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
 }
 
 // -------------------------------------------------------------------------
@@ -281,54 +276,52 @@ template< class _TImage, class _TScalar >
 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 _MarkSphere( const TVertex& idx )
 {
-  /* TODO
-     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;
-
-     TVertex 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 )
-     {
-     TVertex w = mIt.GetIndex( );
-     typename _TImage::PointType p;
-     marks->TransformIndexToPhysicalPoint( w, p );
-     mIt.Set( 1 );
-
-     } // rof
-  */
+  typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
+
+  static const double _eps = std::sqrt( double( TImage::ImageDimension + 1 ) );
+  const TScalarImage* dmap = this->m_DistanceMap->GetOutput( );
+  TMarks* marks = this->GetMarks( );
+  typename TImage::SpacingType spac = dmap->GetSpacing( );
+  typename TImage::RegionType region = dmap->GetRequestedRegion( );
+
+  typename TImage::PointType cnt;
+  dmap->TransformIndexToPhysicalPoint( idx, cnt );
+  double r = double( dmap->GetPixel( idx ) ) * _eps;
+
+  TVertex i0, i1;
+  for( unsigned int d = 0; d < TImage::ImageDimension; ++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 < TImage::ImageDimension; ++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 )
+  {
+    TVertex w = mIt.GetIndex( );
+    typename TImage::PointType p;
+    marks->TransformIndexToPhysicalPoint( w, p );
+    mIt.Set( 1 );
+
+  } // rof
 }
 
 #endif // __fpa__Image__SkeletonFilter__hxx__