一、前言

滤波作为最基础的图像处理手段之一，在图像处理领域占有重要位置，常被用于图像去噪、尺度分解等。从均值滤波到滚动导向滤波，滤波不断朝着精准分离图像中不同尺度信息的方向前进。
我在文中整理了双边滤波、导向滤波、滚动制导滤波三种在图像处理中常见且在论文中经常被使用的滤波方法。这三种滤波较之最基础的均值、高斯滤波有着更加优异的性能和可研究空间。直到现在，一些论文中提出的先进滤波算法仍是以它们为基础进行改进得到的。例如在图像多尺度分解领域中，以这三种滤波和它们的变体组成的图像分解模型几乎占据了基于滤波的图像多尺度分解的半壁江山。至今仍有学者在该方向上发表论文。
希望通过本文的讲述让大家一次性学会这三种最重要的滤波（当然，大家要先懂得滤波的概念并学习过如高斯滤波等一些最基础的滤波），话不多说，上正文！

二、双边滤波(Bilateral filter)

2.1 双边滤波的理论介绍及公式推导

双边滤波由C. Tomasi在1998年提出，是一种经典的非线性空间滤波方法。在滤波器稀疏的制定上，双边滤波同时考虑到了输出像素与邻域内其它像素的欧氏距离和取值的差异，即：同时考虑到了空间域和值域间的差别。如维纳滤波和高斯滤波等只考虑了空间域的滤波方法，在滤波后对边缘信息的保护效果不理想；如α-截尾均值滤波器等只考虑值域的滤波方法，在滤波后图像整体模糊，不能有效的保护细节信息。双边滤波器综合考量了空间域和值域对于滤波产生的影响，因而能达到保持边缘，降噪平滑的效果，是一种良好的边缘保持滤波器。
双边滤波通过基于高斯分布的加权平均方法实现。以图像中具体的像素值求解为里说明双边滤波的实现原理：图像在 ( i , j ) (i,j) (i,j)处经过双边滤波后的输出像素值 g g g依赖于邻域内像素值f的加权组合。
g ( i , j ) = ∑ k , l f ( k , l ) w ( i , j , k , l ) ∑ k , l w ( i , j , k , l ) g(i, j)=\frac{\sum_{k, l} f(k, l) w(i, j, k, l)}{\sum_{k, l} w(i, j, k, l)} g(i,j)=∑k,l​w(i,j,k,l)∑k,l​f(k,l)w(i,j,k,l)​
其中， k , l k,l k,l表示邻域像素的位置，权重系数 w ( i , j , k , l ) w(i,j,k,l) w(i,j,k,l)为空间域核d与值域核r的乘积。空间域核d是指基于高斯函数计算当前点与中心点的欧式距离。
d ( i , j , k , l ) = exp ⁡ ( − ( i − k ) 2 + ( j − l ) 2 2 σ d 2 ) d(i, j, k, l)=\exp \left(-\frac{(i-k)^2+(j-l)^2}{2 \sigma_d^2}\right) d(i,j,k,l)=exp(−2σd2​(i−k)2+(j−l)2​)
值域核 r r r是指基于高斯函数计算当前点与中心点像素值的差的绝对值。
r ( i , j , k , l ) = exp ⁡ ( − ∥ f ( i , j ) − f ( k , l ) ∥ 2 2 σ r 2 ) r(i, j, k, l)=\exp \left(-\frac{\|f(i, j)-f(k, l)\|^2}{2 \sigma_r^2}\right) r(i,j,k,l)=exp(−2σr2​∥f(i,j)−f(k,l)∥2​)
由空间域核与值域核的计算公式可得权重系数的计算公式为：
w ( i , j , k , l ) = exp ⁡ ( − ( i − k ) 2 + ( j − l ) 2 2 σ d 2 − ∥ f ( i , j ) − f ( k , l ) ∥ 2 2 σ r 2 ) w(i, j, k, l)=\exp \left(-\frac{(i-k)^2+(j-l)^2}{2 \sigma_d^2}-\frac{\|f(i, j)-f(k, l)\|^2}{2 \sigma_r^2}\right) w(i,j,k,l)=exp(−2σd2​(i−k)2+(j−l)2​−2σr2​∥f(i,j)−f(k,l)∥2​)

2.2 双边滤波的matlab程序实现

双边滤波实现函数：

%适用于单通道图像的双边滤波程序

