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   // Finish base algorithm
 109   this->Superclass::_AfterGenerateData( );
 110
 111   // Prepare backtracking objects
 112   this->m_EndPoints.clear( );
 113   this->m_BifurcationPoints.clear( );
 114   if( this->m_Candidates.size( ) == 0 )
 115     return;
 116   this->InvokeEvent( TEndEvent( ) );
 117   this->InvokeEvent( TStartBacktrackingEvent( ) );
 118
 119   // Get some input values
 120   const I* input = this->GetInput( );
 121   typename I::SpacingType spac = input->GetSpacing( );
 122   double ms = double( spac[ 0 ] );
 123   for( unsigned int d = 1; d < I::ImageDimension; ++d )
 124     ms =( ms < double( spac[ d ] ) )? double( spac[ d ] ): ms;
 125
 126   // Prepare labels
 127   TLabelImage* label = this->GetLabelImage( );
 128   label->FillBuffer( 0 );
 129
 130   // Object to hold marks
 131   std::set< TVertex, TVertexCompare > tree_marks;
 132
 133   // Object to hold all branches
 134   this->m_Branches.clear( );
 135
 136   // First label
 137   this->m_NumberOfBranches = 1;
 138
 139   // Iterate over the candidates, starting from the candidates that
 140   // are near thin branches
 141   typename _TCandidates::const_reverse_iterator cIt =
 142     this->m_Candidates.rbegin( );
 143   for( ; cIt != this->m_Candidates.rend( ); ++cIt )
 144   {
 145     // If pixel has been already labelled, pass
 146     TVertex v = cIt->second;
 147     if( label->GetPixel( v ) != 0 )
 148       continue;
 149
 150     // Compute nearest start candidate
 151     TVertex endpoint =
 152       this->_MaxInRegion(
 153         input, v,
 154         double( std::sqrt( input->GetPixel( v ) ) ) * double( 1.5 )
 155         );
 156
 157     // Re-check labelling
 158     if( this->_Node( endpoint ).Label == Self::FarLabel )
 159       continue;
 160     if( label->GetPixel( endpoint ) != 0 )
 161       continue;
 162     if( this->m_EndPoints.find( endpoint ) != this->m_EndPoints.end( ) )
 163       continue;
 164     this->m_EndPoints.insert( endpoint );
 165     std::cout << "endpoint " << endpoint << " inserted " << this->m_EndPoints.size( ) << std::endl;
 166
 167     // Get the path all the way to seed
 168     std::vector< TVertex > path;
 169     this->GetMinimumSpanningTree( )->
 170       GetPath( path, endpoint, this->GetSeed( 0 ) );
 171
 172     // Mark branches
 173     bool start = true;
 174     bool change = false;
 175     TVertex last_start = endpoint;
 176     typename std::vector< TVertex >::const_iterator pIt = path.begin( );
 177     for( ; pIt != path.end( ); ++pIt )
 178     {
 179       this->InvokeEvent(
 180         TBacktrackingEvent( *pIt, this->m_NumberOfBranches )
 181         );
 182       if( start )
 183       {
 184         if( tree_marks.find( *pIt ) == tree_marks.end( ) )
 185         {
 186           // Mark a region around current point as visited
 187           tree_marks.insert( *pIt );
 188           this->_Label( *pIt, TLabel( this->m_NumberOfBranches ) );
 189         }
 190         else
 191         {
 192           // A bifurcation point has been reached!
 193           if( *pIt != this->GetSeed( 0 ) )
 194           {
 195             this->m_Branches[ *pIt ][ last_start ] = this->m_NumberOfBranches;
 196             last_start = *pIt;
 197             this->m_BifurcationPoints.insert( *pIt );
 198             std::cout << "bifurcation = " << *pIt << " " << this->GetSeed( 0 ) << std::endl;
 199             this->m_NumberOfBranches++;
 200             start = false;
 201
 202           } // fi
 203
 204         } // fi
 205       }
 206       else
 207       {
 208         if(
 209           std::find(
 210             this->m_BifurcationPoints.begin( ),
 211             this->m_BifurcationPoints.end( ),
 212             *pIt
 213             ) == this->m_BifurcationPoints.end( )
 214           )
 215         {
 216           // Mark a sphere around current point as visited
 217           this->_Label( *pIt, TLabel( this->m_NumberOfBranches ) );
 218           change = true;
 219         }
 220         else
 221         {
 222           // A bifurcation point has been reached!
 223           if( *pIt != this->GetSeed( 0 ) )
 224           {
 225             this->m_Branches[ *pIt ][ last_start ] = this->m_NumberOfBranches;
 226             last_start = *pIt;
 227             this->m_NumberOfBranches++;
 228
 229           } // fi
 230
 231         } // fi
 232
 233       } // fi
 234
 235     } // rof
 236     if( start || change )
 237       this->m_NumberOfBranches++;
 238     this->InvokeEvent( TEndBacktrackingEvent( ) );
 239
 240   } // rof
 241
 242   typename TBranches::const_iterator bIt = this->m_Branches.begin( );
 243   unsigned int leo = 0;
 244   for( ; bIt != this->m_Branches.end( ); ++bIt )
 245   {
 246     typename TBranch::const_iterator brIt = bIt->second.begin( );
 247     for( ; brIt != bIt->second.end( ); ++brIt )
 248     {
 249       std::cout << bIt->first << " " << brIt->first << std::endl;
 250       leo++;
 251     }
 252   }
 253
 254   // Re-enumerate labels
 255   std::map< TLabel, unsigned long > histo;
 256   for(
 257     typename std::set< TVertex, TVertexCompare >::iterator treeIt = tree_marks.begin( );
 258     treeIt != tree_marks.end( );
 259     ++treeIt
 260     )
 261     histo[ label->GetPixel( *treeIt ) ]++;
 262
 263   std::map< TLabel, TLabel > changes;
 264   TLabel last_change = 1;
 265   for( TLabel i = 1; i <= this->m_NumberOfBranches; ++i )
 266   {
 267     if( histo[ i ] != 0 )
 268       changes[ i ] = last_change++;
 269
 270   } // rof
 271   this->m_NumberOfBranches = changes.size( );
 272
 273   std::cout << leo << " - " << this->m_EndPoints.size( ) << " - " << this->m_BifurcationPoints.size( ) << " - " << this->m_NumberOfBranches << std::endl;
 274
 275   /*
 276     for(
 277     typename TTree::iterator treeIt = this->m_FullTree.begin( );
 278     treeIt != this->m_FullTree.end( );
 279     ++treeIt
 280     )
 281     {
 282     TMark old = treeIt->second.second;
 283     treeIt->second.second = changes[ old ];
 284
 285     } // fi
 286
 287     // Construct reduced tree
 288     for(
 289     typename TVertices::const_iterator eIt = this->m_EndPoints.begin( );
 290     eIt != this->m_EndPoints.end( );
 291     ++eIt
 292     )
 293     {
 294     typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt );
 295     if( tIt != this->m_FullTree.end( ) )
 296     {
 297     TMark id = tIt->second.second;
 298     do
 299     {
 300     tIt = this->m_FullTree.find( tIt->second.first );
 301     if( tIt == this->m_FullTree.end( ) )
 302     break;
 303
 304     } while( tIt->second.second == id );
 305     this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id );
 306
 307     } // fi
 308
 309     } // rof
 310
 311     for(
 312     typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( );
 313     bIt != this->m_BifurcationPoints.end( );
 314     ++bIt
 315     )
 316     {
 317     typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt );
 318     if( tIt != this->m_FullTree.end( ) )
 319     {
 320     TMark id = tIt->second.second;
 321     do
 322     {
 323     tIt = this->m_FullTree.find( tIt->second.first );
 324     if( tIt == this->m_FullTree.end( ) )
 325     break;
 326
 327     } while( tIt->second.second == id );
 328     this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id );
 329
 330     } // fi
 331
 332     } // rof
 333   */
 334 }
 335
 336 // -------------------------------------------------------------------------
 337 template< class I, class O >
 338 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 339 _SetResult( const TVertex& v, const _TNode& n )
 340 {
 341   this->Superclass::_SetResult( v, n );
 342   this->m_Candidates.insert( _TCandidate( n.Result, v ) );
 343 }
 344
 345 // -------------------------------------------------------------------------
 346 template< class I, class O >
 347 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 348 _TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >::
 349 _Region( const TVertex& c, const double& r )
 350 {
 351   typename I::ConstPointer input = this->GetInput( );
 352   typename I::SpacingType spac = input->GetSpacing( );
 353   _TRegion region = input->GetLargestPossibleRegion( );
 354   typename I::IndexType idx0 = region.GetIndex( );
 355   typename I::IndexType idx1 = idx0 + region.GetSize( );
 356
 357   // Compute region size and index
 358   typename I::IndexType i0, i1;
 359   _TSize size;
 360   for( unsigned int d = 0; d < I::ImageDimension; ++d )
 361   {
 362     long s = long( std::ceil( r / double( spac[ d ] ) ) ) + 3;
 363     i0[ d ] = c[ d ] - s;
 364     i1[ d ] = c[ d ] + s;
 365
 366     if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
 367     if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
 368     if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
 369     if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
 370     size[ d ] = i1[ d ] - i0[ d ];
 371
 372   }  // rof
 373
 374   // Prepare region and return it
 375   region.SetIndex( i0 );
 376   region.SetSize( size );
 377
 378   return( region );
 379 }
 380
 381 // -------------------------------------------------------------------------
 382 template< class I, class O >
 383 template< class _T >
 384 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
 385 TVertex fpa::Image::DijkstraWithEndPointDetection< I, O >::
 386 _MaxInRegion( const _T* image, const TVertex& v, const double& r )
 387 {
 388   typedef itk::ImageRegionConstIteratorWithIndex< _T > _TIt;
 389
 390   _TIt iIt( image, this->_Region( v, r ) );
 391   iIt.GoToBegin( );
 392   TVertex max_vertex = iIt.GetIndex( );
 393   typename _T::PixelType max_value = iIt.Get( );
 394   for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 395   {
 396     typename _T::PixelType value = iIt.Get( );
 397     if( value > max_value )
 398     {
 399       max_value = value;
 400       max_vertex = iIt.GetIndex( );
 401
 402     } // fi
 403
 404   } // rof
 405   return( max_vertex );
 406 }
 407
 408 // -------------------------------------------------------------------------
 409 template< class I, class O >
 410 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
 411 _Label( const TVertex& v, const TLabel& l )
 412 {
 413   typedef itk::ImageRegionIteratorWithIndex< TLabelImage > _TIt;
 414
 415   double d = std::sqrt( double( this->GetInput( )->GetPixel( v ) ) );
 416   _TRegion region = this->_Region( v, d );
 417   if( region.GetNumberOfPixels( ) > 0 )
 418   {
 419     _TIt lIt( this->GetLabelImage( ), region );
 420     for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
 421       lIt.Set( l );
 422   }
 423   else
 424     this->GetLabelImage( )->SetPixel( v, l );
 425 }
 426
 427 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
 428
 429 // eof - $RCSfile$