高斯与高斯滤波器的卷积

Question

简单的任务..我想用高斯平滑一些向量..这只是一个测试用例，稍后我想将它应用到图像上。

import numpy as np
import scipy.stats
import scipy.ndimage

m = 7  # size of the 'signal'
n = 7  # size of the filter
sgm = 2  # dev for standard distr
weight_conv = np.zeros(2*m*n).reshape(2*n, m)  # Weights for the convolution

input_signal = np.array(range(m))  # input signal..
x1 = np.linspace(-4*sgm, 4*sgm, n)  # x-values for the normal-dstr
input_filter = scipy.stats.norm.pdf(x1, loc=0, scale=sgm)

# create my own weight matrix
for i in range(weight_conv.shape[1]):
    weight_conv[i:(len(input_filter)+i), i] = input_filter

# My own way of calculating the convolution
np.sum(weight_conv * input_signal, axis=1)
# Convolution provided by numpy
np.convolve(input_signal, input_filter)
# Apply the scipy gaussian filter...
scipy.ndimage.filters.gaussian_filter(input_signal, sigma=sgm)
scipy.ndimage.filters.gaussian_filter1d(input_signal, sigma=sgm)

现在我的想法是这些都应该是相似的。我的方法确实产生与 numpy 卷积相似的输出，但 scipy 方法不同......

scipy.ndimage.filters.gaussian_filter(input_signal, sigma=sgm)
array([1, 1, 2, 3, 3, 4, 4])

现在肯定是 scipy 正在做一些不同的事情。但是什么？我不知道。我检查了源代码，似乎他们只是在使用带有高斯核的卷积（这也是我正在做的事情）。但答案并没有加起来......

有人有不同的想法吗？

Answer 1

在@filippo 和SO-question 的帮助下，我能够重建scipy 实现。该方法将信息传播到任一方的方式至关重要。

下面是一个代码 sn-p 显示它们是如何相等的

import numpy as np
import scipy.stats
import scipy.ndimage
import matplotlib.pyplot as plt
np.set_printoptions(linewidth=160)


m_init = 7
# Create any signal here...
input_signal_init = []  
input_signal_init = np.arange(m_init)
input_signal_init = np.random.choice(range(m_init),m_init)
# Convert to float for better results in scipy.
input_signal_init = np.array(input_signal_init).astype(float)

# Simulating method='reflect'
input_signal = np.array([*input_signal_init[::-1], *input_signal_init, *input_signal_init[::-1]])

# Define new length of input signal
m = len(input_signal)  
# Properties of the Gaussian
sgm = 2  # dev for standard distr
radius = 4 * sgm
x = numpy.arange(-radius, radius+1)
n = len(x)
weight_conv = np.zeros(m*(n+m)).reshape(n+m, m)  

# Calculate the gaussian
p = np.polynomial.Polynomial([0, 0, -0.5 / (sgm * sgm)])
input_filter = numpy.exp(p(x), dtype=numpy.double)
input_filter /= input_filter.sum()

# Calculate the filter weights
for i in range(weight_conv.shape[1]):
    weight_conv[i:(len(input_filter)+i), i] = input_filter

# My own way of calculating the convolution
self_conv = np.sum(weight_conv * input_signal, axis=1)[(2*m_init+1):(3*m_init+1)]
# Convolution provided by numpy
numpy_conv = np.convolve(input_signal, input_filter)[(2*m_init+1):(3*m_init+1)]
# Convolution by scipy with method='reflect'
# !! Here we use t[![enter image description here][2]][2]he original 'input_signal_init'
scipy_conv = scipy.ndimage.filters.gaussian_filter(input_signal_init, sigma=sgm)

绘制结果总是让人相信你做得很好......所以

plt.plot(scipy_conv, 'r-')
plt.plot(self_conv, 'bo')
plt.plot(numpy_conv, 'k.-')
plt.show()

给出以下image。

还可以验证，将 scipy 过滤器的mode 设置为'constant' 也会创建相同的卷积。

Answer 2

Scipy 多维高斯滤波器使用更大的内核。默认情况下，内核半径被截断为 4 sigma，在您的情况下，这应该有点类似于 17x17 过滤器。

请参阅_gaussian_kernel1d 了解具体实现。它还使用了几个 1d 可分离的相关性，但这应该没有太大区别。

另一个关键区别在于输出向量的大小和精度。来自 ndimage 文档：

中间数组存储在与输出相同的数据类型中。因此，对于精度较低的输出类型，结果可能不精确，因为存储的中间结果可能精度不足。这可以通过指定更精确的输出类型来防止。

因此，在您的情况下，输出精度受input_signal.dtype 的限制。尝试使用浮点输入数组或不同的数组作为输出。

输出大小和边缘处理有点棘手，不确定是否有办法从 np.convolve 和 scipy.ndimage.gaussian_filter 获得相同的行为