1.3 霧靄圖像清晰化及其實(shí)現(xiàn)

隨著計(jì)算機(jī)技術(shù)、傳感器技術(shù)的飛速發(fā)展，數(shù)字圖像處理在智能交通、偵察勘探、導(dǎo)航定位等領(lǐng)域發(fā)揮的作用日益突出。

由于霧靄天氣的影響，會(huì)使能見(jiàn)度變低，圖像采集設(shè)備采集的圖像一般都會(huì)出現(xiàn)比較嚴(yán)重的退化與失真，會(huì)對(duì)智能交通、航拍測(cè)繪、視頻導(dǎo)航造成嚴(yán)重的影響。因此，通過(guò)現(xiàn)代圖像處理的方法提高霧靄天氣下拍攝圖像的能見(jiàn)度具有重要的作用。

1.3.1 Retinex理論

Retinex（視網(wǎng)膜Retina和大腦皮層Cortex的縮寫(xiě)）理論是一種建立在科學(xué)實(shí)驗(yàn)和科學(xué)分析基礎(chǔ)上的基于人類視覺(jué)系統(tǒng)（Human Visual System）的圖像增強(qiáng)理論。該算法的基本原理模型最早是由Edwin Land（埃德溫·蘭德）于1971年提出的理論，并在顏色恒常性的基礎(chǔ)上提出的一種圖像增強(qiáng)方法。Retinex理論的基本內(nèi)容是：物體的顏色是由物體對(duì)長(zhǎng)波（紅）、中波（綠）和短波（藍(lán)）光線的反射能力決定的，而不是由反射光強(qiáng)度的絕對(duì)值決定的；物體的色彩不受光照非均性的影響，具有一致性，即Retinex理論是以色感一致性（顏色恒常性）為基礎(chǔ)的。

根據(jù)Edwin Land提出的理論，一幅給定的圖像S（x,y）分解成兩幅不同的圖像：反射物體圖像R（x,y）和入射光圖像L（x,y），其原理示意圖如圖1.3.1所示。

對(duì)于給定圖像S中的每個(gè)點(diǎn)（x,y），用公式可以表示為：

實(shí)際上，Retinex理論就是通過(guò)圖像S來(lái)得到物體的反射性質(zhì)R，也就是去除了入射光L的性質(zhì)從而得到物體原本該有的樣子。

1.3.2 基于Retinex理論的圖像增強(qiáng)的基本步驟

步驟一：利用取對(duì)數(shù)的方法將照射光分量和反射光分量分離，即：

S′（x,y）＝r（x,y）+l（x,y）＝log（R（x,y））+log（L（x,y））

步驟二：用高斯模板對(duì)原圖像作卷積，即相當(dāng)于對(duì)原圖像作低通濾波，得到低通濾波后的圖像D（x,y）,F（x,y）表示高斯濾波函數(shù)：

D（x, y）＝S（x, y）*F（x, y）

步驟三：在對(duì)數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增強(qiáng)的圖像G（x,y）：

G（x,y）＝S′（x,y）－log（D（x,y））

步驟四：對(duì)G（x,y）取反對(duì)數(shù)，得到增強(qiáng)后的圖像R（x,y）：

R（x,y）＝exp（G（x,y））

步驟五：對(duì)R（x,y）作對(duì)比度增強(qiáng)，得到最終的結(jié)果圖像。

1.3.3 多尺度Retinex算法

Jobson D．等人提出了多尺度Retinex算法，多尺度算法（MSR）的基本公式是：

其中，Ri（x,y）是Retinex的輸出，i∈R,G,B，表示3個(gè)顏色譜帶，F（x,y）是高斯濾波函數(shù)，Wn表示尺度的權(quán)重因子，N表示使用尺度的個(gè)數(shù)。N＝3，表示彩色圖像，i∈R,G,B;N＝1，表示灰度圖像。從公式中可以看出：MSR算法的特點(diǎn)是能產(chǎn)生包含色調(diào)再現(xiàn)和動(dòng)態(tài)范圍壓縮這兩個(gè)特性的輸出圖像。

