Compute the product of neighborhood for each cell in a matrix with numpy/scipy

Question

I'm trying to implement an image processing algorithm which involves computing the product of 4-adjacent neighborhood for each cell. That is, to compute a new matrix Y for X where y[i, j] = x[i-1, j] * x[i, j-1] * x[i+1, j] * x[i, j+1]. Out-of-bound neighbors should be ignored.

Now I can only think of this approach: use scipy.ndimage.filters.correlate and pass in weights with zeros and one 1 to get four matrices, each containing the neighbor for each cell in a direction, like passing in weight = [[0, 0, 0], [1, 0, 0], [1, 1]] and I get a[i, j] = x[i-1, j], and with other weights I can get b[i, j] = x[i, j-1], c[i, j] = x[i+1, j], d[i, j] = x[i, j+1]. Then I use np.multiply to compute the product of these four matrices.

However, this approach is a bit too slow, and I can't ignore the boundaries. Is there another way to do it with numpy/scipy so I don't have to resort to for loops?

Answer 1

I think this better matches what you're asking for:

import numpy as np

x = np.array([
    [1, 2, 3, 4, 5],
    [6, 7, 8, 9, 1],
    [2, 3, 4, 5, 6],
    [7, 8, 9, 1, 2]
])

y = np.ones_like(x)
y[+1:, :] *= x[:-1, :]
y[:-1, :] *= x[+1:, :]
y[:, +1:] *= x[:, :-1]
y[:, :-1] *= x[:, +1:]

y
#>>> array([[  12,   21,   64,  135,    4],
#>>>        [  14,  288,  756,  160,  270],
#>>>        [ 126,  448, 1080,  216,   10],
#>>>        [  16,  189,   32,   90,    6]])

Note that the first *= can be an assignment if you need extra speed.