自个儿曾经说过优化后的ExpBlur比BoxBlur还要快,互连网上有很多有关的代码和博客对该算法实行教学和优化

    表达:本文全数算法的关联到的优化均指在PC上开始展览的,对于其他构架是或不是适宜未知,请自行试验。

     
SSE图像算法优化连串五:超高速指数模糊算法的完结和优化(10000*一千0在100ms左右落实) 一文中,我已经说过优化后的ExpBlur比BoxBlur还要快,那几个时候自身比较的BoxBlur算法是经过积分图+SSE完结的,笔者在09年其它3个博客账号上早已提供过一篇那些稿子彩色图像高速模糊之懒惰算法,里面也介绍了一种高效的图像模糊算法,这几个算法的实施时间基本也是和半径非亲非故的。在今年的SSE优化学习之路上本身一度也考虑过将该算法使用SSE完成,但立时认为这几个算法逐像素同时逐行都以上下正视的(单纯的逐像素重视算法笔者在指数模糊里有涉嫌怎样用SSE优化),是无力回天用SSE处理的,一向没考虑,直到日前有情侣提出有些基于局地局方差的算法希望能提速,笔者又重新触发灵感,终于将那种算法也落到实处的命令集达成,并且测试速度比积分图快二倍,比解析opencv中BoxFilter的兑现并提议进一步加快的方案(源码共享)这篇小说的快慢快3倍,比opencv的cvSmooth函数快5倍,在一台老旧的I3台式机上拍卖三千*贰仟的灰度图达到了6ms的快慢,本文分享该优化进度并提供灰度版本的优化代码供大家学习和座谈。

   
 那里的高斯模糊使用的是舆论《Recursive implementation of the Gaussian
filter》里描述的递归算法。

      BoxFilter,最经典的一种世界操作,在广大的场子中都具有广大的利用,作为多个很基础的函数,其质量的高低也直接影响着此外有关函数的本性,最特异莫如未来很好的EPF滤波器:GuideFilter。因而其优化的档次和品位是可怜关键的,互连网上有很多相关的代码和博客对该算法进行讲解和优化,提出了诸多O(1)算法,但所谓的0(1)算法也有高低之分,0(1)只是意味着执行时间和有个别参数无关,但笔者的耗费时间如故有分别。比较盛行的就是积累直方图法,其实这么些是老大不可取的,因为稍大的图就会造成溢出,那里大家分析下
opencv的连带代码,看看她是什么样开始展览优化的。

  在彩色图像高速模糊之懒惰算法一文附带的代码中(VB6的代码)是针对性2二位的图像,为了切磋方便,我们先贴出七个人灰度的C++的代码:

           
  图片 1

     
首先找到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 }

     
仔细考察和清楚上述公式,在forward进度中,n是递增的,由此,如若在拓展forward以前,把in数据先全部的赋值给w,然后式子(9a)就能够改为:

   

  以前没觉察到,就那样简单的代码用C写后能生出速度也是很诱人的,两千*两千的图能连成一气39ms,假如在编写翻译选项里勾选编写翻译器的“启用增强指令集:流式处理
SIMD 增加 2 (/arch:SSE2)”,
则系统会对上述全数浮点计算的片段使用相关的SIMD指令优化,如下图所示:

      w[n] = B w[n] +
(b1 w[n-1] + b2 w[n-2] +
b3 w[n-3]) / b0;          ———>    
(1a)

     Box Filter
是一种队列可分其他滤波,由此,opencv也是那样处理的,先举行行方向的滤波,获得中间结果,然后再对中级结果进行列方向的拍卖,获得终极的结果。

                       
  图片 2

 
  在backward进程中,n是递减的,由此在backward前,把w的数码完全的拷贝到out中,则式子(9b)变为:

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

  那一个时候3000*3000的图能做到25ms,,基本上接近小编创新的OPENCV的代码的快慢了。

     out[n] = B out[n] +
(b1 
out[n+1] + b2 out[n+2] +
b3 
out[n+3]) / b0 ;     <———     (1b)

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;
            }
        }
    }
};

  一句话来说述下各函数的效应先。

   
 从编制程序角度看来,backward中的拷贝是一点一滴没有须要的,因而 (1b)能够一直写为:

  那些代码考虑了八个通道以及多种数据类型的处境,为了分析方便大家重写下单通道时的中坚代码。

  IM_GetMirrorPos函数首借使得到某1个任务Pos根据镜像的方法处理时在Length动向的坐标,那里Pos能够为负值,这几个第1是用来取得中期的坐标偏移的。      

           w[n] = B w[n] +
(b1 w[n+1] + b2 w[n+2] +
b3 w[n+3]) / b0 ;               <———    
(1c)

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。

   
 
从进程上考虑,浮点除法非常的慢,由此,大家再一次定义b1 = b1 / b0, b2 = b2 /
b0, b3 = b3 / b0, 最后获得大家使用的递归公式:

  上述代码中RowData是指对某一行像素实行扩充后的像素值,其上涨幅度为Width

  SumofArray_C首若是计量3个数组的有着的因素的和。

           w[n] = B w[n] +
