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( 3 );
  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_FinalTree.clear( );
 103   this->m_EndPoints.clear( );
 104   this->m_BifurcationPoints.clear( );
 105   if( this->m_Candidates.size( ) == 0 )
 106     return;
 107   this->InvokeEvent( TEndEvent( ) );
 108
 109   // Get some input values
 110   const I* input = this->GetInput( );
 111   typename I::SpacingType spac = input->GetSpacing( );
 112   double max_spac = spac[ 0 ];
 113   for( unsigned int d = 1; d < I::ImageDimension; ++d )
 114     max_spac =
 115       ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
 116   max_spac *= double( 3 );
 117
 118   // Prepare an object to hold marks
 119   typename TMarkImage::Pointer marks = this->GetOutputMarkImage( );
 120   marks->FillBuffer( 0 );
 121
 122   // Iterate over the candidates, starting from the candidates that
 123   // are near thin branches
 124   typename _TCandidates::const_reverse_iterator cIt =
 125     this->m_Candidates.rbegin( );
 126   this->m_NumberOfBranches = 0;
 127   for( ; cIt != this->m_Candidates.rend( ); ++cIt )
 128   {
 129     // If pixel hasn't been visited, continue
 130     TVertex v = cIt->second;
 131     if( marks->GetPixel( v ) != 0 )
 132       continue;
 133
 134     // Compute nearest start candidate
 135     _TRegion region = this->_Region( v, max_spac );
 136     itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
 137     iIt.GoToBegin( );
 138     TVertex max_vertex = iIt.GetIndex( );
 139     _TPixel max_value = iIt.Get( );
 140     for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
 141     {
 142       _TPixel value = iIt.Get( );
 143       if( value > max_value )
 144       {
 145         max_value = value;
 146         max_vertex = iIt.GetIndex( );
 147
 148       } // fi
 149
 150     } // rof
 151
 152     // Keep endpoint
 153     if( marks->GetPixel( max_vertex ) != 0 )
 154       continue;
 155     this->m_NumberOfBranches++;
 156     this->m_EndPoints.push_back( max_vertex );
 157
 158     // Construct path (at least the part that hasn't been iterated)
 159     do
 160     {
 161       if(
 162         std::find(
 163           this->m_BifurcationPoints.begin( ),
 164           this->m_BifurcationPoints.end( ),
 165           max_vertex
 166           ) == this->m_BifurcationPoints.end( )
 167         )
 168       {
 169         // Mark a sphere around current point as visited
 170         double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
 171         region = this->_Region( max_vertex, dist * double( 1.1 ) );
 172         itk::ImageRegionIteratorWithIndex< TMarkImage >
 173           mIt( marks, region );
 174         for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 175           mIt.Set( this->m_NumberOfBranches );
 176
 177         // Next vertex in current path
 178         this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
 179         this->m_FinalTree[ max_vertex ] = this->_Parent( max_vertex );
 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       } // fi
 188       max_vertex = this->_Parent( max_vertex );
 189
 190     } while( max_vertex != this->_Parent( max_vertex ) );
 191
 192     /* TODO
 193        bool terminate = false;
 194        do
 195        {
 196        if( this->m_FinalTree.find( max_vertex ) == this->m_FinalTree.end( ) )
 197        {
 198        // Mark a sphere around current point as visited
 199        double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
 200        region = this->_Region( max_vertex, dist * double( 1.25 ) );
 201        itk::ImageRegionIteratorWithIndex< TMarkImage >
 202        mIt( marks, region );
 203        for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 204        mIt.Set( true );
 205
 206        // Next vertex in current path
 207        this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
 208        this->m_FinalTree[ max_vertex ] = this->_Parent( max_vertex );
 209        }
 210        else
 211        {
 212        // A bifurcation point has been reached!
 213        this->m_BifurcationPoints.push_back( max_vertex );
 214        terminate = true;
 215
 216        } // fi
 217        max_vertex = this->_Parent( max_vertex );
 218
 219        } while( max_vertex != this->_Parent( max_vertex ) && !terminate );
 220
 221        if( !terminate )
 222        {
 223        this->m_FinalTree[ max_vertex ] = max_vertex;
 224        this->InvokeEvent( TEndBacktrackingEvent( this->m_NumberOfBranches ) );
 225
 226        } // fi
 227     */
 228
 229   } // rof
 230 }
 231
 232 // -------------------------------------------------------------------------
 233 template< class I, class C >
 234 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 235 _UpdateNeigh( _TNode& nn, const _TNode& n )
 236 {
 237   C nc = this->_Cost( nn.Vertex, n.Vertex );
 238   if( TCost( 0 ) < nc )
 239   {
 240     nc += TCost( 1 );
 241     nn.Cost = n.Cost + ( TCost( 1 ) / std::pow( nc, 4 ) );
 242     nn.Result = nn.Cost;
 243     return( true );
 244   }
 245   else
 246     return( false );
 247 }
 248
 249 // -------------------------------------------------------------------------
 250 template< class I, class C >
 251 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 252 _UpdateResult( _TNode& n )
 253 {
 254   bool ret = this->Superclass::_UpdateResult( n );
 255   if( ret )
 256   {
 257     TCost nc = this->_Cost( n.Vertex, n.Parent );
 258     if( TCost( 0 ) < nc )
 259     {
 260       TCost cc = n.Cost / nc;
 261       this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
 262       /* TODO: debug code
 263          this->GetOutput( )->SetPixel( n.Vertex, cc );
 264       */
 265     } // fi
 266
 267   } // fi
 268   return( ret );
 269 }
 270
 271 // -------------------------------------------------------------------------
 272 template< class I, class C >
 273 typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 274 _TRegion fpa::Image::DijkstraWithSphereBacktracking< I, C >::
 275 _Region( const TVertex& c, const double& r )
 276 {
 277   typename I::ConstPointer input = this->GetInput( );
 278   typename I::SpacingType spac = input->GetSpacing( );
 279   _TRegion region = input->GetLargestPossibleRegion( );
 280   typename I::IndexType idx0 = region.GetIndex( );
 281   typename I::IndexType idx1 = idx0 + region.GetSize( );
 282
 283   // Compute region size and index
 284   typename I::IndexType i0, i1;
 285   _TSize size;
 286   for( unsigned int d = 0; d < I::ImageDimension; ++d )
 287   {
 288     long s = long( std::ceil( r / double( spac[ d ] ) ) );
 289     i0[ d ] = c[ d ] - s;
 290     i1[ d ] = c[ d ] + s;
 291
 292     if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
 293     if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
 294     if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
 295     if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
 296     size[ d ] = i1[ d ] - i0[ d ];
 297
 298   }  // rof
 299
 300   // Prepare region and return it
 301   region.SetIndex( i0 );
 302   region.SetSize( size );
 303   return( region );
 304 }
 305
 306 #endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__
 307
 308 // eof - $RCSfile$