%主程序
function q = gradient_guidedfilter(I, p, eps)

%   GUIDEDFILTER   O() 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)

%   - regularization parameter: eps  

r=;

[hei, wid] = size(I);

N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+)^ 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;  

%weight

epsilon=(0.001*(max(p(:))-min(p(:))))^;

r1=;  

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

mean_I1 = boxfilter(I, r1) ./ N1;

mean_II1 = boxfilter(I.*I, r1) ./ N1;

var_I1 = mean_II1 - mean_I1 .* mean_I1;  

chi_I=sqrt(abs(var_I1.*var_I));

weight=(chi_I+epsilon)/(mean(chi_I(:))+epsilon);       

gamma = (/(mean(chi_I(:))-min(chi_I(:))))*(chi_I-mean(chi_I(:)));

gamma =  - ./( + exp(gamma));  

%result

a = (cov_Ip + (eps./weight).*gamma) ./ (var_I + (eps./weight));

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()

[cpp] view plain copy

function imDst = boxfilter(imSrc, r)  

%   BOXFILTER   O() 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, [*r+, *r+], 'sliding', @sum);

%   - But much faster.  

[hei, wid] = size(imSrc);

imDst = zeros(size(imSrc));  

%cumulative sum over Y axis

imCum = cumsum(imSrc, );

%difference over Y axis

imDst(:r+, :) = imCum(+r:*r+, :);

imDst(r+:hei-r, :) = imCum(*r+:hei, :) - imCum(:hei-*r-, :);

imDst(hei-r+:hei, :) = repmat(imCum(hei, :), [r, ]) - imCum(hei-*r:hei-r-, :);  

%cumulative sum over X axis

imCum = cumsum(imDst, );

%difference over X axis

imDst(:, :r+) = imCum(:, +r:*r+);

imDst(:, r+:wid-r) = imCum(:, *r+:wid) - imCum(:, :wid-*r-);

imDst(:, wid-r+:wid) = repmat(imCum(:, wid), [, r]) - imCum(:, wid-*r:wid-r-);

end

%运行程序

clear

I = double(imread('D:\数字图像处理\研究方向\Filter Smooth\images\tulips.png')) / 255;
% if size(I,3)==3
% I=rgb2gray(I);
% end

p = I;
r=16;
eps = 0.8^2; % try eps=0.1^2, 0.2^2, 0.4^2

q_guide(:,:,1)=guidedfilter(I(:,:,1), p(:,:,1), r, eps);
q_guide(:,:,2)=guidedfilter(I(:,:,2), p(:,:,2), r, eps);
q_guide(:,:,3)=guidedfilter(I(:,:,3), p(:,:,3), r, eps);

q(:,:,1) = gradient_guidedfilter(I(:,:,1), p(:,:,1), eps);
q(:,:,2) = gradient_guidedfilter(I(:,:,2), p(:,:,2), eps);
q(:,:,3) = gradient_guidedfilter(I(:,:,3), p(:,:,3), eps);

figure;imshow([I,q_guide,q]);title('原图，引导滤波，改进引导滤波 eps=0.8^2');