lib/fpa/Image/SkeletonFilter.hxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #ifndef __fpa__Image__SkeletonFilter__hxx__
   7 #define __fpa__Image__SkeletonFilter__hxx__
   8
   9 #include <itkImageRegionIteratorWithIndex.h>
  10 #include <itkMinimumMaximumImageCalculator.h>
  11
  12 #include <fpa/Image/Functors/VertexIdentity.h>
  13 #include <fpa/Base/Functors/InvertValue.h>
  14
  15 // -------------------------------------------------------------------------
  16 template< class _TImage, class _TScalar >
  17 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  18 SetInput( const TScalarImage* image )
  19 {
  20   // Do nothing
  21 }
  22
  23 // -------------------------------------------------------------------------
  24 template< class _TImage, class _TScalar >
  25 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  26 SetInput( const TImage* image )
  27 {
  28   this->m_DistanceMap->SetInput( image );
  29   this->Superclass::SetInput( this->m_DistanceMap->GetOutput( ) );
  30 }
  31
  32 // -------------------------------------------------------------------------
  33 template< class _TImage, class _TScalar >
  34 _TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  35 GetInput( )
  36 {
  37   return( this->m_DistanceMap->GetInput( ) );
  38 }
  39
  40 // -------------------------------------------------------------------------
  41 template< class _TImage, class _TScalar >
  42 const _TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  43 GetInput( ) const
  44 {
  45   return( this->m_DistanceMap->GetInput( ) );
  46 }
  47
  48 // -------------------------------------------------------------------------
  49 template< class _TImage, class _TScalar >
  50 typename fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  51 TSkeleton* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  52 GetSkeleton( )
  53 {
  54   return(
  55     dynamic_cast< TSkeleton* >(
  56       this->itk::ProcessObject::GetOutput( this->m_SkeletonIdx )
  57       )
  58     );
  59 }
  60
  61 // -------------------------------------------------------------------------
  62 template< class _TImage, class _TScalar >
  63 typename fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  64 TMarks* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  65 GetMarks( )
  66 {
  67   return(
  68     dynamic_cast< TMarks* >(
  69       this->itk::ProcessObject::GetOutput( this->m_MarksIdx )
  70       )
  71     );
  72 }
  73
  74 // -------------------------------------------------------------------------
  75 template< class _TImage, class _TScalar >
  76 fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  77 SkeletonFilter( )
  78   : Superclass( ),
  79     m_SeedFromMaximumDistance( false )
  80 {
  81   typedef fpa::Image::Functors::VertexIdentity< TScalarImage, _TScalar > _TVertexFunc;
  82   typedef fpa::Base::Functors::InvertValue< _TScalar, _TScalar > _TValueFunc;
  83
  84   typename _TVertexFunc::Pointer vertex_func = _TVertexFunc::New( );
  85   typename _TValueFunc::Pointer value_func = _TValueFunc::New( );
  86   value_func->SetAlpha( 1 );
  87   value_func->SetBeta( 1 );
  88
  89   this->SetFunctor( vertex_func );
  90   this->SetFunctor( value_func );
  91
  92   this->m_DistanceMap = TDistanceMap::New( );
  93   this->m_DistanceMap->InsideIsPositiveOn( );
  94   this->m_DistanceMap->SquaredDistanceOff( );
  95   this->m_DistanceMap->UseImageSpacingOn( );
  96
  97   unsigned int nOutputs = this->GetNumberOfRequiredOutputs( );
  98   this->SetNumberOfRequiredOutputs( nOutputs + 2 );
  99   this->m_SkeletonIdx = nOutputs;
 100   this->m_MarksIdx = nOutputs + 1;
 101
 102   typename TSkeleton::Pointer skeleton = TSkeleton::New( );
 103   typename TMarks::Pointer marks = TMarks::New( );
 104   this->SetNthOutput( this->m_SkeletonIdx, skeleton.GetPointer( ) );
 105   this->SetNthOutput( this->m_MarksIdx, marks.GetPointer( ) );
 106 }
 107
 108 // -------------------------------------------------------------------------
 109 template< class _TImage, class _TScalar >
 110 fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 111 ~SkeletonFilter( )
 112 {
 113 }
 114
 115 // -------------------------------------------------------------------------
 116 template< class _TImage, class _TScalar >
 117 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 118 GenerateData( )
 119 {
 120   // Prepare distance map
 121   this->m_DistanceMap->Update( );
 122
 123   // Update start seed
 124   TVertex seed;
 125   if( this->m_SeedFromMaximumDistance )
 126   {
 127     typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
 128     typename _TMinMax::Pointer minmax = _TMinMax::New( );
 129     minmax->SetImage( this->m_DistanceMap->GetOutput( ) );
 130     minmax->Compute( );
 131     seed = minmax->GetIndexOfMaximum( );
 132   }
 133   else
 134     seed = *( this->BeginSeeds( ) );
 135   this->m_Seeds.clear( );
 136   this->m_Seeds.insert( seed );
 137
 138   // Prepare skeleton candidates queue
 139   this->m_SkeletonQueue.clear( );
 140
 141   // Go!
 142   this->Superclass::GenerateData( );
 143
 144   // Backtracking
 145   _TAdjacencies A;
 146   std::vector< TVertex > end_points;
 147   this->_EndPoints( end_points, A );
 148   this->_Skeleton( end_points, A );
 149 }
 150
 151 // -------------------------------------------------------------------------
 152 template< class _TImage, class _TScalar >
 153 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 154 _SetOutputValue( const TVertex& vertex, const TOutputValue& value )
 155 {
 156   typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
 157
 158   this->Superclass::_SetOutputValue( vertex, value );
 159   double d = double( this->m_DistanceMap->GetOutput( )->GetPixel( vertex ) );
 160   if( d >= double( 0 ) )
 161   {
 162     // Update skeleton candidates
 163     d += double( 1e-5 );
 164     double v = double( value ) / ( d * d );
 165     this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, vertex ) );
 166
 167   } // fi
 168 }
 169
 170 // -------------------------------------------------------------------------
 171 template< class _TImage, class _TScalar >
 172 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 173 _EndPoints( std::vector< TVertex >& end_points, _TAdjacencies& A )
 174 {
 175   typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
 176
 177   TMarks* marks = this->GetMarks( );
 178   TMST* mst = this->GetMinimumSpanningTree( );
 179   typename TImage::SpacingType spac = marks->GetSpacing( );
 180
 181   // Some values
 182   marks->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
 183   marks->SetRequestedRegion( mst->GetRequestedRegion( ) );
 184   marks->SetBufferedRegion( mst->GetBufferedRegion( ) );
 185   marks->SetSpacing( mst->GetSpacing( ) );
 186   marks->SetOrigin( mst->GetOrigin( ) );
 187   marks->SetDirection( mst->GetDirection( ) );
 188   marks->Allocate( );
 189   marks->FillBuffer( 0 );
 190
 191   // BFS from maximum queue
 192   while( this->m_SkeletonQueue.size( ) > 0 )
 193   {
 194     // Get node
 195     typename _TSkeletonQueue::iterator nIt = this->m_SkeletonQueue.begin( );
 196     _TSkeletonQueueValue n = *nIt;
 197     this->m_SkeletonQueue.erase( nIt );
 198
 199     // Mark it and update end-points
 200     unsigned char m = marks->GetPixel( n.second );
 201     if( m != 0 )
 202       continue;
 203     marks->SetPixel( n.second, 1 );
 204     end_points.push_back( n.second );
 205
 206     // Mark path
 207     TVertex it = n.second;
 208     TVertex p = mst->GetParent( it );
 209     while( it != p )
 210     {
 211       this->_MarkSphere( it );
 212       it = p;
 213       p = mst->GetParent( it );
 214
 215     } // elihw
 216     this->_MarkSphere( it );
 217     A[ n.second ] = it;
 218
 219   } // elihw
 220 }
 221
 222 // -------------------------------------------------------------------------
 223 template< class _TImage, class _TScalar >
 224 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 225 _Skeleton( const std::vector< TVertex >& end_points, _TAdjacencies& A )
 226 {
 227   typedef typename TSkeleton::TPath _TPath;
 228   typedef itk::Image< unsigned long, TImage::ImageDimension > _TTagsImage;
 229
 230   TMST* mst = this->GetMinimumSpanningTree( );
 231   TSkeleton* skeleton = this->GetSkeleton( );
 232
 233   // Tag branches
 234   typename _TTagsImage::Pointer tags = _TTagsImage::New( );
 235   tags->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
 236   tags->SetRequestedRegion( mst->GetRequestedRegion( ) );
 237   tags->SetBufferedRegion( mst->GetBufferedRegion( ) );
 238   tags->Allocate( );
 239   tags->FillBuffer( 0 );
 240   for( TVertex it: end_points )
 241   {
 242     TVertex p = mst->GetParent( it );
 243     while( it != p )
 244     {
 245       tags->SetPixel( it, tags->GetPixel( it ) + 1 );
 246       it = p;
 247       p = mst->GetParent( it );
 248
 249     } // elihw
 250     tags->SetPixel( it, tags->GetPixel( it ) + 1 );
 251
 252   } // rof
 253
 254   // Build paths (branches)
 255   for( TVertex it: end_points )
 256   {
 257     TVertex p = mst->GetParent( it );
 258     TVertex sIdx = it;
 259     typename _TPath::Pointer path = _TPath::New( );
 260     path->SetReferenceImage( mst );
 261     while( it != p )
 262     {
 263       if( tags->GetPixel( sIdx ) != tags->GetPixel( it ) )
 264       {
 265         // Ok, a new branch has been added
 266         path->AddVertex( it );
 267         skeleton->AddBranch( path );
 268
 269         // Update a new starting index
 270         path = _TPath::New( );
 271         path->SetReferenceImage( mst );
 272         sIdx = it;
 273       }
 274       else
 275         path->AddVertex( it );
 276       it = p;
 277       p = mst->GetParent( it );
 278
 279     } // elihw
 280
 281     // Finally, add last branch
 282     path->AddVertex( it );
 283     skeleton->AddBranch( path );
 284
 285   } // rof
 286 }
 287
 288 // -------------------------------------------------------------------------
 289 template< class _TImage, class _TScalar >
 290 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 291 _MarkSphere( const TVertex& idx )
 292 {
 293   typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
 294
 295   static const double _eps = std::sqrt( double( TImage::ImageDimension + 1 ) );
 296   const TScalarImage* dmap = this->m_DistanceMap->GetOutput( );
 297   TMarks* marks = this->GetMarks( );
 298   typename TImage::SpacingType spac = dmap->GetSpacing( );
 299   typename TImage::RegionType region = dmap->GetRequestedRegion( );
 300
 301   typename TImage::PointType cnt;
 302   dmap->TransformIndexToPhysicalPoint( idx, cnt );
 303   double r = double( dmap->GetPixel( idx ) ) * _eps;
 304
 305   TVertex i0, i1;
 306   for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 307   {
 308     long off = long( std::ceil( r / double( spac[ d ] ) ) );
 309     if( off < 3 )
 310       off = 3;
 311     i0[ d ] = idx[ d ] - off;
 312     i1[ d ] = idx[ d ] + off;
 313
 314     if( i0[ d ] < region.GetIndex( )[ d ] )
 315       i0[ d ] = region.GetIndex( )[ d ];
 316
 317     if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
 318       i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
 319
 320   } // rof
 321
 322   typename TImage::SizeType size;
 323   for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 324     size[ d ] = i1[ d ] - i0[ d ] + 1;
 325
 326   typename TImage::RegionType neighRegion;
 327   neighRegion.SetIndex( i0 );
 328   neighRegion.SetSize( size );
 329
 330   _TMarksIt mIt( marks, neighRegion );
 331   for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 332   {
 333     TVertex w = mIt.GetIndex( );
 334     typename TImage::PointType p;
 335     marks->TransformIndexToPhysicalPoint( w, p );
 336     mIt.Set( 1 );
 337
 338   } // rof
 339 }
 340
 341 #endif // __fpa__Image__SkeletonFilter__hxx__
 342
 343 // eof - $RCSfile$