超高速BoxBlur算法的落实和优化,源码分享

    表明:本文全体算法的涉及到的优化均指在PC上扩充的,对于其他构架是不是相符未知,请自行试验。

     
SSE图像算法优化类别五:超高速指数模糊算法的贯彻和优化(一千0*一千0在100ms左右达成) 一文中,笔者曾经说过优化后的ExpBlur比博克斯Blur还要快,今年笔者相比的BoxBlur算法是透过积分图+SSE实现的,小编在09年别的二个博客账号上一度提供过一篇那一个稿子彩色图像高速模糊之懒惰算法,里面也介绍了一种高效的图像模糊算法,这些算法的推行时间基本也是和半径非亲非故的。在二零一两年的SSE优化学习之路上本身早就也思虑过将该算法使用SSE实现,但马上认为那几个算法逐像素同期逐行都以左右注重的(单纯的逐像素信赖算法笔者在指数模糊里有关联怎么样用SSE优化),是力所不及用SSE管理的,一贯没思量,直到日前有朋友提议有个别基于局部局方差的算法希望能提速,笔者又再一次触发灵感,终于将这种算法也促成的一声令下集完成,並且测量试验速度比积分图快二倍,比分析opencv中央博物院克斯Filter的贯彻并建议进一步加速的方案(源码分享)那篇文章的进程快3倍,比opencv的cvSmooth函数快5倍,在一台老旧的I3台式机上拍卖两千*两千的灰度图到达了6ms的快慢,本文分享该优化进程并提供灰度版本的优化代码供我们学习和争辩。

      BoxFilter,最优秀的一种世界操作,在比非常多的场地中都独具广阔的选拔,作为一个很基础的函数,其性质的三六九等也一贯影响着其他相关函数的性质,最卓绝莫如以往很好的EPF滤波器:GuideFilter。因而其优化的水准和水准是相当主要的,网络上有相当多相关的代码和博客对该算法实行教学和优化,提议了好多O(1)算法,但所谓的0(1)算法也可能有上下之分,0(1)只是表示实施时间和有些参数非亲非故,但自己的耗费时间如故有分别。相比较盛行的便是储存直方图法,其实这一个是不行不可取的,因为稍大的图就能够形成溢出,这里大家分析下
opencv的相关代码,看看她是何等实行优化的。

  在彩色图像高速模糊之懒惰算法一文附带的代码中(VB6的代码)是对准二十三个人的图像,为了探讨方便,大家先贴出8位灰度的C++的代码:

     
