lib/fpa/Filters/Image/ExtractSkeleton.hxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5 #ifndef __fpa__Filters__Image__ExtractSkeleton__hxx__
   6 #define __fpa__Filters__Image__ExtractSkeleton__hxx__
   7
   8 #include <itkImage.h>
   9 #include <itkImageRegionIteratorWithIndex.h>
  10 #include <itkMinimumMaximumImageCalculator.h>
  11 #include <fpa/Functors/Dijkstra/Invert.h>
  12
  13 // -------------------------------------------------------------------------
  14 template< class _TInputImage, class _TDistanceMap >
  15 fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  16 _TDijkstra::_TDijkstra( )
  17   : Superclass( )
  18 {
  19   // Prepare weight function
  20   typedef fpa::Functors::Dijkstra::Invert< TScalar > _TWeight;
  21   typename _TWeight::Pointer weight = _TWeight::New( );
  22   weight->SetAlpha( 1 );
  23   weight->SetBeta( 1 );
  24   this->SetWeightFunction( weight );
  25 }
  26
  27 // -------------------------------------------------------------------------
  28 template< class _TInputImage, class _TDistanceMap >
  29 fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  30 _TDijkstra::~_TDijkstra( )
  31 {
  32 }
  33
  34 // -------------------------------------------------------------------------
  35 template< class _TInputImage, class _TDistanceMap >
  36 void fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  37 _TDijkstra::_BeforeGenerateData( )
  38 {
  39   this->Superclass::_BeforeGenerateData( );
  40   this->m_SkeletonQueue.clear( );
  41 }
  42
  43 // -------------------------------------------------------------------------
  44 template< class _TInputImage, class _TDistanceMap >
  45 void fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  46 _TDijkstra::_PostComputeOutputValue( TNode& n )
  47 {
  48   typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
  49
  50   this->Superclass::_PostComputeOutputValue( n );
  51   double d = double( this->GetInput( )->GetPixel( n.Vertex ) );
  52   if( d >= double( 0 ) )
  53   {
  54     // Update skeleton candidates
  55     d += double( 1e-5 );
  56     double v = double( n.Value ) / ( d * d );
  57     this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) );
  58
  59   } // fi
  60 }
  61
  62 // -------------------------------------------------------------------------
  63 template< class _TInputImage, class _TDistanceMap >
  64 itk::ModifiedTimeType
  65 fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  66 GetMTime( ) const
  67 {
  68   itk::ModifiedTimeType t = this->Superclass::GetMTime( );
  69   itk::ModifiedTimeType q = this->m_DistanceMap->GetMTime( );
  70   return( ( q < t )? q: t );
  71 }
  72
  73 // -------------------------------------------------------------------------
  74 template< class _TInputImage, class _TDistanceMap >
  75 void fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  76 SetInput( TInputImage* input )
  77 {
  78   this->Superclass::SetNthInput( 0, input );
  79 }
  80
  81 // -------------------------------------------------------------------------
  82 template< class _TInputImage, class _TDistanceMap >
  83 typename fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  84 TInputImage* fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  85 GetInput( )
  86 {
  87   return( dynamic_cast< TInputImage* >( this->Superclass::GetInput( 0 ) ) );
  88 }
  89
  90 // -------------------------------------------------------------------------
  91 template< class _TInputImage, class _TDistanceMap >
  92 const typename fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  93 TInputImage* fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
  94 GetInput( ) const
  95 {
  96   return(
  97     dynamic_cast< const TInputImage* >( this->Superclass::GetInput( 0 ) )
  98     );
  99 }
 100
 101 // -------------------------------------------------------------------------
 102 template< class _TInputImage, class _TDistanceMap >
 103 typename fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
 104 TSkeleton* fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
 105 GetOutput( )
 106 {
 107   return( dynamic_cast< TSkeleton* >( this->Superclass::GetOutput( 0 ) ) );
 108 }
 109
 110 // -------------------------------------------------------------------------
 111 template< class _TInputImage, class _TDistanceMap >
 112 const typename fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
 113 TSkeleton* fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
 114 GetOutput( ) const
 115 {
 116   return(
 117     dynamic_cast< const TSkeleton* >( this->Superclass::GetOutput( 0 ) )
 118     );
 119 }
 120
 121 // -------------------------------------------------------------------------
 122 template< class _TInputImage, class _TDistanceMap >
 123 fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
 124 ExtractSkeleton( )
 125   : Superclass( ),
 126     m_SeedFromMaximumDistance( false )
 127 {
 128   this->SetNumberOfRequiredInputs( 1 );
 129   this->SetNumberOfRequiredOutputs( 1 );
 130   this->SetNthOutput( 0, TSkeleton::New( ) );
 131
 132   this->m_DistanceMap = TDistanceMap::New( );
 133 }
 134
 135 // -------------------------------------------------------------------------
 136 template< class _TInputImage, class _TDistanceMap >
 137 fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
 138 ~ExtractSkeleton( )
 139 {
 140 }
 141
 142 // -------------------------------------------------------------------------
 143 template< class _TInputImage, class _TDistanceMap >
 144 void fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
 145 GenerateData( )
 146 {
 147   // Update distance map
 148   this->m_DistanceMap->SetInput( this->GetInput( ) );
 149   this->m_DistanceMap->Update( );
 150
 151   // Correct seed
 152   if( this->m_SeedFromMaximumDistance )
 153   {
 154     typedef itk::MinimumMaximumImageCalculator< TOutputImage > _TMinMax;
 155     typename _TMinMax::Pointer minmax = _TMinMax::New( );
 156     minmax->SetImage( this->m_DistanceMap->GetOutput( ) );
 157     minmax->Compute( );
 158     this->m_Seed = minmax->GetIndexOfMaximum( );
 159
 160   } // fi
 161
 162   std::cout << int( this->GetInput( )->GetPixel( this->m_Seed ) ) << std::endl;
 163   std::cout << this->m_DistanceMap->GetOutput( )->GetPixel( this->m_Seed ) << std::endl;
 164
 165   // Compute MST
 166   typename _TDijkstra::Pointer dijkstra = _TDijkstra::New( );
 167   dijkstra->SetInput( this->m_DistanceMap->GetOutput( ) );
 168   dijkstra->AddSeed( this->m_Seed );
 169   dijkstra->Update( );
 170
 171   // Compute end-points
 172   const _TMST* mst = dijkstra->GetMinimumSpanningTree( );
 173   std::vector< TIndex > end_points;
 174   this->_EndPoints(
 175     end_points, this->m_DistanceMap->GetOutput( ),
 176     mst, dijkstra->GetSkeletonQueue( )
 177     );
 178   this->m_EndPoints = end_points;
 179
 180   // Compute symbolic branches
 181   /* TODO
 182      typedef std::map< TIndex, TIndex, typename TIndex::LexicographicCompare > _TTags;
 183      _TTags tags, branches;
 184      typename std::vector< TIndex >::const_iterator eIt = end_points.begin( );
 185      for( ; eIt != end_points.end( ); ++eIt )
 186      {
 187      // Tag path
 188      TIndex it = *eIt;
 189      TIndex p = mst->GetParent( it );
 190      typename _TTags::iterator bIt = tags.end( );
 191      while( it != p && bIt == tags.end( ) )
 192      {
 193      typename _TTags::iterator tIt = tags.find( it );
 194      if( tIt != tags.end( ) )
 195      {
 196      // Ok, a bifurcation point has been found
 197      // branch1: tIt->second <-> it (ok)
 198      branches[ tIt->second ] = it;
 199
 200      // branch2: *eit <-> it (ok)
 201      branches[ *eIt ] = it;
 202
 203      // branch3: it <-> until next bifurcation
 204      bIt = tIt;
 205      }
 206      else
 207      tags[ it ] = *eIt;
 208      it = p;
 209      p = mst->GetParent( it );
 210
 211      } // elihw
 212      if( bIt != tags.end( ) )
 213      {
 214      TIndex pTag = bIt->second;
 215      TIndex nTag = bIt->first;
 216      it = bIt->first;
 217      p = it;
 218      while( tags[ it ] == pTag )
 219      {
 220      tags[ it ] = nTag;
 221      p = it;
 222      it = mst->GetParent( it );
 223
 224      } // elihw
 225      tags[ it ] = nTag;
 226      branches[ bIt->first ] = p;
 227      }
 228      else
 229      {
 230      tags[ it ] = *eIt;
 231      branches[ *eIt ] = it;
 232
 233      } // fi
 234
 235      } // rof
 236
 237      // Fill full branches
 238      typedef typename _TMST::TVertices _TVertices;
 239      typedef typename TSkeleton::TPath _TPath;
 240
 241      TSkeleton* sk = this->GetOutput( );
 242      typename _TTags::const_iterator bIt = branches.begin( );
 243      for( ; bIt != branches.end( ); ++bIt )
 244      {
 245      _TVertices v = mst->GetAxis( bIt->first, bIt->second );
 246      typename _TPath::Pointer path = _TPath::New( );
 247      path->SetReferenceImage( this->GetInput( ) );
 248      typename _TVertices::const_reverse_iterator vIt = v.rbegin( );
 249      for( ; vIt != v.rend( ); ++vIt )
 250      path->AddVertex( *vIt );
 251      sk->AddBranch( path );
 252
 253      } // rof
 254   */
 255 }
 256
 257 // -------------------------------------------------------------------------
 258 template< class _TInputImage, class _TDistanceMap >
 259 template< class _TMarksPointer >
 260 void fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
 261 _MarkSphere(
 262   _TMarksPointer& marks, const TOutputImage* dmap, const TIndex& center
 263   )
 264 {
 265   typedef typename _TMarksPointer::ObjectType          _TMarks;
 266   typedef itk::ImageRegionIteratorWithIndex< _TMarks > _TMarksIt;
 267
 268   static const double _eps = std::sqrt( double( Self::Dimension + 1 ) );
 269   typename _TMarks::SpacingType spac = dmap->GetSpacing( );
 270   typename _TMarks::RegionType region = dmap->GetRequestedRegion( );
 271
 272   typename _TMarks::PointType cnt;
 273   dmap->TransformIndexToPhysicalPoint( center, cnt );
 274   double r = double( dmap->GetPixel( center ) ) * _eps;
 275
 276   TIndex i0, i1;
 277   for( unsigned int d = 0; d < Self::Dimension; ++d )
 278   {
 279     long off = long( std::ceil( r / double( spac[ d ] ) ) );
 280     if( off < 3 )
 281       off = 3;
 282     i0[ d ] = center[ d ] - off;
 283     i1[ d ] = center[ d ] + off;
 284
 285     if( i0[ d ] < region.GetIndex( )[ d ] )
 286       i0[ d ] = region.GetIndex( )[ d ];
 287
 288     if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
 289       i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
 290
 291   } // rof
 292
 293   typename _TMarks::SizeType size;
 294   for( unsigned int d = 0; d < Self::Dimension; ++d )
 295     size[ d ] = i1[ d ] - i0[ d ] + 1;
 296
 297   typename _TMarks::RegionType neighRegion;
 298   neighRegion.SetIndex( i0 );
 299   neighRegion.SetSize( size );
 300
 301   _TMarksIt mIt( marks, neighRegion );
 302   for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 303     mIt.Set( true );
 304 }
 305
 306 // -------------------------------------------------------------------------
 307 template< class _TInputImage, class _TDistanceMap >
 308 void fpa::Filters::Image::ExtractSkeleton< _TInputImage, _TDistanceMap >::
 309 _EndPoints(
 310   std::vector< TIndex >& end_points,
 311   const TOutputImage* dmap,
 312   const _TMST* mst,
 313   const _TSkeletonQueue& queue
 314   )
 315 {
 316   typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
 317
 318   // Some values
 319   typedef itk::Image< bool, Self::Dimension > _TMarks;
 320   typename _TMarks::Pointer marks = _TMarks::New( );
 321   marks->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
 322   marks->SetRequestedRegion( mst->GetRequestedRegion( ) );
 323   marks->SetBufferedRegion( mst->GetBufferedRegion( ) );
 324   marks->SetSpacing( mst->GetSpacing( ) );
 325   marks->SetOrigin( mst->GetOrigin( ) );
 326   marks->SetDirection( mst->GetDirection( ) );
 327   marks->Allocate( );
 328   marks->FillBuffer( false );
 329
 330   // Get candidates in maximum to minimum iteration
 331   typename _TSkeletonQueue::const_reverse_iterator nIt = queue.rbegin( );
 332   for( ; nIt != queue.rend( ); ++nIt )
 333   {
 334     // Mark it and update end-points
 335     if( !( marks->GetPixel( nIt->second ) ) )
 336     {
 337       marks->SetPixel( nIt->second, true );
 338       end_points.push_back( nIt->second );
 339
 340       // Mark path
 341       TIndex it = nIt->second;
 342       TIndex p = mst->GetParent( it );
 343       while( it != p )
 344       {
 345         this->_MarkSphere( marks, dmap, it );
 346         it = p;
 347         p = mst->GetParent( it );
 348
 349       } // elihw
 350       this->_MarkSphere( marks, dmap, it );
 351
 352     } // fi
 353
 354   } // rof
 355 }
 356
 357 #endif // __fpa__Filters__Image__ExtractSkeleton__hxx__
 358
 359 // eof - $RCSfile$