在MSR算法的增強(qiáng)過(guò)程中，圖像可能會(huì)因?yàn)樵黾恿嗽肼暥斐蓪?duì)圖像中的局部區(qū)域色彩失真，使得物體的真正顏色效果不能很好地顯現(xiàn)出來(lái)，從而影響了整體視覺(jué)效果。為了彌補(bǔ)這個(gè)缺點(diǎn)，一般情況下會(huì)采用帶色彩恢復(fù)因子C的多尺度算法（MSRCR）來(lái)解決。帶色彩恢復(fù)因子C的多尺度算法是在多個(gè)固定尺度的基礎(chǔ)上考慮色彩不失真恢復(fù)的結(jié)果，在多尺度Retinex算法過(guò)程中，通過(guò)引入一個(gè)色彩因子C來(lái)彌補(bǔ)由于圖像局部區(qū)域?qū)Ρ榷仍鰪?qiáng)而導(dǎo)致的圖像顏色失真的缺陷，通常情況下所引入的色彩恢復(fù)因子C的表達(dá)式為：

其中，Ci表示第i個(gè)通道的色彩恢復(fù)系數(shù)，它的作用是調(diào)節(jié)3個(gè)通道顏色的比例，f（）表示顏色空間的映射函數(shù)。帶色彩恢復(fù)的多尺度Retinex算法通過(guò)色彩恢復(fù)因子C這個(gè)系數(shù)來(lái)調(diào)整原始圖像中3個(gè)顏色通道之間的比例關(guān)系，從而把相對(duì)有點(diǎn)暗的區(qū)域的信息凸顯出來(lái)，以達(dá)到消除圖像色彩失真缺陷的目的。處理后的圖像局域?qū)Ρ榷鹊靡蕴岣?，而且其亮度與真實(shí)的場(chǎng)景很相似，圖像在人們的視覺(jué)感知下顯得極其逼真。

1.3.4 實(shí)例精講