首先找到opencv的代码的职责,其在\opencv\sources\modules\imgproc\src\smooth.cpp中。

  1 /// <summary>
  2 /// 按照Tile方式进行数据的扩展,得到扩展后在原尺寸中的位置,以0为下标。
  3 /// </summary>
  4 int IM_GetMirrorPos(int Length, int Pos)
  5 {
  6     if (Pos < 0)
  7         return -Pos;
  8     else if (Pos >= Length)
  9         return Length + Length - Pos - 2;
 10     else    
 11         return Pos;
 12 }
 13 
 14 void FillLeftAndRight_Mirror_C(int *Array, int Length, int Radius)
 15 {
 16     for (int X = 0; X < Radius; X++)
 17     {
 18         Array[X] = Array[Radius + Radius - X];
 19         Array[Radius + Length + X] = Array[Radius + Length - X - 2];
 20     }
 21 }
 22 
 23 int SumofArray_C(int *Array, int Length)
 24 {
 25     int Sum = 0;
 26     for (int X = 0; X < Length; X++)
 27     {
 28         Sum += Array[X];
 29     }
 30     return Sum;
 31 }
 32 
 33 /// <summary>
 34 /// 将整数限幅到字节数据类型。
 35 /// </summary>
 36 inline unsigned char IM_ClampToByte(int Value)            //    现代PC还是这样直接写快些
 37 {
 38     if (Value < 0)
 39         return 0;
 40     else if (Value > 255)
 41         return 255;
 42     else
 43         return (unsigned char)Value;
 44     //return ((Value | ((signed int)(255 - Value) >> 31)) & ~((signed int)Value >> 31));
 45 }
 46 
 47 int IM_BoxBlur_C(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius)
 48 {
 49     int Channel = Stride / Width;
 50     if ((Src == NULL) || (Dest == NULL))                    return IM_STATUS_NULLREFRENCE;
 51     if ((Width <= 0) || (Height <= 0) || (Radius <= 0))        return IM_STATUS_INVALIDPARAMETER;
 52     if (Radius < 1)                                            return IM_STATUS_INVALIDPARAMETER;
 53     if ((Channel != 1) && (Channel != 3) && (Channel != 4))    return IM_STATUS_NOTSUPPORTED;
 54 
 55     Radius = IM_Min(IM_Min(Radius, Width - 1), Height - 1);        //    由于镜像的需求,要求半径不能大于宽度或高度-1的数据
 56 
 57     int SampleAmount = (2 * Radius + 1) * (2 * Radius + 1);
 58     float Inv = 1.0 / SampleAmount;
 59 
 60     int *ColValue = (int *)malloc((Width + Radius + Radius) * (Channel == 1 ? Channel : 4) * sizeof(int));
 61     int *ColOffset = (int *)malloc((Height + Radius + Radius) * sizeof(int));
 62     if ((ColValue == NULL) || (ColOffset == NULL))
 63     {
 64         if (ColValue != NULL)    free(ColValue);
 65         if (ColOffset != NULL)    free(ColOffset);
 66         return IM_STATUS_OUTOFMEMORY;
 67     }
 68     for (int Y = 0; Y < Height + Radius + Radius; Y++)
 69         ColOffset[Y] = IM_GetMirrorPos(Height, Y - Radius);
 70 
 71     if (Channel == 1)
 72     {
 73         for (int Y = 0; Y < Height; Y++)
 74         {
 75             unsigned char *LinePD = Dest + Y * Stride;
 76             if (Y == 0)
 77             {
 78                 memset(ColValue + Radius, 0, Width * sizeof(int));
 79                 for (int Z = -Radius; Z <= Radius; Z++)
 80                 {
 81                     unsigned char *LinePS = Src + ColOffset[Z + Radius] * Stride;
 82                     for (int X = 0; X < Width; X++)
 83                     {
 84                         ColValue[X + Radius] += LinePS[X];                                            //    更新列数据
 85                     }
 86                 }
 87             }
 88             else
 89             {
 90                 unsigned char *RowMoveOut = Src + ColOffset[Y - 1] * Stride;                //    即将减去的那一行的首地址    
 91                 unsigned char *RowMoveIn = Src + ColOffset[Y + Radius + Radius] * Stride;    //    即将加上的那一行的首地址
 92                 for (int X = 0; X < Width; X++)
 93                 {
 94                     ColValue[X + Radius] -= RowMoveOut[X] - RowMoveIn[X];                                            //    更新列数据
 95                 }
 96             }
 97             FillLeftAndRight_Mirror_C(ColValue, Width, Radius);                //    镜像填充左右数据
 98             int LastSum = SumofArray_C(ColValue, Radius * 2 + 1);                //    处理每行第一个数据                                
 99             LinePD[0] = IM_ClampToByte(LastSum * Inv);
100             for (int X = 1; X < Width; X++)
101             {
102                 int NewSum = LastSum - ColValue[X - 1] + ColValue[X + Radius + Radius];
103                 LinePD[X] = IM_ClampToByte(NewSum * Inv);
104                 LastSum = NewSum;
105             }
106         }
107     }
108     else if (Channel == 3)
109     {
110 
111     }
112     else if (Channel == 4)
113     {
114 
115     }
116     free(ColValue);
117     free(ColOffset);
118     return IM_STATUS_OK;
119 }

   

  在此之前没发掘到,就这样轻易的代码用C写后能生出速度也是很摄人心魄的,2000*两千的图能幸不辱命39ms,假设在编写翻译选项里勾选编写翻译器的“启用巩固指令集:流式管理SIMD 扩展 2 (/arch:SSE2)”,
则系统会对上述全数浮点计算的一些行使相关的SIMD指令优化,如下图所示:

     Box Filter
是一种队列可分别的滤波,因而,opencv也是如此管理的,先进行行方向的滤波,获得中间结果,然后再对中间结果进行列方向的拍卖,获得终极的结果。

                       
  图片 1

     opencv