b1 w[n-1] + b2 w[n-2] +
b3 w[n-3];          ———>     (a)

  • Radius + Radius, Radius是值滤波的半径, 而代码中Size = 2 * Radius +
    1,LinePD即所谓的中间结果,考虑数据类型能发布的界定,必须选用int类型。

  IM_BoxBlur_C函数内部即为模糊的函数体,选用的优化思路完全和随机半径中值滤波(扩张至百分比滤波器)O(1)时间复杂度算法的法则、完毕及效用是一致的。当半径为昂科雷时,方框模糊的值等于以某点为着力,左右左右各扩大Murano像素的的范围内部存款和储蓄器有像素的汇总,像素总个数为(2*R+1)*(2*CRUISER+1)个,当然大家也足以把她分成(2*R+1)列,每列有(2*昂Cora+1)个因素本例的优化措施大家正是把累加数据分为一列一列的,丰裕利用重复消息来落成速度升高。

        w[n] = B w[n] +
b1 w[n+1] + b2 w[n+2] +
b3 w[n+3] ;             <———      (b)

     
对于每行首个点相当粗略,间接用for计算出行方向的内定半径内的累加值。而其后全体点,用前2个累计值加上新进入的点,然后去除掉移出的哪1个点,获得新的累加值。这些主目的在于广大文献中都有提及,并从未什么样异样之处,那里不多说。

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

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

   
上述公式是一维的高斯模糊总括办法,针对二维的图像,正确的做法正是先对各类图像行进行模糊处理得到中间结果,再对中等结果的每列进行模糊操作,最后获得二维的歪曲结果,当然也能够利用先列后行那样的做法。

   
 依据常规的思索,在列方向的拍卖相应和行方向完全相同,只是处理的趋向的两样、处理的数据源的不比以及最终索要对计量结果多1个归一化的历程而已。事实也是如此,有很多算法都以如此做的,不过出于CPU构架的局地缘由(首假设cache
miss的增多),同样的算法沿列方向处理总是会比沿行方向慢3个程度,化解措施有成都百货上千,例如先对中间结果开始展览转置,然后再依照行方向的平整实行拍卖,处理完后在将数据转置回去。转置进度是有特别急速的处理格局的,借助于SSE以及Cache优化,能促成比原来两重for循环快4倍以上的效益。还有一种方法就正如opencv中列处理进程所示,那就是上边将要重点描述的。

       
在更新了上述累加值后,大家开头拍卖总计均值了,对于每行的首个点,SumofArray_C计算了前2*路虎极光+
一个列的累加值,那几个累加值正是该点周边(2*R+1)*(2*奥迪Q7+1)个因素的累积和,对于一行的此外像素,其实就恍如于行方向列累加值的翻新,减去移出的插足新进的列,如下图右边图所示,102到104行即达成了该进程。

     
此外注意点便是,边缘像素的处理,大家见到在公式(a)或许(b)中,w[n]的结果个别凭借于前八个或者后多个成分的值,而对于边缘地点的值,这几个都以不在合理限定内的,常常的做法是镜像数据恐怕重新边缘数据,实践评释,由于是递归的算法,初步值的例外会将结果直接一而再下去,由此,区其他方法对边缘部分的结果要么有自然的震慑的,那里大家采取的归纳的重新边缘像素值的办法。

     
在opencv的代码中,并没有平昔沿列方向处理,而是继续本着行的势头一行一行处理,作者先贴出我要好翻译的有关纯C的代码举行表达:

           图片 3   
 图片 4

     
由于地点公式中的周密均为浮点类型,因而,总结一般也是对浮点实行的,也便是说一般需求先把图像数据转换为浮点,然后进行高斯模糊,在将结果转换为字节类型的图像,由此,上述进程能够独家用一下多少个函数完结:

    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优化的数据
}

               
CalcGaussCof          
//  总结高斯模糊中采取到的全面
      ConvertBGMurano8U2BGRAF      //  将字节数据转换为浮点数据 
      GaussBlurFromLeftToRight    //  水平方向的前向传来
      GaussBlurFromRightToLeft    //  水平方向的反向传播
      GaussBlurFromTopToBottom  
 //   垂直方向的前向传来
      GaussBlurFromBottomToTop  
 //   垂直方向的反向传播
      ConvertBGRAF2BG福睿斯8U      
 //   将结果转换为字节数据

     
上述代码中定义了叁个ColSum用于保存每行有些地点处于列方向上钦点半径内各中等成分的累加值,对于第二行,做特别处理,其余行的种种成分的处理格局和BaseRowFilter里的处理方式很像,只是在挥洒上有所分化,尤其注意的是对第②行的丰裕时,最终1个成分并从未测算在内,这么些处理技术在上面包车型大巴X循环里有反映,请我们精心回味下。

  用_mm_loadl_epi64加载九个字节数据到内部存储器,并用_mm_cvtepu8_epi16将其更换为拾三人的__m128i变量,前边再把低3个人和高3人的多寡分别转换来30个人的,然后用_mm_add_epi32添加,注意前面小编改换函数用了三种分歧的形式,因为这边的数额绝对都以正数,因而也是能够利用_mm_cvtepi16_epi32和_mm_unpackhi_epi16组合Zero实现。

   大家各样分析之。

   
 上述代码那样做的裨益是,能使得的缩减CPU的cache