function B = Bilater_Gray(A,w,sigma_d,sigma_r)
         %输出参数：
         % A为待滤波图像（double类型，取值在[0，1]）
         % w为滤波窗口的半径（e.g:3*3窗口的w值为1，w=3时的滤波效果较好）
         % sigma_d为定义域（空间域）核的方差，通常设置为3
         % sigma_r为值域核的方差，通常设置为0.1
         %输出参数：
         % B为滤波后的图像
% 预先计算高斯距离权重
[X,Y] = meshgrid(-w:w,-w:w);
%创建核距离矩阵，e.g.
%  [x,y]=meshgrid(-1:1,-1:1)
% 
% x =
% 
%     -1     0     1
%     -1     0     1
%     -1     0     1
% 
% 
% y =
% 
%     -1    -1    -1
%      0     0     0
%      1     1     1
%计算定义域核
G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));


%计算值域核H 并与定义域核G 乘积得到双边权重函数F
dim = size(A);
B = zeros(dim);
for i = 1:dim(1)
   for j = 1:dim(2)

         % 确定作用区域
         iMin = max(i-w,1);
         iMax = min(i+w,dim(1));
         jMin = max(j-w,1);
         jMax = min(j+w,dim(2));
         %定义当前核所作用的区域为(iMin:iMax,jMin:jMax)
         I = A(iMin:iMax,jMin:jMax);%提取该区域的源图像值赋给I

         %计算值域核H.
         H = exp(-(I-A(i,j)).^2/(2*sigma_r^2));

         % Calculate bilateral filter response.
         F = H.*G((iMin:iMax)-i+w+1,(jMin:jMax)-j+w+1);
         %在计算边缘部分的点的时候H的大小会变化，例如在计算第一行第一列的点时，
         %H的大小为4*4，因为7*7的其余部分都在图像外面（没有对应的值）。
         %因此适当的对G进行裁切使得G的大小始终能和H对上
         B(i,j) = sum(F(:).*I(:))/sum(F(:));

   end
end

主程序中调用：

%双边滤波主程序
clear all;
close all;
clc;
% 输入图像
fname   = 'Images\lena.jpg';   %改成你要操作的图像
X = double(rgb2gray(imread(fname)))/255;
% 开始滤波
Y = Bilater_Gray(X,3,3,0.1);

%获取细节层，即：被过滤的部分
Z = X - Y;

%结果显示
figure; imshow(X,[]);
figure; imshow(Y,[]);
figure; imshow(Z,[]);
%因为Z中有负值，所以最终强制以图片的形式显示的时候为灰色（正常现象）

结果展示：

图1 双边滤波效果展示(从左到右依次是：原始图像、滤波后的图像、被过滤的图像细节)

三、导向滤波(Guided Fliter)

3.1 导向滤波的理论介绍及公式推导

