...

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

diff --git a/lib/fpa/Image/EndPointsFilter.hxx b/lib/fpa/Image/EndPointsFilter.hxx
index e2be58ffb1b3d89fe2da7e47c6c5807f43f29fbf..c2e1d1c3d2b8b8108bc2711d5c090a15e4dece16 100644 (file)
--- a/lib/fpa/Image/EndPointsFilter.hxx
+++ b/lib/fpa/Image/EndPointsFilter.hxx
@@ -29,93 +29,108 @@ template< class _TDistanceMap, class _TCostMap >
 void fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
 Compute( )
 {
-  /* TODO
-     typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt;
-     typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt;
-     typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
-     typedef typename _TQueue::value_type _TQueueValue;
-
-     // Create queue
-     _TQueue queue;
-     _TDistMapIt dIt(
-     this->m_DistanceMap, this->m_DistanceMap->GetRequestedRegion( )
-     );
-     _TCostMapIt cIt(
-     this->m_CostMap, this->m_CostMap->GetRequestedRegion( )
-     );
-     dIt.GoToBegin( );
-     cIt.GoToBegin( );
-     for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
-     {
-     double d = double( dIt.Get( ) );
-     if( d > 0 )
-     {
-     double v = double( cIt.Get( ) ) * d;
-     queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
-
-     } // fi
-
-     } // rof
-
-     TIndices marks;
-     while( queue.size( ) > 0 )
-     {
-     auto nIt = queue.begin( );
-     auto n = *nIt;
-     queue.erase( nIt );
-
-     if( marks.find( n.second ) == marks.end( ) )
-     {
-     std::cout << queue.size( ) << " " << n.first << std::endl;
-     marks.insert( n.second );
-     this->m_EndPoints.insert( n.second );
-     auto path = this->m_MST->GetPath( n.second );
-     std::cout << path.size( ) << std::endl;
-     for( auto pIt = path.begin( ); pIt != path.end( ); ++pIt )
-     {
-     double d = double( this->m_DistanceMap->GetPixel( *pIt ) );
-     d = std::sqrt( std::fabs( d ) );
-     typename _TCostMap::PointType center;
-     this->m_CostMap->TransformIndexToPhysicalPoint( *pIt, center );
-
-     std::queue< TIndex > q;
-     TIndices m;
-     q.push( *pIt );
-     while( q.size( ) > 0 )
-     {
-     TIndex idx = q.front( );
-     q.pop( );
-     if( m.find( idx ) != m.end( ) )
-     continue;
-     m.insert( idx );
-     marks.insert( idx );
-     for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x )
-     {
-     for( int y = -1; y <= 1; y += 2 )
-     {
-     TIndex idx2 = idx;
-     idx2[ x ] += y;
-     typename _TCostMap::PointType c;
-     this->m_CostMap->TransformIndexToPhysicalPoint( idx2, c );
-     if( this->m_CostMap->GetRequestedRegion( ).IsInside( idx2 ) )
-     {
-     if( center.EuclideanDistanceTo( c ) <= ( d * 1.5 ) )
-     q.push( idx2 );
-
-     } // fi
-
-     } // rof
-
-     } // rof
-
-     } // elihw
-
-     } // rof
-
-     } // fi
-
-     } // elihw
-  */
+  typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt;
+  typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt;
+  typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
+  typedef typename _TQueue::value_type _TQueueValue;
+  typedef itk::Image< unsigned char, _TCostMap::ImageDimension > _TMarks;
+
+  // Some values
+  typename _TDistanceMap::RegionType region =
+    this->m_DistanceMap->GetRequestedRegion( );
+  typename _TMarks::Pointer marks = _TMarks::New( );
+  marks->SetLargestPossibleRegion( this->m_DistanceMap->GetLargestPossibleRegion( ) );
+  marks->SetRequestedRegion( this->m_DistanceMap->GetRequestedRegion( ) );
+  marks->SetBufferedRegion( this->m_DistanceMap->GetBufferedRegion( ) );
+  marks->SetSpacing( this->m_DistanceMap->GetSpacing( ) );
+  marks->SetOrigin( this->m_DistanceMap->GetOrigin( ) );
+  marks->SetDirection( this->m_DistanceMap->GetDirection( ) );
+  marks->Allocate( );
+  marks->FillBuffer( 0 );
+
+  // Create queue
+  _TQueue queue;
+  _TDistMapIt dIt( this->m_DistanceMap, region );
+  _TCostMapIt cIt( this->m_CostMap, region );
+  dIt.GoToBegin( );
+  cIt.GoToBegin( );
+  for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
+  {
+    double d = double( dIt.Get( ) );
+    if( d > 0 )
+    {
+      double v = double( cIt.Get( ) ) / ( d * d );
+      queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
+
+    } // fi
+
+  } // rof
+
+  // BFS from maximum queue
+  while( queue.size( ) > 0 )
+  {
+    // Get node
+    auto nIt = queue.begin( );
+    auto n = *nIt;
+    queue.erase( nIt );
+
+    unsigned char m = marks->GetPixel( n.second );
+    if( ( m & 0x01 ) == 0x01 )
+      continue;
+
+    // Mark it
+    m |= 0x01;
+    marks->SetPixel( n.second, m );
+    this->m_EndPoints.insert( n.second );
+
+    // Get path
+    typename TMST::TPath::Pointer path;
+    this->m_MST->GetPath( path, n.second );
+    for( unsigned long i = 0; i < path->GetSize( ); ++i )
+    {
+      TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i );
+      typename _TCostMap::PointType cnt;
+      this->m_CostMap->TransformIndexToPhysicalPoint( idx, cnt );
+      double d = double( this->m_DistanceMap->GetPixel( idx ) );
+      d *= double( 2 );
+
+      // Mark sphere
+      std::queue< TIndex > q;
+      q.push( idx );
+      while( q.size( ) > 0 )
+      {
+        TIndex v = q.front( );
+        q.pop( );
+        unsigned char m = marks->GetPixel( v );
+        if( ( m & 0x02 ) == 0x02 )
+          continue;
+        m |= 0x03;
+        marks->SetPixel( v, m );
+
+        for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x )
+        {
+          for( int y = -1; y <= 1; y += 2 )
+          {
+            TIndex w = v;
+            w[ x ] += y;
+            if( region.IsInside( w ) )
+            {
+              typename _TCostMap::PointType p;
+              this->m_CostMap->TransformIndexToPhysicalPoint( w, p );
+              if( cnt.EuclideanDistanceTo( p ) <= d )
+                q.push( w );
+
+            } // fi
+
+          } // rof
+
+        } // rof
+
+      } // elihw
+
+    } // rof
+
+  } // elihw
 }
 
 // -------------------------------------------------------------------------