行方向处理的相关代码如下:

  这年三千*3000的图能成就25ms,,基本上临近小编创新的OPENCV的代码的速度了。

template<typename T, typename ST>
struct RowSum :
        public BaseRowFilter
{
    RowSum( int _ksize, int _anchor ) :
        BaseRowFilter()
    {
        ksize = _ksize;
        anchor = _anchor;
    }

    virtual void operator()(const uchar* src, uchar* dst, int width, int cn)
    {
        const T* S = (const T*)src;
        ST* D = (ST*)dst;
        int i = 0, k, ksz_cn = ksize*cn;

        width = (width - 1)*cn;
        for( k = 0; k < cn; k++, S++, D++ )
        {
            ST s = 0;
            for( i = 0; i < ksz_cn; i += cn )
                s += S[i];
            D[0] = s;
            for( i = 0; i < width; i += cn )
            {
                s += S[i + ksz_cn] - S[i];
                D[i+cn] = s;
            }
        }
    }
};

  轻易的叙说下各函数的效率先。

  这些代码思量了两个通道以及八种数据类型的意况,为了深入分析方便大家重写下单通道时的主干代码。

  IM_GetMirrorPos函数首借使获得某八个岗位Pos根据镜像的格局管理时在Length动向的坐标,这里Pos可以为负值,这一个首固然用来得到前期的坐标偏移的。      

for(Z = 0, Value = 0; Z < Size; Z++)    
    Value += RowData[Z];
LinePD[0] = Value;

for(X = 1; X < Width; X ++)
{
    Value += RowData[X + Size - 1] - RowData[X - 1];    
    LinePD[X] = Value;               
}

  FillLeftAndRight_Mirror_C首倘若用来博取两侧镜像数据的(直接得到,不调用IM_GetMirrorPos函数),举个例子比方Array原始数据为
***abcdefgh*** (Length = 8, Radius =
3),则结果为dcbabcdefghgfe。

  上述代码中RowData是指对某一行像素进行扩大后的像素值,其调幅为Width

  SumofArray_C主纵然计量贰个数组的保有的要素的和。

  • Radius + Radius, Radius是值滤波的半径, 而代码中Size = 2 * Radius +
    1,LinePD即所谓的中档结果,思索数据类型能发挥的界定,必得采用int类型。

  IM_BoxBlur_C函数内部即为模糊的函数体,接纳的优化思路完全和随意半径中值滤波(增添至百分比滤波器)O(1)时间复杂度算法的规律、完成及作用是平等的。当半径为Koleos时,方框模糊的值等于以某点为主干,左右左右各扩张ENCORE像素的的限定内装有像素的综合,像素总个数为(2*R+1)*(2*昂Cora+1)个,当然大家也足以把她分成(2*R+1)列,每列有(2*凯雷德+1)个元素本例的优化措施大家正是把累加数据分为一列一列的,丰富利用重复音讯来完结速度提高。

     
对于每行第多少个点很粗大略,间接用for总计出行方向的钦命半径内的累加值。而其后全部一些,用前一个累计值加上新加入的点,然后去除掉移出的哪二个点,得到新的累加值。那些艺术在好多文献中都有谈起,并从未什么样异样之处,这里相当少说。

  我们定义一个Radius + Width + Radius的内存数据ColValue,用来保存列方向的累加数据,对于第一行数据,需要做特殊处理,也就是用原始的方式计算一行像素所有元素在列方向上的和,
