1.3 雾霭图像清晰化及其实现
随着计算机技术、传感器技术的飞速发展，数字图像处理在智能交通、侦察勘探、导航定位等领域发挥的作用日益突出。
由于雾霭天气的影响，会使能见度变低，图像采集设备采集的图像一般都会出现比较严重的退化与失真，会对智能交通、航拍测绘、视频导航造成严重的影响。因此，通过现代图像处理的方法提高雾霭天气下拍摄图像的能见度具有重要的作用。
1.3.1 Retinex理论
Retinex（视网膜Retina和大脑皮层Cortex的缩写）理论是一种建立在科学实验和科学分析基础上的基于人类视觉系统（Human Visual System）的图像增强理论。该算法的基本原理模型最早是由Edwin Land（埃德温·兰德）于1971年提出的理论，并在颜色恒常性的基础上提出的一种图像增强方法。Retinex理论的基本内容是：物体的颜色是由物体对长波（红）、中波（绿）和短波（蓝）光线的反射能力决定的，而不是由反射光强度的绝对值决定的；物体的色彩不受光照非均性的影响，具有一致性，即Retinex理论是以色感一致性（颜色恒常性）为基础的。
根据Edwin Land提出的理论，一幅给定的图像S（x,y）分解成两幅不同的图像：反射物体图像R（x,y）和入射光图像L（x,y），其原理示意图如图1.3.1所示。
图1.3.1 Retinex理论示意图
对于给定图像S中的每个点（x,y），用公式可以表示为：
实际上，Retinex理论就是通过图像S来得到物体的反射性质R，也就是去除了入射光L的性质从而得到物体原本该有的样子。
1.3.2 基于Retinex理论的图像增强的基本步骤
步骤一：利用取对数的方法将照射光分量和反射光分量分离，即：
S′（x,y）＝r（x,y）+l（x,y）＝log（R（x,y））+log（L（x,y））
步骤二：用高斯模板对原图像作卷积，即相当于对原图像作低通滤波，得到低通滤波后的图像D（x,y）,F（x,y）表示高斯滤波函数：
D（x, y）＝S（x, y）*F（x, y）
步骤三：在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像G（x,y）：
G（x,y）＝S′（x,y）－log（D（x,y））
步骤四：对G（x,y）取反对数，得到增强后的图像R（x,y）：
R（x,y）＝exp（G（x,y））
步骤五：对R（x,y）作对比度增强，得到最终的结果图像。
1.3.3 多尺度Retinex算法
Jobson D．等人提出了多尺度Retinex算法，多尺度算法（MSR）的基本公式是：
其中，Ri（x,y）是Retinex的输出，i∈R,G,B，表示3个颜色谱带，F（x,y）是高斯滤波函数，Wn表示尺度的权重因子，N表示使用尺度的个数。N＝3，表示彩色图像，i∈R,G,B;N＝1，表示灰度图像。从公式中可以看出：MSR算法的特点是能产生包含色调再现和动态范围压缩这两个特性的输出图像。
在MSR算法的增强过程中，图像可能会因为增加了噪声而造成对图像中的局部区域色彩失真，使得物体的真正颜色效果不能很好地显现出来，从而影响了整体视觉效果。为了弥补这个缺点，一般情况下会采用带色彩恢复因子C的多尺度算法（MSRCR）来解决。带色彩恢复因子C的多尺度算法是在多个固定尺度的基础上考虑色彩不失真恢复的结果，在多尺度Retinex算法过程中，通过引入一个色彩因子C来弥补由于图像局部区域对比度增强而导致的图像颜色失真的缺陷，通常情况下所引入的色彩恢复因子C的表达式为：
其中，Ci表示第i个通道的色彩恢复系数，它的作用是调节3个通道颜色的比例，f（）表示颜色空间的映射函数。带色彩恢复的多尺度Retinex算法通过色彩恢复因子C这个系数来调整原始图像中3个颜色通道之间的比例关系，从而把相对有点暗的区域的信息凸显出来，以达到消除图像色彩失真缺陷的目的。处理后的图像局域对比度得以提高，而且其亮度与真实的场景很相似，图像在人们的视觉感知下显得极其逼真。
1.3.4 实例精讲
【例1.3.1】 基于Retinex理论进行雾霭天气增强的MATLAB程序，读者可结合程序及注释对基于Retinex理论进行雾霭图像清晰化的基本原理进行深入分析。
*********************************************************************************
clear;
close all;
% 读入图像
I＝imread(＇wu.png＇);
% 取输入图像的R分量
R＝I(:, :,1);
[N1, M1]＝size(R);
% 对R分量进行数据转换，并对其取对数
R0＝double(R);
Rlog＝log(R0+1);
% 对R分量进行二维傅里叶变换
Rfft2＝fft2(R0);
% 形成高斯滤波函数
sigma＝250;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
% 对高斯滤波函数进行二维傅里叶变换
Ffft＝fft2(double(F));
% 对R分量与高斯滤波函数进行卷积运算
DR0＝Rfft2.*Ffft;
DR＝ifft2(DR0);
% 在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像
DRdouble＝double(DR);
DRlog＝log(DRdouble+1);
Rr＝Rlog－DRlog;
% 取反对数，得到增强后的图像分量
EXPRr＝exp(Rr);
% 对增强后的图像进行对比度拉伸增强
MIN ＝ min(min(EXPRr));
MAX ＝ max(max(EXPRr));
EXPRr ＝(EXPRr－MIN)/(MAX－MIN);
EXPRr＝adapthisteq(EXPRr);
% 取输入图像的G分量
G＝I(:, :,2);
[N1, M1]＝size(G);
% 对G分量进行数据转换，并对其取对数
G0＝double(G);
Glog＝log(G0+1);
% 对G分量进行二维傅里叶变换
Gfft2＝fft2(G0);
% 形成高斯滤波函数
sigma＝250;
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
% 对高斯滤波函数进行二维傅里叶变换
Ffft＝fft2(double(F));
% 对G分量与高斯滤波函数进行卷积运算
DG0＝Gfft2.*Ffft;
DG＝ifft2(DG0);
% 在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像
DGdouble＝double(DG);
DGlog＝log(DGdouble+1);
Gg＝Glog－DGlog;
% 取反对数，得到增强后的图像分量
EXPGg＝exp(Gg);
% 对增强后的图像进行对比度拉伸增强
MIN ＝ min(min(EXPGg));
MAX ＝ max(max(EXPGg));
EXPGg ＝(EXPGg－MIN)/(MAX－MIN);
EXPGg＝adapthisteq(EXPGg);
% 取输入图像的B分量
B＝I(:, :,3);
[N1, M1]＝size(B);
% 对B分量进行数据转换，并对其取对数
B0＝double(B);
Blog＝log(B0+1);
% 对B分量进行二维傅里叶变换
Bfft2＝fft2(B0);
% 形成高斯滤波函数
sigma＝250;
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
% 对高斯滤波函数进行二维傅里叶变换
Ffft＝fft2(double(F));
% 对B分量与高斯滤波函数进行卷积运算
DB0＝Gfft2.*Ffft;
DB＝ifft2(DB0);
% 在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像
DBdouble＝double(DB);
DBlog＝log(DBdouble+1);
Bb＝Blog－DBlog;
EXPBb＝exp(Bb);
% 对增强后的图像进行对比度拉伸增强
MIN ＝ min(min(EXPBb));
MAX ＝ max(max(EXPBb));
EXPBb ＝(EXPBb－MIN)/(MAX－MIN);
EXPBb＝adapthisteq(EXPBb);
% 对增强后的图像R、G、B分量进行融合
I0(:, :,1)＝EXPRr;
I0(:, :,2)＝EXPGg;
I0(:, :,3)＝EXPBb;
% 显示运行结果
subplot(121), imshow(I);
subplot(122), imshow(I0);
*********************************************************************************
该程序的运行结果如图1.3.2所示。
图1.3.2 例1.3.1的运行结果
【例1.3.2】 基于Retinex理论进行雾霭图像清晰化的MATLAB程序，读者可结合程序及注释对基于Retinex理论进行雾霭图像清晰化的基本原理进行进一步分析。
*********************************************************************************
clear;
close all;
I＝imread(＇wu.png＇);
% 分别取输入图像的R、G、B三个分量，并将其转换为双精度型
R＝I(:, :,1);
G＝I(:, :,2);
B＝I(:, :,3);
R0＝double(R);
G0＝double(G);
B0＝double(B);
[N1, M1]＝size(R);
% 对R分量进行对数变换
Rlog＝log(R0+1);
% 对R分量进行二维傅里叶变换
Rfft2＝fft2(R0);
% 形成高斯滤波函数(sigma＝128)
sigma＝128;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
% 对高斯滤波函数进行二维傅里叶变换
Ffft＝fft2(double(F));
% 对R分量与高斯滤波函数进行卷积运算
DR0＝Rfft2.*Ffft;
DR＝ifft2(DR0);
% 在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像
DRdouble＝double(DR);
DRlog＝log(DRdouble+1);
Rr0＝Rlog－DRlog;
% 形成高斯滤波函数(sigma＝256)
sigma＝256;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
% 对高斯滤波函数进行二维傅里叶变换
Ffft＝fft2(double(F));
% 对R分量与高斯滤波函数进行卷积运算
DR0＝Rfft2.*Ffft;
DR＝ifft2(DR0);
% 在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像
DRdouble＝double(DR);
DRlog＝log(DRdouble+1);
Rr1＝Rlog－DRlog;
% 形成高斯滤波函数(sigma＝512)
sigma＝512;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
% 对高斯滤波函数进行二维傅里叶变换
Ffft＝fft2(double(F));
% 对R分量与高斯滤波函数进行卷积运算
DR0＝Rfft2.*Ffft;
DR＝ifft2(DR0);
% 在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像
DRdouble＝double(DR);
DRlog＝log(DRdouble+1);
Rr2＝Rlog－DRlog;
% 对上述三次增强得到的图像取均值，作为最终增强的图像
Rr＝(1/3)*(Rr0+Rr1+Rr2);
% 定义色彩恢复因子C
a＝125;
II＝imadd(R0, G0);
II＝imadd(II, B0);
Ir＝immultiply(R0, a);
C＝imdivide(Ir, II);
C＝log(C+1);
% 将增强后的R分量乘以色彩恢复因子，并对其进行反对数变换
Rr＝immultiply(C, Rr);
EXPRr＝exp(Rr);
% 对增强后的R分量进行灰度拉伸
MIN ＝ min(min(EXPRr));
MAX ＝ max(max(EXPRr));
EXPRr ＝(EXPRr－MIN)/(MAX－MIN);
EXPRr＝adapthisteq(EXPRr)
[N1, M1]＝size(G);
% 对G分量进行处理，步骤与对R分量处理的步骤相同，可仿照R分量处理的步骤进行理解
G0＝double(G);
Glog＝log(G0+1);
Gfft2＝fft2(G0);
sigma＝128;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
Ffft＝fft2(double(F));
DG0＝Gfft2.*Ffft;
DG＝ifft2(DG0);
DGdouble＝double(DG);
DGlog＝log(DGdouble+1);
Gg0＝Glog－DGlog;
sigma＝256;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
Ffft＝fft2(double(F));
DG0＝Gfft2.*Ffft;
DG＝ifft2(DG0);
DGdouble＝double(DG);
DGlog＝log(DGdouble+1);
Gg1＝Glog－DGlog;
sigma＝512;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
Ffft＝fft2(double(F));
DG0＝Gfft2.*Ffft;
DG＝ifft2(DG0);
DGdouble＝double(DG);
DGlog＝log(DGdouble+1);
Gg2＝Glog－DGlog;
Gg＝(1/3)*(Gg0+Gg1+Gg2);
a＝125;
II＝imadd(R0, G0);
II＝imadd(II, B0);
Ir＝immultiply(R0, a);
C＝imdivide(Ir, II);
C＝log(C+1);
Gg＝immultiply(C, Gg);
EXPGg＝exp(Gg);
MIN ＝ min(min(EXPGg));
MAX ＝ max(max(EXPGg));
EXPGg ＝(EXPGg－MIN)/(MAX－MIN);
EXPGg＝adapthisteq(EXPGg);
% 对B分量进行处理，步骤与对R分量处理的步骤相同，可仿照R分量处理的步骤进行理解
[N1, M1]＝size(B);
B0＝double(B);
Blog＝log(B0+1);
Bfft2＝fft2(B0);
sigma＝128;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
Ffft＝fft2(double(F));
DB0＝Bfft2.*Ffft;
DB＝ifft2(DB0);
DBdouble＝double(DB);
DBlog＝log(DBdouble+1);
Bb0＝Blog－DBlog;
sigma＝256;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
Ffft＝fft2(double(F));
DB0＝Bfft2.*Ffft;
DB＝ifft2(DB0);
DBdouble＝double(DB);
DBlog＝log(DBdouble+1);
Bb1＝Blog－DBlog;
sigma＝512;
F ＝ zeros(N1, M1);
for i＝1:N1
for j＝1:M1
F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
end
end
F ＝ F./(sum(F(:)));
Ffft＝fft2(double(F));
DB0＝Rfft2.*Ffft;
DB＝ifft2(DB0);
DBdouble＝double(DB);
DBlog＝log(DBdouble+1);
Bb2＝Blog－DBlog;
Bb＝(1/3)*(Bb0+Bb1+Bb2);
a＝125;
II＝imadd(R0, G0);
II＝imadd(II, B0);
Ir＝immultiply(R0, a);
C＝imdivide(Ir, II);
C＝log(C+1);
Bb＝immultiply(C, Bb);
EXPBb＝exp(Bb);
MIN ＝ min(min(EXPBb));
MAX ＝ max(max(EXPBb));
EXPBb ＝(EXPBb－MIN)/(MAX－MIN);
EXPBb＝adapthisteq(EXPBb);
% 对增强后的图像R、G、B分量进行融合
I0(:, :,1)＝EXPRr;
I0(:, :,2)＝EXPGg;
I0(:, :,3)＝EXPBb;
% 显示运行结果
subplot(121), imshow(I);
subplot(122), imshow(I0);
*********************************************************************************
该程序的运行结果如图1.3.3所示。
图1.3.3 例1.3.2的运行结果