基于Matlab与FPGA的双边滤波算法实现

以下文章来源于疯狂的FPGA，作者CrazyFPGA

前面发过中值、均值、高斯滤波的文章，这些只考虑了位置，并没有考虑相似度。那么双边滤波来了，既考虑了位置，有考虑了相似度，对边缘的保持比前几个好很多，当然实现上也是复杂很多。本文将从原理入手，采用Matlab与FPGA设计实现双边滤波算法。

1.双边波算法的实现

滤波算法的基本思路，就是采用周边像素，加权平均计算一个新的像素，来缓减噪声对当前像素的影响。我们已经介绍了均值滤波，中值滤波，高斯滤波算法。但这几种算法都不够完美，还有继续提升的空间：

1）均值滤波：简单粗暴的将窗口内的像素累加后求均值，将噪声平均化，同时边缘纹理也被抹平了，有模糊的作用，作为入门学习用。

2）中值滤波：采用窗口内中值的方法，有效剔除了异常高亮或过暗的噪声，对椒盐噪声的去除效果比较好，但实际的图像会伴随着边缘纹理，由于只考虑中值，也会将图像细节去除，只是比均值滤波稍微好一点而已。

3）高斯滤波：采用欧氏距离，权重呈正态分布，相比均值/中值更优，因为均值/中值权重未考虑距离因素，而高斯噪声则考虑了噪声相对当前像素，呈高斯分布的特性，效果更佳。但高斯滤波也只是考虑了距离，未考虑边缘纹理，对细节的保护仍然不佳。

那么，既考虑噪声高斯分布特性，由将图像边缘纹理考虑进去的滤波算法应运而生。这里我们提出更进一步的滤波—双边滤波。

我们这里在前面高斯滤波的基础上，追加实现双边滤波，计划采用3*3的窗口，高斯滤波权重直接采用上一节的代码生成，重点讲解如何进行权重的计算，及Matlab&FPGA的实现。

1.1.高斯滤波算法理论

双边滤波是一种非线性滤波器，它可以达到降噪平滑，同时又保持边缘的效果。和其他滤波原理一样，双边滤波也是采用加权平均的方法，用周边像素亮度值的加权平均，来代表某个像素的强度，所用的加权平均也是基于高斯分布。

这里双边滤波的权重，不仅考虑了像素的欧式距离（如高斯滤波），还考虑了像素范围的辐射差异（比如像素与中心像素之间相似程度，也是高斯分布的），结合距离与相似度，最终计算得到最终的权重（距离与相似度的高斯分布）。

引用双边滤波算法相关图解，如下所示，其中为只考虑与当前像素空间距离的权重（space weight），而则为只考虑和当前像素相似度的权重（range weight），最后相乘，得到则为同时考虑了距离与相似度的权重，公式累加后最后归一化，得到最终的权重（space & rang）。

由于双边滤波同时考虑了空间距离和像素相似度的影响，因此尤其在具有边缘梯度的图像中，能够有不错的效果。即在平坦区域，空间距离占优势，在边缘区域，像素间相似度占优势，可以直观的用下面这个图来表示（注明出处）：

上述公式太抽象，我们进一步细化，并重新梳理如下。其中为归一化因子，高斯分布参数由于最终需要归一化，就直接省略了：

从公式可知，对于给定的窗口大小（比如3*3），以及确定的方差、，为常数，可以提前计算好高斯模板，具体的计算方式在高斯滤波中已经给出，以3*3窗口，为例，采用Data_Generate_nxn.m生成模板，如下所示（扩大了1024倍后）：

因此接下来需重点处理的是相似度权重，这也是我们Matlab与FPGA RTL实现的重点了。

1.2.双边滤波Matlab实现

1.2.1.浮点双边滤波实现

首先，以3*3窗口，为例，采用Data_Generate_nxn.m生成3*3高斯模板定点化后的结果，相关代码及生成的数据过程如下所示（扩大了1024倍后）：

最终得到的G3高斯模板为定点化后的结果，因此就不需要再计算，接下来只需要关心相似度权重即可，即的计算。针对灰度图像的权重计算，先计算当前像素的相似度分布权重，再计算最终的双边滤波，相关核心Matlab代码如下所示：

这里套用的仍然是我们滤波窗口处理的Matlab相关代码，以3*3为例，具体解释如下：

1）42行：获取以当前像素为中心的n*n的图像窗口A；

2）44行：计算以当前像素为中心的n*n窗口内的相似度权重；

3）45行：将高斯权重与相似度权重点乘，得到同时考虑距离与相似度的双边滤波权重；

