libs/fpa/Image/EndPointsFilter.hxx

   1 #ifndef __fpa__Image__EndPointsFilter__hxx__
   2 #define __fpa__Image__EndPointsFilter__hxx__
   3
   4 #include <functional>
   5 #include <map>
   6 #include <queue>
   7 #include <itkImageRegionConstIteratorWithIndex.h>
   8
   9 // -------------------------------------------------------------------------
  10 template< class _TDistanceMap, class _TCostMap >
  11 const typename fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  12 TIndices& fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  13 GetBifurcations( ) const
  14 {
  15   return( this->m_Bifurcations );
  16 }
  17
  18 // -------------------------------------------------------------------------
  19 template< class _TDistanceMap, class _TCostMap >
  20 const typename fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  21 TIndices& fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  22 GetEndPoints( ) const
  23 {
  24   return( this->m_EndPoints );
  25 }
  26
  27 // -------------------------------------------------------------------------
  28 template< class _TDistanceMap, class _TCostMap >
  29 void fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  30 Compute( )
  31 {
  32   typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt;
  33   typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt;
  34   typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
  35   typedef typename _TQueue::value_type _TQueueValue;
  36   typedef itk::Image< unsigned char, _TCostMap::ImageDimension > _TMarks;
  37
  38   // Some values
  39   typename _TDistanceMap::RegionType region =
  40     this->m_DistanceMap->GetRequestedRegion( );
  41   typename _TMarks::Pointer marks = _TMarks::New( );
  42   marks->SetLargestPossibleRegion( this->m_DistanceMap->GetLargestPossibleRegion( ) );
  43   marks->SetRequestedRegion( this->m_DistanceMap->GetRequestedRegion( ) );
  44   marks->SetBufferedRegion( this->m_DistanceMap->GetBufferedRegion( ) );
  45   marks->SetSpacing( this->m_DistanceMap->GetSpacing( ) );
  46   marks->SetOrigin( this->m_DistanceMap->GetOrigin( ) );
  47   marks->SetDirection( this->m_DistanceMap->GetDirection( ) );
  48   marks->Allocate( );
  49   marks->FillBuffer( 0 );
  50
  51   // Create queue
  52   _TQueue queue;
  53   _TDistMapIt dIt( this->m_DistanceMap, region );
  54   _TCostMapIt cIt( this->m_CostMap, region );
  55   dIt.GoToBegin( );
  56   cIt.GoToBegin( );
  57   for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
  58   {
  59     double d = double( dIt.Get( ) );
  60     if( d > 0 )
  61     {
  62       double v = double( cIt.Get( ) ) / ( d * d );
  63       queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
  64
  65     } // fi
  66
  67   } // rof
  68
  69   // BFS from maximum queue
  70   while( queue.size( ) > 0 )
  71   {
  72     // Get node
  73     auto nIt = queue.begin( );
  74     auto n = *nIt;
  75     queue.erase( nIt );
  76
  77     unsigned char m = marks->GetPixel( n.second );
  78     if( ( m & 0x01 ) == 0x01 )
  79       continue;
  80
  81     // Mark it
  82     m |= 0x01;
  83     marks->SetPixel( n.second, m );
  84     this->m_EndPoints.insert( n.second );
  85
  86     // Get path
  87     typename TMST::TPath::Pointer path;
  88     this->m_MST->GetPath( path, n.second );
  89     for( unsigned long i = 0; i < path->GetSize( ); ++i )
  90     {
  91       TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i );
  92       typename _TCostMap::PointType cnt;
  93       this->m_CostMap->TransformIndexToPhysicalPoint( idx, cnt );
  94       double d = double( this->m_DistanceMap->GetPixel( idx ) );
  95       d *= double( 2 );
  96
  97       // Mark sphere
  98       std::queue< TIndex > q;
  99       q.push( idx );
 100       while( q.size( ) > 0 )
 101       {
 102         TIndex v = q.front( );
 103         q.pop( );
 104         unsigned char m = marks->GetPixel( v );
 105         if( ( m & 0x02 ) == 0x02 )
 106           continue;
 107         m |= 0x03;
 108         marks->SetPixel( v, m );
 109
 110         for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x )
 111         {
 112           for( int y = -1; y <= 1; y += 2 )
 113           {
 114             TIndex w = v;
 115             w[ x ] += y;
 116             if( region.IsInside( w ) )
 117             {
 118               typename _TCostMap::PointType p;
 119               this->m_CostMap->TransformIndexToPhysicalPoint( w, p );
 120               if( cnt.EuclideanDistanceTo( p ) <= d )
 121                 q.push( w );
 122
 123             } // fi
 124
 125           } // rof
 126
 127         } // rof
 128
 129       } // elihw
 130
 131     } // rof
 132
 133   } // elihw
 134 }
 135
 136 // -------------------------------------------------------------------------
 137 template< class _TDistanceMap, class _TCostMap >
 138 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
 139 EndPointsFilter( )
 140   : Superclass( )
 141 {
 142 }
 143
 144 // -------------------------------------------------------------------------
 145 template< class _TDistanceMap, class _TCostMap >
 146 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
 147 ~EndPointsFilter( )
 148 {
 149 }
 150
 151 #endif // __fpa__Image__EndPointsFilter__hxx__
 152
 153 // eof - $RCSfile$