Given a 2d matrix, for instance:
A = array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
how can I select NxN window around a given element so that the window is filled with an arbitrary (e.g. mean) value if it goes outside the original array?
Example:
neighbors(A, x=0, y=0, N=3)
would yield
array([[ 2.5, 2.5, 2.5],
[ 2.5, 0, 1],
[ 2.5, 4, 5]])
You can pad your array. Padding will extend your array with the desired boundary conditions (see the mode parameter for all possible options):
>>> A = np.array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
>>> N = 5 # 5x5 windows
>>> B = np.pad(A, N//2, mode='reflect')
>>> B
array([[10, 9, 8, 9, 10, 11, 10, 9],
[ 6, 5, 4, 5, 6, 7, 6, 5],
[ 2, 1, 0, 1, 2, 3, 2, 1],
[ 6, 5, 4, 5, 6, 7, 6, 5],
[10, 9, 8, 9, 10, 11, 10, 9],
[14, 13, 12, 13, 14, 15, 14, 13],
[10, 9, 8, 9, 10, 11, 10, 9],
[ 6, 5, 4, 5, 6, 7, 6, 5]])
As you can see, the original array is in the center of the matrix, padded by 2 rows and 2 columns (N//2 = 5//2 = 2, both from lef/right and bottom/top). The padded elements are reflected.
For this new array, you can access to the desired neighbours window by just normally indexing the array:
>>> x = 1; y = 1 # corresponds to number 5 in the original array
>>> B[y:y+N, x:x+N]
array([[ 5, 4, 5, 6, 7],
[ 1, 0, 1, 2, 3],
[ 5, 4, 5, 6, 7],
[ 9, 8, 9, 10, 11],
[13, 12, 13, 14, 15]])
You can choose other padding methods, mean is one of the options.
I think this is not too far off:
def neighbors(arr, x, y, N):
left = max(0, x - N)
right = min(arr.shape[0], x + N)
top = max(0, y - N)
bottom = min(arr.shape[1], y + N)
window = arr[left:right+1,top:bottom+1]
fillval = window.mean()
result = np.empty((2*N+1, 2*N+1))
result[:] = fillval
ll = N - x
tt = N - y
result[ll+left:ll+right+1,tt+top:tt+bottom+1] = window
return result
Here's one way with np.lib.pad and np.ix_ -
Ap = np.lib.pad(A.astype(float),1, 'constant', constant_values=(np.nan,np.nan))
Acut = Ap[np.ix_(np.arange(N)+x,np.arange(N)+y)]
Acut[np.isnan(Acut)] = np.nanmean(Acut)
Sample run -
In [76]: A
Out[76]:
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
In [77]: x=1; y=0; N = 3
In [78]: Ap = np.lib.pad(A.astype(float),1, 'constant', constant_values=(np.nan,np.nan))
...: Acut = Ap[np.ix_(np.arange(N)+x,np.arange(N)+y)]
...: Acut[np.isnan(Acut)] = np.nanmean(Acut)
...:
In [79]: Acut
Out[79]:
array([[ 4.5, 0. , 1. ],
[ 4.5, 4. , 5. ],
[ 4.5, 8. , 9. ]])
scipy.ndimage.generic_filter.