1.概念

高斯分布

圖像濾波之高斯濾波介紹

圖像處理算法|高斯濾波

高斯濾波（Gaussian filter）包含很多種，包括低通、高通、帶通等，在圖像上說的高斯濾波通常是指的高斯模糊（Gaussian Blur），是一種高斯低通濾波。通常這個(gè)算法也可以用來模糊圖像，提供模糊濾鏡。也可以用來過濾自然界的高斯白噪聲。

高斯分布（正態(tài)分布）是一個(gè)常見的連續(xù)概率分布，正態(tài)分布的數(shù)學(xué)期望值或期望值μ muμ等于位置參數(shù)，決定了分布的位置，其方差σ 2 sigma^2σ2的開平方或者標(biāo)準(zhǔn)差σ sigmaσ等于尺度參數(shù)，決定了分布的幅度。正態(tài)分布的概率密度函數(shù)曲線呈鐘形，所以又被稱為鐘形曲線。我們常說的標(biāo)準(zhǔn)正態(tài)分布是位置參數(shù)μ = 0 , 方差 σ 2 = 1 mu=0,方差sigma^2=1μ=0,方差σ2=1的正態(tài)分布。
??若隨機(jī)變量X XX服從一個(gè)位置參數(shù)為μ 、方差為 σ 2 mu、方差為sigma^2μ、方差為σ2的正態(tài)分布，可以記為X N ( μ , σ 2 ) X~N(mu,sigma^2)XN(μ,σ2),其概率密度函數(shù)為：
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)
??高斯濾波器是一種根據(jù)高斯函數(shù)的形狀來選擇權(quán)值的線性平滑濾波器，對(duì)于抑制服從正態(tài)分布的噪聲非常有效，一維零均值高斯函數(shù)為：
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).
??其中，高斯分布參數(shù)（σ sigmaσ）決定了高斯函數(shù)的寬度，一維高斯函數(shù)的圖形如下圖：

??二維高斯高斯分布的函數(shù)為：
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.高斯濾波性質(zhì)

高斯函數(shù)具有五個(gè)重要的性質(zhì)，這些性質(zhì)使得它在早期圖像處理中特別有用．這些性質(zhì)表明，高斯平滑濾波器無論在空間域還是在頻率域都是十分有效的低通濾波器，且在實(shí)際圖像處理中得到了工程人員的有效使用．高斯函數(shù)具有五個(gè)十分重要的性質(zhì)，它們是：
??1.二維高斯函數(shù)具有旋轉(zhuǎn)對(duì)稱性，即濾波器在各個(gè)方向上的平滑程度是相同的．一般來說，一幅圖像的邊緣方向是事先不知道的，因此，在濾波前是無法確定一個(gè)方向上比另一方向上需要更多的平滑．旋轉(zhuǎn)對(duì)稱性意味著高斯平滑濾波器在后續(xù)邊緣檢測中不會(huì)偏向任一方向．
??2.高斯函數(shù)是單值函數(shù)．這表明，高斯濾波器用像素鄰域的加權(quán)均值來代替該點(diǎn)的像素值，而每一鄰域像素點(diǎn)權(quán)值是隨該點(diǎn)與中心點(diǎn)的距離單調(diào)增減的．這一性質(zhì)是很重要的，因?yàn)檫吘壥且环N圖像局部特征，如果平滑運(yùn)算對(duì)離算子中心很遠(yuǎn)的像素點(diǎn)仍然有很大作用，則平滑運(yùn)算會(huì)使圖像失真．
??3.高斯函數(shù)的傅立葉變換頻譜是單瓣的．正如下面所示，這一性質(zhì)是高斯函數(shù)付立葉變換等于高斯函數(shù)本身這一事實(shí)的直接推論．圖像常被不希望的高頻信號(hào)所污染(噪聲和細(xì)紋理)．而所希望的圖像特征（如邊緣），既含有低頻分量，又含有高頻分量．高斯函數(shù)付立葉變換的單瓣意味著平滑圖像不會(huì)被不需要的高頻信號(hào)所污染，同時(shí)保留了大部分所需信號(hào)．
??4.高斯濾波器寬度(決定著平滑程度)是由參數(shù)σ表征的，而且σ和平滑程度的關(guān)系是非常簡單的．σ越大，高斯濾波器的頻帶就越寬，平滑程度就越好．通過調(diào)節(jié)平滑程度參數(shù)σ，可在圖像特征過分模糊(過平滑)與平滑圖像中由于噪聲和細(xì)紋理所引起的過多的不希望突變量(欠平滑)之間取得折衷．
??5.由于高斯函數(shù)的可分離性，較大尺寸的高斯濾波器可以得以有效地實(shí)現(xiàn)．二維高斯函數(shù)卷積可以分兩步來進(jìn)行，首先將圖像與一維高斯函數(shù)進(jìn)行卷積，然后將卷積結(jié)果與方向垂直的相同一維高斯函數(shù)卷積．因此，二維高斯濾波的計(jì)算量隨濾波模板寬度成線性增長而不是成平方增長．

