完整版MATLAB典型去雾算法代码.docx

1、完整版MATLAB典型去雾算法代码本节主要介绍基于Retinex理论的雾霭天气图像增强及其实现1.3.1 Rentinex 理论Retinex (视网膜“ Retina”和大脑皮层“ Cortex”的缩写)理论是一种建立 在科学实验和科学分析基础上的基于人类视觉系统 (Human Visual System)的图像增强理论。该算法的基本原理模型最早是由Edwin Land (埃德温？兰德)于1971 年提出的一种被称为的色彩的理论，并在颜色恒常性的基础上提出的一种图像增 强方法。Retinex理论的基本内容是物体的颜色是由物体对长波 (红)、中波(绿)和短波(蓝)光线的反射能力决定的，而不是由

2、反射光强度的绝对值决定的；物 体的色彩不受光照非均性的影响，具有一致性，即 Ret inex理论是以色感一致性(颜色恒常性)为基础的。根据Edwin Land提出的理论，一幅给定的图像S(x,y)分解成两幅不同的图像：图1.3-1 Retin ex理论示意图对于观察图像S中的每个点(x,y)，用公式可以表示为：S(x,y)=R(x,y) L(x,y) (1.3.1)实际上，Retinex理论就是通过图像S来得到物体的反射性质R，也就是去除了 入射光L的性质从而得到物体原本该有的样子。1.3.2 基于Retinex理论的图像增强的基本步骤 步骤一：利用取对数的方法将照射光分量和反射光分量分离，即

