深度学习FPGA实现基础知识16(图像处理中任意核卷积(matlab中conv2函数)的快速实现)

绝地灬酷狼 2022-05-19 10:36 187阅读 0赞

**需求说明：深度学习FPGA实现知识储备**

**来自：[http://www.2cto.com/kf/201411/355943.html][http_www.2cto.com_kf_201411_355943.html]**

**整理来自：时间的诗**  
卷积其实是图像处理中最基本的操作，我们常见的一些算法比如：均值模糊、高斯模糊、锐化、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][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 = PadConvLine * ConvH;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            int X, Y, IndexD, IndexE, IndexK, ExpStride;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            float *CurKer, Inv, Sum = 
           
            0;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            unsigned 
           
            char *PtExp, *PtDest;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
            
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   TImage *Expand;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   IS_RET Ret = GetPadImage(Src, &Expand, Left, Right, Top, Bottom, Edge);                            
          
           
         
          
         
         
          
          
           
         
          
         
          
          
             
           
            // 得到扩展后的数据，可以提速和方便编程，但是多占用一份内存
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            if (Ret != IS_RET_OK) 
           
            return Ret;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   PtExp = Expand->Scan0; PtDest = Dest->Scan0; ExpStride = Expand->Stride;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            for (X = 
           
            0; X < ConvH * ConvW; X ++) Sum += Conv->Data.F[X];
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   Inv = (Sum == 
           
            0 ? 
           
            1: 
           
            1 / Sum);                                                                        
           
            // 如果卷积举证的和为0，则设置为1
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            float *Conv16 = (
           
            float *)_mm_malloc(PadConvLine * ConvH * 
           
            sizeof(
           
            float), 
           
            16);  
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
           
            // 保存16字节对齐的卷积矩阵，以方便使用SSE
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            float *Kernel = (
           
            float *)_mm_malloc(PadConvLine * ExpHeight * 
           
            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
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   }
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            for (Y = 
           
            0; Y < ExpHeight; Y++)
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   {
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                       IndexE = Y * ExpStride;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                       CurKer = Kernel + Y * PadConvLine;                        
           
            // 计算第一列所有像素将要取样的卷积核数据
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                       
           
            for (X = 
           
            0; X < ConvW; X++)
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                       {
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                           CurKer[X] = PtExp[IndexE++];
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                       }
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   }
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            for (X = 
           
            0 ; X < Width ; X ++)
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   {
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                       
           
            if (X != 
           
            0)                                                
           
            // 如果不是第一列，需要更新卷积核的数据
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                       {
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                           
           
            memcpy(Kernel, Kernel + 
           
            1, (PadConvLine * ExpHeight - 
           
            1) * 
           
            sizeof(
           
            float));    
           
            // 往前移动一个数据
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                           IndexK = ConvW - 
           
            1 ;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                           IndexE = IndexK + X;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                           
           
            for (Y = 
           
            0; Y < ExpHeight; Y++)
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                           {
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                               Kernel[IndexK] = PtExp[IndexE];        
           
            // 只要刷新下一个元素
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                               IndexK += PadConvLine;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                               IndexE += ExpStride;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                           }
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                       }
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                       
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                       CurKer = Kernel;    IndexD = X;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                       
           
            for (Y = 
           
            0; Y < Height; Y ++)                            
           
            // 沿列的方向进行更新
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                       {
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                           PtDest[IndexD] = Clamp((
           
            int)( MultiplySSE(Conv16, CurKer, Length) * Inv + 
           
            0.5));        
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
           
            // 直接把函数放在这里也没有啥提速的，注意改函数不会被内联的
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                           CurKer += PadConvLine;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                           IndexD += Stride;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                       }
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   }
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   _mm_free(Conv16);
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   _mm_free(Kernel);
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
                   FreeImage(Expand);
          
           
         
          
         
         
          
          
           
         
          
         
          
          
                   
           
            return IS_RET_OK;
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
               }
          
           
         
          
         
         
          
          
           
         
          
         
          
          
           
            
           }

对于第一个问题，解决的方式很简答，即用空间换时间，新建一副（Width + ConvW - 1, Height + ConvH -1)大小的图像，然后四周的ConvW及ConvH的像素用边缘的值或者边缘镜像的值填充，正中间的则用原来的图复制过来，这样操作后进行取样时不再原图取样，而在这福扩展的图中取样，就避免了坐标判断等if语句的跳转耗时了，上GetPadImage即实现了改功能。  
第二个问题则需要有一定的实现技巧，我们分配一块PadConvLine \* (Height + ConvH - 1) 大小的内存，然后计算原图第一列像素串联起来的需要卷积的部分的数据，这一部分代码如上述44-52行所示。有了这样的数据，如果需要计算第一列的卷积结果，则很简单了，每跳过一列则把被卷积的数据起点增加PadConvLine个元素，在调用上述MultiplySSE函数获得卷积结果。接着则计算第二列像素的卷积值，此时需要整体更新这一列像素串联起来的需要被卷积的数据，更新也很简单，就是把原来的数据整体向左移动一个像素，这个可以用memcpy快速实现，然后在填充入新进来的那个元素，就ok了，接着就是再次调用MultiplySSE函数，如此重复下去。  
经过编码测试，对于3000\*3000的灰度图，15\*15的核在I5的CPU上的测试平均结果为360ms,比matlab的慢了一半。  
最后说明一点，很多人都说用FFT可以快速的实现卷积，并且是O(1)的，我比较同意后半句，但是前面半句是绝对的有问题的，至少在核小于50\*50时，FFT实现的卷积不会比直接实现块。要知道FFT的计算量其实是很大的。

</div>

[http_www.2cto.com_kf_201411_355943.html]: http://www.2cto.com/kf/201411/355943.html
[http_blog.csdn.net_celerychen2009_article_details_38852105]: http://blog.csdn.net/celerychen2009/article/details/38852105