高级去雾算法与原理--暗原色图技术
原文：http://blog.csdn.net/huixingshao/article/details/42834939
很好的PPT解釋：http://wenku.baidu.com/link?url=zbOrNev-7H0Qn-RczEIUsauPEx4ybYtN1T2DCa-tkTI8S6jyQbO8CVBryqcmrMuY6kMJq2PfVT9vQ7AFXFhEAlmrNGKDuwjqzL5pE2wEJ4_
? ? 在圖像去霧這個領域，幾乎沒有人不知道《Single Image Haze Removal Using Dark Channel Prior》這篇文章，該文是2009年CVPR最佳論文。作者何凱明博士，2007年清華大學畢業，2011年香港中文大學博士畢業，可謂是功力深厚，感嘆于國內一些所謂博士的水平，何這樣的博士才可以真正叫做Doctor。
???? 關于何博士的一些資料和論文，大家可以訪問這里：http://research.microsoft.com/en-us/um/people/kahe/
???? 最開始接觸何的這篇論文是在2011年，說實在的那個時候，只是隨便瀏覽了下，看到里面的soft matting過程比較復雜，并且執行速度非常慢，就沒有什么大的興趣。最近又偶爾拾起，仔細研讀，覺得論文的推理步驟特別清晰，講解很到位。恰好適逢瀏覽到其另外一篇文章《Guided Image Filtering》?，其中提到了可以用導向濾波來代替soft matting的過程，且速度很快，因此，我對去霧的興趣算法又大大提高了。 ?
???? 本文主要上是對《Single Image Haze Removal Using Dark Channel Prior》的翻譯、整理、及部分解釋。如果您的英文水平好，建議看原文可能來的更爽些。
????一、論文思想的簡單描述?
?????首先看看暗通道先驗是什么：
? ?? 在絕大多數非天空的局部區域里，某一些像素總會有至少一個顏色通道具有很低的值。換言之，該區域光強度的最小值是個很小的數。
我們給暗通道一個數學定義，對于任意的輸入圖像J，其暗通道可以用下式表達：
????????????????????????????????????
??????式中Jc表示彩色圖像的每個通道?，Ω(x)表示以像素X為中心的一個窗口。?
?? 式（5）的意義用代碼表達也很簡單，首先求出每個像素RGB分量中的最小值，存入一副和原始圖像大小相同的灰度圖中，然后再對這幅灰度圖進行最小值濾波，濾波的半徑由窗口大小決定，一般有WindowSize = 2?*?Radius + 1;??????????
??????暗通道先驗的理論指出：
????????????????????????????????????????????????????????? ?????????????
?????實際生活中造成暗原色中低通道值主要有三個因素：a)汽車、建筑物和城市中玻璃窗戶的陰影，或者是樹葉、樹與巖石等自然景觀的投影；b)色彩鮮艷的物體或表面，在RGB的三個通道中有些通道的值很低（比如綠色的草地／樹／植物，紅色或黃色的花朵／葉子，或者藍色的水面）；c)顏色較暗的物體或者表面，例如灰暗色的樹干和石頭。總之，自然景物中到處都是陰影或者彩色，這些景物的圖像的暗原色總是很灰暗的。
????? 我們拋開論文中列舉的那些例子，自己從網上找幾幅沒有霧的風景照，看看結果如下：
???? ???
????????
????????
???????????? ???? 一些無霧的圖片??????????????????????????????????????????????????????????????????????????????? ? 其暗通道
在看看一些有霧的圖的暗通道：
??????
??????
一些有霧的圖片?????????????????????????????????????????????????????????????????????????????? ???其暗通道
上述暗通道圖像均使用的窗口大小為15*15，即最小值濾波的半徑為7像素。
????? 由上述幾幅圖像，可以明顯的看到暗通道先驗理論的普遍性。在作者的論文中，統計了5000多副圖像的特征，也都基本符合這個先驗，因此，我們可以認為其實一條定理。
????? 有了這個先驗，接著就需要進行一些數學方面的推導來最終解決問題。
首先，在計算機視覺和計算機圖形中，下述方程所描述的霧圖形成模型被廣泛使用：
???????????????????????????????????????????????? ??
自己注釋：因為光線在傳播過程中有三種分散，分別是透射、散射、吸收。當時這里假設的理論是只有透射、散射，不考慮吸收。并且透射率是隨著景物深度的增加而衰減的。
? 其中，I(X)就是我們現在已經有的圖像（待去霧的圖像），J(x)是我們要恢復的無霧的圖像，A是全球大氣光成分， t(x)為透射率。現在的已知條件就是I(X)，要求目標值J(x),顯然，這是個有無數解的方程，因此，就需要一些先驗了。
將式（1）稍作處理，變形為下式：
????????????????????????????????????????????????? ??
??? 如上所述，上標C表示R/G/B三個通道的意思。
??? 首先假設在每一個窗口內透射率t(x)為常數，定義他為，并且A值已經給定，然后對式（7）兩邊求兩次最小值運算，得到下式：
??????????????????????????????? ??
??? 上式中，J是待求的無霧的圖像，根據前述的暗原色先驗理論有：
???????????????????????????????????????????????
???? 因此，可推導出：
????????????????????????????????????????????????????? ???
??? 把式（10）帶入式（8）中，得到：
?????????????????????????????????????????????????
????這就是透射率的預估值。
??? 在現實生活中，即使是晴天白云，空氣中也存在著一些顆粒，因此，看遠處的物體還是能感覺到霧的影響，另外，霧的存在讓人類感到景深的存在，因此，有必要在去霧的時候保留一定程度的霧，這可以通過在式（11）中引入一個在[0,1] 之間的因子，則式（11）修正為：
?????????????????????????????????????????? ????
???? 本文中所有的測試結果依賴于：? ω=0.95。
???? 上述推論中都是假設全球達氣光A值時已知的，在實際中，我們可以借助于暗通道圖來從有霧圖像中獲取該值。具體步驟如下：
??? 1） 從暗通道圖中按照亮度的大小取前0.1%的像素。
????????? 2） 在這些位置中，在原始有霧圖像I中尋找對應的具有最高亮度的點的值，作為A值（天空處I跟A越相似，則其值越接近1，t就越接近0，所以其會導致偏色）。
???? 到這一步，我們就可以進行無霧圖像的恢復了。由式（1）可知：??J =?( I - A)/t + A??
?????現在I,A,t都已經求得了，因此，完全可以進行J的計算。
???? 當投射圖t 的值很小時，會導致J的值偏大，從而使淂圖像整體向白場過度，因此一般可設置一閾值T0，當t值小于T0時，令t=T0，本文中所有效果圖均以T0=0.1為標準計算。
???? 因此，最終的恢復公式如下：
?????????????????????????????
自己的見解：上面的A的求取是基于有霧的圖像的得到的，作為分母其是不可能為零的。根據(12)得出A值越大這透射率的值越大，則表明圖像霧越少。所以這個算法對對于天空的處理有兩種，第一要么把具有天空的單獨出來不處理。
白色的區域的則回進行
???? 當直接用上述理論進行恢復時，去霧的效果其實也是很明顯的，比如下面一些例子：
???????
????? ?
?????? 有霧圖 ??????????????? 去霧圖
??????注意到第一幅圖的原圖兩個字的周圍明顯有一塊不協調的地方，而第二圖頂部水平方向似乎有一塊沒有進行去霧處理，這些都是由于我們的透射率圖過于粗糙了。
????? 要獲得更為精細的透射率圖，何博士在文章中提出了了soft matting方法，能得到非常細膩的結果。但是他的一個致命的弱點就是速度特慢，不使用于實際使用。在2011年，何博士又除了一片論文，其中提到了導向濾波的方式來獲得較好的透射率圖。該方法的主要過程集中于簡單的方框模糊，而方框模糊有多重和半徑無關的快速算法。因此，算法的實用性特強，關于這個導向濾波算法大家在何博士的網站可以自己去研習下，除了在去霧方面外，還有著其他多方面的應用，這部分本文不多述。
???? 使用了導向濾波后的去霧效果：
?????
??????
??? 使用原始的預估透射率圖 ?????????? ?使用導向濾波后的透射率圖
??????
?? (a) 原圖 ??? (b) 去霧結果圖
??????
?????? ?（c）??? 暗通道圖 （d） 導向圖 （原始圖像的灰度圖）
???? ?
????????? （e)? ?預估透射率圖 (f)???使用導向濾波后的透射率圖
二、各參數對去霧結果的影響
?第一：窗口的大小。這個對結果來說是個關鍵的參數，窗口越大，其包含暗通道的概率越大，暗通道也就越黑。我們不去從理論角度分析，從實踐的效果來看，似乎窗口越大，去霧的效果越不明顯，如下圖所示：
?????? ?
?????????????????????????? ?(a) 原始圖像? (b) 窗口大小=11
??????
?(c) 窗口大小=21 (d) 窗口大小=101?
我的建議是窗口大小在11-51之間，即半徑在5-25之間。
???? 式（12）中的ω具有著明顯的意義，其值越小，去霧效果越不明顯，舉例如下：
??????
??????????????????????????????????? (a) 原始圖像? ??? ?? (b)??? ω=0.5?????????
??????????
? ?????????????????????????????????????(c)??? ω=0.8?????????????????????????????????????????????????????????????????? ?? (d)??? ω=1?
????三：編碼的步驟
如果你仔細的分析了原文的細路，加上適當的參考，編碼其實并不是很困難。
1）根據原始圖像求暗通道，參考代碼如下：
for (Y = 0, DarkPt = DarkChannel; Y < Height; Y++){ImgPt = Scan0 + Y * Stride;for (X = 0; X < Width; X++){Min = *ImgPt;if (Min > *(ImgPt + 1)) Min = *(ImgPt + 1);if (Min > *(ImgPt + 2)) Min = *(ImgPt + 2);*DarkPt = Min;ImgPt += 3;DarkPt++;}}MinFilter(DarkChannel, Width, Height, Radius);????這里需要注意的是MinFilter算法的快速實現，提供一篇論文供有需要的朋友學習：STREAMING MAXIMUM-MINIMUM FILTER USING NO MORE THAN THREE COMPARISONS PER ELEMENT?。這個算法的時間復雜度是O(1)的。
????? 2）按文中所描述的算法自動獲得全球大氣光的值。
?????這里說明一點，原始論文中的A最終是取原始像素中的某一個點的像素，我實際上是取的符合條件的所有點的平均值作為A的值，我這樣做是因為，如果是取一個點，則各通道的A值很有可能全部很接近255，這樣的話會造成處理后的圖像偏色和出現大量色斑。原文作者說這個算法對天空部分不需特備處理，我實際發現該算法對有天空的圖像的效果一般都不好。天空會出現明顯的過渡區域。作為解決方案，我增加了一個參數，最大全球大氣光值，當計算的值大于該值時，則就取該值。?
??? ???? ?
??????????????????????? 原圖???????????????????????????????????????????????????????? 未對A值做限定 最大A值限定為220
?????? 3） 按式（12）計算預估的透射率圖。
在式（12）中，每個通道的數據都需要除以對應的A值，即歸一化，這樣做，還存在一個問題，由于A的選取過程，并不能保證每個像素分量值除以A值后都小于1，從而導致t的值可能小于0，而這是不容許的，原文作者并沒有交代這一點是如何處理的。我在實際的編碼中發現，如果真的這樣做了，其效果也并不是很理想?，因此，我最后的辦法是在式（12）中，不考慮A的計算。
??????? 4）計算導向濾波圖。
這里可以直接用原始的圖像做導向圖，當然也可以用其灰度圖，但是用RGB導向圖在下一步的計算中會占用比較大的時間。
??????? 5）按照《Guided Image Filtering》論文中的公式（5）、（6）、（8）編碼計算獲得精細的透射率圖。
網絡上有這個算法的 matlab代碼可下載的，這里貼部分代碼：
function q = guidedfilter(I, p, r, eps)% GUIDEDFILTER O(1) time implementation of guided filter.%% - guidance image: I (should be a gray-scale/single channel image)% - filtering input image: p (should be a gray-scale/single channel image)% - local window radius: r% - regularization parameter: eps[hei, wid] = size(I);N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.% imwrite(uint8(N), 'N.jpg');% figure,imshow(N,[]),title('N');mean_I = boxfilter(I, r) ./ N;mean_p = boxfilter(p, r) ./ N;mean_Ip = boxfilter(I.*p, r) ./ N;cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.mean_II = boxfilter(I.*I, r) ./ N;var_I = mean_II - mean_I .* mean_I;a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;b = mean_p - a .* mean_I; % Eqn. (6) in the paper;mean_a = boxfilter(a, r) ./ N;mean_b = boxfilter(b, r) ./ N;q = mean_a .* I + mean_b; % Eqn. (8) in the paper;end?????? 由上面的代碼，可見，主要的工作量在均值模糊上，而均值模糊是個很快速的算法，關于均值模糊的優化可參考我以前的文章：彩色圖像高速模糊之懶惰算法。
??????還有一點就是，上述計算需要在[0,1]范圍內進行，也就是說導向圖和預估的透射率圖都必須從[0,255]先映射到[0,1]在進行計算。
????? 關于guidedfilter中的半徑r值，因為在前面進行最小值后暗通道的圖像成一塊一塊的，為了使透射率圖更加精細，建議這個r的取值不小于進行最小值濾波的半徑的4倍，如下圖所示：
?????? ??
（a)? r=最小值濾波半徑的2倍 (b) r=最小值濾波半徑的8倍
??????可以看到，當r比較小的時候，在透射率圖中基本看不到什么細節信息，因此恢復處的圖像邊緣處不明顯。
????? 參數eps的取值也有所講究，他主要是為了防止計算中除以0的錯誤以及為了使得某些計算結果不至于過大，一般建議取值0.001或者更小。
????? 如果使用的彩色RGB圖做導向圖，計算時間上會增加不少，所的到的透射率圖的邊緣會比灰度圖所處理的保留了更多的細節，效果上略微比灰度圖好，如下所示：
??????? ?
(a) 原圖??????????????????????????????????????????????????????????????????????????????? ??? (b)預估的透射率圖
?? ??????
??????????????????? ?? (c)使用灰度圖為導向圖得到的透射率圖 ?????? ??? (d)使用RGB圖為導向圖得到的透射率圖
?????????
???? (e)灰度圖為導向圖對應的去霧效果???? ??????? （f）RGB圖導向圖對應的去霧效果
???????以RGB圖為導向圖的計算中，涉及到3*3部分矩陣求逆的過程，如果用非matlab語言寫，可以先借助于matlab的符號計算功能，以及其中的符號計算命令simple，把計算結果算出來，然后再再其他高級語言中實現。
?? ?? ?(6) 按式（22）進行無霧圖像的恢復。
?四、其他一些去霧效果圖
??????
????? ????
? ????????
????? ????
????????????????????????????? 原圖???????????????????????????????????????????????? 去霧效果圖????????????????????????????????????????????????????? 透射率圖
??????上圖中最后一副圖我連續做了兩次去霧的處理。?
????? 在原文中，有這樣一段話：
?Since the scene radiance is usually not as bright as the atmospheric light, the image after haze removal looks dim. So we increase the exposure of J(x) for
display.
?????意思就是說直接去霧后的圖像會比原始的暗，因此在處理完后需要進行一定的曝光增強，但作者沒有說明其是如何增強的， 因此，這里的圖和他論文的效果有所不同時正常的。一般在去霧處理后再用自動色劑之類的算法增強下會獲得比較滿意的結果，如下圖：
?????? ????? ?
??? 原圖 ? 去霧后 +自動色階
????? 去霧算法目前也有著眾多其他的方式，不過我所接觸的，很多都是以這個為基礎，因此，先弄會這個為研究其他的去霧算法能奠定堅實的基礎。
??????網絡上有一些比較好的暗原色先驗去霧的matlab代碼：比如和本文基本對應的matlab資源：http://files.cnblogs.com/Imageshop/cvpr09defog%28matlab%29.rar
后記：稍微有看了幾篇去霧的文章，基本上都是圍繞著獲得透視率圖做文章，比如一些文章講用聯合雙邊濾波方式獲得精細的透射率，從我個人淺薄的認識中，我覺得去霧已基本上跳不出暗原色這個范疇了。
????? 我對雙邊濾波算法那也做了實驗，發現這個的效果也還行，就是速度慢了很多，雙邊濾波的快速算法其實快不起來的，所以這個的實用性不高，我選擇了一些圖像做比較：
?????????
? ???????????????????????? ?? (a) 原圖?????????????????????????????????????????????????????????????? （b)?? 聯合雙邊濾波去霧圖
?????????
? ???????????????????? （c） 導向濾波獲得透射率圖?????????????????????????????????????????????????? (d)聯合雙邊濾波透射率圖(Sigmad=SigmaR=100)
???????上圖可以很明顯的看出聯合雙邊濾波的透射率圖沒有導向濾波的精細，但比原始的粗糙的透射率圖還是好很多，過渡很光滑，因此，也能獲得不錯的視覺去霧效果。
?????? 聯合雙邊濾波器中的算法是參考了OpenCV中相關函數寫的。
?????和平常一樣，提供一個可供大家測試效果的程序：?基于暗原色先驗的圖像去霧演示程序
????
??? 我分別用VB6和C#做了個程序，兩個程序都已經通過各自的語言方式進行了優化，算法部分編碼是一樣的，C#運行速度大約是VB6的1.8倍。
???? 在處理速度上，比 matalb的快了很多倍，在I3的筆記本電腦上，一副1024*768的圖像去霧時間大約在150ms內（以灰度圖為導向圖）。
?
?五、算法的局限性
?
????? 暗原色先驗是一種統計的結果，是對大量戶外無霧照片(outdoor haze-free images)的統計結果，如果目標場景內在的就和大氣光類似，比如雪地、白色背景墻、大海等，則由于前提條件就不正確，因此一般無法獲得滿意的效果，而對于一般的風景照片這個算法能處理的不錯。
????
?
????2013.8.23 后記補充修正：
?????在后續對該算法的關注中，發現自己在前面做出了一個錯誤的判斷，就是關于式（11）中/A的操作。我在前面說這個除法會引起一些問題，因此，去除了這一步。但是后來的實踐證明正是有了這一步，對于對比度低的圖像才可以獲得很好的去霧高對比度圖。
????前面說的/A操作可能會導致t的值小于0，這種情況就可以把t的值直接設置為0來解決。
??? 還有一個事情就是式（11）嚴格的來說是要對原始圖像的每個通道進行歸一化后，再取每個通道R/G/B值的最小值得到中間圖，然后對這個中間圖進行指定半徑的最小值濾波后，通過11式得到粗糙的透射率圖，那么這樣就需要多不少計算，我在實際中發現如果直接用前面的暗通道圖/A進行操作，兩者的效果區別不明顯，因此，可用這種簡便的方式。
???
?上圖是一副經典的測試圖，雖然取得了比較好的效果，不過似乎馬路那一塊的效果不如一些其他人公開的成果那么好。
???
?????這也是一副比較常見的測試圖，該圖也是緊用去霧獲得結果，未做任何的后處理，同CSDN一個的案例庫：圖像去霧的算法研究中的效果相比，在整幅圖像的對比和和協調性上都要好一個檔次。
???? 再如下圖，也比CSDN那個案例庫中的效果要好很多。
?????
???? 還有：
??
??
???
???總結：我對這種去霧算法的效果還是很滿意的，?效果和速度都還比較合適。
???附件的測試程序已經同步更新。
???2013.10,10 后記補充:
????一直聽說C的效率很高，于是用C實現了該算法，并且編制了一個可供其他語言調用的dll文件，然后對最小值濾波算法又進行了特殊的優化，算法速度有了很大的提高，特別是在用VS2010時，編寫C的代碼可以勾選下圖中的SSE選項和快速（fp:fast），程序會調用SSE一些函數，實現指令級別的并行。而這些C#的編譯器是無法實現的。
?????同樣的算法，相對于原始的C#程序，程序的速度能提高一倍左右，對于800*600的圖像，在我的I3的CPU上平均能能達到20fps的速度（只占用了一個核的CPU資源），因此可以適用于不大于該范圍內的實時圖像處理。
???? 同樣我增加了3個內部可調的參數供大家測試。
???? 相應的DLL函數聲明如下：
???? c風格：
void HazeRemovalUseDarkChannelPrior(unsigned char * Src,unsigned char * Dest,int Width,int Height,int Stride, int Radius ,int GuideRadius, int MaxAtom, float Omega,float Epsilon,float T0 )??? c#調用：
[DllImport("HazeRemoval.dll", CallingConvention = CallingConvention.StdCall, CharSet = CharSet.Unicode, ExactSpelling = true)]private static extern void HazeRemovalUseDarkChannelPrior(byte* Src, byte* Dest, int Width, int Height, int Stride, int Radius,int GuideRadius, int MaxAtom, float Omega, float Epsilon, float T0);????VB6調用：
Private Declare Sub HazeRemovalUseDarkChannelPrior Lib "HazeRemoval.dll" (ByVal Src As Long, ByVal dest As Long, ByVal Width As Long, ByVal Height As Long, ByVal Stride As Long, ByVal Radius As Long, ByVal GuideRadius As Long, ByVal MaxAtom As Long, ByVal Omega As Single, ByVal Epsilon As Single, ByVal T0 As Single)?????
???? 調用實例源代碼下載：http://files.cnblogs.com/Imageshop/HazeRemovalTest.rar
?2013.11,22修正：?
通過聯合雙邊濾波求透射率圖的效果上面的不正確的，進最新的研究結果表明，雙邊濾波確實也能獲得精細的投射圖，比如依舊那上面的測試圖像，不同的SigmaS和SigmaR下獲得透射率效果如下：
?????? ?
SigmaS=50,SigmaR=30 ??? SigmaS=100,SigmaR=30
?
?2014.1.12 重大更新（可實現實時效果）：
? ? ?何的算法效果以及普遍的實用性都比其他的去霧算法要好，而主要的問題就是其速度還是不夠快，有著太多的浮點計算。鑒于此，作者也多次試著對代碼進行深層次的優化，包括SSE處理、并行運行等，但由于算法本身的順序執行，無法全程并行，偶爾一個小函數可以并行，但由于其本身執行就特別快，比如不要5ms，你去用并行算法可能耗時還會大一些。因此，一直沒有什么大的進步，對于一副1024*768的彩圖進行去霧需要90ms，這肯定無法滿足需求。
? ? ?最近，在思考，既然暗通道去霧的透射率圖比其他的算法都來的精細，如果適當的降低一點點其精度，其去霧的效果理論上應該不會有太大的區別，于是我想到了一種方式，即求取透射率的時候不是對原圖進行求取，而是先對原圖進行下采樣，比如縮小為原圖的1/4，計算出小圖的透射率，之后在通過插值的方式的獲取原圖大概的透射率，則應該也可以獲得效果。經過實踐，這種方式大大的提高了執行速度，而且效果和原始的方案基本一致，對于1024*768的圖像大約只需要（I3CPU）30ms了，如果進一步取1/9的縮放，則只需要大約20ms，完全可以滿足工業實時性要求高的場合。
? ? 當然，如果你的縮小系數不是特別大的話，比如縮小為原來的0.5大小，可能兩次縮放所用的耗時還抵消了計算小圖的透射率圖所換來的盈利，因此必須合理選擇這個下采樣率。
? ? 要實現這樣的速度，當然還是需要很高的優化技巧的，這些東西還是有所保留比較好。?
? ? 我做了一個程序，集成了本博客中6種圖像去霧的算法：?圖像去霧綜合版本
用暗通道的算法對一段視頻進行了處理，大家可以在這里看到效果：http://video.sina.com.cn/v/b/124538950-1254492273.html
總結
以上是生活随笔為你收集整理的高级去雾算法与原理--暗原色图技术的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 工作篇 之 解决谷歌地图与 Nested
- 下一篇： iOS模拟器使用双指手势