2020华为杯E题--基于暗通道优先算法的能见度估计模型(附代码)
- 一、E题赛题
- 二、赛题分析与思路
- 三、基于暗通道优先算法的能见度估计模型
- 3.1 暗通道理论
- 3.2 大气物理模型
- 3.3 求解透射率 t
- 3.4 求消光系数
- 3.5求解能见度
- 四、结果分析
- 附
一、E题赛题
2020研究生数学建模赛题链接:https://download.csdn.net/download/qq_35759272/13028941
二、赛题分析与思路
建立不依赖能见度仪观测数据的能见度估计算法,即通过视频/图像中的信息获取场景的能见度。通过查阅文献,能见度可以用消光系数求出,而消光系数可以通过区域的透射率和拍摄距离求得。因此,首先通过暗通道理论获取暗通道图像,再通过大气散射物理模型求得图像得透射率,再通过深度图获取目标物与观测点的距离,再用该距离与透射率求出消光系数,再利用公式求出能见度。
分析思路步骤如图 :
三、基于暗通道优先算法的能见度估计模型
3.1 暗通道理论
在绝大多数非天空的局部区域里,某一些像素总会(至少一个颜色通道)具有很低的值。换言之,该区域光强度的最小值是个很小的数(趋于 0)。
基于上述结论,对于一幅图像 J,暗通道用公式描述为:
暗通道函数:
function Jdark = Idark( I )
% output: Jdark = min(min(r),min(g),min(b));
Wnd = 3;
Ir = I(:,:,1);
Ig = I(:,:,2);
Ib = I(:,:,3);
%% 图像拓展
[m,n,~] = size(I);
Irr = zeros(m+Wnd-1, n+Wnd-1);
Irr((Wnd-1)/2 : m+(Wnd-1)/2-1 , (Wnd-1)/2 : n+(Wnd-1)/2-1 ) = Ir;
Igg = zeros(m+Wnd-1, n+Wnd-1);
Igg((Wnd-1)/2 : m+(Wnd-1)/2-1 , (Wnd-1)/2 : n+(Wnd-1)/2-1 ) = Ig;
Ibb = zeros(m+Wnd-1, n+Wnd-1);
Ibb((Wnd-1)/2 : m+(Wnd-1)/2-1, (Wnd-1)/2 : n+(Wnd-1)/2-1 ) = Ib;
%% 暗通道
for i=1:1:m
for j=1:1:n
Rmin = min(min ( Irr(i:i+Wnd-1, j:j+Wnd-1) ));
Gmin = min(min ( Igg(i:i+Wnd-1, j:j+Wnd-1) ));
Bmin = min(min ( Ibb(i:i+Wnd-1, j:j+Wnd-1) ));
Jdark(i,j) = min(min(Rmin,Gmin),Bmin);
end
end
end
3.2 大气物理模型
如果从物理模型角度对有雾图像进行清晰化处理,需要用到大气散射模型。大气散射
物理模型包含两部分,第一部分称为直接衰减项也称直接传播,第二部分称为大气光照。
其原理图如图 5-17 所示:
通过公式表示如下:
公式中 J 代表的是景物反射光的强度,A 代表全局大气光照强度,t 用来描述光线通过介质
透射到成像设备过程中没有被散射的部分。
透射率 t 可表示为:
公式中β代表大气的散射系数,该式表明景物光线是随着景物深度 d 按指数衰减的。
3.3 求解透射率 t
3.4 求消光系数
根据透射率的定义公式,推导出的大气消光系数公式如下所示:
3.5求解能见度
将消光系数的结果带入能见度公式进行能见度计算,能见度的公式如下所示:
主函数:
clc;
clear all;
close all;
J = imread('4.jpg');
J = double(J);
J = J ./255 ;
figure(1); imshow(J);
%% 求暗通道图像 Jdark = min(min());
Jdark = Idark(J);
figure(2);imshow(Jdark,[]);
%%
% 注意:何凯明使用了soft matting方法对得到的粗透射率Jt进行细化
% 本代码采用梯度导向滤波实现
Jdark = gradient_guidedfilter(Jdark,Jdark, 0.04);
figure(3);imshow(Jdark,[]);
%% 大气物理模型 J = I*t + A*(1-t) 【直接衰减项】+【大气光照】
% 透射率 t与深度的关系 t=exp(-a*depth)
w = 0.95; %雾的保留系数
Jt = 1 - w*Jdark; %求解透射率
%% 求解全局大气光照
% 1.首先对输入的有雾图像I求解其暗通道图像Jdark。
% 2.选择Jdark总像素点个数千分之一(N/1000)个最亮的像素点,记录像素点(x,y)坐标
% 3.根据点的坐标分别在原图像J的三个通道(r,g,b)内找到这些像素点并加和得到(sum_r,sum_g,sum_b).
% 4.Ac=[Ar,Ag,Ab]. 其中Ar=sum_r/N; Ag=sum_g/N; Ab=sum_b/N.
[m,n,~] = size(J);
N = floor( m*n./1000 );
MaxPos = [0,0]; % 初始化
for i=1:1:N
MaxValue = max(max(Jdark));
[x,y] = find(Jdark==MaxValue);
Jdack(Jdark==MaxValue) = 0; %最大值置零,寻找下一次次大值
%检查长度
MaxPos = vertcat(MaxPos,[x,y]);
Cnt = length(MaxPos(1));
if Cnt > N
break;
end
end
MaxPosN = MaxPos(2:N+1,:);
Rsum = 0; Jr = J(:,:,1);
Gsum = 0; Jg = J(:,:,2);
Bsum = 0; Jb = J(:,:,3);
for j=1:1:N
Rsum = Rsum + Jr(MaxPosN(j,1),MaxPosN(j,2));
Gsum = Gsum + Jg(MaxPosN(j,1),MaxPosN(j,2));
Bsum = Bsum + Jb(MaxPosN(j,1),MaxPosN(j,2));
end
Ac = [Rsum/N, Gsum/N, Bsum/N];
基于梯度的指导滤波器
function q = gradient_guidedfilter(I, p, 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)
% - regularization parameter: eps
r=16;
[hei, wid] = size(I);
N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.
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;
%weight
epsilon=(0.001*(max(p(:))-min(p(:))))^2;
r1=1;
N1 = boxfilter(ones(hei, wid), r1); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.
mean_I1 = boxfilter(I, r1) ./ N1;
mean_II1 = boxfilter(I.*I, r1) ./ N1;
var_I1 = mean_II1 - mean_I1 .* mean_I1;
chi_I=sqrt(abs(var_I1.*var_I));
weight=(chi_I+epsilon)/(mean(chi_I(:))+epsilon);
gamma = (4/(mean(chi_I(:))-min(chi_I(:))))*(chi_I-mean(chi_I(:)));
gamma = 1 - 1./(1 + exp(gamma));
%result
a = (cov_Ip + (eps./weight).*gamma) ./ (var_I + (eps./weight));
b = mean_p - a .* mean_I;
mean_a = boxfilter(a, r) ./ N;
mean_b = boxfilter(b, r) ./ N;
q = mean_a .* I + mean_b;
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 X 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
四、结果分析
早上 6:30 到 7:39 高速公路的能见度整体呈递减的趋势,这与机场场景下的能见度变化类似,所以模型的估计结果符合实际。
附
第二题解题——基于AlexNet深度网络的能见度估计模型
解题文章链接:2020研究生数学建模E题–AlexNet深度网络解法(大雾能见度估计与预测)(含代码)
第三题解题——基于辅助车道线的大雾能见度估计模型
(1)对于机场视频可以采用——基于暗通道优先算法的能见度估计模型(即本文)
(2)对于高速公路视频截图数据采用——基于辅助车道线的大雾能见度估计与预测模型
解题文章链接:2020华为杯E题–基于辅助车道线的大雾能见度估计与预测(附代码)
第四题解题——基于灰色预测的大雾能见度预测模型
解题文章链接:2020华为杯E题——基于灰色预测的大雾能见度预测模型(附代码)
本博客参考文章:
【1】基于暗通道优先算法的去雾应用Matlab
【2】基于暗通道优先的单幅图像去雾算法(Matlab/C++)