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