3.高斯濾波的原理與實(shí)現(xiàn)

3.1 高斯模板的生成

高斯濾波既能用來過濾高斯噪音，也可用來做高斯模糊。
??要模糊一張圖像，可以直接用均值濾波來做簡單的模糊，但是這樣做顯然不大合理，因?yàn)閳D像是連續(xù)的，離卷積核中心的點(diǎn)關(guān)系更加密切，越遠(yuǎn)的點(diǎn)關(guān)系越疏遠(yuǎn)，這個(gè)時(shí)候就需要加權(quán)平均。明顯的離中心點(diǎn)越近的像素點(diǎn)權(quán)重越大。而正態(tài)分布顯然是一種可取的權(quán)重分配方式，又由于圖像是二維的，所以需要使用二維的高斯函數(shù)。所以高斯濾波的本質(zhì)是利用高斯函數(shù)來生成高斯核（高斯卷積模板）來對(duì)圖像進(jìn)行卷積操作。
??理論上高斯分布在所有定義域上都有非負(fù)值，這就需要一個(gè)無限大的卷積核，但是實(shí)際上，僅需要取均值的三倍標(biāo)準(zhǔn)差（即3 σ 3sigma3σ）內(nèi)的值，以外的部分去掉即可。
??高斯濾波最重要的就是找到高斯模板然后進(jìn)行卷積，以3X3高斯模板為例，假設(shè)中心點(diǎn)的坐標(biāo)為（0，0），根據(jù)二維高斯函數(shù)g ( x , y ) g(x,y)g(x,y)，還需要設(shè)定σ sigmaσ的值，假定σ = 0.8 sigma=0.8σ=0.8(這個(gè)值不宜過大，否則就會(huì)變成均值濾波)，可以根據(jù)坐標(biāo)值來算出對(duì)應(yīng)的高斯模板。
??假定中心點(diǎn)的坐標(biāo)為（0，0），那么距離它最近的八個(gè)點(diǎn)的坐標(biāo)如下：

??將坐標(biāo)以及設(shè)定的σ sigmaσ?guī)攵S高斯函數(shù)中，可以得到對(duì)應(yīng)點(diǎn)的坐標(biāo)的權(quán)重如下：

為了防止高斯濾波后的圖像偏亮或者偏暗，我們還需要對(duì)圖像進(jìn)行歸一化，這九個(gè)點(diǎn)的權(quán)重和等于0.9125991，因此需要分別對(duì)這九個(gè)數(shù)除以0.9125991，最終得到的高斯模板為：

??有了高斯模板就可以對(duì)圖像進(jìn)行卷積了，但是在FPGA中對(duì)于小數(shù)的運(yùn)算不友好，于是我們將這個(gè)3X3模板擴(kuò)大了16倍，得到了近似的整數(shù)模板，卷積完成后再除以16來做定浮點(diǎn)數(shù)的近似計(jì)算。16倍后的模板高斯卷積模板為：

??3X3的高斯模板生成與定浮點(diǎn)數(shù)的MATLAB算法實(shí)現(xiàn)為：

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實(shí)現(xiàn)

生成卷積模板后的FPGA實(shí)現(xiàn)與仿真參考前面幾章，這里只貼高斯算法的FPGA實(shí)現(xiàn)部分：

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實(shí)現(xiàn)以及驗(yàn)證

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));