图像高斯滤波的原理及FPGA实现思路

1.概念

高斯分布

图像滤波之高斯滤波介绍

图像处理算法|高斯滤波

高斯滤波（Gaussian filter）包含很多种，包括低通、高通、带通等，在图像上说的高斯滤波通常是指的高斯模糊（Gaussian Blur），是一种高斯低通滤波。通常这个算法也可以用来模糊图像，提供模糊滤镜。也可以用来过滤自然界的高斯白噪声。

高斯分布（正态分布）是一个常见的连续概率分布，正态分布的数学期望值或期望值μ muμ等于位置参数，决定了分布的位置，其方差σ 2 sigma^2σ2的开平方或者标准差σ sigmaσ等于尺度参数，决定了分布的幅度。正态分布的概率密度函数曲线呈钟形，所以又被称为钟形曲线。我们常说的标准正态分布是位置参数μ = 0 , 方差 σ 2 = 1 mu=0,方差sigma^2=1μ=0,方差σ2=1的正态分布。
  若随机变量X XX服从一个位置参数为μ 、方差为 σ 2 mu、方差为sigma^2μ、方差为σ2的正态分布，可以记为X N ( μ , σ 2 ) X~N(mu,sigma^2)XN(μ,σ2),其概率密度函数为：
g ( x ) = 1 2 π σ e ( − ( x − μ ) 2 2 σ 2 ) (1) g(x)=frac{1} {sqrt{2pi} sigma }e^{(-frac{{(x-mu)}^{2} }{2sigma^{2}})} ag{1}g(x)=2πσ1e(−2σ2(x−μ)2)(1)
  高斯滤波器是一种根据高斯函数的形状来选择权值的线性平滑滤波器，对于抑制服从正态分布的噪声非常有效，一维零均值高斯函数为：
g ( x ) = 1 2 π σ e − x 2 2 σ 2 (2) g(x)=frac{1} {sqrt{2pi} sigma }e^{-frac{x^{2} }{2sigma^{2}}} ag{2}g(x)=2πσ1e−2σ2x2(2).
  其中，高斯分布参数（σ sigmaσ）决定了高斯函数的宽度，一维高斯函数的图形如下图：

  二维高斯高斯分布的函数为：
g ( x , y ) = 1 2 π σ 2 e − x 2 + y 2 2 σ 2 (3) g(x,y)=frac{1}{2pisigma^2}e^{frac{-x^2+y^2}{2sigma^2}} ag{3}g(x,y)=2πσ21e2σ2−x2+y2(3)
  二维高斯分布的图像为：

2.高斯滤波性质

高斯函数具有五个重要的性质，这些性质使得它在早期图像处理中特别有用．这些性质表明，高斯平滑滤波器无论在空间域还是在频率域都是十分有效的低通滤波器，且在实际图像处理中得到了工程人员的有效使用．高斯函数具有五个十分重要的性质，它们是：
  1.二维高斯函数具有旋转对称性，即滤波器在各个方向上的平滑程度是相同的．一般来说，一幅图像的边缘方向是事先不知道的，因此，在滤波前是无法确定一个方向上比另一方向上需要更多的平滑．旋转对称性意味着高斯平滑滤波器在后续边缘检测中不会偏向任一方向．
  2.高斯函数是单值函数．这表明，高斯滤波器用像素邻域的加权均值来代替该点的像素值，而每一邻域像素点权值是随该点与中心点的距离单调增减的．这一性质是很重要的，因为边缘是一种图像局部特征，如果平滑运算对离算子中心很远的像素点仍然有很大作用，则平滑运算会使图像失真．
  3.高斯函数的傅立叶变换频谱是单瓣的．正如下面所示，这一性质是高斯函数付立叶变换等于高斯函数本身这一事实的直接推论．图像常被不希望的高频信号所污染(噪声和细纹理)．而所希望的图像特征（如边缘），既含有低频分量，又含有高频分量．高斯函数付立叶变换的单瓣意味着平滑图像不会被不需要的高频信号所污染，同时保留了大部分所需信号．
  4.高斯滤波器宽度(决定着平滑程度)是由参数σ表征的，而且σ和平滑程度的关系是非常简单的．σ越大，高斯滤波器的频带就越宽，平滑程度就越好．通过调节平滑程度参数σ，可在图像特征过分模糊(过平滑)与平滑图像中由于噪声和细纹理所引起的过多的不希望突变量(欠平滑)之间取得折衷．
  5.由于高斯函数的可分离性，较大尺寸的高斯滤波器可以得以有效地实现．二维高斯函数卷积可以分两步来进行，首先将图像与一维高斯函数进行卷积，然后将卷积结果与方向垂直的相同一维高斯函数卷积．因此，二维高斯滤波的计算量随滤波模板宽度成线性增长而不是成平方增长．

3.高斯滤波的原理与实现

3.1 高斯模板的生成

高斯滤波既能用来过滤高斯噪音，也可用来做高斯模糊。
  要模糊一张图像，可以直接用均值滤波来做简单的模糊，但是这样做显然不大合理，因为图像是连续的，离卷积核中心的点关系更加密切，越远的点关系越疏远，这个时候就需要加权平均。明显的离中心点越近的像素点权重越大。而正态分布显然是一种可取的权重分配方式，又由于图像是二维的，所以需要使用二维的高斯函数。所以高斯滤波的本质是利用高斯函数来生成高斯核（高斯卷积模板）来对图像进行卷积操作。
  理论上高斯分布在所有定义域上都有非负值，这就需要一个无限大的卷积核，但是实际上，仅需要取均值的三倍标准差（即3 σ 3sigma3σ）内的值，以外的部分去掉即可。
  高斯滤波最重要的就是找到高斯模板然后进行卷积，以3X3高斯模板为例，假设中心点的坐标为（0，0），根据二维高斯函数g ( x , y ) g(x,y)g(x,y)，还需要设定σ sigmaσ的值，假定σ = 0.8 sigma=0.8σ=0.8(这个值不宜过大，否则就会变成均值滤波)，可以根据坐标值来算出对应的高斯模板。
  假定中心点的坐标为（0，0），那么距离它最近的八个点的坐标如下：

  将坐标以及设定的σ sigmaσ带入二维高斯函数中，可以得到对应点的坐标的权重如下：

为了防止高斯滤波后的图像偏亮或者偏暗，我们还需要对图像进行归一化，这九个点的权重和等于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));