卷积运算以及高斯滤波器的构造

卷积的数学定义

在图像分析和图像处理中，卷积(convolution)是一种非常重要的运算。卷积是一个积分运算，其反应的是函数f(x)f(x)f(x)在另一个函数h(x)h(x)h(x)上移动时所叠加的量。函数fff和hhh在有限域[0,t][0,t][0,t]上的一维卷积为：

(f∗h)(t)=∫0tf(τ)h(t−τ)dτ=∫0tf(t−τ)h(τ)dτ(f*h)(t) = \int^{t}_{0}f(\tau)h(t-\tau)\mathrm{d}\tau = \int^{t}_{0}f(t - \tau)h(\tau)\mathrm{d}\tau (f∗h)(t)=∫0t​f(τ)h(t−τ)dτ=∫0t​f(t−τ)h(τ)dτ

需要注意的是，卷积积分的上下限实际为(−∞,+∞)(-\infty, +\infty)(−∞,+∞)，但是此处我们假设负坐标部分的值为0，因此这里可以限定在区间[0,t][0,t][0,t]中。[1]

根据一维卷积的定义，我们可以得到函数f(x,y)f(x,y)f(x,y)在另一个函数h(x,y)h(x,y)h(x,y)上移动时的卷积：

(f∗h)(x,y)=∫−∞+∞f(a,b)h(x−a,y−b)dadb=∫−∞+∞f(x−a,y−b)h(a,b)dadb\begin{aligned} (f*h)(x,y) &= \int^{+\infty}_{-\infty}f(a,b)h(x-a,y-b)\mathrm{d}a\mathrm{d}b \\ &= \int^{+\infty}_{-\infty}f(x-a,y-b)h(a,b)\mathrm{d}a\mathrm{d}b \end{aligned} (f∗h)(x,y)​=∫−∞+∞​f(a,b)h(x−a,y−b)dadb=∫−∞+∞​f(x−a,y−b)h(a,b)dadb​

数字图像分析中的卷积

数字图像是一个二维的离散数据，因此对于数字图像而言，需要将积分运算（∫\int∫)改为加和运算（∑\sum∑），从而得到离散卷积(discrete convolution)。

对于数字图像而言，在图像平面上存在有限的定义域，对于定义域之外的离散卷积结果为0。因此，这个特点并不妨碍我们在数字图像上执行卷积运算。此时，卷积表示的是使用滤波器hhh对图像执行一个线性滤波处理。

我们在图像f(i,j)f(i,j)f(i,j)的一个局部邻域O\mathcal{O}O内计算其所有像素的线性组合，邻域O\mathcal{O}O中的各个像素的权重由h(i,j)h(i,j)h(i,j)来确定。因此，对于图像中的像素f(i,j)f(i,j)f(i,j)，我们在其邻域O\mathcal{O}O内使用滤波器hhh对其执行卷积运算，得到g(i,j)g(i,j)g(i,j)：

g(i,j)=∑(m,n)∈Of(i−m,j−n)h(m,n)g(i,j)=\displaystyle\sum_{(m,n)\in\mathcal{O}}f(i-m,j-n)h(m,n) g(i,j)=(m,n)∈O∑​f(i−m,j−n)h(m,n)

在上式中，称hhh为卷积掩膜(convolution mask)或者滤波器，并且一般会选择具有奇数行和列的矩形邻域O\mathcal{O}O，以方便确定领域的中心。[1]

利用高斯滤波器来模糊图像

通常，图像处理软件会提供“模糊”（blur）滤镜，使图片产生模糊的效果。具体如下图所示：

模糊的本质就是利用像素邻域O\mathcal{O}O内其他像素的线性运算来计算该像素的值，从而提高该像素和邻域内其他像素之间的相关性。

模糊的方式和效果取决于线性运算时的滤波器hhh的选取。hhh的选择有很多，其中有一种就是高斯滤波器，其对应的模糊效果称之为“高斯模糊”。“高斯模糊”的本质就是利用高斯函数生成线性运算时各像素的权重，即利用高斯函数（当μ=0\mu=0μ=0且各方向的σ\sigmaσ相等时）生成滤波器hhh：

G(x,y)=12πσ2e−x2+y22σ2G(x,y) = \frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}} G(x,y)=2πσ21​e−2σ2x2+y2​

hhh中心点的坐标为(0,0)(0,0)(0,0)，则整个hhh的坐标如下所示：

利用G(x,y)G(x,y)G(x,y)的公式，即可计算得到σ\sigmaσ为指定值时的高斯滤波器hhh。具体细节此处不再赘述，具体可以参考：高斯模糊的算法。

不依赖其他库生成高斯滤波器

当μ=0\mu=0μ=0且各方向的σ\sigmaσ相等时，对于h(i,j)h(i,j)h(i,j)我们得到：

h(i,j)=G(i,j)∑i,jG(i,j)=12πσ2e−(i2+j2)/(2σ2)∑i,j12πσ2e−(i2+j2)/(2σ2)=12πσ2e−(i2+j2)/(2σ2)12πσ2∑i,je−(i2+j2)/(2σ2)=e−(i2+j2)/(2σ2)∑i,je−(i2+j2)/(2σ2)\begin{aligned}

