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

Question

鉴于我已f(x,y)加载图像，例如，

\n\n

\n\n

我想计算\xe2\x88\x82/\xe2\x88\x82x \xe2\x88\x82/\xe2\x88\x82y G*f图像 的高斯导数f，其中G是高斯滤波器，*表示卷积。使用 Scipy 可以轻松完成此操作：

\n\n

from scipy.ndimage.filters import gaussian_filter\nimshow(gaussian_filter(g, sigma, order=1))\n

\n\n

这样sigma=50会产生以下结果：

\n\n

\n\n

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

\n\n

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

\n\n

尽管如此，结果看起来还是合理的：

\n\n

\n\n

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

\n\n

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

\n\n

\n\n

现在这是出乎意料的！这个结果没有任何意义，不是吗？

\n\n

作为记录，我还检查了偏微分\xe2\x88\x82/\xe2\x88\x82x G*f和\xe2\x88\x82/\xe2\x88\x82y G*f。然而

\n\n

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

\n\n

看起来很合理

\n\n

\n\n

另一个

\n\n

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

\n\n

才不是：

\n\n

\n\n

这是为什么？

\n

Answer 1

gaussian_filter当 和order均cval非零时，会出现一个错误。具体来说，它在这里：

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

滤波器执行重复的一维卷积，并且每次都传入cval一维滤波器。问题是，如果有任何导数，则应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()

回报

这是预期的结果。（我的原始图像与您的有点不同，并且我使用了不同的可视化工具。）