1.3 霧靄圖像清晰化及其實現

隨著計算機技術、傳感器技術的飛速發展，數字圖像處理在智能交通、偵察勘探、導航定位等領域發揮的作用日益突出。

由于霧靄天氣的影響，會使能見度變低，圖像采集設備采集的圖像一般都會出現比較嚴重的退化與失真，會對智能交通、航拍測繪、視頻導航造成嚴重的影響。因此，通過現代圖像處理的方法提高霧靄天氣下拍攝圖像的能見度具有重要的作用。

1.3.1 Retinex理論

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

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

對于給定圖像S中的每個點（x,y），用公式可以表示為：

實際上，Retinex理論就是通過圖像S來得到物體的反射性質R，也就是去除了入射光L的性質從而得到物體原本該有的樣子。

1.3.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）作對比度增強，得到最終的結果圖像。

1.3.3 多尺度Retinex算法

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

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

在MSR算法的增強過程中，圖像可能會因為增加了噪聲而造成對圖像中的局部區域色彩失真，使得物體的真正顏色效果不能很好地顯現出來，從而影響了整體視覺效果。為了彌補這個缺點，一般情況下會采用帶色彩恢復因子C的多尺度算法（MSRCR）來解決。帶色彩恢復因子C的多尺度算法是在多個固定尺度的基礎上考慮色彩不失真恢復的結果，在多尺度Retinex算法過程中，通過引入一個色彩因子C來彌補由于圖像局部區域對比度增強而導致的圖像顏色失真的缺陷，通常情況下所引入的色彩恢復因子C的表達式為：

其中，Ci表示第i個通道的色彩恢復系數，它的作用是調節3個通道顏色的比例，f（）表示顏色空間的映射函數。帶色彩恢復的多尺度Retinex算法通過色彩恢復因子C這個系數來調整原始圖像中3個顏色通道之間的比例關系，從而把相對有點暗的區域的信息凸顯出來，以達到消除圖像色彩失真缺陷的目的。處理后的圖像局域對比度得以提高，而且其亮度與真實的場景很相似，圖像在人們的視覺感知下顯得極其逼真。

1.3.4 實例精講

【例1.3.1】 基于Retinex理論進行霧靄天氣增強的MATLAB程序，讀者可結合程序及注釋對基于Retinex理論進行霧靄圖像清晰化的基本原理進行深入分析。

    *********************************************************************************
    clear;
    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: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));
    % 對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;
    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));
    % 對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;
    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));
    % 對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);
    *********************************************************************************

該程序的運行結果如圖1.3.2所示。

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

    *********************************************************************************
    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)＝exp(－((i－N1/2)＾2+(j－M1/2)＾2)/(2*sigma*sigma));
            end
    end
    F ＝ 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
            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));
    % 對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);
    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));
    % 對R分量與高斯濾波函數進行卷積運算
    DR0＝Rfft2.*Ffft;
    DR＝ifft2(DR0);
    % 在對數域中，用原圖像減去低通濾波后的圖像，得到高頻增強的圖像
    DRdouble＝double(DR);
    DRlog＝log(DRdouble+1);
    Rr2＝Rlog－DRlog;
    % 對上述三次增強得到的圖像取均值，作為最終增強的圖像
    Rr＝(1/3)*(Rr0+Rr1+Rr2);
    % 定義色彩恢復因子C
    a＝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);
    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);
    % 對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
          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);
    % 對增強后的圖像R、G、B分量進行融合
    I0(:, :,1)＝EXPRr;
    I0(:, :,2)＝EXPGg;
    I0(:, :,3)＝EXPBb;
    % 顯示運行結果
    subplot(121), imshow(I);
    subplot(122), imshow(I0);
    *********************************************************************************

該程序的運行結果如圖1.3.3所示。