详见78至于86行代码,当然这里只计算了中间Width范围内的数据,当不是第一行时,如下图左边所示,新的累加值只需减去移出的哪一行像素值同时加上移入的一行像素值,详见90到96
行代码。

  上面只计算了中间Width范围内的累加值,为了方便后续代码的编写以及使用SSE优化,下面的FillLeftAndRight_Mirror_C主要作用就是填充左边和右边分别填充数据,而且是按照镜像的方式进行填充。

   
 依据正规的妄想,在列方向的拍卖相应和行方向大同小异,只是管理的势头的两样、管理的数据源的两样以及尾声索要对计量结果多一个归一化的经过而已。事实也是如此,有十分多算法都以如此做的,不过出于CPU构架的一些原因(首假若cache
miss的增多),同样的算法沿列方向管理总是会比沿行方向慢七个程度,化解措施有好多,比如先对中级结果开展转置,然后再根据行方向的条条框框举办拍卖,管理完后在将数据转置回去。转置进度是有丰硕高效的管理方式的,借助于SSE以及Cache优化,能落到实处比原本两重for循环快4倍以上的功力。还也是有一种方法就正如opencv中列管理进程所示,这多亏下边就要入眼描述的。

       
在更新了上述累加值后,大家开始拍卖计算均值了,对于每行的首先个点,SumofArray_C总计了前2*Odyssey+
1个列的累加值,这几个累加值正是该点左近(2*R+1)*(2*ENVISION+1)个要素的积攒和,对于一行的别的像素,其实就类似于行方向列累加值的翻新,减去移出的加盟新进的列,如下图左侧图所示,102到104行即实现了该进程。

     
在opencv的代码中,并未直接沿列方向管理,而是继续本着行的偏向一行一行管理,笔者先贴出我要好翻译的关于纯C的代码举办表明:

           图片 2   
 图片 3

    for (Y = 0; Y < Size - 1; Y++)            //    注意没有最后一项哦                      
    {
        int *LinePS = (int *)(Sum->Data + ColPos[Y] * Sum->WidthStep);
        for(X = 0; X < Width; X++)    ColSum[X] += LinePS[X];
    }

    for (Y = 0; Y < Height; Y++)
    {
        unsigned char* LinePD    = Dest->Data + Y * Dest->WidthStep;    
        int *AddPos              = (int*)(Sum->Data + ColPos[Y + Size - 1] * Sum->WidthStep);
        int *SubPos              = (int*)(Sum->Data + ColPos[Y] * Sum->WidthStep);

        for(X = 0; X < Width; X++)
        {
            Value = ColSum[X] + AddPos[X];
            LinePD[X] = Value * Scale;                    
            ColSum[X] = Value - SubPos[X];
        }
    }
  原理基本上就是这样,这个算法占用的额外内存很明显很少,但是不支持In-Place操作。
  为了追求速度,我们把整个过程能用SSE优化的地方都用SSE优化。
  首先是第79至86行的数据,这个很容易优化:

for (int Z = -Radius; Z <= Radius; Z++)
{
    unsigned char *LinePS = Src + ColOffset[Z + Radius] * Stride;
    int BlockSize = 8, Block = Width / BlockSize;
    for (int X = 0; X < Block * BlockSize; X += BlockSize)
    {
        int *DestP = ColValue + X + Radius;
        __m128i Sample = _mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(LinePS + X)));
        _mm_storeu_si128((__m128i *)DestP, _mm_add_epi32(_mm_loadu_si128((__m128i *)DestP), _mm_cvtepi16_epi32(Sample)));
        _mm_storeu_si128((__m128i *)(DestP + 4), _mm_add_epi32(_mm_loadu_si128((__m128i *)(DestP + 4)), _mm_unpackhi_epi16(Sample, _mm_setzero_si128())));
    }
    //  处理剩余不能被SSE优化的数据
}

     
上述代码中定义了一个ColSum用于保存每行有些地方处于列方向上钦点半径内各中等成分的累加值,对于第一行,做特殊管理,别的行的各样成分的处理方式和BaseRowFilter里的管理情势很像,只是在书写上有所分歧,极其注意的是对第一行的增进时,最终贰个因素并不曾测算在内,那几个管理技术在底下的X循环里有展示,请大家留心回味下。

  用_mm_loadl_epi64加载8个字节数据到内部存款和储蓄器,并用_mm_cvtepu8_epi16将其改变为拾四个人的__m128i变量,前面再把低4位和高4位的多寡分别调换来三13人的,然后用_mm_add_epi32抬高,注意后边笔者改换函数用了二种不相同的章程,因为那边的数额相对都以正数,因而也是能够动用_mm_cvtepi16_epi32和_mm_unpackhi_epi16组合Zero实现。

   
 上述代码那样做的好处是,能卓有功效的缩减CPU的cache
miss,可是总的计算量是从未变动的,因而能管用的增高速度。

  再来看看92到95行代码,那个也很简单。

   
 针对PC,在opencv内部,其在列方向上选择了SSE进行进一步的优化,大家贴出其管理uchar类型的代码:

