计算机视觉小实例 No.1 基于直方图优化的图像去雾技术

直方图均衡化

直方图均衡化主要用于增强灰度值动态范围偏小的图像的对比度。
什么叫做动态范围，表示图像中所包含的从“最暗”至“最亮”的范围，动态范围越大，所能表现的层次越丰富，所包含的色彩空间也越广。

直方图

直方图是图像的一种统计表达形式。
直方图可以在一定程度上反映图像的概貌性描述，包括图像的灰度范围、灰度分布、整幅图像的亮度均值和阴暗对比度等

f = imread('coins.png');
g = imhist(f, 25);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), bar(g);

上面这个直方图的代表意义就是将灰度分为25份，这张图片中，每份所占的数量（可以看出来，这个图偏暗）。

直方图均衡化

把原始图像的灰度统计直方图变换为均匀分布的形式，这样就增加了像素灰度值的动态范围，从而达到增强图像整体对比度的效果（就是将直方图的每个小柱子均匀分布，所以也叫直方图均匀化）。

主要步骤

1. 预处理
   输入图像，计算该图像直方图。
2. 灰度变换表
   根据输入图像的直方图计算灰度值变化表。
3. 查表变换
   将预处理得到的直方图使用灰度变换表进行查表变换操作，通过遍历每一个像元（像素点），将原始图像灰度值x放入交换表中，可得到该点的新灰度值。

以上面那个图为例

f = imread('coins.png');
g = imhist(f);
subplot(2, 2, 1), imshow(f);
subplot(2, 2, 2), bar(g);
h = histeq(f); % 全局直方图均衡化
g = imhist(h);
subplot(2, 2, 3), imshow(h);
subplot(2, 2, 4), bar(g);

效果还是很明显的，但是在这里，我们可以看到，图片变得更加亮，显得更加真实，细节更多，但是有部分细节已经消失了，主要是这张图，不应该用这个方法进行处理，这张图应该用对比度增强adapthisteq（这个函数其实是局部直方图均衡化-使用CLAHE算法，也叫限制对比度自适应直方图均衡）进行处理（或者用线性滤波器）。

f = imread('coins.png');
subplot(1, 2, 1), imshow(f), title('原图');
g = adapthisteq(f);
subplot(1, 2, 2), imshow(g), title('增强后');

上面这个图就是使用了对比度增强，有点跑题了- -。

f = imread('天安门.jpg');
r = f(:, :, 1);
g = f(:, :, 2);
b = f(:, :, 3);
g = cat(3, histeq(r), histeq(g), histeq(b));
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(g);

效果明显的看下面这张图

$$对rgb通道进行处理$$

 f = imread('天安门.jpg');
hsi = rgb2hsi(f);
h = hsi (:, :, 1);
s = hsi (:, :, 2);
i = hsi (:, :, 3);
g = cat(3, h, s, histeq(i));
g = hsi2rgb(g);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(g);

$$对hsi通道的i通道进行处理$$
可以看出来，这张图的下半部分，去雾效果就非常明显了（是不是很鲜活，很神奇，其实就是将细节放大），但是上半部分出现了光圈- -，

我们只需要使用adapthisteq局部直方图均衡化即可，直接看效果

f = imread('天安门.jpg');
hsi = rgb2hsi(f);
h = hsi (:, :, 1);
s = hsi (:, :, 2);
i = hsi (:, :, 3);
g = cat(3, h, s, adapthisteq(i));
g = hsi2rgb(g);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(g);

可以看到，光圈没有了- -
adapthisteq这个函数还有很多参数，可以直接help adapthisteq查看帮助文档，不过直接用即可adapthisteq(f)。

不过呢，上面的全局均衡化，局部均衡化，只是简单的进行处理，得到的图像偏暗，存在模糊的地方。
我们可以使用imfilter线性滤波来进行处理

