1 #ifndef __FPA__BASE__FASTMARCHING__HXX__
2 #define __FPA__BASE__FASTMARCHING__HXX__
7 #include <itkExceptionObject.h>
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( );
14 // -------------------------------------------------------------------------
15 template< class V, class C, class VV, class VC, class B >
16 fpa::Base::FastMarching< V, C, VV, VC, B >::
22 // -------------------------------------------------------------------------
23 template< class V, class C, class VV, class VC, class B >
24 fpa::Base::FastMarching< V, C, VV, VC, B >::
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 )
35 if( !( this->_IsMarked( nn.Vertex ) ) )
38 double F = double( this->_Value( nn.Vertex ) );
39 r = this->_SolveEikonal( re, nn, F );
42 nn.Result = TResult( re );
43 nn.Cost = TCost( re );
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 )
57 this->_Neighs( n, neighs );
59 // Compute interest nodes (i.e. the smallest valued in each dimension)
62 while( i < neighs.size( ) )
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 );
75 _TNode( nm, rm, n.FrontId ):
76 _TNode( np, rp, n.FrontId )
79 minNodes.push_back( _TNode( np, rp, n.FrontId ) );
81 minNodes.push_back( _TNode( nm, rm, n.FrontId ) );
84 std::sort( minNodes.begin( ), minNodes.end( ), _TNodeEikonalSort( ) );
87 double al = double( 0 );
88 double be = double( 0 );
89 double ga = double( -1 ) / ( F * F );
90 re = Self::LargestValue;
92 typename _TNodes::const_iterator nIt = minNodes.begin( );
94 while( nIt != minNodes.end( ) && valid )
96 if( !( re < double( nIt->Result ) ) )
98 double d = this->_Norm( n.Vertex, nIt->Vertex );
99 double s = double( 1 ) / ( d * d );
100 double v = double( nIt->Result );
106 double de = ( be * be ) - ( al * ga );
107 valid = !( de < double( 0 ) );
109 re = ( be + std::sqrt( de ) ) / al;
113 nIt = minNodes.end( );
119 #endif // __FPA__BASE__FASTMARCHING__HXX__