視覺式交通監控之能見度改善

全文

(1)第一章簡介 1.1 研究動機「攝影」一直是紀錄事物的一種方便且又有效率好方法。尤其是近年來數位相機、數位攝影機的發展，重量是越做越輕，尺寸也是越來越輕薄短小，也不需要拿到照相館去沖洗相片，直接傳輸到電腦甚至在機器上就可以直接觀賞，更是大大的增加了「攝影」的便利性。也就是因為數位相機、數位攝影機的普及，造成了許許多多相關的研究如雨後春筍般蓬勃的展開。其中，設置在各種交通要道的視覺式交通監控系統，就是一種重要的應用。利用各種設置在路口的交通監視器，來達到車輛的計數、流量的監控、甚至交通號誌的調整等等。「攝影」的環境大致上來說可以分成兩種情況。第一種環境是指室內拍攝。這種拍攝環境通常比較單純，也較少受到一些雜訊的干擾。在這種環境之下的照片品質通常也會比較高，只是所能應用的種類就比較少了。第二種環境是指在室外拍攝。可想而知，這種拍攝環境就比室內拍攝複雜的多了，常常受到外界的影響，而拍出來的品質也會比較差，但是所能應用的種類就相對的多的多了。上述所提到的視覺式交通監控系統，在大多數的情況下皆屬於室外的應用。只要是室外的應用，或多或少都會受到天氣的影響。尤其是起霧、下雨或者是下雪等惡劣的天氣，更會影響照片的能見度，進而對我們的視覺式交通監視系統造成辨識上的困難。許許多多如車輛計數的誤差，交通流量的誤判等，皆因壞天氣的影響而更加嚴重。但實際上，在許多戶外的應用與實驗當中，常常假設相 1.

(2) 機或人眼所接收景物的光不會因為空氣的傳遞而衰退或改變。換句話說，就是假設實驗的環境天氣晴朗且能見度良好，所要觀察的目標物可以很清楚的被看到。在此假設下，相機或人眼所接受的光，完完全全是景物所散發出的光。但此假設，在現實環境中不可能永遠存在。所以，一個可靠的系統，必須要能容忍在不同的環境中(如：起霧、下雨、下雪等)也能確實的運作。就是此動機讓我們想研究如何使戶外的視覺式交通監控系統能容忍各種的天氣環境，並且更加的穩定。. 1.2 研究目的當景物表面發出的光經空中傳遞到達相機時，因為與空氣交互作用，造成亮度與色彩的改變。我們把這些交互作用分成主要三類，分別是：散射(scattering)、吸收(absorption)與放射(emission)。其中，又以散射對光的影響最為嚴重。散射的形成是由於空氣中飄浮著微小的粒子，當光行經空中時與粒子產生作用而形成散射。所以，空氣中不同程度的粒子分布，就會影響最終所看到的圖片影像。圖 1.1 為起霧時從山上鳥瞰高速公路的情景，很明顯的整張圖看起來霧茫茫的，能見度並不佳。圖 1.2 為下雨時在從路旁拍攝高速公路的情景。圖 1.3 為下雨時從天橋拍攝一般道路的情景。. 2.

(3) 圖 1.1 高速公路起霧時. 圖 1.2 高速公路下雨時. 圖 1.3 一般道路下雨時. 天氣之所以有不同的狀態(如：起霧、下雨、下雪等)最主要是由於不同狀態且不同大小的懸浮微粒散布在大氣中所造成的。以下介紹幾個常見的壞天氣主要的成因：見表 1.1 所示，整理出不同狀態與不同大小的懸浮微粒與不同天氣的關係。. 3.

(4) 天氣種類微粒狀態. 半徑 ( µ m). 濃度 (cm3 ). 空氣. 微小顆粒. 10−4. 1019. 霾. 漂浮物質. 10−2 - 1. 103 - 10. 霧. 小水滴. 1 - 10. 100 - 10. 雲. 小水滴. 1 - 10. 300 - 10. 雨. 水滴. 102 - 104 10-2 - 10-5. 表 1.1 各種天氣狀態與粒子的關西. 下列詳細定義不同的天氣狀態：霾：是由散布在大氣之中的漂浮物質所構成的。這些漂浮物質的成分包括火山灰、植物的分泌物、燃燒後的餘燼、灰塵、海鹽等。最初，漂浮在空氣中的只有一些微小的水分子，水分子漸漸的吸附了空氣周遭的灰塵等汙染物而越來越大，最終形成了霾。其質點極其細微，肉眼無法辨識。霾在大氣中多呈乳白色，惟如對遠地明亮之背景，則呈黃色或橘紅色；如對較陰暗之背景，則顯示淡藍色。霾的分子大小比空氣還大，但是比霧小。霧：當大氣中的相對濕度接近於飽和的時候，上述構成霾的漂浮物質會慢慢的凝結成小水滴，最終形成了霧。而我們常常聽到的靄則是介於霾與霧之間的狀態。霧有數種不同類型，主要根據凝固的成因不同分辨。雲：組成的成分大致與霧相同。霧與雲的主要分別是霧接觸地面，而雲則不會。 4.

(5) 因為雲反射和散射所有波段的電磁波，所以雲的顏色成灰度色，雲層比較薄時成白色，但是當它們變得太厚或濃密而使得陽光不能通過的話，它們可以看起來是灰色或黑色的。雨：當雲為水成雲或冰成雲時，雲與雲之間互相碰撞結合使得雲中水滴變大，當水的重力大到上升氣流無法將其「托住」時，水滴下降變形成了雨，實際上，雨在下落時曾經做過數次垂直運動，這是由上升氣流的強與弱有關的。我們看到的影像，就是因為受到大氣中懸浮微粒的影響，而造成亮度或色度的改變。我們將在下一章節詳細的討論當景物表面發出的光經空中傳遞到達相機時，受到懸浮微粒的散射所影響的影像公式。我們最主要的目的，不僅僅是想了解光與空氣是如何的作用。更進一步的，希望藉由公式的推導，回復光在受未干擾時的狀態，也就是景物表面最初所散發的光。. 1.3 文獻探討在影像處理的領域裡，越來越多的研究致力於影像的恢復以及影像的分析。其中，影像的分析是為了獲得影像內容的資訊來幫助我們做各種的應用。例如在文獻[Lev 07]中，利用改變相機鏡頭的孔徑來估計景物的距離，建立景深圖(depth map)。在[Yua 07]中利用兩張圖，一張雜訊圖而另一張模糊的圖來達到去模糊的目的。判斷邊緣的銳利程度也在[Fat 07]中僅需要一張輸入影像就能做到。而最常用到的雜訊的估計在[Liu 06]也只利用一張輸入的影像就可做到。 5.

(6) 改善壞天氣對視覺的影響一直是一個非常具有挑戰的問題。大多數的文獻在解決此問題時，常常會利用各種附加的資訊或者人工的調整來幫助。在[Tan 00] 中假設景物的距離都已經知道，來改善從空中俯瞰地面的影像中，所受到各種天氣造成的影響。在[Sch 01]中利用兩張圖片來去除被偏極化的霧。此兩張圖片必須是同一個場景，並且對應於偏光鏡中不同的角度，再設法預估偏極化的程度 (degree of polarization)來達到除霧的目的。進一步的在[Shw 06]中研究在高頻的空間頻道(spatial-bands)中自動估計利用偏極化的方法除霧所需要的參數。這些方法都成功改善壞天氣對於視覺所造成的影響。在[Sch 07]中更進一步的解決了上述一系列的研究中，所產生共同的缺點─noise在運算的過程中被放大的現象。而在[Nar 03]中也利用復原的影像來估計景物的距離，並且重建了真實景物的結構圖。在[Oak 07]中假設環境光(airlight)為整張圖中一常數，並假設此一常數已經得到，觀察自然的影像發現，像素點亮度的平均值與標準差成比例。在[Tan 08] 中假設環境光在整張影像中為平滑且均勻的前提下，加強整張影像的對比程度來改善壞天氣對視覺所造成的影響。此一方法也強調並不是還原景物原本的樣貌，而是使對比程度增強使景物看起來更加清晰。大氣中雲霧的存在，也影響了衛星照片、高空影像等地面的監控。多光譜影像的拍攝，是由衛星從太空中拍攝地球所反射的六種頻道的光。這種影像很容易因為雲的存在、空氣的汙染甚至於沙塵暴等影響其品質，所以很多的研究致力於 6.

(7) 解決這些影響。Dark-object subtraction[Cha 88]，此方法假設雲霧的分布都是均勻的狀態下，對於每個頻道最黑的物體都減去一個常數，去除雲霧所造成的影響。而此常數利用亮度的直方圖並且以人為的判斷來決定。[Zha 02]中，對於兩個受雲霧影響比較大的頻道中計算haze-optimized transform使其影響減至最低。對比程度的加強也是許多研究的重點。例如線性映射、直方圖平衡和伽瑪更正都常出現在商用的影像處理軟體上。另外還有許多的研究如：Laplacian pyramid [Rah 96]、wavelet decomposition [Lu 94]、 single scale unsharp-mask filter [Wik 07] 和 the multi-scale bilateral filter [Fat 07]等都是利用放大亮度的變化量使對比提高。承如之前所提到的，雲霧對於影像的影響大致可以視為在影像的亮度中加入一環境光常數。所以，如何預估環境光常數就顯得非常重要了。. 7.

(8) 第二章利用偏光鏡回復影像 2.1 影像退化(Degradation)模式模式影像退化. 圖 2.1 成像過程. 見圖 2.1 所示，當景物表面發出的光經空中傳遞到達相機時，我們所看到的影像是由兩種光源所構成的。第一種光來自於景物表面，經空中傳遞衰退後到達相機，這種光我們稱它為 direct light。第二種光是來自於周遭環境(ambient)的光，例如陽光(sunlight)、天光(skylight)、周遭物體反射的光等，在受到大氣中(如空氣、水氣、懸浮)粒子的散射，其中朝向像機所在的方向傳遞的光，我們稱它為 airlight。令 D 表 direct light，定義如下： D = J ( x )t ( x ) ,. (2.1). 其中 x 表任一像點的位置， J ( x ) 為景物表面所散發出的光， z. t = exp[− ∫ β( z′, λ )dz ′] 表大氣的透射率(transmittance of incoherent light)，其與景物 0. 和相機之間的距離 z 以及大氣的散射係數 β 有關，當衰退係數與位置無關時， 8.

