1.主程序
clear all;clc close all kenlRatio = .01;%窗口大小比例 minAtomsLight = 240; %原始论文中的A最终是取原始像素中的某一个点的像素，我实际上是取的符合条件的所有点的平均值作为A的值。 %如果是取一个点，则各通道的A值很有可能全部很接近255，这样的话会造成处理后的图像偏色和出现大量色斑。 %原文作者说这个算法对天空部分不需特别处理，我实际发现该算法对有天空的图像的效果一般都不好。天空会出现明显的过渡区域。 %作为解决方案，我增加了一个参数，最大全球大气光值minAtomsLight，当计算的值大于该值时，则就取该值。 image_name = 'defog.jpg';%图像名 img=imread(image_name);%读取图像 figure,imshow(uint(img)), title('原始图像');%显示原始图像 sz=size(img);%读取图像大小 w=sz(2);%图像宽度 h=sz(1);%图像高度 dc = zeros(h,w);%生成一个与原图像大小相同的零矩阵 for y=1:h for x=1:w dc(y,x) = min(img(y,x,:));%计算每个像素点RGB中最小的值 end end figure,imshow(uint8(dc)), title('Min(R,G,B)');%显示RGB取最小值后的图像 krnlsz = floor(max([3, w*kenlRatio, h*kenlRatio]));%计算窗口大小 dc2 = minfilt2(dc, [krnlsz,krnlsz]);%最小值滤波 dc2(h,w)=0; figure,imshow(uint8(dc2)), title('After Minfilter ');%显示最小值滤波后的图像 t = 255 - dc2;%计算透射率 figure,imshow(uint8(t)),title('t'); t_d=double(t)/255;%归一化 sum(sum(t_d))/(h*w) A = min([minAtomsLight, max(max(dc2))])%计算全球大气光的值 J = zeros(h,w,3);%生成一个与原图像大小相同的三维矩阵 img_d = double(img);%将整型转化为浮点型 J(:,:,1) = (img_d(:,:,1) - (1-t_d)*A)./t_d;%计算去雾后的R通道 J(:,:,2) = (img_d(:,:,2) - (1-t_d)*A)./t_d;%计算去雾后的G通道 J(:,:,3) = (img_d(:,:,3) - (1-t_d)*A)./t_d;%计算去雾后的B通道 figure,imshow(uint8(J)), title('J');%去雾后的图像 % figure,imshow(rgb2gray(uint8(abs(J-img_d)))), title('J-img_d'); % a = sum(sum(rgb2gray(uint8(abs(J-img_d))))) / (h*w) % return; %---------------------------------- r = krnlsz*4;%滤波半径 eps = 10^-6;%调整参数 % filtered = guidedfilter_color(double(img)/255, t_d, r, eps); filtered = guidedfilter(double(rgb2gray(img))/255, t_d, r, eps);%计算导向滤波图（灰度）求得精确的透射率 t_d = filtered;%将原来的透射率修改为新的透射率 figure,imshow(t_d,[]),title('filtered t'); J(:,:,1) = (img_d(:,:,1) - (1-t_d)*A)./t_d;%计算去雾后的R通道 J(:,:,2) = (img_d(:,:,2) - (1-t_d)*A)./t_d;%计算去雾后的G通道 J(:,:,3) = (img_d(:,:,3) - (1-t_d)*A)./t_d;%计算去雾后的B通道 % img_d(1,3,1) figure,imshow(uint8(J)), title('J_guild_filter');%显示图像 %---------------------------------- %imwrite(uint8(J), ['_', image_name])
2.vanherk
function Y = vanherk(X,N,TYPE,varargin) % VANHERK Fast max/min 1D filter % % Y = VANHERK(X,N,TYPE) performs the 1D max/min filtering of the row % vector X using a N-length filter. % The filtering type is defined by TYPE = 'max' or 'min'. This function % uses the van Herk algorithm for min/max filters that demands only 3 % min/max calculations per element, independently of the filter size. % % If X is a 2D matrix, each row will be filtered separately. % % Y = VANHERK(...,'col') performs the filtering on the columns of X. % % Y = VANHERK(...,'shape') returns the subset of the filtering specified % by 'shape' : % 'full' - Returns the full filtering result, % 'same' - (default) Returns the central filter area that is the % same size as X, % 'valid' - Returns only the area where no filter elements are outside % the image. % % X can be uint8 or double. If X is uint8 the processing is quite faster, so % dont't use X as double, unless it is really necessary. % % Initialization [direc, shape] = parse_inputs(varargin{:}); if strcmp(direc,'col') X = X'; end if strcmp(TYPE,'max') maxfilt = 1; elseif strcmp(TYPE,'min') maxfilt = 0; else error([ 'TYPE must be ' char(39) 'max' char(39) ' or ' char(39) 'min' char(39) '.']) end % Correcting X size fixsize = 0; addel = 0; if mod(size(X,2),N) ~= 0 fixsize = 1; addel = N-mod(size(X,2),N); if maxfilt f = [ X zeros(size(X,1), addel) ]; else f = [X repmat(X(:,end),1,addel)]; end else f = X; end lf = size(f,2); lx = size(X,2); clear X % Declaring aux. mat. g = f; h = g; % Filling g & h (aux. mat.) ig = 1:N:size(f,2); ih = ig + N - 1; g(:,ig) = f(:,ig); h(:,ih) = f(:,ih); if maxfilt for i = 2 : N igold = ig; ihold = ih; ig = ig + 1; ih = ih - 1; g(:,ig) = max(f(:,ig),g(:,igold)); h(:,ih) = max(f(:,ih),h(:,ihold)); end else for i = 2 : N igold = ig; ihold = ih; ig = ig + 1; ih = ih - 1; g(:,ig) = min(f(:,ig),g(:,igold)); h(:,ih) = min(f(:,ih),h(:,ihold)); end end clear f % Comparing g & h if strcmp(shape,'full') ig = [ N : 1 : lf ]; ih = [ 1 : 1 : lf-N+1 ]; if fixsize if maxfilt Y = [ g(:,1:N-1) max(g(:,ig), h(:,ih)) h(:,end-N+2:end-addel) ]; else Y = [ g(:,1:N-1) min(g(:,ig), h(:,ih)) h(:,end-N+2:end-addel) ]; end else if maxfilt Y = [ g(:,1:N-1) max(g(:,ig), h(:,ih)) h(:,end-N+2:end) ]; else Y = [ g(:,1:N-1) min(g(:,ig), h(:,ih)) h(:,end-N+2:end) ]; end end elseif strcmp(shape,'same') if fixsize if addel > (N-1)/2 disp('hoi') ig = [ N : 1 : lf - addel + floor((N-1)/2) ]; ih = [ 1 : 1 : lf-N+1 - addel + floor((N-1)/2)]; if maxfilt Y = [ g(:,1+ceil((N-1)/2):N-1) max(g(:,ig), h(:,ih)) ]; else Y = [ g(:,1+ceil((N-1)/2):N-1) min(g(:,ig), h(:,ih)) ]; end else ig = [ N : 1 : lf ]; ih = [ 1 : 1 : lf-N+1 ]; if maxfilt Y = [ g(:,1+ceil((N-1)/2):N-1) max(g(:,ig), h(:,ih)) h(:,lf-N+2:lf-N+1+floor((N-1)/2)-addel) ]; else Y = [ g(:,1+ceil((N-1)/2):N-1) min(g(:,ig), h(:,ih)) h(:,lf-N+2:lf-N+1+floor((N-1)/2)-addel) ]; end end else % not fixsize (addel=0, lf=lx) ig = [ N : 1 : lx ]; ih = [ 1 : 1 : lx-N+1 ]; if maxfilt Y = [ g(:,N-ceil((N-1)/2):N-1) max( g(:,ig), h(:,ih) ) h(:,lx-N+2:lx-N+1+floor((N-1)/2)) ]; else Y = [ g(:,N-ceil((N-1)/2):N-1) min( g(:,ig), h(:,ih) ) h(:,lx-N+2:lx-N+1+floor((N-1)/2)) ]; end end elseif strcmp(shape,'valid') ig = [ N : 1 : lx]; ih = [ 1 : 1: lx-N+1]; if maxfilt Y = [ max( g(:,ig), h(:,ih) ) ]; else Y = [ min( g(:,ig), h(:,ih) ) ]; end end if strcmp(direc,'col') Y = Y'; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [direc, shape] = parse_inputs(varargin) direc = 'lin'; shape = 'same'; flag = [0 0]; % [dir shape] for i = 1 : nargin t = varargin{i}; if strcmp(t,'col') & flag(1) == 0 direc = 'col'; flag(1) = 1; elseif strcmp(t,'full') & flag(2) == 0 shape = 'full'; flag(2) = 1; elseif strcmp(t,'same') & flag(2) == 0 shape = 'same'; flag(2) = 1; elseif strcmp(t,'valid') & flag(2) == 0 shape = 'valid'; flag(2) = 1; else error(['Too many / Unkown parameter : ' t ]) end end
3.导向滤波
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. % imwrite(uint8(N), 'N.jpg'); % figure,imshow(N,[]),title('N'); 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 q = guidedfilter_color(I, p, r, eps) % GUIDEDFILTER_COLOR O(1) time implementation of guided filter using a color image as the guidance. % % - guidance image: I (should be a color (RGB) image) % - filtering input image: p (should be a gray-scale/single channel image) % - local window radius: r % - regularization parameter: eps [hei, wid] = size(p); N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels. mean_I_r = boxfilter(I(:, :, 1), r) ./ N; mean_I_g = boxfilter(I(:, :, 2), r) ./ N; mean_I_b = boxfilter(I(:, :, 3), r) ./ N; mean_p = boxfilter(p, r) ./ N; mean_Ip_r = boxfilter(I(:, :, 1).*p, r) ./ N; mean_Ip_g = boxfilter(I(:, :, 2).*p, r) ./ N; mean_Ip_b = boxfilter(I(:, :, 3).*p, r) ./ N; % covariance of (I, p) in each local patch. cov_Ip_r = mean_Ip_r - mean_I_r .* mean_p; cov_Ip_g = mean_Ip_g - mean_I_g .* mean_p; cov_Ip_b = mean_Ip_b - mean_I_b .* mean_p; % variance of I in each local patch: the matrix Sigma in Eqn (14). % Note the variance in each local patch is a 3x3 symmetric matrix: % rr, rg, rb % Sigma = rg, gg, gb % rb, gb, bb var_I_rr = boxfilter(I(:, :, 1).*I(:, :, 1), r) ./ N - mean_I_r .* mean_I_r; var_I_rg = boxfilter(I(:, :, 1).*I(:, :, 2), r) ./ N - mean_I_r .* mean_I_g; var_I_rb = boxfilter(I(:, :, 1).*I(:, :, 3), r) ./ N - mean_I_r .* mean_I_b; var_I_gg = boxfilter(I(:, :, 2).*I(:, :, 2), r) ./ N - mean_I_g .* mean_I_g; var_I_gb = boxfilter(I(:, :, 2).*I(:, :, 3), r) ./ N - mean_I_g .* mean_I_b; var_I_bb = boxfilter(I(:, :, 3).*I(:, :, 3), r) ./ N - mean_I_b .* mean_I_b; a = zeros(hei, wid, 3); for y=1:hei for x=1:wid Sigma = [var_I_rr(y, x), var_I_rg(y, x), var_I_rb(y, x); var_I_rg(y, x), var_I_gg(y, x), var_I_gb(y, x); var_I_rb(y, x), var_I_gb(y, x), var_I_bb(y, x)]; Sigma = Sigma + eps * eye(3); cov_Ip = [cov_Ip_r(y, x), cov_Ip_g(y, x), cov_Ip_b(y, x)]; a(y, x, :) = cov_Ip * inv(Sigma + eps * eye(3)); % Eqn. (14) in the paper; end end b = mean_p - a(:, :, 1) .* mean_I_r - a(:, :, 2) .* mean_I_g - a(:, :, 3) .* mean_I_b; % Eqn. (15) in the paper; q = (boxfilter(a(:, :, 1), r).* I(:, :, 1)... + boxfilter(a(:, :, 2), r).* I(:, :, 2)... + boxfilter(a(:, :, 3), r).* I(:, :, 3)... + boxfilter(b, r)) ./ N; % Eqn. (16) in the paper; end
4.boxfilter
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
5.最小值滤波
function Y = minfilt2(X,varargin) % MINFILT2 Two-dimensional min filter % % Y = MINFILT2(X,[M N]) performs two-dimensional minimum % filtering on the image X using an M-by-N window. The result % Y contains the minimun value in the M-by-N neighborhood around % each pixel in the original image. % This function uses the van Herk algorithm for min filters. % % Y = MINFILT2(X,M) is the same as Y = MINFILT2(X,[M M]) % % Y = MINFILT2(X) uses a 3-by-3 neighborhood. % % Y = MINFILT2(..., 'shape') returns a subsection of the 2D % filtering specified by 'shape' : % 'full' - Returns the full filtering result, % 'same' - (default) Returns the central filter area that is the % same size as X, % 'valid' - Returns only the area where no filter elements are outside % the image. % % See also : MAXFILT2, VANHERK % % Initialization [S, shape] = parse_inputs(varargin{:}); % filtering Y = vanherk(X,S(1),'min',shape);%用van Herk算法对每一行进行滤波 Y = vanherk(Y,S(2),'min','col',shape);%用van Herk算法对每一列进行滤波 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [S, shape] = parse_inputs(varargin) shape = 'same'; flag = [0 0]; % size shape for i = 1 : nargin t = varargin{i}; if strcmp(t,'full') & flag(2) == 0 shape = 'full'; flag(2) = 1; elseif strcmp(t,'same') & flag(2) == 0 shape = 'same'; flag(2) = 1; elseif strcmp(t,'valid') & flag(2) == 0 shape = 'valid'; flag(2) = 1; elseif flag(1) == 0 S = t; flag(1) = 1; else error(['Too many / Unkown parameter : ' t ]) end end if flag(1) == 0 S = [3 3]; end if length(S) == 1; S(2) = S(1); end if length(S) ~= 2 error('Wrong window size parameter.') end
6.最大值滤波
function Y = maxfilt2(X,varargin) % MAXFILT2 Two-dimensional max filter % % Y = MAXFILT2(X,[M N]) performs two-dimensional maximum % filtering on the image X using an M-by-N window. The result % Y contains the maximun value in the M-by-N neighborhood around % each pixel in the original image. % This function uses the van Herk algorithm for max filters. % % Y = MAXFILT2(X,M) is the same as Y = MAXFILT2(X,[M M]) % % Y = MAXFILT2(X) uses a 3-by-3 neighborhood. % % Y = MAXFILT2(..., 'shape') returns a subsection of the 2D % filtering specified by 'shape' : % 'full' - Returns the full filtering result, % 'same' - (default) Returns the central filter area that is the % same size as X, % 'valid' - Returns only the area where no filter elements are outside % the image. % % See also : MINFILT2, VANHERK % % Initialization [S, shape] = parse_inputs(varargin{:}); % filtering Y = vanherk(X,S(1),'max',shape); Y = vanherk(Y,S(2),'max','col',shape); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [S, shape] = parse_inputs(varargin) shape = 'same'; flag = [0 0]; % size shape for i = 1 : nargin t = varargin{i}; if strcmp(t,'full') & flag(2) == 0 shape = 'full'; flag(2) = 1; elseif strcmp(t,'same') & flag(2) == 0 shape = 'same'; flag(2) = 1; elseif strcmp(t,'valid') & flag(2) == 0 shape = 'valid'; flag(2) = 1; elseif flag(1) == 0 S = t; flag(1) = 1; else error(['Too many / Unkown parameter : ' t ]) end end if flag(1) == 0 S = [3 3]; end if length(S) == 1; S(2) = S(1); end if length(S) ~= 2 error('Wrong window size parameter.') end