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