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     this->m_EndPoints.push_back( max_vertex );
 157
 158     // Construct path (at least the part that hasn't been iterated)
 159     bool start = true;
 160     bool change = false;
 161     do
 162     {
 163       if( start )
 164       {
 165         if( this->m_FullTree.find( max_vertex ) == this->m_FullTree.end( ) )
 166         {
 167           // Mark a sphere around current point as visited
 168           double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
 169           region = this->_Region( max_vertex, dist * double( 1.5 ) );
 170           itk::ImageRegionIteratorWithIndex< TMarkImage >
 171             mIt( marks, region );
 172           for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 173             mIt.Set( this->m_NumberOfBranches );
 174
 175           // Next vertex in current path
 176           this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
 177           this->m_FullTree[ max_vertex ] =
 178             TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
 179           // std::cout << "New: " << this->m_NumberOfBranches << std::endl;
 180         }
 181         else
 182         {
 183           // A bifurcation point has been reached!
 184           this->m_BifurcationPoints.push_back( max_vertex );
 185           this->m_NumberOfBranches++;
 186
 187           this->m_FullTree[ max_vertex ] =
 188             TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
 189           // std::cout << "Bifurcation: " << this->m_NumberOfBranches << std::endl;
 190
 191           start = false;
 192
 193         } // fi
 194       }
 195       else
 196       {
 197         if(
 198           std::find(
 199             this->m_BifurcationPoints.begin( ),
 200             this->m_BifurcationPoints.end( ),
 201             max_vertex
 202             ) == this->m_BifurcationPoints.end( )
 203           )
 204         {
 205           // Mark a sphere around current point as visited
 206           double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
 207           region = this->_Region( max_vertex, dist * double( 1.5 ) );
 208           itk::ImageRegionIteratorWithIndex< TMarkImage >
 209             mIt( marks, region );
 210           for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 211             mIt.Set( this->m_NumberOfBranches );
 212
 213           // Next vertex in current path
 214           this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
 215           this->m_FullTree[ max_vertex ] =
 216             TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
 217           change = true;
 218           // std::cout << "Change: " << this->m_NumberOfBranches << std::endl;
 219         }
 220         else
 221         {
 222           // A bifurcation point has been reached!
 223           // TODO: this->m_BifurcationPoints.push_back( max_vertex );
 224           this->m_NumberOfBranches++;
 225           this->m_FullTree[ max_vertex ] =
 226             TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
 227           // std::cout << "Change bifurcation: " << this->m_NumberOfBranches << std::endl;
 228
 229         } // fi
 230
 231       } // fi
 232       max_vertex = this->_Parent( max_vertex );
 233
 234     } while( max_vertex != this->_Parent( max_vertex ) );
 235     if( start || change )
 236       this->m_NumberOfBranches++;
 237     this->m_FullTree[ max_vertex ] = TTreeNode( max_vertex, this->m_NumberOfBranches );
 238
 239     this->InvokeEvent( TEndBacktrackingEvent( ) );
 240
 241   } // rof
 242
 243   std::map< TMark, unsigned long > histo;
 244   for(
 245     typename TTree::iterator treeIt = this->m_FullTree.begin( );
 246     treeIt != this->m_FullTree.end( );
 247     ++treeIt
 248     )
 249     histo[ treeIt->second.second ]++;
 250
 251   std::map< TMark, TMark > changes;
 252   TMark last_change = 1;
 253   for( TMark i = 1; i <= this->m_NumberOfBranches; ++i )
 254   {
 255     if( histo[ i ] != 0 )
 256       changes[ i ] = last_change++;
 257
 258   } // rof
 259   this->m_NumberOfBranches = changes.size( );
 260
 261   for(
 262     typename TTree::iterator treeIt = this->m_FullTree.begin( );
 263     treeIt != this->m_FullTree.end( );
 264     ++treeIt
 265     )
 266   {
 267     TMark old = treeIt->second.second;
 268     treeIt->second.second = changes[ old ];
 269
 270   } // fi
 271
 272   // Construct reduced tree
 273   for(
 274     typename TVertices::const_iterator eIt = this->m_EndPoints.begin( );
 275     eIt != this->m_EndPoints.end( );
 276     ++eIt
 277     )
 278   {
 279     typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt );
 280     if( tIt != this->m_FullTree.end( ) )
 281     {
 282       TMark id = tIt->second.second;
 283       do
 284       {
 285         tIt = this->m_FullTree.find( tIt->second.first );
 286         if( tIt == this->m_FullTree.end( ) )
 287           break;
 288
 289       } while( tIt->second.second == id );
 290       this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id );
 291
 292     } // fi
 293
 294   } // rof
 295
 296   for(
 297     typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( );
 298     bIt != this->m_BifurcationPoints.end( );
 299     ++bIt
 300     )
 301   {
 302     typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt );
 303     if( tIt != this->m_FullTree.end( ) )
 304     {
 305       TMark id = tIt->second.second;
 306       do
 307       {
 308         tIt = this->m_FullTree.find( tIt->second.first );
 309         if( tIt == this->m_FullTree.end( ) )
 310           break;
 311
 312       } while( tIt->second.second == id );
 313       this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id );
 314
 315     } // fi
 316
 317   } // rof
 318 }
 319
 320 // -------------------------------------------------------------------------
 321 template< class I, class C >
 322 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 323 _UpdateNeigh( _TNode& nn, const _TNode& n )
 324 {
 325   C nc = this->_Cost( nn.Vertex, n.Vertex );
 326   if( TCost( 0 ) < nc )
 327   {
 328     typename I::PointType pnn, pn;
 329     this->GetInput( )->TransformIndexToPhysicalPoint( nn.Vertex, pnn );
 330     this->GetInput( )->TransformIndexToPhysicalPoint( n.Vertex, pn );
 331
 332
 333     nc += TCost( 1 );
 334     nn.Cost = n.Cost + ( TCost( pnn.EuclideanDistanceTo( pn ) ) / std::pow( nc, 4 ) );
 335     nn.Result = nn.Cost;
 336     return( true );
 337   }
 338   else
 339     return( false );
 340 }
 341
 342 // -------------------------------------------------------------------------
 343 template< class I, class C >
 344 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 345 _UpdateResult( _TNode& n )
 346 {
 347   bool ret = this->Superclass::_UpdateResult( n );
 348   if( ret )
 349   {
 350     TCost nc = this->_Cost( n.Vertex, n.Parent );
 351     if( TCost( 0 ) < nc )
 352     {
 353       TCost cc = n.Cost / nc;
 354       this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
 355       /* TODO: debug code
 356          this->GetOutput( )->SetPixel( n.Vertex, cc );
 357       */
 358     } // fi
 359
 360   } // fi
 361   return( ret );
 362 }
 363
 364 // -------------------------------------------------------------------------
 365 template< class I, class C >
 366 typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 367 _TRegion fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 368 _Region( const TVertex& c, const double& r )
 369 {
 370   typename I::ConstPointer input = this->GetInput( );
 371   typename I::SpacingType spac = input->GetSpacing( );
 372   _TRegion region = input->GetLargestPossibleRegion( );
 373   typename I::IndexType idx0 = region.GetIndex( );
 374   typename I::IndexType idx1 = idx0 + region.GetSize( );
 375
 376   // Compute region size and index
 377   typename I::IndexType i0, i1;
 378   _TSize size;
 379   for( unsigned int d = 0; d < I::ImageDimension; ++d )
 380   {
 381     long s = long( std::ceil( r / double( spac[ d ] ) ) );
 382     i0[ d ] = c[ d ] - s;
 383     i1[ d ] = c[ d ] + s;
 384
 385     if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
 386     if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
 387     if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
 388     if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
 389     size[ d ] = i1[ d ] - i0[ d ];
 390
 391   }  // rof
 392
 393   // Prepare region and return it
 394   region.SetIndex( i0 );
 395   region.SetSize( size );
 396   return( region );
 397 }
 398
 399 #endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__
 400
 401 // eof - $RCSfile$