【MATLAB图像处理实用案例详解（1）】—— 基于直方图优化的图像去雾技术

一、背景意义

雾霾天气往往会给人类的生产和生活带来极大不便，也大大增加了交通事故的发生概率。一般而言，在恶劣天气(如雾天、雨天等)条件下，户外景物图像的对比度和颜色会被改变或退化，图像中蕴含的许多特征也会被覆盖或模糊，这会导致某些视觉系统(如电子卡口、门禁监控等)无法正常工作。因此，从在雾霾天气下采集的退化图像中复原和增强景物的细节信息具有重要的现实意义。数字图像处理技术已被广泛应用于科学和工程领域，如地形分类系统、户外监控系统、自动导航系统等。为了保证视觉系统全天候正常工作，就必须使视觉系统适应各种天气状况。
本案例将展开对雾霾天气下的图像清晰化技术的讨论，雾天天气下的图像清晰化技术也有可能对其他恶劣天气下的图像清晰化技术的发展起到促进作用。

二、理论基础

2.1 空域图像増强

图像增强是指按特定的需要突出一幅图像中的某些信息，并同时削弱或去除某些不需要的信息的处理技术。
图像增强的主要作用是相对于原来的图像，处理后的图像能更加有效地满足某些特定应用的要求。

根据图像处理空间的不同，图像增强基本上可分为两大类：频域处理法、空域处理法。

频域处理法的基础是卷积定理，其通过进行某种图像变换（如傅里叶变换、小波变换等）得到频域结果并进行修改的方法来实现对图像的增强处理。
空域处理法是直接对图像中的像素进行处理，一般是以图像灰度映射变换为基础并且根据图像增强的目标来采用所需的映射变换，常见的图像对比度增强、图像的灰度层次优化等处理均属于空域处理法。

2.2 直方图均衡化

直方图是图像的一种统计表达形式。对于一幅灰度图像来说，其灰度统计直方图可以反映该图像中不同灰度级出现的统计情况。一般而言，图像的视觉效果和其直方图有对应关系，通过调整或变换其直方图的形状会对图像显示效果有很大影响。

直方图均衡化主要用于增强灰度值动态范围偏小的图像的对比度。该方法的基本思想是把原始图像的灰度统计直方图变换为均匀分布的形式，这样就增加了像素灰度值的动态范围，从而达到增强图像整体对比度的效果。
直方图均衡化也叫作直方图均匀化，就是把给定图像的直方图分布变换成均匀分布的直方图，是一种较为常用的灰度增强算法。
直方图均衡化概括起来包括以下三个主要步骤：

预处理。输入图像，计算该图像直方图。
灰度变换表。根据输入图像的直方图计算灰度值变换表。
查表变换。执行变换x’=H(x) ，表示对步骤1中得到的直方图使用步骤2 得到的灰度值变换表进行查表变换操作，通过遍历整幅图像的每一个像元，将原始图像灰度值x 放入变换表H(x) 中，可得到变换后的新灰度值x’

图像经直方图均衡化后，将会包含更多的信息量，进而能突出某些图像特征。具有均匀分布直方图的图像，其信息量H最大。

为进行图像去雾实验，采用全局直方图均衡化、局部直方图均衡化算法来进行图像去雾实验，并选择Retinex 增强算法来作为直方图去雾算法的延伸。

三、方法选择

3.1 全局直方图算法

全局直方图均衡算法对于整幅图像的像素使用相同的直方图变换，对于那些像素值分布比较均衡的图像来说，算法的效果很好。如果图像中包括明显比图像其他区域暗或者亮的部分，在这些部分的对比度将得不到有效的增强。
直接采用整体直方图均衡化对于色彩均衡的图像的确也有较好的效果，残留的雾感较前面的方法较浓，但也去除了大部分雾。
效果展示：

全局直方图算法核心代码如下：

function In = RemoveFogByGlobalHisteq(I, flag)if nargin < 2flag = 1;
end
R = I(:,:,1);
G = I(:,:,2);
B = I(:,:,3);
M = histeq(R);
N = histeq(G);
L = histeq(B);
In = cat(3, M, N, L);
if flagfigure;subplot(2, 2, 1); imshow(I); title('原图像', 'FontWeight', 'Bold');subplot(2, 2, 2); imshow(In); title('处理后的图像', 'FontWeight', 'Bold');Q = rgb2gray(I);W = rgb2gray(In);subplot(2, 2, 3); imhist(Q, 64); title('原灰度直方图', 'FontWeight', 'Bold');subplot(2, 2, 4); imhist(W, 64); title('处理后的灰度直方图', 'FontWeight', 'Bold');
end

3.2 局部直方图算法

可以看到，全局直方图均衡化增强只是将原图像的直方图进行了均衡化，未能有效保持原始图像的局部特征，容易出现色彩失真问题。通过局部直方图处理可以在一定层度上保持原有图像的局部特征，提高去雾的效果。所谓的局部直方图处理就是选择固定尺寸的滑动窗口作用于原始图像来进行局部的处理，通过对RGB三通道进行局部直方图均衡化，然后再整合到新的图像中。
效果展示：

局部直方图算法核心代码如下：

