Accessing neighboring cells for numpy array

Question

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]

Answer 1

Please clarify your question

• Is it really intended that one loop iteration depends on the other, as mentioned by @jakevdp in the comments?
• If this is the case, how exactly should be border pixels be handeled? This will affect the whole result due to the dependence from one loop iteration to the others
• Please add a working reference implementation (You are getting an out of bounds error in your reference implementation)

Borders untouched, dependend loop iterations

I don't see any other way than using a compiler in this way. In this example I use Numba, but you can also do quite the same in Cython if this is preverred.

import numpy as np
import numba as nb

@nb.njit(fastmath=True)
def without_borders(arr):
  for x in range(1,arr.shape[0]-1):
    for y in range(1,arr.shape[1]-1):
      # Find 10% of current cell
      reduce_by = arr[x,y] * 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]
  return arr

Timings

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

#non-compiled verson: 6.7s
#compiled version:    6ms (the first call takes about 450ms due to compilation overhead)

This is realy easy to do an gives about a factor of 1000x. Depending on the first 3 Points there might be some more optimizations possible.

Answer 2

No need for loops, avoid the usual python loops, they are very slow. For greater efficiency, rely on numpy's build in matrix operation, "universal" functions, filters, masks and conditions whenever you can. https://realpython.com/numpy-array-programmin For complicated computations vectorization is not too bad see some chart and benchmarks Most efficient way to map function over numpy array (just do not use it for simpler matrix operations, like squaring of cells, build in functions will overperform)

Easy to see that each internal cell would be multiplied on .9 up to 8 times due 8 neighbors (that is reduced by .1), and additionally due to be a central cell, yet it cannot be reduced below .25/.9 = 5/18. For border and corner cell number of decreases fells to 6 and 3 times.

Therefore

x1 = 700  # for debugging use lesser arrays
x2 = 1400

neighbors = 8 # each internal cell has 8 neighbors


for i in range(neighbors):
     view1 = arr[1:-1, 1:-1] # internal cells only
     arr [1:x1, 1:-1] = np.multiply(view1,.9, where = view1 > .25)

arr [1:-1, 1:-1] *= .9 

Borders and corners are be treated in same way with neighbours = 5 and 3 respectively and different views. I guess all three cases can be joined in one formula with complicated where case, yet speed up would be moderate, as borders and corners take a small fraction of all cells.

Here I used a small loop, yet it just 8 repetitions. It should be can get rid of the loop too, using power, log, integer part and max functions, resulting in a bit clumsy, but somewhat faster one-liner, something around

      numpy.multiply( view1, x ** numpy.max( numpy.ceil( (numpy.log (* view1/x... / log(.9)

We can also try another useful technique, vectorization. The vectorization is building a function which then can be applied to all the elements of the array.

For a change, lets preset margins/thresholds to find out exact coefficient to multiply on . Here is what code to look like

n = 8
decrease_by = numpy.logspace(1,N,num=n, base=x, endpoint=False)

margins = decrease_by * .25

# to do : save border rows for further analysis, skip this for simplicity now    
view1 = a [1: -1, 1: -1]

def decrease(x):


    k = numpy.searchsorted(margin, a)
    return x * decrease_by[k]

f = numpy.vectorize(decrease)
f(view1)

Remark 1 One can try use different combinations of approaches, e.g. use precomputed margins with matrix arithmetics rather than vectorization. Perhaps there are even more tricks to slightly speed up each of above solutions or combinations of above.

Remark 2 PyTorch has many similarity with Numpy functionality but can greatly benefit from GPU. If you have a decent GPU consider PyTorch. There were attempt on gpu based numpy (gluon, abandoned gnumpy, minpy) More on gpu's https://stsievert.com/blog/2016/07/01/numpy-gpu/