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
 112   // 3. Build skeleton and keep track of bifurcations
 113   this->_Skeleton( dmap, cmap, mst, ep->Get( ), bi->Get( ), sk );
 114 }
 115
 116 // -------------------------------------------------------------------------
 117 template< class _TDistanceMap, class _TCostMap >
 118 void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
 119 _EndPoints(
 120   const TDistanceMap* dmap,
 121   const TCostMap* cmap,
 122   const TMST* mst,
 123   TIndicesData& end_points
 124   )
 125 {
 126   typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt;
 127   typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt;
 128   typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
 129   typedef typename _TQueue::value_type _TQueueValue;
 130   typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
 131
 132   static const double _eps = std::sqrt( double( TMarks::ImageDimension + 1 ) );
 133
 134   // Some values
 135   auto marks = this->GetMarks( );
 136   marks->SetLargestPossibleRegion( dmap->GetLargestPossibleRegion( ) );
 137   marks->SetRequestedRegion( dmap->GetRequestedRegion( ) );
 138   marks->SetBufferedRegion( dmap->GetBufferedRegion( ) );
 139   marks->SetSpacing( dmap->GetSpacing( ) );
 140   marks->SetOrigin( dmap->GetOrigin( ) );
 141   marks->SetDirection( dmap->GetDirection( ) );
 142   marks->Allocate( );
 143   marks->FillBuffer( 0 );
 144
 145   // Create queue
 146   _TQueue queue;
 147   _TDistMapIt dIt( dmap, dmap->GetRequestedRegion( ) );
 148   _TCostMapIt cIt( cmap, cmap->GetRequestedRegion( ) );
 149   dIt.GoToBegin( );
 150   cIt.GoToBegin( );
 151   for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
 152   {
 153     double d = double( dIt.Get( ) );
 154     if( d >= double( 0 ) )
 155     {
 156       d += double( 1e-5 );
 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   while( queue.size( ) > 0 )
 167   {
 168     // Get node
 169     auto nIt = queue.begin( );
 170     auto n = *nIt;
 171     queue.erase( nIt );
 172
 173     // Mark it and update end-points
 174     unsigned char m = marks->GetPixel( n.second );
 175     if( m != 0 )
 176       continue;
 177     marks->SetPixel( n.second, 1 );
 178     end_points.insert( n.second );
 179
 180     // Mark path
 181     auto spac = marks->GetSpacing( );
 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 r = double( dmap->GetPixel( idx ) ) * _eps;
 190
 191       TIndex i0, i1;
 192       for( unsigned int d = 0; d < TMarks::ImageDimension; ++d )
 193       {
 194         long off = long( std::ceil( r / double( spac[ d ] ) ) );
 195         if( off < 3 )
 196           off = 3;
 197         i0[ d ] = idx[ d ] - off;
 198         i1[ d ] = idx[ d ] + off;
 199
 200         if( i0[ d ] < region.GetIndex( )[ d ] )
 201           i0[ d ] = region.GetIndex( )[ d ];
 202
 203         if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
 204           i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
 205
 206       } // rof
 207
 208       typename TMarks::SizeType size;
 209       for( unsigned int d = 0; d < TMarks::ImageDimension; ++d )
 210         size[ d ] = i1[ d ] - i0[ d ] + 1;
 211
 212       typename TMarks::RegionType neighRegion;
 213       neighRegion.SetIndex( i0 );
 214       neighRegion.SetSize( size );
 215
 216       _TMarksIt mIt( marks, neighRegion );
 217       for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 218       {
 219         TIndex w = mIt.GetIndex( );
 220         typename _TCostMap::PointType p;
 221         marks->TransformIndexToPhysicalPoint( w, p );
 222         mIt.Set( 1 );
 223
 224       } // rof
 225
 226     } // rof
 227
 228   } // elihw
 229 }
 230
 231 // -------------------------------------------------------------------------
 232 template< class _TDistanceMap, class _TCostMap >
 233 void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
 234 _Skeleton(
 235   const TDistanceMap* dmap,
 236   const TCostMap* cmap,
 237   const TMST* mst,
 238   const TIndicesData& end_points,
 239   TIndicesData& bifurcations,
 240   TSkeleton* skeleton
 241   )
 242 {
 243   typedef typename TMST::TPath _TPath;
 244
 245   for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
 246   {
 247     typename _TPath::Pointer path;
 248     mst->GetPath( path, *eIt );
 249     skeleton->AddBranch( path );
 250
 251   } // rof
 252 }
 253
 254 #endif // __fpa__Image__SkeletonFilter__hxx__
 255
 256 // eof - $RCSfile$