基于MATLAB retinex的图像去雾系统

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

1　Rentinex理论

Retinex(视网膜“Ｒｅｔina”和大脑皮层“Cortex”的缩写）理论是一种建立在科学实验和科学分析基础上的基于人类视觉系统（Hｕman Viｓual　System）的图像增强理论。该算法的基本原理模型最早是由Ｅdｗin Ｌand(埃德温•兰德）于1９7１年提出的一种被称为的色彩的理论，并在颜色恒常性的基础上提出的一种图像增强方法。Ｒeｔｉnex 理论的基本内容是物体的颜色是由物体对长波(红)、中波（绿）和短波（蓝)光线的反射能力决定的，而不是由反射光强度的绝对值决定的;物体的色彩不受光照非均性的影响，具有一致性，即Retｉｎex理论是以色感一致性（颜色恒常性)为基础的。

根据Eｄwin Lanｄ提出的理论,一幅给定的图像S(x,y)分解成两幅不同的图像:反射物体图像R(x,y)和入射光图像L(x,y),其原理示意图如图8.3－1所示。

图-1 Retinｅx理论示意图

对于观察图像S中的每个点(ｘ,y）,用公式可以表示为:　ﻫ 　 　 　 　 　 S(ｘ,y）＝R(x，y)×Ｌ(x，y) 　 　　　　 　　 　

（1.3．1）

实际上，Retinｅx理论就是通过图像S来得到物体的反射性质R，也就是去除了入射光Ｌ的性质从而得到物体原本该有的样子。

2 基于Reｔｉｎｅｘ理论的图像增强的基本步骤

步骤一: 利用取对数的方法将照射光分量和反射光分量分离，即:

S'(x， ｙ)=r(ｘ， y)+ｌ(ｘ, ｙ)＝ｌog(Ｒ(ｘ， y))+lｏｇ(Ｌ(x, ｙ))；

步骤二:用高斯模板对原图像做卷积，即相当于对原图像做低通滤波,得到低通滤波后的图像D（x,y）,F(x, y）表示高斯滤波函数：

Ｄ（ｘ， y)＝Ｓ(x, y) *F（x, y)；

步骤三:在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像G（x,　ｙ）:

G(ｘ,y)=S＇(ｘ, y)－log(D(ｘ, ｙ)) ;

步骤四:对G（ｘ,y)取反对数，得到增强后的图像R(x， y):

Ｒ（ｘ， y）=eｘp（G(ｘ, y))；

步骤五：对R（ｘ，y)做对比度增强,得到最终的结果图像。

３ 多尺度Retinｅx算法

4 例程精讲

例程1是基于Retｉnｅx理论进行雾霭天气增强的MATＬAB程序,读者可结合程序及注释对基于Retｉnex理论进行雾霭天气增强的基本原理进行进一步分析，该程序的运行结果如图-2所示。

例程1：

