function q = guidedfilter(G, input, r, eps)
%   GUIDEDFILTER O(1) time implementation of guided filter.
%
%   - guidance image: G (should be a gray-scale/single channel image)
%   - filtering input image: input (should be a gray-scale/single channel image)
%   - local window radius: r
%   - regularization parameter: eps

% The essence of this function is the use of boxfilter. First, the matrix N
% means the weight of every pixel in the whole image. Then the average
% operation is done with the matrix N. This operation is done with the
% whole image, avoiding using a small window to scan the whole image.

[hei, wid] = size(G);
N = boxfilter(ones(hei, wid), r);  % the size of each local patch; N=(2r+1)^2 except for boundary pixels.

mean_G = boxfilter(G, r) ./ N;
mean_input = boxfilter(input, r) ./ N;
mean_Ginput = boxfilter(G.*input, r) ./ N;
cov_Ginput = mean_Ginput - mean_G .* mean_input; % this is the covariance of (G, input) in each local patch.

mean_GG = boxfilter(G.*G, r) ./ N;
var_G = mean_GG - mean_G .* mean_G;

a = cov_Ginput ./ (var_G + eps); % Eqn.(5) in the paper;
b = mean_input - a .* mean_G; % Eqn.(6) in the paper;

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

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