[FIXED] Accessing neighboring cells for numpy array

Issue

How can I access and modify the surrounding 8 cells for a 2D numpy array in an efficient manner?

I have a 2D numpy array like this:

arr = np.random.rand(720, 1440)

For each grid cell, I want to reduce by 10% of the center cell, the surrounding 8 cells (fewer for corner cells), but only if the surrounding cell value exceeds 0.25. I suspect that the only way to do this is using a for loop but would like to see if there are better/faster solutions.

-- EDIT: For loop based soln:

arr = np.random.rand(720, 1440)

for (x, y), value in np.ndenumerate(arr):
    # Find 10% of current cell
    reduce_by = value * 0.1

    # Reduce the nearby 8 cells by 'reduce_by' but only if the cell value exceeds 0.25
    # [0] [1] [2]
    # [3] [*] [5]
    # [6] [7] [8]
    # * refers to current cell

    # cell [0]
    arr[x-1][y+1] = arr[x-1][y+1] * reduce_by if arr[x-1][y+1] > 0.25 else arr[x-1][y+1]

    # cell [1]
    arr[x][y+1] = arr[x][y+1] * reduce_by if arr[x][y+1] > 0.25 else arr[x][y+1]

    # cell [2]
    arr[x+1][y+1] = arr[x+1][y+1] * reduce_by if arr[x+1][y+1] > 0.25 else arr[x+1][y+1]

    # cell [3]
    arr[x-1][y] = arr[x-1][y] * reduce_by if arr[x-1][y] > 0.25 else arr[x-1][y]

    # cell [4] or current cell
    # do nothing

    # cell [5]
    arr[x+1][y] = arr[x+1][y] * reduce_by if arr[x+1][y] > 0.25 else arr[x+1][y]

    # cell [6]
    arr[x-1][y-1] = arr[x-1][y-1] * reduce_by if arr[x-1][y-1] > 0.25 else arr[x-1][y-1]

    # cell [7]
    arr[x][y-1] = arr[x][y-1] * reduce_by if arr[x][y-1] > 0.25 else arr[x][y-1]

    # cell [8]
    arr[x+1][y-1] = arr[x+1][y-1] * reduce_by if arr[x+1][y-1] > 0.25 else arr[x+1][y-1]

Solution

This answer assumes that you really want to do exactly what you wrote in your question. Well, almost exactly, since your code crashes because indices get out of bounds. The easiest way to fix that is to add conditions, like, e.g.,

if x > 0 and y < y_max:
    arr[x-1][y+1] = ...

The reason why the main operation cannot be vectorized using numpy or scipy is that all cells are “reduced” by some neighbor cells that have already been “reduced”. Numpy or scipy would use the unaffected values of the neighbors on each operation. In my other answer I show how to do this with numpy if you are allowed to group operations in 8 steps, each along the direction of one particular neighbor, but each using the unaffected value in that step for that neighbor. As I said, here I presume you have to proceed sequentially.

Before I continue, let me swap x and y in your code. Your array has a typical screen size, where 720 is the height and 1440 the width. Images are usually stored by rows, and the rightmost index in an ndarray is, by default, the one that varies more rapidly, so everything makes sense. It's admittedly counter-intuitive, but the correct indexing is arr[y, x].

The major optimization that can be applied to your code (that cuts execution time from ~9 s to ~3.9 s on my Mac) is not to assign a cell to itself when it's not necessary, coupled with in-place multiplication and with [y, x] instead of [y][x] indexing. Like this:

y_size, x_size = arr.shape
y_max, x_max = y_size - 1, x_size - 1
for (y, x), value in np.ndenumerate(arr):
    reduce_by = value * 0.1
    if y > 0 and x < x_max:
        if arr[y - 1, x + 1] > 0.25: arr[y - 1, x + 1] *= reduce_by
    if x < x_max:
        if arr[y    , x + 1] > 0.25: arr[y    , x + 1] *= reduce_by
    if y < y_max and x < x_max:
        if arr[y + 1, x + 1] > 0.25: arr[y + 1, x + 1] *= reduce_by
    if y > 0:
        if arr[y - 1, x    ] > 0.25: arr[y - 1, x    ] *= reduce_by
    if y < y_max:
        if arr[y + 1, x    ] > 0.25: arr[y + 1, x    ] *= reduce_by
    if y > 0 and x > 0:
        if arr[y - 1, x - 1] > 0.25: arr[y - 1, x - 1] *= reduce_by
    if x > 0:
        if arr[y    , x - 1] > 0.25: arr[y    , x - 1] *= reduce_by
    if y < y_max and x > 0:
        if arr[y + 1, x - 1] > 0.25: arr[y + 1, x - 1] *= reduce_by

The other optimization (that brings execution time further down to ~3.0 s on my Mac) is to avoid the boundary checks by using an array with extra boundary cells. We don't care what value the boundary contains, because it will never be used. Here is the code:

y_size, x_size = arr.shape
arr1 = np.empty((y_size + 2, x_size + 2))
arr1[1:-1, 1:-1] = arr
for y in range(1, y_size + 1):
    for x in range(1, x_size + 1):
        reduce_by = arr1[y, x] * 0.1
        if arr1[y - 1, x + 1] > 0.25: arr1[y - 1, x + 1] *= reduce_by
        if arr1[y    , x + 1] > 0.25: arr1[y    , x + 1] *= reduce_by
        if arr1[y + 1, x + 1] > 0.25: arr1[y + 1, x + 1] *= reduce_by
        if arr1[y - 1, x    ] > 0.25: arr1[y - 1, x    ] *= reduce_by
        if arr1[y + 1, x    ] > 0.25: arr1[y + 1, x    ] *= reduce_by
        if arr1[y - 1, x - 1] > 0.25: arr1[y - 1, x - 1] *= reduce_by
        if arr1[y    , x - 1] > 0.25: arr1[y    , x - 1] *= reduce_by
        if arr1[y + 1, x - 1] > 0.25: arr1[y + 1, x - 1] *= reduce_by
arr = arr1[1:-1, 1:-1]

For the records, if the operations could be vectorized using numpy or scipy, the speed-up with respect to this solution would be at least by a factor of 35 (measured on my Mac).

N.B.: if numpy did operations on array slices sequentially, the following would yield factorials (i.e., products of positive integers up to a number) – but it does not:

>>> import numpy as np
>>> arr = np.arange(1, 11)
>>> arr
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
>>> arr[1:] *= arr[:-1]
>>> arr
array([ 1,  2,  6, 12, 20, 30, 42, 56, 72, 90])

Answered By - Walter Tross