miss,然而总的计算量是不曾改变的,因而能卓有功能的滋长速度。

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

     
 CalcGaussCof,这些很简单,就遵照杂谈中提议的总计公式依据用户输入的参数来计量,当然结合下方面大家关系的除法变乘法的优化,注意,为了继续的一些标注的标题,笔者讲上述公式(a)和(b)中的周密B写为b0了。

   
 针对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优化的数据
void CalcGaussCof(float Radius, float &B0, float &B1, float &B2, float &B3)
{
    float Q, B;
    if (Radius >= 2.5)
        Q = (double)(0.98711 * Radius - 0.96330);                            //    对应论文公式11b
    else if ((Radius >= 0.5) && (Radius < 2.5))
        Q = (double)(3.97156 - 4.14554 * sqrt(1 - 0.26891 * Radius));
    else
        Q = (double)0.1147705018520355224609375;

    B = 1.57825 + 2.44413 * Q + 1.4281 * Q * Q + 0.422205 * Q * Q * Q;        //    对应论文公式8c
    B1 = 2.44413 * Q + 2.85619 * Q * Q + 1.26661 * Q * Q * Q;
    B2 = -1.4281 * Q * Q - 1.26661 * Q * Q * Q;
    B3 = 0.422205 * Q * Q * Q;

    B0 = 1.0 - (B1 + B2 + B3) / B;
    B1 = B1 / B;
    B2 = B2 / B;
    B3 = B3 / B;
}

图片 5View Code

  那里也是3回性处理八个像素,那里自身利用了其余一种转移技术来把陆位转换为十八个人,不过后面包车型地铁Diff因为有正有负,要转换为三十四人就亟须采纳_mm_cvtepi16_epi32来更换,不可能用unpack连串组合函数来促成,因为unpack不会活动符号位,作者在那边吃了一些次亏了。

  由上述代码可知,B0+B1+B2+B3=1,就是归一化的周详,这也是算法能够递归进行的要紧之一。

  代码相比多,作者稍稍精简下(注意,精简后的是不行运营的,只是更鲜明的发挥了思路)。

  接着是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优化的数据
}
void ConvertBGR8U2BGRAF(unsigned char *Src, float *Dest, int Width, int Height, int Stride)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        float *LinePD = Dest + Y * Width * 3;
        for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3)
        {
            LinePD[0] = LinePS[0];    LinePD[1] = LinePS[1];    LinePD[2] = LinePS[2];
        }
    }
}

     
在行方向上,ColSum种种成分的换代互相之间是没有其余关联的,因此,借助于SSE能够兑现三次性的五个成分的翻新,并且上述代码还将首先行的独特总括也用SSE达成了,即便这一部分总结量是老大小的。

  镜像正是倒序的进程,直接使用SSE的shuffle函数很有益达成。

  我们可以品尝本人把里面的X循环再开始展览探访效果。

   
 在具体的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;
}
void GaussBlurFromLeftToRight(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePD = Data + Y * Width * 3;
        float BS1 = LinePD[0], BS2 = LinePD[0], BS3 = LinePD[0];          //  边缘处使用重复像素的方案
        float GS1 = LinePD[1], GS2 = LinePD[1], GS3 = LinePD[1];
        float RS1 = LinePD[2], RS2 = LinePD[2], RS3 = LinePD[2];
        for (int X = 0; X < Width; X++, LinePD += 3)
        {
            LinePD[0] = LinePD[0] * B0 + BS1 * B1 + BS2 * B2 + BS3 * B3;
            LinePD[1] = LinePD[1] * B0 + GS1 * B1 + GS2 * B2 + GS3 * B3;         // 进行顺向迭代
            LinePD[2] = LinePD[2] * B0 + RS1 * B1 + RS2 * B2 + RS3 * B3;
            BS3 = BS2, BS2 = BS1, BS1 = LinePD[0];
            GS3 = GS2, GS2 = GS1, GS1 = LinePD[1];
            RS3 = RS2, RS2 = RS1, RS1 = LinePD[2];
        }
    }
}

     
注意到在种种SSE代码后边,总还有部分C代码,那是因为我们实际上数目大幅并不一定是4的整数倍,未被SSE处理的一部分必须运用C达成,那实在对读者来说是13分好的事务,因为大家能从这一部分C代码中搞精通上述的SSE代码是干什么的。

  使用多少个__m128i
