Comparing value with neighbor elements in numpy

Question

Let's say I have a numpy array

    a b c
A = i j k
    u v w

I want to compare the value central element with some of its eight neighbor elements (along the axis or along the diagonal). Is there any faster way except the nested for loop (it's too slow for big matrix)?

To be more specific, what I want to do is compare value of element with it's neighbors and assign new values.

For example:

    if (j == 1):
      if (j>i) & (j>k):
        j = 999
      else:
        j = 0
    if (j == 2):
      if (j>c) & (j>u):
        j = 999
      else:
        j = 0
   ...

something like this.

Answer 1

Your operation contains lots of conditionals, so the most efficient way to do it in the general case (any kind of conditionals, any kind of operations) is using loops. This could be done efficiently using numba or cython. In special cases, you can implement it using higher level functions in numpy/scipy. I'll show a solution for the specific example you gave, and hopefully you can generalize from there.

Start with some fake data:

A = np.asarray([
    [1, 1, 1, 2, 0],
    [1, 0, 2, 2, 2],
    [0, 2, 0, 1, 0],
    [1, 2, 2, 1, 0],
    [2, 1, 1, 1, 2]
])

We'll find locations in A where various conditions apply.

• 1a) The value is 1
• 1b) The value is greater than its horizontal neighbors
• 2a) The value is 2
• 2b) The value is greater than its diagonal neighbors

Find locations in A where the specified values occur:

cond1a = A == 1
cond2a = A == 2

This gives matrices of boolean values, of the same size as A. The value is true where the condition holds, otherwise false.

Find locations in A where each element has the specified relationships to its neighbors:

# condition 1b: value greater than horizontal neighbors
f1 = np.asarray([[1, 0, 1]])
cond1b = A > scipy.ndimage.maximum_filter(
    A, footprint=f1, mode='constant', cval=-np.inf)

# condition 2b: value greater than diagonal neighbors
f2 = np.asarray([
    [0, 0, 1],
    [0, 0, 0],
    [1, 0, 0]
])
cond2b = A > scipy.ndimage.maximum_filter(
    A, footprint=f2, mode='constant', cval=-np.inf)

As before, this gives matrices of boolean values indicating where the conditions are true. This code uses scipy.ndimage.maximum_filter(). This function iteratively shifts a 'footprint' to be centered over each element of A. The returned value for that position is the maximum of all elements for which the footprint is 1. The mode argument specifies how to treat implicit values outside boundaries of the matrix, where the footprint falls off the edge. Here, we treat them as negative infinity, which is the same as ignoring them (since we're using the max operation).

Set values of the result according to the conditions. The value is 999 if conditions 1a and 1b are both true, or if conditions 2a and 2b are both true. Else, the value is 0.

result = np.zeros(A.shape)
result[(cond1a & cond1b) | (cond2a & cond2b)] = 999

The result is:

[
    [  0,   0,   0,   0,   0],
    [999,   0,   0, 999, 999],
    [  0,   0,   0, 999,   0],
    [  0,   0, 999,   0,   0],
    [  0,   0,   0,   0, 999]
]

You can generalize this approach to other patterns of neighbors by changing the filter footprint. You can generalize to other operations (minimum, median, percentiles, etc.) using other kinds of filters (see scipy.ndimage). For operations that can be expressed as weighted sums, use 2d cross correlation.

This approach should be much faster than looping in python. But, it does perform unnecessary computations (for example, it's only necessary to compute the max when the value is 1 or 2, but we're doing it for all elements). Looping manually would let you avoid these computations. Looping in python would probably be much slower than the code here. But, implementing it in numba or cython would probably be faster because these tools generate compiled code.