《数字图像工程案例》 --图像去雾处理

1、雾霭的形成机理

雾实际上是由悬浮颗粒在大气中的微小液滴构成的气溶胶，常呈现乳白色，其底部位于地球表面，所以也可以看作是接近地面的云。霭其实跟雾区别不大，它的一种解释是轻雾，多呈现灰白色，与雾的颜色十分接近。广义的雾包括雾、霾、沙尘、烟等一切导致视觉效果受限的物理现象。由于雾的存在，户外图像质量降低，如果不处理，往往满足不了相关研究、应用的要求。在雾的影响下，经过物体表面的光被大气中的颗粒物吸收和反射，导致获取的图像质量差，细节模糊、色彩暗淡。

2、设计步骤

（一）基于直方图均衡化的图像去雾算法

1、转换为灰度图像后对图像进行直方图均衡
由于直方图均衡仅限于灰度图像，所以我们将输入的RGB图像先转变为灰度图像，再进行图像增强。简要过程如图2-1所示。

图2-1 RGB转灰度图像进行直方图均衡化过程

（1） 读入图像，将彩色图像进行降维转换成灰度图像；

（2） 对灰度图像的直方图进行均衡化处理；

（3） 输出均衡化后的灰色图像。

接下来我们考虑想要得到彩色图像，于是从不同空间对图像进行均衡化处理。

2.1、在RGB空间对图像进行直方图均衡

（1） 读入图像，将图像存储于RGB空间；

（2） 分别对RGB空间的R、G、B三个分量的灰度直方图进行均衡化处理；

（3） 输出图像。

2.2、在HSV空间对图像进行直方图均衡

（1） 读入图像，将图像由RGB空间转换到HSV空间；

（2） 对HSV空间饱和度和亮度分量（S、V分量）的灰度直方图进行均衡化处理；

（3） 将图像由HSV空间转换到RGB空间并输出。

图2-3 在HSI空间对图像进行直方图均衡

2.3、在YCrCb空间对图像进行直方图均衡

（1） 读入图像，将图像由RGB空间转换到YCrCb空间；

（2） 对YCrCb空间亮度分量（Y分量）的灰度直方图进行均衡化处理；

（3） 将图像由YCrCb空间转换到RGB空间并输出图像。

图2-4 在YCrCb空间对图像进行直方图均衡

3基于暗通道先验的图像去雾算法

根据何凯明的暗通道先验算法原理，大致的流程如图2-6所示。

图2-6 暗通道先验过程图

（1）根据原始有雾图像求暗通道

用for循环求出每个像素RGB分量中的最小值，存入一副和原始图像大小相同的灰度图中，然后再对这副灰度图使用ordfilt2进行最小值滤波，滤波的半径由窗口大小决定。

（2）求解全局大气光照

根据原理应该选择暗通道内图像总像素点个数千分之一个最亮的像素点，并记录这些像素点坐标，再根据这些点的坐标分别在原图像I的三个通道内找到这些像素点并加和。
考虑到这个思路在MATLAB实现中在寻找这些像素点时，会找到不止1000点，因为在图像中有很多值相等的像素点；为了方便计算，在求A时取最亮值的点集中的某一个点。

（3）求透射率t（x）

首先设定原始程序中去雾系数w=0.95，根据公式先对透射率有一个预估；接下来利用导向滤波（Guided Filter）来对透射率图。

以原始图像的其中一个通道为引导图，使得透射率图的细节更加准确。而为了提高计算速度，在导向滤波中使用盒子滤波器（Box filter），其本质是通过下采样减少像素点，计算后再进行上采样恢复到原有的尺寸大小。

（4）求解无雾图像

根据公式

求解出去雾的图像即可。在初始程序中采用

为标准计算。

4、设计结果与分析
（一）基于直方图均衡化的图像去雾算法
直方图均衡化是图像增强中的一种常用方法，该算法以概率论为基础，基于空间不变思想，运用灰度运算来实现直方图的变换，从而达到增强图像的目的，更适合于景物深度变化很小的图像。分为全局直方图均衡化算法和局部直方图均衡化算法（在本次实验中局部均衡方法去雾效果比较差，所以不对该方法进行赘述）。

全局直方图均衡化算法实现简单，对单景深图像的复原效果较好，但对于场景深度多变的图像的复原效果好，但是对于场景深度多变的额图像，这种方法就很难反映图像中局部景物深度的变化。

