【图像去雾】基于暗通道图像去雾matlab源码

 何恺明的暗通道先验(dark channel prior)去雾算法是CV界去雾领域很有名的算法，关于该算法的论文"Single Image Haze Removal Using Dark Channel Prior"一举获得2009年CVPR最佳论文。作者统计了大量的无雾图像，发现一条规律：每一幅图像的每一个像素的RGB三个颜色通道中，总有一个通道的灰度值很低。基于这个几乎可以视作是定理的先验知识，作者提出暗通道先验的去雾算法。

       作者首先介绍去雾模型如下：

如果你不是CV界新手的话，应该对上式倒背如流，在此不再对上式中的各个参量作详细介绍。对于任意一幅输入图像，定义其暗通道的数学表达式为：

其中c表示rgb三通道中的某一通道。上式表示在一幅输入图像中，先取图像中每一个像素的三通道中的灰度值的最小值，得到一幅灰度图像，再在这幅灰度图像中，以每一个像素为中心取一定大小的矩形窗口，取矩形窗口中灰度值最小值代替中心像素灰度值(最小值滤波)，从而得到该雾天图像的暗通道图像。暗通道图像为灰度图像，通过大量统计并观察发现，暗通道图像的灰度值是很低的，所以将整幅暗通道图像中所有像素的灰度值近似为0，即：

        作者在文中假设大气光A为已知量，以便节省篇幅介绍传输函数的求解方法。在此介绍一种简单普遍的求解大气光的方法：对于任意一幅雾天图像，取其暗通道图像灰度值最大的0.1%的像素点对应于原雾天图像的像素位置的每个通道的灰度值的平均值，从而计算出每个通道的大气光值，即大气光值A是一个三元素向量，每一个元素对应于每一个颜色通道。

       对于成像模型，将其归一化，即两边同时除以每个通道的大气光值：

假设在图像中一定大小的矩形窗口内，传输函数的值为定值，对上式两边用最小化算子(minimum operators)作最小化运算：

 由于在矩形区域内为定值，故将其拿出运算符外部。由于场景辐射(scene radiance)是无雾图像，将暗通道先验应用于J，则有：

由于总是正值，则有：

 将上式代入到最小化运算的式子中，即可得到传输函数的估计值为：

为了防止去雾太过彻底，恢复出的景物不自然，应引入参数，重新定义传输函数为：

        对于求解得到的传输函数，应用滤波器对其进行细化，文章中介绍的方法是软抠图的方法，此方法过程复杂，速度缓慢，因此采用导向滤波对传输函数进行滤波。导向滤波的原理此处不再赘述，其伪代码为:

 上述伪代码中，I表示导向图像(guided image),p为输入图像(input image),q为输出图像(output image)，表示均值滤波，r为窗口半径。

