图像处理算法 | 高斯滤波

1. 高斯滤波原理

高斯滤波器的本质是一个二维低通系统，它通常是以一个矩阵的形式出现，比如式（1）这个5×5高斯核（又称为kernel），用它对图像做卷积，图像的高频分量会被打到低频，最终呈现的效果就是图像变模糊了。
$$gaus\_kernel=\left[\begin{array}{ccccc}4& 16& 24& 16& 4\\ 16& 64& 96& 64& 16\\ 24& 96& 144& 96& 24\\ 16& 64& 96& 64& 16\\ 4& 16& 24& 16& 4\end{array}\right]\text{\hspace{1em}(1)}$$
生成它的是一个以e为底的指数函数，通常用exp来标识，如（2）所示：
$$f(x,y)=\frac{1}{2\pi {\sigma }^2}exp(-\frac{x^2+y^2}{2{\sigma }^2})=\frac{1}{2\pi {\sigma }^2}{e}^{(-\frac{x^2+y^2}{2{\sigma }^2})}(2)$$
公式（2）就是常见的二维高斯函数，它是一个和距离有关的函数，以5×5的窗口大小为例，中心点作为坐标原点，x和y分别代表横纵坐标，其坐标如下图1（a）所示，x2+y2是当前坐标到圆心的欧式距离（高斯核的旋转对称性由此而来），x2+y2的结果如下图1（b）所示，离中心点越远，x2+y2的值越大。

考虑到exp里面有个负号，所以公式（2）等价为：
$$f(x,y)=\frac{1}{2\pi {\sigma }^2}exp(-\frac{x^2+y^2}{2{\sigma }^2})=\frac{1}{2\pi {\sigma }^2}(\frac{1}{e}{)}^{(\frac{x^2+y^2}{2{\sigma }^2})}(3)$$
0<1/e<1，所以f(x,y)是个单调减函数，那么x2+y2的值越大，求出的值就会越小，所以高斯核呈现出中心点权重最大，依次向四周递减的形状。

公式里还有一个变量σ（方差），方差的物理意义是描述一组数据的分散程度，方差越小，数据就越聚拢；方差越大，数据就越分散。对高斯公式来说，σ越小，窗口内权重越聚拢，得到的图像就越清晰；σ越大，窗口内权重越分散，得到的图像就越模糊。可以换个角度来理解它，如果窗口内的能量值固定为1，分配的越集中，每个位置能分到的能量就越多，高斯的波峰就会越高；分配的越分散，每个位置能分到的能量就越少，高斯的波峰就会越低，所以高斯的极限情况就是均值，所有的位置分配到了相同的能量。

为了进一步观察窗口大小和σ之间的关系，用matlab写一段代码做个实验，分别配置不同的窗口大小win_m、win_n和sigma，显示出16个gauss kernel的三维分布图，如图2所示：

图2有一个需要注意的地方，以sigma=1这列的三维图为例，直观感受会觉得窗口越小，高斯权重越分散，窗口越大，高斯权重越集中，好像在sigma一样的情况下，5×5窗口对图像的模糊力度似乎要比17×17的更佳？其实不是的，这里要注意三维图的坐标，把sigma=1时5×5和17×17的三维图单独放大一下，如下图3所示，右图的17×17其实也只显示了5×5的范围，超过5×5的部分差不多被抑制为0了，所以理论上这两个窗口的模糊力度应该差不多。

再看一下sigma=10的情况，同样是5×5和17×17的对比，把图放大如下图4所示，这个就更明显了，直观上好像5×5比17×17更分散，但是看坐标轴就能发现，17×17的范围显然更大，而且z轴显示出的归一化权重也说明了这个问题，和5×5相比，17×17的权重更小一些，在窗口能量为1的情况下，窗口越大，每个位置能分到的能量也就越少，所以17×17的波峰要比5×5低。

从以上的分析可以看出，如果选择了一个较大的窗口，要把sigma配置的大一些，才能发挥出大窗口的模糊效果，而且大窗口相对小窗口的调试范围更大，可以满足多种需求。

2. 图像处理效果与对比

更进一步的，用以上16个kernel分别对下面的Lena图进行高斯滤波（代码在附录中，需要可自行下载），做个滤波效果的对比。从图5上看，sigma=10且窗口17×17的模糊效果最好；sigma=1时，受sigma分散度的限制，即使是大一些的窗口，也没能发挥出更模糊的滤波效果，这很符合上面对高斯公式的分析。由此也能看出，如果对kernel的性质掌握的足够，是可以通过配置直接判断图像处理效果的。

进一步解释图像高斯滤波，其本质是图像和高斯核的二维卷积，以5×5高斯核为例，按照从上到下从左到右的顺序，每次取图像中一个5×5的窗口数据与高斯核进行卷积（乘加计算），得到的最终结果就是当前中心点的高斯滤波结果。图像四周的边界上会存在缺数据的情况，一般的做法是复制最边界的图像，如果窗口是5×5，那就是图像的上下左右边界各扩边2次。

