Pixel neighbors in 2d array (image) using Python

Question

I have a numpy array like this:

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

I need to create a function let's call it "neighbors" with the following input parameter:

• x: a numpy 2d array
• (i,j): the index of an element in a 2d array
• d: neighborhood radius

As output I want to get the neighbors of the cell i,j with a given distance d. So if I run

neighbors(im, i, j, d=1) with i = 1 and j = 1 (element value = 5) 

I should get the indices of the following values: [1,2,3,4,6,7,8,9]. I hope I make it clear. Is there any library like scipy which deal with this?

I've done something working but it's a rough solution.

def pixel_neighbours(self, p):

    rows, cols = self.im.shape

    i, j = p[0], p[1]

    rmin = i - 1 if i - 1 >= 0 else 0
    rmax = i + 1 if i + 1 < rows else i

    cmin = j - 1 if j - 1 >= 0 else 0
    cmax = j + 1 if j + 1 < cols else j

    neighbours = []

    for x in xrange(rmin, rmax + 1):
        for y in xrange(cmin, cmax + 1):
            neighbours.append([x, y])
    neighbours.remove([p[0], p[1]])

    return neighbours

How can I improve this?

Caravette

Answer 1

EDIT: ah crap, my answer is just writing im[i-d:i+d+1, j-d:j+d+1].flatten() but written in a incomprehensible way :)

---

The good old sliding window trick may help here:

import numpy as np
from numpy.lib.stride_tricks import as_strided

def sliding_window(arr, window_size):
    """ Construct a sliding window view of the array"""
    arr = np.asarray(arr)
    window_size = int(window_size)
    if arr.ndim != 2:
        raise ValueError("need 2-D input")
    if not (window_size > 0):
        raise ValueError("need a positive window size")
    shape = (arr.shape[0] - window_size + 1,
             arr.shape[1] - window_size + 1,
             window_size, window_size)
    if shape[0] <= 0:
        shape = (1, shape[1], arr.shape[0], shape[3])
    if shape[1] <= 0:
        shape = (shape[0], 1, shape[2], arr.shape[1])
    strides = (arr.shape[1]*arr.itemsize, arr.itemsize,
               arr.shape[1]*arr.itemsize, arr.itemsize)
    return as_strided(arr, shape=shape, strides=strides)

def cell_neighbors(arr, i, j, d):
    """Return d-th neighbors of cell (i, j)"""
    w = sliding_window(arr, 2*d+1)

    ix = np.clip(i - d, 0, w.shape[0]-1)
    jx = np.clip(j - d, 0, w.shape[1]-1)

    i0 = max(0, i - d - ix)
    j0 = max(0, j - d - jx)
    i1 = w.shape[2] - max(0, d - i + ix)
    j1 = w.shape[3] - max(0, d - j + jx)

    return w[ix, jx][i0:i1,j0:j1].ravel()

x = np.arange(8*8).reshape(8, 8)
print x

for d in [1, 2]:
    for p in [(0,0), (0,1), (6,6), (8,8)]:
        print "-- d=%d, %r" % (d, p)
        print cell_neighbors(x, p[0], p[1], d=d)

Didn't do any timings here, but it's possible this version has reasonable performance.

For more info, search the net with phrases "rolling window numpy" or "sliding window numpy".

Iraq

Answer 2

Have a look at scipy.ndimage.generic_filter.

As an example:

import numpy as np
import scipy.ndimage as ndimage

def test_func(values):
    print(values)
    return values.sum()


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

footprint = np.array([[1,1,1],
                      [1,0,1],
                      [1,1,1]])

results = ndimage.generic_filter(x, test_func, footprint=footprint)

By default, it will "reflect" the values at the boundaries. You can control this with the mode keyword argument.

However, if you're wanting to do something like this, there's a good chance that you can express your problem as a convolution of some sort. If so, it will be much faster to break it down into convolutional steps and use more optimized functions (e.g. most of scipy.ndimage).

Cementum

This method is slower than my custom own method using for loops – Snobbish

Answer 3

I don't know about any library functions for this, but you can easily write something like this yourself using the great slicing functionality of numpy:

import numpy as np
def neighbors(im, i, j, d=1):
    n = im[i-d:i+d+1, j-d:j+d+1].flatten()
    # remove the element (i,j)
    n = np.hstack((n[:len(n)//2], n[len(n)//2+1:] ))
    return n

Of course you should do some range checks to avoid out-of-bounds access.

Acevedo

What does b refer to here? – Caustic

b is not defined, as James said – Idell

Answer 4

By using max and min, you handle pixels at the upper and lower bounds:

im[max(i-1,0):min(i+2,i_end), max(j-1,0):min(j+2,j_end)].flatten()

Dimmer

Answer 5

I agree with Joe Kingtons response, just an add to the footprints

import numpy as np
from scipy.ndimage import generate_binary_structure
from scipy.ndimage import iterate_structure
foot = np.array(generate_binary_structure(2, 1),dtype=int)

or for bigger/different footprints for ex.

np.array(iterate_structure(foot , 2),dtype=int)

Oto

Answer 6

from skimage.util import view_as_windows
window_shape = (3, 3, 3)
neighborhoods = view_as_windows(img[200:220], window_shape) 

Here is one D if you need:

>>> A = np.arange(10)
>>> A
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> window_shape = (3,)
>>> B = view_as_windows(A, window_shape)
>>> B.shape
(8, 3)
>>> B
array([[0, 1, 2],
       [1, 2, 3],
       [2, 3, 4],
       [3, 4, 5],
       [4, 5, 6],
       [5, 6, 7],
       [6, 7, 8],
       [7, 8, 9]]

reference :https://scikit-image.org/docs/stable/api/skimage.util.html#skimage.util.view_as_windows

Moody

Answer 7

Possibly use a KDTree in SciPy ?

Slavic

Just for whatever it's worth, for regularly-gridded data, a quadtree is not ideal. Indexing the grid that the data is on gives you a much faster lookup for neighbors. – Cementum

Answer 8

We first init our matrix of interest using numpy

import numpy as np

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

print(x)

[[1 2 3]
 [4 5 6]
 [7 8 9]]

Our neighbors is a function of distance for instance we might be interested in neighbors of distance 2 this tells us how should we pad our matrix x. We choose to pad with zeros but you can pad with whatever you like might be mean,mode,median of a row/column

d = 2

x_padded = np.pad(x,d,mode='constant')

print(x_padded)

[[0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0]
 [0 0 1 2 3 0 0]
 [0 0 4 5 6 0 0]
 [0 0 7 8 9 0 0]
 [0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0]]

We use x_padded matrix to get neighbors of any value in matrix x. Let (i,j) and (s,t) be indexes of x and x_padded respectively. Now we need to translate (i,j) to (s,t) to get neighbors of (i,j)

i,j = 2,1
s,t = 2*d+i+1, 2*d+j+1

window = x_padded[i:s, j:t]

print(window)

[[0 1 2 3 0]
 [0 4 5 6 0]
 [0 7 8 9 0]
 [0 0 0 0 0]
 [0 0 0 0 0]]

Please Note!!! the indexes (i,j) point to any value you wish to get its neighbors in matrix x

One might wish to iterate over each point in matrix x, get its neighbors and do computation using the neighbors for instance in Image Processing, the convolution with a kernel. One might do the following to get neighbors of each pixel in an image x

for i in range(x.shape[0]):
    for j in range(x.shape[1]):
        i,j = 2,1
        s,t = 2*d+i+1, 2*d+j+1
        window = x_padded[i:s, j:t]

Benefit