本文選自CVPR 2016, 文章鏈接Dana Berman, Tali Treibitz, Shai Avidan. Non-Local Image Dehazing
復現源碼見我的Github
無霧圖像和有霧圖像的RGB空間表示
一幅沒有霧霾的圖像可以用幾百個不同的顏色很好的近似。將沒有霧霾的圖像每一個像素值表示為RGB空間的一個點，一幅圖像的所有像素在RGB空間中的位置會發生如下圖(b)所示的聚類現象：
\[\mathbf{\mathit{I}}(x) = t(x)\mathbf{\mathit{J}}(x) + [1 - t(x)]\mathbf{\mathit{A}}(x)$$(1) 每一個聚類包含的像素點分布在整副圖像區域內，他們具有相近的RGB值，但是距離攝像機的距離遠近不同。根據霧霾模型的公式，由於在同一聚類里的像素點分布在遠近不同的位置所以同一聚類里不同的像素點 $t$ 取不同值（ $t$ 是一個只和景物據攝像機距離有關的量），因此無霧霾圖像的聚類在有霧霾圖像中被拉申成直線，稱為霧霾線。直線的一端坐標值為無霧霾圖像的聚類點的坐標，另一端為環境光 $\mathbf{\mathit{A}}$ 。對於霧霾圖像中的所有聚類點，其對應霧霾線相交與一點，該點坐標即為環境光 $\mathbf{\mathit{A}}$ 。 <div align="center"> <img src="http://images2015.cnblogs.com/blog/810956/201705/810956-20170518222948510-1621450583.png" alt="picture_1" /> </div> ## **霧霾去除算法** ### **檢測霧霾線** 在本實驗中，選取[Single image haze removal using dark channel prior](http://cn.bing.com/academic/profile?id=00c685d62fb1b326466854c7302cfa20&encoded=0&v=paper_preview&mkt=zh-cn)的環境光估計算法，對於一幅霧霾圖像計算其暗通道，然后取暗通道值前0.1%的像素點中最亮的像素值作為環境光的估計。\ 然后定義$\mathbf{\mathit{I}}_{A}$如下： $$\mathbf{\mathit{I}}_{A}(x) = \mathbf{\mathit{I}}(x) - \mathbf{\mathit{A}} = t(x)[\mathbf{\mathit{J}}(x) - \mathbf{\mathit{A}}]$$ (2) 上面的公式將圖像在RGB空間中的像素坐標進行了平移變換，環境光的坐標被變換到了原點。將$\mathbf{\mathit{I}}_{A}$表示成球坐標的形式： $$\mathbf{\mathit{I}}_{A}(x) = [r(x),\theta(x),\phi(x)]$$(3) 這樣霧霾圖像的像素坐標就被表示成為球坐標空間中圍繞着環境光（坐標原點）分布的坐標點。我們來觀察(2)式，對於具有相同$\mathbf{\mathit{J}}$和$\mathbf{\mathit{A}}$的景物點，其離攝像機的距離只影響 $t$ 的取值，並且在球坐標系中改變 $t$ 只影響 $r$ 。因此對於無霧霾圖像中的每一個聚類點，其對應的霧霾線在上述變換以后的球坐標系中都具有相同的$\theta,\phi$值。也就是說，具有相同的$\theta,\phi$值的像素點其對應的無霧霾圖像的像素具有近似的值。為了確定哪些像素具有相同的$\theta,\phi$值，需要將圖像根據$\theta,\phi$進行聚類。為此要先要對球面進行等距離剖分，需要注意的是等分$[0,2\pi]\times[0,\pi]$是無法得到均勻剖分的點的，你可以看一下地球儀的經緯剖分結果是否是這樣。理論上說如果想要在一個球面上得到20個以上等距離剖分的點是不可能的，原因按住不表，有興趣的可以搜索伯努利多面體。但是可以通過迭代細分正二十面體的方法來近似得到等距離剖分，具體步驟見[Colorado State University General Circulation Model](http://kiwi.atmos.colostate.edu/rr/groupPIX/ross/ross1/ross1.html)。下面給出C++實現剖分的方法： ```C++ void subdivide(icosahedron& src, polyhedron& dst, int num) { dst = src.i; double r = src.radius; while(dst.vertex_table.size() < num) { std::map<std::vector<int>, int> mid_table; const int size = dst.plane_table.size(); for(int i = 0; i != size; ++i) { auto f = *(dst.plane_table.begin()); dst.plane_table.pop_front(); int mid_idx12, mid_idx13, mid_idx23; if(!is_present(f[0], f[1], mid_table, mid_idx12)) { cv::Point3d pt1 = dst.vertex_table[f[0]], pt2 = dst.vertex_table[f[1]]; dst.vertex_table.emplace_back((pt1.x + pt2.x)/2, (pt1.y + pt2.y)/2, (pt1.z + pt2.z)/2); mid_idx12 = static_cast<int>(dst.vertex_table.size()) - 1; scale2unit(dst.vertex_table[mid_idx12], r); std::vector<int> temp = { f[0], f[1] }; mid_table[temp] = mid_idx12; } if(!is_present(f[0], f[2], mid_table, mid_idx13)) { cv::Point3d pt1 = dst.vertex_table[f[0]], pt2 = dst.vertex_table[f[2]]; dst.vertex_table.emplace_back((pt1.x + pt2.x)/2, (pt1.y + pt2.y)/2, (pt1.z + pt2.z)/2); mid_idx13 = static_cast<int>(dst.vertex_table.size()) - 1; scale2unit(dst.vertex_table[mid_idx13], r); std::vector<int> temp = { f[0], f[2] }; mid_table[temp] = mid_idx13; } if(!is_present(f[1], f[2], mid_table, mid_idx23)) { cv::Point3d pt1 = dst.vertex_table[f[1]], pt2 = dst.vertex_table[f[2]]; dst.vertex_table.emplace_back((pt1.x + pt2.x)/2, (pt1.y + pt2.y)/2, (pt1.z + pt2.z)/2); mid_idx23 = static_cast<int>(dst.vertex_table.size()) - 1; scale2unit(dst.vertex_table[mid_idx23], r); std::vector<int> temp = { f[1], f[2] }; mid_table[temp] = mid_idx23; } std::vector<int> t = { f[0], mid_idx12, mid_idx13 }; dst.plane_table.push_back(t); t[0] = mid_idx23; dst.plane_table.push_back(t); t[2] = f[1]; dst.plane_table.push_back(t); t[1] = mid_idx13; t[2] = f[2]; dst.plane_table.push_back(t); } } } ``` 為了驗證剖分的正確性，將剖分結果用MATLAB繪制出來，如下： <div align="center"> <img src="http://images2015.cnblogs.com/blog/810956/201705/810956-20170518175049369-1063660519.jpg" alt="picture_2“ /> </div> 為了實現快速的查找，對還需對剖分以后的點集建立KD樹（因為我們只在乎每一個像素點的$\theta,\phi$值，且剖分后的點都在球面上具有相同的 $r$ ，所以只需要對剖分后的點的$(\theta,\phi)$坐標建立KD樹)。下面是建樹的C++實現: ```C++ kd_node* build_kdTree(std::vector<cv::Point2d>& sph_table, kd_node* p, std::vector<int>& subset) { kd_node* r = new kd_node; r->parent = p; if(subset.size() == 1) { r->data = subset[0]; r->dimension = 0; r->left = r->right = nullptr; r->is_leaf = true; return r; } std::vector<std::vector<int>> subsets; r->dimension = dimension_choice(sph_table, subset); r->data = split(sph_table, subset, subsets, r->dimension); r->is_leaf = false; r->left = subsets[0].size() != 0 ? build_kdTree(sph_table, r, subsets[0]) : nullptr; r->right = subsets[1].size() != 0 ? build_kdTree(sph_table, r, subsets[1]) : nullptr; return r; } ``` 為了驗證聚類結果， 將聚類結果用不同的顏色顯示如下(因為一共使用了三種顏色，聚類的數目遠大於3，所以同一種顏色並不一定是同一個聚類)： <div align="center"> <img src="http://images2015.cnblogs.com/blog/810956/201705/810956-20170518175212728-434800058.jpg" alt="picture_3“ /> </div> 其對應的原圖如下： <div align="center"> <img src="http://images2015.cnblogs.com/blog/810956/201705/810956-20170518175229885-1867568387.jpg" alt="picture_4“ /> </div> ### **傳輸系數初始值的估計** 對於一條給定的霧霾線，里面的所有像素具有近似相等的$J$和$A$的值, $r(x)$有下式決定： $$r(x) = t(x)\left \| J - A \right \|$$(4) 當$t = 1$對應最大的徑向坐標： $$r_{max} = \left \| J - A \right \|$$(5) 聯合(4)(5)得： $$t(x) = r(x)/r_{max}$$(6) 對於每一條聚類線，由下式估計 $ r_{max} $: $$\hat{r}_{max} = \max_{x\in H}[r(x)]$$(7) 其中$H$表示每一條霧霾線。\ 所以每一個像素點可以估計出對應的傳輸系數為： $$\tilde{t}(x) = r/\hat{r}_{max}$$(8) 本文估計的傳輸系數如下： <div align="center"> <img src="http://images2015.cnblogs.com/blog/810956/201705/810956-20170518175244947-1624899983.jpg" alt="picture_5“ /> </div> ### **傳輸系數正則化** 由於$J$始終取正值，所以根據(1)式可以給出傳輸系數的下界: $$t_{LB}(x) = 1 - \min_{c \in {R,G,B}}{I_{c}(x)/A_{c}}$$ (9) 所以對傳輸系數進行下界約定后有： $$\tilde{t}_{LB}(x) = \max[\tilde{t}(x),t_{LB}(x)]$$(10) 除此之外，上面對傳輸系數的估計只是基於霧霾線的假設。可以預見的是一個聚類里面的像素點數目越少，上面的假設的可信度越小。還有，傳輸系數的估計還要考慮道原始像素的空間相似性，我們有理由相信領域內像素值越接近的兩個像素點其對應的傳輸系數也越接近。所以為了權衡上面的兩個條件，最終的傳輸系數有下面的最小化問題給出： $$\sum_{x}\frac{[\tilde{t}(x) - \tilde{t}_{LB}(x)]^2}{\delta^2(x)} + \lambda\sum_{x}\sum_{y\in N_{x}}\frac{[\tilde{t}(x) - \tilde{t}(y)]^2}{\left \| I(x) - I(y) \right \|^2}$$(11) 其中$\lambda$的作用是用來權衡兩個假設之間的關系，這里面取0.1。$N_{x}$表示 $x$ 像素的四領域。$\delta(x)$是每一條霧霾線得到的傳輸系數估計值的標准差。原文中指出$\delta(x)$扮演着一個重要角色，因為如果一個聚類里面的點數分布隨着$\delta(x)$變大而減小，這里說的是像素點的分布情況而不是點數。 這個最優化問題的解法很簡單，因為是一個很簡單的凸函數，很容易通過求導為零來解，實際上就是求解下面一個線性方程組： $$(\frac{2}{\delta^2_{i}} + 4\lambda \sum_{p}\frac{1}{\left \| I_{i} - I_{p} \right \|^2})\tilde{t}_{i} - 4\lambda\sum_{p}\frac{\tilde{t}_{p}}{\left \| I_{i} - I_{p} \right \|^2} = \frac{2}{\delta^2_{i}}\tilde{t}_{LBi}$$(12) 本文利用的數值解法是GSL庫提供的GMRES算法。正則化結果如下： <div align="center"> <img src="http://images2015.cnblogs.com/blog/810956/201705/810956-20170518175304682-2012579814.jpg" alt="picture_6“ /> </div> ### **去霧** 去霧步驟很簡單，將上面估計出來的傳輸系數和環境光強度帶入下面公式即可達到去霧的目的。 $$\tilde{J} = (I(x) - [1 - \tilde{t}(x)]A)/\tilde{t}(x)\]
注意的是，一定要防止數值的一出，如果單一通道的數值溢出很容易出現色斑。本人前面的工作已經在兩個星期之前完成，但是由於沒有檢測道溢出的問題使這個工作一直道今天才完成
去霧結果如下，這個結果沒有做論文中所說的對比度的拉伸，所以結果有些偏暗：
