超高速指数模糊算法的实现和优化(10000*10000在100ms左右实现)。

Posted 只(挚)爱图像处理

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了超高速指数模糊算法的实现和优化(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

   

     上述界面里的算法都是经过了SSE优化的,最近一直在研究这方面的东西,又心得就会到这里来记录一下。

 

 

  

以上是关于超高速指数模糊算法的实现和优化(10000*10000在100ms左右实现)。的主要内容,如果未能解决你的问题,请参考以下文章

Codecs系列转载:手淘H265编解码算法与工程优化

计算指数函数的算法

高斯模糊算法的全面优化过程分享。

开盘指数预测基于matlab遗传算法优化小波神经网络股票开盘指数预测含Matlab源码 2014期

O效率的表面模糊算法优化。

优化算法之指数加权平均详解