4）46行：归一化双边滤波权重；

5）47行：当前窗口与双边滤波权重卷积，累加后得到最终的结果。

这里44、45行是最重要的地方，最终结果就是完成上面的计算结果。我们验证一下效果如何，如下所示，左图为原图，右图为3*3窗口，sigma_d = 3, sigma_r = 0.1的滤波结果：

最后我们封装function便于调用，再给出完整的代码。如下所示。这里的sigma_d就是，而sigma_r即是。  

% 灰度图像双边滤波算法实现

% I为输入的灰度图像

% n为滤波的窗口大小，为奇数

function B=bilateral_filter_gray(I,n,sigma_d, sigma_r)   

% ---------------------------------------------------

% 仅供function自测使用

% clear all;   close all;  clc;

% I = rgb2gray(imread('../../images/Scart.jpg'));    % 读取jpg图像

% n = 3; sigma_d = 3; sigma_r = 0.1; 

dim = size(I);   %读取图像高度、宽度

w=floor(n/2);   %窗口 [-w, w]

% ---------------------------------------------------

% 计算n*n高斯模板

G1=zeros(n,n);   %n*n高斯模板

for i=-w : w

for j=-w : w

G1(i+w+1, j+w+1) = exp(-(i^2 + j^2)/(2*sigma_d^2)) ;

end

end

% 归一化n*n高斯滤波模板

temp = sum(G1(:));

G2 = G1/temp;

% n*n高斯模板 *1024定点化

G3 = floor(G2*1024);

I = double(I);

% ---------------------------------------------------

% 计算n*n双边滤波模板+ 滤波结果

h = waitbar(0,'Speed of bilateral filter process...');  %创建进度条

B = zeros(dim);

for i=1 : dim(1)

for j=1 : dim(2)

if(idim(1)-w || jdim(2)-w)

B(i,j) = I(i,j);    %边缘像素取原值

else

A = I(i-w:i+w, j-w:j+w);

%              H =  exp( -(I(i,j)-A).^2/(2*sigma_r^2)  ) ;

H = exp( -((A-I(i,j))/255).^2/(2*sigma_r^2))  ;

F = G3.*H;

F2=F/sum(F(:));

B(i,j) = sum(F2(:) .*A(:));

end

end

waitbar(i/dim(1));

end

close(h);   % Close waitbar.

I = uint8(I);

B = uint8(B);

% ---------------------------------------------------

% 仅供function自测使用

% subplot(121);imshow(I);title('【1】原始图像');

% subplot(122);imshow(B);title('【2】双边滤波结果');

我们进一步进行测试，与迁移章中的高斯滤波，在相同窗口，及sigma_d下进行对比，相关代码及结果如下所示（Bilateral_Filter_Test1.m）：

clear all; 

close all;

clc;

% -------------------------------------------------------------------------

% Read PC image to Matlab

% IMG1= imread('../../images/Lenna.jpg');    % 读取jpg图像

IMG1= imread('../../images/Scart.jpg');    % 读取jpg图像

IMG1 = rgb2gray(IMG1);

h = size(IMG1,1);         % 读取图像高度

w = size(IMG1,2);         % 读取图像宽度

imshow(IMG1);title('【1】原图');

% -------------------------------------------------------------------------

figure;

% -------------------------------------------------------------------------

IMG2 = imfilter(IMG1, fspecial('gaussian',[3,3],1), 'replicate');

IMG3 = imfilter(IMG1, fspecial('gaussian',[3,3],2), 'replicate');

IMG4 = imfilter(IMG1, fspecial('gaussian',[5,5],3), 'replicate');

subplot(231);imshow(IMG2);title('【1】高斯滤波3*3, sigma = 1');

subplot(232);imshow(IMG3);title('【2】高斯滤波3*3, sigma = 2');

subplot(233);imshow(IMG4);title('【3】高斯滤波5*5, sigma = 3');

% -------------------------------------------------------------------------

IMG6 = bilateral_filter_gray(IMG1, 3, 1, 0.1); 

IMG7 = bilateral_filter_gray(IMG1, 3, 2, 0.3);   

IMG8 = bilateral_filter_gray(IMG1, 5, 3, 0.8);  

subplot(234);imshow(IMG6);title('【4】双边滤波3*3, sigma = [1, 0.1]');

subplot(235);imshow(IMG7);title('【5】双边滤波3*3, sigma = [2, 0.3]');

subplot(236);imshow(IMG8);title('【6】双边滤波5*5, sigma = [3, 0.8]');

