matlab暗通道影象去霧演算法實現
阿新 • • 發佈：2019-02-09
何凱明的經典去霧演算法，啥都不說啦，直接上程式碼，OK！
main.m
clear
clc
close all
kenlRatio = 0.01;
minAtomsLight = 240;
% image_name = 'test images\21.bmp';
image_name = 'C:\Users\Administrator\Desktop\Evernote\3.png';
img=imread(image_name);
figure,imshow(uint8(img)), title('src');
sz=size(img);
w=sz(2);
h=sz(1);
dc = zeros(h,w);
for guidedfilter.m
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;
endguidedfilter_color.m
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;
endboxfilter.m
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
maxfilt2.m
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
minfilt2.m
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);
Y = vanherk(Y,S(2),'min','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
vanherk.m
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程式執行結果如下：
見賢思齊，一點一點，不斷進步，總會柳暗花明，加油！