Python中的滑动Gabor滤波器

Question

从skimage的gabor滤波器示例中获取，计算图像的gabor滤波器很容易：

import numpy as np

from scipy import ndimage as nd

from skimage import data
from skimage.util import img_as_float
from skimage.filter import gabor_kernel

brick = img_as_float(data.load('brick.png'))

kernel = np.real(gabor_kernel(0.15, theta = 0.5 * np.pi,sigma_x=5, sigma_y=5))

filtered = nd.convolve(brick, kernel, mode='reflect')

mean = filtered.mean()
variance = filtered.var()

brick就是一个numpy数组。假设我有一个5000*5000的numpy数组。我想要实现的是生成两个新的5000*5000的numpy数组，其中像素是15*15窗口中心的gabor滤波器的平均值和方差。

有人能帮我实现这个吗？

编辑

为什么我会被downvote？无论如何，为了澄清，我展示了如何在单个图像上计算gabor滤波器的示例。我只想在一个非常大的图像的小正方形子集上计算gabor滤波器（因此使用滑动窗口）。

- JEquihua

---

你能否解释一下你所说的“gabor滤波器的均值和方差值”的含义？该滤波器似乎是一个31x31像素的线性滤波器。你想要在这个区域（31x31）内得到滚动平均值和平均数，还是想要得到滤波后图像的滚动平均值和方差？ - DrV

Answer 1

我不知道有标准的方法来做到这一点，但您可以直接自己完成。

在卷积中，每个像素是移位Gabor滤波器值乘以图像像素的总和。也就是说，每个卷积像素基本上是均值，只是一个常数归一化因子，所以filtered基本上是您的均值。

方差稍微有些困难，因为它是平方和，当然，您需要在计算总和之前先计算平方。但是，您可以通过预先平方图像和核来轻松完成这项工作。

N = kernel.shape[0]*kernel.shape[1]
mean = nd.convolve(brick, kernel, mode='reflect')/N
var = nd.convolve(brick*brick, kernel*kernel, mode='reflect')/N - mean*mean

- tom10

Answer 2

如果您只想计算图像的滑动平均值（使用所有1的正方形核进行卷积），则快速方法如下：

# fsize is the filter size in pixels
# integrate in the X direction
r_sum = numpy.sum(img[:, :fsize], axis=1)
r_diff =  img[:, fsize:] - img[:, :-fsize]
r_int = numpy.cumsum(numpy.hstack((r_sum.reshape(-1,1), r_diff)), axis=1)

# integrate in the Y direction
c_sum = numpy.sum(r_img[:fsize, :], axis=0)
c_diff = r_img[fsize:, :] - r_img[:-fsize, :]
c_int = numpy.cumsum(numpy.vstack((c_sum, c_diff)), axis=0)

# now we have an array of sums, average can be obtained by division
avg_img = c_int / (f_size * f_size)

该方法返回一张图像，它在两个方向上都比原图小1个像素大小，因此您需要自己处理边框效果。最外层的像素本来就不好，但如果需要，选择正确的边框填充方法是由您决定的。该算法是获得平均值（计算量最少）的最快方法，特别是比numpy.convolve快得多。

如果像上面那样同时对图像及其平方进行平均，则可以使用类似的技巧来计算方差。

npts = fsize * fsize
variance = (rolling_sum(img**2) - rolling_sum(img)/npts) / npts

rolling_sum是一个滑动求和的函数（也就是上面算法中去掉最后一步除法的部分）。因此，只需要两个滑动求和（图像及其平方）即可计算出滚动方差。

（警告：上述代码未经测试，仅用于说明思路。）

- DrV