Edge-Preserving Decomposition-Based Single Image Haze Removal TIP-2015

#标题:Edge-Preserving Decomposition-Based Single Image Haze Removal TIP-2015

基于边缘保留分解的单图像去雾

---

前言

     本论文主要简介单图像去雾算法，不按照先验知识。本文提出一种 新的模型：把简化暗通道分解为
     两层分别为基础层和细节层。

1. 求解大气光–通过四叉树细分的分层搜索方法估计
2. 求解传输图–根据加权引导图像滤波（WGIF），GIF作用是边缘保留的为局部滤波，加权作用消除边缘附近的伪影。

一、Modelling of Haze Images

****1. 大气雾模型

      Xc( p) = Zc( p)t( p) + Ac(1 − t( p)),     (5)
       t( p) = e−αd( p),                                 (6)
       0 ≤ t( p) ≤ 1                                        (8)

Z为去雾图，X为输入有雾图，右侧第一括号部分为放大因子，第二括号部分为细节层。

****2. 雾图像简化暗通道分解模型

Am = min{Ar , Ag, Ab},
Xm ( p) = min{Xr ( p), Xg( p), Xb( p)}

 Im=min(I,[],3);

Am − Xm ( p) = (Am − Zm( p))t( p). (13)

**模型：

Am − JdX ( p) = (Am − JdZ ( p))t( p). (16)

（16）函数为暗通道分解模型分为（）细节层+（）基础层

二、求解

1.估计全局大气光——四叉树细分的分层查找

1. 将输入雾图分为四个矩形块
2. 求每个矩形块中的（均值-方差）
3. 取四个块中最大的值
4. 继续1步骤直到矩形块小于设定阈值
5. 取像素最大的值 ，使􀀁(Xr ( p), Xg( p), Xb( p)) − (255, 255, 255)􀀁最小

代码如下（示例）：

function [Airlight]=AirlightEstimation(img_rgb)
%使用四叉树求大气光，均值减标准差
%estimation of the Global Atmosperic Light A
%img_rgb为输入的图像
% Airlight为求出的大气光A在R/G/B三通道的值是一个A[3]
[hei, wid, ~] = size(img_rgb);

figure, imshow(img_rgb); hold on; %hierarchical search display
line([1 wid], [hei/2 hei/2], 'Color', 'g');
line([wid/2 wid/2], [1 hei], 'Color', 'g');
yoffset = 1;
xoffset = 1;

while(hei * wid > 200)
   % divide the previous image into four equal subblocks
    img_ul = img_rgb(1: hei/2, 1: wid/2, :); 
    img_ur = img_rgb(1:hei/2, wid/2+1:wid, :);
    img_ll = img_rgb(hei/2+1:hei, 1:wid/2, :);
    img_lr = img_rgb(hei/2+1:hei, wid/2+1:wid, :);

    % calculate block mean - std
    dpscore(1) = mean2(img_ul(:,:,1)) - std2(img_ul(:,:,1));
    dpscore(2) = mean2(img_ul(:,:,2)) - std2(img_ul(:,:,2));
    dpscore(3) = mean2(img_ul(:,:,3)) - std2(img_ul(:,:,3));
    afscore(1) = sum(dpscore(:));
    
    dpscore(1) = mean2(img_ur(:,:,1)) - std2(img_ur(:,:,1));
    dpscore(2) = mean2(img_ur(:,:,2)) - std2(img_ur(:,:,2));
    dpscore(3) = mean2(img_ur(:,:,3)) - std2(img_ur(:,:,3));
    afscore(2) = sum(dpscore(:));

    dpscore(1) = mean2(img_ll(:,:,1)) - std2(img_ll(:,:,1));
    dpscore(2) = mean2(img_ll(:,:,2)) - std2(img_ll(:,:,2));
    dpscore(3) = mean2(img_ll(:,:,3)) - std2(img_ll(:,:,3));
    afscore(3) = sum(dpscore(:));
    
    dpscore(1) = mean2(img_lr(:,:,1)) - std2(img_lr(:,:,1));
    dpscore(2) = mean2(img_lr(:,:,2)) - std2(img_lr(:,:,2));
    dpscore(3) = mean2(img_lr(:,:,3)) - std2(img_lr(:,:,3));
    afscore(4) = sum(dpscore(:));
    
    % select the region with the highest score and 
    % divide it further into four smaller regions
    [maxscore, maxind] = max(afscore);%返回数组最大值对应的索引1*1或者1*2或者1*3或1*4
    clear img_rgb;           
    if(maxind == 1)
      img_rgb = img_ul;
      [hei, wid, ~] = size(img_ul);
    elseif(maxind == 2)
      img_rgb = img_ur;
      [hei, wid, ~] = size(img_ur);
      xoffset = xoffset + wid;
    elseif(maxind == 3)
      img_rgb = img_ll;
      [hei, wid, ~] = size(img_ll);
      yoffset = yoffset + hei;
    elseif(maxind == 4)
      img_rgb = img_lr;
      [hei, wid, ~] = size(img_lr);
      yoffset = yoffset + hei;
      xoffset = xoffset + wid;
    end
    line([xoffset, xoffset+wid], [yoffset+hei/2, yoffset+hei/2], 'Color', 'g');
    line([xoffset+wid/2, xoffset+wid/2], [yoffset, yoffset+hei], 'Color', 'g');