int BlockSize = 8, Block = Width / BlockSize;
__m128i Zero = _mm_setzero_si128();
for (int X = 0; X < Block * BlockSize; X += BlockSize)
{
    int *DestP = ColValue + X + Radius;
    __m128i MoveOut = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(RowMoveOut + X)), Zero);
    __m128i MoveIn = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(RowMoveIn + X)), Zero);
    __m128i Diff = _mm_sub_epi16(MoveIn, MoveOut);                        //    注意这个有负数也有正数的,有负数时转换为32位是不能用_mm_unpackxx_epi16体系的函数
    _mm_storeu_si128((__m128i *)DestP, _mm_add_epi32(_mm_loadu_si128((__m128i *)DestP), _mm_cvtepi16_epi32(Diff)));
    _mm_storeu_si128((__m128i *)(DestP + 4), _mm_add_epi32(_mm_loadu_si128((__m128i *)(DestP + 4)), _mm_cvtepi16_epi32(_mm_srli_si128(Diff, 8))));
}
//  处理剩余不能被SSE优化的数据

图片 4View Code

  这里也是贰回性管理8个像素,这里本人利用了别的一种转移本事来把8位调换为16人,不过后边的Diff因为有正有负,要转变为三11个人就不可能不运用_mm_cvtepi16_epi32来转变,不能用unpack种类组合函数来落到实处,因为unpack不会移动符号位,小编在此处吃了一点次亏掉。

  代码相当多,作者不怎么精简下(注意,精简后的是不可运维的,只是更清晰的表明了思路)。

  接着是FillLeftAndRight_Mirror_C函数的优化,改写如下:

virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
{
    int i;
    int* SUM;
    bool haveScale = scale != 1;
    double _scale = scale;
    if( width != (int)sum.size() )
    {
        sum.resize(width);
        sumCount = 0;
    }

    SUM = &sum[0];
    if( sumCount == 0 )
    {
        memset((void*)SUM, 0, width*sizeof(int));
        for( ; sumCount < ksize - 1; sumCount++, src++ )
        {
            const int* Sp = (const int*)src[0];
            i = 0;

            for( ; i <= width-4; i+=4 )
            {
                __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
                __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
                _mm_storeu_si128((__m128i*)(SUM+i),_mm_add_epi32(_sum, _sp));
            }
            for( ; i < width; i++ )
                SUM[i] += Sp[i];
        }
    }
    else
    {
        src += ksize-1;
    }

    for( ; count--; src++ )
    {
        const int* Sp = (const int*)src[0];
        const int* Sm = (const int*)src[1-ksize];
        uchar* D = (uchar*)dst;

        i = 0;
        const __m128 scale4 = _mm_set1_ps((float)_scale);
        for( ; i <= width-8; i+=8 )
        {
            __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
            __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));

            __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
                                         _mm_loadu_si128((const __m128i*)(Sp+i)));
            __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
                                          _mm_loadu_si128((const __m128i*)(Sp+i+4)));

            __m128i _s0T = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
            __m128i _s0T1 = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s01)));

            _s0T = _mm_packs_epi32(_s0T, _s0T1);

            _mm_storel_epi64((__m128i*)(D+i), _mm_packus_epi16(_s0T, _s0T));

            _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
            _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
        }
        for( ; i < width; i++ )
        {
            int s0 = SUM[i] + Sp[i];
            D[i] = saturate_cast<uchar>(s0*_scale);
            SUM[i] = s0 - Sm[i];
        }

        dst += dststep;
    }
}
void FillLeftAndRight_Mirror_SSE(int *Array, int Length, int Radius)
{
    int BlockSize = 4, Block = Radius / BlockSize;
    for (int X = 0; X < Block * BlockSize; X += BlockSize)
    {
        __m128i SrcV1 = _mm_loadu_si128((__m128i *)(Array + Radius + Radius - X - 3));
        __m128i SrcV2 = _mm_loadu_si128((__m128i *)(Array + Radius + Length - X - 5));
        _mm_storeu_si128((__m128i *)(Array + X), _mm_shuffle_epi32(SrcV1, _MM_SHUFFLE(0, 1, 2, 3)));
        _mm_storeu_si128((__m128i *)(Array + Radius + Length + X), _mm_shuffle_epi32(SrcV2, _MM_SHUFFLE(0, 1, 2, 3)));
    }
    //  处理剩余不能被SSE优化的数据
}

     
在行方向上,ColSum每种成分的更新互相之间是未曾另外涉及的,因而,借助于SSE能够达成一次性的七个要素的立异,何况上述代码还将第一行的新鲜总计也用SSE实现了,纵然那部分总括量是老大小的。

  镜像正是倒序的进度,直接动用SSE的shuffle函数很有利实现。

   
 在切切实实的SSE细节上,对于uchar类型,尽管中间结果是用int类型表明的,不过出于SSE没有整形的除法指令(是还是不是有?),由此地点借用了浮点数的乘法SSE指令完结,当然就多了_mm_cvtepi32_ps
