一种可实时处理 O(1)复杂度图像去雾算法的实现。

　　在我博文的一系列的文章，有不少算法都于去雾有关，比如限制对比度自适应直方图均衡化算法原理、实现及效果、局部自适应自动色阶/对比度算法在图像增强上的应用这两个增强算法都有一定的去雾能力，而最直接的就是《Single Image Haze Removal Using Dark Channel Prior》一文中图像去雾算法的原理、实现、效果及其他一文，描述了暗通道去雾这一state-of-the-art algorithms的过程和实现，虽几经优化，对于常用的视频1024*768大小的图片，算法处理部分还是需要70MS的时间（I3 笔记本CPU），因此，这一算法用于实时要求时还有一定的难度，并且优化后的算法基本无法并行，而可并行的算法重复计算大，由于不熟悉GPU方面的理念，不晓得使用不优化的算法靠GPU是否能有多大速度的提升。

为此，我一直在找寻相关的论文，这种找寻的踪迹一般就是看到一篇好论文--》看其参考文献--》再看参考文献的参考文献，这样循环下去。然后有某种机会或巧合，又看到一篇好论文，重复前面的过程，你就会发现很多交集，慢慢的就会有一些好运向你招手。话说我原本只看英文的文献，所以一直忽略了国内的文章，前几日，一个QQ朋友推荐了一篇清华大学的论文，下载后稍微看了下，觉得其描述的结果还是比较吸引人的，于是就实现了下，实时的效果应该说很不错，这里就简单的介绍并推荐给大家。

原始论文下载：基于单幅图像的快速去雾.pdf ，作者刘倩，陈茂银，周东华，感谢他们（她们）。

算法原理没有什么复杂的地方，其实说原理，还不如说经验或实验，因为论文中可以用理论来推导的公式确实不多。不过这也没关系，有用的东西就应该拿来用。

算法的执行流程直接贴用原图的文字来说明吧：

先挑点小瑕疵，比如步骤5中的L有个下标，而其他涉及到L的地方去没有，这叫前后不一致。论文中也还有一些其他的地方有些小错误，所以看啊，即使是清华的文章，在编辑审核这一块还是相当的不严谨。

我们先来看看算法，1、2、3、4步都没有什么说的，第五步是求大气透射率的过程，这里ρ是一个用来调节的参数，当ρ值越大时，结果图像整体越暗，去雾的效果更明显，ρ较小时，图像偏白，有明显的雾气。第六步求全局大气光A值，用了很简单的方式，即求原始图像的RGB所有像素分量的最大值（这个估计99%都为255了）何暗通道的最大值的平均值，并且注意到RGB三个通道用的A值都为同一个数字。这个和我在何凯明的文章的分析也有相似的地方。第七步用了标准的去雾模型来求结果值。

在来看看算法的效率问题，从算法的初步分析来看，算法的效率取决于第3步和第7步，第三步中，使用了均值模糊，目前已经存在大量的O(1)均值模糊算法，可是O(1)只能表示算法的执行速度和参数的大小无关系，并不表示算法就很快。比如基于积分图的模糊算法是广为认知的O(1)算法，但是他也存在很多问题，最严重的就是数据的溢出，当图像较大和偏白时，对图像积分图的累加和存在超出int.Maxvalue所能表达的范围的问题，解决办法就是积分图内的数据全部使用long类型表示，这将导致程序多占用Width*Height*4字节的大小的内存，且在32位系统还流行的情况进一步降低程序的速度（32位系统64位整数的计算速度要比32位整数慢)。积分图的另外一个问题就是计算积分图的过程难以并行化，因为一个像素的积分值是依赖于其前面一系列像素的相关结果值的。另外一种优化方式就是先计算行方向的平均值，然后再计算列方向的值。这种方式在同一行（列）内，算法依旧必行顺序执行，这也是因为前后影响的原因。但是不同行（列）之间的计算是没有任何关系，因此非常适合GPU这种可大规模并行计算的场合，但不适于CPU这种重量级的线程并发（反而会慢）。这种算法如果为了精度会需要一个和原图一样大小，占用字节Width*Height*4字节大小的的中转区用来保存中间计算的结果。在彩色图像高速模糊之懒惰算法一文中，我采用了另外一种处理方法，利用列直方图相关的技术，只需对每个循环的起始位置处的像素做特殊处理，其他位置的利用简单的一加一简即可获得累加和，从而快速的实现模糊，我实际的编码表明，这种方式比其他的方式都要快。但是有一个缺点，不适合于并行计算，不过在CPU上这个很有优势。

