【老生谈算法】matlab实现图像去雾算法——图像去雾

基于matlab的图像去雾算法详细讲解与实现-附matlab实现源代码

本文主要介绍基于Retinex理论的雾霭天气图像增强及其实现。并通过编写两个程序来实现图像的去雾功能。

1 Rentinex理论
Retinex（视网膜“Retina”和大脑皮层“Cortex”的缩写）理论是一种建立在科学实验和科学分析基础上的基于人类视觉系统（Human Visual System）的图像增强理论。该算法的基本原理模型最早是由Edwin Land（埃德温•兰德）于1971年提出的一种被称为的色彩的理论，并在颜色恒常性的基础上提出的一种图像增强方法。Retinex 理论的基本内容是物体的颜色是由物体对长波（红）、中波（绿）和短波（蓝）光线的反射能力决定的，而不是由反射光强度的绝对值决定的；物体的色彩不受光照非均性的影响，具有一致性，即Retinex理论是以色感一致性（颜色恒常性）为基础的。

根据Edwin Land提出的理论，一幅给定的图像S(x,y)分解成两幅不同的图像：反射物体图像R(x,y)和入射光图像L(x,y)，其原理示意图如图8.3-1所示。

图-1 Retinex理论示意图

对于观察图像S中的每个点(x,y)，用公式可以表示为：
S(x,y)=R(x,y)×L(x,y) （1.3.1）
实际上，Retinex理论就是通过图像S来得到物体的反射性质R，也就是去除了入射光L的性质从而得到物体原本该有的样子。

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)做对比度增强，得到最终的结果图像。

3 多尺度Retinex算法
D Jobson等人提出了多尺度Retinex算法，多尺度算法的基本公式是：
其中，是Retinex的输出，表示3个颜色谱带，是高斯滤波函数，表示尺度的权重因子，表示使用尺度的个数，=3，表示彩色图像，。=1，表示灰度图像。从公式中可以看出：MSR算法的特点是能产生包含色调再现和动态范围压缩这两个特性的输出图像。

在MSR算法的增强过程中，图像可能会因为增加了噪声而造成对图像中的局部区域色彩失真，使得物体的真正颜色效果不能很好的显现出来，从而影响了整体视觉效果。为了弥补这个缺点，一般情况下会应用带色彩恢复因子C的多尺度算法（MSRCR）来解决。带色彩恢复因子C的多尺度算法(MSRCR)]是在多个固定尺度的基础上考虑色彩不失真恢复的结果，在多尺度Retinex算法过程中，我们通过引入一个色彩因子C来弥补由于图像局部区域对比度增强而导致图像颜色失真的缺陷，通常情况下所引入的色彩恢复因子C的表达式为

其中,表示第个通道的色彩恢复系数，它的作用是用来调节3个通道颜色的比例，表示的是颜色空间的映射函数。带色彩恢复的多尺度Retinex算法（MSRCR）通过色彩恢复因子C这个系数来调整原始图像中三个颜色通道之间的比例关系，从而通过把相对有点暗的区域的信息凸显出来，以达到消除图像色彩失真的缺陷。处理后的图像局域对比度提高，而且它的亮度与真实的场景很相似，图像在人们视觉感知下显得极其逼真。因此，MSR算法具有较好的颜色再现性、亮度恒常性以及动态范围压缩等特性。

4 例程精讲
例程1是基于Retinex理论进行雾霭天气增强的MATLAB程序，读者可结合程序及注释对基于Retinex理论进行雾霭天气增强的基本原理进行进一步分析，该程序的运行结果如图-2所示。

例程1：

clear;
closeall;
% 读入图像
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);
fori=1:N1
forj=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;
fori=1:N1
forj=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;
fori=1:N1
forj=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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

图-2例程1的运行结果

例程2是基于Retinex理论进行雾霭天气增强的MATLAB程序，读者可结合程序及注释对基于Retinex理论进行雾霭天气增强的基本原理进行进一步分析，该程序的运行结果如图-3所示。

例程2：
****************************************************************************************
clear;
closeall;
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);
fori=1:N1
forj=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);
fori=1:N1
forj=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);
fori=1:N1
forj=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);
fori=1:N1
forj=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);
fori=1:N1
forj=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);
fori=1:N1
forj=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);
fori=1:N1
forj=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);
fori=1:N1
forj=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);
fori=1:N1
forj=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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265

图-3例程2的运行结果