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