finding shortest path (revisited)

[Wolfram Moebius]

Hello all,

I am trying to build on some earlier code to find the shortest path after having computed the travel time (https://github.com/wmoebius/optpath).

I came to the conclusion that local interpolation is better suited than the current approach of generating an 'interpolation object' using scipy.

This is for two reasons: (1) It will be easier to integrate into the module. (2) Anecdotally, I had problems with interpolation creating local minima in arrival times in which path finding got stuck. I hope that more careful numerics can avoid that problem by using what we know about the the matrix of travel times (e.g., that there are no local minima of travel times).

In his book, Sethian states: "We do so by using second order Heun's method to integrate the ordinary differential equation. Starting from an initial point P0, we take one step using Euler's method to compute a temporary value. We then evaluate [nabla] T at this point by bilinear interpolation from the values of the gradient at each corner of the cell which contains the particular. These gradient values at cell corners are found by straightforward upwind differences using neighboring values. Once the gradient is evaluated at the temporary point, Heun's method is used to construct an average gradient. This is used to advance the position of the point, and the process is repeated."

I have two questions:

(1) Does anybody know of an alternative description (but see below)?

(2) Does anybody have a better understanding on 'straightforward upwind differences'? I know the concept is part of  FMM, but I am, at the moment, not well versed in the different approximation schemes that exist.

Thanks,

Wolfram

---

For completness:

In the 1998 paper on triangulated surfaces, Kimmel and Sethian write: "We use second order
Huen’s integration method on the triangulated surface with a switch to a first order scheme at sonic points in the gradient. When integrating within a given triangle, its three neighboring
triangles support the computation and are used to interpolate T as a second order polynomial whose six coefficients are computed from those given T values at the vertices."

[Jason Furtney]

Hi Wolfram,

Can you share an example where your current approach is not working?

It looks like by default the interpolation uses 3rd order polynomials,
maybe it would work better with second order or even linear
interpolation?

Another thought is that the problem is a manifestation of the
propagation of first order error from the narrow-band initialization.
This experimental branch uses a higher order narrow-band
initialization and gives more accurate results.
https://github.com/scikit-fmm/scikit-fmm/tree/new-init

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 view this discussion on the web visit https://groups.google.com/d/msgid/scikit-fmm/CAH-rdbT61WfccX7%2B1Fr5151193d0nZg1KJwhTS6Z6yx3zMmnOw%40mail.gmail.com.

[wolfram...@gmail.com]

Hi Jason,

Apologies for the delay. Took me a while to find the time to get this into a state I can share.

The not so minimal example speed matrix is here: https://www.dropbox.com/s/w4aapqskr1d6ffr/optpathchallenge.npz?dl=0

> # load data

> data = np.load('optpathchallenge.npz')
> sm = data['arr_0']
> phi = data['arr_1']

> # FMM

> tt=skfmm.travel_time(phi,sm,periodic=(1,0))

> # tracing path

> #OptPath.optpath_scipyode(np.arange(sm.shape[1]), np.arange(sm.shape[0]), tt, phi, (10000,500), 2)
> OptPath.optpath_eulerforward(np.arange(sm.shape[1]), np.arange(sm.shape[0]), tt, phi, (10000,500), 2)

Linear interpolation doesn't work in the current implementation, definitely something to look into.

I was thinking of doing this from scratch without the SciPy routines, thus my original message about whether you or anybody else knows whether one strategically only includes a subset of points for interpolation.

I don't think the initialisation error is related to this.

Best,

Wolfram