再观察下第七步。第七步存在两个需要优化的地方，第一，存在除法；第二，有浮点运算。如果直接编码必然会带来性能损失，但是，观察下在第七步的公式中，只有两个自变量，H(X)和L(X),并且自变量的取值都为[0,255]之间的整数，因此，如果事先建议一个查找表，由于这个查找表的计算量只有 256*256次，要远远的小于直接计算的次数，必然能提高程序的速度。256这个数字还有个好处，就是可以用移位来辅助计算查表表的下标，部分参考代码如下所示：

unsigned char * Table = (unsigned char  *) malloc (256*256*sizeof(unsigned char));
for (Y = 0; Y < 256; Y++)
{
    Index=Y<<8;
    for (X = 0; X < 256; X++)
    {
        Value = (Y-X) /(1-X*InvA);
        if (Value > 255)
            Value = 255;
        else if (Value < 0)
            Value = 0;
        Table[Index++]= Value;
    }
}

　　其中InvA = 1 / A;A为全局大气光值。 Value需为double类型的变量。

当我们需要计算F(x)时，查表的方式为 F(X)=Table[(H[X]<<8 )+ L[X]];

实际的效果表明，这样的方式对于1024*768的图，可以提速10ms。

那么对于其他步骤也有很多优化的注意事项，比如计算M(X)中所有元素的平均值M_av这一块，完全没有必要在开一个循环，而是可以在进行步骤3的时候同步进行，大家知道，循环楚了要计算循环体内部的东西外，还要有个循环计数器的更新的，何必浪费这个时间呢。

for (Y = 0, DarkPt = DarkChannel; Y < Height; Y++)
{
    ImgPt = Src + 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;                                        //    三通道的最小值
        Sum    +=    Min;                                        //  累积以方便后面求平均值
        ImgPt += 3;
        DarkPt++;
    }
}
Mean =(double)Sum /(Width*Height*255);

还有比如第6步中，分别求原图RGB三像素最大值以及安通道中的最大值的过程，传统的过程如下代码：

int Max1 =0 ;
for (Y = 0; Y < Height; Y++)
{
    ImgPt = Src + Y * Stride;
    for (X = 0; X < Width; X++)
    {
        if (Max1 < *(ImgPt)) Max1 = *(ImgPt);
        if (Max1 < *(ImgPt + 1)) Max1 = *(ImgPt + 1);
        if (Max1 < *(ImgPt + 2)) Max1 = *(ImgPt + 2);
        ImgPt += 3;
    }
}

　　特别是对于求原图的最大值，实际上很多情况下这个值都为255，因此如果Max1变量已经是255，则循环完全没有必要进行下去了，因此，如果改为下述代码，必然可以减少计算量：

for (Y = 0; Y < Height; Y++)
{
    ImgPt = Src + Y * Stride;
    for (X = 0; X < Width; X++)
    {
        if (Max1 < *(ImgPt)) Max1 = *(ImgPt);
        if (Max1 < *(ImgPt + 1)) Max1 = *(ImgPt + 1);
        if (Max1 < *(ImgPt + 2)) Max1 = *(ImgPt + 2);
        ImgPt += 3;
    }
    if (Max1==255) break;
}

注意，这个break语句必须放在Y循环中，如果放在X循环中，虽然提前退出循环的可能性会增加，但是判断的工作量带来的损失更多。
综合上述优化，我用C++写了个DLL，对于1024*768的图像在我的I3的机器上平均能达到18ms每副图像的计算速度，相当于56fps,只占用了单核的资源，考虑解码、显示等等其他过程所占用的时间，应该是能够靠CPU实现20fps的实时速度的。

在内存占用上，约需要> 3*Width*Height+256*256字节的空间（不包括图像本身的），如果用在连续的视频处理上，这部分内存就不需要频繁的分配和释放，可能也对速度的保证有好处。

程序下载地址： https://files.cnblogs.com/Imageshop/FastHazeRemovalTest.rar

　参数的选取：

　为了获得好的效果，该算法需要选择恰当的参数，我们为此做了一些测试。对于半径参数，我的个人建议是取值不要小于50或图像宽度和高度最大值的的1/20，比如对下面的图像（原图大小不是下图中的大小，这里是为了方便浏览缩小显示的），ρ取1.28时（对于上图中的去雾程序选择为128），半径取不同参数时的效果：