h(i,j)&=\frac{G(i,j)}{\displaystyle\sum_{i,j}{G(i,j)}} \\ &=\frac{\frac{1}{2\pi\sigma^2}e^{-(i^2+j^2)/(2\sigma^2)}}{\displaystyle\sum_{i,j}{\frac{1}{2\pi\sigma^2}e^{-(i^2+j^2)/(2\sigma^2)}}} \\ &=\frac{\frac{1}{2\pi\sigma^2}e^{-(i^2+j^2)/(2\sigma^2)}}{\frac{1}{2\pi\sigma^2}\displaystyle\sum_{i,j}e^{-(i^2+j^2)/(2\sigma^2)}} \\ &=\frac{e^{-(i^2+j^2)/(2\sigma^2)}}{\displaystyle\sum_{i,j}e^{-(i^2+j^2)/(2\sigma^2)}}

\end{aligned} h(i,j)​=i,j∑​G(i,j)G(i,j)​=i,j∑​2πσ21​e−(i2+j2)/(2σ2)2πσ21​e−(i2+j2)/(2σ2)​=2πσ21​i,j∑​e−(i2+j2)/(2σ2)2πσ21​e−(i2+j2)/(2σ2)​=i,j∑​e−(i2+j2)/(2σ2)e−(i2+j2)/(2σ2)​​

使用python得到二维的高斯滤波器的代码如下：

def gaussian_kernel_2d(sigma, width):
    if width == 1:
        return np.ones((1, 1))

kernel_radius = np.floor(width >> 1)
    ax            = np.arange(-kernel_radius, kernel_radius + 1., dtype=np.float32)
    xx, yy        = np.meshgrid(ax, ax)
    kernel        = np.exp(-(xx**2 + yy**2) / (2. * sigma**2))
    return kernel / np.sum(kernel)

详细代码可以参考高斯滤波器的生成和可视化。

滤波器的可分离特性

图像处理中的可分离滤波器(separable filter)可以写成两个更简单的滤波器的乘积。通常，可以把二维卷积运算分离为两个一维滤波器，从而降低图像的卷积运算复杂度。

例如，对于W×HW \times HW×H的图像，如果采用n×nn \times nn×n的二维滤波器，则卷积运算的复杂度为O(M⋅N⋅n2)O(M·N·n^2)O(M⋅N⋅n2)，而如果用两个1×n1 \times n1×n的一维滤波器，则卷积运算的复杂度为O(2M⋅N⋅n)O(2M·N·n)O(2M⋅N⋅n)。[4]

Karas Pavel和Svoboda David在Algorithms for Efficient Computation of Convolution[5]的3. Separable convolution部分指出：

给定一个行向量u→=(u1,u2,...um)\overrightarrow{u}=(u_1,u_2,...u_m)u=(u1​,u2​,...um​)和一个列向量v→T=(v1,v2,...,vn)\overrightarrow{v}^T=(v_1,v_2,...,v_n)vT=(v1​,v2​,...,vn​)，则其二者的卷积定义为：

u→∗v→=v→u→=An×m\begin{aligned} \overrightarrow{u}*\overrightarrow{v}&=\overrightarrow{v}\ \overrightarrow{u} \\ &=\mathbf{A}_{n \times m} \end{aligned}

u∗v​=v u=An×m​​

并且，当且仅当矩阵An×m\mathbf{A}_{n \times m}An×m​的秩r(A)=1r(\mathbf{A})=1r(A)=1时，滤波器才是可分离的滤波器。Gaussian滤波器和Sobel滤波器就是这种可分离的滤波器。

利用这个特性，对于二维高斯滤波器h(x,y)h(x,y)h(x,y)而言，我们将其分离成两个一维滤波器。对图像的二维卷积运算从而简化为对图像的两次一维卷积，具体如下所示：

h(x,y)=(hx∗hy)(x)(f∗h)(x,y)=f∗(hx∗hy)=(f∗hx)∗fy\begin{aligned} h(x,y) &= (h_x * h_y)(x) \\ (f*h)(x,y) &= f*(h_x * h_y) \\ &=(f*h_x)*f_y \end{aligned} h(x,y)(f∗h)(x,y)​=(hx​∗hy​)(x)=f∗(hx​∗hy​)=(f∗hx​)∗fy​​

同时，利用如上的特性，我们还可以利用一维高斯滤波器的结果生成二维高斯滤波器，代码如下所示：

def gaussian_kernel_2d_acc(self):
    if self.width == 1:
        self.kernel = np.ones((1, 1))
        return

self.gaussian_kernel_1d()
    k = self.kernel
    self.kernel = k.reshape(self.width,1) * k.reshape(1,self.width)

更详细的代码可以参考高斯滤波器的生成。

[1]: 图像处理、分析与机器视觉（第4版）
[2]: 如何通俗易懂的解释卷积？
[3]: 高斯模糊的算法
[4]: separable filter
[5]: Karas Pavel and Svoboda David (January 16th 2013). Algorithms for Efficient Computation of Convolution, Design and Architectures for Digital Signal Processing, Gustavo Ruiz and Juan A. Michell, IntechOpen, DOI: 10.5772/51942. Available from: https://www.intechopen.com/books/design-and-architectures-for-digital-signal-processing/algorithms-for-efficient-computation-of-convolution