以及 _mm_cvtps_epi32如此的类型转换的函数。假若有整形除法,那就能够更加好了。

  总计数组的丰裕和优化也便于。

   
 SSE的实现,无非也正是用_mm_loadu_si128加载数据,用_mm_add_epi32, _mm_mul_ps之类的函数举办着力的演算,用_mm_storeu_si128来保存数据,并未怎么特别复杂的一对,注意到思索到数码的布满性,不鲜明都是16字节对齐的,由此在加载和保存是亟需使用u方面包车型大巴函数,其实现在的_mm_loadu_si128和_mm_load_si128在管理速度上曾经没有极度鲜明的分歧了。

int SumofArray_SSE(int *Array, int Length)
{
    int BlockSize = 8, Block = Length / BlockSize;
    __m128i Sum1 = _mm_setzero_si128();
    __m128i Sum2 = _mm_setzero_si128();
    for (int X = 0; X < Block * BlockSize; X += BlockSize)
    {
        Sum1 = _mm_add_epi32(Sum1, _mm_loadu_si128((__m128i *)(Array + X + 0)));
        Sum2 = _mm_add_epi32(Sum2, _mm_loadu_si128((__m128i *)(Array + X + 4)));
    }
    int Sum = SumI32(_mm_add_epi32(Sum1, Sum2));
    //  处理剩余不能被SSE优化的数据
    return Sum;
}

     
注意到在各类SSE代码前边,总还只怕有一对C代码,那是因为大家实际上数据小幅并不一定是4的整几倍,未被SSE管理的一部分必得使用C完成,那实则对读者来讲是老大好的业务,因为我们能从这一部分C代码中搞通晓上述的SSE代码是为什么的。

  使用三个__m128i
变量首要是为着丰硕利用XMM寄放器,个中SumI32函数如下,首假使为了总计__m128i
内七个整数的累加值。

  以上正是opencv的BoxFilter达成的要害代码,要是读者去看具体细节,opencv还应该有针对手提式有线电话机上的Neon优化,其实那一个和SSE的意趣基本是毫无二致的。

//    Convert 4 32-bit integer values to 4 unsigned char values.
void _mm_storesi128_4char(__m128i Src, unsigned char *Dest)
{
    __m128i T = _mm_packs_epi32(Src, Src);
    T = _mm_packus_epi16(T, T);
    *((int*)Dest) = _mm_cvtsi128_si32(T);
}

  那是否还会有革新的空间吗,从自己的体味中,在对垂直方向的拍卖上,应该未有了,不过程度方向呢,
SSE有未有用武之地,经过自己的想想自身感到依旧有的。

  代码不解释。

  在行方向的企图中,那些for循环是人命关天的揣度。

  最终正是100到104行的代码的改变。

