基于matlab的圖像去霧算法詳細講解與實現(xiàn)-附matlab實現(xiàn)源代碼

本文主要介紹基于Retinex理論的霧靄天氣圖像增加及其實現(xiàn)。并通過編寫兩個程序來實現(xiàn)圖像的去霧功能。1Rentinex理論Retinex〔視網(wǎng)膜“Retina〞和大腦皮層“Cortex〞的縮寫〕理論是一種建立在科學試驗和科學分析根底上的基于人類視覺系統(tǒng)〔HumanVisualSystem〕的圖像增加理論。該算法的根本原理模型最早是由EdwinLand〔埃德溫?蘭德〕于1971年提出的一種被稱為的色調(diào)的理論，并在顏色恒常性的根底上提出的一種圖像增加方法。Retinex理論的根本內(nèi)容是物體的顏色是由物體對長波〔紅〕、中波〔綠〕和短波〔藍〕光線的反射實力確定的，而不是由反射光強度的肯定值確定的；物體的色調(diào)不受光照非均性的影響，具有一樣性，即Retinex理論是以色感一樣性〔顏色恒常性〕為根底的。依據(jù)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來得到物體的反射性質(zhì)R，也就是去除了入射光L的性質(zhì)從而得到物體原本該有的樣子。2基于Retinex理論的圖像增加的根本步驟步驟一:利用取對數(shù)的方法將照耀光重量和反射光重量別離，即：S＇(x,y)=r(x,y)+l(x,y)=log(R(x,y))+log(L(x,y))；步驟二：用高斯模板對原圖像做卷積，即相當于對原圖像做低通濾波，得到低通濾波后的圖像D(x,y)，F(xiàn)〔x,y〕表示高斯濾波函數(shù)：D(x,y)=S(x,y)*F〔x,y〕；步驟三：在對數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增加的圖像G〔x,y〕：

G(x,y)=S＇(x,y)-log(D(x,y))

；步驟四：對G(x,y)取反對數(shù)，得到增加后的圖像R(x,y)：

