代码
主程序代码
m1.m
% MATLAB 图形用户界面 (GUI) 用于实现去雾
function varargout = m1(varargin)
% 设置 GUI 为单例模式，确保只能打开一个实例。
gui_Singleton = 1;
% 定义了 GUI 的状态，包括界面名字、是否单例模式等。
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @m1_OpeningFcn, ...
'gui_OutputFcn', @m1_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
% 判断输入参数是否是一个字符串，用于判断 GUI 的回调函数。
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
% 判断是否需要返回输出参数。
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% 在 GUI 可见之前执行的函数。将默认命令行输出作为 GUI 的输出。
function m1_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
guidata(hObject, handles);
% 设置 GUI 的输出函数。
function varargout = m1_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;
% 选择图片导入
function pushbutton1_Callback(hObject, eventdata, handles)
[imgfilename,imgpathname] = uigetfile({'*.jpg;*.png'},'输入图像');
if imgfilename
img = imread([imgpathname imgfilename]);
% 导入原图到第一个坐标系上
axes(handles.axes1);
imshow(img);
title('原图','FontWeight','Bold');
handles.imgfilename=imgfilename;
handles.imgpathname = imgpathname;
handles.img = img;
end
guidata(hObject,handles)
% 当按钮 pushbutton2 被点击时执行的回调函数。该函数使用全局直方图均衡化算法对图像进行去雾处理，
% 并在第二个坐标轴上显示去雾后的图像
function pushbutton2_Callback(hObject, eventdata, handles)
pnum=get(handles.popupmenu1,'value');
switch pnum
case 1
handles.img1=GlobalHistogram(handles.img);%执行全局直方图均衡化去雾
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
axes(handles.axes2);imshow(handles.img1,[]);%结果图显示在axes2
guidata(hObject,handles);
set(handles.edit1,'String','全局直方图算法');%将方法名显示在edit1中
case 2
handles.img1=Localhistogram1(handles.img);%执行局部直方图均衡化去雾
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
axes(handles.axes2);imshow(handles.img1,[]);%结果图显示在axes2
guidata(hObject,handles);
set(handles.edit1,'String','局部直方图算法');%将方法名显示在edit1中
case 3
handles.img1=Retinexsuanfa(handles.img);%执行Retinex单尺度算法去雾
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
axes(handles.axes2);imshow(handles.img1,[]);%结果图显示在axes2
guidata(hObject,handles);
set(handles.edit1,'String','Retinex单尺度算法');%将方法名显示在edit1中
case 4
handles.img1=Retinexsuanfa1(handles.img);%执行Retinex多尺度算法去雾
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
axes(handles.axes2);imshow(handles.img1,[]);%结果图显示在axes2
guidata(hObject,handles);
set(handles.edit1,'String','Retinex多尺度算法');%将方法名显示在edit1中
case 5
method = 'manual';
% Window Size设置为15
wsz = 15;
A = Airlight(handles.img, method, wsz);
% 计算边界约束
wsz = 3;
ts = Boundcon(handles.img, A, 30, 300, wsz);
% 改进传输评估
% 正则化参数，该参数越多，就越接近原始逐片传输
lambda = 2;
% 使用上下文信息
t = CalTransmission(handles.img, ts, lambda, 0.5);
% 去雾
handles.img1 = Dehazefun(handles.img, t, A, 0.85);
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
axes(handles.axes2);imshow(handles.img1,[]);%结果图显示在axes2
guidata(hObject,handles);
set(handles.edit1,'String','Dark_channel算法');%将方法名显示在edit1中
case 6
% estimate_airlight用于估计图像的大气光照（airlight）。它基于一种3*2D霍夫变换方法，其中每个点使用一组固定角度为给定位置进行投票。
A = reshape(estimate_airlight(im2double(handles.img).^(1)),1,1,3);
% 图像去雾
[handles.img1,trans_refined] = non_local_dehazing(handles.img, A, 1 ); % imgdata是输入图片handles是指针，A和gamma=1是它的参数
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
axes(handles.axes2);imshow(handles.img1,[]);%结果图显示在axes2
guidata(hObject,handles);
set(handles.edit1,'String','non_local算法');%将方法名显示在edit1中
otherwise
end
title('去雾后','FontWeight','Bold');
% 显示直方图
function pushbutton3_Callback(hObject, eventdata, handles)
img=rgb2gray(handles.img);
img1=rgb2gray(handles.img1);
axes(handles.axes3);imhist(img);title('原图灰度直方图','fontweight','bold');
axes(handles.axes4);imhist(img1);title('去雾后灰度直方图','fontweight','bold');
% 退出系统
function pushbutton4_Callback(hObject, eventdata, handles)
choice = questdlg('确定要退出系统?', ...
'退出', ...
'确定','取消','取消');
switch choice
case '确定'
close;
case '取消'
return;
end
% --- Executes on selection change in popupmenu1.
% 当用户改变下拉菜单的选择时触发，用来重置相关显示区域，清除之前的结果。
function popupmenu1_Callback(hObject, eventdata, handles)
pnum=get(handles.popupmenu1,'value');
if isequal(pnum==1, ~1)
axes(handles.axes2);
cla reset;
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
set(handles.edit1,'string','')
end
if isequal(pnum==2, ~2)
axes(handles.axes2);
cla reset;
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
set(handles.edit1,'string','')
end
if isequal(pnum==3, ~3)
axes(handles.axes2);
cla reset;
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
set(handles.edit1,'string','')
end
if isequal(pnum==4, ~4)
axes(handles.axes2);
cla reset;
axes(handles.axes3);
cla reset;
axes(handles.axes4);
cla reset;
set(handles.edit1,'string','')
end
% 仅当控件创建时调用，可以用来设置默认属性
function popupmenu1_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
% ispc 函数检查当前操作系统是否为 Windows（PC）。如果是在 Windows 上运行，则继续执行后续代码。
% isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 检查编辑框的背景颜色是否为MATLAB默认的用户界面控件背景颜色。
% set 函数将编辑框的背景颜色设置为白色 ('white')。这主要是为了确保在不同系统上的视觉一致性，因为MATLAB的默认背景颜色可能会随系统主题而变化。
function edit1_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
% 保存图片
% function pushbutton5_Callback(hObject, eventdata, handles)
% imwrite(handles.img1,'D:\课程作业\1.jpg')
function pushbutton5_Callback(hObject, eventdata, handles)
% 构造新的文件名，添加 'dehazing_' 前缀到原始文件名。
[~, name, ext] = fileparts(handles.imgfilename);
newFileName = ['dehazing_' name ext];
% 指定保存路径。这里假设你想要保存到与导入图片相同的目录。
savePath = fullfile(fileparts(handles.imgpathname), newFileName);
% 保存去雾后的图像。
imwrite(handles.img1, savePath);
% 可选：显示保存成功的消息框。
msg = sprintf('图片已保存为: %s', savePath);
uiwait(msgbox(msg, '保存成功', 'help'));
全局直方图均衡化
GlobalHistogram.m
function img1 = GlobalHistogram( img )
% 从输入图像img中分别各通道的像素值。
r=img(:,:,1);
g=img(:,:,2);
b=img(:,:,3);
% 对通道的像素值进行直方图均衡化处理
R=histeq(r);
G=histeq(g);
B=histeq(b);
% 将经过均衡化处理的通道重新组合成一个新的彩色图像
img1=cat(3,R,G,B);
end
局部直方图均衡化
Localhistogram1.m
function img1 = Localhistogram1( img )
% 分离通道
R =img(:, :, 1);
G =img(:, :, 2);
B =img(:, :, 3);
%对RGB分量进行局部直方图均衡化
IR=adapthisteq(R,'numtiles',[2 2],'clipLimit',0.02,'Distribution','rayleigh');
IG=adapthisteq(G,'numtiles',[2 2],'clipLimit',0.02,'Distribution','rayleigh');
IB=adapthisteq(B,'numtiles',[2 2],'clipLimit',0.02,'Distribution','rayleigh');
% 恢复到RGB
img1= cat(3, IR, IG, IB);
end
Retinex单尺度算法
Retinexsuanfa.m
function img = Retinexsuanfa( img )
R = img(:, :, 1);
%取对数将照射光分量和反射光分量分离
R0 = double(R);
Rlog = log(R0+1);
Rfft2 = fft2(R0);
%定义模板大小
[N1, M1] = size(R);
sigma1 =250;
F1= fspecial('gaussian', [N1,M1], sigma1);
Ffft1 = fft2(double(F1));
%高斯模板和原图像作卷积，低通滤波
DR0 = Rfft2.* Ffft1;
DR = ifft2(DR0);
%在对数域中，原始图像减去低通滤波后的图像，得到高频增强的图像
DRlog = log(DR +1);
Rr = Rlog - DRlog;
%取反对数
EXPRr = exp(Rr);
%对增强后的图像对比度拉伸增强
MIN = min(min(EXPRr));
MAX = max(max(EXPRr));
EXPRr1 = (EXPRr - MIN)/(MAX - MIN);
EXPRr2 = adapthisteq(EXPRr1);
G = img(:, :, 2);
G0 = double(G);
Glog = log(G0+1);
Gfft2 = fft2(G0);
%定义模板大小
[N2, M2] = size(G);
sigma2= 250;
F2= fspecial('gaussian', [N2,M2], sigma2);
Ffft2 = fft2(double(F2));
%高斯模板和原图像作卷积，低通滤波
DG0 = Gfft2.* Ffft2;
DG = ifft2(DG0);
DGlog = log(DG +1);
Gg = Glog - DGlog;
EXPGg = exp(Gg);
MIN = min(min(EXPGg));
MAX = max(max(EXPGg));
EXPGg1 = (EXPGg - MIN)/(MAX - MIN);
EXPGg2 = adapthisteq(EXPGg1);
B = img(:, :, 3);
B0 = double(B);
Blog = log(B0+1);
Bfft2 = fft2(B0);
%定义模板大小
[N3, M3] = size(B);
sigma3= 250;
F3 = fspecial('gaussian', [N3,M3], sigma3);
Ffft3 = fft2(double(F3));
%高斯模板和原图像作卷积，低通滤波
DB0 = Bfft2.* Ffft3;
DB = ifft2(DB0);
DBlog = log(DB+1);
Bb = Blog - DBlog;
EXPBb = exp(Bb);
MIN = min(min(EXPBb));
MAX = max(max(EXPBb));
EXPBb1 = (EXPBb - MIN)/(MAX - MIN);
EXPBb2 = adapthisteq(EXPBb1);
%RGB三通道合成一副图像
img= cat(3, EXPRr2, EXPGg2, EXPBb2);
end
Retinex多尺度算法
Retinexsuanfa1
function img = Retinexsuanfa1( img )
R = img(:, :, 1);
%取对数将照射光分量和反射光分量分离
R0 = double(R);
Rlog = log(R0+1);
Rfft2 = fft2(R0);
G = img(:, :, 2);
G0 = double(G);
Glog = log(G0+1);
Gfft2 = fft2(G0);
B = img(:, :, 3);
B0 = double(B);
Blog = log(B0+1);
Bfft2 = fft2(B0);
%色彩恢复因子C
a=2.5;
I1=imadd(R0,G0);
I2=imadd(I1,B0);
Ir=immultiply(R0,a);
C1=imdivide(Ir,I2);
C=log(C1+1);
%定义模板大小
[N1, M1] = size(R);
sigma1= 128;
F1= fspecial('gaussian', [N1,M1], sigma1);
Ffft1 = fft2(double(F1));
%高斯模板和原图像作卷积，低通滤波
DR01 = Rfft2.* Ffft1;
DR1 = ifft2(DR01);
%在对数域中，原始图像减去低通滤波后的图像，得到高频增强的图像
DRlog1 = log(DR1 +1);
Rr1 = Rlog - DRlog1;
sigma2= 256;
F2= fspecial('gaussian', [N1,M1], sigma2);
Ffft2 = fft2(double(F2));
%高斯模板和原图像作卷积，低通滤波
DR02 = Rfft2.* Ffft2;
DR2 = ifft2(DR02);
%在对数域中，原始图像减去低通滤波后的图像，得到高频增强的图像
DRlog2 = log(DR2 +1);
Rr2 = Rlog - DRlog2;
sigma3= 512;
F3 = fspecial('gaussian', [N1,M1], sigma3);
Ffft3 = fft2(double(F3));
%高斯模板和原图像作卷积，低通滤波
DR03 = Rfft2.* Ffft3;
DR3 = ifft2(DR03);
%在对数域中，原始图像减去低通滤波后的图像，得到高频增强的图像
DRlog3 = log(DR3 +1);
Rr3 = Rlog - DRlog3;
Rr=(1/3)*(Rr1+Rr2+Rr3);
Rr4=immultiply(C,Rr);
%取反对数
EXPRr = exp(Rr4);
%对增强后的图像对比度拉伸增强
MIN = min(min(EXPRr));
MAX = max(max(EXPRr));
EXPRr1 = (EXPRr - MIN)/(MAX - MIN);
EXPRr2 = adapthisteq(EXPRr1);
%定义模板大小
[N2, M2] = size(G);
sigma4= 128;
F4= fspecial('gaussian', [N2,M2], sigma4);
Ffft4 = fft2(double(F4));
%高斯模板和原图像作卷积，低通滤波
DG01= Gfft2.* Ffft4;
DG1 = ifft2(DG01);
DGlog1 = log(DG1 +1);
Gg1 = Glog - DGlog1;
sigma5= 256;
F5= fspecial('gaussian', [N2,M2], sigma5);
Ffft5 = fft2(double(F5));
%高斯模板和原图像作卷积，低通滤波
DG02 = Gfft2.* Ffft5;
DG2 = ifft2(DG02);
DGlog2 = log(DG2 +1);
Gg2 = Glog - DGlog2;
sigma6= 512;
F6 = fspecial('gaussian', [N2,M2], sigma6);
Ffft6 = fft2(double(F6));
%高斯模板和原图像作卷积，低通滤波
DG03 = Gfft2.* Ffft6;
DG3 = ifft2(DG03);
DGlog3 = log(DG3 +1);
Gg3 = Glog - DGlog3;
Gg=(1/3)*(Gg1+Gg2+Gg3);
Gg4=immultiply(C,Gg);
EXPGg = exp(Gg4);
MIN = min(min(EXPGg));
MAX = max(max(EXPGg));
EXPGg1 = (EXPGg - MIN)/(MAX - MIN);
EXPGg2 = adapthisteq(EXPGg1);
%定义模板大小
[N3, M3] = size(B);
sigma7= 128;
F7= fspecial('gaussian', [N3,M3], sigma7);
Ffft7 = fft2(double(F7));
%高斯模板和原图像作卷积，低通滤波
DB01 = Bfft2.* Ffft7;
DB1 = ifft2(DB01);
DBlog1 = log(DB1+1);
Bb1 = Blog - DBlog1;
sigma8= 256;
F8= fspecial('gaussian', [N3,M3], sigma8);
Ffft8 = fft2(double(F8));
%高斯模板和原图像作卷积，低通滤波
DB02 = Bfft2.* Ffft8;
DB2 = ifft2(DB02);
DBlog2 = log(DB2+1);
Bb2 = Blog - DBlog2;
sigma9= 512;
F9 = fspecial('gaussian', [N3,M3], sigma9);
Ffft9 = fft2(double(F9));
%高斯模板和原图像作卷积，低通滤波
DB03 = Bfft2.* Ffft9;
DB3 = ifft2(DB03);
DBlog3 = log(DB3+1);
Bb3 = Blog - DBlog3;
Bb=(1/3)*(Bb1+Bb2+Bb3);
Bb4=immultiply(C,Bb);
EXPBb = exp(Bb4);
MIN = min(min(EXPBb));
MAX = max(max(EXPBb));
EXPBb1 = (EXPBb - MIN)/(MAX - MIN);
EXPBb2 = adapthisteq(EXPBb1);
img= cat(3, EXPRr2, EXPGg2, EXPBb2);
end
Dark_channel算法
Airlight.m 估计图像的大气光
function A = Airlight(HazeImg, method, wsz)
% estimating the global airlight
%
if strcmp(method, 'manual')
h = figure, imshow(HazeImg, []);
title('manual airlight estimation: left click to pick a most hazy pixel. ')
[x, y] = ginput(1);
A = HazeImg(round(y), round(x), :);
A = double(A) -1;
A = min(A, 255);
close(h);
elseif strcmp(method, 'he')
A = airlight_he(HazeImg, wsz);
elseif strcmp(method, 'our')
A = airlight_our(HazeImg, wsz);
else
error('parameter error.');
end
%%
function A = airlight_our(HazeImg, wsz)
% estimating A channel by channel separately
for k = 1 : 3
minImg = ordfilt2(double(HazeImg(:, :, k)), 1, ones(wsz), 'symmetric');
A(k) = max(minImg(:));
end
%%
function A = airlight_he(HazeImg, wsz)
% estimating A using He's method
hsv = rgb2hsv(HazeImg);
GrayImg = hsv(:, :, 3);
[nRows, nCols, bt] = size(HazeImg);
% computing dark channel
DarkImg = min(double(HazeImg), [], 3);
DarkImg = ordfilt2(DarkImg, 1, ones(wsz), 'symmetric');
%
topDark = sort(DarkImg(:), 'descend');
idx = round(0.001 * length(topDark));
val = topDark(idx);
id_set = find(DarkImg >= val); % the top 0.1% brightest pixels in the dark channel
BrightPxls = GrayImg(id_set);
iBright = find(BrightPxls >= max(BrightPxls));
id = id_set(iBright); id = id(1);
row = mod(id, nRows);
col = floor(id / nRows) + 1;
% A is a vector
A = HazeImg(row, col, :);
A = double(A(:));
Boundcon.m计算边界约束
function [t_bdcon, t_b] = Boundcon(HazeImg, A, C0, C1, sz)
% patch-wise transmission from boundary constraint
if length(A) == 1
A = A * ones(3, 1);
end
if length(C0) == 1
C0 = C0 * ones(3, 1);
end
if length(C1) == 1
C1 = C1 * ones(3, 1);
end
HazeImg = double(HazeImg);
% pixel-wise boundary
t_r = max((A(1) - HazeImg(:, :, 1)) ./ (A(1) - C0(1)), (HazeImg(:, :, 1) - A(1)) ./ (C1(1) - A(1) ));
t_g = max((A(2) - HazeImg(:, :, 2)) ./ (A(2) - C0(2)), (HazeImg(:, :, 2) - A(2)) ./ (C1(2) - A(2) ));
t_b = max((A(3) - HazeImg(:, :, 3)) ./ (A(3) - C0(3)), (HazeImg(:, :, 3) - A(3)) ./ (C1(3) - A(3) ));
t_b = max(cat(3, t_r, t_g, t_b), [], 3);
t_b = min(t_b, 1);
% minimum filtering
se = strel('square', sz);
t_bdcon = imclose(t_b, se);
CalTransmission使用上下文信息
function t = CalTransmission(HazeImg, t, lambda, param)
% estimating the transmission function
[nRows, nCols] = size(t);
% differential filters bank
% Note: filter must have odd size to ensure the correct boundary
nsz = 3; NUM = nsz * nsz;
d{1} = [5, 5, 5; -3, 0, -3; -3, -3, -3];
d{2} = [-3, 5, 5; -3, 0, 5; -3, -3, -3];
d{3} = [-3, -3, 5; -3, 0, 5; -3, -3, 5];
d{4} = [-3, -3, -3; -3, 0, 5; -3, 5, 5];
d{5} = [5, 5, 5; -3, 0, -3; -3, -3, -3];
d{6} = [-3, -3, -3; 5, 0, -3; 5, 5, -3];
d{7} = [5, -3, -3; 5, 0, -3; 5, -3, -3];
d{8} = [5, 5, -3; 5, 0, -3; -3, -3, -3];
% normalizing filters
num_filters = length(d);
for k = 1 : num_filters
d{k} = d{k} / norm(d{k}(:));
end
% calculating weighting function
for k = 1 : num_filters
WFun{k} = CalWeightFun(HazeImg, d{k}, param);
end
% precomputing constant terms
Tf = fft2(t);
DS = 0;
for k = 1 : num_filters
D{k} = psf2otf(d{k}, [nRows, nCols]);
DS = DS + abs(D{k}).^2;
end
% cyclic looping for refining t
beta = 1; beta_rate = 2 * sqrt(2);
beta_max = 2^8;
Outiter = 0;
while beta < beta_max
% updating parameters
gamma = lambda / beta;
% show the results
Outiter = Outiter + 1;
fprintf('Outer iteration %d; beta %.3g\n', Outiter, beta);
figure(1000), imshow(t, []); title(num2str(Outiter)); pause(0.05);
% fixing t, solving u
DU = 0;
for k = 1 : num_filters
dt{k} = imfilter(t, d{k}, 'circular');
u{k} = max(abs(dt{k}) - WFun{k} / beta / num_filters, 0) .* sign(dt{k});
DU = DU + fft2(imfilter(u{k}, flipud(fliplr(d{k})), 'circular'));
end
% fixing u, solving t;
t = abs(ifft2((gamma * Tf + DU) ./ ( gamma + DS)));
% increasing beta
beta = beta * beta_rate;
end
close(1000);
%%
function WFun = CalWeightFun(HazeImg, D, param)
% parameters setting
sigma = param;
% calculating the weighting function
HazeImg = double(HazeImg) / 255;
% weighting function
method = 'circular';
d_r = imfilter(HazeImg(:, :, 1), D, method);
d_g = imfilter(HazeImg(:, :, 2), D, method);
d_b = imfilter(HazeImg(:, :, 3), D, method);
WFun = exp(-(d_r.^2 + d_g.^2 + d_b.^2) / sigma / 2);
Dehazefun.m去雾
function rImg = Dehazefun(HazeImg, t, A, delta)
% dehaze an image given t and A
t = max(abs(t), 0.0001).^delta;
% extropolation to dehaze
HazeImg = double(HazeImg);
if length(A) == 1
A = A * ones(3, 1);
end
R = (HazeImg(:, :, 1) - A(1)) ./ t + A(1); %R = max(R, 0); R = min(R, 255);
G = (HazeImg(:, :, 2) - A(2)) ./ t + A(2); %G = max(G, 0); G = min(G, 255);
B = (HazeImg(:, :, 3) - A(3)) ./ t + A(3); %B = max(B, 0); B = min(B, 255);
rImg = cat(3, R, G, B) ./ 255;
non_local算法
estimate_airlight.m估计图像的大气光照
function [ Aout ] = estimate_airlight( img, Amin, Amax, N, spacing, K, thres )
%% Verify input params, set defaults when necessary (same as published results)
if ~exist('thres','var') || isempty(thres), thres = 0.01 ; end;
if ~exist('spacing','var') || isempty(spacing), spacing = 0.02 ; end; %1/M in the paper
if ~exist('n_colors','var') || isempty(N), N = 1000 ; end; %number of colors clusters
if ~exist('K','var') || isempty(K), K = 40 ; end; %number of angles
% Define search range for the air-light. The search range is different for each
% color channel. These values were used in all of our experiments.
if ~exist('Amin','var') || isempty(Amin), Amin = [0,0.05,0.1]; end;
if ~exist('Amax','var') || isempty(Amax), Amax = 1; end;
% Air-light search range, accept a scalar if identical for all color channels
if isscalar(Amin), Amin = repmat(Amin,1,3); end
if isscalar(Amax), Amax = repmat(Amax,1,3); end
%% Convert input image to an indexed image
[img_ind, points] = rgb2ind(img, N);
[h,w,~] = size(img);
% Remove empty clusters
idx_in_use = unique(img_ind(:));
idx_to_remove = setdiff(0:(size(points,1)-1),idx_in_use);
points(idx_to_remove+1,:) = [];
img_ind_sequential = zeros(h,w);
for kk = 1:length(idx_in_use)
img_ind_sequential(img_ind==idx_in_use(kk)) = kk;
end
% Now the min value of img_ind_sequential is 1 rather then 0, and the indices
% correspond to points
% Count the occurences if each index - this is the clusters' weight
[points_weight,~] = histcounts(img_ind_sequential(:),size(points,1));
points_weight = points_weight./(h*w);
if ~ismatrix(points), points = reshape(points,[],3); end % verify dim
%% Define arrays of candidate air-light values and angles
angle_list = reshape(linspace(0, pi, K),[],1);
% Use angle_list(1:end-1) since angle_list(end)==pi, which is the same line
% in 2D as since angle_list(1)==0
directions_all = [sin(angle_list(1:end-1)) , cos(angle_list(1:end-1)) ];
% Air-light candidates in each color channel
ArangeR = Amin(1):spacing:Amax(1);
ArangeG = Amin(2):spacing:Amax(2);
ArangeB = Amin(3):spacing:Amax(3);
%% Estimate air-light in each pair of color channels
% Estimate RG
Aall = generate_Avals(ArangeR, ArangeG);
[~, AvoteRG] = vote_2D(points(:,1:2), points_weight, directions_all, Aall, thres );
% Estimate GB
Aall = generate_Avals(ArangeG, ArangeB);
[~, AvoteGB] = vote_2D(points(:,2:3), points_weight, directions_all, Aall, thres );
% Estimate RB
Aall = generate_Avals(ArangeR, ArangeB);
[~, AvoteRB] = vote_2D(points(:,[1,3]), points_weight, directions_all, Aall, thres);
%% Find most probable airlight from marginal probabilities (2D arrays)
% Normalize (otherwise the numbers are quite large)
max_val = max( [max(AvoteRB(:)) , max(AvoteRG(:)) , max(AvoteGB(:)) ]);
AvoteRG2 = AvoteRG./max_val;
AvoteGB2 = AvoteGB./max_val;
AvoteRB2 = AvoteRB./max_val;
% Generate 3D volumes from 3 different 2D arrays
A11 = repmat( reshape(AvoteRG2, length(ArangeG),length(ArangeR))', 1,1,length(ArangeB));
tmp = reshape(AvoteRB2, length(ArangeB),length(ArangeR))';
A22 = repmat(reshape(tmp, length(ArangeR),1,length(ArangeB)) , 1,length(ArangeG),1);
tmp2 = reshape(AvoteGB2, length(ArangeB),length(ArangeG))';
A33 = repmat(reshape(tmp2, 1, length(ArangeG),length(ArangeB)) , length(ArangeR),1,1);
AvoteAll = A11.*A22.*A33;
[~, idx] = max(AvoteAll(:));
[idx_r,idx_g,idx_b] = ind2sub([length(ArangeR),length(ArangeG),length(ArangeB)],idx);
Aout = [ArangeR(idx_r), ArangeG(idx_g), ArangeB(idx_b)];
end % function estimate_airlight_2D
%% Sub functions
function Aall = generate_Avals(Avals1, Avals2)
%Generate a list of air-light candidates of 2-channels, using two lists of
%values in a single channel each
%Aall's length is length(Avals1)*length(Avals2)
Avals1 = reshape(Avals1,[],1);
Avals2 = reshape(Avals2,[],1);
A1 = kron(Avals1, ones(length(Avals2),1));
A2 = kron(ones(length(Avals1),1), Avals2);
Aall = [A1, A2];
end % function generate_Avals
function [Aout, Avote2] = vote_2D(points, points_weight, directions_all, Aall, thres)
n_directions = size(directions_all,1);
accumulator_votes_idx = false(size(Aall,1), size(points,1), n_directions);
for i_point = 1:size(points,1)
for i_direction = 1:n_directions
% save time and ignore irelevant points from the get-go
idx_to_use = find( (Aall(:, 1) > points(i_point, 1)) & (Aall(:, 2) > points(i_point, 2)));
if isempty(idx_to_use), continue; end
% calculate distance between all A options and the line defined by
% i_point and i_direction. If the distance is smaller than a thres,
% increase the cell in accumulator
dist1 = sqrt(sum([Aall(idx_to_use, 1)-points(i_point, 1), Aall(idx_to_use, 2)-points(i_point, 2)].^2,2));
%dist1 = dist1 - min(dist1);
dist1 = dist1./sqrt(2) + 1;
dist = -points(i_point, 1)*directions_all(i_direction,2) + ...
points(i_point, 2)*directions_all(i_direction,1) + ...
Aall(idx_to_use, 1)*directions_all(i_direction,2) - ...
Aall(idx_to_use, 2)*directions_all(i_direction,1);
idx = abs(dist)<2*thres.*dist1;
if ~any(idx), continue; end
idx_full = idx_to_use(idx);
accumulator_votes_idx(idx_full, i_point,i_direction) = true;
end
end
% use only haze-lined that are supported by 2 points or more
accumulator_votes_idx2 = (sum(uint8(accumulator_votes_idx),2))>=2;
accumulator_votes_idx = bsxfun(@and, accumulator_votes_idx ,accumulator_votes_idx2);
accumulator_unique = zeros(size(Aall,1),1);
for iA = 1:size(Aall,1)
idx_to_use = find(Aall(iA, 1) > points(:, 1) & (Aall(iA, 2) > points(:, 2)));
points_dist = sqrt((Aall(iA,1) - points(idx_to_use,1)).^2+(Aall(iA,2) - points(idx_to_use,2)).^2);
points_weight_dist = points_weight(idx_to_use).*(5.*exp(-reshape(points_dist,1,[]))+1);
accumulator_unique(iA) = sum(points_weight_dist(any(accumulator_votes_idx(iA,idx_to_use,:),3)));
end
[~, Aestimate_idx] = max(accumulator_unique);
Aout = Aall(Aestimate_idx,:);
Avote2 = accumulator_unique;
end
wls_optimization优化传输图
function out = wls_optimization(in, data_weight, guidance, lambda)
small_num = 0.00001;
if ~exist('lambda','var') || isempty(lambda), lambda = 0.05; end
[h,w,~] = size(guidance);
k = h*w;
guidance = rgb2gray(guidance);
% Compute affinities between adjacent pixels based on gradients of guidance
dy = diff(guidance, 1, 1);
dy = -lambda./(sum(abs(dy).^2,3) + small_num);
dy = padarray(dy, [1 0], 'post');
dy = dy(:);
dx = diff(guidance, 1, 2);
dx = -lambda./(sum(abs(dx).^2,3) + small_num);
dx = padarray(dx, [0 1], 'post');
dx = dx(:);
% Construct a five-point spatially inhomogeneous Laplacian matrix
B = [dx, dy];
d = [-h,-1];
tmp = spdiags(B,d,k,k);
ea = dx;
we = padarray(dx, h, 'pre'); we = we(1:end-h);
so = dy;
no = padarray(dy, 1, 'pre'); no = no(1:end-1);
D = -(ea+we+so+no);
Asmoothness = tmp + tmp' + spdiags(D, 0, k, k);
% Normalize data weight
data_weight = data_weight - min(data_weight(:)) ;
data_weight = 1.*data_weight./(max(data_weight(:))+small_num);
% Make sure we have a boundary condition for the top line:
% It will be the minimum of the transmission in each column
% With reliability 0.8
reliability_mask = data_weight(1,:) < 0.6; % find missing boundary condition
in_row1 = min( in,[], 1);
data_weight(1,reliability_mask) = 0.8;
in(1,reliability_mask) = in_row1(reliability_mask);
Adata = spdiags(data_weight(:), 0, k, k);
A = Adata + Asmoothness;
b = Adata*in(:);
% Solve
% out = lsqnonneg(A,b);
out = A\b;
out = reshape(out, h, w);
adjust全局线性对比度拉伸
function adj = adjust(img,percen)
if ~exist('percen','var') || isempty(percen), percen=[0.01 0.99]; end;
% linear contrast stretch to [0,1], identical on all colors
minn=min(img(:));
img=img-minn;
img=img./max(img(:));
% limit the change magnitude so the WB would not be drastically changed
contrast_limit = stretchlim(img,percen);
val = 0.2;
contrast_limit(2,:) = max(contrast_limit(2,:), 0.2);
contrast_limit(2,:) = val*contrast_limit(2,:) + (1-val)*max(contrast_limit(2,:), mean(contrast_limit(2,:)));
contrast_limit(1,:) = val*contrast_limit(1,:) + (1-val)*min(contrast_limit(1,:), mean(contrast_limit(1,:)));
adj=imadjust(img,contrast_limit,[],1);
non_local_dehazing去雾
function [img_dehazed, transmission] = non_local_dehazing(img_hazy, air_light, gamma)
%% Validate input
[h,w,n_colors] = size(img_hazy);
if (n_colors ~= 3) % input verification
error(['Non-Local Dehazing reuires an RGB image, while input ',...
'has only ',num2str(n_colors),' dimensions']);
end
if ~exist('air_light','var') || isempty(air_light) || (numel(air_light)~=3)
error('Dehazing on sphere requires an RGB airlight');
end
if ~exist('gamma','var') || isempty(gamma), gamma = 1; end
img_hazy = im2double(img_hazy);
img_hazy_corrected = img_hazy.^gamma; % radiometric correction
%% Find Haze-lines
% Translate the coordinate system to be air_light-centric (Eq. (3))
dist_from_airlight = double(zeros(h,w,n_colors));
for color_idx=1:n_colors
dist_from_airlight(:,:,color_idx) = img_hazy_corrected(:,:,color_idx) - air_light(:,:,color_idx);
end
% Calculate radius (Eq. (5))
radius = sqrt( dist_from_airlight(:,:,1).^2 + dist_from_airlight(:,:,2).^2 +dist_from_airlight(:,:,3).^2 );
% Cluster the pixels to haze-lines
% Use a KD-tree impementation for fast clustering according to their angles
dist_unit_radius = reshape(dist_from_airlight,[h*w,n_colors]);
dist_norm = sqrt(sum(dist_unit_radius.^2,2));
dist_unit_radius = bsxfun(@rdivide, dist_unit_radius, dist_norm);
n_points = 1000;
% load pre-calculated uniform tesselation of the unit-sphere
fid = fopen(['TR',num2str(n_points),'.txt']);
points = cell2mat(textscan(fid,'%f %f %f')) ;
fclose(fid);
mdl = KDTreeSearcher(points);
ind = knnsearch(mdl, dist_unit_radius);
%% Estimating Initial Transmission
% Estimate radius as the maximal radius in each haze-line (Eq. (11))
K = accumarray(ind,radius(:),[n_points,1],@max);
radius_new = reshape( K(ind), h, w);
% Estimate transmission as radii ratio (Eq. (12))
transmission_estimation = radius./radius_new;
% Limit the transmission to the range [trans_min, 1] for numerical stability
trans_min = 0.1;
transmission_estimation = min(max(transmission_estimation, trans_min),1);
%% Regularization
% Apply lower bound from the image (Eqs. (13-14))
trans_lower_bound = 1 - min(bsxfun(@rdivide,img_hazy_corrected,reshape(air_light,1,1,3)) ,[],3);
transmission_estimation = max(transmission_estimation, trans_lower_bound);
% Solve optimization problem (Eq. (15))
% find bin counts for reliability - small bins (#pixels<50) do not comply with
% the model assumptions and should be disregarded
bin_count = accumarray(ind,1,[n_points,1]);
bin_count_map = reshape(bin_count(ind),h,w);
bin_eval_fun = @(x) min(1, x/50);
% Calculate std - this is the data-term weight of Eq. (15)
K_std = accumarray(ind,radius(:),[n_points,1],@std);
radius_std = reshape( K_std(ind), h, w);
radius_eval_fun = @(r) min(1, 3*max(0.001, r-0.1));
radius_reliability = radius_eval_fun(radius_std./max(radius_std(:)));
data_term_weight = bin_eval_fun(bin_count_map).*radius_reliability;
lambda = 0.1;
transmission = wls_optimization(transmission_estimation, data_term_weight, img_hazy, lambda);
%% Dehazing
% (Eq. (16))
img_dehazed = zeros(h,w,n_colors);
leave_haze = 1.06; % leave a bit of haze for a natural look (set to 1 to reduce all haze)
for color_idx = 1:3
img_dehazed(:,:,color_idx) = ( img_hazy_corrected(:,:,color_idx) - ...
(1-leave_haze.*transmission).*air_light(color_idx) )./ max(transmission,trans_min);
end
% Limit each pixel value to the range [0, 1] (avoid numerical problems)
img_dehazed(img_dehazed>1) = 1;
img_dehazed(img_dehazed<0) = 0;
img_dehazed = img_dehazed.^(1/gamma); % radiometric correction
% For display, we perform a global linear contrast stretch on the output,
% clipping 0.5% of the pixel values both in the shadows and in the highlights
adj_percent = [0.005, 0.995];
img_dehazed = adjust(img_dehazed,adj_percent);
img_dehazed = im2uint8(img_dehazed);
end % function non_local_dehazing
non_local算法还需要标识文件TR1000.txt，所有代码文件的链接如下：
通过网盘分享的文件：去雾
链接: https://pan.baidu.com/s/1YjNoFjbQ-TSEelS94d0uTA?pwd=3cbc 提取码: 3cbc
–来自百度网盘超级会员v4的分享
效果展示
主程序的运行效果
全局直方图算法去雾
局部直方图算法去雾
单尺度算法去雾
多尺度算法去雾
Dark_channel算法去雾
non_local算法去雾
扫一扫
136
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
