前面发过中值、均值、高斯滤波的文章，这些只考虑了位置，并没有考虑相似度。那么双边滤波来了，既考虑了位置，有考虑了相似度，对边缘的保持比前几个好很多，当然实现上也是复杂很多。本文将从原理入手，采用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的滤波结果：
，而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(i<w+1 </w+1|| i>dim(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对双边滤波强度的影响程度如何，我们进一步分析。
其次是sigma_r，同时3*3窗口，sigma_d=3下对比，如下所示：
综上，感性的分析，对滤波强度的影响，首先是滤波窗口，其次是sigma_r，最后才是sigma_d。所以固定窗口下，sigma_r有较大的权重，也正是因此，相似度是双边权重的一个重要因素，当相似度高时，则权重收距离的影响更大；反之则收到相似度影响更大，因此也能有效地保护边缘，这就是双边滤波。
1.2.2. 定点双边滤波实现
，这里不难看出，相邻2像素的插值的绝对值，一定∈[0,255]，那么我们可以提前计算好
并定点化，那么就可以用查找表的方法来实现相似度权重的计算了。相关代码如下所示：
其中框框部分较为重要，详解如下：
1）第一个框：实现的是0-255下，
的查找表，同时为了避免出现浮点，将结果扩大了1024倍。
2）第二个框：根据查找表H，实现
的计算，其结果位宽为10bit；
3）第三个框：仍然将F2扩大1024倍后再运算，防止归一化出现浮点；
4）第四个框：将结果缩小1024倍，得到最终8bit的计算结果。
这样采用查找表的方式，我们不仅避免了指数的运算，同时也避免了浮点的运算，那么得到了定点化的
与
，后续的计算都不在话下了。这是FPGA实现时通用的定点化方法，也是本文后续的实现方式。
% 灰度图像双边滤波算法实现
% 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(i<w+1 </w+1|| i>dim(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】双边滤波结果');
最后，经常说双边滤波是一种磨皮算法，即去掉了斑纹又保持了边缘，那我们不妨找一个美女的头像测试一下，如下所示，为带噪声的神奇女侠（盖尔.加朵），在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实现
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的强度为例，高斯权重生成的数据如下所示：
% ----------------------------------------------------------------------
fp_gray = fopen('.\Similary_LUT.v','w');
fprintf(fp_gray,'//Sigma_r = %d\n', sigma_r);
fprintf(fp_gray,'module Similary_LUT\n');
fprintf(fp_gray,'(\n');
fprintf(fp_gray,' input\t\t[7:0]\tPre_Data,\n');
fprintf(fp_gray,' output\treg\t[9:0]\tPost_Data\n');
fprintf(fp_gray,');\n\n');
fprintf(fp_gray,'always@(*)\n');
fprintf(fp_gray,'begin\n');
fprintf(fp_gray,'\tcase(Pre_Data)\n');
% -----------------------------
% 计算相似度指数*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,'\t8''d%03d : Post_Data = 10''d%d; \n',i, Similary_ARRAY(i+1) );
end
fprintf(fp_gray,'\tendcase\n');
fprintf(fp_gray,'end\n');
fprintf(fp_gray,'\nendmodule\n');
fclose(fp_gray);
这里特别需要注意的是，为了最后的映射的数据保证10bit位宽（方便FPGA RTL设计），当Similary_ARRAY(i)=1时，其结果取1023，而非1024。
执行后生成Similar_LUT.v，完成查找表的映射，文件代码如下，该文件可直接用于sigma_r下exp指数的运算。
