图像高斯滤波的原理与FPGA实现思路
图像高斯滤波的原理与FPGA实现思路
高斯滤波是图像处理中常用的一种平滑技术，能够有效去除图像中的噪声。本文将详细介绍高斯滤波的原理，并展示如何在FPGA上实现这一算法。
1. 概念
高斯分布（正态分布）是一个常见的连续概率分布，其数学期望值或期望值μ等于位置参数，决定了分布的位置；其方差σ²的开平方或标准差σ等于尺度参数，决定了分布的幅度。正态分布的概率密度函数曲线呈钟形，因此又被称为钟形曲线。标准正态分布是指位置参数μ=0、方差σ²=1的正态分布。
若随机变量X服从一个位置参数为μ、方差为σ²的正态分布，可以记为X~N(μ,σ²)，其概率密度函数为：
$$
g(x) = \frac{1}{\sqrt{2\pi} \sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}} \tag{1}
$$
高斯滤波器是一种根据高斯函数的形状来选择权值的线性平滑滤波器，对于抑制服从正态分布的噪声非常有效。一维零均值高斯函数为：
$$
g(x) = \frac{1}{\sqrt{2\pi} \sigma} e^{-\frac{x^2}{2\sigma^2}} \tag{2}
$$
其中，高斯分布参数σ决定了高斯函数的宽度。一维高斯函数的图形如下：
二维高斯分布的函数为：
$$
g(x, y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}} \tag{3}
$$
二维高斯分布的图像为：
2. 高斯滤波性质
高斯函数具有五个重要的性质，这些性质使得它在早期图像处理中特别有用：
- 二维高斯函数具有旋转对称性，即滤波器在各个方向上的平滑程度是相同的。
- 高斯函数是单值函数，这意味着高斯滤波器用像素邻域的加权均值来代替该点的像素值，而每一邻域像素点权值是随该点与中心点的距离单调增减的。
- 高斯函数的傅立叶变换频谱是单瓣的，这意味着平滑图像不会被不需要的高频信号所污染，同时保留了大部分所需信号。
- 高斯滤波器宽度是由参数σ表征的，σ越大，高斯滤波器的频带就越宽，平滑程度就越好。
- 由于高斯函数的可分离性，较大尺寸的高斯滤波器可以得以有效地实现。
3. 高斯滤波的原理与实现
3.1 高斯模板的生成
高斯滤波既能用来过滤高斯噪音，也可用来做高斯模糊。要模糊一张图像，可以直接用均值滤波来做简单的模糊，但是这样做显然不大合理，因为图像是连续的，离卷积核中心的点关系更加密切，越远的点关系越疏远，这个时候就需要加权平均。明显的离中心点越近的像素点权重越大。而正态分布显然是一种可取的权重分配方式，又由于图像是二维的，所以需要使用二维的高斯函数。所以高斯滤波的本质是利用高斯函数来生成高斯核（高斯卷积模板）来对图像进行卷积操作。
理论上高斯分布在所有定义域上都有非负值，这就需要一个无限大的卷积核，但是实际上，仅需要取均值的三倍标准差（即3σ）内的值，以外的部分去掉即可。
高斯滤波最重要的就是找到高斯模板然后进行卷积，以3X3高斯模板为例，假设中心点的坐标为（0，0），根据二维高斯函数g(x,y)，还需要设定σ的值，假定σ=0.8（这个值不宜过大，否则就会变成均值滤波），可以根据坐标值来算出对应的高斯模板。
假定中心点的坐标为（0，0），那么距离它最近的八个点的坐标如下：
将坐标以及设定的σ带入二维高斯函数中，可以得到对应点的坐标的权重如下：
为了防止高斯滤波后的图像偏亮或者偏暗，我们还需要对图像进行归一化，这九个点的权重和等于0.9125991，因此需要分别对这九个数除以0.9125991，最终得到的高斯模板为：
有了高斯模板就可以对图像进行卷积了，但是在FPGA中对于小数的运算不友好，于是我们将这个3X3模板扩大了16倍，得到了近似的整数模板，卷积完成后再除以16来做定浮点数的近似计算。16倍后的模板高斯卷积模板为：
3X3的高斯模板生成与定浮点数的MATLAB算法实现为：
clear;
clc;
close all;
sigma = 0.8;
A = exp(-(1+1)/(2*sigma*sigma))/(2*pi*sigma*sigma);
B = exp(-(1+0)/(2*sigma*sigma))/(2*pi*sigma*sigma);
C = exp(-(0+0)/(2*sigma*sigma))/(2*pi*sigma*sigma);
D = A*4 + B*4 + C;
gauss_double = [A,B,A;B,C,B;A,B,A];
gauss_normal = gauss_double / sum(sum(gauss_double));
gauss_integer = floor(gauss_normal/gauss_normal(1,1));
3.2 高斯算法的FPGA实现
生成卷积模板后的FPGA实现与仿真参考前面几章，这里只贴高斯算法的FPGA实现部分：
module gassin_filter#(
parameter DW = 8
)(
input wire clk ,
input wire rst_n ,
input wire matrix_de ,
input wire [DW-1:0] matrix11 ,
input wire [DW-1:0] matrix12 ,
input wire [DW-1:0] matrix13 ,
input wire [DW-1:0] matrix21 ,
input wire [DW-1:0] matrix22 ,
input wire [DW-1:0] matrix23 ,
input wire [DW-1:0] matrix31 ,
input wire [DW-1:0] matrix32 ,
input wire [DW-1:0] matrix33 ,
output wire gassin_data_de ,
output wire [DW-1:0] gassin_data
);
//gassin_filter
// [1 2 1]
// [2 4 2]
// [1 2 1]
reg [1:0] matrix_de_r ;
reg [DW+2:0] one_line ;
reg [DW+2:0] two_line ;
reg [DW+2:0] three_line ;
reg [DW+4:0] sum_matrix ;
always @(posedge clk)begin
if(rst_n == 0)begin
matrix_de_r <= 0;
end
else begin
matrix_de_r <= {matrix_de_r[0],matrix_de};
end
end
always @(posedge clk)begin
if(rst_n==0)begin
one_line <= 0 ;
two_line <= 0 ;
three_line <= 0 ;
end
else if(matrix_de)begin
one_line <= matrix11 + 2*matrix12 + matrix13 ;
two_line <= 2*matrix21 + 4*matrix22 + 2*matrix23 ;
three_line <= matrix31 + 2*matrix32 + matrix33 ;
end
else begin
one_line <= 0 ;
two_line <= 0 ;
three_line <= 0 ;
end
end
always @(posedge clk)begin
if(rst_n == 0)begin
sum_matrix <= 0;
end
else if(matrix_de_r[0])begin
sum_matrix <= one_line + two_line + three_line ;
end
else begin
sum_matrix <= 0;
end
end
assign gassin_data = sum_matrix[DW+4:4] ;
assign gassin_data_de = matrix_de_r[1] ;
endmodule
3.3 高斯算法的MATLAB实现以及验证
clear;
clc;
close all;
sigma = 0.8;
A = exp(-(1+1)/(2*sigma*sigma))/(2*pi*sigma*sigma);
B = exp(-(1+0)/(2*sigma*sigma))/(2*pi*sigma*sigma);
C = exp(-(0+0)/(2*sigma*sigma))/(2*pi*sigma*sigma);
D = A*4 + B*4 + C;
gauss_double = [A,B,A;B,C,B;A,B,A];
gauss_normal = gauss_double / sum(sum(gauss_double));
gauss_integer = floor(gauss_normal/gauss_normal(1,1));
GRAY = imread('../img/gray.bmp');
[row,col] = size(GRAY);
gassin_padding = zeros(row+2,col+2);
gassin_result = zeros(row,col);
for i = 1:row
for j = 1:col
gassin_padding(i+1,j+1) = GRAY(i,j);
end
end
for i = 1:row+2
gassin_padding(i,1) = gassin_padding(i,2);
gassin_padding(i,col+2) = gassin_padding(i,col+1);
end
for i = 1:col+2
gassin_padding(1,i) = gassin_padding(2,i);
gassin_padding(row+2,i) = gassin_padding(row+1,i);
end
for i = 2:row+1
for j = 2:col+1
matrix11 = gassin_padding(i-1,j-1);
matrix12 = gassin_padding(i-1,j);
matrix13 = gassin_padding(i-1,j+1);
matrix21 = gassin_padding(i,j-1);
matrix22 = gassin_padding(i,j);
matrix23 = gassin_padding(i,j+1);
matrix31 = gassin_padding(i+1,j-1);
matrix32 = gassin_padding(i+1,j);
matrix33 = gassin_padding(i+1,j+1);
matrix = [matrix11,matrix12,matrix13;matrix21,matrix22,matrix23;matrix31,matrix32,matrix33];
gassin_mult = matrix.* gauss_integer;
sum_gassin_matrix = sum(sum(gassin_mult()));
gassin_result(i-1,j-1) = sum_gassin_matrix/16;
end
end
a = textread('../data/gassin_filter.txt','%s');
IMdec1 = hex2dec(a);
IM1 = reshape(IMdec1,col,row);
fpga_Y = uint8(IM1)';
b = textread('../data/pre.txt','%s');
subplot(1,3,1)
matlab_Y = uint8(floor(gassin_result));
imshow(matlab_Y),title('MATLAB gassin算法图像');
subplot(1,3,2)
imshow(fpga_Y),title('FPGA gassin算法图像');
subplot(1,3,3)
imshow(GRAY),title('原图像');
sub = matlab_Y - fpga_Y;
min_sub = min(min(sub));
max_sub = max(max(sub));