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   this->m_BifurcationPoints.clear( );
  77   if( this->m_Candidates.size( ) == 0 )
  78     return;
  79   this->InvokeEvent( TEndEvent( ) );
  80
  81   // Get some input values
  82   const I* input = this->GetInput( );
  83   typename I::SpacingType spac = input->GetSpacing( );
  84   double max_spac = spac[ 0 ];
  85   for( unsigned int d = 1; d < I::ImageDimension; ++d )
  86     max_spac =
  87       ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
  88   max_spac *= double( 3 );
  89
  90   // Prepare an object to hold marks
  91   typename _TMarkImage::Pointer marks = _TMarkImage::New( );
  92   marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
  93   marks->SetRequestedRegion( input->GetRequestedRegion( ) );
  94   marks->SetBufferedRegion( input->GetBufferedRegion( ) );
  95   marks->SetOrigin( input->GetOrigin( ) );
  96   marks->SetSpacing( input->GetSpacing( ) );
  97   marks->SetDirection( input->GetDirection( ) );
  98   marks->Allocate( );
  99   marks->FillBuffer( false );
 100
 101   // Iterate over the candidates, starting from the candidates that
 102   // are near thin branches
 103   typename _TCandidates::const_reverse_iterator cIt =
 104     this->m_Candidates.rbegin( );
 105   for( int backId = 0; cIt != this->m_Candidates.rend( ); ++cIt )
 106   {
 107     // If pixel hasn't been visited, continue
 108     TVertex v = cIt->second;
 109     if( marks->GetPixel( v ) )
 110       continue;
 111
 112     // Compute nearest start candidate
 113     _TRegion region = this->_Region( v, max_spac );
 114     itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 115     iIt.GoToBegin( );
 116     TVertex max_vertex = iIt.GetIndex( );
 117     _TPixel max_value = iIt.Get( );
 118     for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 119     {
 120       _TPixel value = iIt.Get( );
 121       if( value > max_value )
 122       {
 123         max_value = value;
 124         max_vertex = iIt.GetIndex( );
 125
 126       } // fi
 127
 128     } // rof
 129
 130     // Keep endpoint
 131     if( marks->GetPixel( max_vertex ) )
 132       continue;
 133     backId++;
 134     this->m_EndPoints.push_back( max_vertex );
 135
 136     // Construct path (at least the part that hasn't been iterated)
 137     bool terminate = false;
 138     do
 139     {
 140       if( this->m_FinalTree.find( max_vertex ) == this->m_FinalTree.end( ) )
 141       {
 142         // Mark a sphere around current point as visited
 143         double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
 144         region = this->_Region( max_vertex, dist * double( 1.25 ) );
 145         itk::ImageRegionIteratorWithIndex< _TMarkImage >
 146           mIt( marks, region );
 147         for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 148           mIt.Set( true );
 149
 150         // Next vertex in current path
 151         this->InvokeEvent( TBacktrackingEvent( max_vertex, backId ) );
 152         this->m_FinalTree[ max_vertex ] = this->_Parent( max_vertex );
 153       }
 154       else
 155       {
 156         // A bifurcation point has been reached!
 157         this->m_BifurcationPoints.push_back( max_vertex );
 158         terminate = true;
 159
 160       } // fi
 161       max_vertex = this->_Parent( max_vertex );
 162
 163     } while( max_vertex != this->_Parent( max_vertex ) && !terminate );
 164
 165     if( !terminate )
 166     {
 167       this->m_FinalTree[ max_vertex ] = max_vertex;
 168       this->InvokeEvent( TEndBacktrackingEvent( backId ) );
 169
 170     } // fi
 171
 172   } // rof
 173 }
 174
 175 // -------------------------------------------------------------------------
 176 template< class I, class C >
 177 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 178 _UpdateNeigh( _TNode& nn, const _TNode& n )
 179 {
 180   C nc = this->_Cost( nn.Vertex, n.Vertex );
 181   if( TCost( 0 ) < nc )
 182   {
 183     nc += TCost( 1 );
 184     nn.Cost = n.Cost + ( TCost( 1 ) / std::pow( nc, 4 ) );
 185     nn.Result = nn.Cost;
 186     return( true );
 187   }
 188   else
 189     return( false );
 190 }
 191
 192 // -------------------------------------------------------------------------
 193 template< class I, class C >
 194 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 195 _UpdateResult( _TNode& n )
 196 {
 197   bool ret = this->Superclass::_UpdateResult( n );
 198   if( ret )
 199   {
 200     TCost nc = this->_Cost( n.Vertex, n.Parent );
 201     if( TCost( 0 ) < nc )
 202     {
 203       TCost cc = n.Cost / nc;
 204       this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
 205       /* TODO: debug code
 206          this->GetOutput( )->SetPixel( n.Vertex, cc );
 207       */
 208     } // fi
 209
 210   } // fi
 211   return( ret );
 212 }
 213
 214 // -------------------------------------------------------------------------
 215 template< class I, class C >
 216 typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 217 _TRegion fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 218 _Region( const TVertex& c, const double& r )
 219 {
 220   typename I::ConstPointer input = this->GetInput( );
 221   typename I::SpacingType spac = input->GetSpacing( );
 222   _TRegion region = input->GetLargestPossibleRegion( );
 223   typename I::IndexType idx0 = region.GetIndex( );
 224   typename I::IndexType idx1 = idx0 + region.GetSize( );
 225
 226   // Compute region size and index
 227   typename I::IndexType i0, i1;
 228   _TSize size;
 229   for( unsigned int d = 0; d < I::ImageDimension; ++d )
 230   {
 231     long s = long( std::ceil( r / double( spac[ d ] ) ) );
 232     i0[ d ] = c[ d ] - s;
 233     i1[ d ] = c[ d ] + s;
 234
 235     if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
 236     if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
 237     if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
 238     if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
 239     size[ d ] = i1[ d ] - i0[ d ];
 240
 241   }  // rof
 242
 243   // Prepare region and return it
 244   region.SetIndex( i0 );
 245   region.SetSize( size );
 246   return( region );
 247 }
 248
 249 #endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__
 250
 251 // eof - $RCSfile$