变量首假若为了充裕利用XMM寄存器,在那之中SumI32函数如下,首即使为着计算__m128i
内八个整数的累加值。

  不多描述,请我们把上述代码和公式(a)对应一下就清楚了。

  以上正是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);
}

     
有了GaussBlurFromLeftToRight的参考代码,GaussBlurFromRightToLeft的代码就不会有怎么样大的不方便了,为了幸免有个别懒人不看小说不考虑,那里作者不贴那段的代码。

  那是还是不是还有创新的长空吧,从作者的认知中,在对垂直方向的拍卖上,应该没有了,不过程度方向呢,
SSE有没有用武之地,经过本身的盘算自个儿以为依然有的。

  代码不解释。

     
那么垂直方向上简单的做只要求改变下循环的主旋律,以及历次指针扩张量更改为Width
* 3即可。

  在行方向的计量中,那几个for循环是关键的乘除。

  最后便是100到104行的代码的转换。

     
那么我们来考虑下垂直方向能否不那样处理啊,指针每一趟扩展Width *
3会促成严重的Cache
miss,特别是对此宽度稍微大点的图像,大家精心考察垂直方向,会发觉依旧能够依照Y
 — X那样的大循环格局也是能够的,代码如下:

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);
}
void GaussBlurFromTopToBottom(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePD3 = Data + (Y + 0) * Width * 3;
        float *LinePD2 = Data + (Y + 1) * Width * 3;
        float *LinePD1 = Data + (Y + 2) * Width * 3;
        float *LinePD0 = Data + (Y + 3) * Width * 3;
        for (int X = 0; X < Width; X++, LinePD0 += 3, LinePD1 += 3, LinePD2 += 3, LinePD3 += 3)
        {
            LinePD0[0] = LinePD0[0] * B0 + LinePD1[0] * B1 + LinePD2[0] * B2 + LinePD3[0] * B3;
            LinePD0[1] = LinePD0[1] * B0 + LinePD1[1] * B1 + LinePD2[1] * B2 + LinePD3[1] * B3;
            LinePD0[2] = LinePD0[2] * B0 + LinePD1[2] * B1 + LinePD2[2] * B2 + LinePD3[2] * B3;
        }
    }
}

     
 自己觉得固然那里的操作是内外信赖的,全局无法并行化,不过观望这一行代码:

  在此之前觉得那么些算法难以使用SSE优化,重要困难就在此地,可是在就学了Opencv的积分图技术时,这些标题就一举成功了,进一步分析还发现Opencv的代码写的并不圆满,还有创新的空中,见上述代码,使用两回对同样数据移位就能够得到充分,由P3
P2 P1 P0变为P3+P2+P1+P0 P2+P1+P0 P1+P0 P0。大家稍事分析一下。          
  

  正是说大家不是一下子就把列方向的前向传来举办完,而是每便只举办一个像素的传遍,当一行全体像素都实行完了列方向的前向传来后,在切换来下一行,这样处理就幸免了Cache
miss出现的情景。

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就应当是:

   
 注意到列方向的边缘地方,为了便利代码的拍卖,大家在高度方向上上下分别扩张了1个行的像素,当举行完全中学间的得力行的行方向前向和反向传来后,依据前述的重新边缘像素的艺术填充上下那空出的三行数据。

      其中RowData[X + Size – 1] –
RowData[X – 1];
并不是左右信赖的,是能够相互的,因而只要提前急速的估算出装有的差值,那么提速的空中还比较大,即要求超前计算出上边包车型客车函数:

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

   
 同样的道理,GaussBlurFromBottomToTop的代码可由读者本身补充进去。

 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)); 这一句为什么要这样写,恐怕还是读者自己体会比较好,似乎不好用文字表达。

   
 最终的ConvertBGRAF2BG帕杰罗8U也很简单,就是3个for循环:

  那里用SSE达成则卓殊不难了:

   最终贰个_mm_storesi128_4char是本人自个儿定义的2个将二个__m128i里面包车型客车多少个叁10个人整数转换为字节数并保留到内存中的函数,详见附属类小部件代码。

