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