例程2是基于Retinex理论进行雾霭天气增强的MATLAB程序，读者可结合程序及注释对基于Retinex理论进行雾霭天气增强的基本原理进行进一步分析，该程序的运行结果如图-3所示。 例程2：
**************************************************************************************** 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(:)));