卷积与滤波器

前言

开完题了,论文的审稿意见也下来了。不出意料,有两个评委给了 No,一个给了 Yes。值得开心的是他们都说这是一个新颖的想法,不足还是自己的写作方面,没有把模型说清楚,而且实验做的也不够充分,接下来可能就是要慢慢改了。循环神经网络部分的实验内容也总结完了,虽然还是模模糊糊,但是这些东西大体上在我看来也没那么高深莫测,倒是一些细节方面确实挺难的。

卷积运算

在学习卷积神经网络之前需要再好好巩固一下什么是卷积运算,然后顺便把大四学习的滤波器的内容总结一下。知乎上有不少通俗易懂的解释可以参考通俗易懂地理解卷积,其实卷积是之前在优化算法中提到的移动平均的推广。

定义

$f$ 和 $g$ 的卷积写为 $(f*g)(n)$,其离散的定义为:
$$
(f*g)[n]=\sum_{\tau=-\infty}^\infty f[\tau]g[n-\tau]
$$
连续的定义为:

$$
(f*g)(n)=\int_{-\infty}^\infty f(\tau)g(n-\tau)d\tau
$$
其中 $n=\tau+(n-\tau)$,借鉴一下马同学在知乎中的例子,离散卷积的应用场景:两个普通骰子点数加起来等于四的概率。
$$
f[1]g[3]+f[2]g[2]+f[3]g[1]
$$
符合卷积的定义,写成标准形式就是:
$$
(f*g)[4]=\sum_{m=1}^3f[4-m]g[m]
$$
连续卷积的应用场景:追踪飞船的位置 $f(t)$,由于噪声的影响,我们需要对得到的结果进行加权平均,时间上越近的测量结果越相关,相关函数为 $g(t)$,所以对飞船的位置的估计为:
$$
\int_0^T f(t)g(T-t)dt
$$
这就是连续的卷积,也就是将信号 $f(t)$ 和翻转平移后的信号 $g(t)$ 进行积分。在卷积神经网络中,卷积的第一个参数 $f$ 通常叫做输入,第二个参数 $g$ 通常叫做核函数,输出有时候叫特征映射。

计算离散卷积

计算离散卷积 $f[n]*g[n]$ 一般有三种方法:

  1. 直接计算
  2. 快速傅里叶变换
  3. 分段卷积

第一种方法就是直接根据定义来计算,第二种和第三种方法则都用到了快速傅里叶变换的知识,第三种方法是先将信号分成一小段一小段再以第二种方法来计算,时间复杂度会小一点,这里就只介绍第二种方法。三种方法的时间复杂度可以参考维基百科

快速傅里叶变换

两个离散信号在时域(time domain)做卷积相当于这两个信号在频域(frequence domain)相乘。

$$
y[t]=f[t]*g[t]\leftrightarrow Y[f]=F[f]G[f]
$$

根据定义不难证明以上公式,由于乘法比较简单,因此在计算卷积的时候先将信号从时域转成频域,然后计算 $Y[f]$ ,在将计算结果用傅里叶变换的逆变换转回时域。

滤波器

卷积运算通常用于图像处理领域中的滤波操作,首先先来看看数字图像处理领域的女神 Lena,有兴趣的可以点击查看原图。由于这张图包含了各种细节、平滑区域、阴影和纹理,因此广泛用于图像处理。

这是一张 400*400 三通道的图片,为了更方便演示算法,我们通常在灰度图上进行操作(Matlab 中可以使用 rgb2gray() 函数获取图像的灰度值)。那么如何把灰度图和波联系起来呢?我们取灰度图的第一行像素,然后根据灰度值画出曲线就能得到波:

图像的频率反映了图像的像素灰度在空间中变化的情况,是灰度在平面空间上的梯度。

图中曲线波动较大代表灰度图明暗差距大的地方,即高频成分较强,低频成分较弱;而波动较小代表灰度图明暗差距较小的地方,即低频成分较强,高频成分较弱。

这段话虽然能理解,但是学了傅里叶变换之后感觉懵懵的,总是觉得这一整个波不就是由同样的正弦波组成的吗?后来才醒悟过来,所谓的高频低频是相对的。先举个例子两张图片,一面墙壁的图像和国际象棋棋盘,前者灰度值分布平坦,其低频成分就较强,而高频成分较弱;而后者具有快速空间变化的图像来说,其高频成分会相对较强,低频则较弱。那么在一张图片中怎么理解呢?那就是将需要比较的两个地方分别截取出来,然后进行周期性延拓,因此波动较大的地方截取出来周期性延拓,进行傅里叶变换后,得到一系列正弦波,高频的正弦波就比较强(相比较于时域比较平坦的区域)。

