Getting number of neighbours for a 3d numpy arrays

[Nobody]

Some common ways to handle the boundary condition:

1. Generate clamped indices, test for validity and substitute invalid
entries with an "identity" element. E.g.

ijk = np.mgrid[:a,:b,:c]
i,j,k = ijk
i0,j0,k0 = np.maximum(0,ijk-1)
i1,j1,k1 = np.minimum(np.array([[[[a,b,c]]]]).T-1,ijk+1)

n1 = (i>0  ) & myarray[i0,j,k]
n2 = (i<a-1) & myarray[i1,j,k]
n3 = (j>0  ) & myarray[i,j0,k]
n4 = (j<b-1) & myarray[i,j1,k]
n5 = (k>0  ) & myarray[i,j,k0]
n6 = (k<c-1) & myarray[i,j,k1]

2. Similar, but use np.roll() instead of clamped indices:

n1 = (i>0  ) & np.roll(myarray, 1,axis=0)
n2 = (i<a-1) & np.roll(myarray,-1,axis=0)
n3 = (j>0  ) & np.roll(myarray, 1,axis=1)
n4 = (j<b-1) & np.roll(myarray,-1,axis=1)
n5 = (k>0  ) & np.roll(myarray, 1,axis=2)
n6 = (k<c-1) & np.roll(myarray,-1,axis=2)

3. Split each axis into 3 ranges: start (has no predecessor), middle, end
(has no successor). Handle the permutations separately. So for 3
dimensions, you get 3x3x3=27 cases:

1 central block: (a-2)x(b-2)x(c-2)
3x2 = 6 faces: 2 each of 1x(b-2)x(c-2), (a-2)x1x(c-2), (a-2)x(b-2)x1
3x4 = 12 edges: 4 each of (a-2)x1x1, 1x(b-2)x1, 1x1x(c-2)
8 corners: 1x1x1

More complicated code-wise, but the sub-arrays are views rather than
copies. If the array is particularly large, the efficiency gain may
justify the complexity.