论文阅读：Guided Image Filtering, ECCV2010
Paper: Guided Image Filtering, ECCV2010
http://research.microsoft.com/en-us/um/people/kahe/eccv10/index.html
The Guided Image Filter is known as an edge-preserving smoothing filter. Compared to the traditional bilateral filter, it also retains the characteristics of edge preservation and non-iterative computation, while maintaining gradient information and ensuring that its computational complexity remains independent of the filter size.
经典的双边滤波算法如下图：
滤波权重系数由两个部分构成：反映空间位置关系的高斯函数与反映邻域颜色差异的高斯函数并存。该系统具有edge-preserving特性然而其算法计算复杂度为O(r²)
Guided filter原理如下，其中两张输入图像（p和I）可以相同，也可以不同：
将guided filter表达成类似于双边滤波的形式，如下图：
可以看出, 滤波权重系数W_{ij}完全取决于guided image. 对于每个包含像素i和j的局部窗口, 计算该窗口内所有像素的均值与方差. 然后通过前述公式计算权重系数. 邻域像素权重系数的具体分布情况如下图所示:
可以看出, 滤波权重系数W_{ij}完全取决于guided image. 对于每个包含像素i和j的局部窗口, 计算该窗口内所有像素的均值与方差. 然后通过前述公式计算权重系数. 邻域像素权重系数的具体分布情况如下图所示:
具体实现可参照作者的Matlab代码：
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;
endboxfilter可采用类似于快速积分图的方法，算法复杂度为O(1):
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