lib/fpa/Image/EndPointsFilter.hxx

   1 #ifndef __fpa__Image__EndPointsFilter__hxx__
   2 #define __fpa__Image__EndPointsFilter__hxx__
   3
   4 #include <map>
   5 #include <queue>
   6 #include <itkImageRegionConstIteratorWithIndex.h>
   7
   8 // -------------------------------------------------------------------------
   9 template< class _TDistanceMap, class _TCostMap >
  10 const typename fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  11 TIndices& fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  12 GetBifurcations( ) const
  13 {
  14   return( this->m_Bifurcations );
  15 }
  16
  17 // -------------------------------------------------------------------------
  18 template< class _TDistanceMap, class _TCostMap >
  19 const typename fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  20 TIndices& fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  21 GetEndPoints( ) const
  22 {
  23   return( this->m_EndPoints );
  24 }
  25
  26 // -------------------------------------------------------------------------
  27 template< class _TDistanceMap, class _TCostMap >
  28 void fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
  29 Compute( )
  30 {
  31   typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt;
  32   typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt;
  33   typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
  34   typedef typename _TQueue::value_type _TQueueValue;
  35
  36   // Create queue
  37   _TQueue queue;
  38   _TDistMapIt dIt(
  39     this->m_DistanceMap, this->m_DistanceMap->GetRequestedRegion( )
  40     );
  41   _TCostMapIt cIt(
  42     this->m_CostMap, this->m_CostMap->GetRequestedRegion( )
  43     );
  44   dIt.GoToBegin( );
  45   cIt.GoToBegin( );
  46   for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
  47   {
  48     double d = double( dIt.Get( ) );
  49     if( d > 0 )
  50     {
  51       double v = double( cIt.Get( ) ) * d;
  52       queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
  53
  54     } // fi
  55
  56   } // rof
  57
  58   TIndices marks;
  59   while( queue.size( ) > 0 )
  60   {
  61     auto nIt = queue.begin( );
  62     auto n = *nIt;
  63     queue.erase( nIt );
  64
  65     if( marks.find( n.second ) == marks.end( ) )
  66     {
  67       std::cout << queue.size( ) << " " << n.first << std::endl;
  68       marks.insert( n.second );
  69       this->m_EndPoints.insert( n.second );
  70       auto path = this->m_MST->GetPath( n.second );
  71       std::cout << path.size( ) << std::endl;
  72       for( auto pIt = path.begin( ); pIt != path.end( ); ++pIt )
  73       {
  74         double d = double( this->m_DistanceMap->GetPixel( *pIt ) );
  75         d = std::sqrt( std::fabs( d ) );
  76         typename _TCostMap::PointType center;
  77         this->m_CostMap->TransformIndexToPhysicalPoint( *pIt, center );
  78
  79         std::queue< TIndex > q;
  80         TIndices m;
  81         q.push( *pIt );
  82         while( q.size( ) > 0 )
  83         {
  84           TIndex idx = q.front( );
  85           q.pop( );
  86           if( m.find( idx ) != m.end( ) )
  87             continue;
  88           m.insert( idx );
  89           marks.insert( idx );
  90           for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x )
  91           {
  92             for( int y = -1; y <= 1; y += 2 )
  93             {
  94               TIndex idx2 = idx;
  95               idx2[ x ] += y;
  96               typename _TCostMap::PointType c;
  97               this->m_CostMap->TransformIndexToPhysicalPoint( idx2, c );
  98               if( this->m_CostMap->GetRequestedRegion( ).IsInside( idx2 ) )
  99               {
 100                 if( center.EuclideanDistanceTo( c ) <= ( d * 1.5 ) )
 101                   q.push( idx2 );
 102
 103               } // fi
 104
 105             } // rof
 106
 107           } // rof
 108
 109         } // elihw
 110
 111       } // rof
 112
 113     } // fi
 114
 115   } // elihw
 116 }
 117
 118 // -------------------------------------------------------------------------
 119 template< class _TDistanceMap, class _TCostMap >
 120 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
 121 EndPointsFilter( )
 122   : Superclass( )
 123 {
 124 }
 125
 126 // -------------------------------------------------------------------------
 127 template< class _TDistanceMap, class _TCostMap >
 128 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
 129 ~EndPointsFilter( )
 130 {
 131 }
 132
 133 #endif // __fpa__Image__EndPointsFilter__hxx__
 134
 135 // eof - $RCSfile$