cleaｒ；
clｏsｅ　alｌ;
％ 读入图像
I=iｍｒeaｄ('wu.ｐng')；
% 取输入图像的R分量
Ｒ=I(:,：,1);
[N1，M１]=sｉze（Ｒ);
%　对R分量进行数据转换,并对其取对数
R０=douｂle(R）;
Rlog=ｌog(R0+１)；
% 对R分量进行二维傅里叶变换
Rfｆt2=ffｔ2（R０）;
% 形成高斯滤波函数
sigma=２５０;
F = zｅros(N1,M１);
fｏｒ ｉ=1:N1
       for　j=1:M1
   　 　　 Ｆ(i，ｊ）=exp(-((i-N1/2)^2+(j-M1／２）^2)/(2*ｓigma*ｓigmａ));
　  　   eｎd
end
F = Ｆ./(sum（F(:)))；
%　对高斯滤波函数进行二维傅里叶变换
Fffｔ=fｆt2(douｂle(Ｆ)）;
% 对Ｒ分量与高斯滤波函数进行卷积运算
ＤR0=Rfft２.*Fｆｆt;
DR=iｆｆｔ２(DＲ0）;
% 在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像
DRdｏｕｂle=douｂｌｅ（DR)；
DRｌｏｇ=lｏg（DRdouble+1);
Rr=Rlog-ＤRｌｏg；
% 取反对数，得到增强后的图像分量
EXＰRｒ=exp(Rｒ);
％ 对增强后的图像进行对比度拉伸增强
ＭＩN =　mｉn（min(EXＰRr）)；
MAX = ｍａｘ（maｘ(EXPRr)）；
EXＰRｒ ＝ (ＥＸＰRr-MIN)/(MAＸ-MIN)；
EXＰRr=ａdaｐtｈisｔｅq(EXＰRr)；

% 取输入图像的G分量
Ｇ=I(：,:,2);
［N1,M1]=size(G);
%　对G分量进行数据转换,并对其取对数
G0＝dｏuｂlｅ(Ｇ）;
Glｏg＝log(G0+1）;
% 对G分量进行二维傅里叶变换
Gfft2=ｆft2(G0）;
％　形成高斯滤波函数
siｇｍa=250；
ｆoｒ i=1:Ｎ1
   　  　ｆor j=1：M1
    　 　 F(ｉ，ｊ)=exp(-(（ｉ-N1/2）＾2+(j-M1/2)^2)/(2*sigma*sigmａ）);
  　 　　 end
eｎd
F ＝ Ｆ./(sum(Ｆ（:)));
％ 对高斯滤波函数进行二维傅里叶变换
Fｆｆt=ffｔ２(doublｅ(F));
％　对G分量与高斯滤波函数进行卷积运算
DＧ０=Ｇfｆt2.*Ffｆt；
DＧ=iｆfｔ2(DG0）;
％ 在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像
DGdｏuｂle＝douｂle（ＤG）；
DGｌｏg=loｇ(DGdｏuble+1);
Gg=Ｇlｏｇ-DGloｇ;
% 取反对数,得到增强后的图像分量
EXPGg=eｘp（Gｇ）;
%　对增强后的图像进行对比度拉伸增强
MIN =　ｍin(miｎ（ＥXPＧg));
ＭAＸ ＝　mａｘ(max(EXPGg));
EXPＧg　= (ＥXPGg-MＩN)/（MＡX-ＭＩＮ);
EXPGｇ=ａdapｔｈisteｑ(EXPGg);

%　取输入图像的B分量
B=I(:,:,3）;
[N1,Ｍ1]=size(Ｂ);
％　对Ｂ分量进行数据转换,并对其取对数
B０=dｏｕｂle(Ｂ）；
Bloｇ=ｌog(Ｂ0+1);
％　对Ｂ分量进行二维傅里叶变换
Bffｔ２＝fft2（B0）;
％ 形成高斯滤波函数
sigma=25０;
for　i＝1:N１
 　    　for j=１:M1
  　    　F（ｉ,ｊ)=exｐ(-((i－N１/2)＾2+(j-M1/2）＾2)/(2*siｇｍa*ｓiｇma));
       end
end
F =　Ｆ./（sum(F(：)）);
% 对高斯滤波函数进行二维傅里叶变换
Ffft=ｆft2(doubｌe(Ｆ）);
％　对B分量与高斯滤波函数进行卷积运算
DＢ０＝Ｇfft２.*Ffft;
DB＝ifft2(DB０);
%　在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像
DBdｏｕble=doｕble(ＤB）;
DＢlｏg＝ｌoｇ(DBｄouｂｌe+1);
Bb＝Ｂlog-ＤBｌog；
EXPＢｂ=exp(Bb);
%　对增强后的图像进行对比度拉伸增强
MIＮ　= mｉn（ｍin（EXPＢb))；
MAＸ　= max(maｘ(ＥXＰBｂ));
EXＰＢb = (EXＰBb－ＭIN)/(MAX-ＭIＮ);
ＥＸＰＢb=adaｐthｉsteq(EＸＰBｂ）;
% 对增强后的图像R、G、Ｂ分量进行融合
I0（:,:，1)=EＸPＲr；
Ｉ0(:,：,２)＝EXPGｇ；
I0(:，:,3)=EXＰBb;
%　显示运行结果
suｂplot(12１）,imｓhｏw(Ｉ)；
ｓubｐlｏｔ（12２)，ｉmshｏw（I０);
*****＊***＊****＊**＊*＊***＊****＊**＊*****＊*＊**＊＊＊＊＊*****＊＊*****＊**＊*******＊************＊***＊

图-2 例程１的运行结果

例程2是基于Reｔineｘ理论进行雾霭天气增强的ＭATＬAB程序,读者可结合程序及注释对基于Reｔinｅx理论进行雾霭天气增强的基本原理进行进一步分析,该程序的运行结果如图－3所示。

例程2:

****＊********＊＊******＊**＊********＊*＊************＊＊＊＊*＊*****＊******＊*＊***＊*＊*****＊**＊*＊*＊
ｃlear;
cloｓｅ aｌl；
Ｉ=imrｅad（'wｕ.png＇);
％ 分别取输入图像的R、Ｇ、B三个分量，并将其转换为双精度型
R=I(:,:，1)；
G=Ｉ（:，:,2);
B=I(:,:,３);
Ｒ0=doublｅ(Ｒ);
Ｇ０＝doubｌe(G）；
Ｂ0=ｄoublｅ（B）；
 
[N1,Ｍ1]＝sizｅ（R);

% 对R分量进行对数变换
Rlog=log（R0+１);
％ 对R分量进行二维傅里叶变换
Rfft2=fft2(R０）;
% 形成高斯滤波函数(sigmａ=１２8)
sigma=128;
F =　zerｏｓ(N1,M1）;
foｒ i=1：Ｎ1
  　   　for j＝1:M１
　  　    F（ｉ,j)＝exp(－(（i－N1／2)^2+(j-M1/2)^2）/(2*sigｍa*siｇｍa));
  　 　  end
ｅnd
F = F./（sｕｍ(Ｆ（:）)）;
% 对高斯滤波函数进行二维傅里叶变换
Ｆfft＝ｆｆt2(doｕbｌe(Ｆ));
％ 对Ｒ分量与高斯滤波函数进行卷积运算
ＤＲ0=Rfft2.*Fffｔ;
DR=ifft2(ＤR0)；
％　在对数域中，用原图像减去低通滤波后的图像,得到高频增强的图像
DＲｄouble=ｄouｂｌe(DR);
DRloｇ=log(ＤRdｏuble＋1)；
Ｒr0=Rlｏｇ-DＲｌoｇ;

％　形成高斯滤波函数（siｇma=256）
ｓigｍa=２56;
F = zerｏｓ(N１,M１);
for ｉ＝1:Ｎ1
  　 　　 for j=１:Ｍ1
　       F(i，j）=exp(-（(i-N1/2)^2＋(ｊ－Ｍ１/2)^2)/（2*sigｍa*sｉgmａ)）;
  　    end
ｅnd
F = F./(sum（F(:）));
% 对高斯滤波函数进行二维傅里叶变换
Ｆffｔ=ｆft2(dｏuble（F));
% 对Ｒ分量与高斯滤波函数进行卷积运算
DR0＝Rffｔ2.＊Fｆfｔ；
DR＝iｆft２(DR0)；
％ 在对数域中,用原图像减去低通滤波后的图像，得到高频增强的图像
DＲｄouble=doubｌe(DR);
DRlog=log(DRdｏublｅ+1);
Rr１=Rlog-DRｌｏｇ；

% 形成高斯滤波函数(sigma＝５1２）
sigmａ＝512;
 F ＝ zｅros（N１,Ｍ1);
for i＝1:Ｎ１
　　     fｏr　j=1:M1
 　     　F（i,j）=exp(-（(i－N１／２)^2+（j－M1／２)^2)/(2＊sigma*sｉgma));
  　 　　 end
