图像处理中任意核卷积(matlab中conv2函数)的快速实现。 - c++编程基础

　　卷积其实是图像处理中最基本的操作，我们常见的一些算法比如：均值模糊、高斯模糊、锐化、Sobel、拉普拉斯、prewitt边缘检测等等一些和领域相关的算法，都可以通过卷积算法实现。只不过由于这些算法的卷积矩阵的特殊性，一般不会直接实现它，而是通过一些优化的手段让计算量变小。但是有些情况下卷积矩阵的元素值无甚规律或者有特殊要求，无法通过常规手段优化，这个时候只能通过原始的方式实现。因此，如何快速的实现图像的任意卷积矩阵操作也有必要做适当的研究。

目前，通过友人共享或自己搜索找到的一片关于任意核算法优化的文章有： Reshuf?ing: A Fast Algorithm for Filtering with Arbitrary Kernels，改文章称能够提高原始程序速度的40%左右，但是原始的程序是如何写的也还不明白。

在matlab中有几个函数都与图像卷积有关，比如imfilter就可以实现卷积，或者 conv2也行，他们的速度都是相当快的，比如3000*3000的灰度图，卷积矩阵大小为15*15，在I5的CPU上运行时间只要170ms左右，相当的给力。

在Celery的博客中，也提到了他的优化后的conv2和matlab相当甚至快于matlab，详见http://blog.csdn.net/celerychen2009/article/details/38852105。

由于matlab的代码中使用到了IPL库进行加速，目前我写的Conv2函数还无法做到和其相当，对于任何核速度约为matlab的一半。

简单的记录下我在做卷积过程中用到的优化吧。

原始的卷积的实现需要四重循环，简单的表达如下：

for (Y = 0; Y < Height; Y++)
{
    for (X = 0; X < Width; X++)
    {
        Index = .....;
        Sum = 0;
        for (XX = 0; XX < ConvW; XX++)
        {
            for (YY = 0; YY < ConvH; YY++)
            {
                Index1 = ..... ;
                Index2 = ..... ;
                Sum += Conv[Index1] * Pixel[Index2];
            }
        }
        Dest[Index] = Sum / Weight;
    }
}

　　当卷积矩阵较大时，计算量将会很大，而且由于程序中的内存访问很频繁，cache miss现象比较严重，因此效率极为低下。

我的优化方法主要包括以下几个方面：

一：使用SSE进行乘法计算，由于SSE可以一次性进行4个单精度浮点数的计算，因此可以有明显的速度提升。

二：通过适当的处理方式，对每个取样点周边的卷积矩阵内的元素进行集中，使得每移动一个像素点不会需要从内存中进行大量的搜索工作。

具体来说实现过程如下：

1、为了使用SSE的优势，首先将卷积矩阵进行调整，调整卷积矩阵一行的元素个数，使其为不小于原始值的4的整数倍，并且让新的卷积矩阵的内存布局符合SSE相关函数的16字节对齐的要求。

　　实现代码如下：

float *Conv16 = (float *)_mm_malloc(PadConvLine * ConvH * sizeof(float), 16);                        //    保存16字节对齐的卷积矩阵，以方便使用SSE
                                                
for(Y = 0; Y < ConvH; Y++) 
{
    memcpy (Conv16 + Y * PadConvLine, Conv->Data.F + Y * ConvW , ConvW * sizeof(float));            //    复制卷积矩阵的数据
    memset(Conv16 + Y * PadConvLine + ConvW, 0, (PadConvLine - ConvW) * sizeof(float));                //    把冗余部分的卷积数据设置为0
}

其中PadConvLine = Pad4(ConvW) 以及Pad4的原型为： #define Pad4(bits) (((bits) + 3) / 4 * 4)；

注意_mm_malloc函数分配的内存中的值是随机值，对于扩展的部分一定要填充0，否则就会破坏卷积的结果。

那么如果我们也同时获得了需要被卷积的部分数据的话（卷积核肯定和卷积矩阵一样大小，且也应该是16字节对齐的），可以用如下的SSE的代码进行乘法计算：

