This is great, thank you for the PR.

I think one day we should implement a proper convolution instead of those matrix products (or even circular convolution with Fourier as an option).

you can add you name at the top of the file as Author . I will do the merge in the next days. I am still waiting for a few things but I will do a new release shortly.

Agreed a proper convolution does sound useful in some circumstances. Interesting to note a simple usage of scipy.ndimage.filter.gaussian_filter with this kernel is ~4X slower on the example images. I would guess matrix products are simply better optimized at this scale.

def K(x):
    # sigma to match previous behavior
    sigma = [np.sqrt(reg / 2) * dim for dim in x.shape]
    x = gaussian_filter(x, sigma=sigma, mode="constant")
    return x