导向滤波是由何凯明等人在2010年提出来的一种非线性空间滤波算法。它在继承了双边滤波良好的边缘保留特性的同时，克服了双边滤波在主要边缘附近存在梯度变形的问题。
可以将导向滤波的过程看作是一个普通的线性平移变换滤波模型。滤波模型的输入量为：输入图像 p p p，引导图像 I I I，经过滤波模型可得到滤波后的图像 q q q。在这一过程中，对于 i i i位置像素点的滤波操作可以被表达成一个加权平均：
q i = ∑ j W i j ( I ) p j q_i=\sum_j W_{i j}(I) p_j qi​=j∑​Wij​(I)pj​
其中，i和j分别表示像素的下标。滤波核 W i j W_{i j} Wij​只与引导图像I相关。滤波器 W i j ( I ) W_{ij} (I) Wij​(I)与图像 p p p呈线性相关。引导图像 I I I是预先设定的，可以与 p p p保持一致。当 I = p I=p I=p时，导向滤波将退化为双边滤波。
假设输出图像 q q q是引导图像I的一个局部线性变换，即：在窗口 ω k ω_k ωk​上， q q q和 I I I线性相关：
q i = a k I i + b k , ∀ i ∈ ω k q_i=a_k I_i+b_k, \quad \forall i \in \omega_k qi​=ak​Ii​+bk​,∀i∈ωk​
其中 ω k ω_k ωk​是一个半径为r的方形局部化窗口，当 r r r为一个确定的常数值时， a k a_k ak​和 b k b_k bk​也将是一组确定的常系数。这保证了在一个局部区域里，如果引导图像 I I I存在边缘，输出图像 q q q也会保持引导图像的边缘特性。因为，对于边缘附近的像素点而言存在着 ∇ q = a ∇ I ∇q=a∇I ∇q=a∇I。对于这个局部线性变换模型而言，求解出系数 a a a和 b b b即可得到输出图像 q q q。为了求解系数 a a a和 b b b，需要约束输入图像 p p p，这里将输入图像被滤波器过滤的部分记为 n n n。
q i = p i − n i q_i=p_i-n_i qi​=pi​−ni​
为了在滤波的过程中极大程度的保留原始图像中的梯度信息，需要在保证线性模型的成立的前提下最小化 p p p和 q q q之间的差异。即：
argmin ⁡ ∑ i ∈ ω k ( q i − p i ) 2 argmin ⁡ ∑ i ∈ ω k ( a k I i + b k − p i ) 2 \begin{gathered} \operatorname{argmin} \sum_{i \in \omega_k}\left(q_i-p_i\right)^2 \\ \operatorname{argmin} \sum_{i \in \omega_k}\left(a_k I_i+b_k-p_i\right)^2 \end{gathered} argmini∈ωk​∑​(qi​−pi​)2argmini∈ωk​∑​(ak​Ii​+bk​−pi​)2​
在解决上式的最优化问题时，为防止 a k a_k ak​过大引入一个正则化参数 ϵ ϵ ϵ，可得最小化滤波窗口的损失函数为：
E ( a k , b k ) = ∑ i ∈ ω k ( ( a k I i + b k − p i ) 2 + ϵ a k 2 ) E\left(a_k, b_k\right)=\sum_{i \in \omega_k}\left(\left(a_k I_i+b_k-p_i\right)^2+\epsilon a_k^2\right) E(ak​,bk​)=i∈ωk​∑​((ak​Ii​+bk​−pi​)2+ϵak2​)
通过对参数求偏导的方法可以求解 a k a_k ak​和 b k b_k bk​。
a k = 1 ∣ ω ∣ ∑ i ϵ ω k I i p i − μ k p ˉ k σ k 2 + ϵ b k = p ˉ k − a k μ k \begin{gathered} a_k=\frac{\frac{1}{|\omega|} \sum_{i \epsilon \omega_k} I_i p_i-\mu_k \bar{p}_k}{\sigma_k^2+\epsilon} \\ b_k=\bar{p}_k-a_k \mu_k \end{gathered} ak​=σk2​+ϵ∣ω∣1​∑iϵωk​​Ii​pi​−μk​pˉ​k​​bk​=pˉ​k​−ak​μk​​
其中， μ k μ_k μk​和 σ k 2 σ_k^2 σk2​分别代表引导图像 I I I在窗口 ω k ω_k ωk​中的平均值和方差， p ˉ k \bar{p} _k pˉ​k​表示输入图像p在窗口内像素的均值， ∣ ω ∣ |ω| ∣ω∣代表窗口中的像素个数。

3.2 导向滤波matlab代码实现

导向滤波实现函数：

%导向滤波实现函数，适用于单通道图像
function q = Guided_filter(I, p, r, eps)
         %   时间复杂度为O(1).
         %输入参数：
         %   引导图像 I (double类型，取值为[0，1])
         %   待滤波图像 p (double类型，取值为[0，1])
         %   滤波窗口半径 r
         %   正则化参数: eps（可以取：0.2^2）
         %输出参数：
         %   滤波后的图像q

[hei, wid] = size(I);
N = boxfilter(ones(hei, wid), r); %boxfilter是一个求窗口内所有元素的和的程序(程序实现方式见下一代码块)
                                  %因此N代表了图像I在半径为r的窗口内的像素个数

mean_I = boxfilter(I, r) ./ N;
mean_p = boxfilter(p, r) ./ N;
mean_Ip = boxfilter(I.*p, r) ./ N;
cov_Ip = mean_Ip - mean_I .* mean_p; % 求局部区域内的协方差，即ak表达式的分子部分

mean_II = boxfilter(I.*I, r) ./ N;
var_I = mean_II - mean_I .* mean_I;

a = cov_Ip ./ (var_I + eps); % 求ak;
b = mean_p - a .* mean_I;    % 求bk;

mean_a = boxfilter(a, r) ./ N;
mean_b = boxfilter(b, r) ./ N;

q = mean_a .* I + mean_b;    %求滤波结果q;
end

boxfilter函数：