end
 rectangle('Position', [xoffset, yoffset, wid, hei], 'FaceColor', 'r', 'EdgeColor', 'g');

% choose the color vector which minimizes the distance as the atmospheric light
nDiff = abs(img_rgb - ones(hei, wid, 3));               
nDistance = nDiff(:,:,1) + nDiff(:,:,2) + nDiff(:,:,3);  %here abs-sum can be replaced by square-sum
[rowind, colind] = find(nDistance ==min(nDistance(:)));
row=rowind(1);
col=colind(1);

Airlight(1)=img_rgb(row,col,1);
Airlight(2)=img_rgb(row,col,2);
Airlight(3)=img_rgb(row,col,3);
end

2.估计传输图——利用WGIF

1. 求引导图G
   G( p) = Am − Xm( p). （17）

2. 求传输图
   
   按照WGIF的成本函数进行计算变量：
   

3. 结合（17）与（18）可以变为以下公式
   

传输图

代码如下（示例）：

%加权引导滤波求传输图
N1=m*n/255;                   %图像像素总值归一化处理
N = boxfilter(ones(m, n),60); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.
    
mean_ImJm = boxfilter(Im.*Jm, 60) ./ N;  %求均值
mean_Im = boxfilter(Im,60) ./ N;         %求均值
mean_Jm = boxfilter(Jm, 60) ./ N;        %求均值
% this is the covariance of (ImJm, Im*Jm) in each local patch.
cov_ImJm = mean_ImJm - mean_Im.* mean_Jm; %求相关

mean_ImIm = boxfilter(Im.*Im,60) ./ N; 
var_I = mean_ImIm - mean_Im .* mean_Im;   %求Im的方差

lamda=1;                                 %可变的变量
L = 1;                                   % dynamic range
epsilon=(0.001*L)^2;            
 
c=var_I+epsilon*ones(m,n);
%求加权引导滤波的权重
gamma0=N1*(var_I+epsilon*ones(m,n))/255*sum(c(:));
gamma = imgaussfilt(gamma0, 2);  %防止块效应，将权值进行高斯滤波处理，2为标准差

a1=cov_ImJm ./(var_I+lamda./gamma);  %求a
b1=mean_Jm -a1.*mean_Im;             %求b
 
mean_a1= boxfilter(a1, 60) ./ N;  %求均值
mean_b1 = boxfilter(b1, 60) ./ N;  %求均值
 
 fi=mean_a1.*Im+mean_b1;     %求fi也就是Am-t(p)

 t1=ones(m,n)-fi/Am;       %估算出传输图

3.估计去雾图——利用公式（9）

4.结果分析及改进

由于第3章中求出的结果中天空区域的噪音被放大，故需要把公式9中的放大因素进行限制，因为放大因素在天空区域很大。论文提到说在传输图上加一个非负的补偿项根据输入雾图的雾图级设置的，具体代码如下：

haze_level = 0.0625;   % 0(light), 0.03125(normal), 0.0625(heavy)雾图级别;
t_adjusted =t1 .^(1+haze_level);   % Eqn(28)，调整传输图

amplify_factor0 = max(0.1, 1./t_adjusted - 1);   %限制放大因子
amplify_factor1 = min(amplify_factor0 ,5);       %限制放大因子

结果图：

---

总结

总结：
以上就是今天要讲的内容，因为第一次发，可能做的不是很好，还有论文中一些代码部分参考其他人代码进行的改进，有些代码是自己想的（注意：由于浏览的时候没有留下记录故他人的贡献在这里未进行引）。