[FIXED] Obtain the indices of neighboring elements within a radius

Issue

I have a NxN matrix and I want to find indices (or at least count) of non-zero elements within a radius R for each element in matrix. The matrix will be large so I plan to do the calculations using Numpy, and I need the indices because this matrix is linked to another with data-storage (which is a python list [[]] of objects).

The matrix looks like this:

where Green is the selected cell, Yellow - the surrounding area.

The main problem is the moving stencil, either in loop or vectorized.

I tried the following code which gives the correct coordinates (all within radius and excluding center) relative to the point, but I dont know how to apply it to the numpy array:

r = 2
neighbor_coords = []
for i in list(itertools.product(range(-r, r+1), repeat=2)):
    if any(i) and np.sqrt(i[0]**2 + i[1]**2) <= r:
        neighbor_coords.append(i)

>>> [(-2, 0), (-1, -1), (-1, 0), (-1, 1), (0, -2), (0, -1), (0, 1), (0, 2), (1, -1), (1, 0), (1, 1), (2, 0)]

Expected output is the indices or count of neighboring elements.

Solution

This might get you part way there. You can use the convolve function in scipy to convolve your matrix with a kernel that you generate. This below will count the non-zero entries within the reach of the kernel. Note that I transformed the input matrix to ones so that it was summing the non-zero locations, not the values.

In [74]: import numpy as np                                                     

In [75]: from scipy.ndimage.filters import *                                    

In [76]: a                                                                      
Out[76]: 
array([[1, 2, 0, 0],
       [1, 1, 0, 0],
       [0, 0, 0, 0],
       [5, 5, 5, 5]])

In [77]: a_ones = np.where(a>0, 1, 0)                                           

In [78]: a_ones                                                                 
Out[78]: 
array([[1, 1, 0, 0],
       [1, 1, 0, 0],
       [0, 0, 0, 0],
       [1, 1, 1, 1]])

In [79]: k                                                                      
Out[79]: 
array([[0, 1, 0],
       [1, 0, 1],
       [0, 1, 0]])

In [80]: convolve(a_ones, k, mode='constant')                                   
Out[80]: 
array([[2, 2, 1, 0],
       [2, 2, 1, 0],
       [2, 2, 1, 1],
       [1, 2, 2, 1]])

In [81]:        

Answered By - AirSquid