滤波器根据滤波的目的可以分为两种,通过低频的滤波器称为低通滤波器,通过高频的滤波器称为高通滤波器。实现滤波的方式也分为两种,空间滤波和频域滤波。

空间滤波器

空间滤波器指对图像一个邻域内的像素执行预定义的操作,滤波产生一个新像素,它的坐标等于邻域中心的坐标,滤波器的中心访问图像的每一个像素后就生成了滤波后的图像。通常会将滤波的结果存在新的图像中,避免改变图像内容的同时进行滤波操作。通常使用奇数尺寸 $3\times 3$ 的滤波器,对一张 $4\times 4$的图像进行滤波,第一行第一个像素和第三行第二个像素的滤波过程如下图所示:

图像边界的像素点没有那么多邻域,通常需要零填充(如下图所示)或者其他操作,然后对中间的像素进行滤波操作。

根据上面提到的预定义的操作,空间滤波分为线性空间滤波和非线性空间滤波。非线性空间滤波例如统计排序滤波器(中值、最大值和最小值)和自适应滤波器等等。举个例子,最大值滤波器就是输出邻域内的像素的灰度值的最大值,如下图所示:

$7=max(4, 7, 7, 2, 2, 5, 2, 4, 3)$,就是说新的图像的灰度值为原图像邻域内像素的灰度值的最大值。线性空间滤波器有一个和邻域对应的滤波器系数矩阵 $\omega$,在滤波过程中,每个邻域内的像素的灰度值分别与对应位置的系数相乘,最后再相加,如下图所示:

对于图像 $f$ 中任意一点 $(x, y)$,滤波器的响应 $g(x, y)$ 是滤波器系数与该滤波器所包围的像素的灰度值的乘积之和:
$$
g(x, y)=\omega(-1, -1)f(x-1, y-1)+\omega(-1, 0)f(x-1, y)…\omega(1, 1)f(x+1, y+1)
$$
对于奇数尺寸 $(2a+1)\times (2b+1)$ 的滤波器,有:
$$
g(x, y)=\sum_{s=-a}^a\sum_{t=-b}^b\omega(s, t)f(x+s, y+t)
$$
这个操作叫做相关,有点类似于二维卷积操作,不过卷积需要滤波器翻转:
$$
\omega(x, y)*f(x, y)=\sum_{s=-a}^a\sum_{t=-b}^b\omega(s, t)f(x-s, y-t)
$$
所以如果滤波器是对称的,那么相关和卷积将得到相同的结果,我们将翻转后的滤波器系数和图像进行卷积即可实现线性空间滤波操作。

边缘提取

大四选修数字图像处理的实验三就要求实现四种(Sobel 算子、Prewitt 算子、Roberts 算子和 Marr 算子)边缘提取函数(空间高通滤波),代码也都托管在 Github。不同的边缘提取函数的区别就是算子的不同,即滤波器系数矩阵不同,这里只详细介绍一下 Sobel 算子:

用这两个算子与原图像做卷积运算可以分别得到纵向和横向灰度值的梯度(即变化率),像素两边的灰度值差距越大,卷积的结果的绝对值也就越大,因此可以得到图像的边缘:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
[img, map] = imread('Lena.jpg');
img = rgb2gray(img);

subplot(131), imshow(img);

hx = [-1 -2 -1; 0 0 0; 1 2 1];
hy = hx';
Gx = conv2(im2double(img), hx);
Gy = conv2(im2double(img), hy);
sobel = abs(Gx) + abs(Gy);
subplot(132), imshow(sobel);

threshold = max(max(sobel)) * 0.18;
sobel(sobel >= threshold) = 1;
sobel(sobel < threshold) = 0;

subplot(133), imshow(sobel);

得到横向和纵向的梯度后,使用以下公式计算灰度的大小(注意是矩阵元素的平方 sqrt(Gx.^2+Gy.^2)):
$$
G=\sqrt{G_x^2+G_y^2}
$$
但是为了提高运算效率,通常使用不开方的近似值,即使会损失一定的精度:
$$
|G|=|G_x|+|G_y|
$$
最后可以进行阈值处理,将灰度图变成二值图像,让检测的边缘更加明显。其他的算子如下所示:

  • Roberts 算子:hx = [-1 -1 -1; 0 0 0; 1 1 1]; hy = hx';
  • Prewitt 算子:hx = [-1 0; 0 1]; hy = [0 -1; 1 0];

频域滤波器

实验二要求实现三种(理想、布特沃斯和高斯滤波)低通滤波器和高通滤波器,由于时域上的卷积操作等于频域上的乘积操作,因此线性空间滤波和频域滤波一一对应,但是频域滤波不能进行非线性滤波。