采用局部增强方法可以使图像的每一个区域的对比度都能得到很大的改善，但是同时这种方法也将灰度变化缓慢的区域，如天空等也进行了误增强，从而导致景物影像因许多噪声的干扰而表现得不自然。

原始有雾图像如图3-1所示。

(a) 有雾图像1

(b) 有雾图像2

图3-1 有雾图像

对全局均衡结果图的效果进行分析（程序aft.m）：

对有雾图像RGB通道进行均衡化处理后融合的图像效果最好，清晰度较高，且颜色较符合地物颜色，但天空处有较多的块状噪声，且树叶产生了较大的畸变，结果如图3-1。

图3-1 对RGB通道进行均衡化处理结果

图3-2 对RGB通道进行均衡化结果及其直方图前后对比

将只对HSI空间的亮度分量[I通道]和对HSI空间的亮度、饱和度分量[I、S通道]进行均衡化，通道进行的均衡化处理的图像颜色对比度明显发生的变化，第二种颜色过于鲜艳而发生畸变，第一种较RGB均衡方式的颜色对比度较差，且树叶处畸变较差结果如图3-3。

图3-3 对HSI通道进行均衡化处理结果

图3-4 对HSI通道进行均衡化前后结果对比及其直方图

对YCrCb空间的亮度分量[Y通道]进行均衡处理，所得图像在天空处处理的较差，有大片空白，但是对于树叶的处理较其他几种方式都好，结果如图3-5。

图3-5 对YCrCb通道进行均衡化处理结果

图3-6 对YCrCb通道进行均衡化结果前后对比及其直方图
当考虑将RGB图像转变为灰度图像在进行直方图均衡化的结果如图3-7所示，

图3-7 对灰度进行均衡化结果前后对比及其直方图

图3-8 对灰度进行均衡化结果前后对比及其直方图

如果忽略灰度图直方图均衡化的视觉效果问题，在去雾效果方面灰度变换的结果还是比较好的。

接着再使用有雾图像2的直方图均衡化，结果如图3-9所示，

图3-8 对灰度进行均衡化结果前后对比及其直方图

由于有雾图像2在前景处没有雾的影响，在整体图像增强后虽然也起到了一定去雾的目的，但是图像效果不太好。

总的来说，直方图均衡化的方法对灰度变化不大的树叶和天空处理效果较差，但是对于建筑物等灰度变化较大的区域效果还是不错的。

直方图均衡化通过使用累积函数对灰度值进行“调整”以实现对比度的增强。具体说来，即把原图像对应的灰度直方图从比较集中的灰度区间通过一定的转换变成在全部灰度范围内的均匀分布。就是对图像进行非线性拉伸，重新分配图像像素值，使一定灰度范围内的像素数量大致相同。

通过使用累积函数对灰度值进行“调整”以实现对比度的增强，一定程度上使暗的地方更暗、亮的地方更亮，从而在视觉上达到去雾的效果。

但是上面的去雾效果仅仅局限于灰度图片，而不能直接对彩色图像进行处理。经过上述直方图处理后，达到均衡化图像的目的，但是图像会在原图基础上失真，即整幅图像的颜色发生改变，有时会严重偏离去雾的初衷。

（二）基于暗通道先验的图像去雾算法

1、实验过程结果图

（1）根据原始有雾图像求暗通道

图4-1 有雾图像1的暗通道图像

图4-2 有雾图像2的暗通道图像

通过两幅图的暗通道图像尤其是有雾图像2中可以看出在近距离的树叶处非常的暗，通过该暗通道图像去估算大气光照和透射率，建立模型就可以很好地在去雾的同时尽可能保留无雾区域的特征细节。

（2）求解全局大气光照
求解有图像1和2的大气光照数A：

（3）求透射率

图4-3 图像1的透射率图

图4-4 图像2的透射率图

通过两幅图像透射率图可以看出利用导向滤波函数求得的透射率图像较为精细，也符合大气物理退化模型中的情况。

（4）求解无雾图像

图4-5 图像1的去雾图像

有雾图像1的去雾结果相较于使用直方图均衡化中效果最好的对于RGB通道进行直方图均衡化处理结果来说，去雾的程度更大，颜色更加真实，特征细节处的保持的也较好。去雾后的图像会比原始的暗，因此在处理完后最好进行一定的曝光增强，但由于视觉上的效果还不错，我们选择不进行其他操作以免影响去雾效果的观察。

