去雾算法C++实现

去雾算法代码实现

前言：暑假闲着没啥事儿就乘着这个机会好好学习学习算法原理吧。某虽不才，虽然实现效果没有大佬们的好，但是我的代码通俗易懂，并且完全开源（C++），欢迎大家前来相互学习，批评指正。

原理

本文对原理不再赘述，了解去雾算法原理请去往：原理

根据何凯明博士在论文中给出的公式，我们可以推导出最终待求图像J(x)的计算公式：
J(x) = [I(x) - A] / t(x) + A
其中I(x)就是原图像，A是大气光成分，在计算时需要BGR三个通道分开进行；t(x)是透射率，它的计算也是最麻烦的，同样需要三通道分开进行。

显然，我们需要求得A，t(x)两个参数。

步骤

1.计算暗通道图像

暗通道的计算比较简单，首先把BGR三个颜色通道分离；随后，就原图中每个像素点而言，比较出该像素点三个颜色通道灰度值中的最小值，然后将该最小值作为像素值存入暗通道图像的对应位置中。最后别忘了滤波。
代码实现入下：

darkImg = Mat::zeros(img.rows, img.cols, CV_8UC1);//创建空矩阵
    for (int i = 0; i<img.rows; i++)
    {

        for (int j = 0; j<img.cols; j++)
        {
            darkImg.at<uint8_t>(i, j) = min
                    (
                        min(img.at<Vec<uint8_t, 3>>(i, j)[0], img.at<Vec<uint8_t, 3>>(i, j)[1]),
                        min(img.at<Vec<uint8_t, 3>>(i, j)[0], img.at<Vec<uint8_t, 3>>(i, j)[2])
                    );
            if(darkImg.at<uint8_t>(i, j) == 255)
                x++;
        }
    }

    blur(darkImg, darkImgOut, Point(3,3));//均值滤波

暗通道效果图：

2.计算A值

A值的计算需要分通道进行，也就是说我们需要计算三个A值，我们设他们分别为：AB，AG，AR，以AB为例：
根据算法原理，我们首先需要在暗通道图像中寻找亮度值（灰度值）排名在前0.1%的像素点，然后再根据这些像素点的位置到原图中寻找亮度最大的像素点，并将其作为该通道的A值。

由于本人对OpenCV库函数了解还不透彻，所以不知道有没有函数可以较为方便的实现像素点筛选，于是只能自己实现（知道的同学可以评论区留言哦）。

我使用了简单的循环实现了这一功能：由于前0.1%的像素值一定很接近255，那么我们可以从255开始对暗通道图进行遍历，每降低一个灰度值，就统计一下像素点个数，如果超过了总像素点个数的0.1%就停止循环，同时把当前的灰度值作为阈值。
那么问题来了，如何用暗通道图中的这些前0.1%的像素点到原图的各个颜色通道中去寻找灰度最大的点呢？还记得那个刚刚求出的阈值吗？我们再来一次遍历，不过这次是在颜色通道中进行遍历，用if语句进行判断，条件就是暗通道中当前点的像素值是否大于阈值，如果是，那么我们就把颜色通道中的该位置的像素值装入事先定义好的数组中，如果不是就跳过进行下一个。由于我们只需比较大小，不用知道最大像素值点的位置，所以比较简单。遍历完之后我们得到一个一维数组，最后调用 求数组的最大值函数就可得出该颜色通道的A值了。

代码实现如下：

    int cou = 0, pix = 255;
    
    for(pix = 255; pix>=0; pix--)//确定灰度值（pix）为多少以上的点入选
    {
        for (int i = 0; i<img.rows; i++)
        {
            for (int j = 0; j<img.cols; j++)
            {
                if((darkImgOut.at<uint8_t>(i, j)>=pix)&&(cou<=count))
                {cou++;}
                else if(cou>count)
                {flag0 = 1;}
            }
        }
        if(flag0)//count满退出循环
            break;
    }

    Mat brightB, brightG, brightR;//各通道亮度图像
    int map[3][10000] = {{0}};//开辟一段内存

    ChannelSplit(img)[0].copyTo(brightB);
    ChannelSplit(img)[1].copyTo(brightG);
    ChannelSplit(img)[2].copyTo(brightR);

    int AB, AG, AR;//各通道A值
    int n = 0;

    if(x<count)
    {
        for (int i = 0; i<img.rows; i++)//获得索引图
        {
            for (int j = 0; j<img.cols; j++)
            {
                if(darkImgOut.at<uint8_t>(i, j) >= pix)//满足的原图像素值装入map
                {
                    map[0][n] = brightB.at<uint8_t>(i, j);
                    map[1][n] = brightG.at<uint8_t>(i, j);
                    map[2][n] = brightR.at<uint8_t>(i, j);
                    n++;
                }
            }
        }
        AB = (*max_element(map[0], map[0] + count));//求得AB
        AG = (*max_element(map[1], map[1] + count));//求得AG
        AR = (*max_element(map[2], map[2] + count));//求得AR
    }
    else if(x>=count)
    {
        for (int i = 0; i<img.rows; i++)//获得索引图
        {
            for (int j = 0; j<img.cols; j++)
            {
                if(darkImgOut.at<uint8_t>(i, j) == 255)
                {
                    map[0][n] = brightB.at<uint8_t>(i, j);//把亮度图中对应暗图前0.1%的点的像素值放入map数组中
                    map[1][n] = brightG.at<uint8_t>(i, j);
                    map[2][n] = brightR.at<uint8_t>(i, j);
                    n++;
                }
            }
        }
        AB = (*max_element(map[0], map[0] + count));//对map数组排序求得AB
        AG = (*max_element(map[1], map[1] + count));
        AR = (*max_element(map[2], map[2] + count));
    }

