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;
 340     mst->GetPath( path, v, seed );
 341
 342     // Backtracking to find endpoints and bifurcations
 343     bool adding_new_points = true;
 344     typename TVertices::const_iterator pIt = path.begin( );
 345     for( ; pIt != path.end( ) && adding_new_points; ++pIt )
 346     {
 347       this->InvokeEvent( TBacktrackingEvent( *pIt, 1 ) );
 348       if( tree_marks.find( *pIt ) == tree_marks.end( ) )
 349       {
 350         // Mark current point as a tree point
 351         tree_marks.insert( *pIt );
 352
 353         // Mark a region around current point as visited
 354         vr = std::sqrt( double( input->GetPixel( *pIt ) ) );
 355         _TRegion region = this->_Region( *pIt, vr );
 356         if( region.GetNumberOfPixels( ) > 0 )
 357         {
 358           itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 359           for( iIt.GoToBegin( ); !iIt.IsAtEnd( ); ++iIt )
 360             pixel_marks[ iIt.GetIndex( ) ] = true;
 361         }
 362         else
 363           pixel_marks[ *pIt ] = true;
 364       }
 365       else
 366       {
 367         // A bifurcation point has been reached!
 368         if( *pIt != seed )
 369         {
 370           bifurcations->Insert( *pIt );
 371           adding_new_points = false;
 372
 373         } // fi
 374
 375       } // fi
 376
 377     } // rof
 378     this->InvokeEvent( TEndBacktrackingEvent( ) );
 379
 380   } // rof
 381 }
 382
 383 // -------------------------------------------------------------------------
 384 template< class I, class O >
 385 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 386 _FindBranches( )
 387 {
 388   // Configure and obtain inputs
 389   TVertex seed = this->GetSeed( 0 );
 390   const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
 391   const TUniqueVertices* endpoints = this->GetEndPoints( );
 392   const TUniqueVertices* bifurcations = this->GetBifurcations( );
 393   TBranches* branches = this->GetBranches( );
 394   branches->Clear( );
 395
 396   // Reconstruct pixels
 397   typename TUniqueVertices::ConstIterator eIt = endpoints->Begin( );
 398   for( ; eIt != endpoints->End( ); ++eIt )
 399   {
 400     // Get the path all the way to global seed
 401     TVertices path;
 402     mst->GetPath( path, *eIt, seed );
 403
 404     TVertex start_vertex = *eIt;
 405     typename TVertices::const_iterator pIt = path.begin( );
 406     for( ; pIt != path.end( ); ++pIt )
 407     {
 408       if( bifurcations->Find( *pIt ) )
 409       {
 410         branches->SetValue( start_vertex, *pIt, 1 );
 411         start_vertex = *pIt;
 412
 413       } // fi
 414
 415     } // rof
 416
 417     // Finish with branches to global seed
 418     branches->SetValue( start_vertex, seed, 1 );
 419
 420   } // rof
 421 }
 422
 423 // -------------------------------------------------------------------------
 424 template< class I, class O >
 425 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 426 _LabelAll( )
 427 {
 428   // Configure and obtain inputs
 429   const I* input = this->GetInput( );
 430   const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
 431   TBranches* branches = this->GetBranches( );
 432   TLabelImage* label = this->GetLabelImage( );
 433   label->FillBuffer( 0 );
 434
 435   // Label branches
 436   typename TBranches::Iterator bIt = branches->Begin( );
 437   TLabel actual_label = 0;
 438   for( ; bIt != branches->End( ); ++bIt )
 439   {
 440     typename TBranches::RowIterator brIt = branches->Begin( bIt );
 441     for( ; brIt != branches->End( bIt ); ++brIt )
 442     {
 443       actual_label++;
 444       brIt->second = actual_label;
 445
 446       TVertices path;
 447       mst->GetPath( path, bIt->first, brIt->first );
 448       typename TVertices::const_iterator pIt = path.begin( );
 449       for( ; pIt != path.end( ); ++pIt )
 450       {
 451         double d = std::sqrt( double( input->GetPixel( *pIt ) ) );
 452         _TRegion region = this->_Region( *pIt, d );
 453         itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 454         itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
 455         iIt.GoToBegin( );
 456         lIt.GoToBegin( );
 457         for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
 458         {
 459           // Mask real input image
 460           if(
 461             !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) &&
 462             lIt.Get( ) == TLabel( 0 )
 463             )
 464             lIt.Set( actual_label );
 465
 466         } // rof
 467
 468       } // rof
 469
 470     } // rof
 471
 472   } // rof
 473   this->m_NumberOfBranches = actual_label;
 474
 475   // Label bifurcations
 476   /* TODO: Is this necessary?
 477      typename TUniqueVertices::const_iterator bifIt =
 478      this->m_Bifurcations.begin( );
 479      for( ; bifIt != this->m_Bifurcations.end( ); ++bifIt )
 480      {
 481      actual_label++;
 482      double d = std::sqrt( double( input->GetPixel( *bifIt ) ) );
 483      _TRegion region = this->_Region( *bifIt, d );
 484      itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 485      itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
 486      iIt.GoToBegin( );
 487      lIt.GoToBegin( );
 488      for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
 489      {
 490      // Mask real input image
 491      if( !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) )
 492      lIt.Set( actual_label );
 493
 494      } // rof
 495
 496      } // rof
 497   */
 498 }
 499
 500 // -------------------------------------------------------------------------
 501 template< class I, class O >
 502 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 503 _TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >::
 504 _Region( const TVertex& c, const double& r )
 505 {
 506   typename I::ConstPointer input = this->GetInput( );
 507   typename I::SpacingType spac = input->GetSpacing( );
 508   _TRegion region = input->GetLargestPossibleRegion( );
 509   typename I::IndexType idx0 = region.GetIndex( );
 510   typename I::IndexType idx1 = idx0 + region.GetSize( );
 511
 512   // Compute region size and index
 513   typename I::IndexType i0, i1;
 514   _TSize size;
 515   for( unsigned int d = 0; d < I::ImageDimension; ++d )
 516   {
 517     // NOTE: 3 is a minimum neighborhood size
 518     long s =
 519       long( std::ceil( r / double( spac[ d ] ) ) ) +
 520       long( this->m_SafetyNeighborhoodSize );
 521     i0[ d ] = c[ d ] - s;
 522     i1[ d ] = c[ d ] + s;
 523
 524     if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
 525     if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
 526     if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
 527     if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
 528     size[ d ] = i1[ d ] - i0[ d ];
 529
 530   }  // rof
 531
 532   // Prepare region and return it
 533   region.SetIndex( i0 );
 534   region.SetSize( size );
 535
 536   return( region );
 537 }
 538
 539 // -------------------------------------------------------------------------
 540 template< class I, class O >
 541 template< class _T >
 542 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 543 TVertex fpa::Image::DijkstraWithEndPointDetection< I, O >::
 544 _MaxInRegion( const _T* image, const TVertex& v, const double& r )
 545 {
 546   typedef itk::ImageRegionConstIteratorWithIndex< _T > _TIt;
 547
 548   _TIt iIt( image, this->_Region( v, r ) );
 549   iIt.GoToBegin( );
 550   TVertex max_vertex = iIt.GetIndex( );
 551   typename _T::PixelType max_value = iIt.Get( );
 552   for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 553   {
 554     typename _T::PixelType value = iIt.Get( );
 555     if( value > max_value )
 556     {
 557       max_value = value;
 558       max_vertex = iIt.GetIndex( );
 559
 560     } // fi
 561
 562   } // rof
 563   return( max_vertex );
 564 }
 565
 566 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
 567
 568 // eof - $RCSfile$