burg算法的matlab代码实现_导向滤波算法及其matlab代码实现

导向滤波同样是一种平滑滤波算法，其与最小二乘滤波和双边滤波相比，同样是一种具有边缘保持的功能的图形滤波算法，可以用于处理图形噪点较多的图像，而且此种滤波算法与最小二乘滤波和双边滤波相比，有其独特的特点，让我们来首先简述一下其滤波的基本原理：

滤波过程示意图

这是来源于Kaiming He（何恺明教授）其文章中的原理图，其中的p为输入图像，其中的q为输出图像，I为输入的引导图像，其中引导图像经过一个双边滤波核与输入图像进行结合，最终获得输出的图像。

导向滤波其基本的算法原理可以见下图：

算法原理

其中输出图像q有：

其中的输出图像q又有如下表达：

我们可以看到的是Wij为权值，这个权值限定了图像q的输出，而这个权值的来源是引导图像I，所以W是一个与I有关的数值，且其为一个滤波核，其中i与j为图像的像素坐标引导，由此表达式我们可以直观的看出，输出图像和输入图像p，存在的是一种线性关系。

在导向滤波中，核心的假设在于，在一个以k为核心的处理窗口中，输出图像的q是一种线性的模式，也即：

其中的ak和bk为窗口中的常系数，同时我们需要设定一个正方形的窗格，此窗格的半径为r。同时ak和bk的来源，我们都通过上文得知，其由输入的图像来进行限制，我们认为，输出的图像是输入图像减去噪声而形成的，故有：

其中的ni为图像中的噪声，而我们需要将此图像中的噪声降低到最低，所以需要找到一种比较好的方法来达到这一点：

我们可以看到上式为一个关于ak和bk的像素损失函数，同时我们又有一个调整参数，我们将调整参数设定为e，这个调整参数限定了ak。在何凯明教授之后的文中提到，此函数是一种脊回归模型，我们可以获取ak和bk其显式解为：

其中uk为引导图I中的窗口中的像素均值，而

为窗口中的像素方差，而下式的e

则为调整参数，而w为窗口的像素点的个数，而下式中的pk为：

那么当我们将此应用于具体的图像中时，对于一个N*N的处理窗格，一个像素点就会被包含N次，也即我们需要利用w，来进行取均值的处理：

同时我们注意到ak和bk，所以我们可以将上式改写为：

这样我们就建立了一个，从p，再依据I，最后到q的一种映射的过程，这也就是导向滤波整体算法原理，而当我们的导向图I和输入图片为同一张图片时，此种滤波算法就具有边缘保持的能力。

我们总结一下此种算法即：其输入端有4个输入：

输入图像p，引导图像I，处理窗口半径r，调整参数e。

而其输出为一张图像：

输出图像q。

其具体的算法实现过程：

我们利用matlab进行对原始的噪点较多的图形进行处理，首先我们是可以得知，此处我们输入图像和引导图像为同一张图像，此处我们设置三种窗口的规格

r0=2;r1=3;r2=4；

而对于调整参数，我们同样设置三种：

Eps0=0.01，eps1=0.04，eps2=0.09；

故对于此滤波器的调用方式为：

subplot(3,3,1);

gfout1_1=guidedfilter(gfI,gfp,r0,eps0);

imshow(gfout1_1);

title('r=2,eps=0.01');

lena女神

Guidedfilter.h:
function q = guidedfilter(I, p, r, eps)
 
%   - guidance image: I (should be a gray-scale/single channel image)
%   - filtering input image: p (should be a gray-scale/single channel image)
%   - local window radius: r
%   - regularization parameter: eps
 
[hei, wid] = size(I);
N = boxfilter(ones(hei, wid), r); 
 
mean_I = boxfilter(I, r) ./ N;
mean_p = boxfilter(p, r) ./ N;
mean_Ip = boxfilter(I.*p, r) ./ N;
% this is the covariance of (I, p) in each local patch.
cov_Ip = mean_Ip - mean_I .* mean_p; 
 
mean_II = boxfilter(I.*I, r) ./ N;
var_I = mean_II - mean_I .* mean_I;
 
a = cov_Ip ./ (var_I + eps); 
b = mean_p - a .* mean_I; 
 
mean_a = boxfilter(a, r) ./ N;
mean_b = boxfilter(b, r) ./ N;
 
q = mean_a .* I + mean_b; 
end

boxfilter.h:
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

算法总体评价：导向滤波作为三大边缘保持滤波器之一，其平滑效果和边缘保持能力值得肯定，但是其同时因为其算法原因会产生光晕（halo）此点何恺明教授也在其文章中详细说明，其与另外两种边缘保持滤波算法相比，其一大特点就在于其输入，其输入为一张输入图像和一张引导图像，所以这点是一个非常鲜明的特点，同时需注意的是，要是其存在边缘保持的滤波效果，需要使引导图像和输入图像为同一张图像，而且其处理窗口和损失参数的设置对其滤波效果的影响较大，此点可详见上文。

参考文献

____________________________________________________________________________________________

• Guided Image Filtering, by Kaiming He, Jian Sun, and Xiaoou Tang, in TPAMI 2013.

SANGHUSUN

2020.02.10