end
Ｆ　=　F./（sum(Ｆ(:)）);
% 对高斯滤波函数进行二维傅里叶变换
Ffft=ffｔ2(ｄｏuble(F));
% 对R分量与高斯滤波函数进行卷积运算
DR0=Rffｔ2.＊Ｆffｔ;
ＤR=ifft2(ＤR0);
％ 在对数域中,用原图像减去低通滤波后的图像,得到高频增强的图像
DRdouble=doｕblｅ(ＤR);
DRｌog=log(DRdouble＋1);
Rｒ2=Rlｏｇ-DRlog;

% 对上述三次增强得到的图像取均值作为最终增强的图像
Ｒr=(1/３)*（Ｒr0＋Rr1+Rr２）;
 
％ 定义色彩恢复因子C
a＝125;
II＝imａdd(R0，G0);
IＩ=iｍadd（II,B0);
Iｒ=iｍｍuｌtｉplｙ(R０，a);
C=iｍdivide(Ir,II）;
Ｃ=loｇ(C+1）；
 
%　将增强后的R分量乘以色彩恢复因子,并对其进行反对数变换
Rr=immｕltiplｙ(C，Ｒｒ）；
ＥＸＰＲr=exp(Rｒ）;

％　对增强后的R分量进行灰度拉伸 
ＭIN =　min（min(EXPRr));
MAX　＝ maｘ(max(EXPＲr));
ＥXPＲr =　(EXPＲｒ-MIＮ)／(MAX-MIN)；
ＥXＰRr＝ａdapthisteｑ（EXPＲｒ)；
 
[Ｎ1,M1］=ｓiｚe（G）；
% 对G分量进行处理，步骤与对Ｒ分量处理的步骤相同,请读者仿照R分量处理的步骤进行理解。
Ｇ0=ｄoｕblｅ(G)；
Glog＝log(G0+1）；
 

Gffｔ2=ｆft2(G0）;
ｓigma=12８;
Ｆ ＝ zerｏｓ(Ｎ1,M1）;
fｏr i＝1:N1
   　  　for ｊ=1:Ｍ1
　 　  　　　F(ｉ，j)=exp(－((ｉ-N1/2)^２+(j-M1/２)^2)/(2*sigmａ*sｉgma));
       ｅnd
enｄ
F = Ｆ.／(ｓuｍ（F(：）));
 
Fffｔ=ｆfｔ2(ｄoublｅ(F));
 
DG0=Gfft2.＊Ｆfft;
ＤG＝ifft２（DＧ0);
 
DGdｏuｂlｅ=double（DG);
DGlｏg=log(DGdouble+1);
Gg0=Glog-DGloｇ；
 
sｉgma=２56;
 
Ｆ ＝ ｚeros（N1，Ｍ1);
for i=１:N1
   　　  for　j=1:M1
     　  F(i,j）=eｘp(-((ｉ-N1/2)^２+(j-M1/2）^2)／(2＊sｉgma＊ｓｉgmａ）);
 　   　　ｅnd
eｎd
Ｆ = F./(sｕm(F（:)）);
 
Ffft=fｆt2(dｏｕblｅ（F));
 
DG0＝Gｆft2.*Ffft;
DＧ＝ｉfｆt2(DＧ0);
 
DＧｄouｂle=double(ＤG);
DGlog=lｏｇ(DGdｏuｂle＋1)；
Gg１＝Gloｇ－DGｌog;
 
 
sigma=512；
 
F　= zerｏs(N１,Ｍ1);
foｒ i＝1:N1
　     　for　j=１:M1
   　    F(i,j)＝exp（-(（i-N1/２）^2+（j－M1/２）＾２）/（2*sigmａ＊sigma));
 　 　　 　ｅnd
ｅnd
F　=　F．/(ｓum(F（：）))；
 
Ffft=ffｔ２(doｕｂlｅ（F));
 
DＧ0=Gfft2．＊Ffft；
DＧ＝iffｔ2(DＧ0);
 
DGdｏuble=douｂle(DＧ)；
DＧloｇ=loｇ(DGｄoubｌe＋１);
Ｇg2=Ｇlog-DGlog；
 
Gg＝（1/3）*（Ｇg0+Gg１+Gg2）;
 
 
a=125;
IＩ=imadd（R0，G0）;
ＩI＝imadd（II，B0)；
Ｉr=ｉｍmuｌtiply(R0,a）;
 
Ｇg=immｕlｔiply(Ｃ,Ｇg)；
 
 
EXPGg=ｅxp(Ｇg)；
 
MIN ＝　min（min（EXＰGg));
MAX　= max(max(EXPGｇ));
EXＰGg = (ＥXPGg-MＩN)／(MＡX－MＩＮ);
EXPGｇ=ａdapthisｔeq(ＥXPGｇ)；
 
% 对B分量进行处理,步骤与对R分量处理的步骤相同,请读者仿照R分量处理的步骤进行理解。
 [N1,M1]=sizｅ(Ｂ);
 
B0=doublｅ(Ｂ)；
Ｂloｇ=log（B0+1);
 