【例1.3.1】 基于Retinex理論進(jìn)行霧靄天氣增強(qiáng)的MATLAB程序，讀者可結(jié)合程序及注釋對(duì)基于Retinex理論進(jìn)行霧靄圖像清晰化的基本原理進(jìn)行深入分析。

    *********************************************************************************
    clear;
    close all;
    % 讀入圖像
    I＝imread(＇wu.png＇);
    % 取輸入圖像的R分量
    R＝I(:, :,1);
    [N1, M1]＝size(R);
    % 對(duì)R分量進(jìn)行數(shù)據(jù)轉(zhuǎn)換，并對(duì)其取對(duì)數(shù)
    R0＝double(R);
    Rlog＝log(R0+1);
    % 對(duì)R分量進(jìn)行二維傅里葉變換
    Rfft2＝fft2(R0);
    % 形成高斯濾波函數(shù)
    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));
            end
    end
    F ＝ F./(sum(F(:)));
    % 對(duì)高斯濾波函數(shù)進(jìn)行二維傅里葉變換
    Ffft＝fft2(double(F));
    % 對(duì)R分量與高斯濾波函數(shù)進(jìn)行卷積運(yùn)算
    DR0＝Rfft2.*Ffft;
    DR＝ifft2(DR0);
    % 在對(duì)數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增強(qiáng)的圖像
    DRdouble＝double(DR);
    DRlog＝log(DRdouble+1);
    Rr＝Rlog－DRlog;
    % 取反對(duì)數(shù)，得到增強(qiáng)后的圖像分量
    EXPRr＝exp(Rr);
    % 對(duì)增強(qiáng)后的圖像進(jìn)行對(duì)比度拉伸增強(qiáng)
    MIN ＝ min(min(EXPRr));
    MAX ＝ max(max(EXPRr));
    EXPRr ＝(EXPRr－MIN)/(MAX－MIN);
    EXPRr＝adapthisteq(EXPRr);
    % 取輸入圖像的G分量
    G＝I(:, :,2);
    [N1, M1]＝size(G);
    % 對(duì)G分量進(jìn)行數(shù)據(jù)轉(zhuǎn)換，并對(duì)其取對(duì)數(shù)
    G0＝double(G);
    Glog＝log(G0+1);
    % 對(duì)G分量進(jìn)行二維傅里葉變換
    Gfft2＝fft2(G0);
    % 形成高斯濾波函數(shù)
    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));
           end
    end
    F ＝ F./(sum(F(:)));
    % 對(duì)高斯濾波函數(shù)進(jìn)行二維傅里葉變換
    Ffft＝fft2(double(F));
    % 對(duì)G分量與高斯濾波函數(shù)進(jìn)行卷積運(yùn)算
    DG0＝Gfft2.*Ffft;
    DG＝ifft2(DG0);
    % 在對(duì)數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增強(qiáng)的圖像
    DGdouble＝double(DG);
    DGlog＝log(DGdouble+1);
    Gg＝Glog－DGlog;
    % 取反對(duì)數(shù)，得到增強(qiáng)后的圖像分量
    EXPGg＝exp(Gg);
    % 對(duì)增強(qiáng)后的圖像進(jìn)行對(duì)比度拉伸增強(qiáng)
    MIN ＝ min(min(EXPGg));
    MAX ＝ max(max(EXPGg));
    EXPGg ＝(EXPGg－MIN)/(MAX－MIN);
    EXPGg＝adapthisteq(EXPGg);
    % 取輸入圖像的B分量
    B＝I(:, :,3);
    [N1, M1]＝size(B);
    % 對(duì)B分量進(jìn)行數(shù)據(jù)轉(zhuǎn)換，并對(duì)其取對(duì)數(shù)
    B0＝double(B);
    Blog＝log(B0+1);
    % 對(duì)B分量進(jìn)行二維傅里葉變換
    Bfft2＝fft2(B0);
    % 形成高斯濾波函數(shù)
    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));
           end
    end
    F ＝ F./(sum(F(:)));
    % 對(duì)高斯濾波函數(shù)進(jìn)行二維傅里葉變換
    Ffft＝fft2(double(F));
    % 對(duì)B分量與高斯濾波函數(shù)進(jìn)行卷積運(yùn)算
    DB0＝Gfft2.*Ffft;
    DB＝ifft2(DB0);
    % 在對(duì)數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增強(qiáng)的圖像
    DBdouble＝double(DB);
    DBlog＝log(DBdouble+1);
    Bb＝Blog－DBlog;
    EXPBb＝exp(Bb);
    % 對(duì)增強(qiáng)后的圖像進(jìn)行對(duì)比度拉伸增強(qiáng)
    MIN ＝ min(min(EXPBb));
    MAX ＝ max(max(EXPBb));
    EXPBb ＝(EXPBb－MIN)/(MAX－MIN);
    EXPBb＝adapthisteq(EXPBb);
     % 對(duì)增強(qiáng)后的圖像R、G、B分量進(jìn)行融合
    I0(:, :,1)＝EXPRr;
    I0(:, :,2)＝EXPGg;
    I0(:, :,3)＝EXPBb;
    % 顯示運(yùn)行結(jié)果
    subplot(121), imshow(I);
    subplot(122), imshow(I0);
    *********************************************************************************

該程序的運(yùn)行結(jié)果如圖1.3.2所示。