图4-6 图像2的去雾图像

图4-6更加能够突出暗通道先验算法去雾的优点，对于无雾区域的处理将区别于有雾区域，可以减小本身不需要去雾图像区域的改变。

2、影响暗通道先验算法去雾效果的参数研究
计算的标准参数值如下表所示：
1 去雾系数ω 0.95
2 滤波窗口大小（形成暗通道的滤波半径）r 9
3 透射率窗口大小（导向半径）R 80
4 透射率下限值 0.1
5 eps

（1）去雾系数ω

粗透射率为

越大，去雾程度越大，图像对比度显得越暗；反之越小，去雾越不彻底，图像越泛白。一般取

5、对比两种去雾算法

直方图均衡化

暗通道先验

优点	
对去雾图像没有要求	不仅在物理上有效，而且能够处理包括在雾浓度很大情况下的远距离物体；可以减小无雾区域的变样

缺点	去雾后的图像或为灰色图像或产生颜色失真，去雾效果不稳定	得到的去雾图像比较暗，无法处理大面积天空或类似区域


虽然原始图像受雾影响模糊不清、颜色不真实的图像，但经过暗通道先验算法处理后在细节上更清晰，颜色更真实，视觉效果更令人满意。这种算法适
合以上各种场景的图像去雾，无论雾的分布均匀与否，雾的浓度薄或浓，场景中深度变坏如何，都能有效地去除图像中雾，有效地提高了图像的清晰度。
从而很大程度上提高了雾天的能见度。但是该算法存在一定的局限性，暗原色先验是一种统计的结果，是对大量户外无雾照片的统计结果，如果目标
场景内在的就和大气光类似，比如雪地、白色背景墙、大海等，则由于前提条件就不正确，因此一般无法获得满意的效果，而对于一般的风景照片这个
算法能处理的不错。另外，由于景物退化与场景深度呈非线性关系，由此带来的最大问题是很难保证建立的景物退化模型的正确性和宽适性。目前，大
多数的图像复原方法都建立在大气散射模型的基础上，并受到了此类模型的限制。

附代码

1、aft.m
1.%利用直方图均衡化去雾
2.
3.%===========================================================
4.
5.close all;
6.
7.clear all;
8.
9.wu=imread('D:\学习\2023\专业\数字图像工程\图片1.png');
10.
11.%--------------------------------------------------------------------------
12.
13.%对饱和度与亮度进行直方图均衡处理
14.
15.%h：色度
16.
17.%s：饱和度
18.
19.%v：亮度
20.
21.hsv=rgb2hsv(wu);
22.
23.wu1=rgb2gray(wu);
24.
25.h=hsv(:,:,1);
26.
27.s=hsv(:,:,2);
28.
29.v=hsv(:,:,3);
30.
31.S=histeq(s);
32.
33.V=histeq(v);
34.
35.result_hsv=hsv2rgb(h,S,V);
36.
37.%--------------------------------------------------------------------------
38.
39.% 对RGB每个通道进行直方图均衡处理
40.
41.% R：红色
42.
43.% G：绿色
44.
45.% B：蓝色
46.
47.r=wu(:,:,1);
48.
49.g=wu(:,:,2);
50.
51.b=wu(:,:,3);
52.
53.R=histeq(r);
54.
55.G=histeq(g);
56.
57.B=histeq(b);
58.
59.result_rgb=cat(3,R,G,B);%将RGB三个分量放入一个矩阵
60.
61.result_rgb1=rgb2gray(result_rgb);
62.
63.%--------------------------------------------------------------------------
64.
65.% 对YCbCr的亮度进行直方图均衡处理
66.
67.%y：颜色的亮度，
68.
69.%Cb：蓝色的浓度偏移成分
70.
71.%Cr：红色的浓度偏移成分
72.
73.ycbcr=rgb2ycbcr(wu);
74.
75.y=ycbcr(:,:,1);
76.
77.cb=ycbcr(:,:,2);
78.
79.cr=ycbcr(:,:,3);
80.
81.Y=histeq(y);
82.
83.result_ycbcr=ycbcr2rgb(cat(3,Y,cb,cr));%%将RGB三个分量放入一个矩阵
84.
85.result_ycbcr1=rgb2gray(result_ycbcr);
86.
87.%将彩色图像转换为灰度图像后，再做直方图均衡处理
88.
89.whole=rgb2gray(wu);
90.
91.W=histeq(whole);
92.
93.%--------------------------------------------------------------------------
94.
95.% 结果显示
96.
97.figure(1)
98.
99.subplot(2,2,1),imshow(wu),title('原始图像')
100.
101.subplot(2,2,2),imshow(result_hsv),title('HSV')
102.
103.subplot(2,2,3),imhist(wu1);
104.
105.subplot(2,2,4),imhist(result_hsv);
106.
107.figure(2)
108.
109.subplot(2,2,1),imshow(wu),title('原始图像')
110.
111.subplot(2,2,2),imshow(result_rgb),title('RGB')
112.
113.subplot(2,2,3),imhist(wu1);
114.
115.subplot(2,2,4),imhist(result_rgb1);
116.
117.figure(3)
118.
119.subplot(2,2,1),imshow(wu),title('原始图像')
120.
121.subplot(2,2,2),imshow(result_ycbcr),title('YCbCr')
122.
123.subplot(2,2,3),imhist(wu1);
124.
125.subplot(2,2,4),imhist(result_ycbcr1);
126.
127.figure(4)
128.
129.subplot(2,2,1),imshow(wu),title('原始图像')
130.
131.subplot(2,2,2),imshow(W),title('灰度图像直方图均衡化')
132.
133.subplot(2,2,3),imhist(wu1);
134.
135.subplot(2,2,4),imhist(W);