for(X = 1; X < Width; X ++)
{
    Value += RowData[X + Size - 1] - RowData[X - 1];    
    LinePD[X] = Value;               
}
int BlockSize = 4, Block = (Width - 1) / BlockSize;
__m128i OldSum = _mm_set1_epi32(LastSum);
__m128 Inv128 = _mm_set1_ps(Inv);
for (int X = 1; X < Block * BlockSize + 1; X += BlockSize)
{
    __m128i ColValueOut = _mm_loadu_si128((__m128i *)(ColValue + X - 1));
    __m128i ColValueIn = _mm_loadu_si128((__m128i *)(ColValue + X + Radius + Radius));
    __m128i ColValueDiff = _mm_sub_epi32(ColValueIn, ColValueOut);                            //    P3 P2 P1 P0                                                
    __m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, 4));        //    P3+P2 P2+P1 P1+P0 P0
    __m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, 8));                 //    P3+P2+P1+P0 P2+P1+P0 P1+P0 P0
    __m128i NewSum = _mm_add_epi32(OldSum, Value);
    OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(3, 3, 3, 3));                              //    重新赋值为最新值
    __m128 Mean = _mm_mul_ps(_mm_cvtepi32_ps(NewSum), Inv128);
    _mm_storesi128_4char(_mm_cvtps_epi32(Mean), LinePD + X);
}

     
 自身认为尽管这里的操作是内外正视的,全局不能够并行化,不过观察这一行代码:

  在此以前感到这些算法难以使用SSE优化,重要难题就在那边,但是在上学了Opencv的积分图手艺时,这几个难点就一蹴而就了,进一步深入分析还发掘Opencv的代码写的并不到家,还也可能有革新的半空中,见上述代码,使用四次对同一数据移位就能够获取充分,由P3
P2 P1 P0变为P3+P2+P1+P0 P2+P1+P0 P1+P0 P0。大家有个别解析一下。          
  

Value += RowData[X + Size - 1] - RowData[X - 1];    

  __m128i ColValueDiff =
_mm_sub_epi32(ColValueIn, ColValueOut);
那句代码获得了移入和移出种类的差值,我们计为;  P3 P2 P1 P0
(高位到未有)由于算法的增进性情,假使说OldSum的原始值为A3 A3 A3
A3,那么新的NewSum就应有是:

      其中RowData[X + Size – 1] –
RowData[X – 1];
并非左右信赖的,是足以并行的,因而即使提前急忙的测算出全部的差值,那么提速的空间还非常大,即要求超前总括出下边包车型地铁函数:

    A3+P3+P2+P1+P0  A3+P2+P1+P0  A3+P1+P0  A3+P0;

 for(X = 0; X < (Width - 1); X++)
     Diff[X] = AddPos[X] - SubPos[X];
__m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, 4)); 这句的_mm_slli_si128得到中间结果 P2 P1 P0 0, 再和P3 P2 P1 P0相加得到

    Value_Temp =  P3+P2   P2+P1  P1+P0   P0

__m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, 8));这句的_mm_slli_si128得到中间结果P1+P0 P0 0 0,再和P3+P2 P2+P1 P1+P0  P0相加得到:

    Value = P3+P2+P1+P0   P2+P1+P0   P1+P0   P0
