I have an array of dimensions 100x100x1000. For every pixel of the array, I want to perform the following operation:
A = np.random.rand(100, 100, 1000)
values = [[0, 1, 2, 3],
[0, 2, 4, 5],
[3, 0, 1, 2],
[4, 3, 3, 0]]
vmax = np.amax(values)
for i in range(A.shape[0]-2):
for j in range(A.shape[1]-2):
for v in values:
for t in A.shape[2] - vmax:
result[i,j] = A[i+1, j, t+v[0]] * A[i, j+1, t+v[1]] * A[i+2, j+1, t+v[2]] * A[i+1, j+2, t+v[3]]
I wrote the code as inefficiently as I could just to make it clear that it's basically 4 nested loops for every single pixel. The most efficient way I've found is by indexing and multiplying:
A = np.random.rand(100, 100, 1000)
values = [[0, 1, 2, 3],
[0, 2, 4, 5],
[3, 0, 1, 2],
[4, 3, 3, 0]]
result = np.zeros((98, 98))
for v in values:
result += A[-1:1, :-2, v[0]:-vmax+v[0]] * A[:-2, 1:-1, v[1]:-vmax+v[1]] * A[-2:, 1:-1, v[2]:-vmax+v[2]] * A[-1:1, 2:, v[3]:-vmax+v[3]]
However, even this method can be extremely time-consuming for long combinations of "values" or for larger arrays "A". Are there any other approaches that could speed this up? My knowledge of multiprocessing is quite limited so I don't know if it would yield faster results than numpy's multiplications. Or since the operation is always the same, is there a more efficient way to vectorise the operation?