2、antongd.m
1.%暗通道去雾算法
2.
3.%===========================================================
4.
5.%暗原色去雾算法是建立在户外自然场景暗通道优先法则的基础上的去雾方法
6.
7.%其实就是解一个方程：I(x)=J(x)t(x)+A(1-t(x))
8.
9.%其中I(x)是受到雾气污染的图像,J(x)是我们需要求的去雾后的图像
10.
11.%t(x)是天空中云层的透射分布率，A是天空的亮度
12.
13.%===========================================================
14.
15.close all;
16.
17.clear all;
18.
19.I = double('D:\学习\2023\专业\数字图像工程\图片1.png')/255;
20.
21.[h,w,c] = size(I); %获取图像大小
22.
23.w0 = 0.95; %去雾系数
24.
25.dehaze = zeros(h,w,c); %初始化结果图像
26.
27.win_dark = zeros(h,w); %初始化暗影通道图像
28.
29.%--------------------------------------------------------------------------
30.
31.%计算暗通道
32.
33.%将三个通道中最暗的值赋给win_dark(i,j),使得三维图变成了二维图
34.
35.for i=1:h
36.
37.for j=1:w
38.
39.win_dark(i,j) = min(I(i,j,:));
40.
41.end
42.
43.end
44.
45.win_dark = ordfilt2(win_dark,1,ones(9,9),'symmetric'); %9*9最小值滤波
46.
47.%--------------------------------------------------------------------------
48.
49.%计算大气亮度A
50.
51.dark_channel = win_dark;
52.
53.A = max(max(dark_channel));
54.
55.[i0,j0] = find(dark_channel==A);
56.
57.i = i0(1);
58.
59.j = j0(1);
60.
61.A = mean(I(i,j,:));
62.
63.%--------------------------------------------------------------------------
64.
65.%计算透射率t(x)
66.
67.transmission = 1-w0*win_dark/A; %透射率预估
68.
69.gray_I = I(:,:,3); %这里gray_I可以是RGB图像中任何一个通道
70.
71.p = transmission; %透射率图
72.
73.r = 50;
74.
75.eps = 10^-3;
76.
77.transmission_filter = guidedfilter(gray_I,p,r,eps);
78.
79.%用guided filter对trasmission map做soft matting获得更为精细的透射率图
80.
81.t0=0.1; %设置阈值
82.
83.t1 = max(t0,transmission_filter);
84.
85.%--------------------------------------------------------------------------
86.
87.for i=1:c
88.
89.for j=1:h
90.
91.for l=1:w
92.
93.dehaze(j,l,i)=(I(j,l,i)-A)/t1(j,l)+A; %无雾图像恢复
94.
95.end
96.
97.end
98.
99.end
100.
101.%--------------------------------------------------------------------------
102.
103.%显示图像
104.
105.figure
106.
107.subplot(121)
108.
109.imshow(win_dark);title('暗通道')
110.
111.subplot(122)
112.
113.imshow(t1);title('透射率图')
114.
115.figure
116.
117.subplot(121)
118.
119.imshow(I);title('去雾前')
120.
121.subplot(122)
122.
123.imshow(dehaze);title('去雾后')

