lib/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
  37   // Create queue
  38   _TQueue queue;
  39   _TDistMapIt dIt(
  40     this->m_DistanceMap, this->m_DistanceMap->GetRequestedRegion( )
  41     );
  42   _TCostMapIt cIt(
  43     this->m_CostMap, this->m_CostMap->GetRequestedRegion( )
  44     );
  45   dIt.GoToBegin( );
  46   cIt.GoToBegin( );
  47   for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
  48   {
  49     double d = double( dIt.Get( ) );
  50     if( d > 0 )
  51     {
  52       double v = double( cIt.Get( ) ) * d;
  53       queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
  54
  55     } // fi
  56
  57   } // rof
  58
  59   TIndices marks;
  60   while( queue.size( ) > 0 )
  61   {
  62     auto nIt = queue.begin( );
  63     auto n = *nIt;
  64     queue.erase( nIt );
  65
  66     if( marks.find( n.second ) == marks.end( ) )
  67     {
  68       std::cout << queue.size( ) << " " << n.first << std::endl;
  69       marks.insert( n.second );
  70       this->m_EndPoints.insert( n.second );
  71       auto path = this->m_MST->GetPath( n.second );
  72       std::cout << path.size( ) << std::endl;
  73       for( auto pIt = path.begin( ); pIt != path.end( ); ++pIt )
  74       {
  75         double d = double( this->m_DistanceMap->GetPixel( *pIt ) );
  76         d = std::sqrt( std::fabs( d ) );
  77         typename _TCostMap::PointType center;
  78         this->m_CostMap->TransformIndexToPhysicalPoint( *pIt, center );
  79
  80         std::queue< TIndex > q;
  81         TIndices m;
  82         q.push( *pIt );
  83         while( q.size( ) > 0 )
  84         {
  85           TIndex idx = q.front( );
  86           q.pop( );
  87           if( m.find( idx ) != m.end( ) )
  88             continue;
  89           m.insert( idx );
  90           marks.insert( idx );
  91           for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x )
  92           {
  93             for( int y = -1; y <= 1; y += 2 )
  94             {
  95               TIndex idx2 = idx;
  96               idx2[ x ] += y;
  97               typename _TCostMap::PointType c;
  98               this->m_CostMap->TransformIndexToPhysicalPoint( idx2, c );
  99               if( this->m_CostMap->GetRequestedRegion( ).IsInside( idx2 ) )
 100               {
 101                 if( center.EuclideanDistanceTo( c ) <= ( d * 1.5 ) )
 102                   q.push( idx2 );
 103
 104               } // fi
 105
 106             } // rof
 107
 108           } // rof
 109
 110         } // elihw
 111
 112       } // rof
 113
 114     } // fi
 115
 116   } // elihw
 117 }
 118
 119 // -------------------------------------------------------------------------
 120 template< class _TDistanceMap, class _TCostMap >
 121 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
 122 EndPointsFilter( )
 123   : Superclass( )
 124 {
 125 }
 126
 127 // -------------------------------------------------------------------------
 128 template< class _TDistanceMap, class _TCostMap >
 129 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
 130 ~EndPointsFilter( )
 131 {
 132 }
 133
 134 #endif // __fpa__Image__EndPointsFilter__hxx__
 135
 136 // eof - $RCSfile$