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