Separable filter on numpy array

Question

Say I have a numpy array a, and I want to create a new array, b such that b[i, j] is a function of, say:

a[i-1, j-1], a[i-1, j  ], a[i-1, j+1],
a[i  , j-1], a[i  , j  ], a[i  , j+1],
a[i+1, j-1], a[i+1, j  ], a[i+1, j+1]

What would be the fastest way to do this?

As this is a separable filter, is there any way to run this in multiple threads? (not processes, because I would have to copy the data back)

Or is writing C code to bypass the GIL mandatory?

Partial solutions (like assuming the function is linear) are welcome too.

Answer 1

An idealized numpy way of working with a sliding window like this is to construct a 4D array

C.shape = (N,M,3,3)

where

C[i,j,:,:] = np.array([a[i-1, j-1], a[i-1, j  ], a[i-1, j+1],
                       a[i  , j-1], a[i  , j  ], a[i  , j+1],
                       a[i+1, j-1], a[i+1, j  ], a[i+1, j+1]])

and write your function do some sort of reduction on the last 2 dimensions. sum or mean would be typical, e.g.

B = C.sum(axis=(2,3))

Other SO questions show how to use np.lib.stride_tricks.as_strided to construct such an array. But with only a 3x3 subarray, it might be just as fast to do something like

C = np.zeros((N,M,3,3))
C[:,:,0,0] = a[:-1,:-1]
etc.

(or use hstack and vstack to the same effect).

But a nice thing (or maybe not so nice) about the strided approach is that it doesn't involve copy any data of a - it is just a view.

As to splitting the job into pieces, I can imagine using slices of C (on the 1st 2 dimensions), e.g.

 C[0:100,0:100,:,:].sum(axis=(2,3))