void ConvertBGRAF2BGR8U(float *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePS = Src + Y * Width * 3;
        unsigned char *LinePD = Dest + Y * Stride;
        for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3)
        {
            LinePD[0] = LinePS[0];    LinePD[1] = LinePS[1];    LinePD[2] = LinePS[2];
        }
    }
}
        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];

  至于2三个人图像的优化,是个相比狼狈的景况,因为SSE3遍性处理6个三11位数,而23位各类像素只有二个轻重,那种气象相似依旧把她恢弘到七个轻重,比如说ColValue就足以改成4大路的,然后累积和也供给处理成4通道的,那自然须求一定的迷你淫技,那里本人不想多谈,留点东西给自身。当然也得以设想先把2四个人分解到一个灰度内部存款和储蓄器,然后采取灰度的算法进行测算,前面在合成到2几人,这么些解释和合成的历程也是足以采用SSE加快的,假设你结合八线程,还能把三个灰度进程的处理并行化,那样除了多占用内部存款和储蓄器外,速度比至二级处理22位要快(因为平昔处清理计算法不可能并行的,前后注重的来头)。此外正是最终在盘算列累积求平均值的长河也变得进一步自然了,不会产出灰度那种要在__mm128i里边开始展览添加的历程,而是径直得三个SSE变量累加。

   在使得的限定内,上述浮点总计的结果不会超出byte所能表达的范围,因而也不要求特其余拓展Clamp操作。

  和列方向的SSE处理代码分裂的是,那里加载的是uchar数据,因而加载的函数就截然分歧,处理的方法也有分别,上边多少个SSE函数各位查查MSDN就能了然其意思,如故很有意味的。

  还说一点,未来多数的CPU都帮忙AVX256了,还足以行使AVX进一步加速,就像代码该起来也不是很难,有趣味的对象可以自身尝试。

     
 最终正是局地内部存储器分配和自由的代码了,如下所示:

  经过这么的处理,经过测试发现,速度能够比opencv的版本在拉长百分之三十,也是卓殊的大悲大喜。

  能够说,近期以此速度已经大致达到了CPU的终端了,不过测试过IPP的速度,就像比那些还要快点,不消除他动用了AVX,也不免除他采纳多核的财富。

void GaussBlur(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius)
{
    float B0, B1, B2, B3;
    float *Buffer = (float *)malloc(Width * (Height + 6) * sizeof(float) * 3);

    CalcGaussCof(Radius, B0, B1, B2, B3);
    ConvertBGR8U2BGRAF(Src, Buffer + 3 * Width * 3, Width, Height, Stride);

    GaussBlurFromLeftToRight(Buffer + 3 * Width * 3, Width, Height, B0, B1, B2, B3);
    GaussBlurFromRightToLeft(Buffer + 3 * Width * 3, Width, Height, B0, B1, B2, B3);        //    如果启用多线程,建议把这个函数写到GaussBlurFromLeftToRight的for X循环里,因为这样就可以减少线程并发时的阻力

    memcpy(Buffer + 0 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float));
    memcpy(Buffer + 1 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float));
    memcpy(Buffer + 2 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float));

    GaussBlurFromTopToBottom(Buffer, Width, Height, B0, B1, B2, B3);

    memcpy(Buffer + (Height + 3) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float));
    memcpy(Buffer + (Height + 4) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float));
    memcpy(Buffer + (Height + 5) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float));

    GaussBlurFromBottomToTop(Buffer, Width, Height, B0, B1, B2, B3);

    ConvertBGRAF2BGR8U(Buffer + 3 * Width * 3, Dest, Width, Height, Stride);

    free(Buffer);
}

     
再有的优化只怕提速有限,比如把剩余的一部分for循环内部分成四路等等。

  那些的优化对于BoxBlur来说是第贰的一步,不过更首要的是她能够行使在不少场所,比如图像的一部分均方差计算,也能够动用类似的技艺拓展加快,两幅图像大的一些平方差也是能够如此优化的,后续作者会不难的谈下他在那上边加速的利用。

  正如上所述,分配了Height +
6行的内部存款和储蓄器区域,主借使为着便利垂直方向的处理,在实践GaussBlurFromTopToBottom从前根据重复边缘的尺度复制3行,然后在GaussBlurFromBottomToTop从前在复制底部边缘的3行像素。

   
 在笔者的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

     
至此,3个简易而又飞速的高斯模糊就大旨到位了,代码数量也不多,也未曾啥难度,不明白为何许四个人搞不定。

     
在不利用八线程(就算本算法非凡适合二十三四线程总计),不选取GPU,只使用单线程\CPU举办计算的事态下,私家认为眼下自身这一个代码是很难有质的超过的,从有些地方讲,代码中的用于总括的岁月占据的小时比从内部存款和储蓄器等待数据的小时大概还要短,而就像也从未更好的算法能进一步压缩内部存款和储蓄器数据的访问量。

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

     
依照本人的测试,上述措施代码在一台I5-6300HQ
2.30GHZ的记录簿上针对一副三千*三千的2二个人图像的处理时间大体必要370ms,假如在C++的编写翻译选项的代码生成选项里的启用增强指令集选取–>流式处理SIMD扩张2(/arch:sse2),则编译后的程序大约必要220ms的小运。

     
自身在此间约请各位高手对当前自个儿优化的这么些代码实行更进一步的优化,希望高人不要客气。

图片 6

     
大家尝试选择系统的片段财富进一步升高速度,首先大家想到了SSE优化,关于那上边英特尔有一篇官方的篇章和代码:

  运营界面:

 

        IIR Gaussian Blur Filter
Implementation using Intel® Advanced Vector
Extensions
 [PDF
513KB]

图片 7

  糟糕意思,图太小,速度为0ms了。

     source
code: gaussian_blur.cpp [36KB]

  

