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算法

　　D Jobsｏn等人提出了多尺度Ｒetｉneｘ算法,多尺度算法的基本公式是：

　 其中,

是Rｅtineｘ的输出，

表示3个颜色谱带，

是高斯滤波函数,

表示尺度的权重因子,

表示使用尺度的个数，

＝3，表示彩色图像,

。

＝1，表示灰度图像。从公式中可以看出：ＭSR算法的特点是能产生包含色调再现和动态范围压缩这两个特性的输出图像。

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

其中，

表示第个通道的色彩恢复系数，它的作用是用来调节3个通道颜色的比例,

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

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　例程２的运行结果