lib/fpa/Image/SkeletonFilter.hxx

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