Numerical issue?

[jkfurtney]

Probably related to #18 and or #15
From: Weilun Sun

import skfmm
import numpy as np

phi = np.array([[ 1, 1, -1, -1, -1, -1, 1],
                         [ 1, 1, 1, 1, 1, 1, 1],
                         [ 1, 1, 1, 1, 1, 1, 3],
                         [ 1, 1, 1, 3, 3, 3, 3],
                         [ 3, 3, 3, 1, 3, 3, 3]])

speed = np.array([[ 1. , 1. , 1, 1, 1, 1, 1. ],
                             [ 1. , 1. , 0.2, 0.5, 0.997, 0.678, 1. ],
                             [ 1. , 1. , 1. , 1. , 1. , 1. , 1. ],
                             [ 1. , 1. , 1. , 1. , 1. , 1. , 1. ],
                             [ 1. , 1. , 1. , 1. , 1. , 1. , 1. ]])

skfmm.travel_time(phi,speed)

Try this one?
The input looks legit to me. zero boundary is clear and speed is all positive.
If I change that 0.2 to 0.5 (something larger) in speed, there's no error. So I'm guessing it could be some numerical issues.

[dsmccormick8492]

@jkfurtney have you been able to address these apparent numerical issues? I have been trying the code with some of my pseudo-speed fields and it will fail with travel_time(phi, speed, order=2), but will work with order=1.

Thanks

[jkfurtney]

Thanks for the message and the reminder. It has been in the back of my mind to get to this. Do you have any new examples of this bug? Thanks again and I will try to get to it soon.

[dsmccormick8492]

@jkfurtney Attached I have written some code using your FMM library. I put some wrappers around the distance() and travel_time() functions as I'm using it specifically for evaluating on 2D and 3D grids related to fluid flow topics where masking for non-communication cell values is important to treat and essential to the geodesic distance calculation. It's probably ugly code as I'm not a professional programmer.

I have two test cases in the main() block:

1. An synthetic grid using sin() functions to create some spatial variability in the grids. phi and speed are identical and there are no masked values. This grid works fine for the order=2 travel_time() calculation. In the binary case, the (geodesic) distance() is equivalent to the travel_time() results for contour propagation.

2. I import a 2D binary image (of a river delta). In the wrapper routines for distance() and travel_time() I mask the grids for zero values (non-river). I set as single zero-valued source point for the distance() and travel_time() calculations. In this case, in travel_time(), the order=2 fails and I have to fall back to order=1.

I have other examples with continuous grids that could be used to illustrate this inability to use the order=2 calculations.

I hope this is clear.

Please contact me if you have any further questions. And thanks again for your work on this. It's a great contribution.

Cheers, David

David McCormick

fmm_travel_time_test.zip

[kjexini]

We have also had lots of issues with the "Negative discriminant..." error. According to the PhD thesis by Hanna Källén p. 20, the travel time update should updated according to the quadratic equation (eq. 2.11 in the thesis) if it is solvable, but otherwise another equation should be used (eq. 2.14), which basically is that the new time should be 1/speed + min(neighbors' travel time). Could this be of help?

[jkfurtney]

Thanks for the message and for the link to the thesis, I will check this out. Could you provide an example of the "Negative discriminant error" please?

[kjexini]

The example from Weilun Sun above gives the same error I'm referring to: Negative discriminant in time marcher quadratic.

[jkfurtney]

OK, when I implement equation 2.14 I get what looks like a reasonable answer.

import skfmm
import numpy as np

phi = np.array([[ 1, 1, -1, -1, -1, -1, 1],
                [ 1, 1,  1,  1,  1,  1, 1],
                [ 1, 1,  1,  1,  1,  1, 3],
                [ 1, 1,  1,  3,  3,  3, 3],
                [ 3, 3,  3,  1,  3,  3, 3]])

speed = np.array([[ 1. , 1. , 1,   1,   1,     1,     1 ],
                  [ 1. , 1. , 0.2, 0.5, 0.997, 0.678, 1. ],
                  [ 1. , 1. , 1,   1,   1,     1,     1. ],
                  [ 1. , 1. , 1,   1,   1,     1,     1. ],
                  [ 1. , 1. , 1,   1,   1,     1,     1. ]])

skfmm.travel_time(phi,speed)

Out[1]:
array([[1.5, 0.5, 0.35355339, 0.5, 0.5,        0.35355339, 0.5  ],
       [2.5, 1.5, 2.5,        1,   0.50150451, 0.73746313, 1.5  ],
       [3.5, 2.5, 3,          2,   1.50150451, 1.73746313, 2.5  ],
       [4.5, 3.5, 4,          3,   2.50150451, 2.73746313, 3.5  ],
       [5.5, 4.5, 5,          4,   3.50150451, 3.73746313, 4.5  ]])

[jkfurtney]

Let me read the thesis a little more carefully and double check a few things. I want to double check that I am choosing the solution correctly and that that this generalizes to higher dimensions.

[jkfurtney]

For now here is a branch with eq 2.14: https://github.com/scikit-fmm/scikit-fmm/tree/travel-time-fix

[kjexini]

That was fast! =)

