数字图像处理（二）：图像去雾

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

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

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

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

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

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