I am trying to use scipy.ndimage.correlate to replicate the output of IDL convol() function. The IDL function only calculates elements where there is full overlap between the input and the kernel.
So, for example if
input = np.array([[1,2,3],[4,5,6],[7,8,9]]]
kernel = np.array([[0,0,0],[0,1,0],[0,0,0]])
scipy.ndimage.correlate(input, kernel) = [[0,0,0],[0,5,0],[0,0,0]]
Is this possible?
In SciPy, if you want to control the behavior of how scipy.ndimage.correlate() extends the input
argument beyond the edges of the input, you can use the mode
argument. This unfortunately doesn't help you, because there is no mode argument that results in setting the output to zero when an edge is encountered.
That being said, I have two ideas for how you could solve this.
Step one, determine the places where your correlation is valid. You could do this by convolving a uniform filter of the same size as your kernel over an input of all zeros. This should be zero everywhere. Then, you tell it to fill in the edges with a non-zero edge, and you find all the places where the kernel overlaps with some edge value. You check if the output is zero. If it is, then the kernel is completely contained within the input.
Step two, you can use this array to set the areas where the kernel overlaps with an edge to zero.
Here's a full example.
import numpy as np
import scipy
input = np.array([[1,2,3],[4,5,6],[7,8,9]])
kernel = np.array([[0,0,0],[0,1,0],[0,0,0]])
def idl_convol(input_array, kernel):
validity_input = np.zeros_like(input_array, dtype='float32')
defined = scipy.ndimage.uniform_filter(validity_input, kernel.shape, mode='constant', cval=1) == 0
return defined * scipy.ndimage.correlate(input_array, kernel)
print(idl_convol(input, kernel))
Output:
[[0 0 0]
[0 5 0]
[0 0 0]]
I also tried an approach which computes this by setting the edges of the array to zero. On large arrays, this is much faster than idea #1.
import numpy as np
import scipy
input = np.array([[1,2,3],[4,5,6],[7,8,9]])
kernel = np.array([[0,0,0],[0,1,0],[0,0,0]])
def idl_convol_trim(input_array, kernel):
kernel = np.asarray(kernel)
output = scipy.ndimage.correlate(input_array, kernel)
trim_edges(output, kernel)
return output
def trim_edges(output, kernel):
for axis_pos, axis_length in enumerate(kernel.shape):
assert axis_length % 2 == 1, "Kernel must have odd size"
zeros_to_add = axis_length // 2
if zeros_to_add == 0:
continue
output_axis = np.moveaxis(output, axis_pos, 0)
output_axis[:zeros_to_add] = 0
output_axis[-zeros_to_add:] = 0
print(idl_convol_trim(input, kernel))
Timing data for 1000x1000 input:
Find valid via convolution
14.3 ms ± 177 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Mask invalid via indexing
6.38 ms ± 25.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
The downside of this approach is 1) it's more complex, and 2) it's harder to modify to make changes like an even-sized kernel or non-centered kernel.
Note for people who are not using IDL: these code examples use correlate()
, not convolve()
. This is because IDL defines convol()
as correlation. If you actually want convolution, use scipy.ndimage.convolve().