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

成像模型：

I(x)=S(x)t(x)+(1−t(x))VI(x) = S(x)t(x) + (1 - t(x))V

1

A method for estimating the atmospheric light in the atmospheric scattering model, denoted as V here

Imaging model:

I(x)=S(x)t(x)+(1−t(x))VI(x) = S(x)t(x) + (1 - t(x))V

1

方法简介

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

  1. 分层操作。

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

    image-20210207200553809

  2. 计算特征。

    对每一层的每一块分别计算下面三个特征。

    • 特征 1：

      μi,jdGB=1∣Ωi,j∣∑x∈Ωi,jCi,jdGB(x)\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)

      f1i,j=1−μi,jdGB∈[0,1]f_{1_{i, j}}=1-\mu_{i, j}^{d^{G B}} \in[0,1]

      其中 Ci,jdGB(x)C_{i,j}^{d^{G B}}(x) 表示仅使用 绿色和蓝色通道计算的暗通道图。∣Ωi,j∣\left|\Omega_{i, j}\right| 表示第 j 层 第 i 个区域的像素点个数。

    • 特征 2 ：

      σi,jc=1∣Ωi,j∣−1∑x∈Ωi,j(μi,jc−Ii,jc(x))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}}

      f2i,j=σi,jR+σi,jG+σi,jB3∈[0,1]f_{2_{i, j}}=\frac{\sigma_{i, j}^{R}+\sigma_{i, j}^{G}+\sigma_{i, j}^{B}}{3} \in[0,1]

      其中 μi,jc\mu_{i, j}^{c} 代表 第 j 层 第 i 个区域的像素均值，I 为原图。

    • 特征 3 ：

      f3i,j=1∣Ωi,j∣∑x∈Ωi,j(Gi,j′(x))∈[0,1]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]

      其中 Gi,j′(⋅)G_{i, j}^{\prime}( · ) 表示 第 j 层 第 i 个区域 的灰度图的梯度图。

  3. 计算均值。

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

    image-20210207203350713

  4. 计算融合特征图。

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

    Fi=(f1i+f2i+f3i)3F_i = \frac{(f_{1_i} + f_{2_i} + f_{3_i})}{3}

  5. 区域选择。

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

Method Overview

The process of estimating V is a vertical hierarchical search technique, which mainly consists of the following steps:

  1. Layering operation.

    Layering is not actual image layering; instead, the image is divided into blocks using segmentation at different scales to compute features. The i-th layer divides the image into 4i4^i blocks.

    image-20210207200553809

  2. Computing features.

    For each block of each layer, the following three features are computed.

    • Feature 1:

      μi,jdGB=1∣Ωi,j∣∑x∈Ωi,jCi,jdGB(x)\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)

      f1i,j=1−μi,jdGB∈[0,1]f_{1_{i, j}}=1-\mu_{i, j}^{d^{G B}} \in[0,1]

      where Ci,jdGB(x)C_{i,j}^{d^{G B}}(x) denotes the dark channel image computed using only the green and blue channels, and ∣Ωi,j∣\left|\Omega_{i, j}\right| denotes the number of pixels in the i-th region of the j-th layer.

    • Feature 2:

      σi,jc=1∣Ωi,j∣−1∑x∈Ωi,j(μi,jc−Ii,jc(x))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}}

      f2i,j=σi,jR+σi,jG+σi,jB3∈[0,1]f_{2_{i, j}}=\frac{\sigma_{i, j}^{R}+\sigma_{i, j}^{G}+\sigma_{i, j}^{B}}{3} \in[0,1]

      where μi,jc\mu_{i, j}^{c} represents the pixel mean of the i-th region of the j-th layer, and I is the original image.

    • Feature 3:

      f3i,j=1∣Ωi,j∣∑x∈Ωi,j(Gi,j′(x))∈[0,1]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]

      where Gi,j′(⋅)G_{i, j}^{\prime}( · ) denotes the gradient map of the grayscale image of the i-th region of the j-th layer.

  3. Computing the mean.

    For each region of the last layer, compute the mean over all layers covering that region. This is done for each feature, as shown in the figure:

    image-20210207203350713

  4. Computing the fused feature map.

    After computing the mean for each feature in the above way, three feature maps are obtained. The three feature maps are normalized, then added and averaged to obtain the final feature map.

    Fi=(f1i+f2i+f3i)3F_i = \frac{(f_{1_i} + f_{2_i} + f_{3_i})}{3}

  5. Region selection.

    The final region selection is the region corresponding to the minimum value on the feature map F. Next, the brightest pixel is selected on the green-blue dark channel image of the original image in that region, and the R G B values of that pixel in the original image are returned.

Matlab 实现

绿蓝暗通道

Matlab Implementation

Green-Blue Dark Channel

1
2
3
4
5
6
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

Feature 1

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

特征 2

Feature 2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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

Feature 3

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

分层的同时计算每层特征

Computing Features of Each Layer While Layering

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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

归一化

Normalization

1
2
3
4
5
function [outputs] = Normalization(arg)
    MAX = max(arg(:));
    MIN = min(arg(:));
    outputs = (arg - MIN)./ (MAX - MIN);
end

标记函数

Marking Function

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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

主函数

Main Function

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
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

结果图

先来一些比较好的

5

2

5

5

5

再来一些不好的

364

5

5

5

5

Results

First, Some Better Ones

5

2

5

5

5

Then Some Worse Ones

364

5

5

5

5