(9) β( z ′, λ ) 可令為一常數，亦即 β( z ′, λ ) = β ；t 於是可寫成 z. z. 0. 0. t = exp[− ∫ βdz′] = exp[−β ∫ dz ′] = exp[−βz ] 。. 圖 2.2 Airlight 的計算. 有關 airlight 的估計見圖 2.2 所示，其乃累積由觀察者至景物一表面積所成的角錐狀體積內之環境光，代表了相機與景物之間的 airlight。令 z 為觀察者與景物之間的距離， d ω 為角錐體的 solid angle， dV 則表示在距觀察者 z′ 處的一極小的體積，其計算式為 dV = d ω z′2 dz ′ 。在 dV 內的環境光因散射而被我們觀察到的. intensity 可以用下列式子表示： dI ( z ′, λ ) = dVkβ(λ ) = dω z ′2 dz ′kβ(λ ). (2.2). 其中 β(λ ) 表對波長為 λ 的光之衰退係數，令其與波長無關亦即 β(λ ) = β ， k 表單位體積的環境光。因此(2.2)式可以表示成： dI ( z ′) = dVkβ = d ω z ′2 dz′kβ dI ( z′) 在傳到觀察者時，會衰退為 dE ( z ′) = 9. dI ( z′)e-βz′ ， z ′2. (2.3). (2.4).

(10) 此為觀察者實際所接受到的 irradiance。而由 dV 來至觀察者的 radiance 即為 dL ( z′) =. dE ( z ′) dI ( z ′)e-βz ′ = dω d ω z ′2. 把(2.3)式代入(2.5)式 dL ( z′) =. (2.5) dI ( z ′)e-βz ′ d ω z ′2 dz′kβe-βz′ = = kβe-βz′ dz′，將其對 z′ 從 0 2 2 dω z′ dω z′. 積分到 z 得到 L( z ) = k (1 − t ) 其中 t = e. -βz. (2.6). 。如果景物在地平線上無窮遠處時 z = ∞ ， t = 0 ，我們得到. Lh (∞) = k 。(2.6)式變成 L( z ) = Lh (∞)(1 − t ) 。令 Alight 表 Airlight 即依此定義如下 Alight = A[1 − t ( x )] ,. (2.7). 其中 A 為一常數向量，代表在無窮遠處尚未因距離衰退之單位體積環境光之彩色向量。前述影像乃是由 direct light ((2.1)式) 和 airlight ((2.7)式) 這兩種光源相加而成 I ( x ) = t ( x)J ( x ) + (1 − t ( x )) A .. (2.8). 上式顯示 direct light 和 airlight 是以線性相加，其中 t ( x)J ( x ) 說明景物表面散發出的光 J ( x ) 隨距離衰減，而 (1 − t ( x )) A 說明只有 A 不被衰退的部份才能到達觀察者。本研究的目的是要由已知的影像， I ( x ) 求得在沒有大氣衰退和 airlight 影響下的景物表面散發出的光 J ( x ) 。. 2.2*光的偏振性就先前對光的認識而言，光是具有二相性的，從整體上來看，波的性質較為顯著，但從細微的角度來看，粒子性則顯得較為重要。 10.

(11) 一般光源含有大量的發光原子或分子，它們所發出的光的振動方向彼此毫不相關。但由於發光原子和分子的數量相當龐大，因而在與光波傳播方向垂直的平面上，與任何一個振動方向相聯繫的原子或分子數目基本上是相等的。換言之，這樣的光源在一切可振動方向上的光振動機率是相等的。從能量的角度上來看，光振動平面內任何方向上都具有相同的能量，或者說光振動的功率密度都相同，這種振動方向漫無規律的光波就是自然光。如下圖 2.3 所示，圖中光進行的方向與紙面垂直，我們利用箭號表示振動的方向；若振動的方向沒有偏好(亦即所有方向都會均勻地出現)，就說這道光是未偏極化未偏極化的。未偏極化. 圖 2.3 自然光示意圖. 由於光是電磁波的一種，到目前為止我們僅考慮在固定平面上的電場向量；而此電場向量所在之平面也稱為電場振動平面；而此時在固定平面振動傳播之光波稱為偏振光。若光波的波向量總是沿著一條直線作反覆振動，則稱為線偏振光 11.

(12) (linearly polarized light)。線偏振光的振動面是固定的平面，故也稱為平面偏振光。如下圖 2.4 為一極端的例子，前者只有水平的振動方向、而後者只有垂直的振動方向。. 圖 2.4 線偏振光示意圖. 若線偏振光與自然光相摻雜，則為部分偏振光(partially polarized light)，這時發出沿線偏振方向光振動的原子或分子比發出沿其他方向光振動的原子或分子要多，因而這個方向光振動的功率密度比其他方向的光振動的功率密度要大。通常用偏振度(degree of polarization)p 來量度線偏振的程度與一切振動相同，光振動也可分解成兩個方向相垂直的振動。自然光中的每一個光振動都可分解成這兩個互相垂直的光振動。由於自然光中各方向上的光振動都相等，因而它們在這二方向上的分振動之和也相等，即自然光可用功率密度相等、振動方向互相垂直的兩個線偏振光來表示。且每個線偏振光的功率密度都 12.

(13) 是自然光功率密度的一半大自然中不少光線都是偏振光，例如玻璃表面的反映、花草的反光、水氣的反射等等都是例子之一，在攝影角度而言，這些光線都令照片色彩減弱，可免則免。偏光鏡以一種平衡排列的結晶體構成，這種結構只讓與晶體平衡的光波通過，向其他角度震動的光波就一律被阻擋，如下面示意圖：. 圖 2.5 偏光鏡示意圖. 本研究就是利用光的偏振性與偏光鏡的使用來求得在沒有大氣衰退和. airlight 影響下的景物表面散發出的光 J ( x ) 。 *註：此章節資料來自於. http://ykuo.ncue.edu.tw/report/045-Polarization%20of%20Light.doc 2.3 影像公式在 2.1 節裡有提到，airlight 是因為受到大氣中懸浮粒子的散射，而朝向相機 13.

(14) 的方向傳遞的光。在這散射折射的過程中光常常會被偏極化。因此我們假設. airlight 為部分偏振光，並且利用偏光鏡的調整來估計 airlight 的量。當轉動偏光鏡時，會因為偏光鏡濾掉其他方向的偏振光而使亮度有所改變。當偏光鏡轉到某一角度時，airlight 對亮度的貢獻最少，我們稱此狀態為最佳狀態(best state)。此時偏光鏡濾掉了大多數被偏極化的 airlight 表 A(x ) min 。而當偏光鏡轉到某另一角度時，使 airlight 的亮度最強表 A(x ) max 。通常 A(x )min 與 A(x ) max 其偏光鏡所相差的角度為 90 度。所以((2.7)式)中的 airlight 可以表示為：. Alight = A(x )min + A(x )max. (2.9). 假設 direct light ((2.1)式)未被偏極化，並且能量平均分散於上述兩種狀態上。因此透過偏光鏡被觀察到的 intensity 可以用下列式子表示：. I (x )min = D(x ) / 2 + A(x )min , I (x )max = D (x ) / 2 + A(x )max. (2.10). 通常仿照可見度的定義用偏振度(degree of polarization)p 來量度偏振的程度。因此 airlight 的偏振度定義如下：. p (x ) ≡ ( A(x ) max − A(x )min ) / Alight. (2.11). 在狹窄的視野裡 p 的變化並不大，所以假設 p 為一常數，且 0 ≤ p ≤ 1 。所以 A(x )min = Alight (1- p )/2， A(x ) max = Alight (1+ p )/2。因此((2.9)式)又可以表示成：. I (x )min = D(x ) / 2 + Alight (1- p )/2，I (x )max = D (x ) / 2 + Alight (1+ p )/2 並且，可以很簡單的從((2.8)式)和((2.12)式)得到： Iˆ (x ) total = I (x ) max + I (x ) min 14. (2.12).

(15) 2.4 恢復影像根據上面的影像公式，如果要恢復原本的影像 J ( x ) 。首先必須先估計. airlight。方法如下： Aˆ light = ( I (x )max − I (x )min ) / p. 再來，必須估計大氣的透射率： tˆ( x ) = 1 − Aˆ light /A. (2.13) (2.14). 最後，再利用上面求得的 airlight 和大氣透射率恢復原本的影像： Jˆ (x )=(I (x ) total -Aˆ light )/tˆ( x ). (2.15). 上述的方法有兩個問題存在：第一，在估計 airlight 時必須要先知道常數 p 的值。第二，在估計 t 時必須要先知道向量 A 的值。下面的章節就在介紹我們如何得到這兩個參數。. 2.5 估計 p 本章節利用了 ICA 的概念，使 airlight 與 direct light 分開，進而從中找出我們需要的參數 p。 2.5.1 Facilitating 線性的 ICA 根據((2.12)式)，把得到的兩張影像寫成下列的方程式：.  I max   Alight   min  = M    D  I   (1 + p ) / 2 1/ 2  其中， M=    (1 − p ) / 2 1/ 2 . (2.16). (2.17). 15.

