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. -- https://mail.python.org/mailman/listinfo/python-list