【短道速滑十】非局部均值滤波的指令集优化和加速（针对5*5的搜索特例，可达到单核1080P灰度图 28ms/帧的速度）。

非局部均值滤波(Non Local Means）作为三大最常提起来的去燥和滤波算法之一（双边滤波、非局部均值、BM3D），也是有着很多的论文作为研究和比较的对象，但是也是有着致命的缺点，速度慢，严重的影响了算法的应用范围。目前在已有的文献中尚未看到在不对算法的本质原理上做更改的情况下，能取得实时的效果，本文呢，也不求得到这个目的，只是对现有的开放的资源上来取得更进一步的提升。

　　标准的NL-Means算法中，一般有三个参数，搜索半径SearchRadius，块半径PatchRadius，以及一个决定平滑程度的高斯函数参数Delta。在百度上能够搜索到的大部分文章所描述的提速算法都是使用积分图来提升NL-Means的速度，这也是目前来说唯一比较靠谱的优化技术，通过积分图，可以做到算法和块半径PatchRadius的大小基本无关，和Delta也无关，和SearchRadius成平方关系。

　　因此，在我们很多的严重的噪音图像中，SearchRadius至少需要取到7以上（涉及15*15= 225个领域范围）才有明显的效果，因此，这就相当于要计算225次全图的某种计算，即使是每次这种计算只需要1ms（通常，任何图像处理算法无法超越同等内存大小的memcpy的大小的），也需要225ms的，因此，确实比较感尴。

　　通过多线程方式可以适当对这个过程进行加速，毕竟每个像素点的处理相对来说还是独立的，但是，这个加速也收到物理核心的限制，就是8核的机器，满利用，也无法达到8倍的加速效果。

话说回来，这么大的计算量，用GPU也都是很吃力的。

　　目前使用积分图来记性NL-Means算法的比较好的文章也是大家常看的还是这一篇：

　　　　　非局部均值滤波（NL-means）算法的积分图加速原理与C++实现

　　其实积分图一直有个问题，可能很多搞图像的人都没有注意到，或者说这个问题可能对某些算法的影响还不是很大，不足有对大家注意到或者关注到，即积分图存在着一下2个方面的问题和缺点：

　　1、当图像较大时，积分图无法使用int类型来保存，我们必须选择能够容纳更大数据范围的数据类型来存储。

　　如果我们保存的是一副字节图像的积分图，考虑极端情况，每个图像的值都是255，则积分图像最多只可保存uint.Maxvalue / 255 = 16843009个像素，大约4100*4100大小的图像，一般来说，这个规模对于实际的应用来说是足够了。

　　但是我们在实际中用到的很多情况，不是直接的图像积分图，还常用的有平方积分图，即求图像平方值后的积分图，这个时候每一个元素的最大值达到了65025，极限情况下uint类型只可保存uint.Maxvalue / 255 =66051个元素的总和，只大概是指256*256个图像的大小了，已经远远的无法满足实际的应用需求。

　　因此，为了实现这个要求，我们必须选择能够容纳更大数据范围的数据类型，这里有三个选择：long long(int64) / float / double

　　第一个long long类型，即64位的整形数据，这个数据的表达范围已经完全够我们在图像中使用积分图了，而且保存的数据是非常准确的，特别是对于图像方面的数据来说（都是整形的），但是有个致命的问题：速度相当相当的慢，特别是同样的计算和int类型比较的话，那真的不是一个档次上的。我一直不太理解，现在大部分都是64位系统，为什么对64位的数据的支持还是这么的弱。而且我们看大部分指令集优化的函数对64位整形的支持都比较少。因此，非常不建议使用这个类型。

　　第二个float类型。如果使用这个类型，保存的数据范围是没有什么大的问题的，我们在网络上看到的文章大部分也是使用这个类型来保存结果的。但是，我在实践中多次遇到用float类型得不到正确的结果的问题，后来发现核心的原因是float的计算精度严重不足，特别是对于积分图这种连续的加法的计算，累计误差会越来越严重，当计算量足够大时，就会出现明显的误差瑕疵。因此，这个数据类型从本质上来说，对积分图是不够安全的。

　　关于这一点，实际上已经有不少作者注意到了，我在博文：SSE图像算法优化系列十四：局部均方差及局部平方差算法的优化中也有提及。

　　第三个是double类型，这个类型也是64位的，因此，数据范围毫无问题，计算精度经过测试，也是没有什么问题的，而且，编译器和指令集的支持和优化做的都还很不错，因此，个人认为这个数据类型是用来保存积分图最为合适的类型，但是有一个不友好的特点，计算速度慢，而且指令集对其能加速的空间有限。

　　2、积分图虽然能做到某些算法和参数无关，但是其并不是最佳的最速度

　　使用积分图技术，首先是要分配积分图占用的那部分额外的内存（而且是相当客观的内存），其次，积分图本身的计算也是需要一定时间的。还有，积分图的计算必须对边界部分做特别的判断，这个当参数较大时，计算量有的时候还是相当可观的。

　　还是回到我们的非局部均值滤波上吧。上面说了这么多，意思就是虽然非局部均值可以用积分图去优化，但是还是不是很好，那有没有更好的更快的实现呢，其实在我的下面两篇博客里就已经有了相关的技术。

　　一是： SSE图像算法优化系列十三：超高速BoxBlur算法的实现和优化（Opencv的速度的五倍）

二是： SSE图像算法优化系列十四：局部均方差及局部平方差算法的优化

　　因为非局部均值里使用的积分图是差的平方的积分图，因此我们只要对上述参考博文一里面的参与积分的数据稍微换一下即可，一个简单的C++代码如下所示：

//    SearchRadius    搜索的领域半径
//    PatchRadius      计算相似度的块的半径
//    Delta     高斯平滑的参数
int IM_NLM_Denoising(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int SearchRadius, int PatchRadius, float Delta)
{
    int Status = IM_STATUS_OK;
    int Channel = Stride / Width;
    if (Src == NULL)                        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))        return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1))                        return IM_STATUS_INVALIDPARAMETER;    

    int ExpandW = Width + SearchRadius + PatchRadius + SearchRadius + PatchRadius;
    int ExpandH = Height + SearchRadius + PatchRadius + SearchRadius + PatchRadius;
    int DiffW = Width + PatchRadius + PatchRadius;
    int DiffH = Height + PatchRadius + PatchRadius;

    float *Weight = (float *)calloc(Width * Height , sizeof(float));
    float *Sum = (float *)calloc(Width * Height , sizeof(float));
    int *ColValue = (int *)malloc((Width + PatchRadius + PatchRadius) * sizeof(int));
    unsigned char *Expand = (unsigned char *)malloc(ExpandH * ExpandW * sizeof(unsigned char));
    if ((Weight == NULL) || (Sum == NULL) || (ColValue == NULL) || (Expand == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    Status = IM_GetExpandImage(Src, Expand, Width, Height, Stride, ExpandW, SearchRadius + PatchRadius, SearchRadius + PatchRadius, SearchRadius + PatchRadius, SearchRadius + PatchRadius, IM_EDGE_MIRROR);
    if (Status != IM_STATUS_OK)    goto FreeMemory;

    int Area = (2 * PatchRadius + 1) * (2 * PatchRadius + 1);
    float Inv = -1.0f / Area / Delta / Delta;


    for (int YY = -SearchRadius; YY <= SearchRadius; YY++)
    {
        for (int XX = -SearchRadius; XX <= SearchRadius; XX++)
        {
            for (int Y = 0; Y < Height; Y++)
            {
                float *LinePS = Sum + Y * Width;
                float *LinePW = Weight + Y * Width;
                unsigned char *LinePE = Expand + (Y + YY + SearchRadius + PatchRadius) * ExpandW + XX + SearchRadius + PatchRadius;
                if (Y == 0)
                {
                    memset(ColValue, 0, DiffW * sizeof(int));
                    for (int Z = -PatchRadius; Z <= PatchRadius; Z++)
                    {
                        unsigned char *LineP1 = Expand + (Z + PatchRadius + YY + SearchRadius) * ExpandW + XX + SearchRadius;
                        unsigned char *LineP2 = Expand + (Z + PatchRadius + SearchRadius) * ExpandW + SearchRadius;
                        for (int X = 0; X < DiffW; X++)
                        {
                            int Value = LineP2[X] - LineP1[X];
                            ColValue[X] += Value * Value;
                        }
                    }
                }
                else
                {
                    unsigned char *LineOut1 = Expand + (Y - 1 + YY + SearchRadius) * ExpandW + XX + SearchRadius;
                    unsigned char *LineOut2 = Expand + (Y - 1 + SearchRadius) * ExpandW + SearchRadius;
                    unsigned char *LineIn1 = Expand + (Y + PatchRadius + PatchRadius + YY + SearchRadius) * ExpandW + XX + SearchRadius;
                    unsigned char *LineIn2 = Expand + (Y + PatchRadius + PatchRadius + SearchRadius) * ExpandW + SearchRadius;
                    for (int X = 0; X < DiffW; X++)
                    {
                        int Out = LineOut2[X] - LineOut1[X];
                        int In = LineIn2[X] - LineIn1[X];
                        ColValue[X] -= Out * Out - In * In;                                    //    更新列数据
                    }
                }
                int SumA = IM_SumofArray(ColValue, PatchRadius * 2 + 1);                //    处理每行第一个数据    
                float W = IM_Fexp(SumA * Inv);
                LinePW[0] += W;
                LinePS[0] += W * LinePE[0];
                for (int X = 1; X < Width; X++)
                {
                    SumA = SumA - ColValue[X - 1] + ColValue[X + PatchRadius + PatchRadius];
                    float W = IM_Fexp(SumA * Inv);
                    LinePW[X] += W;
                    LinePS[X] += W * LinePE[X];
                }
            }
        }
    }
    for (int Y = 0; Y < Height; Y++)
    {
        int Index = Y * Width;
        unsigned char *LinePD = Dest + Y * Stride;
        for (int X = 0; X < Width; X++)
        {
            LinePD[X] = Sum[Index + X] / (Weight[Index + X]);
        }
    }
FreeMemory:
    if (Expand != NULL)            free(Expand);
    if (Weight != NULL)            free(Weight);
    if (Sum != NULL)               free(Sum);
    if (ColValue != NULL)         free(ColValue);
    return Status;
}

　　对于常用的搜索半径为10，块半径为3, 500 * 500的灰度图耗时大概是500ms，相当的慢啊。

　　通过我一系列博文里的资料，可以知道上面的循环内的代码其实是很容易进行指令集优化的，基本上我那个方框模糊的优化是同一个技巧，经过指令集优化后，500*500的灰度图的耗时大概在200ms左右，如果加上线程技术，可以优化到75ms左右，这个时间和非局部均值滤波（NL-means）算法的CUDA优化加速文章里提的CUDA的优化速度基本已经差不多了。

　　另外，在搜索半径较小时，一种可行的优化方式是进行行列分离的卷积，即先计算中心行的结果，然后已这个结果为原始数据，在计算列方向的卷积，注意不是同时计算中心列和中心行的累加，而是有顺序的，这样就可以利用到周边领域所有的信息，这个时候计算量就会大为下降，计算的耗时也可以明显提高。这种优化方式必须测试是否对去燥的效果有很大的影响，因为他有可能会产生较为明显的水平或垂直线条效果。

　　另外，如果搜索半径较大，还可以尝试在上述基础上再进行45度和135度两个方向的卷积，以便抵消这种线条效果，同样提速也还是很明显的。

　　作为一个特例，有些情况下我们可能需要搜索半径为2，块大小也为2的非局部均值滤波，这种尺寸的滤波对于高强度的高斯噪音是没有什么去燥效果的，但是对于小规模的噪音，比如一般视频会议中的噪音还是有较强的抑制作用，因此，也还是有应用场景的，但是这种场景对算法的速度提出了极高的要求，如果考虑流畅性，给这个算法的处理不易超过20ms，而现在的视频流越来越高清了，因此，对这个算法的优化处理就必须做的更到位。

　　针对这个特例，我又做了一些优化，首先，因为是小半径，所以可以使用行列分离的算法，即先计算如下区域的合成结果：

　　计算得到中间的结果后再通过中间结果计算下述区域的合成结果作为最终值：

　　这里面的优化技巧有很多很多。

　　我这里提几个想法供参考：　　

　　1、中心点C的处的权重不需要计算，因为必然为1。

　　2、在进行水平计算时，去除掉中线后，只有4个点了，PP/P/N/NN，我们可以通过一些组合手段把他们的权重以及权重乘以值的量一次性计算出来，这样就可以直接把中间值给搞定了，如此做的好处是，不用多次保存累加的权重值和乘以值，也就是说舍弃掉了前面代码里的Weight和Sum变量，这对算法速度的提高也是有极大的贡献的，垂直的计算也是一样的。

3、那个exp函数是个比较耗时的函数，在http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/中一个比较快速的函数，虽然精度一般，但是经过测试，可以满足本算法的需求，但是其计算量小很多。

　　4、去掉中心点后，无论在垂直或者水平方向都只有4个数据，这为SIMD指令的某些读取和优化提供了无尚的便利和方便，而有些数据真的是巧合加天工设计，比如水平方向处理时，我们需要一次性处理4个点（出去中间那个），而我们用_mm_loadl_epi64恰好可以读取8个字节，这个时候，读取的8个字节如下：

　　　　　　　　A0 　　A1 　　A2 　　A3 　　A4 　　A5 　　A6 　　A7

　　正好可以组成4对5个组合 :

　　　　　　　　　　A0　　A1　　A2　　A3　　A4

　　　　　　　　　　A1 　　A2 　 A3 　　A4 　　A5

　　　　　　　　　　A2 　　A3 　A4 　　A5 　　A6

　　　　　　　　　　A3 　　A4 　A5 　　A6 A7

　　如果多一个字节，都不好处理了（主要是不好读取），多么完美的事情啊。

　　通过多种优化方式后，对于常见的1080P视频（1920*1080），处理其Y分量（灰度）耗时能做到28ms了，如果使用2个线程或4个线程，则可以完全满足实时的需求，如果是720P的视频，则单核也能到20ms。