lib/fpa/Image/SkeletonFilter.hxx

   1 #ifndef __fpa__Image__SkeletonFilter__hxx__
   2 #define __fpa__Image__SkeletonFilter__hxx__
   3
   4 #include <queue>
   5 #include <itkImageRegionConstIteratorWithIndex.h>
   6 #include <itkImageRegionIteratorWithIndex.h>
   7
   8 // -------------------------------------------------------------------------
   9 #define fpa_Image_SkeletonFilter_InputMacro( i_n, i_t, i_i )            \
  10   template< class _TDistanceMap, class _TCostMap >                      \
  11   typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::     \
  12   i_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::         \
  13   Get##i_n( )                                                           \
  14   {                                                                     \
  15     return(                                                             \
  16       dynamic_cast< i_t* >( this->Superclass::GetInput( i_i ) )         \
  17       );                                                                \
  18   }                                                                     \
  19   template< class _TDistanceMap, class _TCostMap >                      \
  20   const                                                                 \
  21   typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::     \
  22   i_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::         \
  23   Get##i_n( ) const                                                     \
  24   {                                                                     \
  25     return(                                                             \
  26       dynamic_cast< const i_t* >( this->Superclass::GetInput( i_i ) )   \
  27       );                                                                \
  28   }                                                                     \
  29   template< class _TDistanceMap, class _TCostMap >                      \
  30   void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::         \
  31   Set##i_n( i_t* input )                                                \
  32   {                                                                     \
  33     this->Superclass::SetNthInput( i_i, input );                        \
  34   }
  35
  36 fpa_Image_SkeletonFilter_InputMacro( DistanceMap, TDistanceMap, 0 );
  37 fpa_Image_SkeletonFilter_InputMacro( CostMap, TCostMap, 1 );
  38 fpa_Image_SkeletonFilter_InputMacro( MinimumSpanningTree, TMST, 2 );
  39
  40 // -------------------------------------------------------------------------
  41 #define fpa_Image_SkeletonFilter_OutputMacro( o_n, o_t, o_i )           \
  42   template< class _TDistanceMap, class _TCostMap >                      \
  43   typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::     \
  44   o_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::         \
  45   Get##o_n( )                                                           \
  46   {                                                                     \
  47     return(                                                             \
  48       dynamic_cast< o_t* >( this->Superclass::GetOutput( o_i ) )        \
  49       );                                                                \
  50   }                                                                     \
  51   template< class _TDistanceMap, class _TCostMap >                      \
  52   const typename                                                        \
  53   fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::              \
  54   o_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::         \
  55   Get##o_n( ) const                                                     \
  56   {                                                                     \
  57     return(                                                             \
  58       dynamic_cast< const o_t* >( this->Superclass::GetOutput( o_i ) )  \
  59       );                                                                \
  60   }
  61
  62 fpa_Image_SkeletonFilter_OutputMacro( Skeleton, TSkeleton, 0 );
  63 fpa_Image_SkeletonFilter_OutputMacro( EndPoints, TIndices, 1 );
  64 fpa_Image_SkeletonFilter_OutputMacro( Bifurcations, TIndices, 2 );
  65 fpa_Image_SkeletonFilter_OutputMacro( Marks, TMarks, 3 );
  66
  67 // -------------------------------------------------------------------------
  68 template< class _TDistanceMap, class _TCostMap >
  69 fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
  70 SkeletonFilter( )
  71   : Superclass( )
  72 {
  73   this->SetNumberOfRequiredInputs( 3 );
  74   this->SetNumberOfRequiredOutputs( 3 );
  75
  76   typename TIndices::Pointer end_points = TIndices::New( );
  77   typename TIndices::Pointer bifurcations = TIndices::New( );
  78   typename TSkeleton::Pointer skeleton = TSkeleton::New( );
  79   typename TMarks::Pointer marks = TMarks::New( );
  80   this->SetNthOutput( 0, skeleton.GetPointer( ) );
  81   this->SetNthOutput( 1, end_points.GetPointer( ) );
  82   this->SetNthOutput( 2, bifurcations.GetPointer( ) );
  83   this->SetNthOutput( 3, marks.GetPointer( ) );
  84 }
  85
  86 // -------------------------------------------------------------------------
  87 template< class _TDistanceMap, class _TCostMap >
  88 fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
  89 ~SkeletonFilter( )
  90 {
  91 }
  92
  93 // -------------------------------------------------------------------------
  94 template< class _TDistanceMap, class _TCostMap >
  95 void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
  96 GenerateData( )
  97 {
  98   // 0. I/O objects
  99   const TDistanceMap* dmap = this->GetDistanceMap( );
 100   const TCostMap* cmap = this->GetCostMap( );
 101   const TMST* mst = this->GetMinimumSpanningTree( );
 102   TIndices* ep = this->GetEndPoints( );
 103   TIndices* bi = this->GetBifurcations( );
 104   TSkeleton* sk = this->GetSkeleton( );
 105
 106   // 1. Check input correspondance
 107   // TODO
 108
 109   // 2. Detect end-points
 110   this->_EndPoints( dmap, cmap, mst, ep->Get( ) );
 111   std::cout << "endpoints" << std::endl;
 112
 113   // 3. Build skeleton and keep track of bifurcations
 114   this->_Skeleton( dmap, cmap, mst, ep->Get( ), bi->Get( ), sk );
 115 }
 116
 117 // -------------------------------------------------------------------------
 118 template< class _TDistanceMap, class _TCostMap >
 119 void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
 120 _EndPoints(
 121   const TDistanceMap* dmap,
 122   const TCostMap* cmap,
 123   const TMST* mst,
 124   TIndicesData& end_points
 125   )
 126 {
 127   typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt;
 128   typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt;
 129   typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
 130   typedef typename _TQueue::value_type _TQueueValue;
 131   typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
 132
 133   static const double _eps = std::sqrt( double( TMarks::ImageDimension + 1 ) );
 134
 135   // Some values
 136   auto marks = this->GetMarks( );
 137   marks->SetLargestPossibleRegion( dmap->GetLargestPossibleRegion( ) );
 138   marks->SetRequestedRegion( dmap->GetRequestedRegion( ) );
 139   marks->SetBufferedRegion( dmap->GetBufferedRegion( ) );
 140   marks->SetSpacing( dmap->GetSpacing( ) );
 141   marks->SetOrigin( dmap->GetOrigin( ) );
 142   marks->SetDirection( dmap->GetDirection( ) );
 143   marks->Allocate( );
 144   marks->FillBuffer( 0 );
 145
 146   // Create queue
 147   _TQueue queue;
 148   _TDistMapIt dIt( dmap, dmap->GetRequestedRegion( ) );
 149   _TCostMapIt cIt( cmap, cmap->GetRequestedRegion( ) );
 150   dIt.GoToBegin( );
 151   cIt.GoToBegin( );
 152   for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
 153   {
 154     double d = double( dIt.Get( ) );
 155     if( d > double( 0 ) )
 156     {
 157       double v = double( cIt.Get( ) ) / ( d * d );
 158       queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
 159
 160     } // fi
 161
 162   } // rof
 163
 164   // BFS from maximum queue
 165   auto region = dmap->GetRequestedRegion( );
 166   double init_v = queue.begin( )->first;
 167   while( queue.size( ) > 0 )
 168   {
 169     // Get node
 170     auto nIt = queue.begin( );
 171     auto n = *nIt;
 172     queue.erase( nIt );
 173
 174     unsigned char m = marks->GetPixel( n.second );
 175     if( m != 0 )
 176       continue;
 177
 178     std::cout << queue.size( ) << std::endl;
 179
 180     // Mark it and update end-points
 181     m |= 0x01;
 182     marks->SetPixel( n.second, m );
 183     end_points.insert( n.second );
 184
 185     // Mark path
 186     auto spac = marks->GetSpacing( );
 187     typename TMST::TPath::Pointer path;
 188     mst->GetPath( path, n.second );
 189     for( unsigned long i = 0; i < path->GetSize( ); ++i )
 190     {
 191       TIndex idx = path->GetVertex( i );
 192       typename _TCostMap::PointType cnt;
 193       cmap->TransformIndexToPhysicalPoint( idx, cnt );
 194       double r = double( dmap->GetPixel( idx ) ) * _eps;
 195
 196       TIndex i0, i1;
 197       for( unsigned int d = 0; d < TMarks::ImageDimension; ++d )
 198       {
 199         long off = long( std::ceil( r / double( spac[ d ] ) ) );
 200         if( off == 0 )
 201           off = 1;
 202         i0[ d ] = idx[ d ] - off;
 203         i1[ d ] = idx[ d ] + off;
 204
 205         if( i0[ d ] < region.GetIndex( )[ d ] )
 206           i0[ d ] = region.GetIndex( )[ d ];
 207
 208         if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
 209           i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
 210
 211       } // rof
 212
 213       typename TMarks::SizeType size;
 214       for( unsigned int d = 0; d < TMarks::ImageDimension; ++d )
 215         size[ d ] = i1[ d ] - i0[ d ] + 1;
 216
 217       typename TMarks::RegionType neighRegion;
 218       neighRegion.SetIndex( i0 );
 219       neighRegion.SetSize( size );
 220
 221       _TMarksIt mIt( marks, neighRegion );
 222       for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 223       {
 224         TIndex w = mIt.GetIndex( );
 225         typename _TCostMap::PointType p;
 226         marks->TransformIndexToPhysicalPoint( w, p );
 227         if( cnt.EuclideanDistanceTo( p ) <= r )
 228           mIt.Set( mIt.Get( ) | 0x02 );
 229
 230       } // rof
 231
 232     } // rof
 233
 234   } // elihw
 235 }
 236
 237 // -------------------------------------------------------------------------
 238 template< class _TDistanceMap, class _TCostMap >
 239 void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
 240 _Skeleton(
 241   const TDistanceMap* dmap,
 242   const TCostMap* cmap,
 243   const TMST* mst,
 244   const TIndicesData& end_points,
 245   TIndicesData& bifurcations,
 246   TSkeleton* skeleton
 247   )
 248 {
 249   typedef typename TMST::TPath _TPath;
 250
 251   for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
 252   {
 253     typename _TPath::Pointer path;
 254     mst->GetPath( path, *eIt );
 255     skeleton->AddBranch( path );
 256
 257   } // rof
 258 }
 259
 260 #endif // __fpa__Image__SkeletonFilter__hxx__
 261
 262 // eof - $RCSfile$