lib/fpa/Image/DijkstraWithSphereBacktracking.hxx

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