f = imread('天安门.jpg');
hsi = rgb2hsi(f);
h = hsi (:, :, 1);
s = hsi (:, :, 2);
i = hsi (:, :, 3);
g = cat(3, h, s, adapthisteq(imfilter(i, fspecial('average'))));
g = hsi2rgb(g);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(g);

效果不是很明显，这个时候我们就可以使用Retinex增强处理的MSR。

function in = retinex(f)
% =========================
% 使用Retinex增强算法来对图像进行增强
% f = retinex(f);
% 注意传入的f必须是rgb格式的
% =========================

% 先提取rgb通道
fr = f(:, :, 1);
fg = f(:, :, 2);
fb = f(:, :, 3);

% 将数据归一化
mr = mat2gray(im2double(fr));
mg = mat2gray(im2double(fg));
mb = mat2gray(im2double(fb));

alpha = 15 ^ 2 / 2;

n = 125;

n1 = floor((n + 1) / 2);

for i = 1 : n
        for j = 1 : n
            b(i, j) = exp(-((i - n1) ^ 2 + (j - n1) ^ 2) / (4 * alpha)) / (pi * alpha);
        end
end

% 卷积滤波
nr1 = imfilter(mr, b, 'conv', 'replicate');
ng1 = imfilter(mg, b, 'conv', 'replicate');
nb1 = imfilter(mb, b, 'conv', 'replicate');

ur1 = log(nr1);
ug1 = log(ng1);
ub1 = log(nb1);

tr1 = log(mr);
tg1 = log(mg);
tb1 = log(mb);

yr1 = (tr1 - ur1) / 3;
yg1 = (tg1 - ug1) / 3;
yb1 = (tb1 - ub1) / 3;

alpha = 80 ^ 2 / 2;

% n = 31;

for i = 1 : n
        for j = 1 : n
            a(i, j) = exp(-((i - n1) ^ 2 + (j - n1) ^ 2) / (4 * alpha)) / (6 * pi * alpha);
        end
end

% 卷积滤波
nr2 = imfilter(mr, a, 'conv', 'replicate');
ng2 = imfilter(mg, a, 'conv', 'replicate');
nb2 = imfilter(mb, a, 'conv', 'replicate');

ur2 = log(nr2);
ug2 = log(ng2);
ub2 = log(nb2);

tr2 = log(mr);
tg2 = log(mg);
tb2 = log(mb);

yr2 = (tr2 - ur2) / 3;
yg2 = (tg2 - ug2) / 3;
yb2 = (tb2 - ub2) / 3;

imshow(yr2);

alpha = 120 ^ 2 / 2;

% n = 501;

for i = 1 : n
        for j = 1 : n
            e(i, j) = exp(-((i - n1) ^ 2 + (j - n1) ^ 2) / (4 * alpha)) / (4 * pi * alpha);
        end
end

% 卷积滤波
nr3 = imfilter(mr, e, 'conv', 'replicate');
ng3 = imfilter(mg, e, 'conv', 'replicate');
nb3 = imfilter(mb, e, 'conv', 'replicate');

ur3 = log(nr3);
ug3 = log(ng3);
ub3 = log(nb3);

tr3 = log(mr);
tg3 = log(mg);
tb3 = log(mb);

yr3 = (tr3 - ur3) / 3;
yg3 = (tg3 - ug3) / 3;
yb3 = (tb3 - ub3) / 3;

dr = yr1 + yr2 + yr3;
dg = yg1 + yg2 + yg3;
db = yb1 + yb2 + yb3;

cr = im2uint8(dr);
cg = im2uint8(dg);
cb = im2uint8(db);

in = cat(3, cr, cg, cb);

end

上面这个是Retinex的MSR函数

f = imread('天安门.jpg');
subplot(1, 2, 1), imshow(f), title('原图');
subplot(1, 2, 2), imshow(retinex(f)), title('MSR处理后');

怎么样，是不是超级明显- -，如果嫌太亮，可以将滤大小调高。

参考文献
1. 《MATLAB计算机视觉与深度学习实战》 刘衍琦、詹福宇、蒋献文、周华英. 电子工业出版社2017年6月出版