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 {
 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   const I* input = this->GetInput( );
 207   typename I::SpacingType spac = input->GetSpacing( );
 208   double max_spac = double( spac[ 0 ] );
 209   for( unsigned int d = 1; d < I::ImageDimension; ++d )
 210     max_spac =
 211       ( double( spac[ d ] ) > max_spac )? double( spac[ d ] ): max_spac;
 212   max_spac *= double( 3 );
 213
 214   // Correct seeds
 215   for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s )
 216   {
 217     TVertex seed = this->m_SeedVertices[ s ];
 218     _TNode n = this->m_Seeds[ seed ];
 219     _TRegion region = this->_Region( seed, max_spac );
 220     itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 221
 222     iIt.GoToBegin( );
 223     _TPixel max_value = iIt.Get( );
 224     for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 225     {
 226       if( iIt.Get( ) > max_value )
 227       {
 228         this->m_SeedVertices[ s ] = iIt.GetIndex( );
 229         max_value = iIt.Get( );
 230
 231       } // fi
 232
 233     } // rof
 234     this->m_Seeds.erase( seed );
 235     n.Parent = this->m_SeedVertices[ s ];
 236     this->m_Seeds[ this->m_SeedVertices[ s ] ] = n;
 237
 238   } // rof
 239
 240   // End initialization
 241   this->Superclass::_BeforeGenerateData( );
 242   this->m_Candidates.clear( );
 243 }
 244
 245 // -------------------------------------------------------------------------
 246 template< class I, class O >
 247 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 248 _AfterGenerateData( )
 249 {
 250   // Finish base algorithm
 251   this->Superclass::_AfterGenerateData( );
 252
 253   // Check if backtracking has any sense
 254   if( this->m_Candidates.size( ) == 0 )
 255     return;
 256   this->InvokeEvent( TEndEvent( ) );
 257   this->InvokeEvent( TStartBacktrackingEvent( ) );
 258
 259   // Detect endpoints and bifurcations
 260   this->_EndPointsAndBifurcations( );
 261
 262   // Find branches
 263   this->_FindBranches( );
 264
 265   // Label pixels and branches
 266   this->_LabelAll( );
 267 }
 268
 269 // -------------------------------------------------------------------------
 270 template< class I, class O >
 271 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 272 _SetResult( const TVertex& v, const _TNode& n )
 273 {
 274   this->Superclass::_SetResult( v, n );
 275   this->m_Candidates.insert( _TCandidate( n.Result, v ) );
 276 }
 277
 278 // -------------------------------------------------------------------------
 279 template< class I, class O >
 280 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 281 _EndPointsAndBifurcations( )
 282 {
 283   // Prepare backtracking objects
 284   TUniqueVertices* endpoints = this->GetEndPoints( );
 285   TUniqueVertices* bifurcations = this->GetBifurcations( );
 286   endpoints->Clear( );
 287   bifurcations->Clear( );
 288
 289   // Get some input values
 290   TVertex seed = this->GetSeed( 0 );
 291   const I* input = this->GetInput( );
 292   const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
 293
 294   // Prepare a pixel marks structure
 295   typedef std::map< TVertex, bool, TVertexCompare > _TPixelMarks;
 296   _TPixelMarks pixel_marks;
 297
 298   // Object to hold tree-nodes marks
 299   typedef std::set< TVertex, TVertexCompare > _TTreeMarks;
 300   _TTreeMarks tree_marks;
 301
 302   // Iterate over the candidates, starting from the candidates that
 303   // are near thin branches (high accumulated cost)
 304   typename _TCandidates::const_reverse_iterator cIt =
 305     this->m_Candidates.rbegin( );
 306   for( ; cIt != this->m_Candidates.rend( ); ++cIt )
 307   {
 308     // If pixel has already been visited, pass
 309     TVertex v = cIt->second;
 310     if( pixel_marks[ v ] )
 311       continue;
 312
 313     // Correct it to nearest start candidate (high distance value)
 314     double vr = std::sqrt( double( input->GetPixel( v ) ) );
 315     v = this->_MaxInRegion( input, v, vr * double( 1.5 ) );
 316
 317     // Now, check for real marking conditions
 318     //   1. Has it been visited by dijkstra?
 319     if( this->_Node( v ).Label == Self::FarLabel )
 320       continue;
 321     //   2. Is it already marked?
 322     if( pixel_marks[ v ] )
 323       continue;
 324     //   3. Is it already an endpoint?
 325     if( endpoints->Find( v ) )
 326       continue;
 327     //   4. Ok, it is completely new!
 328     endpoints->Insert( v );
 329
 330     // Get the path all the way to global seed
 331     TVertices path;
 332     mst->GetPath( path, v, seed );
 333
 334     // Backtracking to find endpoints and bifurcations
 335     bool adding_new_points = true;
 336     typename TVertices::const_iterator pIt = path.begin( );
 337     for( ; pIt != path.end( ) && adding_new_points; ++pIt )
 338     {
 339       this->InvokeEvent( TBacktrackingEvent( *pIt, 1 ) );
 340       if( tree_marks.find( *pIt ) == tree_marks.end( ) )
 341       {
 342         // Mark current point as a tree point
 343         tree_marks.insert( *pIt );
 344
 345         // Mark a region around current point as visited
 346         vr = std::sqrt( double( input->GetPixel( *pIt ) ) );
 347         _TRegion region = this->_Region( *pIt, vr );
 348         if( region.GetNumberOfPixels( ) > 0 )
 349         {
 350           itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 351           for( iIt.GoToBegin( ); !iIt.IsAtEnd( ); ++iIt )
 352             pixel_marks[ iIt.GetIndex( ) ] = true;
 353         }
 354         else
 355           pixel_marks[ *pIt ] = true;
 356       }
 357       else
 358       {
 359         // A bifurcation point has been reached!
 360         if( *pIt != seed )
 361         {
 362           bifurcations->Insert( *pIt );
 363           adding_new_points = false;
 364
 365         } // fi
 366
 367       } // fi
 368
 369     } // rof
 370     this->InvokeEvent( TEndBacktrackingEvent( ) );
 371
 372   } // rof
 373 }
 374
 375 // -------------------------------------------------------------------------
 376 template< class I, class O >
 377 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 378 _FindBranches( )
 379 {
 380   // Configure and obtain inputs
 381   TVertex seed = this->GetSeed( 0 );
 382   const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
 383   const TUniqueVertices* endpoints = this->GetEndPoints( );
 384   const TUniqueVertices* bifurcations = this->GetBifurcations( );
 385   TBranches* branches = this->GetBranches( );
 386   branches->Clear( );
 387
 388   // Reconstruct pixels
 389   typename TUniqueVertices::ConstIterator eIt = endpoints->Begin( );
 390   for( ; eIt != endpoints->End( ); ++eIt )
 391   {
 392     // Get the path all the way to global seed
 393     TVertices path;
 394     mst->GetPath( path, *eIt, seed );
 395
 396     TVertex start_vertex = *eIt;
 397     typename TVertices::const_iterator pIt = path.begin( );
 398     for( ; pIt != path.end( ); ++pIt )
 399     {
 400       if( bifurcations->Find( *pIt ) )
 401       {
 402         branches->SetValue( start_vertex, *pIt, 1 );
 403         start_vertex = *pIt;
 404
 405       } // fi
 406
 407     } // rof
 408
 409     // Finish with branches to global seed
 410     branches->SetValue( start_vertex, seed, 1 );
 411
 412   } // rof
 413 }
 414
 415 // -------------------------------------------------------------------------
 416 template< class I, class O >
 417 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 418 _LabelAll( )
 419 {
 420   // Configure and obtain inputs
 421   const I* input = this->GetInput( );
 422   const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
 423   TBranches* branches = this->GetBranches( );
 424   TLabelImage* label = this->GetLabelImage( );
 425   label->FillBuffer( 0 );
 426
 427   // Label branches
 428   typename TBranches::Iterator bIt = branches->Begin( );
 429   TLabel actual_label = 0;
 430   for( ; bIt != branches->End( ); ++bIt )
 431   {
 432     typename TBranches::RowIterator brIt = branches->Begin( bIt );
 433     for( ; brIt != branches->End( bIt ); ++brIt )
 434     {
 435       actual_label++;
 436       brIt->second = actual_label;
 437
 438       TVertices path;
 439       mst->GetPath( path, bIt->first, brIt->first );
 440       typename TVertices::const_iterator pIt = path.begin( );
 441       for( ; pIt != path.end( ); ++pIt )
 442       {
 443         double d = std::sqrt( double( input->GetPixel( *pIt ) ) );
 444         _TRegion region = this->_Region( *pIt, d );
 445         itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 446         itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
 447         iIt.GoToBegin( );
 448         lIt.GoToBegin( );
 449         for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
 450         {
 451           // Mask real input image
 452           if(
 453             !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) &&
 454             lIt.Get( ) == TLabel( 0 )
 455             )
 456             lIt.Set( actual_label );
 457
 458         } // rof
 459
 460       } // rof
 461
 462     } // rof
 463
 464   } // rof
 465   this->m_NumberOfBranches = actual_label;
 466
 467   // Label bifurcations
 468   /* TODO: Is this necessary?
 469      typename TUniqueVertices::const_iterator bifIt =
 470      this->m_Bifurcations.begin( );
 471      for( ; bifIt != this->m_Bifurcations.end( ); ++bifIt )
 472      {
 473      actual_label++;
 474      double d = std::sqrt( double( input->GetPixel( *bifIt ) ) );
 475      _TRegion region = this->_Region( *bifIt, d );
 476      itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 477      itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
 478      iIt.GoToBegin( );
 479      lIt.GoToBegin( );
 480      for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
 481      {
 482      // Mask real input image
 483      if( !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) )
 484      lIt.Set( actual_label );
 485
 486      } // rof
 487
 488      } // rof
 489   */
 490 }
 491
 492 // -------------------------------------------------------------------------
 493 template< class I, class O >
 494 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 495 _TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >::
 496 _Region( const TVertex& c, const double& r )
 497 {
 498   typename I::ConstPointer input = this->GetInput( );
 499   typename I::SpacingType spac = input->GetSpacing( );
 500   _TRegion region = input->GetLargestPossibleRegion( );
 501   typename I::IndexType idx0 = region.GetIndex( );
 502   typename I::IndexType idx1 = idx0 + region.GetSize( );
 503
 504   // Compute region size and index
 505   typename I::IndexType i0, i1;
 506   _TSize size;
 507   for( unsigned int d = 0; d < I::ImageDimension; ++d )
 508   {
 509     // NOTE: 3 is a minimum neighborhood size
 510     long s = long( std::ceil( r / double( spac[ d ] ) ) ) + 3;
 511     i0[ d ] = c[ d ] - s;
 512     i1[ d ] = c[ d ] + s;
 513
 514     if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
 515     if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
 516     if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
 517     if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
 518     size[ d ] = i1[ d ] - i0[ d ];
 519
 520   }  // rof
 521
 522   // Prepare region and return it
 523   region.SetIndex( i0 );
 524   region.SetSize( size );
 525
 526   return( region );
 527 }
 528
 529 // -------------------------------------------------------------------------
 530 template< class I, class O >
 531 template< class _T >
 532 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 533 TVertex fpa::Image::DijkstraWithEndPointDetection< I, O >::
 534 _MaxInRegion( const _T* image, const TVertex& v, const double& r )
 535 {
 536   typedef itk::ImageRegionConstIteratorWithIndex< _T > _TIt;
 537
 538   _TIt iIt( image, this->_Region( v, r ) );
 539   iIt.GoToBegin( );
 540   TVertex max_vertex = iIt.GetIndex( );
 541   typename _T::PixelType max_value = iIt.Get( );
 542   for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 543   {
 544     typename _T::PixelType value = iIt.Get( );
 545     if( value > max_value )
 546     {
 547       max_value = value;
 548       max_vertex = iIt.GetIndex( );
 549
 550     } // fi
 551
 552   } // rof
 553   return( max_vertex );
 554 }
 555
 556 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
 557
 558 // eof - $RCSfile$