(16)  Aˆ   I max  因此，要估計的部分就可以寫成  light  = W  min   Dˆ  I . (2.18). −1/ p   1/ p 其中， W=   ( p − 1) / p ( p + 1) / p . (2.19). ICA 最主要的觀念為：找出一個矩陣 W 產生出「好的」 Aˆ light 與 Dˆ 。為了使用 ICA，必須假設 airlight 與 direct light 為互相統計獨立(statistical independent)的，但 airlight 與 direct light 皆與距離 z 有關。Airlight 會隨著距離增加而變大而 direct. light 則會隨著距離增加而遞減。所以上述假設在實際的情況下並不成立。非常幸運的，統計相依(statistical dependence)並不會發生在所有的影像區域。事實上，上述相依的情形只會發生在 very low 的空間頻率上。因此，必須把原本的影像轉換到更適合的表示方法，並且排除低頻的部分。在這裡，我們使用小波轉換(wavelet transformation)，定義如下： Dc ( x ) = w{D( x )}. (2.20). 其中 c 代表小波的頻道， w 代表小波線性轉換。同樣的把 Alight ， Aˆ light ， Dˆ ， I max 與 I min 都分別執行小波轉換並且得到：Alightc，Aˆ lightc，Dˆ c，I max c 與 I min c。因此((2.18) 式)經過轉換後得到：  Aˆ lightc   I c max    = W  min  ˆ  Ic   Dc . (2.21). 2.5.2 Haar 離散小波轉換小波轉換為一種將資料由空間域轉變成頻率域方法。其中 Haar 小波轉換為眾多小波轉換中最簡單的一種。Haar 函數所使用之離散小波轉換使用的方法是將影像的所有像素視為分別是各自獨立的數值，並對這些數值做相加、相減的運 16.

(17) 算，以求得這張影像的對於頻率上的分頻。其中相加後的值會越來越大，值越大也就代表該值越重要，此部份就是低頻的部份。相反地，相減後的值代表著像素與像素之間的差距，如同我們做邊緣偵測時的方法一樣，當遇到影像中的物體邊緣時，像素間的差值就會很大；所以在 Haar 函數離散小波轉換中相減的部份，會明顯的強調出物體的邊緣，此部份是高頻的部份。Haar 離散小波轉換大致可以分為兩步驟，接下來我們將以 4 × 4 的影像為例詳細的說明 Haar 離散小波轉換的過程：若給予一大小為 4 × 4 的影像 B，且 B 為. a. b. c. d. e. f. g. h. i. j. k. l. m. n. o. p. 圖 2.6 4 × 4 的影像. 則第一步驟水平的小波轉換的結果為：. 17.

(18) a′. b′. c′. d′. e′. f′. g′. h′. i′. j′. k′. l′. m′. n′. o′. p′. 圖 2.7 水平小波轉換後的影像. 其中， a′ =. a+b c+d a−b c−d e+ f g+h ， b′ = ， c′ = ， d′ = ， e′ = ， f′= ， 2 2 2 2 2 2. g′ =. e− f g −h i+ j k +l i− j k −l m+n ， h′ = ， i′ = ， j′ = ，k′ = ，l′ = ， m′ = ， 2 2 2 2 2 2 2. n′ =. o+ p m−n o− p ， o′ = ， p′ = 。而第一步驟水平的小波轉換可以把影像分為 2 2 2. 高頻以及低頻兩部分如下圖所示：. 圖 2.8 水平分頻後影像頻率的分布狀況. 其中，L 為低頻的部分而 H 為高頻的部分。之後再將第一步驟所產生的結果做垂直的分頻。則第二步驟垂直的小波轉換結果為： 18.

(19) a′′. b′′. c′′. d ′′. e′′. f ′′. g ′′. h′′. i′′. j′′. k ′′. l ′′. m′′. n′′. o′′. p′′. 圖 2.9 垂直小波轉換後的結果其中，a′′ =. a′ + e′ b′ + f ′ c′ + g ′ d ′ + h′ i′ + m′ j ′ + n′ ，b′′ = ，c′′ = ，d ′′ = ，e′′ = ， f ′′ = ， 2 2 2 2 2 2. g ′′ =. k ′ + o′ l ′ + p′ a′ − e′ b′ − f ′ c′ − g ′ d ′ − h′ ， h′′ = ， i′′ = ， j′′ = ， k ′′ = ， l ′′ = ， 2 2 2 2 2 2. m′′ =. i′ − m′ j ′ − n′ k ′ − o′ l ′ − p′ ， n′′ = ， o′′ = ， p′′ = 。在L頻率中，相加的集合因 2 2 2 2. 為是由低頻中又計算出來的低頻，以LL來代表；而相減的集合是低頻中的高頻，以LH來表示。在H頻率中，相加的集合為高頻中的低頻，以HL來表示；相減的集合是高頻中的高頻，以HH來表示。如下圖所示：. 圖 2.10 垂直分頻後影像頻率的分布狀況離散小波轉換不僅僅只能做一階轉換而已，還能再從第一階轉換結果後的LL子 19.

(20) 頻帶繼續做轉換，出來的結果稱為第二階離散小波轉換，並且還可以繼續以第二階離散小波轉換再從其LL子頻帶做下去，但是通常最多只做到三階而已。在實驗裡，我們就是利用三階的Haar小波離散轉換並且去除最低頻的部分. LLLL頻率來估計我們的p值，下圖為 I max 執行完三次小波轉換後的結果：. 圖2.11 Haar小波轉換後的結果. 其中最左上角的區塊頻率為LLLL頻率，就是要排除的頻域。 2.5.3 最佳化的標準統計相依(Statistical dependency)的程度可以用mutual information. (MI)來測量。 Aˆ lightc 與 Dˆ c 的 MI 可以表示成： I ( Aˆ lightc , Dˆ c ) = H Aˆ. lightc. 其中， H Aˆ. lightc. + H Dˆ − H Aˆ c. (2.22 ). ˆ. lightc , Dc. 與 H Dˆ 為 Aˆ lightc 與 Dˆ c 的marginal entropies。 H Aˆ. ˆ. lightc , Dc. c. 為它們的joint. entropy。如果 Aˆ lightc 與 Dˆ c 越獨立，則它們彼此的mutual information(MI)越小。我 20.

(21) 們再看矩陣W((2.19)式)，它的結構是由標度(scale)p決定。代表 A% lightc 是由 Aˆ lightc 根據p而估計的： A% lightc = I c max − I c min. (2.23). 同樣的代表 D% c = w1 I c max + w2 I c min. (2.24). 也是由 Dˆ c 根據標度 p 而估計的，其中 w1 ≡ ( p − 1) ， w2 ≡ ( p + 1). (2.25). % 1 因此分開 A% lightc 與 D% c 的矩陣為： W= w  1. (2.26). −1 w2 . 為了把 Aˆ lightc 與 Dˆ c 分開，同樣的也是把 A% lightc 與 D% c 分開必須使 A% lightc 與 D% c 的統計相依(statistical dependency)最低。因此必須使 A% lightc 與 D% c 的 mutual information 最小， I ( A% lightc , D% c ) = H A% lightc + H D% c − H A% lightc , D% c. (2.27). % − H max min ，因此((2.27)式)又可以改其中 A% lightc 與 D% c 的 joint entropy 為 log det(W) I ,I c. c. % − H max min 寫成： I ( A% lightc , D% c ) = H A% lightc + H D% c − log det(W) I ,I c. (2.28). c. % 不相依，所以我們把它忽略。而 H % 也為一常數， H I max , I min 為一常數並且與 W Alightc c. c. 同樣的也把它忽略。現在最佳化的問題被簡化為：. min{H D% c − log w2 + w1 } w1, w2. (2.29). 2.5.4 藉由 Probability Model 來簡化本章節主要的目的為利用機率密度函數probability density. function (PDF)來簡化最佳化的問題。其中常被用在影像中的 PDF 模型為 generalized Laplacian： ρ PDF(D% c )=c( ρ ) exp(− D% c ). (2.30). 其中， ρ ∈ (0, 2) 為分布的參數。 c( ρ ) 為一正規化的常數。 D% c 的 entropy 根據[Bof 21.

(22) 01]和[Zib 01]又可以定義為：. {. H D% c = ε − log PDF(D% c ). }. (2.31). 其中 ε 代表期望值。再把((2.30)式)帶入((2.31)式)中，得到 1 Hˆ D% c = C ( ρ ) + N. ∑ D% ( x ). ρ. (2.32). c. x. 其中 N 為影像中像素點的個數， C ( ρ ) = log[C ( ρ )] 。 C ( ρ ) 與 D% c 互相獨立，因此也被我們忽略。最後我們要解的最佳化問題變成 min{− log w2 + w1 + w1, w2. 1 N. ∑ D% ( x) c. ρ. }. (2.33). x. 再用 ρ = 1 來簡化((2.33)式)最後得到： min{− log w2 + w1 + w1, w2. = min{− log w2 + w1 + w1w2. 1 N. ∑ D% ( x) }. 1 N. ∑ wI. c. x max 1 c. ( x ) + w2 I cmin ( x ) }. x. 又， w1 ≡ p − 1 , w2 ≡ p +1 ≡ w1 + 2 所以， 1 N. min{− log 2 w1 + 2 + w1w2. = min{− log 2 p + p. 1 N. ∑ wI. max 1 c. ( x ) + ( w1 + 2) I cmin ( x ) }. x. ∑ ( p − 1) I. max c. ( x ) + ( p + 1) I cmin ( x ) }. (2.34). x. 2.5.5 用平均的方式決定 p 理論上來說，偏振度(degree of polarization)p 應該與不同的小波頻率 c 獨立。在不同頻率 c 所估計出來的 p 理應相同。但實則不然，不同頻率所預估出的 pˆ 可能會不一樣。造成此種現象的原因為：有些頻率遵守我們再 2.5.1 節裡獨立的假設，但有些頻率卻沒有。所以根據實驗的結果，我們對影像作三階 Haar 小波離 22.

(23) 散轉換，並對每個 HH 頻率都估計出 pˆ c ，再以平均的方式求的我們最後的 pˆ 。. 2.6 估計 A 本章節描述我們如何估計另一個重要的參數 A 。首先，先利用[Chu 06]的方法估計出像素銳利的程度。並且選出兩個或以上最銳利並且 radiance 相近的像素點 (xk , yk ) ， k =1,2,....... 。很容易從((2.13)式)((2.14)式)((2.15)式)中得到：. Iˆ total (k ) = J + SAˆ (k ). (2.35). 其中， Aˆ (k ) 可以利用上面章節的方法求得，並且定義 S 如下： S ≡ (1 − J / A). (2.36). 因此，每個被選定的像素點 Iˆ total (k ) 與 Aˆ (k ) 皆為已知。很容易的從((2.35)式)中看出此方程式為一條以 S 為斜率的直線方程式。我們利用之前已選定的兩個或以上的像素點，利用回歸直線的方法描繪出此直線，並且求得 J 與 S。再利用下式： A = J /(1 − S ). 求得我們所需要的 A 值。. 23. (2.37).

