How to Get the Neighboring Elements in a Numpy Array With Taking Boundaries into Account

how to get the neighboring elements in a numpy array with taking boundaries into account?

import numpy as np
a = np.array([0,1,2,3,4,5,6,7,8,9])
num_neighbor=3

for index in range(len(a)):
left = a[:index][-num_neighbor:]
right= a[index+1:num_neighbor+index+1]
print(index,left,right)

yields

(0, array([], dtype=int32), array([1, 2, 3]))
(1, array([0]), array([2, 3, 4]))
(2, array([0, 1]), array([3, 4, 5]))
(3, array([0, 1, 2]), array([4, 5, 6]))
(4, array([1, 2, 3]), array([5, 6, 7]))
(5, array([2, 3, 4]), array([6, 7, 8]))
(6, array([3, 4, 5]), array([7, 8, 9]))
(7, array([4, 5, 6]), array([8, 9]))
(8, array([5, 6, 7]), array([9]))
(9, array([6, 7, 8]), array([], dtype=int32))

The reason why a[index-num_neighbor:index] does not work when index<num_neighbor is because of slicing rules #3 and #4:

Given s[i:j]:

If i or j is negative, the index is relative to the end of the string:
len(s) + i or len(s) + j is substituted.

The slice of s from i to j is defined as the sequence of items with
index k such that i <= k < j. If i or j is greater than len(s), use
len(s). If i is omitted or None, use 0. If j is omitted or None, use
len(s). If i is greater than or equal to j, the slice is empty.

So when index=1, then a[index-num_neighbor:index] = a[-2:1] = a[10-2:1] = a[8:1] = [].

Python Numpy Array geht values of neighbours

You can take advantage of python indexing wrapping around for negative indices.

def wrap_nb(x,i,j):
return x[np.ix_(*((z-1, z, z+1-S) for z,S in zip((i,j), x.shape)))].ravel()

This requires i and j to be nonnegative and less than the shape of x.

If that is not guaranteed:

def wrap_nb(x,i,j):
return x[np.ix_(*(np.r_[z-1:z+2]%S for z,S in zip((i,j), x.shape)))].ravel()

Examples:

>>> wrap_nb(x,1,-2)
array([ 2, 3, 4, 6, 7, 8, 10, 11, 12])
>>> wrap_nb(x,0,-1)
array([15, 16, 13, 3, 4, 1, 7, 8, 5])
>>> wrap_nb(x,0,0)
array([16, 13, 14, 4, 1, 2, 8, 5, 6])

How to create an array of neighbors out from each element in an numpy 2d array

You can use np.roll() to shift the array and generate an array of the neighbors. This has the added benefit of "rolling" the boundaries (ie the right side of the rightmost element is the leftmost element).

In regards to your example, here is a possible implementation using np.roll although I doubt it is the cleanest:

left=np.roll(self.grid,1,axis=1)
right=np.roll(self.grid,-1,axis=1)
down=np.roll(self.grid,-1,axis=0)
up=np.roll(self.grid,1,axis=0)
leftbottom=np.roll(left,-1,axis=0)
rightbottom=np.roll(right,-1,axis=0)
lefttop=np.roll(left,1,axis=0)
righttop=np.roll(right,1,axis=0)

neighbors2 = [np.array(i) for i in zip(*(np.ndarray.flatten(i) \
for i in (left, right, down, up, leftbottom, rightbottom, lefttop, righttop)))]

Accessing neighboring cells for numpy array

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])

Find neighbors in a matrix?

It might be hard in other languages but in Python this is quite easy. Here is a function that can do what you asked for:

def neighbors(radius, row_number, column_number):
return [[a[i][j] if i >= 0 and i < len(a) and j >= 0 and j < len(a[0]) else 0
for j in range(column_number-1-radius, column_number+radius)]
for i in range(row_number-1-radius, row_number+radius)]

Here is a 2D list:

 a = [[ 11,  21,  31,  41,  51,  61,  71],
[ 12, 22, 32, 42, 52, 62, 72],
[ 13, 23, 33, 43, 53, 63, 73],
[ 14, 24, 34, 44, 54, 64, 74],
[ 15, 25, 35, 45, 55, 65, 75],
[ 16, 26, 36, 46, 56, 66, 76],
[ 17, 27, 37, 47, 57, 67, 77]]

See List comprehensions.

