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