1、本文主要介紹基于Retinex理論的霧靄天氣圖像增強及其實現。并通過編寫兩個程序來實現圖像的去霧功能。1 Rentinex理論Retinex（視網膜“Retina”和大腦皮層“Cortex”的縮寫）理論是一種建立在科學實驗和科學分析基礎上的基于人類視覺系統（Human Visual System）的圖像增強理論。該算法的基本原理模型最早是由Edwin Land（埃德溫蘭德）于1971年提出的一種被稱為的色彩的理論，并在顏色恒常性的基礎上提出的一種圖像增強方法。Retinex 理論的基本內容是物體的顏色是由物體對長波（紅）、中波（綠）和短波（藍）光線的反射能力決定的，而不是由反射光強度的絕對值決

2、定的；物體的色彩不受光照非均性的影響，具有一致性，即Retinex理論是以色感一致性（顏色恒常性）為基礎的。根據Edwin Land提出的理論，一幅給定的圖像S(x,y)分解成兩幅不同的圖像：反射物體圖像R(x,y)和入射光圖像L(x,y)，其原理示意圖如圖8.3-1所示。圖-1 Retinex理論示意圖對于觀察圖像S中的每個點(x,y)，用公式可以表示為： S(x,y)=R(x,y)×L(x,y) （1.3.1）實際上，Retinex理論就是通過圖像S來得到物體的反射性質R，也就是去除了入射光L的性質從而得到物體原本該有的樣子。2 基于Retinex理論的圖像增強的基本步驟步驟一:

3、 利用取對數的方法將照射光分量和反射光分量分離，即：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

4、)；步驟五：對R(x,y)做對比度增強，得到最終的結果圖像。3 多尺度Retinex算法 D Jobson等人提出了多尺度Retinex算法，多尺度算法的基本公式是： 其中，是Retinex的輸出，表示3個顏色譜帶，是高斯濾波函數，表示尺度的權重因子，表示使用尺度的個數，=3，表示彩色圖像，。=1，表示灰度圖像。從公式中可以看出：MSR算法的特點是能產生包含色調再現和動態范圍壓縮這兩個特性的輸出圖像。 在MSR算法的增強過程中，圖像可能會因為增加了噪聲而造成對圖像中的局部區域色彩失真，使得物體的真正顏色效果不能很好的顯現出來，從而影響了整體視覺效果。為了彌補這個缺點，一般情況下會應用帶色彩恢復

5、因子C的多尺度算法（MSRCR）來解決。帶色彩恢復因子C的多尺度算法(MSRCR)是在多個固定尺度的基礎上考慮色彩不失真恢復的結果，在多尺度Retinex算法過程中，我們通過引入一個色彩因子C來彌補由于圖像局部區域對比度增強而導致圖像顏色失真的缺陷，通常情況下所引入的色彩恢復因子C的表達式為其中,表示第個通道的色彩恢復系數，它的作用是用來調節3個通道顏色的比例，表示的是顏色空間的映射函數。帶色彩恢復的多尺度Retinex算法（MSRCR）通過色彩恢復因子C這個系數來調整原始圖像中三個顏色通道之間的比例關系，從而通過把相對有點暗的區域的信息凸顯出來，以達到消除圖像色彩失真的缺陷。處理后的圖像局域

6、對比度提高，而且它的亮度與真實的場景很相似，圖像在人們視覺感知下顯得極其逼真。因此，MSR算法具有較好的顏色再現性、亮度恒常性以及動態范圍壓縮等特性。4 例程精講例程1是基于Retinex理論進行霧靄天氣增強的MATLAB程序，讀者可結合程序及注釋對基于Retinex理論進行霧靄天氣增強的基本原理進行進一步分析，該程序的運行結果如圖-2所示。例程1：clear;close all;% 讀入圖像I=imread('wu.png');% 取輸入圖像的R分量R=I(:,:,1);N1,M1=size(R);% 對R分量進行數據轉換,并對其取對數R0=double(R);Rlog=lo

7、g(R0+1);% 對R分量進行二維傅里葉變換Rfft2=fft2(R0);% 形成高斯濾波函數sigma=250;F = zeros(N1,M1);for i=1:N1 for j=1:M1 F(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=do

8、uble(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(

9、G0);% 形成高斯濾波函數sigma=250;for i=1:N1 for j=1:M1 F(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;% 取反對數

10、，得到增強后的圖像分量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;for i=1:N1 for j=1:M1 F(i,j

11、)=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

12、= 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所示。例

13、程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);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:N1 for j=1:M1 F(i,j)=ex

14、p(-(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);for i=1:N1 f

15、or j=1:M1 F(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

16、(N1,M1);for i=1:N1 for j=1:M1 F(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;% 對上述三次增強得到的圖像取均值作為最終

17、增強的圖像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(EX

18、PRr); N1,M1=size(G);% 對G分量進行處理，步驟與對R分量處理的步驟相同，請讀者仿照R分量處理的步驟進行理解。G0=double(G);Glog=log(G0+1); Gfft2=fft2(G0);sigma=128;F = zeros(N1,M1);for i=1:N1 for j=1:M1 F(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(D

19、G);DGlog=log(DGdouble+1);Gg0=Glog-DGlog; sigma=256; F = zeros(N1,M1);for i=1:N1 for j=1:M1 F(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,M

20、1);for i=1:N1 for j=1:M1 F(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=imdivid

21、e(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);for i=1:N1 fo

22、r j=1:M1 F(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);for i=1:N1 for j=1:M1 F(i,j)=exp(-(i-N1/2)2+(j-M1/2)2)/(2*sigma*sigma); endendF = F./(s