如上图所示，第一行为高斯滤波，第二行为双边滤波，不难发现，双边滤波对边缘的保护成都更好。那么边缘滤波的3个参数：n、sigma_d、sigma_r对双边滤波强度的影响程度如何，我们进一步分析。

相关代码详见Bilateral_Filter_Test2.m首先是sigma_d，如下所示，可见滤波强度类似，sigma_d对双边权重的影响并不大。  

其次是sigma_r，同时3*3窗口，sigma_d=3下对比，如下所示：

最后进行滤波窗口系数n的测试，磨皮的效果出来了，可见对滤波强度的影响最大。  

综上，感性的分析，对滤波强度的影响，首先是滤波窗口，其次是sigma_r，最后才是sigma_d。所以固定窗口下，sigma_r有较大的权重，也正是因此，相似度是双边权重的一个重要因素，当相似度高时，则权重收距离的影响更大；反之则收到相似度影响更大，因此也能有效地保护边缘，这就是双边滤波。

1.2.2. 定点双边滤波实现

我们的目标是FPGA实现，但上述结Matlab计算中，主要是相似度权重计算中，仍然有浮点的参与。主要涉及公式，这里不难看出，相邻2像素的插值的绝对值，一定∈[0,255]，那么我们可以提前计算好并定点化，那么就可以用查找表的方法来实现相似度权重的计算了。相关代码如下所示：  

其中框框部分较为重要，详解如下：

1）第一个框：实现的是0-255下，的查找表，同时为了避免出现浮点，将结果扩大了1024倍。

2）第二个框：根据查找表H，实现的计算，其结果位宽为10bit；

3）第三个框：仍然将F2扩大1024倍后再运算，防止归一化出现浮点；

4）第四个框：将结果缩小1024倍，得到最终8bit的计算结果。

这样采用查找表的方式，我们不仅避免了指数的运算，同时也避免了浮点的运算，那么得到了定点化的与，后续的计算都不在话下了。这是FPGA实现时通用的定点化方法，也是本文后续的实现方式。

再次封装我们的函数，bilateral_filter_gray_INT(I,n,sigma_d, sigma_r)，增加“_INT”与bilateral_filter_gray区分，Matlab代码如下所示：  

% 灰度图像双边滤波算法实现

% I为输入的灰度图像

% n为滤波的窗口大小，为奇数

function B=bilateral_filter_gray_INT(I,n,sigma_d, sigma_r)   

% clear all;   close all;  clc;

% I = rgb2gray(imread('../../images/Scart.jpg'));    % 读取jpg图像

% n = 3; sigma_d = 3; sigma_r = 0.1; 

dim = size(I);   %读取图像高度、宽度

w=floor(n/2);   %窗口 [-w, w]

% ---------------------------------------------------

% 计算n*n高斯模板

G1=zeros(n,n);   %n*n高斯模板

for i=-w : w

for j=-w : w

G1(i+w+1, j+w+1) = exp(-(i^2 + j^2)/(2*sigma_d^2)) ;

end

end

% 归一化n*n高斯滤波模板

temp = sum(G1(:));

G2 = G1/temp;

% n*n高斯模板 *1024定点化

G3 = floor(G2*1024);

% ---------------------------------------------------

% 计算相似度指数*1024结果

% H = zeros(1, 256);

for i=0 : 255

H(i+1) = exp( -(i/255)^2/(2*sigma_r^2));

end

H = floor(H *1024);

I = double(I);

% ---------------------------------------------------

% 计算n*n双边滤波模板+ 滤波结果

h = waitbar(0,'Speed of bilateral filter process...');  %创建进度条

B = zeros(dim);

for i=1 : dim(1)

for j=1 : dim(2)

if(idim(1)-w || jdim(2)-w)

B(i,j) = I(i,j);    %边缘像素取原值

else

A = I(i-w:i+w, j-w:j+w);

%              H =  exp( -(I(i,j)-A).^2/(2*sigma_r^2)  ) ;

F1 = reshape(H(abs(A-I(i,j))+1), n, n);    %计算相似度权重(10bit)

F2 = F1 * G3;                                       %计算双边权重（20bit）

F3=F2*1024/sum(F2(:));                        %归一化双边滤波权重（扩大1024）

B(i,j) = sum(F3(:) .*A(:))/1024;                %卷积后得到最终累加的结果（缩小1024）

end

end

waitbar(i/dim(1));

end

close(h);   % Close waitbar.

I = uint8(I);

B = uint8(B);