I haven't gone through the math for higher dimensions, but eq. 2.14 makes at least intuitive sense I think.

[kjexini]

I have tried out the fix, and it solves the bug, but for cases where the bug did not occur I get very different results compared to before. I have attached some debug data. When I run the following snippet with the master branch I get output 404 (segmented volume size), but when I run it with the fix code I get output 1521.

debug_file = debug_dir / '0047-2020-01-16 16:36:41.228766-cropped.pkl'
debug_data = pickle.loads(debug_file.read_bytes())
travel_time = skfmm.travel_time(**debug_data)
np.sum(travel_time < 3)

0047-2020-01-16 16:36:41.228766-cropped.pkl.zip

[jkfurtney]

OK thanks, it looks like there is more to the story than just added eq. 21.4 I will investigate this soon.

It would be really nice if we had some exact solutions for travel time for a set of cases where the discriminant transitions from positive to negative, can you think of anything?

[kjexini]

What I meant with the debug example is that there were some unintentional side effects from the fix. The example did not have any negative discriminant (since it worked before the fix). After consulting a C++ expert (my husband), he concluded that the problem is that the variable 'pass' is not assigned when the discriminant is positive (and then is either false or a random value). So he suggested adding *pass = true; on line 93 in travel_time_marcher.cpp, and also adding *pass = true; in distance_marcher.cpp e.g. on line 93. This solved my problem - the segmentation size goes back to what is was before the "Negative Discrimant"-fix.

Regarding exact solutions, I guess one could do it by hand, I don't have time unfortunately. Another idea is to get another implementation of fast march to compare to.

[jkfurtney]

OK got it, uninitialized variable error! I will check things over and get an updated release out soon. Thanks again.

[kjexini]

After thinking through this problem a bit more I see two issues:

1. The segmentations should not have been larger when using 2.14 compared to before (but this was the case in my debug example above). My intuition is that when you use 2.14 the "front" is approaching from one dimension only. This can only lead to larger travel times compared to when it is approaching from many dimensions. My best guess is that using 2.14 in the fix fails to take into account spacing different to 1 (dx option).

2. If you have more than two dimensions, I don't think anymore it is entirely sufficient to use 2.14 as a fallback (but I don't think it is too bad of an approximation). For example, if data is three-dimensional, failing to solve the quadratic equation means that you don't approach the target point from the full three dimensions. But then one should check if you can approach it from two dimensions before one retorts to 2.14, since 2.14 means approaching from only one dimension.

I would like to understand the equations in scikit-fmm better to investigate this a bit more, do you have a reference with their derivation?

[jkfurtney]

Thanks for the note and I agree that there are more issues. In the thesis, eq. 2.10 and beyond are only for first-order upwind derivatives and I have not checked to what extent the material that follows applies to second order. There is another issue (#18) that is not related to the negative discriminant, it may be related to mixing first- and second-order during the initialization.

The reference I worked with is "Level Set Methods and Fast Marching Methods" by Sethian. The approach is the same as in the thesis except second order upwind derivatives are used. Substituting these into eq 2.8 should give the polynomial coefficients the module uses in the point update. I have been meaning to write out this derivation fully and I will get to it when I can.

[f-fanni]

Hi, I added a pull request with a more generic fix for the second order update in case the discriminant is less than 0. All the examples that I found in this thread seems to pass, and the modification are quite minimal, but I don't have any exact solution to compare the current results.

[jkfurtney]

OK, better late than never. Thanks to @f-fanni we now have fix for this in master. I will ask around for testing and if everything looks good this coming week I will put a release out.