Bｆft2＝fft２(B0);
sｉgma＝１28;
 
F ＝ zeros（Ｎ1,M１)；
fｏr i=１:N1
 　     ｆor ｊ＝１：M１
       　F(ｉ,ｊ)=exp（-((ｉ-N1/2)^2＋(j-M1/２）＾2)/(2*sｉgmａ*sｉｇｍa));
  　 　 　eｎd
end
F　=　F．/(suｍ(F(:)));
 
Ｆfft=fｆｔ2(doubｌe(F));
 
DB０=Bfft2.*Ffｆｔ;
DB=ｉfft2(DＢ０);
 
DＢdoubｌe=ｄouｂle(DB)；
DBｌoｇ=lｏg(DBｄouble+1)；
Bb0=Blog-DBloｇ;
 
sigｍa=256;
 
F =　zerｏs（N１，M1）;
for i=1:N1
    　　 ｆor j=1：M１
 　  　　　 Ｆ（i,ｊ)=eｘp(-(（i-N１／2）＾2+(ｊ-M１/2)＾2)/（２*sigma*ｓigｍａ));
       end
eｎd
F ＝ Ｆ./(ｓuｍ(Ｆ(:)))；
 
Ffft=ffｔ2(dｏuble（F));
 
ＤB0=Ｂfft2．＊Ｆffｔ;
ＤB=ｉfft2(DB0);
 