好简单的过程。

  OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(3, 3, 3, 3)); 这一句为什么要这样写,恐怕还是读者自己体会比较好,似乎不好用文字表达。

  这里用SSE达成则非常轻巧了:

   最终三个_mm_storesi128_4char是自家自身定义的二个将1个__m128i里面包车型地铁4个34个人整数调换为字节数并保留到内存中的函数,详见附属类小部件代码。

        unsigned char *AddPos = RowData + Size * Channel;
        unsigned char *SubPos = RowData;
        X = 0;                    //    注意这个赋值在下面的循环外部,这可以避免当Width<8时第二个for循环循环变量未初始化            
        __m128i Zero = _mm_setzero_si128();
        for (; X <= (Width - 1) * Channel - 8; X += 8)
        {
            __m128i Add = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i const *)(AddPos + X)), Zero);        
            __m128i Sub = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i const *)(SubPos + X)), Zero);        
            _mm_store_si128((__m128i *)(Diff + X + 0), _mm_sub_epi32(_mm_unpacklo_epi16(Add, Zero), _mm_unpacklo_epi16(Sub, Zero)));        //    由于采用了_aligned_malloc函数分配内存,可是使用_mm_store_si128
            _mm_store_si128((__m128i *)(Diff + X + 4), _mm_sub_epi32(_mm_unpackhi_epi16(Add, Zero), _mm_unpackhi_epi16(Sub, Zero)));
        }
        for(; X < (Width - 1) * Channel; X++)
            Diff[X] = AddPos[X] - SubPos[X];

  至于29个人图像的优化,是个相比较难堪的情境,因为SSE叁遍性管理4个三十一人数,而贰16个人各类像素唯有3个轻重,这种情景相似依然把她恢弘到4个轻重,举个例子说ColValue就足以改成4通道的,然后储存和也亟需管理成4坦途的,这本来必要一定的小巧淫技,这里自身不想多谈,留点东西给自身。当然也得以怀念先把贰十二个人分解到3个灰度内存,然后选拔灰度的算法举办测算,前边在合成到22人,那些解释和合成的进度也是能够运用SSE加快的,假使您结合多线程,还足以把3个灰度进度的拍卖并行化,那样除了多占用内部存款和储蓄器外,速度比至二级管理贰十四个人要快(因为间接处清理计算法无法并行的,前后信任的来头)。别的正是最后在计算列积攒求平均值的经过也变得特别自然了,不会出现灰度这种要在__mm128i在那之中实行增添的进度,而是径直得八个SSE变量累加。

  和列方向的SSE管理代码分化的是,这里加载的是uchar数据,因而加载的函数就天壤之别,管理的主意也会有分别,上边多少个SSE函数各位查查MSDN就能够知道其意思,依然很有深意的。

  还说一点,未来半数以上的CPU都扶助AVX256了,还能利用AVX进一步加快,就如代码该起来亦非很难,有乐趣的爱人可以和谐尝试。

  经过这么的拍卖,经过测量检验发现,速度能够比opencv的本子在升高百分之二十,也是十分的悲喜。

  能够说,方今以此速度已经大约达到了CPU的巅峰了,不过测量检验过IPP的速度,仿佛比这一个还要快点,不免除他利用了AVX,也不清除他选择多核的能源。

     
再有的优化大概提速有限,譬喻把剩余的一些for循环内部分成四路等等。

  这几个的优化对于BoxBlur来说是不可或缺的一步,不过更注重的是她能够动用在大多场馆,例如图像的一些均方差总计,也得以使用类似的能力拓宽加快,两幅图像大的有些平方差也是能够那样优化的,后续作者会轻易的谈下他在那方面加速的施用。

   
 在本人的I5台式机中,使用上述算法对1024*1024的三通道彩色图像进行半径为5的测量试验,举行玖拾捌回的测算纯算法部分的耗费时间为800ms,假设是纯C版本大概为1800ms;当半径为200时,SSE版本约为950ms,
C版本约壹玖零零ms;当半径为400时,SSE版本约为一千ms,
C版本约2100ms;可见,半径增大,耗费时间稍有增添,这重大是出于算法中有一部分起初化的乘除和半径有关,但是这几个都以不重大的。

  源代码下载:https://files.cnblogs.com/files/Imageshop/FastBlur.rar

     
在不选用八线程(就算本算法特别适合三十二线程总结),不行使GPU,只利用单线程\CPU举办测算的场合下,村办感到这段日子自己这一个代码是很难有质的超越的,从有些方面讲,代码中的用于总结的年华吞没的年华比从内部存款和储蓄器等待数据的流年恐怕还要短,而如同也未有越来越好的算法能更为回降内部存款和储蓄器数据的访谈量。

  彩色图工程:https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

     
本身在这里诚邀各位高手对近期本人优化的这么些代码举行更上一层楼的优化,希望高人不要客气。

图片 5

  运营分界面:

 

图片 6

  不佳意思,图太小,速度为0ms了。

  

图片 7

  本文的代码是针对性常用的图像数据开展的优化管理,在广大场合下,必要对int恐怕float类型进行拍卖,比方GuideFilter,如若读者知道了本文解读的代码的准则,更改代码以便适应他们则是一件特别轻便的作业,若是您不会,那么也请你不用给笔者留言,也许本身得以有偿提供,因为从精神上讲笔者快乐提供渔,并不是鱼,渔具已经送给您了,你却找笔者要鱼,那只可以…………..。

  

   
  本文源代码下载地址(VS二〇〇八编写制定):Box
Filter
 

  

图片 8

 

 

 

 

 

 ****************************我:
laviewpbt   时间: 二零一五.12.17    联系QQ:  33184777
转发请保留本行音信**********************