```

%-------------------------------------- clc; clear; close all;

%% -----------图像去雾算法---------------- %% 加载图片 img = imread('foggybench.jpg'); figure;imshow(img);title('雾图'); %% 去雾函数 Deimg = anyuanse(img); figure;imshow(Deimg);title('去雾图'); %% 输出结果，分辨率300dpi并保存为tiff图片 imwrite(Deimg,'1.tiff','tiff','Resolution',300); function R = anyuanse(m_img) % 原始图像

windark=zeros(imgsize ,1);

for cc=1:imgsize windark(cc)=1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 windark=reshape(windark,h,w); %计算分块深度darkchannel for j=1+winsize:w-winsize for i=winsize+1:h-winsize mposmin = min(I(i,j,:)); for n=j-winsize:j+winsize
for m=i-winsize:i+winsize if(windark(m,n)>mposmin) windark(m,n)=mposmin; end end end

end

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 %图像透射率预处理，深度图反相 for cc=1:imgsize windark(cc)=1-windark(cc); end %%%%%%%%%%%%%%%%%图像软抠图开始%%%%%%%%%%%%%%%%%%%%% %选定精确dark value坐标 winb = zeros(img_size,1);

for ci=1:h for cj=1:w if(rem(ci-8,15)<1) %没有余数？ if(rem(cj-8,15)<1) winb(ci*w+cj)=windark(ci*w+cj); end end

end

end

%显示分块darkchannel nebsize = 9; winsize = 1; epsilon = 0.0000001; %指定矩阵形状 indsM=reshape([1:imgsize],h,w); %计算矩阵L tlen = imgsizeneb_size^2; row_inds=zeros(tlen ,1); col_inds=zeros(tlen,1); vals=zeros(tlen,1); len=0; for j=1+win_size:w-win_size for i=win_size+1:h-win_size if(rem(ci-8,15)<1) if(rem(cj-8,15)<1) continue; end end win_inds=indsM(i-win_size:i+win_size,j-win_size:j+win_size); win_inds=win_inds(:);%列显示 winI=I(i-win_size:i+win_size,j-win_size:j+win_size,:); winI=reshape(winI,neb_size,c); %三个通道被拉平成为一个二维矩阵 39 winmu=mean(winI,1)'; %求每一列的均值 如果第二个参数为2 则为求每一行的均值 //矩阵变向量 winvar=inv(winI'winI/neb_size-win_muwinmu' +epsilon/nebsize*eye(c)); %求方差

winI=winI-repmat(win_mu',neb_size,1);%求离差
  tvals=(1+winI*win_var*winI')/neb_size;% 求论文所指的矩阵L

  row_inds(1+len:neb_size^2+len)=reshape(repmat(win_inds,1,neb_size),...
                                         neb_size^2,1);
  col_inds(1+len:neb_size^2+len)=reshape(repmat(win_inds',neb_size,1),...
                                         neb_size^2,1);
  vals(1+len:neb_size^2+len)=tvals(:);
  len=len+neb_size^2;
end

end

vals=vals(1:len); rowinds=rowinds(1:len); colinds=colinds(1:len); %创建稀疏矩阵 A=sparse(rowinds,colinds,vals,imgsize,imgsize); %求行的总和 sumA为列向量 sumA=sum(A,2); %spdiags(sumA(:),0,imgsize,imgsize) 创建imgsize大小的稀疏矩阵其元素是sumA中的列元素放在由0指定的对角线位置上。 A=spdiags(sumA(:),0,imgsize,img_size)-A;

%创建稀疏矩阵 D=spdiags(winb(:),0,imgsize,imgsize); lambda=1; x=(A+lambda*D)(lambda*winb(:).*win_b(:)); %%%%%%%%%%%%%%%%%%%%%%%%%软图像抠图结束%%%%%%%%%%%%%%%55

%去掉0-1范围以外的数 alpha=max(min(reshape(x,h,w),1),0);%图像透射率

A=220/255; %大气光没有去计算 %去雾

for i=1:c for j=1:h for l=1:w dehaze(j,l,i)=(I(j,l,i)-A)/alpha(j,l)+A; end end end R = dehaze;function R = anyuanse(mimg) % 原始图像 I=double(mimg)/255;

% 获取图像大小 [h,w,c]=size(I); winsize = 7; imgsize=wh; dehaze=zeros(img_sizec,1); dehaze=reshape(dehaze,h,w,c);

windark=zeros(imgsize ,1);

for cc=1:imgsize windark(cc)=1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 windark=reshape(windark,h,w); %计算分块深度darkchannel for j=1+winsize:w-winsize for i=winsize+1:h-winsize mposmin = min(I(i,j,:)); for n=j-winsize:j+winsize
for m=i-winsize:i+winsize if(windark(m,n)>mposmin) windark(m,n)=mposmin; end end end

end

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 %图像透射率预处理，深度图反相 for cc=1:imgsize windark(cc)=1-windark(cc); end %%%%%%%%%%%%%%%%%图像软抠图开始%%%%%%%%%%%%%%%%%%%%% %选定精确dark value坐标 winb = zeros(img_size,1);

for ci=1:h for cj=1:w if(rem(ci-8,15)<1) %没有余数？ if(rem(cj-8,15)<1) winb(ci*w+cj)=windark(ci*w+cj); end end

end

end

%显示分块darkchannel nebsize = 9; winsize = 1; epsilon = 0.0000001; %指定矩阵形状 indsM=reshape([1:imgsize],h,w); %计算矩阵L tlen = imgsizeneb_size^2; row_inds=zeros(tlen ,1); col_inds=zeros(tlen,1); vals=zeros(tlen,1); len=0; for j=1+win_size:w-win_size for i=win_size+1:h-win_size if(rem(ci-8,15)<1) if(rem(cj-8,15)<1) continue; end end win_inds=indsM(i-win_size:i+win_size,j-win_size:j+win_size); win_inds=win_inds(:);%列显示 winI=I(i-win_size:i+win_size,j-win_size:j+win_size,:); winI=reshape(winI,neb_size,c); %三个通道被拉平成为一个二维矩阵 39 winmu=mean(winI,1)'; %求每一列的均值 如果第二个参数为2 则为求每一行的均值 //矩阵变向量 winvar=inv(winI'winI/neb_size-win_muwinmu' +epsilon/nebsize*eye(c)); %求方差

winI=winI-repmat(win_mu',neb_size,1);%求离差
  tvals=(1+winI*win_var*winI')/neb_size;% 求论文所指的矩阵L

  row_inds(1+len:neb_size^2+len)=reshape(repmat(win_inds,1,neb_size),...
                                         neb_size^2,1);
  col_inds(1+len:neb_size^2+len)=reshape(repmat(win_inds',neb_size,1),...
                                         neb_size^2,1);
  vals(1+len:neb_size^2+len)=tvals(:);
  len=len+neb_size^2;
end

end

R = dehaze; ```