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)\]
注意的是，必定要防止数值的一出，若是单一通道的数值溢出很容易出现色斑。本人前面的工做已经在两个星期以前完成，可是因为没有检测道溢出的问题使这个工做一直道今天才完成
去雾结果以下，这个结果没有作论文中所说的对比度的拉伸，因此结果有些偏暗：

这些工做是本人研一代陪用来打发时间的，有错误欢迎指正，谢谢！