图片 8

     
笔者只是简短的看了下那篇文章,感觉他里面用到的计算公式和Deriche滤波器的很像,和本文参考的Recursive
implementation
不太雷同,并且其SSE代码对能处理的图还有不少限量,对本身这里的参考意义相当的小。

  本文的代码是本着常用的图像数据开始展览的优化处理,在广大地方下,需求对int恐怕float类型实行拍卖,比如GuideFilter,如若读者知道了本文解读的代码的法则,更改代码以便适应他们则是一件万分简单的政工,借使您不会,那么也请您不用给本身留言,只怕本身能够有偿提供,因为从精神上讲小编喜爱提供渔,而不是鱼,渔具已经送给你了,你却找作者要鱼,那只好…………..。

  

     
大家先看下要旨的计量的SSE优化,注意到  GaussBlurFromLeftToRight
的代码中,
其基本的一个钱打二16个结部分是多少个乘法,然而她唯有三个乘法总计,假若能够凑成四行,那么就能够充裕利用SSE的批量计量效率了,约等于假若能扩大多个大路,使得GaussBlurFromLeftToRight变为如下方式:

   
  本文源代码下载地址(VS二〇一〇编纂):Box
Filter
 

  

void GaussBlurFromLeftToRight(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePD = Data + Y * Width * 4;
        float BS1 = LinePD[0], BS2 = LinePD[0], BS3 = LinePD[0];                //  边缘处使用重复像素的方案
        float GS1 = LinePD[1], GS2 = LinePD[1], GS3 = LinePD[1];
        float RS1 = LinePD[2], RS2 = LinePD[2], RS3 = LinePD[2];
        float AS1 = LinePD[3], AS2 = LinePD[3], AS3 = LinePD[3];
        for (int X = 0; X < Width; X++, LinePD += 4)
        {
            LinePD[0] = LinePD[0] * B0 + BS1 * B1 + BS2 * B2 + BS3 * B3;
            LinePD[1] = LinePD[1] * B0 + GS1 * B1 + GS2 * B2 + GS3 * B3;         // 进行顺向迭代
            LinePD[2] = LinePD[2] * B0 + RS1 * B1 + RS2 * B2 + RS3 * B3;
            LinePD[3] = LinePD[3] * B0 + AS1 * B1 + AS2 * B2 + AS3 * B3;
            BS3 = BS2, BS2 = BS1, BS1 = LinePD[0];
            GS3 = GS2, GS2 = GS1, GS1 = LinePD[1];
            RS3 = RS2, RS2 = RS1, RS1 = LinePD[2];
            AS3 = AS2, AS2 = AS1, AS1 = LinePD[3];
        }
    }
}

图片 9

 

  则很不难就把上述代码转换到SSE能够规范处理的代码了。

 

 

  而对于Y方向的代码,你细心观望会发现,无论是及通道的图,天然的就能够动用SSE进行处理,详见上面包车型大巴代码。

 

 

  好,大家还是1个三个的来分析:

 ****************************作者:
laviewpbt   时间: 二零一四.12.17    联系QQ:  33184777
转载请保留本行音讯**********************

 

   第一个函数
CalcGaussCof 无需实行任何的优化。

   

 

     
首个函数 ConvertBG奥迪Q58U2BGRAF根据上述分析需求重新写,因为急需追加一个通路,新的大道的值填0可能其余值都能够,但提议填0,那对某个SSE函数很有用,小编把这几个函数的SSE实现共享一下:

void ConvertBGR8U2BGRAF_SSE(unsigned char *Src, float *Dest, int Width, int Height, int Stride)
{
    const int BlockSize = 4;
    int Block = (Width - 2) / BlockSize;
    __m128i Mask = _mm_setr_epi8(0, 1, 2, -1, 3, 4, 5, -1, 6, 7, 8, -1, 9, 10, 11, -1);            //    Mask为-1的地方会自动设置数据为0
    __m128i Zero = _mm_setzero_si128();
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        float *LinePD = Dest + Y * Width * 4;
        int X = 0;
        for (; X < Block * BlockSize; X += BlockSize, LinePS += BlockSize * 3, LinePD += BlockSize * 4)
        {
            __m128i SrcV = _mm_shuffle_epi8(_mm_loadu_si128((const __m128i*)LinePS), Mask);        //    提取了16个字节,但是因为24位的,我们要将其变为32位的,所以只能提取出其中的前12个字节
            __m128i Src16L = _mm_unpacklo_epi8(SrcV, Zero);
            __m128i Src16H = _mm_unpackhi_epi8(SrcV, Zero);
            _mm_store_ps(LinePD + 0, _mm_cvtepi32_ps(_mm_unpacklo_epi16(Src16L, Zero)));        //    分配内存时已经是16字节对齐了,然后每行又是4的倍数的浮点数,因此,每行都是16字节对齐的
            _mm_store_ps(LinePD + 4, _mm_cvtepi32_ps(_mm_unpackhi_epi16(Src16L, Zero)));        //    _mm_stream_ps是否快点?
            _mm_store_ps(LinePD + 8, _mm_cvtepi32_ps(_mm_unpacklo_epi16(Src16H, Zero)));
            _mm_store_ps(LinePD + 12, _mm_cvtepi32_ps(_mm_unpackhi_epi16(Src16H, Zero)));
        }
        for (; X < Width; X++, LinePS += 3, LinePD += 4)
        {
            LinePD[0] = LinePS[0];    LinePD[1] = LinePS[1];    LinePD[2] = LinePS[2];    LinePD[3] = 0;        //    Alpha最好设置为0,虽然在下面的CofB0等SSE常量中通过设置ALPHA对应的系数为0,可以使得计算结果也为0,但是不是最合理的
        }
    }
}

  稍作解释,_mm_loadu_si128二遍性加载拾陆个字节的数量到SSE寄存器中,对于二十四位图像,17个字节里带有了5

  • 1 /
    二个像素的新闻,而大家的靶子是把那个数据转换为4通道的音讯,由此,大家只可以三遍性的领取处16/4=多少个像素的值,那借助于_mm_shuffle_epi8函数和正好的Mask来达成,_mm_unpacklo_epi8/_mm_unpackhi_epi八分别领取了SrcV的高八位和低陆位的八个字节数据并将它们转换为七个16进制数保存到Src16L和Src16H中,
    而_mm_unpacklo_epi16/_mm_unpackhi_epi16则越是把十十一位数据扩张到三11人整形数据,最后经过_mm_cvtepi32_ps函数把整形数据转换为浮点型。

  或者有人注意到了 int Block = (Width – 2) /