%这是一个通用的求半径为r的滤波窗口内所有元素的和的函数。
%也可以用matlab自带的colfilt来实现
function imDst = boxfilter(imSrc, r)

%   BOXFILTER   O(1) time box filtering using cumulative sum
%
%   - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));
%   - Running time independent of r; 
%   - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);
%   - But much faster.

[hei, wid] = size(imSrc);
imDst = zeros(size(imSrc));

%cumulative sum over Y axis
imCum = cumsum(imSrc, 1);
%difference over Y axis
imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);

%cumulative sum over X axis
imCum = cumsum(imDst, 2);
%difference over Y axis
imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
end

主程序调用：

%用于测试导向滤波的程序
close all;
fname_I   = 'Images\lena.jpg';   %改成你要操作的图像
fname_P   = 'Images\lena.jpg';   %改成你要操作的图像
%这里将引导图像设置成了原始图像，滤波效果会有所减弱
I = double(rgb2gray(imread(fname)))/255;
P = double(rgb2gray(imread(fname)))/255;

r = 2;
eps = 0.1^2;

q = Guided_filter(I, p, r, eps);

n = p - q;

%结果显示
figure; imshow(p,[]);  %原始图像
figure; imshow(q,[]);  %导向滤波后的图像
figure; imshow(n,[]);  %导向滤波过滤掉的图像的信息

图2 经过导向滤波过滤后的图像(图像展示顺序与图1的一致)

四、滚动导向滤波(RollingGuidedFilter)

4.1 滚动导向滤波的理论介绍及公式推导

