lib/fpa/Image/DijkstraWithEndPointDetection.hxx

   1 #ifndef __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
   2 #define __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
   3
   4 #include <vector>
   5 #include <itkImageRegionConstIteratorWithIndex.h>
   6
   7 // -------------------------------------------------------------------------
   8 template< class I, class O >
   9 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
  10 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
  11 GetLabelImage( )
  12 {
  13   return(
  14     dynamic_cast< TLabelImage* >(
  15       this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
  16       )
  17     );
  18 }
  19
  20 // -------------------------------------------------------------------------
  21 template< class I, class O >
  22 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
  23 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
  24 GetLabelImage( ) const
  25 {
  26   return(
  27     dynamic_cast< const TLabelImage* >(
  28       this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
  29       )
  30     );
  31 }
  32
  33 // -------------------------------------------------------------------------
  34 template< class I, class O >
  35 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
  36 GraftLabelImage( itk::DataObject* obj )
  37 {
  38   TLabelImage* lbl =
  39     dynamic_cast< TLabelImage* >(
  40       this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
  41       );
  42   if( lbl != NULL )
  43     this->GraftNthOutput( this->m_LabelImageIndex, lbl );
  44 }
  45
  46 // -------------------------------------------------------------------------
  47 template< class I, class O >
  48 fpa::Image::DijkstraWithEndPointDetection< I, O >::
  49 DijkstraWithEndPointDetection( )
  50   : Superclass( )
  51 {
  52   this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( );
  53   this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 1 );
  54   this->itk::ProcessObject::SetNthOutput(
  55     this->m_LabelImageIndex, TLabelImage::New( )
  56     );
  57 }
  58
  59 // -------------------------------------------------------------------------
  60 template< class I, class O >
  61 fpa::Image::DijkstraWithEndPointDetection< I, O >::
  62 ~DijkstraWithEndPointDetection( )
  63 {
  64 }
  65
  66 // -------------------------------------------------------------------------
  67 template< class I, class O >
  68 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
  69 _BeforeGenerateData( )
  70 {
  71   // Correct seeds
  72   /* TODO
  73      for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
  74      {
  75      _TNode seed = this->m_Seeds[ s ];
  76      _TRegion region = this->_Region( seed.Vertex, max_spac );
  77      itk::ImageRegionConstIteratorWithIndex< I > iIt( img, region );
  78
  79      iIt.GoToBegin( );
  80      _TPixel max_value = iIt.Get( );
  81      for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
  82      {
  83      if( iIt.Get( ) > max_value )
  84      {
  85      seed.Vertex = iIt.GetIndex( );
  86      seed.Parent = seed.Vertex;
  87      max_value = iIt.Get( );
  88
  89      } // fi
  90
  91      } // rof
  92      this->m_Seeds[ s ] = seed;
  93
  94      } // rof
  95   */
  96
  97   // End initialization
  98   this->Superclass::_BeforeGenerateData( );
  99   this->m_Candidates.clear( );
 100 }
 101
 102 // -------------------------------------------------------------------------
 103 template< class I, class O >
 104 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 105 _AfterGenerateData( )
 106 {
 107   this->Superclass::_AfterGenerateData( );
 108
 109   // Finish base algorithm
 110   /* TODO
 111      this->m_FullTree.clear( );
 112      this->m_ReducedTree.clear( );
 113   */
 114   this->m_EndPoints.clear( );
 115   this->m_BifurcationPoints.clear( );
 116   if( this->m_Candidates.size( ) == 0 )
 117     return;
 118   this->InvokeEvent( TEndEvent( ) );
 119   this->InvokeEvent( TStartBacktrackingEvent( ) );
 120
 121   // Get some input values
 122   const I* input = this->GetInput( );
 123   typename I::SpacingType spac = input->GetSpacing( );
 124   double max_spac = spac[ 0 ];
 125   for( unsigned int d = 1; d < I::ImageDimension; ++d )
 126     max_spac =
 127       ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
 128   max_spac *= double( 3 );
 129   TLabelImage* label = this->GetLabelImage( );
 130   label->FillBuffer( 0 );
 131
 132   // Prepare an object to hold marks
 133   /* TODO
 134      typename TMarkImage::Pointer marks = this->GetOutputMarkImage( );
 135      marks->FillBuffer( 0 );
 136   */
 137
 138   // Iterate over the candidates, starting from the candidates that
 139   // are near thin branches
 140   typename _TCandidates::const_reverse_iterator cIt =
 141     this->m_Candidates.rbegin( );
 142   this->m_NumberOfBranches = 1;
 143   std::map< TLabel, std::pair< TVertex, TVertex > > branches;
 144   for( ; cIt != this->m_Candidates.rend( ); ++cIt )
 145   {
 146     // If pixel hasn't been visited, continue
 147     TVertex v = cIt->second;
 148     if( label->GetPixel( v ) != 0 )
 149       continue;
 150
 151     // Compute nearest start candidate
 152     _TRegion region = this->_Region( v, max_spac );
 153     itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 154     iIt.GoToBegin( );
 155     TVertex max_vertex = iIt.GetIndex( );
 156     _TPixel max_value = iIt.Get( );
 157     for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 158     {
 159       _TPixel value = iIt.Get( );
 160       if( value > max_value )
 161       {
 162         max_value = value;
 163         max_vertex = iIt.GetIndex( );
 164
 165       } // fi
 166
 167     } // rof
 168
 169     // Keep endpoint
 170     if( label->GetPixel( max_vertex ) != 0 )
 171       continue;
 172     this->m_EndPoints.push_back( max_vertex );
 173
 174     // Get the path all the way to seed
 175     std::vector< TVertex > path;
 176     this->GetMinimumSpanningTree( )->
 177       GetPath( path, max_vertex, this->GetSeed( 0 ) );
 178
 179     // Mark branches
 180     bool start = true;
 181     bool change = false;
 182     TVertex last_start = max_vertex;
 183     typename std::vector< TVertex >::const_iterator pIt = path.begin( );
 184     for( ; pIt != path.end( ); ++pIt )
 185     {
 186       this->InvokeEvent(
 187         TBacktrackingEvent( *pIt, this->m_NumberOfBranches )
 188         );
 189       if( start )
 190       {
 191         TLabel lbl = label->GetPixel( *pIt );
 192         if( lbl == 0 || lbl == this->m_NumberOfBranches )
 193         {
 194           // Mark a sphere around current point as visited
 195           double dist = std::sqrt( double( input->GetPixel( *pIt ) ) );
 196           region = this->_Region( max_vertex, dist * double( 1.1 ) );
 197           itk::ImageRegionIteratorWithIndex< TLabelImage >
 198             lIt( label, region );
 199           for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
 200           {
 201             if( lIt.Get( ) == 0 )
 202               lIt.Set( this->m_NumberOfBranches );
 203
 204           } // rof
 205
 206           // Next vertex in current path
 207           // TODO: this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
 208           /*
 209             this->m_FullTree[ max_vertex ] =
 210             TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
 211             std::cout << "New: " << this->m_NumberOfBranches << std::endl;
 212           */
 213         }
 214         else
 215         {
 216           // A bifurcation point has been reached!
 217           branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt );
 218           std::cout << "bif " << this->m_NumberOfBranches << std::endl;
 219           last_start = *pIt;
 220           this->m_BifurcationPoints.push_back( *pIt );
 221           this->m_NumberOfBranches++;
 222           start = false;
 223
 224         } // fi
 225       }
 226       else
 227       {
 228         if(
 229           std::find(
 230             this->m_BifurcationPoints.begin( ),
 231             this->m_BifurcationPoints.end( ),
 232             *pIt
 233             ) == this->m_BifurcationPoints.end( )
 234           )
 235         {
 236           // Mark a sphere around current point as visited
 237           double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
 238           region = this->_Region( max_vertex, dist * double( 1.1 ) );
 239           itk::ImageRegionIteratorWithIndex< TLabelImage >
 240             lIt( label, region );
 241           for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
 242             lIt.Set( this->m_NumberOfBranches );
 243
 244           // Next vertex in current path
 245           /* TODO
 246              this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
 247              this->m_FullTree[ max_vertex ] =
 248              TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
 249           */
 250           change = true;
 251           // std::cout << "Change: " << this->m_NumberOfBranches << std::endl;
 252         }
 253         else
 254         {
 255           // A bifurcation point has been reached!
 256           // TODO: this->m_BifurcationPoints.push_back( max_vertex );
 257           branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt );
 258           std::cout << "change " << this->m_NumberOfBranches << std::endl;
 259
 260           last_start = *pIt;
 261           this->m_NumberOfBranches++;
 262           /* TODO
 263              this->m_FullTree[ max_vertex ] =
 264              TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
 265           */
 266           // std::cout << "Change bifurcation: " << this->m_NumberOfBranches << std::endl;
 267
 268         } // fi
 269
 270       } // fi
 271
 272     } // rof
 273     if( start || change )
 274       this->m_NumberOfBranches++;
 275     this->InvokeEvent( TEndBacktrackingEvent( ) );
 276
 277
 278
 279     /*
 280       this->InvokeEvent( TEndBacktrackingEvent( ) );
 281     */
 282
 283   } // rof
 284
 285   std::cout << this->m_NumberOfBranches << " " << branches.size( ) << std::endl;
 286   std::cout << this->m_EndPoints.size( ) << " " << this->m_BifurcationPoints.size( ) << std::endl;
 287
 288
 289   // Re-enumerate labels
 290   /*
 291   std::map< TLabel, unsigned long > histo;
 292   for(
 293     typename TTree::iterator treeIt = this->m_FullTree.begin( );
 294     treeIt != this->m_FullTree.end( );
 295     ++treeIt
 296     )
 297     histo[ treeIt->second.second ]++;
 298
 299   std::map< TMark, TMark > changes;
 300   TMark last_change = 1;
 301   for( TMark i = 1; i <= this->m_NumberOfBranches; ++i )
 302   {
 303     if( histo[ i ] != 0 )
 304       changes[ i ] = last_change++;
 305
 306   } // rof
 307   this->m_NumberOfBranches = changes.size( );
 308   */
 309
 310   /*
 311     for(
 312     typename TTree::iterator treeIt = this->m_FullTree.begin( );
 313     treeIt != this->m_FullTree.end( );
 314     ++treeIt
 315     )
 316     {
 317     TMark old = treeIt->second.second;
 318     treeIt->second.second = changes[ old ];
 319
 320     } // fi
 321
 322     // Construct reduced tree
 323     for(
 324     typename TVertices::const_iterator eIt = this->m_EndPoints.begin( );
 325     eIt != this->m_EndPoints.end( );
 326     ++eIt
 327     )
 328     {
 329     typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt );
 330     if( tIt != this->m_FullTree.end( ) )
 331     {
 332     TMark id = tIt->second.second;
 333     do
 334     {
 335     tIt = this->m_FullTree.find( tIt->second.first );
 336     if( tIt == this->m_FullTree.end( ) )
 337     break;
 338
 339     } while( tIt->second.second == id );
 340     this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id );
 341
 342     } // fi
 343
 344     } // rof
 345
 346     for(
 347     typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( );
 348     bIt != this->m_BifurcationPoints.end( );
 349     ++bIt
 350     )
 351     {
 352     typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt );
 353     if( tIt != this->m_FullTree.end( ) )
 354     {
 355     TMark id = tIt->second.second;
 356     do
 357     {
 358     tIt = this->m_FullTree.find( tIt->second.first );
 359     if( tIt == this->m_FullTree.end( ) )
 360     break;
 361
 362     } while( tIt->second.second == id );
 363     this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id );
 364
 365     } // fi
 366
 367     } // rof
 368   */
 369 }
 370
 371 // -------------------------------------------------------------------------
 372 template< class I, class O >
 373 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 374 _SetResult( const TVertex& v, const _TNode& n )
 375 {
 376   this->Superclass::_SetResult( v, n );
 377
 378   TResult vv = TResult( this->_VertexValue( v ) );
 379   if( TResult( 0 ) < vv )
 380     this->m_Candidates.insert( _TCandidate( n.Result / vv, v ) );
 381 }
 382
 383 // -------------------------------------------------------------------------
 384 template< class I, class O >
 385 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 386 _TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >::
 387 _Region( const TVertex& c, const double& r )
 388 {
 389   typename I::ConstPointer input = this->GetInput( );
 390   typename I::SpacingType spac = input->GetSpacing( );
 391   _TRegion region = input->GetLargestPossibleRegion( );
 392   typename I::IndexType idx0 = region.GetIndex( );
 393   typename I::IndexType idx1 = idx0 + region.GetSize( );
 394
 395   // Compute region size and index
 396   typename I::IndexType i0, i1;
 397   _TSize size;
 398   for( unsigned int d = 0; d < I::ImageDimension; ++d )
 399   {
 400     long s = long( std::ceil( r / double( spac[ d ] ) ) );
 401     i0[ d ] = c[ d ] - s;
 402     i1[ d ] = c[ d ] + s;
 403
 404     if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
 405     if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
 406     if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
 407     if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
 408     size[ d ] = i1[ d ] - i0[ d ];
 409
 410   }  // rof
 411
 412   // Prepare region and return it
 413   region.SetIndex( i0 );
 414   region.SetSize( size );
 415   return( region );
 416 }
 417
 418 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
 419
 420 // eof - $RCSfile$