DBdouble＝doubｌe（ＤB）；
DＢｌog=log(ＤBdouble＋1);
Bb1=Bｌog－DBｌog;
 
 
sｉgma＝5１２;
 
F = zｅｒoｓ(N１,M1）;
fｏr i=１:N1
 　　    for　j=1：M１
 　  　  　Ｆ(ｉ，j)=exp（-((ｉ-N1/2)^2+(j-M1／２)＾2)/(2*ｓｉgｍa*sigｍa)）;
    　  end
ｅnd
F = F.／(sum(Ｆ（：））);
 
Ｆfｆt＝ffｔ2(double(F));
 
ＤＢ０=Rfft2.*Fｆｆt;
DB=iffｔ2(ＤB０);
 
 
Bｂ＝immｕｌtｉpｌy(C,Bb)；
EXＰBb＝exp(Ｂb);
 
ＭＩN ＝ min(ｍin(EXPBｂ))；
MAX = max(maｘ(ＥXPBb））；
EＸPBｂ　= (EXＰBb－MIＮ）/(MAＸ-ＭIN)；
EXPＢb=adaｐtｈistｅq(EXPBb）;

%　对增强后的图像Ｒ、G、B分量进行融合
I０(:，:,１）=EXPRｒ；
Ｉ０(:,:，2)=EXPGｇ;
I0(：,：,３)=EXPBb;
 
％　显示运行结果
subpｌot（121),imshow（I)；
subｐlｏｔ（１2２)，imｓhow(I0）;
***＊***＊*＊***＊**＊*＊****＊*********＊＊＊****＊*****＊＊＊*＊＊********＊********＊**＊***************

图-3　例程２的运行结果