Non-Local Image Dehazing 復現

時間 2019-11-07

標籤 non local image dehazing 復現简体版

原文原文鏈接

本文選自CVPR 2016, 文章連接Dana Berman, Tali Treibitz, Shai Avidan. Non-Local Image Dehazing 復現源碼見個人Githubhtml

無霧圖像和有霧圖像的RGB空間表示

一幅沒有霧霾的圖像能夠用幾百個不一樣的顏色很好的近似。將沒有霧霾的圖像每個像素值表示爲RGB空間的一個點，一幅圖像的全部像素在RGB空間中的位置會發生以下圖(b)所示的聚類現象：node

\[\mathbf{\mathit{I}}(x) = t(x)\mathbf{\mathit{J}}(x) + [1 - t(x)]\mathbf{\mathit{A}}(x)\](1)c++

每個聚類包含的像素點分佈在整副圖像區域內，他們具備相近的RGB值，可是距離攝像機的距離遠近不一樣。根據霧霾模型的公式，因爲在同一聚類裏的像素點分佈在遠近不一樣的位置因此同一聚類裏不一樣的像素點 $t$ 取不一樣值（ $t$ 是一個只和景物據攝像機距離有關的量），所以無霧霾圖像的聚類在有霧霾圖像中被拉申成直線，稱爲霧霾線。直線的一端座標值爲無霧霾圖像的聚類點的座標，另外一端爲環境光 $\mathbf{\mathit{A}}$ 。對於霧霾圖像中的全部聚類點，其對應霧霾線相交與一點，該點座標即爲環境光 $\mathbf{\mathit{A}}$ 。git

霧霾去除算法

檢測霧霾線

在本實驗中，選取Single image haze removal using dark channel prior的環境光估計算法，對於一幅霧霾圖像計算其暗通道，而後取暗通道值前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。下面給出C++實現剖分的方法：github

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繪製出來，以下：算法

<img src="http://images2015.cnblogs.com/blog/810956/201705/810956-20170518175049369-1063660519.jpg" alt="picture_2「 />
ide

爲了實現快速的查找，對還需對剖分之後的點集創建KD樹（由於咱們只在意每個像素點的$\theta,\phi$值，且剖分後的點都在球面上具備相同的 $r$ ，因此只須要對剖分後的點的$(\theta,\phi)$座標創建KD樹)。下面是建樹的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，因此同一種顏色並不必定是同一個聚類)：優化

<img src="http://images2015.cnblogs.com/blog/810956/201705/810956-20170518175212728-434800058.jpg" alt="picture_3「 />
ui

其對應的原圖以下：

傳輸係數初始值的估計

對於一條給定的霧霾線，裏面的全部像素具備近似相等的$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)
本文估計的傳輸係數以下：

傳輸係數正則化

因爲$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算法。正則化結果以下：

去霧

去霧步驟很簡單，將上面估計出來的傳輸係數和環境光強度帶入下面公式便可達到去霧的目的。
\[\tilde{J} = (I(x) - [1 - \tilde{t}(x)]A)/\tilde{t}(x)\]
注意的是，必定要防止數值的一出，若是單一通道的數值溢出很容易出現色斑。本人前面的工做已經在兩個星期以前完成，可是因爲沒有檢測道溢出的問題使這個工做一直道今天才完成
去霧結果以下，這個結果沒有作論文中所說的對比度的拉伸，因此結果有些偏暗：

這些工做是本人研一代陪用來打發時間的，有錯誤歡迎指正，謝謝！

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。