Updated missing "and" in the solution - pls review

How to find all neighbour values near the edge in array?

This is a connected component labeling problem. You could use scipy.ndimage to identify the connected components, check which slices of the found objects contain 0 as a starting point and use them to fill the new array:

from scipy import ndimage

# labels the connected components with a different digit
x_components, _ = ndimage.measurements.label(a, np.ones((3, 3)))
# returns slices with the bounding boxes
bboxes = ndimage.measurements.find_objects(x_components)
# fills a new array with 1 on those slices
b = np.zeros_like(a)
for bbox in s:
if bbox[0].start == 0:
b[bbox] = a[bbox]

print(b)

array([[0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
[0, 0, 0, 1, 1, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 1, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

Python get get average of neighbours in matrix with na value

Shot #1

This assumes you are looking to get sliding windowed average values in an input array with a window of 3 x 3 and considering only the north-west-east-south neighborhood elements.

For such a case, signal.convolve2d with an appropriate kernel could be used. At the end, you need to divide those summations by the number of ones in kernel, i.e. kernel.sum() as only those contributed to the summations. Here's the implementation -

import numpy as np
from scipy import signal

# Inputs
a = [[1,2,3],[3,4,5],[5,6,7],[4,8,9]]

# Convert to numpy array
arr = np.asarray(a,float)

# Define kernel for convolution
kernel = np.array([[0,1,0],
[1,0,1],
[0,1,0]])

# Perform 2D convolution with input data and kernel
out = signal.convolve2d(arr, kernel, boundary='wrap', mode='same')/kernel.sum()

Shot #2

This makes the same assumptions as in shot #1, except that we are looking to find average values in a neighborhood of only zero elements with the intention to replace them with those average values.

Approach #1: Here's one way to do it using a manual selective convolution approach -

import numpy as np

# Convert to numpy array
arr = np.asarray(a,float)

# Pad around the input array to take care of boundary conditions
arr_pad = np.lib.pad(arr, (1,1), 'wrap')

R,C = np.where(arr==0) # Row, column indices for zero elements in input array
N = arr_pad.shape[1] # Number of rows in input array

offset = np.array([-N, -1, 1, N])
idx = np.ravel_multi_index((R+1,C+1),arr_pad.shape)[:,None] + offset

arr_out = arr.copy()
arr_out[R,C] = arr_pad.ravel()[idx].sum(1)/4

Sample input, output -

In [587]: arr
Out[587]:
array([[ 4., 0., 3., 3., 3., 1., 3.],
[ 2., 4., 0., 0., 4., 2., 1.],
[ 0., 1., 1., 0., 1., 4., 3.],
[ 0., 3., 0., 2., 3., 0., 1.]])

In [588]: arr_out
Out[588]:
array([[ 4. , 3.5 , 3. , 3. , 3. , 1. , 3. ],
[ 2. , 4. , 2. , 1.75, 4. , 2. , 1. ],
[ 1.5 , 1. , 1. , 1. , 1. , 4. , 3. ],
[ 2. , 3. , 2.25, 2. , 3. , 2.25, 1. ]])

To take care of the boundary conditions, there are other options for padding. Look at numpy.pad for more info.

Approach #2: This would be a modified version of convolution based approach listed earlier in Shot #1. This is same as that earlier approach, except that at the end, we selectively replace
the zero elements with the convolution output. Here's the code -

import numpy as np
from scipy import signal

# Inputs
a = [[1,2,3],[3,4,5],[5,6,7],[4,8,9]]

# Convert to numpy array
arr = np.asarray(a,float)

# Define kernel for convolution
kernel = np.array([[0,1,0],
[1,0,1],
[0,1,0]])

# Perform 2D convolution with input data and kernel
conv_out = signal.convolve2d(arr, kernel, boundary='wrap', mode='same')/kernel.sum()

# Initialize output array as a copy of input array
arr_out = arr.copy()

# Setup a mask of zero elements in input array and
# replace those in output array with the convolution output
mask = arr==0
arr_out[mask] = conv_out[mask]

Remarks: Approach #1 would be the preferred way when you have fewer number of zero elements in input array, otherwise go with Approach #2.

Numpy: fast calculations considering items' neighbors and their position inside the array

EDIT I have kept my original answer at the bottom. As Paul points out in the comments, the original answer didn't really answer the OP's question, and could be more easily achieved with an ndimage filter. The following much more cumbersome function should do the right thing. It takes two arrays, a and c, and returns the windowed minimum of a and the values in c at the positions of the windowed minimums in a:

def neighbor_min(a, c):
ac = np.concatenate((a[None], c[None]))
rows, cols = ac.shape[1:]
ret = np.empty_like(ac)

# Fill in the center
win_ac = as_strided(ac, shape=(2, rows-2, cols, 3),
strides=ac.strides+ac.strides[1:2])
win_ac = win_ac[np.ogrid[:2, :rows-2, :cols] +
[np.argmin(win_ac[0], axis=2)]]
win_ac = as_strided(win_ac, shape=(2, rows-2, cols-2, 3),
strides=win_ac.strides+win_ac.strides[2:3])
ret[:, 1:-1, 1:-1] = win_ac[np.ogrid[:2, :rows-2, :cols-2] +
[np.argmin(win_ac[0], axis=2)]]

# Fill the top, bottom, left and right borders
win_ac = as_strided(ac[:, :2, :], shape=(2, 2, cols-2, 3),
strides=ac.strides+ac.strides[2:3])
win_ac = win_ac[np.ogrid[:2, :2, :cols-2] +
[np.argmin(win_ac[0], axis=2)]]
ret[:, 0, 1:-1] = win_ac[:, np.argmin(win_ac[0], axis=0),
np.ogrid[:cols-2]]
win_ac = as_strided(ac[:, -2:, :], shape=(2, 2, cols-2, 3),
strides=ac.strides+ac.strides[2:3])
win_ac = win_ac[np.ogrid[:2, :2, :cols-2] +
[np.argmin(win_ac[0], axis=2)]]
ret[:, -1, 1:-1] = win_ac[:, np.argmin(win_ac[0], axis=0),
np.ogrid[:cols-2]]
win_ac = as_strided(ac[:, :, :2], shape=(2, rows-2, 2, 3),
strides=ac.strides+ac.strides[1:2])
win_ac = win_ac[np.ogrid[:2, :rows-2, :2] +
[np.argmin(win_ac[0], axis=2)]]
ret[:, 1:-1, 0] = win_ac[:, np.ogrid[:rows-2],
np.argmin(win_ac[0], axis=1)]
win_ac = as_strided(ac[:, :, -2:], shape=(2, rows-2, 2, 3),
strides=ac.strides+ac.strides[1:2])
win_ac = win_ac[np.ogrid[:2, :rows-2, :2] +
[np.argmin(win_ac[0], axis=2)]]
ret[:, 1:-1, -1] = win_ac[:, np.ogrid[:rows-2],
np.argmin(win_ac[0], axis=1)]
# Fill the corners
win_ac = ac[:, :2, :2]
win_ac = win_ac[:, np.ogrid[:2],
np.argmin(win_ac[0], axis=-1)]
ret[:, 0, 0] = win_ac[:, np.argmin(win_ac[0], axis=-1)]
win_ac = ac[:, :2, -2:]
win_ac = win_ac[:, np.ogrid[:2],
np.argmin(win_ac[0], axis=-1)]
ret[:, 0, -1] = win_ac[:, np.argmin(win_ac[0], axis=-1)]
win_ac = ac[:, -2:, -2:]
win_ac = win_ac[:, np.ogrid[:2],
np.argmin(win_ac[0], axis=-1)]
ret[:, -1, -1] = win_ac[:, np.argmin(win_ac[0], axis=-1)]
win_ac = ac[:, -2:, :2]
win_ac = win_ac[:, np.ogrid[:2],
np.argmin(win_ac[0], axis=-1)]
ret[:, -1, 0] = win_ac[:, np.argmin(win_ac[0], axis=-1)]

return ret

The return is a (2, rows, cols) array that can be unpacked into the two arrays:

>>> a = np.random.randint(100, size=(5,5))
>>> c = np.random.randint(100, size=(5,5))
>>> a
array([[42, 54, 18, 88, 26],
[80, 65, 83, 31, 4],
[51, 52, 18, 88, 52],
[ 1, 70, 5, 0, 89],
[47, 34, 27, 67, 68]])
>>> c
array([[94, 94, 29, 6, 76],
[81, 47, 67, 21, 26],
[44, 92, 20, 32, 90],
[81, 25, 32, 68, 25],
[49, 43, 71, 79, 77]])
>>> neighbor_min(a, c)
array([[[42, 18, 18, 4, 4],
[42, 18, 18, 4, 4],
[ 1, 1, 0, 0, 0],
[ 1, 1, 0, 0, 0],
[ 1, 1, 0, 0, 0]],

[[94, 29, 29, 26, 26],
[94, 29, 29, 26, 26],
[81, 81, 68, 68, 68],
[81, 81, 68, 68, 68],
[81, 81, 68, 68, 68]]])

The OP's case could then be solved as:

def bd_from_ac(a, c):
b,d = neighbor_min(a, c)
return a*b, d

And while there is a serious performance hit, it is pretty fast still:

In [3]: a = np.random.rand(1000, 1000)

In [4]: c = np.random.rand(1000, 1000)

In [5]: %timeit bd_from_ac(a, c)
1 loops, best of 3: 570 ms per loop

You are not really using the coordinates of the minimum neighboring element for anything else than fetching it, so you may as well skip that part and create a min_neighbor function. If you don't want to resort to cython for fast looping, you are going to have to go with rolling window views, such as outlined in Paul's link. This will typically convert your (m, n) array into a (m-2, n-2, 3, 3) view of the same data, and you would then apply np.min over the last two axes.

Unfortunately you have to apply it one axis at a time, so you will have to create a (m-2, n-2, 3) copy of your data. Fortunately, you can compute the minimum in two steps, first windowing and minimizing along one axis, then along the other, and obtain the same result. So at most you are going to have intermediate storage the size of your input. If needed, you could even reuse the output array as intermediate storage and avoid memory allocations, but that is left as exercise...

The following function does that. It is kind of lengthy because it has to deal not only with the central area, but also with the special cases of the four edges and four corners. Other than that it is a pretty compact implementation:

def neighbor_min(a):
rows, cols = a.shape
ret = np.empty_like(a)

# Fill in the center
win_a = as_strided(a, shape=(m-2, n, 3),
strides=a.strides+a.strides[:1])
win_a = win_a.min(axis=2)
win_a = as_strided(win_a, shape=(m-2, n-2, 3),
strides=win_a.strides+win_a.strides[1:])
ret[1:-1, 1:-1] = win_a.min(axis=2)

# Fill the top, bottom, left and right borders
win_a = as_strided(a[:2, :], shape=(2, cols-2, 3),
strides=a.strides+a.strides[1:])
ret[0, 1:-1] = win_a.min(axis=2).min(axis=0)
win_a = as_strided(a[-2:, :], shape=(2, cols-2, 3),
strides=a.strides+a.strides[1:])
ret[-1, 1:-1] = win_a.min(axis=2).min(axis=0)
win_a = as_strided(a[:, :2], shape=(rows-2, 2, 3),
strides=a.strides+a.strides[:1])
ret[1:-1, 0] = win_a.min(axis=2).min(axis=1)
win_a = as_strided(a[:, -2:], shape=(rows-2, 2, 3),
strides=a.strides+a.strides[:1])
ret[1:-1, -1] = win_a.min(axis=2).min(axis=1)

# Fill the corners
ret[0, 0] = a[:2, :2].min()
ret[0, -1] = a[:2, -2:].min()
ret[-1, -1] = a[-2:, -2:].min()
ret[-1, 0] = a[-2:, :2].min()

return ret

You can now do things like:

>>> a = np.random.randint(10, size=(5, 5))
>>> a
array([[0, 3, 1, 8, 9],
[7, 2, 7, 5, 7],
[4, 2, 6, 1, 9],
[2, 8, 1, 2, 3],
[7, 7, 6, 8, 0]])
>>> neighbor_min(a)
array([[0, 0, 1, 1, 5],
[0, 0, 1, 1, 1],
[2, 1, 1, 1, 1],
[2, 1, 1, 0, 0],
[2, 1, 1, 0, 0]])

And your original question can be solved as:

def bd_from_ac(a, c):
return a*neighbor_min(a), neighbor_min(c)

As a performance benchmark:

In [2]: m, n = 1000, 1000

In [3]: a = np.random.rand(m, n)

In [4]: c = np.random.rand(m, n)

In [5]: %timeit bd_from_ac(a, c)
1 loops, best of 3: 123 ms per loop


Related Topics



Leave a reply



Submit