function In = RemoveFogByLocalHisteq(I, flag)g1 = GetLocalHisteq(I(:, :, 1));
g2 = GetLocalHisteq(I(:, :, 2));
g3 = GetLocalHisteq(I(:, :, 3));In = cat(3, g1, g2, g3);if flagfigure;subplot(2, 2, 1); imshow(I); title('原图像', 'FontWeight', 'Bold');subplot(2, 2, 2); imshow(In); title('处理后的图像', 'FontWeight', 'Bold');Q = rgb2gray(I);W = rgb2gray(In);subplot(2, 2, 3); imhist(Q, 64); title('原灰度直方图', 'FontWeight', 'Bold');subplot(2, 2, 4); imhist(W, 64); title('处理后的灰度直方图', 'FontWeight', 'Bold');
endfunction g = GetLocalHisteq(I)
% 对灰度图像，进行局部直方图均衡化
% 输入参数：
%  I——图像矩阵
% 输出参数：
%  g——结果图像
% 调用库函数adapthisteq，执行局部均衡化增强
g = adapthisteq(I,'clipLimit',0.02,'Distribution','rayleigh');

3.3 Retinex算法

Retinex理论主要包含了两个方面的内容：物体的颜色是由物体对长波、中波和短波光线的反射能力决定的，而不是由反射光强度的绝对值决定的；物体的色彩不受光照非均匀性的影响，具有一致性。人眼感知物体的亮度取决于环境的照明和物体表面对照射光的反射，其数学表达式为：I(x,y)=L(x,y)∗R(x,y)I(x,y)=L(x,y)*R(x,y)I(x,y)=L(x,y)∗R(x,y)，其中I(x,y)I(x,y)I(x,y)为摄像机接收到的图像信息；L(x,y)L(x,y)L(x,y)代表环境光的照射分量；R(x,y)R(x,y)R(x,y)表示携带图像细节信息的目标物体反射分量。
对该式取对数有
Log(R(x,y))=Log(I(x,y))−Log(L(x,y))Log(R(x,y))=Log(I(x,y))-Log(L(x,y))Log(R(x,y))=Log(I(x,y))−Log(L(x,y))；
Retinex算法的提出者指出这个L(x,y)可以通过对图像数据I(x,y)进行高斯模糊得到。通过对数操作减去这个L(x,y)就可以得到去雾后的增强图像。
实现步骤为：
（1）输入原始图像数据I(x,y)，高斯模糊的半径；
（2）计算原始图像按指定尺寸进行模糊后的图像L(x,y)；
（3）将Log[R(x,y)]量化为0到255范围的像素值，作为最终的输出。

量化的方法有很多，最简单的例子是计算出Log(R(x,y))的最大值Max和最小值Min，然后对每一个值Value，进行线性量化(线性拉伸/映射到全域)。
效果展示：

Retinex算法核心代码如下：

function In = RemoveFogByRetinex(f, flag)
% 用Retinex的MSR实现图像去雾
% 输入参数：
%  f——图像矩阵
%  flag——显示标记
% 输出参数：
%  In——结果图像if nargin < 2flag = 1;
end
%提取图像的R、G、B分量
fr = f(:, :, 1);
fg = f(:, :, 2);
fb = f(:, :, 3);
%数据类型归一化
mr = mat2gray(im2double(fr));
mg = mat2gray(im2double(fg));
mb = mat2gray(im2double(fb));
%定义alpha参数
alpha = 1200;
%定义模板大小
n = 128;
%计算中心
n1 = floor((n+1)/2);
for i = 1:nfor 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;
%定义beta参数
beta = 55;
%定义模板大小
x = 32;
%计算中心
x1 = floor((n+1)/2);
for i = 1:nfor j = 1:n%高斯函数a(i,j)  = exp(-((i-n1)^2+(j-n1)^2)/(4*beta))/(6*pi*beta);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;
%定义eta参数
eta = 13944.5;
%定义模板大小
l = 500;
%计算中心
l1 = floor((n+1)/2);
for i = 1:nfor j = 1:n%高斯函数e(i,j)  = exp(-((i-n1)^2+(j-n1)^2)/(4*eta))/(4*pi*eta);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);
%结果显示
if flagfigure;subplot(2, 2, 1); imshow(f); title('原图像', 'FontWeight', 'Bold');subplot(2, 2, 2); imshow(In); title('处理后的图像', 'FontWeight', 'Bold');% 灰度化，用于计算直方图Q = rgb2gray(f);M = rgb2gray(In);subplot(2, 2, 3); imhist(Q, 64); title('原灰度直方图', 'FontWeight', 'Bold');subplot(2, 2, 4); imhist(M, 64); title('处理后的灰度直方图', 'FontWeight', 'Bold');
end

四、效果演示

GUI原始界面：

载入图像后的GUI界面：

五、完整代码

ok，以上便是全部内容了，相信你完全可以根据源代码进行图像去雾的GUI设计，当然，如果需要源代码，可以参考资源：
基于直方图优化的图像去雾技术(不带GUI界面)：https://download.csdn.net/download/didi_ya/87373523
基于直方图优化的图像去雾技术(带GUI界面)：https://download.csdn.net/download/didi_ya/87373532
，如果遇到任何运行问题，请私信博主，看到后会免费答疑解惑，当然也可跟你你的需求修改GUI界面。