lib/fpa/Base/FastMarching.hxx

   1 #ifndef __FPA__BASE__FASTMARCHING__HXX__
   2 #define __FPA__BASE__FASTMARCHING__HXX__
   3
   4 #include <cmath>
   5 #include <limits>
   6
   7 #include <itkExceptionObject.h>
   8
   9 // ---------------------------------------------------------------------
  10 template< class V, class C, class VV, class VC, class B >
  11 const double fpa::Base::FastMarching< V, C, VV, VC, B >::
  12 LargestValue = std::numeric_limits< double >::max( );
  13
  14 // -------------------------------------------------------------------------
  15 template< class V, class C, class VV, class VC, class B >
  16 fpa::Base::FastMarching< V, C, VV, VC, B >::
  17 FastMarching( )
  18   : Superclass( )
  19 {
  20 }
  21
  22 // -------------------------------------------------------------------------
  23 template< class V, class C, class VV, class VC, class B >
  24 fpa::Base::FastMarching< V, C, VV, VC, B >::
  25 ~FastMarching( )
  26 {
  27 }
  28
  29 // -------------------------------------------------------------------------
  30 template< class V, class C, class VV, class VC, class B >
  31 bool fpa::Base::FastMarching< V, C, VV, VC, B >::
  32 _UpdateNeigh( _TNode& nn, const _TNode& n )
  33 {
  34   bool r = true;
  35   if( !( this->_IsMarked( nn.Vertex ) ) )
  36   {
  37     double re;
  38     double F = double( this->_Value( nn.Vertex ) );
  39     r = this->_SolveEikonal( re, nn, F );
  40     if( r )
  41     {
  42       nn.Result = TResult( re );
  43       nn.Cost = TCost( re );
  44
  45     } // fi
  46
  47   } // fi
  48   return( r );
  49 }
  50
  51 // -------------------------------------------------------------------------
  52 template< class V, class C, class VV, class VC, class B >
  53 bool fpa::Base::FastMarching< V, C, VV, VC, B >::
  54 _SolveEikonal( double& re, const _TNode& n, const double& F )
  55 {
  56   _TNodes neighs;
  57   this->_Neighs( n, neighs );
  58
  59   // Compute interest nodes (i.e. the smallest valued in each dimension)
  60   _TNodes minNodes;
  61   unsigned int i = 0;
  62   while( i < neighs.size( ) )
  63   {
  64     V nm = neighs[ i++ ].Vertex;
  65     V np = neighs[ i++ ].Vertex;
  66     bool mm = this->_IsMarked( nm );
  67     bool mp = this->_IsMarked( np );
  68     TResult rm = this->_Result( nm );
  69     TResult rp = this->_Result( np );
  70
  71     _TNode Ti;
  72     if( mm && mp )
  73       minNodes.push_back(
  74         ( rm < rp )?
  75         _TNode( nm, rm, n.FrontId ):
  76         _TNode( np, rp, n.FrontId )
  77         );
  78     else if( !mm && mp )
  79       minNodes.push_back( _TNode( np, rp, n.FrontId ) );
  80     else if( mm && !mp )
  81       minNodes.push_back( _TNode( nm, rm, n.FrontId ) );
  82
  83   } // elihw
  84   std::sort( minNodes.begin( ), minNodes.end( ), _TNodeEikonalSort( ) );
  85
  86   // Compute solution
  87   double al = double( 0 );
  88   double be = double( 0 );
  89   double ga = double( -1 ) / ( F * F );
  90   re = Self::LargestValue;
  91
  92   typename _TNodes::const_iterator nIt = minNodes.begin( );
  93   bool valid = true;
  94   while( nIt != minNodes.end( ) && valid )
  95   {
  96     if( !( re < double( nIt->Result ) ) )
  97     {
  98       double d = this->_Norm( n.Vertex, nIt->Vertex );
  99       double s = double( 1 ) / ( d * d );
 100       double v = double( nIt->Result );
 101
 102       al += s;
 103       be += v * s;
 104       ga += v * v * s;
 105
 106       double de = ( be * be ) - ( al * ga );
 107       valid = !( de < double( 0 ) );
 108       if( valid )
 109         re = ( be + std::sqrt( de ) ) / al;
 110       nIt++;
 111     }
 112     else
 113       nIt = minNodes.end( );
 114
 115   } // elihw
 116   return( valid );
 117 }
 118
 119 #endif // __FPA__BASE__FASTMARCHING__HXX__
 120
 121 // eof - $RCSfile$