SSE图像算法优化系列五:超高速指数模糊算法的实现和优化(10000*10000在100ms左右兑现)。

     
今日我们来花点时间重新谈谈一个模糊算法,一个至上简单不过又顶尖牛逼的算法,无论在职能上或者速度上都得以和Boxblur,
stackblur或者是Gaussblur想媲美,效果上,比Boxblur来的更平整,和Gaussblur相似,速度上,经过自己的优化,在PC端比她们多少个都要快一大截,而且着力不需占用额外的内存,实在是一个绝好的算法。

     
算法的着力并不是自己想开仍然发明的,是一个情侣在github上挖掘到的,率属于Cairo这些2D图形库的开源代码,详见:

     
 https://github.com/rubencarneiro/NUX/blob/master/NuxGraphics/CairoGraphics.cpp

       大家称之为Exponential
blur(指数模糊)。

     
 提炼出这一个模糊的基本部分算法紧若是下面那些函数:

   static inline void _blurinner(guchar* pixel, gint* zR,gint* zG,gint* zB, gint* zA, gint alpha,gint aprec,gint zprec)
  {
    gint R, G, B, A;
    R = *pixel;
    G = *(pixel + 1);
    B = *(pixel + 2);
    A = *(pixel + 3);

    *zR += (alpha * ((R << zprec) - *zR)) >> aprec;
    *zG += (alpha * ((G << zprec) - *zG)) >> aprec;
    *zB += (alpha * ((B << zprec) - *zB)) >> aprec;
    *zA += (alpha * ((A << zprec) - *zA)) >> aprec;

    *pixel       = *zR >> zprec;
    *(pixel + 1) = *zG >> zprec;
    *(pixel + 2) = *zB >> zprec;
    *(pixel + 3) = *zA >> zprec;
  }

  
 其中Pixel就是咱们要处理的像素,zR,zG,zB,zA是外部传入的一个累加值,alpha、aprec、zprec是由模糊半径radius生成的一对稳定的周详。

  类似于高斯模糊或者StackBlur,算法也是属于行列分离的,行方向上,先用上述办法从左向右统计,然后在从右向左,接着进行列方向处理,先从上向下,然后在从下向上。当然,行列的乘除也足以扭转。需要注意的是,每一步都是用事先处理过的数量举办的。

  在源代码中用以下五个函数实现以下过程:

     水平反向的处理:

  static inline void _blurrow(guchar* pixels,
                               gint    width,
                               gint    /* height */,  // TODO: This seems very strange. Why is height not used as it is in _blurcol() ?
                               gint    channels,
                               gint    line,
                               gint    alpha,
                               gint    aprec,
                               gint    zprec)
  {
    gint    zR;
    gint    zG;
    gint    zB;
    gint    zA;
    gint    index;
    guchar* scanline;

    scanline = &(pixels[line * width * channels]);

    zR = *scanline << zprec;
    zG = *(scanline + 1) << zprec;
    zB = *(scanline + 2) << zprec;
    zA = *(scanline + 3) << zprec;

    for (index = 0; index < width; index ++)
      _blurinner(&scanline[index * channels], &zR, &zG, &zB, &zA, alpha, aprec,
      zprec);

    for (index = width - 2; index >= 0; index--)
      _blurinner(&scanline[index * channels], &zR, &zG, &zB, &zA, alpha, aprec,
      zprec);
  }

  垂直方向的拍卖:

