lib/fpa/Image/DijkstraWithEndPointDetection.hxx

   1 #ifndef __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
   2 #define __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
   3
   4 #include <vector>
   5 #include <set>
   6 #include <itkImageRegionConstIteratorWithIndex.h>
   7 #include <itkImageRegionIteratorWithIndex.h>
   8
   9 // -------------------------------------------------------------------------
  10 template< class I, class O >
  11 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
  12 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
  13 GetLabelImage( )
  14 {
  15   return(
  16     dynamic_cast< TLabelImage* >(
  17       this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
  18       )
  19     );
  20 }
  21
  22 // -------------------------------------------------------------------------
  23 template< class I, class O >
  24 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
  25 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
  26 GetLabelImage( ) const
  27 {
  28   return(
  29     dynamic_cast< const TLabelImage* >(
  30       this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
  31       )
  32     );
  33 }
  34
  35 // -------------------------------------------------------------------------
  36 template< class I, class O >
  37 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
  38 GraftLabelImage( itk::DataObject* obj )
  39 {
  40   TLabelImage* lbl =
  41     dynamic_cast< TLabelImage* >(
  42       this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
  43       );
  44   if( lbl != NULL )
  45     this->GraftNthOutput( this->m_LabelImageIndex, lbl );
  46 }
  47
  48 // -------------------------------------------------------------------------
  49 template< class I, class O >
  50 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
  51 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
  52 GetEndPoints( )
  53 {
  54   return(
  55     dynamic_cast< TUniqueVertices* >(
  56       this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex )
  57       )
  58     );
  59 }
  60
  61 // -------------------------------------------------------------------------
  62 template< class I, class O >
  63 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
  64 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
  65 GetEndPoints( ) const
  66 {
  67   return(
  68     dynamic_cast< const TUniqueVertices* >(
  69       this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex )
  70       )
  71     );
  72 }
  73
  74 // -------------------------------------------------------------------------
  75 template< class I, class O >
  76 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
  77 GraftEndPoints( itk::DataObject* obj )
  78 {
  79   TUniqueVertices* lbl =
  80     dynamic_cast< TUniqueVertices* >(
  81       this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex )
  82       );
  83   if( lbl != NULL )
  84     this->GraftNthOutput( this->m_EndPointsIndex, lbl );
  85 }
  86
  87 // -------------------------------------------------------------------------
  88 template< class I, class O >
  89 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
  90 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
  91 GetBifurcations( )
  92 {
  93   return(
  94     dynamic_cast< TUniqueVertices* >(
  95       this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex )
  96       )
  97     );
  98 }
  99
 100 // -------------------------------------------------------------------------
 101 template< class I, class O >
 102 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 103 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
 104 GetBifurcations( ) const
 105 {
 106   return(
 107     dynamic_cast< const TUniqueVertices* >(
 108       this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex )
 109       )
 110     );
 111 }
 112
 113 // -------------------------------------------------------------------------
 114 template< class I, class O >
 115 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 116 GraftBifurcations( itk::DataObject* obj )
 117 {
 118   TUniqueVertices* lbl =
 119     dynamic_cast< TUniqueVertices* >(
 120       this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex )
 121       );
 122   if( lbl != NULL )
 123     this->GraftNthOutput( this->m_BifurcationsIndex, lbl );
 124 }
 125
 126
 127 // -------------------------------------------------------------------------
 128 template< class I, class O >
 129 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 130 TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >::
 131 GetBranches( )
 132 {
 133   return(
 134     dynamic_cast< TBranches* >(
 135       this->itk::ProcessObject::GetOutput( this->m_BranchesIndex )
 136       )
 137     );
 138 }
 139
 140 // -------------------------------------------------------------------------
 141 template< class I, class O >
 142 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 143 TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >::
 144 GetBranches( ) const
 145 {
 146   return(
 147     dynamic_cast< const TBranches* >(
 148       this->itk::ProcessObject::GetOutput( this->m_BranchesIndex )
 149       )
 150     );
 151 }
 152
 153 // -------------------------------------------------------------------------
 154 template< class I, class O >
 155 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 156 GraftBranches( itk::DataObject* obj )
 157 {
 158   TBranches* lbl =
 159     dynamic_cast< TBranches* >(
 160       this->itk::ProcessObject::GetOutput( this->m_BranchesIndex )
 161       );
 162   if( lbl != NULL )
 163     this->GraftNthOutput( this->m_BranchesIndex, lbl );
 164 }
 165
 166 // -------------------------------------------------------------------------
 167 template< class I, class O >
 168 fpa::Image::DijkstraWithEndPointDetection< I, O >::
 169 DijkstraWithEndPointDetection( )
 170   : Superclass( ),
 171     m_CorrectSeeds( true ),
 172     m_CorrectEndPoints( true ),
 173     m_SafetyNeighborhoodSize( 3 )
 174 {
 175   this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( );
 176   this->m_EndPointsIndex = this->m_LabelImageIndex + 1;
 177   this->m_BifurcationsIndex = this->m_LabelImageIndex + 2;
 178   this->m_BranchesIndex = this->m_LabelImageIndex + 3;
 179   this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 4 );
 180   this->itk::ProcessObject::SetNthOutput(
 181     this->m_LabelImageIndex, TLabelImage::New( )
 182     );
 183   this->itk::ProcessObject::SetNthOutput(
 184     this->m_EndPointsIndex, TUniqueVertices::New( )
 185     );
 186   this->itk::ProcessObject::SetNthOutput(
 187     this->m_BifurcationsIndex, TUniqueVertices::New( )
 188     );
 189   this->itk::ProcessObject::SetNthOutput(
 190     this->m_BranchesIndex, TBranches::New( )
 191     );
 192 }
 193
 194 // -------------------------------------------------------------------------
 195 template< class I, class O >
 196 fpa::Image::DijkstraWithEndPointDetection< I, O >::
 197 ~DijkstraWithEndPointDetection( )
 198 {
 199 }
 200
 201 // -------------------------------------------------------------------------
 202 template< class I, class O >
 203 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 204 _BeforeGenerateData( )
 205 {
 206   if( this->m_CorrectSeeds )
 207   {
 208     const I* input = this->GetInput( );
 209     typename I::SpacingType spac = input->GetSpacing( );
 210     double max_spac = double( spac[ 0 ] );
 211     for( unsigned int d = 1; d < I::ImageDimension; ++d )
 212       max_spac =
 213         ( double( spac[ d ] ) > max_spac )? double( spac[ d ] ): max_spac;
 214     max_spac *= double( 3 );
 215
 216     // Correct seeds
 217     for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s )
 218     {
 219       TVertex seed = this->m_SeedVertices[ s ];
 220       _TNode n = this->m_Seeds[ seed ];
 221       _TRegion region = this->_Region( seed, max_spac );
 222       itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 223
 224       iIt.GoToBegin( );
 225       _TPixel max_value = iIt.Get( );
 226       for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 227       {
 228         if( iIt.Get( ) > max_value )
 229         {
 230           this->m_SeedVertices[ s ] = iIt.GetIndex( );
 231           max_value = iIt.Get( );
 232
 233         } // fi
 234
 235       } // rof
 236       this->m_Seeds.erase( seed );
 237       n.Parent = this->m_SeedVertices[ s ];
 238       this->m_Seeds[ this->m_SeedVertices[ s ] ] = n;
 239
 240     } // rof
 241
 242   } // fi
 243
 244   // End initialization
 245   this->Superclass::_BeforeGenerateData( );
 246   this->m_Candidates.clear( );
 247 }
 248
 249 // -------------------------------------------------------------------------
 250 template< class I, class O >
 251 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 252 _AfterGenerateData( )
 253 {
 254   // Finish base algorithm
 255   this->Superclass::_AfterGenerateData( );
 256
 257   // Check if backtracking has any sense
 258   if( this->m_Candidates.size( ) == 0 )
 259     return;
 260   this->InvokeEvent( TEndEvent( ) );
 261   this->InvokeEvent( TStartBacktrackingEvent( ) );
 262
 263   // Detect endpoints and bifurcations
 264   this->_EndPointsAndBifurcations( );
 265
 266   // Find branches
 267   this->_FindBranches( );
 268
 269   // Label pixels and branches
 270   this->_LabelAll( );
 271 }
 272
 273 // -------------------------------------------------------------------------
 274 template< class I, class O >
 275 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 276 _SetResult( const TVertex& v, const _TNode& n )
 277 {
 278   this->Superclass::_SetResult( v, n );
 279   this->m_Candidates.insert( _TCandidate( n.Result, v ) );
 280 }
 281
 282 // -------------------------------------------------------------------------
 283 template< class I, class O >
 284 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 285 _EndPointsAndBifurcations( )
 286 {
 287   // Prepare backtracking objects
 288   TUniqueVertices* endpoints = this->GetEndPoints( );
 289   TUniqueVertices* bifurcations = this->GetBifurcations( );
 290   endpoints->Clear( );
 291   bifurcations->Clear( );
 292
 293   // Get some input values
 294   TVertex seed = this->GetSeed( 0 );
 295   const I* input = this->GetInput( );
 296   const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
 297
 298   // Prepare a pixel marks structure
 299   typedef std::map< TVertex, bool, TVertexCompare > _TPixelMarks;
 300   _TPixelMarks pixel_marks;
 301
 302   // Object to hold tree-nodes marks
 303   typedef std::set< TVertex, TVertexCompare > _TTreeMarks;
 304   _TTreeMarks tree_marks;
 305
 306   // Iterate over the candidates, starting from the candidates that
 307   // are near thin branches (high accumulated cost)
 308   typename _TCandidates::const_reverse_iterator cIt =
 309     this->m_Candidates.rbegin( );
 310   for( ; cIt != this->m_Candidates.rend( ); ++cIt )
 311   {
 312     // If pixel has already been visited, pass
 313     TVertex v = cIt->second;
 314     if( pixel_marks[ v ] )
 315       continue;
 316
 317     // Correct it to nearest start candidate (high distance value)
 318     double vr = std::sqrt( double( input->GetPixel( v ) ) );
 319     if( this->m_CorrectEndPoints )
 320       v = this->_MaxInRegion( input, v, vr * double( 1.5 ) );
 321
 322     // Now, check for real marking conditions
 323     //   1. Has it been visited by dijkstra?
 324     if( this->_Node( v ).Label == Self::FarLabel )
 325       continue;
 326     //   2. Is it already marked?
 327     if( pixel_marks[ v ] )
 328       continue;
 329     //   3. Is it already an endpoint?
 330     if( endpoints->Find( v ) )
 331       continue;
 332     //   4. Ok, it is completely new!
 333     endpoints->Insert( v );
 334
 335     // Get the path all the way to global seed
 336     TVertices path = mst->GetPath( v, seed );
 337
 338     // Backtracking to find endpoints and bifurcations
 339     bool adding_new_points = true;
 340     typename TVertices::const_iterator pIt = path.begin( );
 341     for( ; pIt != path.end( ) && adding_new_points; ++pIt )
 342     {
 343       this->InvokeEvent( TBacktrackingEvent( *pIt, 1 ) );
 344       if( tree_marks.find( *pIt ) == tree_marks.end( ) )
 345       {
 346         // Mark current point as a tree point
 347         tree_marks.insert( *pIt );
 348
 349         // Mark a region around current point as visited
 350         vr = std::sqrt( double( input->GetPixel( *pIt ) ) );
 351         _TRegion region = this->_Region( *pIt, vr );
 352         if( region.GetNumberOfPixels( ) > 0 )
 353         {
 354           itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 355           for( iIt.GoToBegin( ); !iIt.IsAtEnd( ); ++iIt )
 356             pixel_marks[ iIt.GetIndex( ) ] = true;
 357         }
 358         else
 359           pixel_marks[ *pIt ] = true;
 360       }
 361       else
 362       {
 363         // A bifurcation point has been reached!
 364         if( *pIt != seed )
 365         {
 366           bifurcations->Insert( *pIt );
 367           adding_new_points = false;
 368
 369         } // fi
 370
 371       } // fi
 372
 373     } // rof
 374     this->InvokeEvent( TEndBacktrackingEvent( ) );
 375
 376   } // rof
 377 }
 378
 379 // -------------------------------------------------------------------------
 380 template< class I, class O >
 381 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 382 _FindBranches( )
 383 {
 384   // Configure and obtain inputs
 385   TVertex seed = this->GetSeed( 0 );
 386   const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
 387   const TUniqueVertices* endpoints = this->GetEndPoints( );
 388   const TUniqueVertices* bifurcations = this->GetBifurcations( );
 389   TBranches* branches = this->GetBranches( );
 390   branches->Clear( );
 391
 392   // Reconstruct pixels
 393   typename TUniqueVertices::ConstIterator eIt = endpoints->Begin( );
 394   for( ; eIt != endpoints->End( ); ++eIt )
 395   {
 396     // Get the path all the way to global seed
 397     TVertices path = mst->GetPath( *eIt, seed );
 398
 399     TVertex start_vertex = *eIt;
 400     typename TVertices::const_iterator pIt = path.begin( );
 401     for( ; pIt != path.end( ); ++pIt )
 402     {
 403       if( bifurcations->Find( *pIt ) )
 404       {
 405         branches->SetValue( start_vertex, *pIt, 1 );
 406         start_vertex = *pIt;
 407
 408       } // fi
 409
 410     } // rof
 411
 412     // Finish with branches to global seed
 413     branches->SetValue( start_vertex, seed, 1 );
 414
 415   } // rof
 416 }
 417
 418 // -------------------------------------------------------------------------
 419 template< class I, class O >
 420 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 421 _LabelAll( )
 422 {
 423   // Configure and obtain inputs
 424   const I* input = this->GetInput( );
 425   const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
 426   const TUniqueVertices* bifurcations = this->GetBifurcations( );
 427   TBranches* branches = this->GetBranches( );
 428   TLabelImage* label = this->GetLabelImage( );
 429   label->FillBuffer( 0 );
 430
 431   // Label branches
 432   typename TBranches::Iterator bIt = branches->Begin( );
 433   TLabel actual_label = 0;
 434   for( ; bIt != branches->End( ); ++bIt )
 435   {
 436     typename TBranches::RowIterator brIt = branches->Begin( bIt );
 437     for( ; brIt != branches->End( bIt ); ++brIt )
 438     {
 439       actual_label++;
 440       brIt->second = actual_label;
 441
 442       TVertices path = mst->GetPath( bIt->first, brIt->first );
 443       typename TVertices::const_iterator pIt = path.begin( );
 444       for( ; pIt != path.end( ); ++pIt )
 445       {
 446         double d = std::sqrt( double( input->GetPixel( *pIt ) ) );
 447         _TRegion region = this->_Region( *pIt, d );
 448         itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 449         itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
 450         iIt.GoToBegin( );
 451         lIt.GoToBegin( );
 452         for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
 453         {
 454           // Mask real input image
 455           if(
 456             !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) &&
 457             lIt.Get( ) == TLabel( 0 )
 458             )
 459             lIt.Set( actual_label );
 460
 461         } // rof
 462
 463       } // rof
 464
 465     } // rof
 466
 467   } // rof
 468   this->m_NumberOfBranches = actual_label;
 469
 470   // Label bifurcations
 471   typename TUniqueVertices::ConstIterator bifIt =
 472     bifurcations->Begin( );
 473   for( ; bifIt != bifurcations->End( ); ++bifIt )
 474   {
 475     actual_label++;
 476     double d = std::sqrt( double( input->GetPixel( *bifIt ) ) );
 477     _TRegion region = this->_Region( *bifIt, d * 1.5 );
 478     itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 479     itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
 480     iIt.GoToBegin( );
 481     lIt.GoToBegin( );
 482     for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
 483     {
 484       // Mask real input image
 485       if( !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) )
 486         lIt.Set( actual_label );
 487
 488     } // rof
 489
 490   } // rof
 491 }
 492
 493 // -------------------------------------------------------------------------
 494 template< class I, class O >
 495 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 496 _TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >::
 497 _Region( const TVertex& c, const double& r )
 498 {
 499   typename I::ConstPointer input = this->GetInput( );
 500   typename I::SpacingType spac = input->GetSpacing( );
 501   _TRegion region = input->GetLargestPossibleRegion( );
 502   typename I::IndexType idx0 = region.GetIndex( );
 503   typename I::IndexType idx1 = idx0 + region.GetSize( );
 504
 505   // Compute region size and index
 506   typename I::IndexType i0, i1;
 507   _TSize size;
 508   for( unsigned int d = 0; d < I::ImageDimension; ++d )
 509   {
 510     // NOTE: 3 is a minimum neighborhood size
 511     long s =
 512       long( std::ceil( r / double( spac[ d ] ) ) ) +
 513       long( this->m_SafetyNeighborhoodSize );
 514     i0[ d ] = c[ d ] - s;
 515     i1[ d ] = c[ d ] + s;
 516
 517     if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
 518     if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
 519     if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
 520     if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
 521     size[ d ] = i1[ d ] - i0[ d ];
 522
 523   }  // rof
 524
 525   // Prepare region and return it
 526   region.SetIndex( i0 );
 527   region.SetSize( size );
 528
 529   return( region );
 530 }
 531
 532 // -------------------------------------------------------------------------
 533 template< class I, class O >
 534 template< class _T >
 535 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 536 TVertex fpa::Image::DijkstraWithEndPointDetection< I, O >::
 537 _MaxInRegion( const _T* image, const TVertex& v, const double& r )
 538 {
 539   typedef itk::ImageRegionConstIteratorWithIndex< _T > _TIt;
 540
 541   _TIt iIt( image, this->_Region( v, r ) );
 542   iIt.GoToBegin( );
 543   TVertex max_vertex = iIt.GetIndex( );
 544   typename _T::PixelType max_value = iIt.Get( );
 545   for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 546   {
 547     typename _T::PixelType value = iIt.Get( );
 548     if( value > max_value )
 549     {
 550       max_value = value;
 551       max_vertex = iIt.GetIndex( );
 552
 553     } // fi
 554
 555   } // rof
 556   return( max_vertex );
 557 }
 558
 559 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
 560
 561 // eof - $RCSfile$