Optimizing implementation of a rotating mask in numpy/scipy

I believe you will have a difficult time optimizing significantly using Python + scipy. However, I was able to make a small improvement by using as_strided to generate rot_filters directly, rather than through boolean indexing. This is based on a very simple n-dimensional windows function. (I wrote it to solve this problem before I realized that a 2d convolution function exists in scipy.) The following code provides a modest 10% speedup on my machine; see below for an explanation of how it works:

import numpy as np
from scipy import ndimage
from numpy.lib import stride_tricks

# pass in `as_strided` as a default arg to save a global lookup
def rotation_matrix2(section, _as_strided=stride_tricks.as_strided):
    section = section.reshape(5, 5)  # sqrt(section.size), sqrt(section.size)
    windows_shape = (3, 3, 3, 3)     # 5 - 3 + 1, 5 - 3 + 1, 3, 3
    windows_strides = section.strides + section.strides
    windows = _as_strided(section, windows_shape, windows_strides)
    rot_filters = windows.reshape(9, 9)
    return rot_filters[rot_filters.std(1).argmin(),:].mean()

def get_rotation_smooth(im, _rm=rotation_matrix2, **kwargs):
    return ndimage.filters.generic_filter(im, _rm, size=5, **kwargs)

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    from scipy.misc import lena
    im = lena()
    im2 = get_rotation_smooth(im)
    #plt.gray()      # Uncomment these lines for
    #plt.imshow(im2) # demo purposes.
    #plt.show()

The above function rotation_matrix2 is equivalent to the following two functions (which together are actually a bit slower than your original function, because windows is more generalized). This does exactly what your original code does -- creates 9 3x3 windows into a 5x5 array, and then reshapes them into a 9x9 array for processing.

def windows(a, w, _as_strided=stride_tricks.as_strided):
    windows_shape = tuple(sa - sw + 1 for sa, sw in zip(a.shape, w))
    windows_shape += w
    windows_strides = a.strides + a.strides
    return _as_strided(a, windows_shape, windows_strides)

def rotation_matrix1(section, _windows=windows):
    rot_filters = windows(section.reshape(5, 5), (3, 3)).reshape(9, 9)
    return rot_filters[rot_filters.std(1).argmin(),:].mean()

windows works with arrays of any dimension, as long as the window has the same number of dimensions. Here's a breakdown of how it works:

    windows_shape = tuple(sa - sw + 1 for sa, sw in zip(a.shape, w))

We can think of the windows array as a n-d array of n-d arrays. The shape of the outer n-d array is dictated by the window's degrees of freedom within the larger array; in every dimension, the number of positions the window can take is equal to the length of the larger array minus the length of the window plus one. In this case, we have a 3x3 window into a 5x5 array, so the outer 2-d array is a 3x3 array.

    windows_shape += w

The shape of the inner n-d array is the same as the shape of the window itself. In our case, this is again a 3x3 array.

Now for the strides. We have to define the strides for the outer n-d array, and for the inner n-d array. But it turns out that they're the same! After all, the window moves through the larger array in just the same way that an individual index moves through the array, right?

    windows_strides = a.strides + a.strides

Now we have all the information we need to create the windows:

    return _as_strided(a, windows_shape, windows_strides)