Segmentation fault when anisotropic weight < sqrt(1/2)

[brianapps]

Reproduction steps

When I run this code.

import edt
import numpy as np
d = np.array([[True, True], [True, False]])
print(edt.edt(d, [0.5, 0.5]))

The process exits with a segmentation fault

Segmentation fault (core dumped)

I'd expect the script to produce these distances

[[0.70710677 0.5       ]
 [0.5        0.        ]]

The same issue occurs with the equivalent C++ code.

bool buffer[] = {true, true, true, false};
float* out = edt::edt(buffer, 2, 2, 0.5, 0.5); 

Discussion

When debugging I noticed situations in squared_edt_1d_parabolic where the variable k becomes negative and memory outside various buffers is accessed.

With the above data set, the routine attempts to calculate the intercept of two parabolas using this equation (grabbed from the README).

Let w2 = w * w
Let factor1 = (r - q) * w2
Let factor2 = r + q

Let s = (f(r) - f(q) + factor1 * factor2) / (2 * factor1);

The parameters are w=0.5, r=1, q=0, f(r)=0.25 and f(q)=INFINITY. This gives s at -6.80e38 but this exceeds the limit of a float so s ends up set to -INFINITY.

When it reaches the while block

while (s <= ranges[k]) {
    k--;
    factor1 = (i - v[k]) * w2;
    factor2 =  i + v[k];
    s = (ff[i] - ff[v[k]] + factor1 * factor2) / (2.0 * factor1);
}

This condition is true for k=0 (because ranges[0] is initialised to -INFINITY). k then becomes -1 and the code ends up accessing values outside various buffers. A segmentation fault tends to follow at some point.

I believe that with anisotropic weights > sqrt(1/2), s will always be greater than -INFINITY and the under-run condition never occurs. My gut feeling is adding a guard on k (in two places) e.g.

while (k > 0 && s <= ranges[k]) {
    ...
}

Will fix the issue but I don't know the code well enough to be sure.

Workaround

Scale the anisotropic weights so they are at least one and divide the results by this scale. e.g.

import edt
import numpy as np
d = np.array([[True, True], [True, False]])
spacing = np.array([0.3, 0.5])
scale = 1.0 / spacing.min()
print(edt.edt(d, spacing * scale) / scale)

[william-silversmith]

Thank you so much for this detailed report! I'll be looking into fixing this shortly.

[william-silversmith]

I looked into this a bit more and it looks like what ends up happening is with the guard enabled, the ranges variable ends up looking like:

[ -inf, -inf, +inf ]

Since the render loop then looks for the parabola vertex just below +inf, this works out even if it's odd. This neatly resolves the problem I had handling infinities and finally makes this code more generalized. I was able to remove two additional scans over the array that converted to and from infinite and finite representations which leads to a small speedup.

Thanks for your help!

[brianapps]

Thanks a lot for doing this, I really appreciate your hard work putting this library together. I've checkout out your branch and run it over a much larger problem case -- it works and everything looks good to me. Pleased to hear this led to improvements elsewhere.

[william-silversmith]

For some reason, removing the tofinite/toinfinite passes didn't work so well on Windows. I don't have access to a copy to experiment so I added it back. I'll make a release in a bit.

[william-silversmith]

You can get the newest version 2.0.3. Thanks again for your detailed report!