频谱图

傅里叶变换可以将信号从时域变换成频域,在数字图行处理中时域体现为空间上的分布,因此傅里叶变换可以将图像的灰度分布函数变成图像的频率分布函数。在做傅里叶变换的时候相当于是先将图片做周期性延拓,然后再进行傅里叶变换,如下图所示,通常我们只取一个周期分析(即红色框内),越亮表示该频率的幅值越大。

由于图像是一个实数信号,可以在数学上严格证明,它的傅里叶变换后的频域信号关于 0 频轴对称。所以上图右边的红框中四个角为低频,中间为高频,而且低频幅值较大。为了便于频域的滤波和频谱的分析,常常在变换之前进行频谱的中心化。如下图所示,离中心越近,频率越低,离中心越远,频率越高,类似于二维的坐标轴。

1
2
3
4
5
[img, map] = imread('Lena.jpg');
img = rgb2gray(img);
subplot(131), imshow(img);
subplot(132), f = fft2(img), imshow(log(1+abs(f)), []);
subplot(133), f = fftshift(f), imshow(log(1+abs(f)), []);

由于幅度值范围很大,所以在显示频谱图之前要取对数处理;imshow(I,[])对取对数后的值进行归一化,自动拉伸动态范围,增大对比度; fftshift() 将低频调整到中间,高频在外围。

理想低通滤波器

这里就只在详细介绍理想低通滤波器,它由下面的函数确定:
$$
\begin{align*}
H(u, v)=\begin{cases}
1 & D(u, v) \leq D_0 \\
0 & D(u, v) > D_0
\end{cases}
\end{align*}
$$
其中 $D_0$ 是正常数,$D(u, v)$ 是频域中的点 $(u, v)$ 和中心的距离,即
$$
D(u, v)=\sqrt{(u-P/2)^2+(v-Q/2)^2}
$$
$P, Q$ 是频域图的尺寸。“理想”的意思是在以 0 频率为中心,半径为 $D_0$ 的圆内,所有频率无衰减地通过,圆外地所有频率则被完全过滤。滤波代码和效果如下所示:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
[img, map] = imread('Lena.jpg');
img = rgb2gray(img);
D0 = 32

subplot(221), imshow(img);
subplot(222), f = fftshift(fft2(img)), imshow(1+log(abs(f)), []);

[P, Q] = size(f);
H = zeros(P, Q);
for u = 1:P
for v = 1:Q
D = sqrt((u - round(P/2))^2 + (v - round(Q/2))^2);
if D <= D0
H(u, v) = 1;
else
H(u, v) = 0;
end;
end;
end;
f = f .* H;

img = uint8(real(ifft2(ifftshift(f))));
subplot(223), imshow(img);
subplot(224), f = fftshift(fft2(img)), imshow(log(abs(f)), []);

figure, surf(1:P, 1:Q, H);
axis([1 m 1 n 0 1]);

对于一幅图像来说,低频就是边缘以内的内容,也就是图像的大部分信息,即图像的大致概貌和轮廓,是图像的近似信息。高频指图像边缘的灰度值变化快,图像的细节处也是属于灰度值急剧变化的区域,正是因为灰度值的急剧变化,才会出现细节。另外噪点也是因为像素点灰度值明显不一样了,即高频部分,因此噪声在高频。因此低通滤波器可以去噪,高通滤波器可以去模糊,各种滤波器对应的公式如下所示,n 阶布特沃斯滤波器的阶数趋于正无穷时,就是理想滤波器,对频率的截止也就会更加尖锐。

滤波器函数
理想低通滤波器$\begin{equation}H(u, v)=\begin{cases}1 & D(u, v) \leq D_0 \\ 0 & D(u, v) > D_0\end{cases}\end{equation}$
n 阶布特沃斯低通滤波器$H(u,v)=\frac{1}{1+\Big[D(u, v)/D_0\Big]^{2n}}$
高斯低通滤波器$H(u, v)=e^{-D^2(u, v)/2D_0^2}$
理想高通滤波器$\begin{equation}H(u, v)=\begin{cases}0 & D(u, v) \leq D_0 \\ 1 & D(u, v) > D_0\end{cases}\end{equation}$
n 阶布特沃斯高通滤波器$H(u,v)=\frac{1}{1+\Big[D_0/D(u, v)\Big]^{2n}}$
高斯高通滤波器$H(u, v)=1-e^{-D^2(u, v)/2D_0^2}$

参考文献

  1. ViatorSun. 图像傅里叶变换的频谱图. https://blog.csdn.net/ViatorSun/article/details/82387854
  2. 冈萨雷斯. 数字图像处理
疏影横斜水清浅,暗香浮动月黄昏