【例1.3.2】 基于Retinex理論進(jìn)行霧靄圖像清晰化的MATLAB程序，讀者可結(jié)合程序及注釋對(duì)基于Retinex理論進(jìn)行霧靄圖像清晰化的基本原理進(jìn)行進(jìn)一步分析。

    *********************************************************************************
    clear;
    close all;
    I＝imread(＇wu.png＇);
    % 分別取輸入圖像的R、G、B三個(gè)分量，并將其轉(zhuǎn)換為雙精度型
    R＝I(:, :,1);
    G＝I(:, :,2);
    B＝I(:, :,3);
    R0＝double(R);
    G0＝double(G);
    B0＝double(B);
    [N1, M1]＝size(R);
    % 對(duì)R分量進(jìn)行對(duì)數(shù)變換
    Rlog＝log(R0+1);
    % 對(duì)R分量進(jìn)行二維傅里葉變換
    Rfft2＝fft2(R0);
    % 形成高斯濾波函數(shù)(sigma＝128)
    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));
            end
    end
    F ＝ F./(sum(F(:)));
    % 對(duì)高斯濾波函數(shù)進(jìn)行二維傅里葉變換
    Ffft＝fft2(double(F));
    % 對(duì)R分量與高斯濾波函數(shù)進(jìn)行卷積運(yùn)算
    DR0＝Rfft2.*Ffft;
    DR＝ifft2(DR0);
    % 在對(duì)數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增強(qiáng)的圖像
    DRdouble＝double(DR);
    DRlog＝log(DRdouble+1);
    Rr0＝Rlog－DRlog;
    % 形成高斯濾波函數(shù)(sigma＝256)
    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));
            end
    end
    F ＝ F./(sum(F(:)));
    % 對(duì)高斯濾波函數(shù)進(jìn)行二維傅里葉變換
    Ffft＝fft2(double(F));
    % 對(duì)R分量與高斯濾波函數(shù)進(jìn)行卷積運(yùn)算
    DR0＝Rfft2.*Ffft;
    DR＝ifft2(DR0);
    % 在對(duì)數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增強(qiáng)的圖像
    DRdouble＝double(DR);
    DRlog＝log(DRdouble+1);
    Rr1＝Rlog－DRlog;
    % 形成高斯濾波函數(shù)(sigma＝512)
    sigma＝512;
    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));
          end
    end
    F ＝ F./(sum(F(:)));
    % 對(duì)高斯濾波函數(shù)進(jìn)行二維傅里葉變換
    Ffft＝fft2(double(F));
    % 對(duì)R分量與高斯濾波函數(shù)進(jìn)行卷積運(yùn)算
    DR0＝Rfft2.*Ffft;
    DR＝ifft2(DR0);
    % 在對(duì)數(shù)域中，用原圖像減去低通濾波后的圖像，得到高頻增強(qiáng)的圖像
    DRdouble＝double(DR);
    DRlog＝log(DRdouble+1);
    Rr2＝Rlog－DRlog;
    % 對(duì)上述三次增強(qiáng)得到的圖像取均值，作為最終增強(qiáng)的圖像
    Rr＝(1/3)*(Rr0+Rr1+Rr2);
    % 定義色彩恢復(fù)因子C
    a＝125;
    II＝imadd(R0, G0);
    II＝imadd(II, B0);
    Ir＝immultiply(R0, a);
    C＝imdivide(Ir, II);
    C＝log(C+1);
    % 將增強(qiáng)后的R分量乘以色彩恢復(fù)因子，并對(duì)其進(jìn)行反對(duì)數(shù)變換
    Rr＝immultiply(C, Rr);
    EXPRr＝exp(Rr);
    % 對(duì)增強(qiáng)后的R分量進(jìn)行灰度拉伸
    MIN ＝ min(min(EXPRr));
    MAX ＝ max(max(EXPRr));
    EXPRr ＝(EXPRr－MIN)/(MAX－MIN);
    EXPRr＝adapthisteq(EXPRr)
    [N1, M1]＝size(G);
    % 對(duì)G分量進(jìn)行處理，步驟與對(duì)R分量處理的步驟相同，可仿照R分量處理的步驟進(jìn)行理解
    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));
            end
    end
    F ＝ 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:N1
            for j＝1:M1
            F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
            end
    end
    F ＝ 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:N1
          for j＝1:M1
            F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
          end
    end
    F ＝ 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);
    % 對(duì)B分量進(jìn)行處理，步驟與對(duì)R分量處理的步驟相同，可仿照R分量處理的步驟進(jìn)行理解
    [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:M1
            F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
          end
    end
    F ＝ 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));
          end
    end
    F ＝ 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:N1
          for j＝1:M1
            F(i, j)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
          end
    end
    F ＝ 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);
    % 對(duì)增強(qiáng)后的圖像R、G、B分量進(jìn)行融合
    I0(:, :,1)＝EXPRr;
    I0(:, :,2)＝EXPGg;
    I0(:, :,3)＝EXPBb;
    % 顯示運(yùn)行結(jié)果
    subplot(121), imshow(I);
    subplot(122), imshow(I0);
    *********************************************************************************

該程序的運(yùn)行結(jié)果如圖1.3.3所示。