static inline void _blurcol(guchar* pixels,
                               gint    width,
                               gint    height,
                               gint    channels,
                               gint    x,
                               gint    alpha,
                               gint    aprec,
                               gint    zprec)
  {
    gint zR;
    gint zG;
    gint zB;
    gint zA;
    gint index;
    guchar* ptr;

    ptr = pixels;

    ptr += x * channels;

    zR = *((guchar*) ptr    ) << zprec;
    zG = *((guchar*) ptr + 1) << zprec;
    zB = *((guchar*) ptr + 2) << zprec;
    zA = *((guchar*) ptr + 3) << zprec;

    for (index = width; index < (height - 1) * width; index += width)
      _blurinner((guchar*) &ptr[index * channels], &zR, &zG, &zB, &zA, alpha,
      aprec, zprec);

    for (index = (height - 2) * width; index >= 0; index -= width)
      _blurinner((guchar*) &ptr[index * channels], &zR, &zG, &zB, &zA, alpha,
      aprec, zprec);
  }

  最后的混淆算法如下所示:

  // pixels   image-data
  // width    image-width
  // height   image-height
  // channels image-channels
  // in-place blur of image 'img' with kernel of approximate radius 'radius'
  // blurs with two sided exponential impulse response
  // aprec = precision of alpha parameter in fixed-point format 0.aprec
  // zprec = precision of state parameters zR,zG,zB and zA in fp format 8.zprecb

  void _expblur(guchar* pixels,
                 gint    width,
                 gint    height,
                 gint    channels,
                 gint    radius,
                 gint    aprec,
                 gint    zprec)
  {
    gint alpha;
    gint row = 0;
    gint col = 0;

    if (radius < 1)
      return;

    // calculate the alpha such that 90% of 
    // the kernel is within the radius.
    // (Kernel extends to infinity)
    alpha = (gint) ((1 << aprec) * (1.0f - expf(-2.3f / (radius + 1.f))));

    for (; row < height; row++)
      _blurrow(pixels, width, height, channels, row, alpha, aprec, zprec);

    for (; col < width; col++)
      _blurcol(pixels, width, height, channels, col, alpha, aprec, zprec);

    return;
  }

  作为一个突出的施用,或者说尽量缩小参数,常用的aprec取值为16,Zprec
取值为7。

   
 回顾下代码,全体过程中除了alpha参数的揣摸涉及到了浮点,其他一些都是整形的乘法和移动操作,因而可以设想,速度相应不慢,而且十分适合于手机端处理。同时注意到_blurrow和_blurcol函数循环分明互相之间是单独的,能够利用多线程并行处理,可是那么些代码重假诺专注于算法的公布,并从未过多的设想更好的频率。

   
 此外一些,很明确,算法的耗时是和Radius参数没有其它关联的,也就是说这也是个O(1)算法。

  我们有点对上述代码做个简化处理,对于灰度图,水平方向的代码可以表达如下:

for (int Y = 0; Y < Height; Y++)
{
    byte *LinePS = Src + Y * Stride;
    byte *LinePD = Dest + Y * Stride;
    int Sum = LinePS[0] << Zprec;
    for (int X = 0; X < Width; X++)      //  从左往右
    {
        Sum += (Alpha * ((LinePS[X] << Zprec) - Sum)) >> Aprec;
        LinePD[X] = Sum >> Zprec;
    }
    for (int X = Width - 1; X >= 0; X--)   //  从右到左
    {
        Sum += (Alpha * ((LinePD[X] << Zprec) - Sum)) >> Aprec;
        LinePD[X] = Sum >> Zprec;
    }
}

  在 高斯模糊算法的通盘优化过程分享(一) 中我们探索过垂直方向处理算法一般不宜直接写,而应当用一个暂时的行缓存举行处理,这样列方向的灰度图的拍卖代码类似上边的:

int *Buffer = (int *)malloc(Width * sizeof(int));
for (int X = 0; X < Width; X++)        Buffer[X] = Src[X] << Zprec;
for (int Y = 0; Y < Height; Y++)
{
    byte *LinePS = Src + Y * Stride;
    byte *LinePD = Dest + Y * Stride;
    for (int X = 0; X < Width; X++)        //  从上到下
    {
        Buffer[X] += (Alpha * ((LinePS[X] << Zprec) - Buffer[X])) >> Aprec;
        LinePD[X] = Buffer[X] >> Zprec;
    }
}
for (int Y = Height - 1; Y >= 0; Y--)      //  从下到上
{
    byte *LinePD = Dest + Y * Stride;
    for (int X = 0; X < Width; X++)
    {
        Buffer[X] += (Alpha * ((LinePD[X] << Zprec) - Buffer[X])) >> Aprec;
        LinePD[X] = Buffer[X] >> Zprec;
    }
}
free(Buffer);

   修改为上述后,测试一个3000*2000的8位灰度图,耗时大约52ms(未使用多线程的),和平时的C语言实现的Boxblur时间基本上。

     
 除了线程外,这些时间是不是还有立异的上空啊,大家先来看望列方向的优化。

   在列方向的  for (int X = 0; X <
Width; X++)
循环内,我们注意到针对Buffer的各种元素的拍卖都是独立和同一的,很醒目这样的经过是很容易拔取SIMD指令优化的,可是循环体内部有一些是unsigned
char类型的多少,为运用SIMD指令,需要更换为int类型较为有利,而最终保存时又需要重新处理为unsigned
char类型的,这种往返转换的耗时和其它计量的涨潮能否来拉动意义呢,我们开展了代码的编制,比如:

    for (int X = 0; X < Width; X++)        //  从上到下
    {
        Buffer[X] += (Alpha * ((LinePS[X] << Zprec) - Buffer[X])) >> Aprec;
        LinePD[X] = Buffer[X] >> Zprec;
    }

  这段代码可以用如下的SIMD指令代替:

