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