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

本文主要介绍基于Retinex理论的雾霭天气图像增加及其实现。并通过编写两个程序来实现图像的去雾功能。1Rentinex理论Retinex〔视网膜“Retina〞和大脑皮层“Cortex〞的缩写〕理论是一种建立在科学试验和科学分析根底上的基于人类视觉系统〔HumanVisualSystem〕的图像增加理论。该算法的根本原理模型最早是由EdwinLand〔埃德温•兰德〕于1971年提出的一种被称为的色调的理论，并在颜色恒常性的根底上提出的一种图像增加方法。Retinex理论的根本内容是物体的颜色是由物体对长波〔红〕、中波〔绿〕和短波〔蓝〕光线的反射实力确定的，而不是由反射光强度的肯定值确定的；物体的色调不受光照非均性的影响，具有一样性，即Retinex理论是以色感一样性〔颜色恒常性〕为根底的。依据EdwinLand提出的理论，一幅给定的图像S(x,y)分解成两幅不同的图像：反射物体图像R(x,y)和入射光图像L(x,y)，其原理示意图如图8.3-1所示。图-1Retinex理论示意图对于视察图像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算法DJobson等人提出了多尺度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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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);****************************************************************************************图-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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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);%定义色调复原因子Ca=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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=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:N1forj=1:M1F(i,j)=exp(-((i-N1/2)^2+(j-M1/2)^2)/(2*sigma*sigma));endendF=F./(sum(F(:)));Ffft=fft2(double(F));DB0=Bfft2.*Ffft;DB=ifft