3、：Sz (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)=Sz (x, y)-log(D(x, y);步骤四：对G(x,y)取反对数，得到增强后的图像R(x, y):R(x, y)=exp(G(x, y);步骤五：对R(x,y)做对比度增强，得到最终的结果图像。133多尺

4、度Retinex算法D Jobson等人提出了多尺度Retinex算法，多尺度算法的基本公式是：NRi(x,y) Wn log Ii(x, y) log Fn(x,y) h(x,y)n 1其中，(x,y )是Ret inex的输出，i R,G,B表示3个颜色谱带，F (x, y)是 高斯滤波函数，Wn表示尺度的权重因子，N表示使用尺度的个数，N =3，表示 彩色图像，i R,G,B。N =1，表示灰度图像。从公式中可以看出： MSR算法的 特点是能产生包含色调再现和动态范围压缩这两个特性的输出图像。在MSR算法的增强过程中，图像可能会因为增加了噪声而造成对图像中 的局部区域色彩失真，使得物体的

5、真正颜色效果不能很好的显现出来， 从而影响了整体视觉效果。为了弥补这个缺点，一般情况下会应用带色彩恢复因子 C的多尺度算法(MSRCR)来解决。带色彩恢复因子C的多尺度算法(MSRCR)是在 多个固定尺度的基础上考虑色彩不失真恢复的结果，在多尺度 Ret inex算法过程中，我们通过引入一个色彩因子 C来弥补由于图像局部区域对比度增强而导致图像颜色失真的缺陷，通常情况下所引入的色彩恢复因子 C的表达式为RMSRCRj (x, y) Ci(x, y)RMSR(x, y) Ci (x, y) fli(x,y) f Ni(x,y) j! I j (x, y)其中，Ci表示第个通道的色彩恢复系数，它的

6、作用是用来调节 3个通道颜色的比例，f(?)表示的是颜色空间的映射函数。带色彩恢复的多尺度 Ret inex算法(MSRCR)通过色彩恢复因子C这个系数来调整原始图像中三个颜色通道之间 的比例关系，从而通过把相对有点暗的区域的信息凸显出来， 以达到消除图像色彩失真的缺陷。处理后的图像局域对比度提高，而且它的亮度与真实的场景很相 似，图像在人们视觉感知下显得极其逼真。因此， MSR算法具有较好的颜色再 现性、亮度恒常性以及动态范围压缩等特性。1.3.4例程精讲例程1.3.1是基于Retinex理论进行雾霭天气增强的 MATLAB程序，读者可 结合程序及注释对基于Ret inex理论进行雾霭天气增

7、强的基本原理进行进一步分 析，该程序的运行结果如图1.3-2所示。例程1.3.1clear;close all ;% 读入图像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);for i=1:N1for j=1:M1F(i,j)=exp(-(i-N1/2)A2+(j-M1/2)A2)/(2*sigma*sig

8、ma);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 =

9、 (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;for i=1:N1for j=1:M1 F(i,j)=exp(-(i-N1/2)A2+(j-M1/2)A2)/(2*sigma*sigma);endendF = F./(sum(F(:);% 对高斯滤波函数进行二维傅里叶变换Ffft=fft2

10、(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);%取输入图

11、像的B分量B=I(:,:,3);N1,M1=size(B);%对B分量进行数据转换，并对其取对数B0=double(B);Blog=log(B0+1);%对B分量进行二维傅里叶变换Bfft2=fft2(B0);% 形成高斯滤波函数sigma=250;for i=1:N1for j=1:M1 F(i,j)=exp(-(i-N1/2)A2+(j-M1/2)A2)/(2*sigma*sigma);endendF = F./(sum(F(:);% 对高斯滤波函数进行二维傅里叶变换Ffft=fft2(double(F);%对B分量与高斯滤波函数进行卷积运算DB0=Gfft2.*Ffft;DB=ifft2

12、(DB0);%在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像DBdouble=double(DB);DBIog=log(DBdouble+1);Bb=Blog-DBIog;EXPBb=exp(Bb);%对增强后的图像进行对比度拉伸增强MIN = min( mi n(EXPBb);MAX = max(max(EXPBb);EXPBb = (EXPBb-MIN)/(MAX-MIN);EXPBb=ad apthisteq(EXPBb);%对增强后的图像R、G、B分量进行融合IO(:,:,1)=EXPRr;IO(:,:,2)=EXPGg;IO(:,:,3)=EXPBb;%显示运行结果su

13、bplot(121),imshow(l);subplot(122),imshow(IO);*1.3-2 例程1.3.1的运行结果例程1.3.2是基于Retinex理论进行雾霭天气增强的 MATLAB程序，读者可 结合程序及注释对基于Ret inex理论进行雾霭天气增强的基本原理进行进一步分 析，该程序的运行结果如图1.3-3所示。例程1.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);

14、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:N1for j=1:M1 F(i,j)=exp(-(i-N1/2F2+(j-M1/2)A2)/(2*sigma*sigma);endendF = F./(sum(F(:);% 对高斯滤波函数进行二维傅里叶变换 Ffft=fft2(double(F);%对R分量与高斯滤波函数进行卷积运算DR0=Rfft2.*Ffft;

15、DR=ifft2(DR0);% 在对数域中，用原图像减去低通滤波后的图像，得到高频增强的图像 DRdouble=double(DR);DRlog=log(DRdouble+1);Rr0=Rlog-DRlog;% 形成高斯滤波函数( sigma=256 )sigma=256;F = zeros(N1,M1);for i=1:N1for j=1:M1 F(i,j)=exp(-(i-N1/2)A2+(j-M1/2)A2)/(2*sigma*sigma);endendF = F./(sum(F(:);% 对高斯滤波函数进行二维傅里叶变换 Ffft=fft2(double(F);%对R分量与高斯滤波函数

16、进行卷积运算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:N1for j=1:M1F(i,j)=exp(-(i-N1/2F2+(j-M1/2)A2)/(2*sigma*sigma);endendF = F./(sum(F(:);% 对高斯滤波函数进行二维傅里叶变换Ffft=fft2(do

17、uble(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);% 将增强后

18、的 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);R分量处理的步骤进行理解%对G分量进行处理，步骤与对 R分量处理的步骤相同，请读者仿照G0=double(G);Glog=log(G0+1);Gfft2=fft2(G0);sigma=128;F = zeros(N1,M1

19、);for i=1:N1for j=1:M1F(i,j)=exp(-(i-N1/2F2+(j-M1/2)A2)/(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);for i=1:N1for j=1:M1F(i,j)=exp(-(i-N1/2)A2+(j-M1/2)A2)/(2*sigma*sigma);

20、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);for i=1:N1for j=1:M1F(i,j)=exp(-(i-N1/2F2+(j-M1/2)A2)/(2*sigma*sigma);endendF = F./(sum(F(:);Ffft=fft2(double(F);DG0=Gfft2.*Ffft;DG=ifft2(DG0);

21、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);R分

22、量处理的步骤进行理解%对B分量进行处理，步骤与对 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:M1F(i,j)=exp(-(i-N1/2F2+(j-M1/2)A2)/(2*sigma*sigma);endendF = F./(sum(F(:);Ffft=fft2(double(F);DB0=Bfft2.*Ffft;DB=ifft2(DB0);DBdouble=double(DB);DBlog=log(

23、DBdouble+1);Bb0=Blog-DBlog;sigma=256;F = zeros(N1,M1);for i=1:N1 for j=1:M1F(i,j)=exp(-(i-N1/2)A2+(j-M1/2)A2)/(2*sigma*sigma);endendF = 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);for i=1:N1for j=1

24、:M1F(i,j)=exp(-(i-N1/2F2+(j-M1/2)A2)/(2*sigma*sigma);endendF = 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(l); subplot(122),imshow(IO);*1.3-3例程1.3.2的运行结果