Re: questions on scikit-fmm

[Jason Furtney]

Dear Wolfram,

Thanks for the reply. What you describe should be possible.

The path finding method on the wiki reconstructs the arrival time
function using cubic interpolation so the paths curve on a sub-pixel
scale. As you point out this is relatively slow. The Dijkstra
algorithm gives discrete paths from node to node. If you want to
enumerate all the discrete paths you should be able to do this much
faster.

Something like this may be all you need:

import numpy as np
import pylab as plt
from skfmm import distance
phi = np.ones((15,15))
phi[7,7] = -1

d = distance(phi)
gx, gy = np.gradient(d, edge_order=2)

angle = np.arctan2(gx,gy)/2.0/np.pi + 1/16.0
angle[angle<0] +=1
parent_code = (angle*8).astype(int)
print parent_code

[[5 5 5 5 5 6 6 6 6 6 7 7 7 7 7]
[5 5 5 5 5 5 6 6 6 7 7 7 7 7 7]
[5 5 5 5 5 5 6 6 6 7 7 7 7 7 7]
[5 5 5 5 5 5 5 6 7 7 7 7 7 7 7]
[5 5 5 5 5 5 5 6 7 7 7 7 7 7 7]
[4 5 5 5 5 5 5 6 7 7 7 7 7 7 0]
[4 4 4 5 5 5 5 6 7 7 7 7 0 0 0]
[4 4 4 4 4 4 4 0 0 0 0 0 0 0 0]
[4 4 4 3 3 3 3 2 1 1 1 1 0 0 0]
[4 3 3 3 3 3 3 2 1 1 1 1 1 1 0]
[3 3 3 3 3 3 3 2 1 1 1 1 1 1 1]
[3 3 3 3 3 3 3 2 1 1 1 1 1 1 1]
[3 3 3 3 3 3 2 2 2 1 1 1 1 1 1]
[3 3 3 3 3 3 2 2 2 1 1 1 1 1 1]
[3 3 3 3 3 2 2 2 2 2 1 1 1 1 1]]

Here each point has a code which gives the direction to the nearest cell.

It should be possible to create a new marching algorithm which returns
the parent cell of each cell as you describe. Have a look at the
implementation of the extension velocity class. The virtual function
finalizePoint() is called when a point is frozen. In this function you
could assign the parent cell. The method does not track an explicit
parent cell but you could find one.

Thanks,
Jason




On Wed, Jun 15, 2016 at 4:57 AM, Wolfram Moebius <w.mo...@tue.nl> wrote:
> Dear Jason,
>
> Thanks for your super quick reply!
>
>> 1) A speed of 0 results in a singularity (divide by zero) in the point
>> update. In the constructor of travel_time_marcher.h there is a check
>> that the speed in a cell is greater than the machine epsilon. If it is
>> less it masks that cell. Maybe we should update the documentation to
>> make this more clear.
> OK, so, yes, I could just set the speed to 0, but it will result in masking the cell. The cleaner version therefore is to mask the cells in phi. Is that what you are saying?
>
>> 2) The method does not calculate individual paths but you can
>> calculate those from the travel time field. Here is an example:
>> https://github.com/scikit-fmm/scikit-fmm/wiki/optimal_path
>>
>> Lei Pan has recently been working on making this faster for
>> calculating many paths.
>>
>> Here is the mailing list thread:
>>
>> https://groups.google.com/forum/?utm_medium=email&utm_source=footer#!msg/scikit-fmm/XvV_KhCUN5c/nrJSQasCIwAJ
>>
>> And here is his code:
>>
>> https://github.com/panlei7/shortest_path
>>
>> We were/are considering adding fast path finding as a feature of
>> scikit-fmm. I stalled out on this discussion, I will ping Lei Pan and
>> cc you.
> Thanks! I looked into Lei Pan’s documents and wrote a little wrapper to apply your code from the wiki to my problem. It does work (at least it looks reasonable to me so far). However, after revisiting Dijkstra’s algorithm and the description of the FMM method, I’m left with the conceptual question of whether it wouldn’t be better to find the lines while computing the arrival time matrix. Although I like the idea of using the arrival time matrix and going back, it seems we are applying numerics to a problem already solved numerically. Lei Pan’s examples (to the extent I understand them at the moment) illustrate that this is not ideal.
>
> In the Dijkstra algorithm it shouldn’t be hard to find paths by some book keeping: For each node, when assigning an arrival time (tentative or final) also write down the node which is its parent. In addition to the arrival time matrix you then get a ‘parent matrix’. (I used this for the plot I mentioned in my last email). Is that not possible with FMM because we can’t identify a parent?
>
>> 3) I will have a look at the paper and let you know if I have any
>> ideas. It should be possible to define some set of approximately
>> fastest paths by introducing some randomness in the algorithm that
>> enumerates the paths. The algorithm looks at the gradient of the
>> travel time function and takes small steps in the down gradient
>> direction. You could perturb the down gradient direction slightly and
>> you would get a set of paths. You would have to be careful that the
>> path does not wander into the obstacles where the travel time gradient
>> is undefined.
> Thanks.
>
> Best regards from the Netherlands,
>
> Wolfram

