Veiling Light Estimation for Underwater Dehazing

注意：这篇文章上次更新于1394天前，文章内容可能已经过时。

求大气散射模型中大气光的一种方法，这里记为V

成像模型：

$I(x) = S(x)t(x) + (1 - t(x))V$

点击查看原文

方法简介

求 V 的过程是一种垂直分层的搜索技术，主要有以下几个步骤：

分层操作。

分层并不是真的图像分层，而是用不同尺度的分割方式对图像进行分块计算特征。第 i 层把图像分为 $4^i$ 块。
计算特征。

对每一层的每一块分别计算下面三个特征。
- 特征 1：
  
  $\mu_{i, j}^{d^{G B}}=\frac{1}{\left|\Omega_{i, j}\right|} \sum_{x \in \Omega_{i, j}} C_{i, j}^{d^{G B}}(x)$
  
  $f_{1_{i, j}}=1-\mu_{i, j}^{d^{G B}} \in[0,1]$
  
  其中 $C_{i,j}^{d^{G B}}(x)$ 表示仅使用绿色和蓝色通道计算的暗通道图。 $\left|\Omega_{i, j}\right|$ 表示第 j 层第 i 个区域的像素点个数。
- 特征 2 ：
  
  $\sigma_{i, j}^{c}=\sqrt{\frac{1}{\left|\Omega_{i, j}\right|-1} \sum_{x \in \Omega_{i, j}}\left(\mu_{i, j}^{c}-I_{i, j}^{c}(x)\right)^{2}}$
  
  $f_{2_{i, j}}=\frac{\sigma_{i, j}^{R}+\sigma_{i, j}^{G}+\sigma_{i, j}^{B}}{3} \in[0,1]$
  
  其中 $\mu_{i, j}^{c}$ 代表第 j 层第 i 个区域的像素均值，I 为原图。
- 特征 3 ：
  
  $f_{3_{i, j}}=\frac{1}{\left|\Omega_{i, j}\right|} \sum_{x \in \Omega_{i, j}}\left(G_{i, j}^{\prime}(x)\right) \quad \in[0,1]$
  
  其中 $G_{i, j}^{\prime}( · )$ 表示第 j 层第 i 个区域的灰度图的梯度图。
计算均值。

对最后一层的每一个区域，分别计算覆盖该区域的所有层的均值。对每个特征都进行此计算，计算方式如图所示：
计算融合特征图。

对每一个特征按上述方式计算均值后即可得到 3 个特征图，对三个特征图进行归一化操作，再相加取均值得到最终特征图。

$F_i = \frac{(f_{1_i} + f_{2_i} + f_{3_i})}{3}$
区域选择。

最终的区域选择为特征图 F 上最小值所对应的区域，接下来在该区域原图的绿蓝暗通道图上选择最亮的像素，返回原图中该像素的 R G B 数值。

Matlab 实现

绿蓝暗通道

function out = gb_dark_channel(input,x)
% x 为滤波窗口大小
G = input(:,:,2);
B = input(:,:,3);
out = ordfilt2(min(G,B),1,ones(x,x),'symmetric');
end

特征 1

function out = f1(input)
out = 1 - mean2(gb_dark_channel(input,15));
end

特征 2

function out = f2(input)
% 此公式和 Matlab 中的 std2 函数的计算结果相同
R = input(:,:,1);
G = input(:,:,2);
B = input(:,:,3);
[m,n,~] = size(input);
u_R = mean2(R) * ones(m,n);
u_G = mean2(G) * ones(m,n);
u_B = mean2(B) * ones(m,n);

sigma_R = (1/(m*n-1) * sum(sum((u_R - R).^2))).^0.5; % sigma_R = std2(R)
sigma_G = (1/(m*n-1) * sum(sum((u_G - G).^2))).^0.5;
sigma_B = (1/(m*n-1) * sum(sum((u_B - B).^2))).^0.5;

out = (sigma_R + sigma_G + sigma_B) / 3;
end

特征 3

function out = f3(input)
I_grey = rgb2gray(input);
[Fx,Fy] = gradient(I_grey);
grad = (Fx .^ 2 + Fy .^ 2).^0.5;
out = mean2(grad);
end

分层的同时计算每层特征

function out = fen_ceng(input,J)
[m,n,~] = size(input);
distx = floor(m / (2^J));
disty = floor(n / (2^J));
out.outimg = cell(2^J,2^J); % 实际上不需要保存图像数据，只保存特征就可以了。
for i = 1:2^J
    for j = 1:2^J
        if i == 2^J 
            out.outimg{i,j} = input((i-1)*distx+1:end,(j-1)*disty+1:j*disty,:);
            out.f1{i,j} = f1(out.outimg{i,j});
            out.f2{i,j} = f2(out.outimg{i,j});
            out.f3{i,j} = f3(out.outimg{i,j});
            continue;
        end
        if j == 2^J
            out.outimg{i,j} = input((i-1)*distx+1:i*distx,(j-1)*disty+1:end,:);
            out.f1{i,j} = f1(out.outimg{i,j});
            out.f2{i,j} = f2(out.outimg{i,j});
            out.f3{i,j} = f3(out.outimg{i,j});
            continue;
        end
        out.outimg{i,j} = input((i-1)*distx+1:i*distx,(j-1)*disty+1:j*disty,:);
        out.f1{i,j} = f1(out.outimg{i,j});
        out.f2{i,j} = f2(out.outimg{i,j});
        out.f3{i,j} = f3(out.outimg{i,j});
    end
