注意:这篇文章上次更新于1302天前,文章内容可能已经过时。
求大气散射模型中大气光的一种方法,这里记为V
成像模型:
方法简介
求 V 的过程是一种垂直分层的搜索技术,主要有以下几个步骤:
-
分层操作。
分层并不是真的图像分层,而是用不同尺度的分割方式对图像进行分块计算特征。第 i 层把图像分为 块。
-
计算特征。
对每一层的每一块分别计算下面三个特征。
-
特征 1:
其中 {i,j}{d{G B}}(x)$ 表示仅使用 绿色和蓝色通道计算的暗通道图。 表示第 j 层 第 i 个区域的像素点个数。
-
特征 2 :
其中 代表 第 j 层 第 i 个区域的像素均值,I 为原图。
-
特征 3 :
其中 {i, j}^{\prime}( · )$ 表示 第 j 层 第 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