[Wolfram Moebius]

Dear Jason,

Sorry for taking so long to get back to you and thank you for the reply.

I did some more thinking and learned more about the FMM algorithm. Though my suggestion might be fast, it seems not to be accurate, I think yours is better. I did some low-level “profiling” and came to the following conclusion: Generating the coefficients for the interpolation is not fast, but is only required once. The actual Runge Kutta part is surprisingly slow, though - in C it should be much faster (from all my simulation experience). Maybe leave RectBivariateSpline
in Python and do the actual path determination (with some more parameters) in C. If you are interested in including this to the module, I’d help (though I can’t do that from scratch alone). Otherwise, I’ll write my own little program for personal use.

Best, Wolfram

[Jason Furtney]

Dear Wolfram,

Thanks for the suggestion. Fortunately, this was relatively easy as
indexing is all done in one place. I added an optional keyword
argument 'periodic' that can be a bool or a sequence. The change is in
this branch:

https://github.com/scikit-fmm/scikit-fmm/tree/periodic

Could you try this out? I did minimal testing, but all the existing
tests still pass. Could you write some tests and/or update the
documentation to reflect this change?

Below is an example:

from skfmm import distance, travel_time
import numpy as np

print distance((1,-1,1,1,1,1,1,1))
print distance((1,-1,1,1,1,1,1,1), periodic=True)

#np.array([ 0.5, -0.5, 0.5, 1.5, 2.5, 3.5, 2.5, 1.5])

phi = np.ones((6,6))
phi[0,1]=-1

print distance(phi, periodic=False)
print
print distance(phi, periodic=True)
print
print distance(phi, periodic=(1,0))
print
print distance(phi, periodic=(0,1))
print
print travel_time(phi, speed=np.ones_like(phi), periodic=True)

Thanks,
Jason


On Sun, Jul 3, 2016 at 10:48 AM, Wolfram Moebius
<wolfram...@gmail.com> wrote:
> Hi Jason,
>
> I have another suggestion. I wasn’t sure how to communicate that so I’m writing an email instead of opening an issue. It’s a feature request: periodic boundary conditions in one or two directions. I can help (i.e., send you a suggested modification of the code), but I currently can’t review all the code carefully so I’d rely on your judgement.
>
> Best,
>
> Wolfram

[Jason Furtney]

Dear Wolfram,

Thanks for the reply on this. Yes, I am interested in getting this
capability into scikit-fmm so please keep me updated about your
progress on this.

It would be a good idea to establish a path finding problem which has
an exact solution so we can test the algorithm. Something like the
brachistochrone problem could be used as a test case. (Although the
brachistochrone problem defines the starting velocity as zero which
would be a problem.) Once we figure out an accurate method it should
be relatively easy to make the code fast.

I will be unavailable from today until Monday July 26th. I can help
out when I am back.

Thanks,
Jason

On Sun, Jul 3, 2016 at 10:48 AM, Wolfram Moebius
<wolfram...@gmail.com> wrote:

[Wolfram Moebius]

Dear Jason,

I took your code as a template and wrote my own version. The results look promising and the code could be a starting point for making it faster.

We can use the fact that what we’re doing is closely related to optics. I have one or two examples for which I have computed the path analytically for a nontrivial speed matrix.

Ping me once you’re back?

Best, Wolfram

[Jason Furtney]

Dear Wolfram,

I created a path finding example with an exact solution:
https://gist.github.com/jkfurtney/32514ded53d99a5869c8883277a698ca

