Guided Image Filtering

在图像滤波中，人们最希望的就是可以将图像中的噪声过滤掉的同时，能够让边缘尽可能的保持。噪声属于高频信号，而边缘其实也是一种高频信号，所以一般的滤波器，比如高斯模糊，均值模糊，都是一种低通滤波器，能够将高频信号做平滑处理，这些kernel在将噪声滤掉的同时，也会将边缘模糊掉，所以保边滤波器的设计一直是图像滤波的重点。

导向滤波是非常有名的一种保边滤波器，与传统的高斯滤波，均值滤波独立于图像的内容不同，导向滤波的 kernel 基于一个 guided image, 简单来说，就是给定一个 guided image I, 一个 input image p, 我们能够得到一个 output image q .

qi=∑jWijpj

i,j 表示像素的索引，Wij 是一个滤波系数，由 guided image I 决定，而与输入图像 p 无关。双边滤波器就是类似这样的一种滤波器，我们可以定义:

Wbfij(I)=1Kiexp(−|xi−xj|2σ2s)exp(−|Ii−Ij|2σ2r)

当 I 和 p 一样的时候，上式就是最原始的双边滤波器的表达式，xi,xj 表示像素的坐标，Ki 是一个归一化的系数，σs 控制像素空间的相似性， σr 控制像素值的相似性。

接下来，我们可以定义导向滤波器以及它的核函数，具体的定义形式如下:

qi=akIi+bk,∀i∈wk

这个表达式意味着，在一个局部区域 wk，导向分量 Ii 和输出 qi 是层线性关系的，
上式保证 了 ▽q=a▽I，意味着输出 q 和 导向图 I 具有同样的边界性质，
为了求解线性系数 ak,bk, 我们可以定义如下的能量函数：

E(ak,bk)=∑i∈wk((akIi+bk−pi)2+ϵa2k)

一般来说，ϵa2k 是一个正则项，以防止 ak 太大，最后利用线性回归，我们可以得到：

ak=1|w|∑i∈wkIipi−μkpk¯σ2k+ϵ

bk=pk¯−akμk

μk,σ2k 是导向图 I 在一个局部区域 wk 的均值和方差，|w| 是局部区域的像素个数，pk¯=1|w|∑i∈wkpi 是输入图 p 在局部区域的均值。

我们可以把这个线性模型应用到整张图像的所有局部区域，但是，由于一个像素点 i 可以同时属于很多不同的局部区域，而且每个局部区域计算出来的 qi 是不一样的，一个简单有效的方法，就是对含有像素点 i 的所有局部区域计算一个 线性模型，然后取平均值：

qi=1|w|∑k:i∈wk(akIi+bk)=ai¯Ii+bi¯

ai¯=1|w|∑k∈wiak

bi¯=1|w|∑k∈wibk

最后奉上代码：

function q = guidedfilter(I, p, r, eps)
%   GUIDEDFILTER   O(1) time implementation of guided filter.
%
%   - 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); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.

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; % this is the covariance of (I, p) in each local patch.

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

a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
b = mean_p - a .* mean_I; % Eqn. (6) in the paper;

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

q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
end

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

可以看到，guided image filtering 的代码实现还是很简单的，就是基于box filter。