Bug in distance computation ? about scikit-fmm HOT 1 CLOSED

Dear Abed,

Thanks for the message. When I run this I get a different answer:

array([[ 1.70710678,  1.        ,  1.70710678],
       [ 1.        ,  0.        ,  1.        ],
       [ 1.70710678,  1.        ,  1.70710678]])

What version are you using?

To address your question: I think this is an error inherent in the method. The good news is that the error is corrected as interface propagates outward. The maximum error in the method occurs along the diagonals and is greatest on the zero level set. Because of the way the gradient operator is discretized the diagonal points can only "see" the points to the N,S,E and W. They cannot "see" the central point. This is corrected as the interface marches out because points start to "see" the curvature of the interface via the N,S,E and W points.

The following shows that as the interface marches out the relative error decreases:

phi = np.ones((31, 31))
phi[15, 15] = 0
d = skfmm.distance(phi)
print d[-1,-1], d[0,0], d[-1,0], d[0,-1]
print (2*15.5**2 )**.5

21.4520926896 21.4520926896 21.4520926896 21.4520926896
21.9203102168

phi = np.ones((301, 301))
phi[150, 150] = 0
d = skfmm.distance(phi)
print d[-1,-1], d[0,0], d[-1,0], d[0,-1]
print (2*150.5**2 )**.5

212.348795302 212.348795302 212.348795302 212.348795302
212.839141137

There are some ways to improve this error. One method is to use a gradient stencil which accounts for the diagonal directions (See: M. Sabry Hassouna and Aly A. Farag, "Multistencils Fast Marching Methods: A Highly Accurate Solution to the Eikonal Equation on Cartesian Domains," IEEE Transaction on Pattern Analysis and Machine Intelligence PAMI, vol. 29, no. 9, pp. 1-12, Sep 2007.)

Another option is to reconstruct the zero level-set using multidimensional cubic interpolation to to better initialize the distance of the adjacent points. There is a new-init branch where we experimented with this approach: https://github.com/scikit-fmm/scikit-fmm/tree/new-init There is also some discussion of this on the mailing list https://groups.google.com/forum/#!topic/scikit-fmm/labxaSVAUqA

from scikit-fmm.