3、guidedfilter.m
1.%guided filter（导向滤波函数）
2.
3.%===========================================================
4.
5.%引导图：I（灰度图/单通道图像）
6.
7.%输入图像：p（灰度图/单通道图像）
8.
9.%本地窗口半径：r
10.
11.%正规化参数：eps
12.
13.%==========================================================================
14.
15.function q = guidefilter(I,p,r,eps)
16.
17.[hei,wid] = size(I);
18.
19.N = boxfilter(ones(hei,wid),r); % the size of each local patch;
20.
21.%N=(2r+1)^2 except for boundary pixels.
22.
23.mean_I = boxfilter(I,r)./N;
24.
25.mean_p = boxfilter(p,r)./N;
26.
27.mean_Ip = boxfilter(I.*p,r)./N;
28.
29.cov_Ip = mean_Ip-mean_I.*mean_p; %协方差
30.
31.mean_II = boxfilter(I.*I,r)./N;
32.
33.var_I = mean_II-mean_I.*mean_I; %方差
34.
35.a = cov_Ip./(var_I+eps);
36.
37.b = mean_p-a.*mean_I;
38.
39.mean_a = boxfilter(a,r)./N;
40.
41.mean_b = boxfilter(b,r)./N;
42.
43.q = mean_a.*I+mean_b;
44.
45.end

4、boxfilter.m
1.%boxfilter（盒子滤波函数）
2.
3.%===========================================================
4.
5.%输入：imSrc
6.
7.%给定的滑动窗口大小：r
8.
9.%===========================================================
10.
11.function imDst = boxfilter(imSrc,r)
12.
13.[hei,wid] = size(imSrc);
14.
15.imDst = zeros(size(imSrc)); %初始化
16.
17.%--------------------------------------------------------------------------
18.
19.imCum = cumsum(imSrc,1); %y轴累计求和，将每一行的数值一次累加到下一行
20.
21.%y轴方向差异
22.
23.imDst(1:r+1,:) = imCum(1+r:2*r+1,:);
24.
25.%将imCum中的1+r到2*r+1行数据复制到imDst函数中的1到r+1行
26.
27.imDst(r+2:hei-r,:) = imCum(2*r+2:hei,:)-imCum(1:hei-2*r-1,:);
28.
29.%将imCum中的2*r+2到hei行数据依次减去1到hei-2*r-1行的值
30.
31.%把最终结果存放到函数imDst的r+2到hei-r行
32.
33.imDst(hei-r+1:hei,:) = repmat(imCum(hei,:),[r,1])-imCum(hei-2*r:hei-r-1,:);
34.
35.%使用repmat函数将imCum变为r*1的格式
36.
37.%然后减去函数imCum中从hei-2*r到hei-r-1行的数值
38.
39.%得到的结果放入imDst的hei-r+1到hei行位置
40.
41.%--------------------------------------------------------------------------
42.
43.imCum = cumsum(imDst,2); %x轴累计求和，将每一列的数值一次累加到下一列
44.
45.%y轴方向差异
46.
47.imDst(:, 1:r+1) = imCum(:,1+r:2*r+1);
48.
49.%将imCum中的1+r到2*r+1列数据复制到imDst函数中的1到r+1列
50.
51.imDst(:, r+2:wid-r) = imCum(:,2*r+2:wid)-imCum(:,1:wid-2*r-1);
52.
53.%将imCum中的2*r+2到hei列数据依次减去1到hei-2*r-1列的值
54.
55.%把最终结果存放到函数imDst的r+2到hei-r列
56.
57.imDst(:, wid-r+1:wid) = repmat(imCum(:,wid),[1,r])-imCum(:,wid-2*r:wid-r-1);
58.
59.%使用repmat函数将imCum变为r*1的格式
60.
61.%然后减去函数imCum中从hei-2*r到hei-r-1列的数值
62.
63.%得到的结果放入imDst的hei-r+1到hei列位置
64.
65.end