图像去雾算法

74 阅读 0 评论 49 点赞

我是靠谱客的博主文静电话，最近开发中收集的这篇文章主要介绍图像去雾算法，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

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```
```xml
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