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

最新推荐文章于 2024-04-25 00:10:33 发布

原创最新推荐文章于 2024-04-25 00:10:33 发布 · 50 阅读

CC 4.0 BY-SA版权

文章标签：

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

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

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

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

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

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

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

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

【图像去雾】基于暗通道图像去雾matlab源码_图像处理_04

假设在图像中一定大小的矩形窗口【图像去雾】基于暗通道图像去雾matlab源码_图像处理_05 内，传输函数的值为定值，对上式两边用最小化算子(minimum operators)作最小化运算：

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

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

【图像去雾】基于暗通道图像去雾matlab源码_图像处理_10

由于【图像去雾】基于暗通道图像去雾matlab源码_matlab_11 总是正值，则有：

【图像去雾】基于暗通道图像去雾matlab源码_图像处理_12

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

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

为了防止去雾太过彻底，恢复出的景物不自然，应引入参数【图像去雾】基于暗通道图像去雾matlab源码_图像处理_14 ，重新定义传输函数为：

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

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

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

上述伪代码中，I表示导向图像(guided image),p为输入图像(input image),q为输出图像(output image)，【图像去雾】基于暗通道图像去雾matlab源码_图像处理_17 表示均值滤波，r为窗口半径。

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

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

win_dark=zeros(img_size ,1);

for cc=1:img_size
   win_dark(cc)=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
win_dark=reshape(win_dark,h,w);
%计算分块深度darkchannel
 for j=1+win_size:w-win_size
    for i=win_size+1:h-win_size
        m_pos_min = min(I(i,j,:));
        for n=j-win_size:j+win_size    
            for m=i-win_size:i+win_size
                if(win_dark(m,n)>m_pos_min)
                    win_dark(m,n)=m_pos_min;
                end
            end
        end
       
    end
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
 %图像透射率预处理，深度图反相
 for cc=1:img_size
   win_dark(cc)=1-win_dark(cc);
 end
 %%%%%%%%%%%%%%%%%图像软抠图开始%%%%%%%%%%%%%%%%%%%%%
 %选定精确dark value坐标
win_b = zeros(img_size,1);
 
for ci=1:h
    for cj=1:w
        if(rem(ci-8,15)<1)     %没有余数？
            if(rem(cj-8,15)<1)
                win_b(ci*w+cj)=win_dark(ci*w+cj);
            end
        end
       
    end
end
 
%显示分块darkchannel
neb_size = 9;
win_size = 1;
epsilon = 0.0000001;
%指定矩阵形状
indsM=reshape([1:img_size],h,w);
%计算矩阵L
tlen = img_size*neb_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); %三个通道被拉平成为一个二维矩阵 3*9
      win_mu=mean(winI,1)';  %求每一列的均值 如果第二个参数为2 则为求每一行的均值  //矩阵变向量
      win_var=inv(winI'*winI/neb_size-win_mu*win_mu' +epsilon/neb_size*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);
row_inds=row_inds(1:len);
col_inds=col_inds(1:len);
%创建稀疏矩阵
A=sparse(row_inds,col_inds,vals,img_size,img_size);
%求行的总和 sumA为列向量
sumA=sum(A,2);
%spdiags(sumA(:),0,img_size,img_size) 创建img_size大小的稀疏矩阵其元素是sumA中的列元素放在由0指定的对角线位置上。
A=spdiags(sumA(:),0,img_size,img_size)-A;


  %创建稀疏矩阵
  D=spdiags(win_b(:),0,img_size,img_size);
  lambda=1;
  x=(A+lambda*D)\(lambda*win_b(:).*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(m_img)
% 原始图像
I=double(m_img)/255;

% 获取图像大小
[h,w,c]=size(I);
win_size = 7;
img_size=w*h;
dehaze=zeros(img_size*c,1);
dehaze=reshape(dehaze,h,w,c);

win_dark=zeros(img_size ,1);

for cc=1:img_size
   win_dark(cc)=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
win_dark=reshape(win_dark,h,w);
%计算分块深度darkchannel
 for j=1+win_size:w-win_size
    for i=win_size+1:h-win_size
        m_pos_min = min(I(i,j,:));
        for n=j-win_size:j+win_size    
            for m=i-win_size:i+win_size
                if(win_dark(m,n)>m_pos_min)
                    win_dark(m,n)=m_pos_min;
                end
            end
        end
       
    end
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
 %图像透射率预处理，深度图反相
 for cc=1:img_size
   win_dark(cc)=1-win_dark(cc);
 end
 %%%%%%%%%%%%%%%%%图像软抠图开始%%%%%%%%%%%%%%%%%%%%%
 %选定精确dark value坐标
win_b = zeros(img_size,1);
 
for ci=1:h
    for cj=1:w
        if(rem(ci-8,15)<1)     %没有余数？
            if(rem(cj-8,15)<1)
                win_b(ci*w+cj)=win_dark(ci*w+cj);
            end
        end
       
    end
end
 
%显示分块darkchannel
neb_size = 9;
win_size = 1;
epsilon = 0.0000001;
%指定矩阵形状
indsM=reshape([1:img_size],h,w);
%计算矩阵L
tlen = img_size*neb_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); %三个通道被拉平成为一个二维矩阵 3*9
      win_mu=mean(winI,1)';  %求每一列的均值 如果第二个参数为2 则为求每一行的均值  //矩阵变向量
      win_var=inv(winI'*winI/neb_size-win_mu*win_mu' +epsilon/neb_size*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;