end
end

归一化

function [outputs] = Normalization(arg)
    MAX = max(arg(:));
    MIN = min(arg(:));
    outputs = (arg - MIN)./ (MAX - MIN);
end

标记函数

function out = tag(input,x,y,distx,disty)
% 在输入图像的 x，y 位置上画一个 大小为 2倍distx * 2倍disty 大小的红框 
I = im2double(input);

I((x-1)*distx+1:(x+1)*distx,(y-1)*disty+1,1) = 1.0;
I((x-1)*distx+1:(x+1)*distx,(y+1)*disty,1) = 1.0;
I((x-1)*distx+1,(y-1)*disty+1:(y+1)*disty,1) = 1.0;
I((x+1)*distx,(y-1)*disty+1:(y+1)*disty,1) = 1.0;

I((x-1)*distx+1:(x+1)*distx,(y-1)*disty+1,2) = 0;
I((x-1)*distx+1:(x+1)*distx,(y+1)*disty,2) = 0;
I((x-1)*distx+1,(y-1)*disty+1:(y+1)*disty,2) = 0;
I((x+1)*distx,(y-1)*disty+1:(y+1)*disty,2) = 0;

I((x-1)*distx+1:(x+1)*distx,(y-1)*disty+1,3) = 0;
I((x-1)*distx+1:(x+1)*distx,(y+1)*disty,3) =0;
I((x-1)*distx+1,(y-1)*disty+1:(y+1)*disty,3) = 0;
I((x+1)*distx,(y-1)*disty+1:(y+1)*disty,3) = 0;

out = I;
end

主函数

function out = Veiling_light_estimation(input,str)

I = im2double(input);
[m,n,~] = size(input);

J = 6; % 最大层数

for i = 1:J
    layer{i} = fen_ceng(I,i);
end
f1_avg = zeros(2^J,2^J);
f2_avg = zeros(2^J,2^J);
f3_avg = zeros(2^J,2^J);

for i = 1:2^J
    for j = 1:2^J
        tmpf1 = 0;
        tmpf2 = 0;
        tmpf3 = 0;
        for k = J:-1:1
            tmpf1 = tmpf1 + layer{1,k}.f1{max(round(i/(2^(J-k))),1),max(round(j/(2^(J-k))),1)};
            tmpf2 = tmpf2 + layer{1,k}.f2{max(round(i/(2^(J-k))),1),max(round(j/(2^(J-k))),1)};
            tmpf3 = tmpf3 + layer{1,k}.f3{max(round(i/(2^(J-k))),1),max(round(j/(2^(J-k))),1)};
        end
        f1_avg(i,j) = tmpf1 / J;
        f2_avg(i,j) = tmpf2 / J;
        f3_avg(i,j) = tmpf3 / J;
    end
end
f1_alpha = Normalization(f1_avg);
f2_beta = Normalization(f2_avg);
f3_gamma = Normalization(f3_avg);

F = (f1_alpha + f2_beta + f3_gamma) ./ 3;

minF = min(min(F));
[x1,y1] = find(minF == F);
x = floor(median(x1));
y = floor(median(y1));
F(x,y) = 1.0;
% figure;imshow(F); title('论文中的Figure 1 (d)')%  figure 1 (d)


distx = floor(m / (2^J));
disty = floor(n / (2^J));
% figure;imshow(I); 
I_tag = tag(I,x,y,distx,disty);
imtool(I_tag);
% rectangle('Position', [(y-1)*disty+1,(x-1)*distx+1 , disty, distx ],'EdgeColor', 'r');
imwrite(I_tag,[str '.png']);
if x == 2^J && y == 2^J
    patch = I((x-1)*distx+1:end,(y-1)*disty+1:end,:);
elseif x == 2^J
    patch = I((x-1)*distx+1:end,(y-1)*disty+1:y*disty,:);
elseif y == 2^J
    patch = I((x-1)*distx+1:x*distx,(y-1)*disty+1:end,:);
else
    patch = I((x-1)*distx+1:x*distx,(y-1)*disty+1:y*disty,:);
end

patch_gb_dark_channel = gb_dark_channel(patch,3);
[xx,yy] = find(max(max(patch_gb_dark_channel)) == patch_gb_dark_channel);

r = mean2(patch(xx,yy,1));
g = mean2(patch(xx,yy,2));
b = mean2(patch(xx,yy,3));

out = [r,g,b];
end

结果图

先来一些比较好的

再来一些不好的

364