高斯核给中心点像素分配最大权重，周围的权重离中心点的距离越远其值越小，对于图像细节来说，可以更好的被保留下来，而不是像均值滤波那样被一视同仁的全部抹平。另外，高斯卷积核是实现尺度变换的唯一变换核，并且是唯一的线性核，这个后面如果有需要可以拿出来再详细解释一下。

3. 高斯kernel定点化

用公式（2）生成的5×5高斯核是一堆浮点数据，假设生成的高斯kernel如式（4）所示（已归一化），sigma配置为1.082。
$$gaus\_kernel=\left[\begin{array}{ccccc}0.0046& 0.0166& 0.0255& 0.0166& 0.0046\\ 0.0166& 0.0598& 0.0917& 0.0598& 0.0166\\ 0.0255& 0.0917& 0.1406& 0.0917& 0.0255\\ 0.0166& 0.0598& 0.0917& 0.0598& 0.0166\\ 0.0046& 0.0166& 0.0255& 0.0166& 0.0046\end{array}\right]\text{\hspace{1em}(4)}$$
但是硬件是不能直接处理浮点数据的，需要对式（4）做个定点化，也就是给它乘以一个整数，然后四舍五入取整。为了保证输出图像和输入图像的位宽一致，最后这个乘进来的整数是要在滤波后除掉的（归一化），所以这个整数取2的n次方最合适，此时滤波后只要对结果右移n位就可以了。选择n=10，也就是2^10=1024对gaus_kernel做定点，就可以得到式（1）的结果。

另外二维高斯核还能拆成两个一维高斯核相乘的样子，如式（5）所示，这是个非常优良的性质，在RTL硬件实现的时候可以节省功耗和面积，先写在这里，后面再详细介绍。
$$gaus\_kernel=\left[\begin{array}{ccccc}4& 16& 24& 16& 4\\ 16& 64& 96& 64& 16\\ 24& 96& 144& 96& 24\\ 16& 64& 96& 64& 16\\ 4& 16& 24& 16& 4\end{array}\right]=\left[\begin{array}{c}2\\ 8\\ 12\\ 8\\ 2\end{array}\right]\times \left[\begin{array}{ccccc}2& 8& 12& 8& 2\end{array}\right]\text{\hspace{1em}(5)}$$
附录：

clc;
close all;

fileName      = 'lena';
fmat          = '.bmp';

input_img   = imread([fileName, fmat]);
src         = double(input_img)/256;
[m,n,k]     = size(src);
% figure,imshow(src),title('src');

%% 高斯滤波效果对比
win_m = 8;
win_n = 8;

figure
sigma = 1;
[x,y] = meshgrid(-win_m:1:win_m,-win_n:1:win_n);
gaus_ker  = (1/2*pi*sigma^2)*exp(-(x.^2+y.^2)/(2*sigma^2));
gaus_ker0 = gaus_ker/sum(gaus_ker(:));  % 归一化
surf(gaus_ker0),title('sigma = 1');

figure
sigma = 2;
[x,y] = meshgrid(-win_m:1:win_m,-win_n:1:win_n);
gaus_ker  = (1/2*pi*sigma^2)*exp(-(x.^2+y.^2)/(2*sigma^2));
gaus_ker1 = gaus_ker/sum(gaus_ker(:));  % 归一化
surf(gaus_ker1),title('sigma = 2');

figure
sigma = 5;
[x,y] = meshgrid(-win_m:1:win_m,-win_n:1:win_n);
gaus_ker  = (1/2*pi*sigma^2)*exp(-(x.^2+y.^2)/(2*sigma^2));
gaus_ker2 = gaus_ker/sum(gaus_ker(:));  % 归一化
surf(gaus_ker2),title('sigma = 5');

figure
sigma = 10;
[x,y] = meshgrid(-win_m:1:win_m,-win_n:1:win_n);
gaus_ker  = (1/2*pi*sigma^2)*exp(-(x.^2+y.^2)/(2*sigma^2));
gaus_ker3 = gaus_ker/sum(gaus_ker(:));  % 归一化
surf(gaus_ker3),title('sigma = 10');

gaus_flt0 = imfilter(src, gaus_ker0, 'rep');
gaus_flt1 = imfilter(src, gaus_ker1, 'rep');
gaus_flt2 = imfilter(src, gaus_ker2, 'rep');
gaus_flt3 = imfilter(src, gaus_ker3, 'rep');

figure,imshow(src);
figure,imshow([gaus_flt0,gaus_flt1,gaus_flt2,gaus_flt3]);
imwrite([gaus_flt0,gaus_flt1,gaus_flt2,gaus_flt3],'lena_17.jpg');