python - Filter array, store adjacency information -

let's have 2d array of (n, n) shape:

import numpy np my_array = np.random.random((n, n)) 

now want computations on "cells" of array, instance ones inside central part of array. avoid doing computations on cells i'm not interested in, here create boolean mask, in spirit:

my_mask = np.zeros_like(my_array, bool) my_mask[40:61,40:61] = true my_array[my_mask] = some_twisted_computations(my_array[my_mask]) 

but if some_twisted_computations() involves values of neighboring cells if inside mask? performance-wise, idea create "adjacency array" (len(my_mask), 4) shape, storing index of 4-connected neighbor cells in flat my_array[mask] array use in some_twisted_computations()? if yes, efficient options computing such adjacency array? should switch lower-level langage/other data structures?

my real-worlds arrays shapes around (1000,1000,1000), mask concerns small subset (~100000) of these values , of rather complex geometry. hope questions make sense...

edit: dirty , slow solution i've worked out:

wall = mask  = 0  top_neighbors = [] down_neighbors = [] left_neighbors = [] right_neighbors = [] indices = []  index, val in np.ndenumerate(wall):     if not val:         continue     indices += [index]     if wall[index[0] + 1, index[1]]:         down_neighbors += [(index[0] + 1, index[1])]     else:         down_neighbors += [i]     if wall[index[0] - 1, index[1]]:         top_neighbors += [(index[0] - 1, index[1])]     else:         top_neighbors += [i]     if wall[index[0], index[1] - 1]:         left_neighbors += [(index[0], index[1] - 1)]     else:         left_neighbors += [i]     if wall[index[0], index[1] + 1]:         right_neighbors += [(index[0], index[1] + 1)]     else:         right_neighbors += [i]     += 1   top_neighbors = [i if type(i) int else indices.index(i) in top_neighbors] down_neighbors = [i if type(i) int else indices.index(i) in down_neighbors] left_neighbors = [i if type(i) int else indices.index(i) in left_neighbors] right_neighbors = [i if type(i) int else indices.index(i) in right_neighbors] 

the best answer depend on nature of computations want do. example, if can expressed summations on neighboring pixels, np.convolve or scipy.signal.fftconvolve can nice solution.

for specific question of efficiently generating arrays of neighbor indices, might try this:

x = np.random.rand(100, 100) mask = x > 0.9  i, j = np.where(mask)  i_neighbors = i[:, np.newaxis] + [0, 0, -1, 1] j_neighbors = j[:, np.newaxis] + [-1, 1, 0, 0]  # need edge cases # best choice depend on application # here we'll change out-of-bounds neighbors # central point itself. i_neighbors = np.clip(i_neighbors, 0, 99) j_neighbors = np.clip(j_neighbors, 0, 99)  # compute vectorized result on neighbors # concrete example, here we'll standard deviation result = x[i_neighbors, j_neighbors].std(axis=1) 

the result array of values corresponding masked region, containing standard deviation of neighboring values. approach work whatever specific problem have in mind!

---

edit: given edited question above, here's how response can adapted generate arrays of indices in vectorized manner:

x = np.random.rand(100, 100) mask = x > -0.9  i, j = np.where(mask)  i_neighbors = i[:, np.newaxis] + [0, 0, -1, 1] j_neighbors = j[:, np.newaxis] + [-1, 1, 0, 0] i_neighbors = np.clip(i_neighbors, 0, 99) j_neighbors = np.clip(j_neighbors, 0, 99)  indices = np.zeros(x.shape, dtype=int) indices[mask] = np.arange(len(i))  neighbor_in_mask = mask[i_neighbors, j_neighbors]  neighbors = np.where(neighbor_in_mask,                      indices[i_neighbors, j_neighbors],                      np.arange(len(i))[:, none])  left_indices, right_indices, top_indices, bottom_indices = neighbors.t