BlockSize;
这一行代码,为啥要-2操作呢,那也是本人在反复测试发现先后连接出现内部存储器错误时才发现到的,因为_mm_loadu_si1283回性加载了5

  • 1 / 二个像素的新闻,当在处理最终一行像素时(其余行不会),要是Block
    取Width/BlockSize,
    则很有恐怕拜会了大于像素范围内的内存,而-2不是-1就是因为分外额外的33.33%像素起的效果。

  接着看档次方向的前向传来的SSE方案:

void GaussBlurFromLeftToRight_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    const  __m128 CofB0 = _mm_set_ps(0, B0, B0, B0);
    const  __m128 CofB1 = _mm_set_ps(0, B1, B1, B1);
    const  __m128 CofB2 = _mm_set_ps(0, B2, B2, B2);
    const  __m128 CofB3 = _mm_set_ps(0, B3, B3, B3);
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePD = Data + Y * Width * 4;
        __m128 V1 = _mm_set_ps(LinePD[3], LinePD[2], LinePD[1], LinePD[0]);
        __m128 V2 = V1, V3 = V1;
        for (int X = 0; X < Width; X++, LinePD += 4)            //    还有一种写法不需要这种V3/V2/V1递归赋值,一次性计算3个值,详见 D:\程序设计\正在研究\Core\ShockFilter\Convert里的高斯模糊,但速度上没啥区别
        {
            __m128 V0 = _mm_load_ps(LinePD);
            __m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1));
            __m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3));
            __m128 V = _mm_add_ps(V01, V23);
            V3 = V2;    V2 = V1;    V1 = V;
            _mm_store_ps(LinePD, V);
        }
    }
}

  和上边的4大路的GaussBlurFromLeftToRight_SSE相比较,你会意识差不离语法上从未有过其余的区分,实在是太不难了,注意本身并未用_mm_storeu_ps,而是间接运用_mm_store_ps,那是因为,第②,分配Data内部存款和储蓄器时,小编使用了_mm_malloc分配了16字节对齐的内部存款和储蓄器,而Data每行成分的个数又都以4的倍数,由此,每行早先地方处的内部存款和储蓄器也是16字节对齐的,所以直接_mm_store_ps完全是可以的。

   
 同理,在笔直方向的前向传来的SSE优化代码就更直白了:

void GaussBlurFromTopToBottom_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    const  __m128 CofB0 = _mm_set_ps(0, B0, B0, B0);
    const  __m128 CofB1 = _mm_set_ps(0, B1, B1, B1);
    const  __m128 CofB2 = _mm_set_ps(0, B2, B2, B2);
    const  __m128 CofB3 = _mm_set_ps(0, B3, B3, B3);
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePS3 = Data + (Y + 0) * Width * 4;
        float *LinePS2 = Data + (Y + 1) * Width * 4;
        float *LinePS1 = Data + (Y + 2) * Width * 4;
        float *LinePS0 = Data + (Y + 3) * Width * 4;
        for (int X = 0; X < Width * 4; X += 4)
        {
            __m128 V3 = _mm_load_ps(LinePS3 + X);
            __m128 V2 = _mm_load_ps(LinePS2 + X);
            __m128 V1 = _mm_load_ps(LinePS1 + X);
            __m128 V0 = _mm_load_ps(LinePS0 + X);
            __m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1));
            __m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3));
            _mm_store_ps(LinePS0 + X, _mm_add_ps(V01, V23));
        }
    }
}

  对地方的代码不想做任何解释了。

  最有难度的应有算是ConvertBGRAF2BGEvoque8U的SSE版本了,由于一些原因,笔者不太情愿分享这些函数的代码,有趣味的仇人能够参考opencv的关于落实。

   
 最终的SSE版本高斯模糊的最首要代码如下:

void GaussBlur_SSE(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius)
{
    float B0, B1, B2, B3;
    float *Buffer = (float *)_mm_malloc(Width * (Height + 6) * sizeof(float) * 4, 16);

    CalcGaussCof(Radius, B0, B1, B2, B3);
    ConvertBGR8U2BGRAF_SSE(Src, Buffer + 3 * Width * 4, Width, Height, Stride);

    GaussBlurFromLeftToRight_SSE(Buffer + 3 * Width * 4, Width, Height, B0, B1, B2, B3);        //    在SSE版本中,这两个函数占用的时间比下面两个要多,不过C语言版本也是一样的
    GaussBlurFromRightToLeft_SSE(Buffer + 3 * Width * 4, Width, Height, B0, B1, B2, B3);        //    如果启用多线程,建议把这个函数写到GaussBlurFromLeftToRight的for X循环里,因为这样就可以减少线程并发时的阻力

    memcpy(Buffer + 0 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));
    memcpy(Buffer + 1 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));
    memcpy(Buffer + 2 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));

    GaussBlurFromTopToBottom_SSE(Buffer, Width, Height, B0, B1, B2, B3);

    memcpy(Buffer + (Height + 3) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float));
    memcpy(Buffer + (Height + 4) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float));
    memcpy(Buffer + (Height + 5) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float));

    GaussBlurFromBottomToTop_SSE(Buffer, Width, Height, B0, B1, B2, B3);

    ConvertBGRAF2BGR8U_SSE(Buffer + 3 * Width * 4, Dest, Width, Height, Stride);

    _mm_free(Buffer);
}

  对于同样的三千*两千的彩色图像,SSE版本的代码耗费时间只有145ms的耗费时间,相对于一般性的C代码有约2.5倍左右的提速,这也未可厚非,因为我们在执行SSE的代码时时多处理了三个大路的总括量的,可是同编写翻译器自个儿的SSE优化220ms,只有1.5倍的提速了,那注明编写翻译器的SSE优化能力可能万分强的。

   
 进一步的优化就是本人上面的注释掉的opemp相关代码,在ConvertBG奥迪Q58U2BGRAF /
GaussBlurFromLeftToRight / GaussBlurFromRightToLeft / ConvertBGRAF2BGCRUISER8U
 伍个函数中,直接助长不难的#pragma omp parallel
for就足以了,至于GaussBlurFromTopToBottom_SSE、
GaussBlurFromBottomToTop_SSE则由于上下行之间数据的相关性,是心有余而力不足完成并行总结的,然而测试表示他们的耗费时间要比水平方向的少很多。

   
比如我们钦赐openmp使用三个线程,在上述机器上(四核的),纯C版本能优化到252ms,而纯SSE的只好增强到100ms左右,编写翻译器本身的SSE优化速度大致是150ms,基本上依旧维持同叁个级其他百分比。

   对于灰度图像,很惋惜,上述的程度方向上的优化措施就随便为力了,一种办法就是把4行灰度像素混洗成一行,不过那一个操作不太方便用SSE达成,其它一种就是把水平方向的数据先转置,然后采取垂直方向的SSE优化代码处理,实现在转置回去,最终对转置的数目再一次开始展览垂直方向SSE优化,当然转置的长河是足以凭借SSE的代码实现的,然则须要分外的一份内部存储器,速度上恐怕和一般的C相比较就不会有那么多的升官了,那一个待有时光了再去测试。

   
 前左右后写那么些大致也花了半个月的年华,写完了总体流程,在倒过来看,真的是万分的简约,有的时候正是那样。

   
 本文并没有提供完整的可以直接执行的漫天代码,需者自勤,提供一段C#的调用工程供有趣味的情人测试和比对(未使用OPENMP版本)。

   
 http://files.cnblogs.com/files/Imageshop/GaussBLur_Sample.rar

图片 10

     后记:后来测试发现,当半径参数较大时,无论是C版本仍然SSE版本都晤面世界形势部有失水准的败笔,感觉像是溢出了,后来意识重庆大学是当半径变大后,B0参数变得不大,以至于用float类型的数目来处理递归已经黔驴技穷保险丰裕的精度了,消除的艺术是利用double类型,那对于C语言版本的话是秒秒的事务,而对于SSE版本则是较大的横祸,double时换到AVX大概改动量十分的小,然而AVX的普及度以及…..,可是,一般景象下,半径不超越75时结果都依旧不错的,那对于绝当先八分之四的使用来说是十足的,半径75时,整个图像已经大概没有其他的底细了,再大,差别也相当的小了。

   
 消除上述难点3个得力的方案就是运用Deriche滤波器,我用该滤波器的float版本对大半径是不汇合世难点的,并且也有连锁SSE参考代码。

 图片 11

   后续小说:高斯模糊算法的健全优化进程分享(二)。

 

相关文章