If you have some other exact solutions let me know. These cases can
help find a good algorithm.

Here are a few thoughts:

How can we better calculate the ending point of the path? It would be
nice to have the final point exactly on the zero contour.

We are using 4th order Runge-Kutta integration which is expensive and
maybe overkill?

The step size in the path integration is set to the grid spacing. Is
this optimal? I tried decreasing this step size and the solution
diverged unexpectedly.

How to compare these numerical solutions to the exact path? Maybe
comparing both x- and y- components vs arc length would be better?

With a grid spacing of 1e-2 the error is 2.4435e-4 and with a grid
spacing of 1e-3 the error is 2.589e-6 so the convergence looks good.
Both the interpolation and integration are higher order than the
solution. Maybe most of the error comes from the travel_time solution?

For smoothly varying speed functions it should be possible to make the
travel_time code in scikit-fmm more accurate. Currently, the speed
function is taken as a constant in each cell. When a point is updated
the local speed is taken as a constant despite the fact that some of
the path to the interface could be in another cell. I will give this
more thought.

Before diving to deep into the above there is a more important source
of error in scikit-fmm. The initially frozen points (the points
adjacent the zero contour) have their distance established by linear
interpolation. These first-order errors are propagated by the
subsequent second-order marching algorithm. I have a proof-of-concept
to resolve this in this branch
https://github.com/scikit-fmm/scikit-fmm/tree/new-init
This branch uses bi-cubic interpolation to reconstruct the front
during the initialization. Some effort is needed to bring this branch
into master. There are some corner cases to resolve and I do not want
to add the scipy module as a dependency to scikit-fmm.

Thanks,
Jason



On Sun, Jul 10, 2016 at 7:43 AM, Wolfram Moebius

> --
> 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/FE925E53-7D0D-44FC-93E2-FA59B9A349C4%40gmail.com.
> For more options, visit https://groups.google.com/d/optout.

[Wolfram Moebius]

Dear Jason,

I just now managed to fork and clone and look at what would be helpful.

I managed to install your code for periodic boundary conditions and so far it looks good. I’ll write an example for periodic boundary conditions (which doubles as a test) and adjust the documentation. I’ll let you know and make a pull request.

Best, Wolfram

[Wolfram Moebius]

Dear Jason,

> I created a path finding example with an exact solution:
> https://gist.github.com/jkfurtney/32514ded53d99a5869c8883277a698ca
>
> If you have some other exact solutions let me know. These cases can
> help find a good algorithm.

Yes. I also have an exact solution in mind where one has to go around a masked region. I’ll write a piece of code after having worked on the periodic boundary condition aspect.

> Here are a few thoughts:
>
> How can we better calculate the ending point of the path? It would be
> nice to have the final point exactly on the zero contour.

Don’t know at the moment.

> We are using 4th order Runge-Kutta integration which is expensive and
> maybe overkill?

Possibly, see below.

> The step size in the path integration is set to the grid spacing. Is
> this optimal? I tried decreasing this step size and the solution
> diverged unexpectedly.

The divergence seems surprising. Independently, I thought we could use a solver with an adaptive step size, basically outsourcing the problem because essentially we’re solving an ODE, right? Python should provide us with this?!

> How to compare these numerical solutions to the exact path? Maybe
> comparing both x- and y- components vs arc length would be better?

Yes. I also like to plot differences between the numerics and the exact solution in addition / instead of the actual plot.

Before we go to more details below, here are two other suggestions (see the non-working example https://gist.github.com/wmoebius/272f9e5a16745ea97d1c19f500717b30):
- I use the fact that RectBivariateSpline can actually give you derivatives. I.e., no need to first take the gradient and then interpolate.
- We could use RectBivariateSpline to get the direction, but use the actual input (speed or velocity matrix) to get the magnitude of the gradient? You’re not doing this, do you?

> With a grid spacing of 1e-2 the error is 2.4435e-4 and with a grid
> spacing of 1e-3 the error is 2.589e-6 so the convergence looks good.
> Both the interpolation and integration are higher order than the
> solution. Maybe most of the error comes from the travel_time solution?

:)