滚动导向滤波是总结了双边滤波和导向滤波的经验后提出的一种新型的“两步滤波”，为了更好的说明它的设计出发点和它与前两种滤波的区别。我将在图像多尺度分解领域对其展开描述。
在图像滤波中（图像多尺度分解领域），双边滤波和导向滤波很好的发挥了边缘保持的特性，通过这两种滤波，图像的边缘特征已经能得到很大程度的保留。但是在实际的滤波过程中，单纯通过正向滤波过滤的过程仍然存在着错误识别和滤波过度的问题。简单来说就是将一些小尺度的边缘特征当作细节纹理过滤出去。
一个理想的滤波器是能够精确识别所有细节纹理然后将他们全部过滤，在这一过程不会存过滤过度或不足的情况。然而小尺度的细节纹理和边缘特征在像素构造上本身并没有太大区别，这就导致仅仅依靠正向操作的滤波器很难做到精准识别和选择性过滤。小尺度结构的保留问题亟需解决。
Qi Zhang等人在2014年提出了滚动导向滤波器，通过寻找尺度感知的方式极大程度的解决了这一问题。在这一问题的实际解决中，作者将这新型滤波器的设计追求定在过滤掉规定尺度下的所有细节纹理。
为了更好的理解滚动导向滤波的实现原理，先介绍细节纹理在图像中或在像素关系中处于的地位：一般来说图像中某一地方的像素与周围像素存在很大的像素值差异，那么存在这种差异的地方叫做震荡（像素值的高低），按照差异的大小分为高对比度/大结构/大震荡，低对比度/小结构/小震荡。震荡只是一种像素表现形式，按照图像的特征来看，震荡主要是图像纹理细节和边缘的反映。而纹理主要是小结构信息（因为长度短，通常所有的信息都能在一个滤波窗口中存下），边缘主要是大结构信息（比较长，一个滤波窗口不能完全存下，但在窗口中所有存在和体现）。
在确定了细节纹理和边缘信息的定义及二者在滤波窗口中的实际像素体现的区别后，现在的问题就变成了如何过滤小结构（细节纹理），但是要尽最大程度的不损伤大结构信息（边缘信息）。从单步处理来看，现有的滤波器存在的问题就是在滤波窗口内无法分辨出什么是细节，什么是边缘，也就是无法针对性的只过滤细节保留边缘。以滤波原理最简单的高斯滤波为例，在规定好滤波参数后，高斯滤波会将滤波窗口内所有符合过滤要求的振荡一并滤除。在这些振荡中有些是需要过滤的细节纹理，有些则是需要保留的边缘信息。
导向滤波要做的是只过滤所规定尺度的细节信息，不过滤边缘信息。传统的正向滤波（减法操作，从源图像中减去滤波器所认为的需要过滤的信息）很难实现，导向滤波的提出者们便另辟蹊径，提出了“现减后加”的理论模式。通过实际滤波结果，高斯滤波虽然将所有的振荡信息全部过滤，但对于大尺度（大于滤波上限）的（真正的）边缘信息高斯滤波只是过滤了一部分，也就是让它们变得模糊起来。那么可以先通过高斯滤波将细节纹理和特征信息全部过滤（减法），然后再用某种手段把被过滤的边缘信息加回来（加法），这样就实现了只过滤所规定尺度的细节信息，不过滤边缘信息的目标。
在加回边缘信息时就用到了高斯对于过滤不同尺度的边缘信息的不同表现这一现象，也就高斯滤波仅仅将大尺度边缘信息变得模糊，而不能彻底过滤。滚动导向滤波承袭导向滤波的设计思路，在滤波的过程中设置引导图像来完成自身的“加法”操作。经过高斯过滤后的引导图像中大尺度的边缘信息只是有所减弱，但还是存在的，用它去引导滤波核可以实现重点保护尺度信息（引导滤波核去保留边缘信息）。
总的来说只要是能在第一次（也是唯一的一次）高斯滤波中没有被彻底抹除的信息经过加边缘信息的操作最终都能被保存。而被保存（加回）的条件就是振荡的尺度要大于滤波尺度（这里讲的尺度可以理解为长度，纵深，是二维的占有面积的程度，和通常讲的特征强度不同）。传统的正向操作的滤波器是只要信息的振荡值大，不论尺度多大都能保存下来，但是RGF不同，RGF滤波的唯一条件就是尺度。因此RGF具有尺度感知特性，可以很好的探查小尺度边缘信息和细节纹理，不同于以往的任何滤波器。
下面给出滚动导向滤波的理论推导：
第一步：利用高斯滤波拆除小结构，实现细节信息和边缘结构的全部过滤。高斯滤波器为：
G ( p ) = 1 K p ∑ q ∈ N ( p ) exp ⁡ ( − ∥ p − q ∥ 2 2 σ s 2 ) I ( q ) K p = ∑ q ∈ N ( p ) exp ⁡ ( − ∥ p − q ∥ 2 2 σ s 2 ) \begin{aligned} &G(p)=\frac{1}{K_p} \sum_{q \in N(p)} \exp \left(-\frac{\|p-q\|^2}{2 \sigma_s^2}\right) I(q) \\ &K_p=\sum_{q \in N(p)} \exp \left(-\frac{\|p-q\|^2}{2 \sigma_s^2}\right) \end{aligned} ​G(p)=Kp​1​q∈N(p)∑​exp(−2σs2​∥p−q∥2​)I(q)Kp​=q∈N(p)∑​exp(−2σs2​∥p−q∥2​)​
其中输入图像为 I I I,输出图像为 G G G， p p p和 q q q表示对应的图像像素坐标， σ s σ_s σs​为标准差， N ( p ) N(p) N(p)表示以像素点 p p p为中心的滤波窗口。
第二步：迭代恢复边缘信息。以双边滤波核为滤波手段，通过迭代恢复边缘信息。迭代恢复的过程用公式表达为：
J t + 1 ( p ) = 1 K p ∑ q ∈ N ( p ) exp ⁡ ( − ∥ p − q ∥ 2 2 σ s 2 − ∥ J t ( p ) − J t ( q ) ∥ 2 2 σ r 2 ) I ( q ) K p = ∑ q ∈ N ( p ) exp ⁡ ( − ∥ p − q ∥ 2 2 σ s 2 − ∥ J t ( p ) − J t ( q ) ∥ 2 2 σ r 2 ) \begin{gathered} J^{t+1}(p)=\frac{1}{K_p} \sum_{q \in N(p)} \exp \left(-\frac{\|p-q\|^2}{2 \sigma_s^2}-\frac{\left\|J^t(p)-J^t(q)\right\|^2}{2 \sigma_r^2}\right) I(q) \\ K_p=\sum_{q \in N(p)} \exp \left(-\frac{\|p-q\|^2}{2 \sigma_s^2}-\frac{\left\|J^t(p)-J^t(q)\right\|^2}{2 \sigma_r^2}\right) \end{gathered} Jt+1(p)=Kp​1​q∈N(p)∑​exp(−2σs2​∥p−q∥2​−2σr2​∥Jt(p)−Jt(q)∥2​)I(q)Kp​=q∈N(p)∑​exp(−2σs2​∥p−q∥2​−2σr2​∥Jt(p)−Jt(q)∥2​)​
其中， J t + 1 J^{t+1} Jt+1表示为第