使用 Scipy 进行高斯滤波的意外行为

Question

假设我加载了一个图像f(x,y)，例如，

我想计算图像f 的高斯导数∂/∂x ∂/∂y G*f，其中G 是高斯滤波器，* 表示卷积。使用 Scipy 很容易做到这一点：

from scipy.ndimage.filters import gaussian_filter
imshow(gaussian_filter(g, sigma, order=1))

使用sigma=50 会产生以下结果：

现在，出于应用原因，我需要使用 mode='constant' 进行计算：

imshow(gaussian_filter(g, sigma, order=1, mode='constant', cval=0))

不过，结果看起来还是合理的：

但是，请注意我图像的背景强度是1，而不是0。因此，使用cval=1应该是合理的：

imshow(gaussian_filter(g, sigma, order=1, mode='constant', cval=1))

现在这是出乎意料的！这个结果没有意义吧？

为了记录，我还检查了偏微分 ∂/∂x G*f 和 ∂/∂y G*f。而

imshow(gaussian_filter(g, sigma, order=[0, 1], mode='constant', cval=1)

看起来很合理

另一个

imshow(gaussian_filter(g, sigma, order=[1, 0], mode='constant', cval=1)

没有：

这是为什么呢？

Answer 1

gaussian_filter 中存在一个错误，当order 和cval 均非零时会出现该错误。 Specifically, it's here：

for axis, sigma, order, mode in axes:
    gaussian_filter1d(input, sigma, axis, order, output, mode, cval, truncate)
    input = output

filter执行重复的1d卷积，每次都传入cval到1d filter。问题是，如果有任何导数，那么 cval 应该设置为 0，因为任何常数的导数都是零。这就是为什么order=[1, 0] 的结果是错误的，而order=[0, 1] 的结果却不是。未经测试（没有 SciPy 开发环境），我认为以下是正确的：

for axis, sigma, order, mode in axes:
    gaussian_filter1d(input, sigma, axis, order, output, mode, cval, truncate)
    if order > 0: 
        cval = 0.0
    input = output

解决方法

非零cval 可以通过在过滤之前从图像中减去它来模拟（并且只有在阶数为零时才在过滤后加回来）。示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage.filters import gaussian_filter

g = np.ones((500, 500))
g[200:300, 200:300] = 2
sigma = 50
cval = 1
gf = gaussian_filter(g-cval, sigma, order=1, mode='constant')
plt.matshow(gf)
plt.show()

返回

这是预期的结果。 （我的原图和你的有点不一样，我用了不同的可视化工具。）