Another feature request

[Adrian Butscher]

Hi Jason,

Thanks again for implementing my "narrow band" feature request.  We've been using it for a few weeks now and it seems quite stable -- and it's very useful since our calculations are that much faster!  We'll be submitting a paper in two weeks which describes an algorithm that uses skfmm.  We'd like to acknowledge and/or cite your work.  I've looked on the web page (briefly) and haven't found instructions for doing so.  Do you have a preference for this?  

I'm also writing because I have another feature request, if possible.  For our purpose, we would like to have access to the vector field grad(F) where F is the computed SDF.  We've written our own first-order-accurate upwind gradient operator to do this, but it doesn't manage to achieve || grad(F) || = 1 everywhere and the output is rather noisy .  There are other upwind schemes out there, and I think your default is a second-order-accurate one.  It would be nice to get grad(F) from your code, since that would presumably minimize these issues.  Do you compute grad(F) explicitly in your code anywhere, or do you only work with the equation || grad(F) || = a(x) where a is the speed function?  If the former, would it be possible to expose to the user the ability to get grad(F)?

Thanks again for your help.

Adrian

[Jason Furtney]

Dear Adrian,

Thanks for the email and glad to hear you calculations are working
out. For acknowledgment please mention that you are using Python and
scikit-fmm.

Regarding the gradient of the signed distance function:

As points are updated the gradient is locally estimated via
second-order upwind derivatives. This gradient is not retained. As I
understand it, if you also want a second-order accurate gradient
additional work is required during the point update. The augmented
fast marching method was developed for this purpose. This method
marches out the signed distance function and its derivatives. See
http://arxiv.org/abs/0905.3409 and http://arxiv.org/abs/1111.6903

You could try using the numpy.gradient function to calculate the
gradient of the returned SDF but I think this will also be noisy.
Calculating the gradient from the scalar field is amplifying the
error.

In the current scikit-fmm implementation the calculation of the
distance to the points adjacent to the zero contour is only
first-order accurate. These first order errors are then propagated.
The first step in implementing the augmented method is to make this
initialization step second order. I started on a branch here:

https://github.com/scikit-fmm/scikit-fmm/tree/bicubic-init

to use cubic interpolation to reconstruct the zero contour and get a
second-order accurate distance to the initially frozen points. (This
was originally described by Chopp, 2001 "Some Improvements of the Fast
Marching Method")

I have this working in two-dimensions as a proof-of-concept but more
testing is needed. Have a look at the file skt.py in this branch for
an example. Currently, the bi-cubic interpolation is not correct for
cells on the outer boundary of the array. If the zero contour does not
intersect the edge of the domain this should not be problem. I have
the same thing partially working in 3D.

You could try this branch and see if the results are any better, give
the initorder=2 keyword argument to the distance function. If you have
an example of what you are trying to do I would be happy to have a
look. I am tied up this week but may have some time next week to have
another crack at the augmented method.

Thanks,
Jason

> --
> You received this message because you are subscribed to the Google Groups
> "scikit-fmm" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to scikit-fmm+...@googlegroups.com.
> To post to this group, send email to sciki...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/scikit-fmm/CABPx5hE_EGnDCx0T5TNQ_b2ZxABLAXT%3Dyc4GD%3Dc6_D6OWfwzRA%40mail.gmail.com.
> For more options, visit https://groups.google.com/d/optout.

[Jason Furtney]

Dear Adrian,

In this github branch:
https://github.com/scikit-fmm/scikit-fmm/tree/bicubic-init

I have the second-order frozen point initialization working in 2D. I
tested a small example and this change makes the magnitude of the
gradient of the SDF less noisy. See the attached image and
https://github.com/scikit-fmm/scikit-fmm/blob/bicubic-init/grad_mag_test.py

This branch is quite old unfortunately, I will rebase it or create a
new branch soon.

A few issues need to be addressed before releasing this feature in a
new version:

* This branch uses scipy.optimize.fsolve to do part of the
initialization. For practical reason, I do not want scikit-fmm to
depend on SciPy. To get around this I plan to implement some kind of
solver based on Newton's method to do what fsolve is currently doing.

* All of the new init code in this branch in written in Python so it
is probably pretty slow.

* More testing is needed.

* A little more effort is needed to get this working correctly in 3D.

It would be great if you could test this out on your problem.

Thanks,
Jason