lib/fpa/Image/SkeletonFilter.hxx

   1 #ifndef __fpa__Image__SkeletonFilter__hxx__
   2 #define __fpa__Image__SkeletonFilter__hxx__
   3
   4 #include <fpa/Base/Functors/Inverse.h>
   5 #include <fpa/Image/Functors/SimpleNeighborhood.h>
   6 #include <itkImageRegionIteratorWithIndex.h>
   7
   8 // -------------------------------------------------------------------------
   9 #define fpa_Image_SkeletonFilter_OutputMacro( o_n, o_t )                \
  10   template< class _TImage >                                             \
  11   typename fpa::Image::SkeletonFilter< _TImage >::                      \
  12   o_t* fpa::Image::SkeletonFilter< _TImage >::                          \
  13   Get##o_n( )                                                           \
  14   {                                                                     \
  15     return(                                                             \
  16       dynamic_cast< o_t* >(                                             \
  17         this->itk::ProcessObject::GetOutput( this->m_##o_n##Idx )       \
  18         )                                                               \
  19       );                                                                \
  20   }
  21
  22 fpa_Image_SkeletonFilter_OutputMacro( Skeleton, TSkeleton );
  23 fpa_Image_SkeletonFilter_OutputMacro( Marks, TMarks );
  24
  25 // -------------------------------------------------------------------------
  26 template< class _TImage >
  27 fpa::Image::SkeletonFilter< _TImage >::
  28 SkeletonFilter( )
  29   : Superclass( )
  30 {
  31   typedef typename _TImage::PixelType _TPixel;
  32   typedef fpa::Image::Functors::SimpleNeighborhood< _TImage::ImageDimension > _TNeighFunc;
  33   typedef fpa::Base::Functors::Inverse< _TPixel, _TPixel > _TInvFunc;
  34
  35   unsigned int nOutputs = this->GetNumberOfRequiredOutputs( );
  36   this->SetNumberOfRequiredOutputs( nOutputs + 2 );
  37   this->m_SkeletonIdx = nOutputs;
  38   this->m_MarksIdx = nOutputs + 1;
  39
  40   typename TSkeleton::Pointer skeleton = TSkeleton::New( );
  41   typename TMarks::Pointer marks = TMarks::New( );
  42   this->SetNthOutput( this->m_SkeletonIdx, skeleton.GetPointer( ) );
  43   this->SetNthOutput( this->m_MarksIdx, marks.GetPointer( ) );
  44
  45   typename _TNeighFunc::Pointer nfunc = _TNeighFunc::New( );
  46   nfunc->SetOrder( 2 );
  47   this->SetNeighborhoodFunction( nfunc );
  48
  49   typename _TInvFunc::Pointer ifunc = _TInvFunc::New( );
  50   this->SetConversionFunction( ifunc );
  51 }
  52
  53 // -------------------------------------------------------------------------
  54 template< class _TImage >
  55 fpa::Image::SkeletonFilter< _TImage >::
  56 ~SkeletonFilter( )
  57 {
  58 }
  59
  60 // -------------------------------------------------------------------------
  61 template< class _TImage >
  62 void fpa::Image::SkeletonFilter< _TImage >::
  63 _BeforeGenerateData( )
  64 {
  65   this->Superclass::_BeforeGenerateData( );
  66   this->m_SkeletonQueue.clear( );
  67 }
  68
  69 // -------------------------------------------------------------------------
  70 template< class _TImage >
  71 void fpa::Image::SkeletonFilter< _TImage >::
  72 _UpdateResult( const _TQueueNode& n )
  73 {
  74   typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
  75
  76   this->Superclass::_UpdateResult( n );
  77   double d = double( this->GetInput( )->GetPixel( n.Vertex ) );
  78   if( d >= double( 0 ) )
  79   {
  80     // Update skeleton candidates
  81     d += double( 1e-5 );
  82     double v = double( n.Result ) / ( d * d );
  83     this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) );
  84
  85   } // fi
  86 }
  87
  88 // -------------------------------------------------------------------------
  89 template< class _TImage >
  90 void fpa::Image::SkeletonFilter< _TImage >::
  91 _AfterGenerateData( )
  92 {
  93   this->Superclass::_AfterGenerateData( );
  94
  95   _TAdjacencies A;
  96   std::vector< TIndex > end_points;
  97   this->_EndPoints( end_points, A );
  98   this->_Skeleton( end_points, A );
  99 }
 100
 101 // -------------------------------------------------------------------------
 102 template< class _TImage >
 103 void fpa::Image::SkeletonFilter< _TImage >::
 104 _EndPoints( std::vector< TIndex >& end_points, _TAdjacencies& A )
 105 {
 106   auto marks = this->GetMarks( );
 107   auto mst = this->GetMinimumSpanningTree( );
 108   auto spac = marks->GetSpacing( );
 109
 110   // Some values
 111   marks->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
 112   marks->SetRequestedRegion( mst->GetRequestedRegion( ) );
 113   marks->SetBufferedRegion( mst->GetBufferedRegion( ) );
 114   marks->SetSpacing( mst->GetSpacing( ) );
 115   marks->SetOrigin( mst->GetOrigin( ) );
 116   marks->SetDirection( mst->GetDirection( ) );
 117   marks->Allocate( );
 118   marks->FillBuffer( 0 );
 119
 120   // BFS from maximum queue
 121   while( this->m_SkeletonQueue.size( ) > 0 )
 122   {
 123     // Get node
 124     auto nIt = this->m_SkeletonQueue.begin( );
 125     auto n = *nIt;
 126     this->m_SkeletonQueue.erase( nIt );
 127
 128     // Mark it and update end-points
 129     unsigned char m = marks->GetPixel( n.second );
 130     if( m != 0 )
 131       continue;
 132     marks->SetPixel( n.second, 1 );
 133     end_points.push_back( n.second );
 134
 135     // Mark path
 136     TIndex it = n.second;
 137     TIndex p = mst->GetParent( it );
 138     while( it != p )
 139     {
 140       this->_MarkSphere( it );
 141       it = p;
 142       p = mst->GetParent( it );
 143
 144     } // elihw
 145     this->_MarkSphere( it );
 146     A[ n.second ] = it;
 147
 148   } // elihw
 149 }
 150
 151 // -------------------------------------------------------------------------
 152 template< class _TImage >
 153 void fpa::Image::SkeletonFilter< _TImage >::
 154 _Skeleton( const std::vector< TIndex >& end_points, _TAdjacencies& A )
 155 {
 156   typedef itk::Image< unsigned long, Self::Dimension > _TTagsImage;
 157   typedef typename TMST::TPath _TPath;
 158
 159   auto mst = this->GetMinimumSpanningTree( );
 160   auto skeleton = this->GetSkeleton( );
 161
 162   // Tag branches
 163   typename _TTagsImage::Pointer tags = _TTagsImage::New( );
 164   tags->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
 165   tags->SetRequestedRegion( mst->GetRequestedRegion( ) );
 166   tags->SetBufferedRegion( mst->GetBufferedRegion( ) );
 167   tags->Allocate( );
 168   tags->FillBuffer( 0 );
 169   for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
 170   {
 171     TIndex it = *eIt;
 172     TIndex p = mst->GetParent( it );
 173     while( it != p )
 174     {
 175       tags->SetPixel( it, tags->GetPixel( it ) + 1 );
 176       it = p;
 177       p = mst->GetParent( it );
 178
 179     } // elihw
 180     tags->SetPixel( it, tags->GetPixel( it ) + 1 );
 181
 182   } // rof
 183
 184   // Build paths (branches)
 185   for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
 186   {
 187     TIndex it = *eIt;
 188     TIndex p = mst->GetParent( it );
 189     TIndex sIdx = *eIt;
 190     typename _TPath::Pointer path = _TPath::New( );
 191     path->SetReferenceImage( mst );
 192     while( it != p )
 193     {
 194       if( tags->GetPixel( sIdx ) != tags->GetPixel( it ) )
 195       {
 196         // Ok, a new branch has been added
 197         path->AddVertex( it );
 198         skeleton->AddBranch( path );
 199
 200         // Update a new starting index
 201         path = _TPath::New( );
 202         path->SetReferenceImage( mst );
 203         sIdx = it;
 204       }
 205       else
 206         path->AddVertex( it );
 207       it = p;
 208       p = mst->GetParent( it );
 209
 210     } // elihw
 211
 212     // Finally, add last branch
 213     path->AddVertex( it );
 214     skeleton->AddBranch( path );
 215
 216   } // rof
 217 }
 218
 219 // -------------------------------------------------------------------------
 220 template< class _TImage >
 221 void fpa::Image::SkeletonFilter< _TImage >::
 222 _MarkSphere( const TIndex& idx )
 223 {
 224   typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
 225
 226   static const double _eps = std::sqrt( double( Self::Dimension + 1 ) );
 227   auto input = this->GetInput( );
 228   auto marks = this->GetMarks( );
 229   auto spac = input->GetSpacing( );
 230   auto region = input->GetRequestedRegion( );
 231
 232   typename _TImage::PointType cnt;
 233   input->TransformIndexToPhysicalPoint( idx, cnt );
 234   double r = double( input->GetPixel( idx ) ) * _eps;
 235
 236   TIndex i0, i1;
 237   for( unsigned int d = 0; d < Self::Dimension; ++d )
 238   {
 239     long off = long( std::ceil( r / double( spac[ d ] ) ) );
 240     if( off < 3 )
 241       off = 3;
 242     i0[ d ] = idx[ d ] - off;
 243     i1[ d ] = idx[ d ] + off;
 244
 245     if( i0[ d ] < region.GetIndex( )[ d ] )
 246       i0[ d ] = region.GetIndex( )[ d ];
 247
 248     if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
 249       i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
 250
 251   } // rof
 252
 253   typename _TImage::SizeType size;
 254   for( unsigned int d = 0; d < Self::Dimension; ++d )
 255     size[ d ] = i1[ d ] - i0[ d ] + 1;
 256
 257   typename _TImage::RegionType neighRegion;
 258   neighRegion.SetIndex( i0 );
 259   neighRegion.SetSize( size );
 260
 261   _TMarksIt mIt( marks, neighRegion );
 262   for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 263   {
 264     TIndex w = mIt.GetIndex( );
 265     typename _TImage::PointType p;
 266     marks->TransformIndexToPhysicalPoint( w, p );
 267     mIt.Set( 1 );
 268
 269   } // rof
 270 }
 271
 272 #endif // __fpa__Image__SkeletonFilter__hxx__
 273
 274 // eof - $RCSfile$