R(x,y)=exp(G(x,y))；步驟五：對R(x,y)做比照度增加，得到最終的結(jié)果圖像。3多尺度Retinex算法DJobson等人提出了多尺度Retinex算法，多尺度算法的根本公式是：其中，是Retinex的輸出，表示3個顏色譜帶，是高斯濾波函數(shù)，表示尺度的權(quán)重因子，表示運用尺度的個數(shù)，=3，表示彩色圖像，。=1，表示灰度圖像。從公式中可以看出：MSR算法的特點是能產(chǎn)生包含色調(diào)再現(xiàn)和動態(tài)范圍壓縮這兩個特性的輸出圖像。在MSR算法的增加過程中，圖像可能會因為增加了噪聲而造成對圖像中的局部區(qū)域色調(diào)失真，使得物體的真正顏色效果不能很好的顯現(xiàn)出來，從而影響了整體視覺效果。為了彌補這個缺點，一般狀況下會應用帶色調(diào)復原因子C的多尺度算法〔MSRCR〕來解決。帶色調(diào)復原因子C的多尺度算法(MSRCR)]是在多個固定尺度的根底上考慮色調(diào)不失真復原的結(jié)果，在多尺度Retinex算法過程中，我們通過引入一個色調(diào)因子C來彌補由于圖像局部區(qū)域比照度增加而導致圖像顏色失真的缺陷，通常狀況下所引入的色調(diào)復原因子C的表達式為其中,表示第個通道的色調(diào)復原系數(shù)，它的作用是用來調(diào)整3個通道顏色的比例，表示的是顏色空間的映射函數(shù)。帶色調(diào)復原的多尺度Retinex算法〔MSRCR〕通過色調(diào)復原因子C這個系數(shù)來調(diào)整原始圖像中三個顏色通道之間的比例關(guān)系，從而通過把相對有點暗的區(qū)域的信息凸顯出來，以到達消退圖像色調(diào)失真的缺陷。處理后的圖像局域比照度提高，而且它的亮度與真實的場景很相像，圖像在人們視覺感知下顯得極其逼真。因此，MSR算法具有較好的顏色再現(xiàn)性、亮度恒常性以及動態(tài)范圍壓縮等特性。4例程精講例程1是基于Retinex理論進展霧靄天氣增加的MATLAB程序，讀者可結(jié)合程序及注釋對基于Retinex理論進展霧靄天氣增加的根本原理進展進一步分析，該程序的運行結(jié)果如圖-2所示。例程1：clear;closeall;%讀入圖像I=imread('wu.png');%取輸入圖像的R重量R=I(:,:,1);[N1,M1]=size(R);%對R重量進展數(shù)據(jù)轉(zhuǎn)換,并對其取對數(shù)R0=double(R);Rlog=log(R0+1);%對R重量進展二維傅里葉變換Rfft2=fft2(R0);%形成高斯濾波函數(shù)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(:)));%對高斯濾波函數(shù)進展二維傅里葉變換Ffft=fft2(double(F));%對R重量與高斯濾波函數(shù)進展卷積運算DR0=Rfft2.*Ffft;DR=ifft2(DR0);%在對數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增加的圖像DRdouble=double(DR);DRlog=log(DRdouble+1);Rr=Rlog-DRlog;%取反對數(shù)，得到增加后的圖像重量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重量進展數(shù)據(jù)轉(zhuǎn)換,并對其取對數(shù)G0=double(G);Glog=log(G0+1);%對G重量進展二維傅里葉變換Gfft2=fft2(G0);%形成高斯濾波函數(shù)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(:)));%對高斯濾波函數(shù)進展二維傅里葉變換Ffft=fft2(double(F));%對G重量與高斯濾波函數(shù)進展卷積運算DG0=Gfft2.*Ffft;DG=ifft2(DG0);%在對數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增加的圖像DGdouble=double(DG);DGlog=log(DGdouble+1);Gg=Glog-DGlog;%取反對數(shù)，得到增加后的圖像重量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重量進展數(shù)據(jù)轉(zhuǎn)換,并對其取對數(shù)B0=double(B);Blog=log(B0+1);%對B重量進展二維傅里葉變換Bfft2=fft2(B0);%形成高斯濾波函數(shù)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(:)));%對高斯濾波函數(shù)進展二維傅里葉變換Ffft=fft2(double(F));%對B重量與高斯濾波函數(shù)進展卷積運算DB0=Gfft2.*Ffft;DB=ifft2(DB0);%在對數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增加的圖像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;%顯示運行結(jié)果subplot(121),imshow(I);subplot(122),imshow(I0);****************************************************************************************圖-2例程1的運行結(jié)果例程2是基于Retinex理論進展霧靄天氣增加的MATLAB程序，讀者可結(jié)合程序及注釋對基于Retinex理論進展霧靄天氣增加的根本原理進展進一步分析，該程序的運行結(jié)果如圖-3所示。例程2：****************************************************************************************clear;closeall;I=imread('wu.png');%分別取輸入圖像的R、G、B三個重量，并將其轉(zhuǎn)換為雙精度型R=I(:,:,1);G=I(:,:,2);B=I(:,:,3);R0=double(R);G0=double(G);B0=double(B);[N1,M1]=size(R);%對R重量進展對數(shù)變換Rlog=log(R0+1);%對R重量進展二維傅里葉變換Rfft2=fft2(R0);%形成高斯濾波函數(shù)〔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(:)));%對高斯濾波函數(shù)進展二維傅里葉變換Ffft=fft2(double(F));%對R重量與高斯濾波函數(shù)進展卷積運算DR0=Rfft2.*Ffft;DR=ifft2(DR0);%在對數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增加的圖像DRdouble=double(DR);DRlog=log(DRdouble+1);Rr0=Rlog-DRlog;%形成高斯濾波函數(shù)〔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(:)));%對高斯濾波函數(shù)進展二維傅里葉變換Ffft=fft2(double(F));%對R重量與高斯濾波函數(shù)進展卷積運算DR0=Rfft2.*Ffft;DR=ifft2(DR0);%在對數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增加的圖像DRdouble=double(DR);DRlog=log(DRdouble+1);Rr1=Rlog-DRlog;%形成高斯濾波函數(shù)〔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(:)));%對高斯濾波函數(shù)進展二維傅里葉變換Ffft=fft2(double(F));%對R重量與高斯濾波函數(shù)進展卷積運算DR0=Rfft2.*Ffft;DR=ifft2(DR0);%在對數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增加的圖像DRdouble=double(DR);DRlog=log(DRdouble+1);Rr2=Rlog-DRlog;%對上述三次增加得到的圖像取均值作為最終增加的圖像Rr=(1/3)*(Rr0+Rr1+Rr2);%定義色調(diào)復原因子Ca=125;II=imadd(R0,G0);II=imadd(II,B0);Ir=immultiply(R0,a);C=imdivide(Ir,II);C=log(C+1);%將增加后的R重量乘以色調(diào)復原因子，并對其進展反對數(shù)變換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