> For smoothly varying speed functions it should be possible to make the
> travel_time code in scikit-fmm more accurate. Currently, the speed
> function is taken as a constant in each cell. When a point is updated
> the local speed is taken as a constant despite the fact that some of
> the path to the interface could be in another cell. I will give this
> more thought.

The FMM people should have thought about this deeply, no?

> Before diving to deep into the above there is a more important source
> of error in scikit-fmm. The initially frozen points (the points
> adjacent the zero contour) have their distance established by linear
> interpolation. These first-order errors are propagated by the
> subsequent second-order marching algorithm. I have a proof-of-concept
> to resolve this in this branch
> https://github.com/scikit-fmm/scikit-fmm/tree/new-init
> This branch uses bi-cubic interpolation to reconstruct the front
> during the initialization. Some effort is needed to bring this branch
> into master. There are some corner cases to resolve and I do not want
> to add the scipy module as a dependency to scikit-fmm.

I need to think about this.

Best, Wolfram

[Wolfram Moebius]

Hi Jason,

Let me reply to my own email as a follow-up.

>> I created a path finding example with an exact solution:
>> https://gist.github.com/jkfurtney/32514ded53d99a5869c8883277a698ca
>>
>> If you have some other exact solutions let me know. These cases can
>> help find a good algorithm.
> Yes. I also have an exact solution in mind where one has to go around a masked region. I’ll write a piece of code after having worked on the periodic boundary condition aspect.
>
>> Here are a few thoughts:
>>
>> How can we better calculate the ending point of the path? It would be
>> nice to have the final point exactly on the zero contour.
> Don’t know at the moment.

I think we could get phi at each point while computing the path and stop once phi changes sign. Then we could either find the final point by intersection of the path with the zero contour of phi (based on curve intersection routines) or we could just do bisection in the ODE solution until we’re sufficiently close to the contour.

>> We are using 4th order Runge-Kutta integration which is expensive and
>> maybe overkill?
> Possibly, see below.
>
>> The step size in the path integration is set to the grid spacing. Is
>> this optimal? I tried decreasing this step size and the solution
>> diverged unexpectedly.
> The divergence seems surprising. Independently, I thought we could use a solver with an adaptive step size, basically outsourcing the problem because essentially we’re solving an ODE, right? Python should provide us with this?!

Do you mind using scipy.integrate.ode?

>> How to compare these numerical solutions to the exact path? Maybe
>> comparing both x- and y- components vs arc length would be better?
> Yes. I also like to plot differences between the numerics and the exact solution in addition / instead of the actual plot.
>
> Before we go to more details below, here are two other suggestions (see the non-working example https://gist.github.com/wmoebius/272f9e5a16745ea97d1c19f500717b30):
> - I use the fact that RectBivariateSpline can actually give you derivatives. I.e., no need to first take the gradient and then interpolate.
> - We could use RectBivariateSpline to get the direction, but use the actual input (speed or velocity matrix) to get the magnitude of the gradient? You’re not doing this, do you?

I can draft something based on the above if you think it’s a good idea. If yes via gist or a pull request?

>> With a grid spacing of 1e-2 the error is 2.4435e-4 and with a grid
>> spacing of 1e-3 the error is 2.589e-6 so the convergence looks good.
>> Both the interpolation and integration are higher order than the
>> solution. Maybe most of the error comes from the travel_time solution?
> :)
>
>> For smoothly varying speed functions it should be possible to make the
>> travel_time code in scikit-fmm more accurate. Currently, the speed
>> function is taken as a constant in each cell. When a point is updated
>> the local speed is taken as a constant despite the fact that some of
>> the path to the interface could be in another cell. I will give this
>> more thought.
> The FMM people should have thought about this deeply, no?
>
>> Before diving to deep into the above there is a more important source
>> of error in scikit-fmm. The initially frozen points (the points
>> adjacent the zero contour) have their distance established by linear
>> interpolation. These first-order errors are propagated by the
>> subsequent second-order marching algorithm. I have a proof-of-concept
>> to resolve this in this branch
>> https://github.com/scikit-fmm/scikit-fmm/tree/new-init
>> This branch uses bi-cubic interpolation to reconstruct the front
>> during the initialization. Some effort is needed to bring this branch
>> into master. There are some corner cases to resolve and I do not want
>> to add the scipy module as a dependency to scikit-fmm.
> I need to think about this.

How can I help with this?

Best, Wolfram

[Jason Furtney]