% subplot(121);imshow(I);title('【1】原始图像');

% subplot(122);imshow(B);title('【2】双边滤波结果');

调用如上函数，与bilteral_Filter_gray函数进行对比，显示结果一致，验证了我们的改进结果，是符合预期的。  

最后，经常说双边滤波是一种磨皮算法，即去掉了斑纹又保持了边缘，那我们不妨找一个美女的头像测试一下，如下所示，为带噪声的神奇女侠（盖尔.加朵），在3*3窗口，sigma_d=3， sigma_r=0.1下的双边滤波效果。可见噪声被去除了，皮肤也变得光滑了，而且边缘纹理得以保留了。

这里处理的是彩色图像的双边滤波，需要RGB三通道，而我们的bilateral_filter_gray_INT()为灰度处理的函数，因此分3次分别处理了3个通道的数据，相关代码如下所示。当然也可以另行设计支持彩色图像的双边滤波函数，这个请读者自己努力，加油！

clear all; 

close all;

clc;

% -------------------------------------------------------------------------

% Read PC image to Matlab

IMG1= imread('../../images/girl.jpg');    % 读取jpg图像

% -------------------------------------------------------------------------

subplot(121);imshow(IMG1);title('【1】原图');

IMG2(:,:,1) = bilateral_filter_gray_INT(IMG1(:,:,1), 3, 3, 0.1); 

IMG2(:,:,2) = bilateral_filter_gray_INT(IMG1(:,:,2), 3, 3, 0.1); 

IMG2(:,:,3) = bilateral_filter_gray_INT(IMG1(:,:,3), 3, 3, 0.1); 

subplot(122);imshow(IMG2);title('【2】双边滤波3*3, sigma = [3, 0.1]');

1.3.双边滤波FPGA实现

我们整理上一节中高斯模板与相似度查找表的生成代码，文件名：Data_Generate_nxn_BF.m，Matlab代码如下所示。

clear all; 

close all;

clc;

sigma_d = 3; sigma_r=0.1;

n=3;

w=floor(n/2);

% ---------------------------------------------------

% 计算n*n高斯模板

G1=zeros(n,n);

for i=-w : w

for j=-w : w

G1(i+w+1, j+w+1) = exp(-(i^2 + j^2)/(2*sigma_d^2)) ;

end

end

% 归一化n*n高斯模板

temp = sum(sum(G1));

G2 = G1/temp;

% n*n高斯模板 *1024定点化

G3 = floor(G2*1024);

% ---------------------------------------------------

% 计算相似度指数*1024结果

H = zeros(1, 256);

for i=0 : 255

H(i+1) = exp( -(i/255)^2/(2*sigma_r^2));

end

H = floor(H *1024);

这里采用3*3的窗口，sigma_d=3, sigma_r=0.3的强度为例，高斯权重生成的数据如下所示：

而相似度查找表直接写入verilog文件，相关代码如下所示（详见Data_Generate_nxn_BF.m）：  

% ----------------------------------------------------------------------

fp_gray = fopen('.Similary_LUT.v','w');

fprintf(fp_gray,'//Sigma_r = %d ', sigma_r);

fprintf(fp_gray,'module Similary_LUT ');

fprintf(fp_gray,'( ');

fprintf(fp_gray,'   input [7:0] Pre_Data, ');

fprintf(fp_gray,'   output reg [9:0] Post_Data ');

fprintf(fp_gray,'); ');

fprintf(fp_gray,'always@(*) ');

fprintf(fp_gray,'begin ');

fprintf(fp_gray,' case(Pre_Data) ');

% -----------------------------

% 计算相似度指数*1024结果

Similary_ARRAY = zeros(1,256);

for i = 0 : 255

temp = exp( -(i/255)^2/(2*sigma_r^2));

if(temp == 1)   %i=0

Similary_ARRAY(i+1) = 1023;

else

Similary_ARRAY(i+1) = floor(temp*1023);

end

fprintf(fp_gray,' 8''d%03d : Post_Data = 10''d%d; ',i, Similary_ARRAY(i+1) );

end

fprintf(fp_gray,' endcase ');

fprintf(fp_gray,'end ');

fprintf(fp_gray,' endmodule ');  

fclose(fp_gray);

这里特别需要注意的是，为了最后的映射的数据保证10bit位宽（方便FPGA RTL设计），当Similary_ARRAY(i)=1时，其结果取1023，而非1024。

执行后生成Similar_LUT.v，完成查找表的映射，文件代码如下，该文件可直接用于sigma_r下exp指数的运算。