int X = 0;
for (X = 0; X < Width - 8; X += 8)            
{
    //    将8个字节数存入到2个XMM寄存器中
    //    方案1:使用SSE4新增的_mm_cvtepu8_epi32的函数,优点是两行是独立的
    __m128i Dst1 = _mm_cvtepu8_epi32(_mm_cvtsi32_si128((*(int *)(LinePD + X + 0))));    //    _mm_cvtsi32_si128把参数中的32位整形数据移到XMM寄存器的最低32位,其他为清0。
    __m128i Dst2 = _mm_cvtepu8_epi32(_mm_cvtsi32_si128((*(int *)(LinePD + X + 4))));    //    _mm_cvtepu8_epi32将低32位的整形数的4个字节直接扩展到XMM的4个32位中去。
    __m128i Buf1 = _mm_loadu_si128((__m128i *)(Buffer + X + 0));
    __m128i Buf2 = _mm_loadu_si128((__m128i *)(Buffer + X + 4));
    Buf1 = _mm_add_epi32(_mm_srai_epi32(_mm_mullo_epi32(_mm_sub_epi32(_mm_slli_epi32(Dst1, Zprec), Buf1), Alpha128), Aprec), Buf1);
    Buf2 = _mm_add_epi32(_mm_srai_epi32(_mm_mullo_epi32(_mm_sub_epi32(_mm_slli_epi32(Dst2, Zprec), Buf2), Alpha128), Aprec), Buf2);
    _mm_storeu_si128((__m128i *)(Buffer + X + 0), Buf1);
    _mm_storeu_si128((__m128i *)(Buffer + X + 4), Buf2);
    _mm_storel_epi64((__m128i *)(LinePD + X), _mm_packus_epi16(_mm_packs_epi32(_mm_srai_epi32(Buf1, Zprec), _mm_srai_epi32(Buf2, Zprec)), Zero));
}
for (; X < Width; X++)        
{
    Buffer[X] += (Alpha * ((LinePD[X] << Zprec) - Buffer[X])) >> Aprec;
    LinePD[X] = Buffer[X] >> Zprec;
}

  原来的三四行代码一下子改为了几十行的代码,会不会变慢呢,其实不用担心,SIMD真的很强劲,测试的结果是3000*2000的图耗时下降到42ms左右,而且垂直方向的耗时占比有原来的60%下降到了35%左右,现在的主题就是水平方向的耗时了。

   
 当图像不是灰度格局时,对于垂直方向的处理和灰度不会有分别,这是因为,只需要扩张循环的尺寸就足以了。

   
 我们再来看看水平方向的优化,当图像是ARGB情势时,也就是原作者的代码,总结过程每隔多少个字节就会再也,这种特点当然也合乎SIMD指令,但是为了方便,必须得先将字节数据先转移为int类型的一个缓冲区中,之后从左到右的总括可以用如下的代码实现:

void ExpFromLeftToRight_OneLine_SSE(int *Data, int Length, int Radius, int Aprec, int Zprec, int Alpha)
{
    int *LinePD = Data;
    __m128i A = _mm_set1_epi32(Alpha);
    __m128i S1 = _mm_slli_epi32(_mm_load_si128((__m128i *)(LinePD)), Zprec);
    for (int X = 0; X < Length; X++, LinePD += 4)
    {
        S1 = _mm_add_epi32(S1, _mm_srai_epi32(_mm_mullo_epi32(_mm_sub_epi32(_mm_slli_epi32(_mm_load_si128((__m128i *)(LinePD)), Zprec), S1), A), Aprec));
        _mm_store_si128((__m128i *)(LinePD), _mm_srai_epi32(S1, Zprec));
    }
}

  在测算完成后结果也会在那一个int类型的缓冲区中,之后再用SSE函数转换为int类型的。

     