Dear Wolfram,

Thanks for the note on this. Regarding using scipy.integrate.ode: My
philosophy for scikit-fmm has been provide an application-agnostic and
light-weight library with few dependencies. Adding a dependency on
SciPy adds around 500 MB of binaries to an install and requires a
FORTRAN compiler to build. scikit-fmm may be used on embedded devices
where this could be prohibitive.

For these reasons, I think we should spin the path finding capability
off into a new project which uses scikit-fmm as a dependency. This
gives us a lot of flexibility. It would be great if you could take the
lead on this. I can help with the development, testing, documentation
and the Python packaging. If this sounds good to you could you create
a new repo? You can give me push access or we can work with PRs. Can
you come up with a project name? Just check that the name is available
on PyPI. We can start with a few scripts in the repo and structure it
into a Python module later.

Yes, your scheme based on the ODE solvers in SciPy sounds good. I
don't use the magnitude of the gradient because the path can be found
with only the direction information. I am sure there are other ways to
find paths. We could also consider using the sundials ODE solvers.

Regarding the new-init branch: this is not as issue when the zero
contour of phi is axis-aligned like in our Brachistochrone curve
example (https://gist.github.com/jkfurtney/32514ded53d99a5869c8883277a698ca).
I am (hopefully) almost finished with the new-init branch.

There are currently no tests against analytic solutions for
skfmm.travel_time() where the speed function is smoothly varying.
There may be some ways to improve the accuracy of travel_time(); but,
without an exact solution we are flying blind. I think the
Brachistochrone curve example can be used for this. I found an
expression for the travel time along the Brachistochrone curve itself
but this is of little use for verification. It should be possible to
find an expression for travel time as a function of x and y. I will
work on this.

Thanks,
Jason

[Wolfram Moebius]

Dear Jason,

Thanks for the note on this. Regarding using scipy.integrate.ode: My
philosophy for scikit-fmm has been provide an application-agnostic and
light-weight library with few dependencies. Adding a dependency on
SciPy adds around 500 MB of binaries to an install and requires a
FORTRAN compiler to build. scikit-fmm may be used on embedded devices
where this could be prohibitive.

For these reasons, I think we should spin the path finding capability
off into a new project which uses scikit-fmm as a dependency. This
gives us a lot of flexibility. It would be great if you could take the
lead on this. I can help with the development, testing, documentation
and the Python packaging. If this sounds good to you could you create
a new repo? You can give me push access or we can work with PRs. Can
you come up with a project name? Just check that the name is available
on PyPI. We can start with a few scripts in the repo and structure it
into a Python module later.

Sounds good! If you think it’s useful, we could than add a lightweight-version of it (e.g., with Runge Kutta) to sciki-fmm later. I made a new repository and added you as a collaborator. Once it’s a little bit more polished, I’d give it a proper name etc… Let’s shift the scientific discussion to issues in the new project?

Yes, your scheme based on the ODE solvers in SciPy sounds good. I
don't use the magnitude of the gradient because the path can be found
with only the direction information. I am sure there are other ways to
find paths. We could also consider using the sundials ODE solvers.

OK

Regarding the new-init branch: this is not as issue when the zero
contour of phi is axis-aligned like in our Brachistochrone curve
example (https://gist.github.com/jkfurtney/32514ded53d99a5869c8883277a698ca).
I am (hopefully) almost finished with the new-init branch.

OK

There are currently no tests against analytic solutions for
skfmm.travel_time() where the speed function is smoothly varying.
There may be some ways to improve the accuracy of travel_time(); but,
without an exact solution we are flying blind. I think the
Brachistochrone curve example can be used for this. I found an
expression for the travel time along the Brachistochrone curve itself
but this is of little use for verification. It should be possible to
find an expression for travel time as a function of x and y. I will
work on this.

That’s overlapping with the new spin-off. I’ll keep thinking of a good example.

On a different note: Did you ever think about adding curvature to scikit-fmm (see http://suttond.github.io/scikit-geodesic/)?

Best, Wolfram

[Makki Karim]

Dear John,

I would like to thank all the "skfmm" contributors.

I've recently used the powerful skfmm distance function to realize this work for estimating curvatures for implicit surfaces: https://ieeexplore.ieee.org/abstract/document/9540688

Kind regards from France,

Karim