(24) 第三章利用單張影像回復影像 3.1 利用偏光鏡回復影像的缺點雖然透過旋轉偏光鏡，所得到的輸入影像 I (x )min 與 I (x ) max 確實可以利用公式的堆導，求得在沒有大氣衰退和 airlight 影響下的景物表面散發出的光 J ( x ) 。但是上述方法卻有下列幾個缺點：第一：所需要的輸入影像為 I (x )min 與 I (x ) max ，顧名思義需要兩張或以上的影像資訊才可以工作。如果我們只得到其中一張影像，並無合適的方法可以推測出另外一張影像，也就是上一章節的方法就不適用。第二：在決定 A(x )min 時，是利用人為的判斷。當轉動偏光鏡使影像亮度最小時就認為 airlight 對亮度的貢獻最少，稱此狀態為 A(x )min 。而 A(x ) max 發生於與. A(x )min 所在的偏光鏡角度相差 90 度的方向上。但是人的判斷多少會因為主觀的認知或者環境的影響而有所誤差，我們並沒有辦法在轉動偏光鏡時制訂一套標準判斷何種相對的亮度為最小、何種為最大。第三，此方法並也不適用於動態的影像。例如當我們在拍攝交通監視系統時，「交通」之所謂為「交通」就是因為車輛會隨著時間的變化而移動。當轉動偏光鏡至某一角度為 A(x )min 時拍攝一張，再轉至另一角度 A(x ) max 時，因為轉動偏光鏡需要些許時間，這時間當中車輛已經移動若干距離造成 I (x )min 與 I (x ) max 兩張影像為不同的場景使 airlight 的估計錯誤。所以車輛移動的區域會因為. airlight 估計的錯誤而不能回復原本真實的影像。如下圖 3.1(a)為 I (x )min，圖 3.1(b) 24.

(25) 為 I (x ) max 很明顯的因為車輛的移動造成兩張影像些許的不同使 airlight 的估計錯誤，造成回復後的影像(圖 3.1(c))紅色圈圈等地方錯誤。. (a). (b). (c) 圖 3.1 動態的影像結果(a)輸入影像 I (x ). min. (b)輸入影像 I (x ). max. (c)回復後的影像. 第四，玻璃或是水面上的反射光通常是偏極光，與我們假設 direct light ((2.1) 式)未被偏極化的前提互相違背。所以如果景物有反射光的部分，例如路面積水、車輛玻璃、湖泊等會造成 airlight 估計的錯誤而無法還原原本的景致。如下圖 3.2(a) 為我們輸入的影像 I (x ) max ，而圖 3.2(b)為回復後的影像。很明顯的因為車輛玻璃 25.

(26) 的反光造成 airlight 估計的錯誤，使我們不能回復原本的影像。. (a). (b). 圖 3.2 影像中有反射光的結果(a)輸入的影像 I (x ) max (b)回復後的影像. 第五，在堆導公式並求得 J ( x ) 的過程中，如果有雜訊的存在，會因為一步步的過程而被放大，造成最後的結果會有較大的雜訊存在。探討上述的缺點發現，之所以有這些問題的存在，不外乎是因為「偏光鏡」、「需要多張影像」所造成的。尤其是應用於戶外的交通監視系統，需要多張影像的輸入尤其不可行。所以我們這章節嘗試著再次利用 ICA 的觀念，致力於探討如何利用「單張輸入影像」、「不需要任何額外的輔助」來回復我們原本的影像。. 3.2 反射率為常數( 的影像反射率為常數 Constant Albedo)的影像影像退化的模式如同 2.1 節所述一模一樣，影像乃是由 direct light ((2.1)式) 和 airlight ((2.7)式) 這兩種光源相加而成，在此就不再贅述。 26.

(27) 景物之所以能被看見，主要是因光打到景物的表面，部分被反射傳到相機，所以每個像點，其值含有一未知的物體反射率向量(R)，其為由 RGB 三個頻道的反射係數所構成。由於不同的像點 x 對應不同的物體表面點，R 應為 x 的函數；但因像點的個數很多，R(x)的估計不易。然而我們可以合理的假設一影像點以及其附近的點所對應的物體表面，具有相同的 R，於是我們可以把 J(x)表為 J ( x ) = l ( x ) R ，其中 l(x)為一純量描述投射在景物表面的光。(2.8)式於是改寫成 I ( x ) = t ( x )l ( x ) R + (1 − t ( x )) A .. (3.1). 我們原本的目的要找 J ( x ) ，即變成要找 l ( x ) 和 R 中的三個反射係數。首先我們把向量 R 拆成兩個分量，第一個分量是 R 與 airlight A 平行的部份. R // = R, A ⋅ A / A ， , 表內積；第二個部份是 R 與 A 垂直的部份 R ⊥ = R − R // 。於是(3.1)式可以寫成：. I ( x ) = t ( x )l ( x ) R + (1 − t ( x )) A = t ( x )l ( x )( R ⊥ + R // ) + (1 − t ( x )) A = t ( x )l ( x ) R ⊥ ( R ⊥ / R ⊥ + R, A / R ⊥ ⋅ A/ A ) + (1 − t ( x )) A , = t ( x )l ′( x )( R ⊥ / R ⊥ + η A/ A ) + (1 − t ( x )) A 其中 l ′( x ) = l ( x ) R ⊥. ， η = R, A / R ⊥ 。. 接下來我們把 I(x)投射到 airlight A 的方向 I A ( x ) = I ( x ), A/ A = t (x )l ′( x )η + (1 − t ( x )) A. 而 I(x)在 A 垂直的方向為. I A⊥ ( x ) =. 2. I (x ) − I A2 ( x ). 27. (3.2).