计算投射函数t(x)

首先需要了解窗口的概念：原理

根据何凯明博士的论文我们可得t(x)的计算公式：
t(x) = 1 - w * min { min { I(y) / A } }
其中最里层的min函数定义域是BGR三个颜色通道，因此我们同样需要对三个通道分开计算。最外层的min函数定义域是窗口中的像素点。窗口最好不要太大，为了便于计算，我把窗口大小设为9。w是一个修正参数，我将其赋值为0.95.

在代码实现中，首先，我使用for循环进行最里层min函数的计算，随后进行最外层min函数的计算。最后得到t(x)后再根据博客开头给出的公式分通道计算JB(x)、JG(x)和JR(x)，然后合并三通道即可。代码实现如下：

    int WindowSize = 9;//定义窗口大小

    Mat t = Mat::zeros(img.rows, img.cols, CV_32FC1);//padding
    Mat tCopy;

    t.copyTo(tCopy);
    //先把通道/A最小值放入tCopy
    for(int i=0; i<img.rows; i++)
    {
        for(int j=0; j<img.cols; j++)
        {
            tCopy.at<float>(i, j) = min
            (
                  min((float)(img.at<Vec<uint8_t, 3>>(i, j)[0])/AB, (float)(img.at<Vec<uint8_t, 3>>(i, j)[1])/AG),
                  min((float)(img.at<Vec<uint8_t, 3>>(i, j)[0])/AB, (float)(img.at<Vec<uint8_t, 3>>(i, j)[2])/AR)
            );
        }
    }

    int m;
    float temp[9] = {0};
    float u = 0.75;//修正系数

    //计算窗口最小值
    for(int i=1; i<img.rows-1; i++)
    {
        for(int j=1; j<img.cols-1; j++)
        {
            for(m=0;m<9;m++)
            {
                temp[m] = tCopy.at<float>(i, j);
            }
            t.at<float>(i, j) = 1 - u * (*min_element(temp, temp + WindowSize));//找到窗口中的最小值

            if(t.at<float>(i, j)<0.1)//对t值进行处理以避免图像整体向白场过度
                t.at<float>(i, j) = 0.1;
        }
    }

    Mat J;
    Mat JB = Mat::zeros(img.rows, img.cols, CV_8UC1);
    Mat JG = Mat::zeros(img.rows, img.cols, CV_8UC1);
    Mat JR = Mat::zeros(img.rows, img.cols, CV_8UC1);

    for(int i=0; i<img.rows; i++)
    {
        for(int j=0; j<img.cols; j++)
        {
            JB.at<uint8_t>(i, j) =
                    1.0 * (img.at<Vec<uint8_t, 3>>(i, j)[0] - JB.at<uint8_t>(i, j)) / t.at<float>(i, j) + AB;
            JG.at<uint8_t>(i, j) =
                    1.0 * (img.at<Vec<uint8_t, 3>>(i, j)[1] - JG.at<uint8_t>(i, j)) / t.at<float>(i, j) + AG;
            JR.at<uint8_t>(i, j) =
                    1.0 * (img.at<Vec<uint8_t, 3>>(i, j)[2] - JR.at<uint8_t>(i, j)) / t.at<float>(i, j) + AR;
            //printf("%d  ",JG.at<uint8_t>(i, j));
        }
    }

    Mat channels[3];
    channels[0] = JB;
    channels[1] = JG;
    channels[2] = JR;

    merge(channels, 3, J);//通道合并

    return J;

效果展示

before：

after：

感觉效果雀食一般般，但是学到知识最重要^ _ ^