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 fpa::Image::DijkstraWithSphereBacktracking< I, C >::
  11 DijkstraWithSphereBacktracking( )
  12   : Superclass( )
  13 {
  14 }
  15
  16 // -------------------------------------------------------------------------
  17 template< class I, class C >
  18 fpa::Image::DijkstraWithSphereBacktracking< I, C >::
  19 ~DijkstraWithSphereBacktracking( )
  20 {
  21 }
  22
  23 // -------------------------------------------------------------------------
  24 template< class I, class C >
  25 void fpa::Image::DijkstraWithSphereBacktracking< I, C >::
  26 _BeforeMainLoop( )
  27 {
  28   const I* img = this->GetInput( );
  29   typename I::SpacingType spac = img->GetSpacing( );
  30   double max_spac = spac[ 0 ];
  31   for( unsigned int d = 1; d < I::ImageDimension; ++d )
  32     max_spac =
  33       ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
  34   max_spac *= double( 3 );
  35
  36   // Correct seeds
  37   for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
  38   {
  39     _TNode seed = this->m_Seeds[ s ];
  40     _TRegion region = this->_Region( seed.Vertex, max_spac );
  41     itk::ImageRegionConstIteratorWithIndex< I > iIt( img, region );
  42
  43     iIt.GoToBegin( );
  44     _TPixel max_value = iIt.Get( );
  45     for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
  46     {
  47       if( iIt.Get( ) > max_value )
  48       {
  49         seed.Vertex = iIt.GetIndex( );
  50         seed.Parent = seed.Vertex;
  51         max_value = iIt.Get( );
  52
  53       } // fi
  54
  55     } // rof
  56     this->m_Seeds[ s ] = seed;
  57
  58   } // rof
  59
  60   // End initialization
  61   this->Superclass::_BeforeMainLoop( );
  62   this->m_Candidates.clear( );
  63 }
  64
  65 // -------------------------------------------------------------------------
  66 template< class I, class C >
  67 void fpa::Image::DijkstraWithSphereBacktracking< I, C >::
  68 _AfterMainLoop( )
  69 {
  70   typedef itk::Image< bool, I::ImageDimension > _TMarkImage;
  71
  72   // Finish base algorithm
  73   this->Superclass::_AfterMainLoop( );
  74   this->m_FinalTree.clear( );
  75   this->m_EndPoints.clear( );
  76   if( this->m_Candidates.size( ) == 0 )
  77     return;
  78   this->InvokeEvent( TEndEvent( ) );
  79
  80   // Get input values
  81   const I* input = this->GetInput( );
  82   typename I::SpacingType spac = input->GetSpacing( );
  83   double max_spac = spac[ 0 ];
  84   for( unsigned int d = 1; d < I::ImageDimension; ++d )
  85     max_spac =
  86       ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
  87   max_spac *= double( 3 );
  88
  89   // Prepare an object to hold marks
  90   typename _TMarkImage::Pointer marks = _TMarkImage::New( );
  91   marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
  92   marks->SetRequestedRegion( input->GetRequestedRegion( ) );
  93   marks->SetBufferedRegion( input->GetBufferedRegion( ) );
  94   marks->SetOrigin( input->GetOrigin( ) );
  95   marks->SetSpacing( input->GetSpacing( ) );
  96   marks->SetDirection( input->GetDirection( ) );
  97   marks->Allocate( );
  98   marks->FillBuffer( false );
  99
 100   // Iterate over the candidates, starting from the candidates that
 101   // are near thin branches
 102   typename _TCandidates::const_reverse_iterator cIt =
 103     this->m_Candidates.rbegin( );
 104   for( int backId = 0; cIt != this->m_Candidates.rend( ); ++cIt )
 105   {
 106     // If pixel hasn't been visited, continue
 107     TVertex v = cIt->second;
 108     if( marks->GetPixel( v ) )
 109       continue;
 110
 111     // Compute nearest start candidate
 112     _TRegion region = this->_Region( v, max_spac );
 113     itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 114     iIt.GoToBegin( );
 115     TVertex max_vertex = iIt.GetIndex( );
 116     _TPixel max_value = iIt.Get( );
 117     for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 118     {
 119       _TPixel value = iIt.Get( );
 120       if( value > max_value )
 121       {
 122         max_value = value;
 123         max_vertex = iIt.GetIndex( );
 124
 125       } // fi
 126
 127     } // rof
 128
 129     // Keep endpoint
 130     if( marks->GetPixel( max_vertex ) )
 131       continue;
 132     backId++;
 133     this->m_EndPoints.push_back( max_vertex );
 134
 135     // Construct path (at least the part that hasn't been iterated)
 136     do
 137     {
 138       if( this->m_FinalTree.find( max_vertex ) == this->m_FinalTree.end( ) )
 139       {
 140         // Mark a sphere around current point as visited
 141         double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
 142         region = this->_Region( max_vertex, dist * double( 1.25 ) );
 143         itk::ImageRegionIteratorWithIndex< _TMarkImage >
 144           mIt( marks, region );
 145         for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 146           mIt.Set( true );
 147
 148         // Next vertex in current path
 149         this->InvokeEvent( TBacktrackingEvent( max_vertex, backId ) );
 150         this->m_FinalTree[ max_vertex ] = this->_Parent( max_vertex );
 151
 152       } // fi
 153       max_vertex = this->_Parent( max_vertex );
 154
 155     } while( max_vertex != this->_Parent( max_vertex ) );
 156     this->m_FinalTree[ max_vertex ] = max_vertex;
 157     this->InvokeEvent( TEndBacktrackingEvent( backId ) );
 158
 159   } // rof
 160 }
 161
 162 // -------------------------------------------------------------------------
 163 template< class I, class C >
 164 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 165 _UpdateNeigh( _TNode& nn, const _TNode& n )
 166 {
 167   C nc = this->_Cost( nn.Vertex, n.Vertex );
 168   if( TCost( 0 ) < nc )
 169   {
 170     nn.Cost = n.Cost + ( TCost( 1 ) / nc );
 171     nn.Result = nn.Cost;
 172     return( true );
 173   }
 174   else
 175     return( false );
 176 }
 177
 178 // -------------------------------------------------------------------------
 179 template< class I, class C >
 180 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 181 _UpdateResult( _TNode& n )
 182 {
 183   bool ret = this->Superclass::_UpdateResult( n );
 184   if( ret )
 185   {
 186     TCost nc = this->_Cost( n.Vertex, n.Parent );
 187     if( TCost( 0 ) < nc )
 188     {
 189       TCost cc = n.Cost / nc;
 190       this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
 191       /* TODO: debug code
 192          this->GetOutput( )->SetPixel( n.Vertex, cc );
 193       */
 194     } // fi
 195
 196   } // fi
 197   return( ret );
 198 }
 199
 200 // -------------------------------------------------------------------------
 201 template< class I, class C >
 202 typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 203 _TRegion fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 204 _Region( const TVertex& c, const double& r )
 205 {
 206   typename I::ConstPointer input = this->GetInput( );
 207   typename I::SpacingType spac = input->GetSpacing( );
 208   _TRegion region = input->GetLargestPossibleRegion( );
 209   typename I::IndexType idx0 = region.GetIndex( );
 210   typename I::IndexType idx1 = idx0 + region.GetSize( );
 211
 212   // Compute region size and index
 213   typename I::IndexType i0, i1;
 214   _TSize size;
 215   for( unsigned int d = 0; d < I::ImageDimension; ++d )
 216   {
 217     long s = long( std::ceil( r / double( spac[ d ] ) ) );
 218     i0[ d ] = c[ d ] - s;
 219     i1[ d ] = c[ d ] + s;
 220
 221     if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
 222     if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
 223     if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
 224     if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
 225     size[ d ] = i1[ d ] - i0[ d ];
 226
 227   }  // rof
 228
 229   // Prepare region and return it
 230   region.SetIndex( i0 );
 231   region.SetSize( size );
 232   return( region );
 233 }
 234
 235 #endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__
 236
 237 // eof - $RCSfile$