(28) 2. =. t ( x )l ′( x )( R ⊥ / R ⊥ + η A/ A ) + (1 − t ( x )) A − (t (x )l ′( x )η + (1 − t ( x )) A )2 (t ( x )l ′( x )( R ⊥ / R ⊥ + η A/ A ) + (1 − t ( x )) A) ⋅. = (t ( x )l ′( x )( R ⊥ / R ⊥ + η A/ A ) + (1 − t ( x )) A) − (t (x )l ′( x )η + (1 − t ( x )) A )2 (t 2 ( x )l ′2 ( x )+t 2 ( x )l ′2 ( x )η 2 +2t ( x )l ′( x )η (1 − t ( x )) A + 2. = (1 − t ( x ))2 A ) − (t 2 ( x )l ′2 ( x )η 2 +2t ( x )l ′( x )η (1 − t ( x )) A +. (3.3). 2. (1 − t ( x ))2 A ) = (t 2 ( x )l ′2 ( x ) = t ( x )l ′( x ). 將(3.3)式代入(3.2)式. I A ( x ) = t (x )l ′( x )η + (1 − t ( x )) A = η I A⊥ ( x ) + (1 − t ( x )) A 我們可整理出 t ( x ) = 1 − ( I A ( x ) − η I A⊥ ( x )) / A 。. (3.4). 令集合 Ω 包含我們正在考慮的影像點及其附近的點。假設 airlight A 為已知，因 I(x)為輸入影像，其 I A ( x ) 和 I A⊥ ( x ) 可知。則對於 Ω 中任一點 x，我們有兩條方程式式(3.2)和(3.3)，其中有三個未知數 l ( x ) ， t ( x ) 和η 。以下我們討論如何求η [Hyv 00]，因為知道η 後，由(3.4)式可以得知透射率 t ( x ) ，再將其代入(2.8) 式，即可得到我們想要的 J ( x ) 。由於大氣的透射率 t 是與景物的距離和大氣的衰退係數有關，而 l(x)為投射在景物表面的光，l 和 t 彼此是不相關的。統計學上，稱 l 和 t 是 uncorrelated，它們的共變異數(covariance) CΩ (l , t ) = 0 。共變異數的定義如下. CΩ ( f , g ) = Ω. −1. ∑ (f (x ) − E. Ω. (f ))( g ( x ) − EΩ ( g )) ，. (3.5). x∈Ω. 其中平均值的定義為 EΩ (f ) = Ω. −1. ∑ f (x) 。 x∈Ω. 28. (3.6).

(29) 由式(3.3) (l ′( x ))−1 =. t( x) ，將式(3.4)的 t ( x ) 代入， I A⊥ ( x ). (l ′( x ))−1 = [1 − ( I A ( x ) − η I A⊥ ( x )) / A ] / I A⊥ ( x ) =. 1 − I A ( x ) / A η I A⊥ ( x ) / A 1 − I A ( x ) / A η + = + I A⊥ ( x ) I A⊥ ( x ) I A⊥ ( x ) A. (3.7). A − I A ( x) η + I A⊥ ( x ) A A. =. 令 h( x ) = ( A − I A ( x )) / I A⊥ ( x ) ，則. (l ′( x ))−1 =. A − I A ( x) η + = ( h( x ) + η ) / A I A⊥ ( x ) A A −1. (3.8)前述因 l 和 t. 不相關，l ′( x ) = l ( x ) R ⊥ ，(l ′( x )) 與 t ( x ) 也是不相關；再由(3.8)式， h( x ) 與 t ( x ) 也不相關。由(3.2)和(3.3)式， I A ( x ) = η I A⊥ ( x ) + (1 − t ( x )) A 。 CΩ ( I A ( x ), h( x )) = Ω. −1. ∑ (I. A. ( x ) − EΩ (I A ( x )))(h( x ) − EΩ (h( x ))). x∈Ω. =Ω. −1. ∑ (η I x∈Ω. A⊥. ( x ) + (1 − t ( x )) A −. EΩ (η I A⊥ ( x ) + (1 − t ( x )) A ))(h( x ) − EΩ (h( x ))) =Ω. −1. ∑ (η I x∈Ω. A⊥. ( x ) − η EΩ (I A⊥ ( x )) + (1 − t ( x )) A − EΩ ((1 − t ( x )) A ))(h( x ) − EΩ (h( x ))). =η Ω. −1. ∑ (I x∈Ω. Ω. −1. A⊥. ( x ) − EΩ (I A⊥ ( x )))(h( x ) − EΩ (h( x ))) +. ∑ (1 − t ( x) − E. Ω. (1 − t ( x ))) A (h( x ) − EΩ (h( x ))). x∈Ω. 29.

(30) Q Ω. −1. ∑ (1 − t ( x ) − E. Ω. (1 − t ( x ))) A (h( x ) − EΩ (h( x ))). x∈Ω. = A Ω. −1. ∑ (1 − t ( x) − 1 + E. Ω. (t ( x )))(h( x ) − EΩ (h( x ))). x∈Ω. = A Ω. −1. ∑ (E. Ω. (t ( x )) − t ( x ))(h( x ) − EΩ (h( x ))). x∈Ω. = A Ω. −1. ∑ (E. Ω. (t ( x ))h( x ) − t ( x )h( x ) −. x∈Ω. EΩ (t ( x )) EΩ (h( x )) + t ( x ) EΩ (h( x ))) = A (Ω. −1. ∑E. Ω. (t ( x ))h( x ) − Ω. x∈Ω. −1. ∑ t ( x ) h( x ) − x∈Ω. Ω. −1. ∑E. Ω. (t ( x )) EΩ (h( x )) + Ω. x∈Ω. −1. ∑ t ( x) E. Ω. (h( x ))). x∈Ω. = A ( EΩ (t ( x )) EΩ (h( x )) − EΩ (t ( x )h( x )) − EΩ (t ( x )) EΩ (h( x )) + EΩ (t ( x )) EΩ (h( x )) = A (− EΩ (t ( x )h( x )) + EΩ (t ( x )) EΩ (h( x )). 因為 h( x ) 與 t ( x ) 不相關， EΩ (t ( x )h( x )) = EΩ (t ( x )) EΩ (h( x )) ，上式為 0，因此 CΩ ( I A ( x ), h( x )) = η CΩ ( I A⊥ ( x ), h( x )) + 0 = η CΩ ( I A⊥ ( x ), h( x )) ，得. η=. CΩ ( I A ( x ), h) CΩ ( I A⊥ ( x ), h). (3.9). 上式要成立，須分母 CΩ ( I A ( x ), h) 不為零。我們將在下一章節探討如何確保分母不為零。. 3.3 反射率不為常數的影像前面我們假設整張影像的像點，皆對應同一性質的物體表面，故假設整張影像所含物體的反射率為一常數，亦即R(x) = R，這也意味著我們算出來的η 也是一個常數。但是一般的影像包含多種性質的景物表面，其反射率並不是一個常數。但一般來說，一點與其周圍的像點，其對應物體的反射率不會相差太多。我們利用此特性，首先計算haze-free像點 I ′ = I − AI A / A (有可能是負值)，再把它 30.

(31) 由RGB轉換到YUV色彩空間。令. θ ( x ) = tan −1 (U ( x ) / V ( x ))，. (3.10). 其中 U = −0.14713* R − 0.28886* G + 0.436* B V = 0.615* R − 0.51499* G − 0.10001* B 。. 上述公式把像素點的色彩以角度值來表示。而以. d (θ1 ,θ 2 ) = min{ θ1 − θ 2 , 2π − θ1 − θ 2 }. (3.11). 來表示兩彩色間的差異。定義一權重共變異數(covariance)[Andrews et al. 1972]為. Cx ( f , g) = 其中 E x ( f ) = Wx =. 1 Wx. 1 Wx. ∑. ∑ ( f ( y) − E. x. ( f ))( g ( y ) − E x ( g )) w( x , y ). (3.12). y∈Ω x. f ( y ) w( x , y ) ， w( x, y ) = exp(−d (θ ( x ), θ ( y )) 2 / σ θ2 ) ，. y∈Ω x. ∑ w( x, y ) 。 Ω. x. 代表以像素點x為中心，加上其周圍鄰近的像素點，再去除. y∈Ω x. 集合 β 中的點後，所成的集合。在我們的實驗中設 σ θ = 1/8，並用24 × 24來當作視窗的大小。視窗的尺寸可以視雜訊的嚴重程度來做放大或者是縮小。最後，再把(13)(14)式中的 EΩ 和 CΩ 替換成上面的 C x 和 E x 。對於每一像點 x ∉ β ，我們計算其η ( x ) ，再由(3.4)式可估計其透射率 tˆ( x ) ，再將其代入(2.8)式，即可得到我們想要的 J ( x ) 。前述所估算的 tˆ( x ) 難免會含有誤差；又在估算像點的透射率時，我們忽略了集合 β 中的像點。為了彌補這些缺失，我們對 tˆ( x ) 作一平滑(smoothing)過程。由於透射率與景物的距離和大氣的衰退係數有關，所以透射率具有空間上的規律 31.

(32) 性。利用此特性，我們以一Gauss Markov random field (GMRF) model [Pér 98]來平滑 tˆ( x ) 。由(GMRF) model，我們得一客觀方程式 ˆ. P(t ( x )) ∝ ∏ e − (t ( x ) −t ( x )). 2. / σ t2 ( x ). x∉β. ∏e. − ( t ( x ) −tˆ ( y )) 2 /( I A ( x ) − I A ( y )) 2 / σ s2. (3.13). y∈N x. 其中集合 N x 含像素點 x 的四個 first-order 鄰近點， σ s = 1/ 5 。我們的目的即在尋找一 t ( x ) 使 P(t ( x )) 最大化。首先將上式取對數 log. log P (t ) = ∑ −(t ( x ) − t$( x )) 2 / σ t2 + x∉β. ∑ −(t ( x ) − tˆ( y)). 2. /( I A ( x ) − I A ( y )) 2 / σ s2. y∈N x. 將上式對 t ( x ) 微分並令它等於 0，. −2(t ( x ) − $t ( x )) / σ t2 +. ∑ −2(t ( x ) − tˆ( y )) /( I. A. ( x ) − I A ( y )) 2 / σ s2 = 0 .. y∈N x. i, 當 x ∉ β ： ⇒ −2(t ( x ) − tˆ( x )) / σ t2 + [−2(t ( x ) − tˆ( yup )) /( I A ( x ) − I A ( yup )) 2 / σ s2 ] + [−2(t ( x ) − tˆ( ydown )) /( I A ( x ) − I A ( ydown ))2 / σ s2 ] + [−2(t ( x ) − tˆ( yright )) /( I A ( x ) − I A ( yright )) 2 / σ s2 ] + [−2(t ( x ) − tˆ( yleft )) /( I A ( x ) − I A ( yleft ))2 / σ s2 ] = 0. 令 σ t2 = a . ⇒. ( I A ( x ) − I A ( yup )) 2 / σ s2 = b,. ( I A ( x ) − I A ( ydown )) 2 / σ s2 = c,. ( I A ( x ) − I A ( yleft ))2 / σ s2 = d , ( I A ( x ) − I A ( yright )) 2 / σ s2 = e.. (3.14)式可寫成 −2(t ( x ) − $t ( x )) / a + [−2(t ( x ) − tˆ( yup )) / b] + [−2(t ( x ) − tˆ( ydown )) / c] + [−2(t ( x ) − tˆ( yright )) / d ] + [−2(t ( x ) − tˆ( yleft )) / e] = 0. 32. (3.14).

(33) ⇒ [−2bcde(t ( x ) − t$( x ))] + [−2acde(t ( x ) − t ( yup ))] + [−2abde(t ( x ) − t ( ydown ))] + [−2abce(t ( x ) − t ( yright ))] + [−2abcd (t ( x ) − t ( yleft ))] = 0. ⇒ −2bcdet ( x ) + 2bcdetˆ(x )-2acdet ( x ) + 2acdetˆ( yup ) − 2abdet ( x ) +2abdetˆ(ydown )-2abcet ( x ) + 2abcetˆ( yright ) − 2abcdt ( x ) +2abcdtˆ(yleft ) ⇒ 2bcdet$( x ) + 2acdetˆ( yup ) + 2abdetˆ( ydown ) + 2abcetˆ( yright ) + 2abcdtˆ( yleft ) = 2bcdet ( x ) + 2acdet ( x ) + 2abdet ( x ) + 2abcet ( x ) + 2abcdt ( x ) = 2t ( x )(bcde + acde + abde + abce + abcd ). t ( x) =. bcdet$( x ) + acdetˆ( yup ) + abdetˆ( ydown ) + abcetˆ( yright ) + abcdtˆ( yleft ) (bcde + acde + abde + abce + abcd ). ii, 當 x ∈ β ：. ∑ −2(t ( x ) − tˆ( y)) /( I. A. ( x ) − I A ( y )) 2 / σ s2 = 0. y∈N x. ⇒ [−2(t ( x ) − tˆ( yup )) /( I A ( x ) − I A ( yup ))2 / σ s2 ] + [−2(t ( x ) − tˆ( ydown )) /( I A ( x ) − I A ( ydown ))2 / σ s2 ] + [−2(t ( x ) − tˆ( yright )) /( I A ( x ) − I A ( yright )) 2 / σ s2 ] + [−2(t ( x ) − tˆ( yleft )) /( I A ( x ) − I A ( yleft ))2 / σ s2 ] = 0. ⇒ [−2(t ( x ) − tˆ( yup )) / b] + [−2(t ( x ) − tˆ( ydown )) / c] + [−2(t ( x ) − tˆ( yright )) / d ] + [−2(t ( x ) − tˆ( yleft )) / e] = 0. ⇒ [−2cde(t ( x ) − t ( yup ))] + [−2bde(t ( x ) − t ( ydown ))] + [−2bce(t ( x ) − t ( yright ))] + [−2bcd (t ( x ) − t ( yleft ))] = 0 ⇒ 2dcet ( yup ) + 2dbet ( ydown ) + 2bcet ( yright ) + 2bcdt ( yleft ) = 2dcet ( x ) + 2dbet ( x ) + 2bcet ( x ) + 2bcdt ( x ) = 2t ( x )(cde + bde + bce + bcd ). ⇒ t( x) =. dcet ( yup ) + dbet ( ydown ) + bcet ( yright ) + bcdt ( yleft ) cde + bde + bce + bcd. 33.

(34) 3.4 估計Airlight 估計 Color A 本節討論如何估計airlight color vector A。先前我們曾論述，在局部範圍內，景物表面光的強弱度l和大氣透射率t，兩者是不相關的。根據此觀念，我們要來尋找一個A，使l和t之間的相關性最小。考慮一小視窗(在實驗裡使用24 × 24的視窗，並且選定此視窗為霧最濃的區域)，假設在此範圍內景物表面的反射率是常數(constant albedo)R。我們使用steepest descent method來估計A；一開始猜測一個. A，然後不斷的更新A直到 C ( J , t ) 2 達到最小時的A即為所求。上述的方法，因為是在三維度的空間中搜尋，在計算上非常複雜且耗時。為了降低時間與計算的複雜度，我們嘗試降低維度。根據第(9)式，在視窗內所有像素點x的I(x)皆在R與A所展成的空間，這是一個二維度的空間。首先在影像中找兩個或以上的小區域 R1 和 R2 ...，每個區域內的反射率R為相同，但不同的區域有不同的反射率。對每一區域 Ri , i = 1, 2 ,…我們將其像點的RGB值表在RGB空間，然後應用PCA (principal component analysis)取共變異(covariance)最活躍的兩個軸 v1i 和 v2i 。首先計算像點的 covariance matrix ， v1i 和 v2i 為此矩陣兩個最大的. eigenvalues所對應的eigenvectors。令 Pi 為由 v1i 和 v2i 所展成的平面。由於A必需同時存在於平面 P1 和 P2 ，A的方向可由此兩平面的交集決定。上述的問題可以轉換成尋找一向量使其投影在 P1 和 P2 的分量和為最大，亦即. max ∑ A, v1i A. 2. 2. + A, v2i ,. 2. A =1. i. 又根據Lagrange-multipliers rule，上式的解其實在求3 × 3矩陣， 34. (3.15).

(35) ∑ v (v ) i i 1. i T 1. + v2i (v2i )T ，最大的eigenvalue所對應的eigenvector。在決定了A的方向之. 後，我們即可使用前述steepest descent method來估計 A ；一開始猜測一個A，然後不斷的更新A直到 C ( J , t ) 2 達到最小。. 3.5 區域圖本章節的目的是為了把輸入的影像 I ，性質相近的區塊劃分成不重疊的區域而形成區域圖 S ，為了自動選擇在3.4節裡估計A所需的小區域 Ri , i = 1, 2 ,…。其最主要的概念是採用分割和合併(divide-and-merge)的方法來切割出我們的區域圖S。我們利用像素點之間的相似程度，來幫助我們切割出更好的區域圖。首先，我們把輸入的影像I沿著水平和垂直的方向平均分割成四個子圖 ri , i = 1, 2,3, 4 ，並且用 Θ i 分別代表四個子圖的特徵向量。對於每個子圖 ri ，我們令. Θi = ( ACi , ATi ) ，其中 ACi 表每個子圖 ri 中像素點的平均色度特徵，而 ATi 表每個子圖 ri 中像素點的平均紋理特徵。不同的應用，對於色度特徵有不同的定義。我們利用不變色度特徵(invariant. chromatic characteristics)來定義我們的色度特徵 {H , C }。其中 H = Eλ / Eλλ 表物質的色彩。E表反射的能量，且其下標的符號代表不同的波長 λ 所對應的E。E (λ , x ) 表示對於像素點 x 所反射的能量，定義如下：. E (λ , x ) = i (λ , x )[(1 − ρ ( x )) 2 R(λ , x ) + ρ ( x )]. (3.16). 其中， x 表每個影像所對應的像素點， λ 表波長， i (λ , x ) 為影像的亮度， ρ ( x ) 為 35.

(36) Fresnel表面反射係數， R(λ , x ) 表物體的反射係數。而 C = Eλ / E 為正規化的顏色。 E 、 Eλ 和 Eλλ 可以經由RGB的線性組合來估計，.  E   0.06 0.63 0.27   R   E  =  0.3 0.04 −0.35 G   λ     Eλλ   0.34 −0.6 0.17   B . (3.17). 接下來，我們利用steerable filter來估計紋理特徵。此方法主要是利用steerable 的兩個filter G2 和 H 2 來估計紋理特徵。其中 G2 為圓對稱高斯(circularly symmetric. Gaussian)的二次微分，H 2 為 G2 的Hilbert轉換。而 G2 和 H 2 被旋轉的情形(以 Gθ 2 和 H θ 2 表示)可以利用有限的基本濾波( basis filers)的線性組合來逼近。對於每個像素 p方向的能量定義如下： E2θ ( p ) = [G2θ ( p )]2 + [ H 2θ ( p )]2. (3.18). 我們可以利用phase angle φ ( p ) 在某方向所具有最大的方向能量來定義像素p的紋理特徵： φ ( p ) = arg max E2θ ( p ). (3.19). θ. 所以，特徵向量就可以定義為： Θ = ( AC , AT ) = ( H / ρ H , C / ρC , φ / ρφ ). (3.20). 其中 ρ 為權重常數。最後，兩個子圖 ri 與 rj 之間的相似程度就可以定義為：. ζ (i, j ) = Θι Θ j / Θι Θ j. (3.21). 我們所採用的方法為分割和合併(divide-and-merge)，顧名思義第一步驟必須對原圖分割。當我們所分割的四個子圖之間的相似程度夠接近就停止分割。不然就分別對上述四個子圖再繼續分割成四等份，並且再判斷其相似程度。如此遞迴下去直到所分割的四個子圖相似程度皆大於某個值才停止。如下圖3.6為原圖而 36.

(37) 圖3.7為分割後的情形。第二步驟為合併，首先先檢查子圖附近之子圖的相似程度。如果相似程度夠大，就把兩子圖合併，如此遞迴下去直到無法再合併為止。合併後的子圖必須再重新計算特徵向量。下圖3.3(c)為計算紋理所需要的邊緣圖，圖3.3(d)為最終的區域圖。. (a). (b). (c). (d). 圖3.3區域圖之步驟(a)輸入影像(b)分割後的影像c)邊緣圖d)區域圖. 3.6 雜訊的估計任何的測量都有雜訊的存在，它在我們除霧的過程中造成不小的影響。雜訊 37.

(38) 會使估計的η 值不準確，而 η 值在我們計算影像 J 時扮演著重要的角色。在這裡我們假設雜訊的分佈是白色高斯(white Gaussian)。令 I = I% + ζ 。其中 I% 表. noise-free 影像， ζ 是高斯雜訊，其平均值為零，變異數是 σ 2 。雜訊會經由一系列的計算而放大甚至擴散；從輸入的影像開始，經由 η 的計算，再到值 t 的計算。最後得到影像 J。這些計算基本上都是線性的(linearity)，我們因此可以假設存在最後影像 J 的雜訊，其分佈還是高斯。為了使雜訊得以順利估計，我們列出一些雜訊必須滿足的條件，並且解釋為什麼。首先每個彩色頻道 I i ， i = r , g , b ,其雜訊的變異數 σ 2i 必須小於 I i /10 ，也要小於 I A⊥ / 6，為要滿足 signal-to-noise 的要求，確保雜訊不會過大而影響了輸出影像的品質。又雜訊 ζ c R ′ 的變異數，必須要小於 (CΩ ( I R′ , h) / 5) 2 。其中， ζ c R ′ 是出現在 CΩ ( I R′ , h) 內的雜訊。這裡要注意一點： Ω 不能太小。如果 Ω 太小，我們沒有足夠的像素點來平均，使 ζ c R ′ 的變異數降低而滿足我們的要求，必定違反第二種情況。直覺上，訂定這種情況是為了確保 E Ω 和 CΩ 在統計上的誤差不會太大。根據上面的討論，我們並不假設這些情況一定會被滿足。相反的，因為那些違反上述情況的像素點，其所隱藏的資訊並不可靠，所以在做計算時，我們把這些點都忽略而不予理會。那些被忽略的像素點集合為 β 。當最後在求 t 值時，其實是針對 x ∉ β 的像素點來做計算。輸入一張影像時，必須估計出適當的 σ 來描述雜訊分布的情形。我們把這些留給使用者自行去設定。在我們的實驗中， σ 維持在1/100 – 1/200之間。最近， 38.

(39) [Liu 06]提出了一個利用單張的影像，估計出雜訊分布的方法。. 39.

(40) 第四章實驗結果實驗的影像分為兩種類型，一種是同一場景拍攝單張影像，相機採用的是. Panasonic DMC-FZ50GT 的一般鏡頭拍攝；另外一種是同一場景拍攝兩張影像，在一般鏡頭前面加裝偏光鏡並旋轉至某角度拍攝。而偏光鏡採用的是 B+W. CIRCULAR-POL。圖 4.1(a)為相機裝置偏光鏡前，而圖 4.1(b)為相機裝置偏光鏡前後，圖 4.3(c)為偏光鏡。. (a). (b). (C) 圖 4.1 實驗器材 40.

(41) 實驗的場景選擇了在台灣常見的天氣型態：起霧和下雨。並且針對各種的天氣型態拍攝了交通影像與非交通影像。再針對交通影像拍攝了一般道路與高速公路兩種類型。分類如下圖所示：. 圖 4.2 場景分類示意圖影像能見度的估計是參考[Hec 02]和[Sch 05]利用對比程度來定義如下：. C(I ) =. 1 [ I c ( x, y ) − I c ]2 ∑ ∑ N ℵ= R ,G , B x, y. ∑. (4.1). Ic. c = R ,G , B. 其中， c 代表彩色的頻道，N 表影像的像素點個數， I c = (1/ N )∑ I c ( x, y ) 代表 x, y. 各頻道亮度的平均值。我們利用此公式來計算處理前與處理後影像的能見度，進而輔助我們比較能見度的差異。. 4.1 偏光鏡為基礎 41.

(42) 本實驗在這部分所使用的影像為鏡頭前加裝偏光鏡並旋轉至若干角度所拍攝。影像原本的解析度為 1600X1200，在經過 sub-sampling 程序後，所處理的影像解析度為 320X240。使用 Pentium Duo Core2 為中央處理器的桌上型電腦作影像處理，處理時間不到 1 秒。首先，如圖 4.3 為非交通影像且起霧的時候，其中圖 4.3(a) 為 I (x )max ，圖. 4.3(b)為 I (x )min ，靠近山區的建築物能見度並不佳。為了求得參數 p 各對兩張影像作 Haar 小波離散轉換，如圖 4.3(c)為 4.3(a)作完三階小波離散轉換後的結果。利用 ICA 求得 p 與演算法求得 A 後就可以利用以上的資訊回復景物原本的景致，如圖 4.3(d)。. (a). (b). (c). (d). 圖 4.3 非交通起霧場景(a) I (x )max (b) I (x )min (c)三階 Haar 小波轉換(d)回復的影像 42.

(43) 其中，下圖 4.4(a)為上圖 4.3(b)在遠處建築物附近紅色框框的放大圖。而圖 4.4(b) 為圖 4.3(d) 在遠處建築物附近紅色框框的放大圖，明顯的 4.4(b)清楚很多。但因為 4.4(b)天空的部分全部變為白色造成對比程度下降，故此例子並無法用上述能見度的估計方法來計算其能見度。. (a). (b). 圖 4.4 放大影像(a)為圖 4.3(b)之放大圖(b)為圖 4.3d 之放大圖如下圖 4.5 在非交通且下雨的情況下，遠處的建築物因為雨的散射使得對比度減少，在處理前的能見度 C ( I ) ≈ 0.27。處理過後如圖 4.5(c)能見度 C ( I ) ≈ 0.30 。. (a). (b). 圖 4.5 非交通且下雨場景(a) I (x )max (b)回復的影像 43.

(44) 在一般道路起霧的場景中，為了模擬交通監視系統通常架設在路口至高點處的情形，我們把相機設置在天橋上並且由上往下的拍攝方式來取得一般道路起霧的交通場景。也為了使起霧對視覺的影響更為顯著，我們通常利用伸縮鏡頭，把場景拉至較遠處。如下圖 4.6 為一般道路中起霧的情形，處理前的能見度 C ( I ) ≈ 0.23 。圖 4.6(c)為處理後的影像。整個影像的能見度明顯提高，處理後的. 能見度 C ( I ) ≈ 0.32 。而圖 4.6(d)為 4.6(a)在遠方起霧比較明顯之處的放大圖。. (a). (b). (c). (d). 圖 4.6 一般道路起霧場景(a) I (x )max (b) I (x )min (c)回復的影像(d)為 4.6(a)紅色區塊放大圖 44.

(45) 在一般道路下雨的場景中，相機同樣設置在天橋上並且由上往下的拍攝方式來取得一般道路下雨的交通場景。值得注意的是，下雨時路面會出現積水的情形如圖 4.7(a)路緣的部分，積水會反射光線並且伴隨偏極化的情形。使得我們原本認定景物原本散發出的光沒有被偏極化的假設因為積水的反射光而不成立，造成. airlight 估計的錯誤。如圖 4.7(c)路緣的部分因 airlight 估計的錯誤而錯誤。圖 4.7(a) 的能見度 C ( I ) ≈ 0.30 ，圖 4.7(c)處理後影像的能見度 C ( I ) ≈ 0.34 。. (a). (b). (c) 圖 4.7 一般道路下雨場景(a) I (x )max (b) I (x )min (c)回復的影像. 45.

(46) 在高速公路起霧的場景中，相機通常架設在路肩或者架設在高速公路旁的小山坡上。為了使起霧對視覺的影響更為顯著，如同我們在拍攝一般道路時，把場景拉至較遠處。下圖 4.8(a)為架設在路肩時起霧的情景( I (x )max )，處理前的能見度 C ( I ) ≈ 0.26 。4.8(b)為處理過後的影像，能見度 C ( I ) ≈ 0.30 。圖 4.9(a)為架設在路旁小山坡的情景( I (x )max )，能見度 C ( I ) ≈ 0.26 。圖 4.9(b)為處理過後的影像，能見度 C ( I ) ≈ 0.34 ，確實比處理前好。. (a). (b). 圖 4.8 高速公路起霧場景一(a)輸入影像(b) 回復的影像. (a). (b). 圖 4.9 高速公路起霧場景二(a)輸入影像(b) 回復的影像 46.

(47) 在高速公路下雨的場景中，相機架設在高速公路旁的小山坡上。與在一般道路的情形一樣，雖然高速公路積水的機率比較小，但如果發生會就會造成 airlight 估計的錯誤使得回復的影像錯誤。如 4.10(a)與 4.11(a)皆為處理前高速公路下雨時的景象( I (x )max )。能見度各為 C ( I ) ≈ 0.30 和 C ( I ) ≈ 0.45 。4.10(b)與 4.11(b)為處理後的影像，前者能見度 C ( I ) ≈ 0.31 ，因為天空的部分皆變成白色，造成對比程度的下降使整張影像能見度並沒有明顯提升。後者也因為此原因造成能見度與處理前的影像一樣。但遠處的建築物附近可以明顯看出能見度比處理前為佳。. (a). (b). 圖 4.10 高速公路下雨場景一(a)輸入影像 (b) 回復的影像. (a). (b). 圖 4.11 高速公路下雨場景二(a)輸入影像(b) 回復的影像 47.

(48) 4.2 單張影像為基礎單張影像為基礎本實驗在這部分所使用的影像為單張影像。與前一章節一樣，影像原本的解析度為 1600X1200，在經過 sub-sampling 程序後，所處理的影像解析度為. 320X240。使用 Pentium Duo Core2 為中央處理器的桌上型電腦作影像處理，處理時間約 20 秒。實驗場景如圖 4.2 所分類，我們針對各種情況皆作測試。首先如圖 4.12(a)為非交通影像且下雨時的場景。因為水氣的散射，使得遠方建築物的能見度明顯降低，此時能見度 C ( I ) ≈ 0.20 。為了估計雜訊和 A 我們先對原影像求區域圖如 4.12(b)所示，區域圖中相同顏色的區塊我們認為是同一性質的區域。求得 A 後就可以利用上一章節的推導求得景物原本的影像 J(x) 如圖 4.12(c)所示，此時能見度 C ( I ) ≈ 0.36，提升約 1.8 倍。而圖 4.12(d)為圖 4.12(a) 在紅色方框的放大圖；圖 4.1(e)為圖 4.12(c) 在紅色方框的放大圖。由此比較後更可以看出處理後的景物能見度明顯改善。. (a). (b). 48.

(49) (b). (d). (e) 圖 4.12 非交通影像下雨場景(a)輸入影像(b)區域圖(c)處理結果(d)處理前遠方放大圖(e)處理後遠方放大圖下圖 4.13(a)為非交通且起霧的場景，處理前的能見度 C ( I ) ≈ 0.28 ，處理後的能見度 C ( I ) ≈ 0.35 。能見度提升約 1.25 倍. (a). (b). 圖 4.13 非交通起霧場景(a)輸入影像(b)處理後的影像 49.

(50) 在一般道路起霧的場景中，因為輸入的影像為單張影像，並不會同上一章節所出現的問題一樣，受到車輛移動的影響造成恢復影像的錯誤。如圖 4.14(a)、. 4.15(a)為一般道路起霧的影像，能見度各為 C ( I ) ≈ 0.22 和 C ( I ) ≈ 0.21 。4.14(b)和 4.15(b)為處理後的影像，改善後的能見度各為 C ( I ) ≈ 0.29 和 C ( I ) ≈ 0.27 ，個別提升了約 1.3 倍和 1.28 倍。. (a). (b). 圖 4.14 一般道路起霧場景一(a)輸入影像(b)處理後的影像. (a). (b). 圖 4.15 一般道路起霧場景二(a)輸入影像(b)處理後的影像 50.

(51) 在一般道路下雨的場景中，此方法並不會受到路面積水的影響而造成錯誤。如圖 4.16(a)、4.17(a)為一般道路下雨的影像，能見度各為 C ( I ) ≈ 0.19 和 C ( I ) ≈ 0.30 。4.16(b)和 4.17(b)為處理後的影像，改善後的能見度各為 C ( I ) ≈ 0.25. 和 C ( I ) ≈ 0.39 ，個別提升了約 1.31 倍和 1.3 倍。. (a). (b). 圖 4.16 一般道路下雨場景一(a)輸入影像(b)處理後的影像. (a). (b). 圖 4.17 一般道路下雨場景二(a)輸入影像(b)處理後的影像. 51.

(52) 在高速公路起霧時的場景，如下圖 4.18(a)為在高速公路旁小山坡上拍攝的影像，此時能見度 C ( I ) ≈ 0.21 。圖 4.18(b)為處理過後的影像，此時能見度 C ( I ) ≈ 0.34，提升約 1.61 倍。圖 4.19(a)為在路肩拍攝的影像，能見度 C ( I ) ≈ 0.26。. 圖 4.19(b)為處理過後的影像，能見度 C ( I ) ≈ 0.33 ，提升約 1.26 倍。. (a). (b). 圖 4.18 高速公路起霧場景一(a)輸入影像(b)處理後的影像. (a). (b). 圖 4.19 高速公路起霧場景二(a)輸入影像(b)處理後的影像. 52.

(53) 在高速公路下雨時的場景，如下圖 4.20(a)為在高速公路旁小山坡上拍攝的影像，此時能見度 C ( I ) ≈ 0.22 。圖 4.20(b)為處理過後的影像，此時能見度 C ( I ) ≈ 0.43 ，提升約 1.95 倍。圖 4.21(a)為另一在高速公路旁小山坡上拍攝的影. 像，能見度 C ( I ) ≈ 0.36 。圖 4.21(b)為處理過後的影像，能見度 C ( I ) ≈ 0.45 ，提升約 1.25 倍。. (a). (b). 圖 4.20 高速公路下雨場景一(a)輸入影像(b)處理後的影像. (a). (b). 圖 4.21 高速公路下雨場景二(a)輸入影像(b)處理後的影像 53.

(54) 4.3 交通影像之能見度統計本實驗在此部分為針對各種交通場景統計其影像改善前後能見度的差異，並以折線圖表示。如下圖為一般道路下雨時能見度統計圖，橫軸為影像編號，在此場景總共張數有 13 張影像。縱軸為能見度。其中，藍色實線為處理前各個影像之能見度，而紅色實線則為處理後各個影像之能見度。我們以同樣影像來做比較，明顯的處理後的影像能見度皆有一定程度的改善。藍色虛線則為此場景所有處理前影像能見度之平均，而紅色虛線則為此場景所有處理後影像能見度之平均。由下表得知影像處理後能見度之平均遠大於影像處理前能見度的平均。. 圖 4.22 一般道路下雨時能見度統計圖. 54.

(55) 下圖則為一般道路起霧時能見度統計圖，此場景總共用 13 張影像來統計：. 圖 4.23 一般道路起霧時能見度統計圖. 下圖為高速公路下雨時能見度統計圖，此場景總共用了 11 張影像來統計：. 圖 4.24 高速公路下雨時能見度統計圖. 55.

(56) 下圖為高速公路起霧時能見度統計圖，此場景總共用了 14 張影像來統計：. 圖 4.25 高速公路起霧時能見度統計圖. 由上述統計資料得知，處理後能見度皆有一定程度的提升。. 56.

(57) 第五章結論與未來工作 5.1 結論近年來隨著科技的進步，數位攝影器材日趨成熟，使得產品價格越來越能被大多數人所接受。電腦的普及使得這些數位攝影器材除了操作簡易外，也很方便我們連接電腦作進一步的後續處理。如此豐富且方便的影像資訊，也造就了許多系統利用影像來協助我們做各式各樣的應用。但許多戶外場景的應用常假設在天氣狀況良好時才有較高的辨識率或正確率，對於各種不同類型的天氣影像容忍度卻非常低。因此本研究提出了只需要單張輸入影像就可以去除我們常見的天氣型態所造成的影響。使各種系統對不同的天氣型態容忍度提高。恢復景物原本的影像，無可避免的可以分為兩大問題。分別為估計 airlight 與估計透射率 t 。本研究藉由兩種不同的方法分別得到 airlight 的值與透射率 t 值，再藉由影像退化模式的公式推估出景物在未受干擾前最原始的亮度。第一種方法先假設 airlight 為偏極化的光，再利用偏光鏡所獲得的資訊與 ICA 的觀念求出我們所需的參數進而估計出 airlight。此時 airlight 的值裡面包含著透射率 t 的資訊，利用公式的推導求得透射率 t 值。此種方法雖可以成功的回復景物原本的影像，但還是有需要兩張以及不適用於交通影像等缺點。因此我們嘗試著用另一種方法只需要單張影像輸入並改善此種方法的缺點。第二種方法我們先把景物散發出的光拆解成反射係數 R 與投射在景物表面的光 l 兩部分。再利用 ICA 的觀念─大氣的透射率 t 是與景物的距離和大氣的衰 57.

(58) 退係數有關，而 l 為投射在景物表面的光，l 和 t 彼此是不相關的。利用此觀念我們分別求出 t 與 A 值，再利用 t 與 A 值求得 airlight。此種方法成功的使能見度增加，並且因為只需要單張輸入影像，也不需為任何其他額外的資訊與人為的操作，所以適用於交通等隨時在變動的場景。. 5.2 未來方向由上述可知，A 的估計對於此方法是最為重要的。但我們在估計 A 時僅對霧最濃的 24X24 區塊作計算，造成如果此區塊變化很大時所估計出的 A 並不盡正確。在未來的工作當中我們需要使估計出的 A 更加的準確。雖然本研究能確實改善各種天氣影響下影像的能見度，但計算的時間約需. 30 秒才能完成，目前還無法做到 real-time 處理。未來必須想辦法改進計算時間，使其可以做到 real-time。如此一來就並不只可以裝置在路口的監視系統，甚至可以裝置在飛機車輛上增加飛行時或行車的安全。使影像能見度的增加並不只可以應用在不同天氣狀態之下。就我們所知，在水底拍攝的影像也是因為受到水分子散射的影響使能見度變差，如圖 5.1 所示。但水底與空氣中不同的是，在水裡因為景物的光透過水傳遞使得整個色調偏藍色。我們未來可以針對在水底的應用上使影像恢復原本的面貌。. 58.

(59) 圖 5.1 水底所拍攝的影像以偏光鏡為基礎的影像恢復之所以不適用於交通影像的原因為，需要兩張影像的輸入且因調整偏光鏡需要時間所已使兩張輸入的影像為不同的景致。未來如果設置一裝置能夠自動且同時擷取兩張影像，則上述問題將可解決。. 59.

(60) 參考文獻: 參考文獻: [Bof 01]P. Bofill and M. Zibulevsky. “Underdetermined blind source separation using sparse representations.” Signal Processing, 81:2353–62, 2001. [ Coz 97]F. Cozman and E. Kroktov. “Depth from scattering.” In Proc.IEEE CVPR, pages 801–806, 1997. [Cha 88]P. S. Chavez. 1988.” An improved dark-object subtraction technique for atmonspheric scattering correction of multispectral data.” Remote Sensing of Environment 24, 450–479. [Chu 06] Y. C. Chung, S. L. Chang, J. M. Wang, and S. W. Chen, “An Edge Analysis Based Blur Measure for Image Processing Applications” National Taiwan Normal Universal. [Fat 07]R. Fattal. 2007. “Image upsampling via imposed edge statistics.” ACM SIGGRAPH 26, 3, 95. [Fat 08]R. Fattal.” Single image dehazing.” to appear in SIGGRAPH,2008. [Hyv 00]A. Hyvrinen, and E. Oja. 2000. “Independent component analysis: Algorithms and applications.” Neural Networks 13, 411–430. [Hec 02]E. Hecht, Optics, pp. 562-628. fourth ed., Addison Wesley, 2002. [Liu 06]C. Liu, W. T. Freeman, R. Szeliski, and S. B. Kang.2006. “Noise estimation from a single image.” In Proceedings of IEEE CVPR, 901–908. [Lev 07]A. Levin, R. Fergus, F. Durand, and W. T. Freeman. 2007. “Image and depth from a conventional camera with a coded aperture.” ACM Transaction on Graphics 26, 3, 70. [Lu 94]J. Lu, and D. M. H. Jr. 1994. “Contrast enhancement via multiscale gradient transformation.” In IEEE International Conference on Image Processing, 482–486. [Nar 00]S. G. Narasimhan, and S. K. Nayar. 2000. “Chromatic framework for vision in bad weather.” In Proceedings of IEEE CVPR, 598–605. 60.

(61) [Nay 99]S. K. Nayar, and S. G. Narasimhan. 1999. “Vision in bad weather.” In Proceedings of IEEE CVPR, 820. [Nar 02]S. G. Narasimhan and S. K. Nayar. “Vision and the atmosphere” Int. J. Comp. Vis., 48:233–254, 2002. [Oak 07]J. P. Oakley, and H. Bu. 2007. “Correction of simple contrast loss in color images.” IEEE Transactions on Image Processing 16, 2, 511–522. [Pér 98]P. Pé rez. 1998.”Markov random fields and images.” In CWI Quar- terly, vol. 11, 413–437 [Ros 99] M. V. Rossum, and T. Nieuwenhuizen. 1999. “Multiple scattering of classical waves: microscopy, mesoscopy and diffusion.”vol. 71, 313–371. [Rah 96]Z. Rahman, D. Jobson, and G. Woodell. 1996. “Multiscale retinex for color image enhancement”. [Sch 07]Y. Y. Schechner, and Y. Averbuch. 2007. “Regularized image recovery in scattering media.” IEEE Transactions on Pattern Analysis and Machine Intelligence 29, 9, 1655–1660. [Sch 07]Y. Y. Schechner, S. G. Narasimhan, and S. K. Nayar. “Polarization-based vision through haze.” App. Opt., 42:511– 525, 2003. [Sch 01] Y. Y. Schechner, S. G. Narasimhan, and S. K. Nayar. 2001. “Instant dehazing of images using polarization.” 325–332. [Shw 06]S. Shwartz, E. Namer, and Y. Y. Schechner. 2006. “Blind haze separation.” In Proceedings of IEEE CVPR, 1984–1991. [Tan 00]K. Tan, and J. P. Oakley. 2000.” Enhancement of color images in poor visibility conditions.” Proceedings of International Conference on Image Processing 2, 788–791. [Tan 08]R. T. Tan. 2008. “Visibility in bad weather from a single image.” Proceedings of IEEE CVPR.. 61.

(62) [Wik 07]WIKIPEDIA, 2007. Unsharp masking — wikipedia, the free encyclopedia. [Yua 07]L. Yuan, J. Sun, L. Quan, and H. Y. Shum. 2007. “Image deblurring with blurred/noisy image pairs.” ACMTransactions on Graphics 26, 3, 1. [Zib 01]M. Zibulevsky and B. A. Pearlmutter. “Blind source separation by sparse decomposition in a signal dictionary.” Neural Computations, 13:863–882, 2001. [Zha 02]Y. Zhang, B. Guindon, and J. Cihlar. 2002. “An image transform to characterize and compensate for spatial variations in thin cloud contamination of landsat images.” Remote Sensing of Environment 82 (October), 173–187.. 62.

(63)