前后一次那类别型的转移的SSE实现速度卓殊快,实现之后的提速也充裕肯定,对3000*2000的32位图像耗时大约由150ms降低到50ms,提速很扎眼。

     
但是对于24位如何做呢,他的乘除过程是3个字节重复的,无法直接采取SIMD的这种优化的不二法门的,同高斯模糊算法的完善优化过程分享(一) 一文类似,大家也是可以把24位的图像补一个Alpha通道然后再转换来int类型的缓冲区中,所以问题化解。

     
最难的是灰度图,因为灰度图的盘算过程是单字节重复的,正如上述代码所示,24位补一位的代价是多1个因素的计量,不过SIMD能四遍性统计4个整形的算法,由此依然很划算的,假设灰度也如此玩,SIMD的涨潮和浪费的测算句完全抵消了,而且还扩大了转移时间,肯定是不确切的,不过我们得以扭转思路,一行内各样要素之间的乘除是连续的,不过只要本身把连续4行的数据混搭为一行,混搭成类似32位这种数据格式,不就是能一向行使32位的算法了吧,最后再拆除回去就OK了。

     比例来说,四行灰度数据如下

     A1 A2 A3 A4 A5 A6 A7……

     B1 B2 B3 B4 B5 B6 B7……

     C1 C2 C3 C4 C5 C6 C7……

     D1 D2 D3 D4 D5 D6 D7……

  混搭为:

     A1 B1 C1 D1 A2 B2 C2 D2 A3 B3 C3 D3
A4 B4 C4 D4 A5 B5 C5 D5 A6 B6 C6 D6 A7 B7 C7 D7………

  
假如直白动用普通C语言混搭,这些过程或者非常耗时的,当然也非得的用SSE实现,大家假使条分缕析看过自己图像转置的SSE优化(帮助8位、24位、32位),提速4-6倍一文的代码,那么些过程实现也很容易。

     有的时候思路真的很重点。

     
在展开了地点的优化后,我曾自己满意过一段时间,因为她的时刻已经在自然水准上超越了SSE优化版本的Boxblur,可是俗话说,处处留心皆学问、开卷有益。当某一天自己留意到aprec的值为16抬高>>aprec这个操作时,大家脑海中就崩出了一个很好的SSE指令:_mm_mulhi_epi16,你们看,一个int类型右移16位不就是取int类型的高16位吗,而在移16位的先头就是个乘法,也就是要开展(a*b)>>16,这个和_mm_mulhi_epi16指令的意思完全一致。 

  可是采纳_mm_mulhi_epi16限令前,我们应该认可下这场景能无法满足数码范围的要求,我们看看需要优化的这句代码

       (Alpha * ((LinePD[X] <<
Zprec) – Buffer[X])) >> Aprec

   
 经过测试,只有radius小于2时,这么些alpha会超出short能发挥的上限,而(LinePD[X]
<< Zprec) –
Buffer[X])这句中LinePD[X]范围是[0,255],Zprec为7,两者相乘的界定不会超越32767,而Buffer[X]是个递归的量,只要第几遍不领先32767,前面就不会超过,由此双方的差也不会小于short能发挥的下限。所以说假若radius大于2,这多少个算式完全符合_mm_mulhi_epi16下令的要求。

   
 
由于_mm_mulhi_epi16两回性可以拍卖8个short类型,其他相应的SSE指令也同时改变为16位的话,理论上又会比用32位的SSE指令快一倍,更为紧要的是,我们最初的int缓冲区也应当改为short类型的缓冲区,对于这种自己耗时就不太多的算法,LOAD和Store指令的耗时是非常值得注意,使用short类型时这一个和内存打交道的效用又联合提升了。

   
 值得注意的是改为16位后,无论是32位、24位仍旧灰度的,写入到缓冲区的多少格式都会有有关的变更(其实仍然有这么些过多技巧我这里没有公布的)。

     最终:3000*2000的灰度图的履行时间为
7-8ms,提升了7倍左右。

   
 本文不享受最后优化的代码,请各位参考本文有关思路自行实现。

     一个测试相比工程:

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

   图片 1

     上述界面里的算法都是经过了SSE优化的,如今间接在研讨这上头的事物,又心得就会到此处来记录一下。

 图片 2