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   // Some values
  38   typename _TDistanceMap::RegionType region =
  39     this->m_DistanceMap->GetRequestedRegion( );
  40
  41   // Create queue
  42   _TQueue queue;
  43   _TDistMapIt dIt( this->m_DistanceMap, region );
  44   _TCostMapIt cIt( this->m_CostMap, region );
  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   // BFS from maximum queue
  60   TIndices marks;
  61   while( queue.size( ) > 0 )
  62   {
  63     // Get node
  64     auto nIt = queue.begin( );
  65     auto n = *nIt;
  66     queue.erase( nIt );
  67     if( marks.find( n.second ) != marks.end( ) )
  68       continue;
  69
  70     // Mark it
  71     marks.insert( n.second );
  72     this->m_EndPoints.insert( n.second );
  73
  74     // Get path
  75     typename TMST::TPath::Pointer path;
  76     this->m_MST->GetPath( path, n.second );
  77     for( unsigned long i = 0; i < path->GetSize( ); ++i )
  78     {
  79       typename TMST::TPath::TContinuousIndex cidx = path->GetVertex( i );
  80       typename _TCostMap::PointType cnt;
  81       this->m_CostMap->TransformContinuousIndexToPhysicalPoint( cidx, cnt );
  82       TIndex idx;
  83       this->m_CostMap->TransformPhysicalPointToIndex( cnt, idx );
  84       double d = double( this->m_DistanceMap->GetPixel( idx ) );
  85
  86       /* TODO
  87          TIndex idx;
  88          for( unsigned int d = 0; d < _TCostMap::ImageDimension; ++d )
  89          idx[ d ] = cidx[ d ];
  90       */
  91
  92     } // rof
  93
  94     // TODO: temporary
  95     queue.clear( );
  96
  97   } // elihw
  98
  99   /* TODO
 100      while( queue.size( ) > 0 )
 101      {
 102      marks.insert( n.second );
 103      this->m_EndPoints.insert( n.second );
 104      auto path = this->m_MST->GetPath( n.second );
 105      std::cout << path.size( ) << std::endl;
 106      for( auto pIt = path.begin( ); pIt != path.end( ); ++pIt )
 107      {
 108      double d = double( this->m_DistanceMap->GetPixel( *pIt ) );
 109      d = std::sqrt( std::fabs( d ) );
 110      typename _TCostMap::PointType center;
 111
 112
 113      std::queue< TIndex > q;
 114      TIndices m;
 115      q.push( *pIt );
 116      while( q.size( ) > 0 )
 117      {
 118      TIndex idx = q.front( );
 119      q.pop( );
 120      if( m.find( idx ) != m.end( ) )
 121      continue;
 122      m.insert( idx );
 123      marks.insert( idx );
 124      for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x )
 125      {
 126      for( int y = -1; y <= 1; y += 2 )
 127      {
 128      TIndex idx2 = idx;
 129      idx2[ x ] += y;
 130      typename _TCostMap::PointType c;
 131      this->m_CostMap->TransformIndexToPhysicalPoint( idx2, c );
 132      if( this->m_CostMap->GetRequestedRegion( ).IsInside( idx2 ) )
 133      {
 134      if( center.EuclideanDistanceTo( c ) <= ( d * 1.5 ) )
 135      q.push( idx2 );
 136
 137      } // fi
 138
 139      } // rof
 140
 141      } // rof
 142
 143      } // elihw
 144
 145      } // rof
 146
 147      } // fi
 148
 149      } // elihw
 150   */
 151 }
 152
 153 // -------------------------------------------------------------------------
 154 template< class _TDistanceMap, class _TCostMap >
 155 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
 156 EndPointsFilter( )
 157   : Superclass( )
 158 {
 159 }
 160
 161 // -------------------------------------------------------------------------
 162 template< class _TDistanceMap, class _TCostMap >
 163 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
 164 ~EndPointsFilter( )
 165 {
 166 }
 167
 168 #endif // __fpa__Image__EndPointsFilter__hxx__
 169
 170 // eof - $RCSfile$