ndimage routines behave badly when output has memory overlap with inputs

[skenop]

Maybe the problem is following
...
tmp = grey_erosion(tmp, size, footprint, structure, output, mode, cval, origin)
tmp and input is the same now
so
numpy.subtract(tmp, input, out=tmp)
gives zero

Thank you

[ilayn]

@yuracpp Do you mind adding a reproducable code snippet that we can copy and paste to see the problem for ourselves?

It is not easy for us to understand anything from this.

[skenop]

import numpy as np
import scipy.ndimage as ndi

a = np.ones((3,3))
a[1,1] = 0
ndi.black_tophat(a, size=3)
array([[0., 0., 0.],
[0., 1., 0.],
[0., 0., 0.]])
ndi.black_tophat(a, size=3, output=a)
array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])

[akahard2dj]

In this part, when output variable gets in grey_erosion(), input data is always changed.
But, problem was the case of input = output.
I have designed the solution by comparing input and output.

the results are always same.

    if numpy.array_equal(input, output):
        tmp = grey_erosion(tmp, size, footprint, structure, None, mode,
                          cval, origin)
    else:
        tmp = grey_erosion(tmp, size, footprint, structure, output, mode,
                           cval, origin)

issue reproduce code

import numpy as np
import scipy.ndimage as ndi

a = np.ones((3,3))
a[1,1] = 0
a=ndi.black_tophat(a, size=3)
print('result1')
print(a)

a = np.ones((3,3))
a[1,1] = 0
a=ndi.black_tophat(a, size=3, output=a)
print('result2')
print(a)

a = np.ones((3,3))
a[1,1] = 0
b=np.ones((3,3))
a=ndi.black_tophat(a, size=3, output=b)
print('result3')
print(a)
print(b)

results

result1
[[0. 0. 0.]
 [0. 1. 0.]
 [0. 0. 0.]]

result2
[[0. 0. 0.]
 [0. 1. 0.]
 [0. 0. 0.]]

result3
[[0. 0. 0.]
 [0. 1. 0.]
 [0. 0. 0.]]
[[0. 0. 0.]
 [0. 1. 0.]
 [0. 0. 0.]]

[skenop]

thanks, the output is correct with this fix.
BTW white_tophat has symmetric issue

a = np.zeros((3,3))
a[1,1] = 1
print(ndi.white_tophat(a, size=3))
[[0. 0. 0.]
[0. 1. 0.]
[0. 0. 0.]]
print(ndi.white_tophat(a, size=3, output=a))
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]

[akahard2dj]

Yes, white_tophat() function has same problem.
I will fix them and PR.

[pv]

The problem is that these routines do not check if the output array overlaps in memory with the inputs.

As mentioned in #9578, same issue is in rank_filter and likely many other functions in ndimage.

Since the ndimage code is not written to enable inplace operations, the fix has to be to make copies along the lines of if np.may_share_memory(output, input): input = input.copy() --- or some equivalent --- in these functions. Possibly, it is simpler to do this in the C routines, in case it could be done in a single place in them.