float MultiplySSE(float *Kernel, float *Conv, int Length)
{
    int Block;  
    const float *Data;                        // 将SSE变量上的多个数值合并时所用指针.
    float Sum = 0;
    if (Length > 16)                        //    可以进行四次SSE计算，测试表明，这个还是快些的    
    {
        const int BlockWidth = 4 * 4;        // 块宽. SSE寄存器能一次处理4个float，然后循环展开4次.
        Block = Length / BlockWidth;        // 块数.    
        float *KernelP = Kernel, *ConvP = Conv;                // SSE批量处理时所用的指针.
        
        __m128 Sum0 = _mm_setzero_ps();         // 求和变量。SSE赋初值0
        __m128 Sum1 = _mm_setzero_ps();
        __m128 Sum2 = _mm_setzero_ps();
        __m128 Sum3 = _mm_setzero_ps();

        for(int I = 0; I < Block; I++)
        {
            Sum0 = _mm_add_ps(Sum0, _mm_mul_ps(_mm_load_ps(KernelP), _mm_load_ps(ConvP)));                    // SSE单精浮点紧缩加法
            Sum1 = _mm_add_ps(Sum1, _mm_mul_ps(_mm_load_ps(KernelP + 4), _mm_load_ps(ConvP + 4)));
            Sum2 = _mm_add_ps(Sum2, _mm_mul_ps(_mm_load_ps(KernelP + 8), _mm_load_ps(ConvP + 8)));
            Sum3 = _mm_add_ps(Sum3, _mm_mul_ps(_mm_load_ps(KernelP + 12), _mm_load_ps(ConvP + 12)));
            KernelP += BlockWidth;
            ConvP += BlockWidth;
        }
        
        Sum0 = _mm_add_ps(Sum0, Sum1);    // 两两合并(0~1).
        Sum2 = _mm_add_ps(Sum2, Sum3);    // 两两合并(2~3).
        Sum0 = _mm_add_ps(Sum0, Sum2);    // 两两合并(0~2).

        Data = (const float *)&Sum0;
        Sum = Data[0] + Data[1] + Data[2] + Data[3];

        Length = Length - Block * BlockWidth;            // 剩余数量.
    }
    if (Length != 0)
    {
        const int BlockWidth = 4;                        //    程序已经保证了数量必然是4的倍数 
        Block = Length / BlockWidth;        
        float *KernelP = Kernel, *ConvP = Conv;                
        __m128 Sum0 = _mm_setzero_ps();        

        for(int I = 0; I < Block; I++)
        {
            Sum0 = _mm_add_ps(Sum0, _mm_mul_ps(_mm_load_ps(KernelP), _mm_load_ps(ConvP)));        
            KernelP += BlockWidth;
            ConvP += BlockWidth;
        }

        Data = (const float *)&Sum0;
        Sum += Data[0] + Data[1] + Data[2] + Data[3];
    }
    return Sum;
}

　　当卷积矩阵（扩充后）的元素数量大于16时，我们采用了4路并行的SSE乘法实现，我在I3的CPU上测试时，2路SSE和4路SSE已经没有啥大的区别了，而在I5的CPU上则4路还是有较为明显的提高，因此采用4路SSE同时运行。当然1路SSE肯定还是比2路慢。另外，如果元素的数量少于16或者大于16但不能被16整除，那么余下的部分由于先前的扩充，剩余元素数量也肯定是4的倍数，因此可以用单路的SSE实现。这也是编码上的技巧。

2、前面提到了需要被卷积的部分数据，这部分如何快速的获取呢。观察最原始的4重循环，其内部的2重即为获取需要被卷积的部分，但是这里其实有很多问题。第一：由于卷积取样时必然有部分取样点的坐标在原始图像的有效范围外，因此必须进行判断，耗时。第二：同样为了使用SSE,也必须把取样的数据放在和扩充的卷积矩阵一样大小的内存中。这里我先贴出我的代码在进行解释具体的实现：

IS_RET __stdcall Conv2(TImage *Src, TMatrix *Conv, TImage *Dest, EdgeMode Edge)
{
    if (Src == NULL || Dest == NULL || Conv == NULL) return IS_RET_ERR_PARA;
    if (Src->Width != Dest->Width || Src->Height != Dest->Height || Src->BitCount != Dest->BitCount || Src->Stride != Dest->Stride) return IS_RET_ERR_PARA;
    if (Src->Scan0 == NULL || Dest->Scan0 == NULL || Conv->Data.F == NULL) return IS_RET_ERR_MEM;
    if (Conv->Width < 1 || Conv->Height < 1) return IS_RET_ERR_PARA;
    
    int Width = Src->Width, Height = Src->Height, Stride = Src->Stride;
    int ConvW = Conv->Width, ConvH = Conv->Height;
    unsigned char *PtSrc = Src->Scan0, *PtDest = Dest->Scan0;


    if (Src->BitCount == 24)
    {

    }
    else
    {
        int Left = ConvW / 2, Top = ConvH / 2, Right = ConvW - Left - 1, Bottom = ConvH - Top - 1, ExpHeight = Height + ConvH - 1;        //    注意核中心那个元素不用扩展，比如核的宽度为3，则只要左右各扩展一个像素就可以了
        int PadConvLine = Pad4(ConvW), Length = PadConvLin

图像处理中任意核卷积(matlab中conv2函数)的快速实现。(一)