语言SSE图像算法优化系列八:自然饱和度(Vibrance)算法的模拟实现及其SSE优化(附源码,可看做SSE图像入门,Vibrance算法也只是用来简单的肤色调整)。

     
今天我们来花点时间再谈谈一个歪曲算法,一个特级简单不过同时超级牛逼的算法,无论以职能上还是快上且足以与Boxblur,
stackblur或者是Gaussblur想媲美,效果达到,比Boxblur来的更平整,和Gaussblur相似,速度直达,经过自家之优化,在PC端比他们三个都使赶早一格外截,而且着力未需占用额外的内存,实在是一个断好之算法。

  Vibrance这个单词搜索翻译一般震荡,抖动或者是嘹亮、活力,但是官方的词汇里还常有未出现了当饱和度这个词,也未理解这之Adobe中文翻译人员怎么会如此处理。但是咱看PS对是功能的说明:

     
算法的中心并无是自家想到要发明的,是一个恋人以github上凿到之,率属于Cairo这个2D图形库的开源代码,详见:

       Vibrance: Adjusts the
saturation so that clipping is minimized as colors approach full
saturation. This setting changes the saturation of all lower-saturated
colors with less effect on the higher-saturated colors. Vibrance also
prevents skin tones from becoming oversaturated.

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

     
 确实是跟饱和度有关的,这样敞亮中文的翻反而却合理,那么只能怪Adobe的开发者为什么叫此意义于个名让Vibrance了。

       我们称之为Exponential
blur(指数模糊)。

     
 闲话不多说了,其实当饱和度也是近日几个本子的PS才面世的作用,在调试有些图片的下会产生不易的职能,也可看做简单的肤色调整的一个算法,比如下面这号小姐,用自饱和度即可以为它们失血过多,也得于他肤色红晕。

     
 提炼出这模糊的核心组成部分算法主要是下边这函数:

     
 语言 1   
 语言 2   
 语言 3

   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生成的片段稳的系数。

     
 那么这算法的内在是何等促成之吧,我从没仔细的失去研究他,但是以开源软件PhotoDemon-master(开源地址:https://github.com/tannerhelland/PhotoDemon,visual
basic
6.0的创作,我的最为爱)提供了一个小相像之效力,我们贴起他针对转移效果的一部分注释:

  类似于高斯歪曲或者StackBlur,算法为是属行列分离之,行方向上,先用上述方法从错误为右侧计算,然后于自右边为左,接着进行排方向处理,先从上向下,然后在从下向上。当然,行列的计为得以掉。需要小心的凡,每一样步都是为此事先处理过的数量进行的。

'***************************************************************************
'Vibrance Adjustment Tool
'Copyright 2013-2017 by Audioglider
'Created: 26/June/13
'Last updated: 24/August/13
'Last update: added command bar
'
'Many thanks to talented contributer Audioglider for creating this tool.
'
'Vibrance is similar to saturation, but slightly smarter, more subtle. The algorithm attempts to provide a greater boost
' to colors that are less saturated, while performing a smaller adjustment to already saturated colors.
'
'Positive values indicate "more vibrance", while negative values indicate "less vibrance"
'
'All source code in this file is licensed under a modified BSD license. This means you may use the code in your own
' projects IF you provide attribution. For more information, please visit http://photodemon.org/about/license/
'
'***************************************************************************

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

 其中的叙说和PS官方文档的讲述有类似之处。

     水平反向的处理:

   我们在贴起他的中坚代码:

  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);
  }
     For x = initX To finalX
        quickVal = x * qvDepth
        For y = initY To finalY
            'Get the source pixel color values
            r = ImageData(quickVal + 2, y)
            g = ImageData(quickVal + 1, y)
            b = ImageData(quickVal, y)

            'Calculate the gray value using the look-up table
            avgVal = grayLookUp(r + g + b)
            maxVal = Max3Int(r, g, b)

            'Get adjusted average
            amtVal = ((Abs(maxVal - avgVal) / 127) * vibranceAdjustment)

            If r <> maxVal Then
                r = r + (maxVal - r) * amtVal
            End If
            If g <> maxVal Then
                g = g + (maxVal - g) * amtVal
            End If
            If b <> maxVal Then
                b = b + (maxVal - b) * amtVal
            End If

            'Clamp values to [0,255] range
            If r < 0 Then r = 0
            If r > 255 Then r = 255
            If g < 0 Then g = 0
            If g > 255 Then g = 255
            If b < 0 Then b = 0
            If b > 255 Then b = 255

            ImageData(quickVal + 2, y) = r
            ImageData(quickVal + 1, y) = g
            ImageData(quickVal, y) = b
        Next
    Next

  垂直方向的拍卖:

  很简短的算法,先要来每个像素RGB分量的最为老价值与平均值,然后求两者之异,之后据悉输入调节量求来调整量。

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

     
 VB的语法有些人或未熟识,我多少开点转翻译成C的代码如下:

  最终的歪曲算法如下所示:

    float VibranceAdjustment = -0.01 * Adjustment;        //       Reverse the vibrance input; this way, positive values make the image more vibrant.  Negative values make it less vibrant.
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char * LinePS = Src + Y * Stride;
        unsigned char * LinePD = Dest + Y * Stride;
        for (int X = 0; X < Width; X++)
        {
            int Blue = LinePS[0],    Green = LinePS[1],    Red = LinePS[2];
            int Avg = (Blue + Green + Green + Red) >> 2;
            int Max = IM_Max(Blue, IM_Max(Green, Red));
            float AmtVal = (abs(Max - Avg) / 127.0f) * VibranceAdjustment;                        //    Get adjusted average
            if (Blue != Max)    Blue += (Max - Blue) * AmtVal;
            if (Green != Max)    Green += (Max - Green) * AmtVal;
            if (Red != Max)    Red += (Max - Red) * AmtVal;
            LinePD[0] = IM_ClampToByte(Blue);
            LinePD[1] = IM_ClampToByte(Green);
            LinePD[2] = IM_ClampToByte(Red);
            LinePS += 3;
            LinePD += 3;
        }
    }
  // 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;
  }

  这个的结果以及PS的是于接近的,最起码趋势是格外相近的,但是细节要无平等,不过可以判的是趋势是指向的,如果你势必要是复制PS的结果,我建议乃花点时间改中的有的常数或者计算方法省。应该能够生出得,国内都有人索出来了。

  作为一个突出的利用,或者说尽量减少参数,常用之aprec取值为16,Zprec
取值为7。

     
我们重点讲下这个算法的优化及其SSE实现,特别是SSE版本代码是本文的基本点。

   
 回顾下代码,整体过程遭到除alpha参数的计涉及到了浮点,其他一些还是整形的乘法和倒操作,因此可想像,速度应该不慢,而且非常适合于手机端处理。同时注意到_blurrow和_blurcol函数循环明显相互之间是独立的,可以行使基本上线程并行处理,但是这个代码主要是把注于算法的抒发,并没了多的设想再好的效率。

      第一步优化,去解不必要计算和除法,很显,这同样句子是本段代码中耗时较为强烈的片

   
 另外一些,很醒目,算法的耗时凡是同Radius参数没有其他涉及的,也就是说这为是个O(1)算法。

        float AmtVal = (abs(Max –
Avg) / 127.0f) * VibranceAdjustment;     

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

  /127.0f可以优化为乘法,同时注意VibranceAdjustment在里不移,可以把他们成及循环的太外层,即改变吧:

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;
    }
}
      float VibranceAdjustment = -0.01 * Adjustment / 127.0f;

  在 高斯模糊算法的宏观优化过程分享(一) 中我们追究了垂直方向处理算法一般不宜直接写,而应当据此一个现之行缓存进行拍卖,这样列方向的灰度图的处理代码类似下面的:

  再理会abs里的参数, Max –
Avg,这生必不可少取绝对值吗,最深价值难道会比平均值小,浪费时间,最后移呢:

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);
      float AmtVal = (Max - Avg) * VibranceAdjustment;

   修改也上述后,测试一个3000*2000之8号灰度图,耗时约52ms(未采取多线程的),和通常的C语言实现的Boxblur时间大多。

    这是浮点版本的概括优化,如果未勾选编译器的SSE优化,直接下FPU,对于一副3000*2000之24各项图像耗时在I5的等同华机械上运行用时大概70毫秒,但当下不是至关重要。

     
 除了线程外,这个日子是不是还有改进之长空为,我们先来探望排方向的优化。

  我们来设想某些近似和固化优化。

   以列方向的  for (int X = 0; X <
Width; X++)
循环内,我们注意到对Buffer的每个元素的处理还是单独和同一之,很引人注目这样的历程是老大容易采取SIMD指令优化的,但是循环体内部发生一对凡unsigned
char类型的数目,为运用SIMD指令,需要转移为int类型较为有利,而结尾保存时又要重新处理为unsigned
char类型的,这种往返换的耗时和其他计量的涨潮能否来带效益为,我们进行了代码的编纂,比如:

     
 第一咱拿/127改吗/128,这基本未影响效果,同时Adjustment默认的限制为[-100,100],把它们为线性扩大一点,比如扩大1.28倍增,扩大至[-128,128],这样以结尾我们一次性移位,减少中间的损失,大概的代码如下:

    for (int X = 0; X < Width; X++)        //  从上到下
    {
        Buffer[X] += (Alpha * ((LinePS[X] << Zprec) - Buffer[X])) >> Aprec;
        LinePD[X] = Buffer[X] >> Zprec;
    }
int IM_VibranceI(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Adjustment)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                    return IM_STATUS_INVALIDPARAMETER;
    if (Channel != 3)                                    return IM_STATUS_INVALIDPARAMETER;

    Adjustment = -IM_ClampI(Adjustment, -100, 100) * 1.28;        //       Reverse the vibrance input; this way, positive values make the image more vibrant.  Negative values make it less vibrant.
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        unsigned char *LinePD = Dest + Y * Stride;
        for (int X = 0; X < Width; X++)
        {
            int Blue, Green, Red, Max;
            Blue = LinePS[0];    Green = LinePS[1];    Red = LinePS[2];
            int Avg = (Blue + Green + Green + Red) >> 2;
            if (Blue > Green)
                Max = Blue;
            else
                Max = Green;
            if (Red > Max)
                Max = Red;
            int AmtVal = (Max - Avg) * Adjustment;                                //    Get adjusted average
            if (Blue != Max)    Blue += (((Max - Blue) * AmtVal) >> 14);
            if (Green != Max)    Green += (((Max - Green) * AmtVal) >> 14);
            if (Red != Max)        Red += (((Max - Red) * AmtVal) >> 14);
            LinePD[0] = IM_ClampToByte(Blue);
            LinePD[1] = IM_ClampToByte(Green);
            LinePD[2] = IM_ClampToByte(Red);
            LinePS += 3;
            LinePD += 3;
        }
    }
    return IM_STATUS_OK;
}

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

  这样优化后,同样大小的图像算法用时35毫秒,效果以及浮点版本的骨干无啥区别。

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

     
 最后我们着重来讲讲SSE版本的优化。

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

  
对于这种单像素点、和领域无关的图像算法,为了能以SSE提高程序速度,一个为主的步骤就是是将每颜色分量分离为独立的连续的变量,对于24员图像,我们了解图像在内存中的布局也:

   
 当图像不是灰度模式时,对于直方向的拍卖和灰度不会见有分,这是为,只需要充实循环的长短就得了。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
B1 G1 R1 B2 G2 R2 B3 G3 R3 B4 G4 R4 B5 G5 R5 B6 G6 R6 B7 G7 R7 B8 G8 R8 B9 G9 R9 B10 G10 R10 B11 G11 R11 B12 G12 R12 B13 G13 R13 B14 G14 R14 B15 G15 R15 B16 G16 R16

   
 我们重新来看看水平方向的优化,当图像是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,提速很显。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
B1 B2 B3 B4 B4 B5 B7 B8 B9 B10 B11 B12 B13 B14 B15 B16 G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12 G13 G14 G15 G16 R1 R2 R3 R4 R5 R6 R7 R8 R9 R10 R11 R12 R13 R14 R15 R16

     
但是对24各类怎么收拾为,他的算计过程是3个字节重复的,无法直接以SIMD的这种优化的措施的,同高斯模糊算法的周全优化过程分享(一) 一温软类,我们吧是足以管24各的图像上一个Alpha通道然后重新换到int类型的缓冲区中,所以问题迎刃而解。

 

     
最为难之是灰度图,因为灰度图的计量过程是单字节重复的,正使上述代码所示,24各项上一各之代价是基本上1单因素的盘算,但是SIMD能一次性计算4独整形的算法,因此要生划算的,如果灰度也这样玩,SIMD的涨价与浪费之计句完全抵消了,而且还增加了变时,肯定是无得当的,但是咱得转变思路,一行内相继要素中的算计是连接的,但是一旦我拿连续4行的数额混搭为同实行,混搭成类似32各那种数据格式,不纵能直接行使32个之算法了吧,最后重复拆除回去就是OK了。

     

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

   
 处理了晚我们又如把他们恢复到原之BGR布局。

     A1 A2 A3 A4 A5 A6 A7……

   
 为了实现这个意义,我参考了采石工大侠的关于代码,分享如下:

     B1 B2 B3 B4 B5 B6 B7……

     我们先贴下代码:

     C1 C2 C3 C4 C5 C6 C7……

    Src1 = _mm_loadu_si128((__m128i *)(LinePS + 0));
    Src2 = _mm_loadu_si128((__m128i *)(LinePS + 16));
    Src3 = _mm_loadu_si128((__m128i *)(LinePS + 32));

    Blue8 = _mm_shuffle_epi8(Src1, _mm_setr_epi8(0, 3, 6, 9, 12, 15, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1));
    Blue8 = _mm_or_si128(Blue8, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, 2, 5, 8, 11, 14, -1, -1, -1, -1, -1)));
    Blue8 = _mm_or_si128(Blue8, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 4, 7, 10, 13)));

    Green8 = _mm_shuffle_epi8(Src1, _mm_setr_epi8(1, 4, 7, 10, 13, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1));
    Green8 = _mm_or_si128(Green8, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-1, -1, -1, -1, -1, 0, 3, 6, 9, 12, 15, -1, -1, -1, -1, -1)));
    Green8 = _mm_or_si128(Green8, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 2, 5, 8, 11, 14)));

    Red8 = _mm_shuffle_epi8(Src1, _mm_setr_epi8(2, 5, 8, 11, 14, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1));
    Red8 = _mm_or_si128(Red8, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-1, -1, -1, -1, -1, 1, 4, 7, 10, 13, -1, -1, -1, -1, -1, -1)));
    Red8 = _mm_or_si128(Red8, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 3, 6, 9, 12, 15)));

     D1 D2 D3 D4 D5 D6 D7……

     
首先,一次性加载48个图像数据到内存,正好停于三只__m128i变量中,同时另外一个不胜好之政工虽是48刚能被3理除,也就是说我们整体的加载了16个24员像素,这样尽管不见面产出断层,只表示下面48独像素可以与今天之48独像素使用同一的办法开展处理。

  混搭为:

     
如齐代码,则Src1负保留在:

     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………

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
B1 G1 R1 B2 G2 R2 B3 G3 R3 B4 G4 R4 B5 G5 R5 B6

  
如果直接以普通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限令的意完全一致。 

      Src2蒙保存在:

  但是以_mm_mulhi_epi16命前,我们应有认同下本场景能不能够满足数量范围的需要,我们省用优化的那句代码

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
G6 R6 B7 G7 R7 B8 G8 R8 B9 G9 R9 B10 G10 R10 B11 G11

       (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类型时是和内存打交道的效率又一道提高了。

  Src3负之数则也:

   
 值得注意的凡移吗16各项后,无论是32各、24各还是灰度的,写副到缓冲区底数据格式都见面发连锁的更改(其实还是产生不少多技我此没有发挥的)。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
R11 B12 G12 R12 B13 G13 R13 B14 G14 R14 B15 G15 R15 B16 G16 R16

     最终:3000*2000之灰度图的行时呢
7-8ms,提高了7倍左右。

 

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

 

     一个测试于工程:

    为了上我们的目的,我们即将动用SSE中强有力的shuffle指令了,如果能将shuffle指令运用的出神入化,可以落很多格外有趣的效用,有如鸠摩智的有点无相功一样,可以催动拈花指发、袈裟服魔攻等等,成就世间能同自我鸠摩智打成平成的无几只人一致的丰功伟绩。哈哈,说远了。

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

   
 简单的敞亮shuffle指令,就是以__m128i变量内之相继数据论指定的逐一进行再次布置,当然这布阵不必然要统统采用原来的多寡,也堪再次某些数据,或者某些位置多按部就班,比如在实行下这长长的指令

   语言 4

    Blue8 = _mm_shuffle_epi8(Src1,
_mm_setr_epi8(0, 3, 6, 9, 12, 15, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1));

     上述界面里之算法都是经了SSE优化的,最近径直于研讨这点的事物,又心得就见面交这里来记录转。

   Blue8中的数量也:

 语言 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
B1 B2 B3 B4 B5 B6 0 0 0 0 0  0  0  0  0

 

_mm_setr_epi8指令的参数顺序可能更适合于我们常用的从左到右的理解习惯,其中的某个参数如果不在0和15之间时,则对应位置的数据就是被设置为0。

 
 可以看看进展上述操作后Blue8的签6个字节已经入我们的需要了。

   于羁押代码的生一致词:

        Blue8 = _mm_or_si128(Blue8, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, 2, 5, 8, 11, 14, -1, -1, -1, -1, -1)));

  这句之继半部分和前的好像,只是内的常数不同,由_mm_shuffle_epi8(Src2,
_mm_setr_epi8(-1, -1, -1, -1, -1, -1, 2, 5, 8, 11, 14, -1, -1, -1,
-1, -1))得到的现数据为:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0 0 0 0 0 0 B7 B8 B9 B10 B11 0 0 0 0 0

 

 

   
 如果把这个临时结果及事先的Blue8进行或者操作还是直接进行加操作,新的Blue8变量则也:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 0 0 0 0 0

 

 

     
最后就无异于句子和Blue8相关的代码为:

Blue8 = _mm_or_si128(Blue8, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 4, 7, 10, 13)));

  后面的shuffle临时之获的变量为:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0 0 0 0 0 0 0 0 0 0 0 B12 B13 B14 B15 B16

 

 

   
 再次与前的Blue8结果进行或者操作得到终极之结果:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12 B13 B14 B15 B16

 

 

   
 对于Green和Red分量,处理的措施及手续是一模一样的,只是出于位置不同,每次进行shuffle操作的常数有所不同,但原理完全平等。

  如果掌握了是因为BGRBGRBGR
—》变为了BBBGGGRRR这样的模式之规律后,那么由BBBGGGRRR–>变为BGRBGRBGR的道理就异常浅显了,这里不赘述,直接贴起代码:

    Dest1 = _mm_shuffle_epi8(Blue8, _mm_setr_epi8(0, -1, -1, 1, -1, -1, 2, -1, -1, 3, -1, -1, 4, -1, -1, 5));
    Dest1 = _mm_or_si128(Dest1, _mm_shuffle_epi8(Green8, _mm_setr_epi8(-1, 0, -1, -1, 1, -1, -1, 2, -1, -1, 3, -1, -1, 4, -1, -1)));
    Dest1 = _mm_or_si128(Dest1, _mm_shuffle_epi8(Red8, _mm_setr_epi8(-1, -1, 0, -1, -1, 1, -1, -1, 2, -1, -1, 3, -1, -1, 4, -1)));

    Dest2 = _mm_shuffle_epi8(Blue8, _mm_setr_epi8(-1, -1, 6, -1, -1, 7, -1, -1, 8, -1, -1, 9, -1, -1, 10, -1));
    Dest2 = _mm_or_si128(Dest2, _mm_shuffle_epi8(Green8, _mm_setr_epi8(5, -1, -1, 6, -1, -1, 7, -1, -1, 8, -1, -1, 9, -1, -1, 10)));
    Dest2 = _mm_or_si128(Dest2, _mm_shuffle_epi8(Red8, _mm_setr_epi8(-1, 5, -1, -1, 6, -1, -1, 7, -1, -1, 8, -1, -1, 9, -1, -1)));

    Dest3 = _mm_shuffle_epi8(Blue8, _mm_setr_epi8(-1, 11, -1, -1, 12, -1, -1, 13, -1, -1, 14, -1, -1, 15, -1, -1));
    Dest3 = _mm_or_si128(Dest3, _mm_shuffle_epi8(Green8, _mm_setr_epi8(-1, -1, 11, -1, -1, 12, -1, -1, 13, -1, -1, 14, -1, -1, 15, -1)));
    Dest3 = _mm_or_si128(Dest3, _mm_shuffle_epi8(Red8, _mm_setr_epi8(10, -1, -1, 11, -1, -1, 12, -1, -1, 13, -1, -1, 14, -1, -1, 15)));

  核心要这些常数的选取。

     
以上是拍卖的率先步,看上去是代码很多,实际上他们之推行时杀快的,3000*2000之觊觎是拆分和合并过程也就算盖2ms。

     
当然是因为字节数据列的表达范围非常少,除了少有的几只片的操作会针对字节类型直接处理外,比如本例的丘RGB的Max值,就得直接用底的SIMD指令实现:

Max8 = _mm_max_epu8(_mm_max_epu8(Blue8, Green8), Red8);

     
很其他基本上计都是力不从心直接以如此的界定外进行了,因此就发必要将数据类型扩展,比如扩展及short类型或者int/float类型。

     
在SSE里展开这样的操作为是非常简单的,SSE提供了大量的数据类型转换的函数和下令,比如有byte扩展及short,则足以据此_mm_unpacklo_epi8和_mm_unpackhi_epi8配合zero来实现:

BL16 = _mm_unpacklo_epi8(Blue8, Zero);
BH16 = _mm_unpackhi_epi8(Blue8, Zero);

  其中

Zero = _mm_setzero_si128();

   很有趣的操作,比如_mm_unpacklo_epi8凡将鲜单__m128i之低8各项交错布置形成一个初的128各数据,如果中间一个参数为0,则就是将另外一个参数的低8独字节无损的扩展为16位了,以上述BL16啊条例,其里面布局也:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
B1 0 B2 0 B3 0 B3 0 B4 0 B5 0 B6 0 B7 0

 

 

  如果我们用开展以int范围外开展测算,则还得更扩大,此时得使_mm_unpackhi_epi16/_mm_unpacklo_epi16相当zero继续展开扩展,这样一个Blue8变量需要4个__m128i
int范围的数量来表达。

     
好,说道这里,我们继承羁押咱们C语言里之马上词:

  int Avg = (Blue + Green + Green + Red) >> 2;

  可以见见,这里的计是无法再byte范围外到位的,中间的Blue

  • Green + Green +
    Red在多数场面下还见面超出255如绝对小于255*4,,因此我们得扩大数据到16号,按上述办法,对Blue8\Green8\Red8\Max8展开扩张,如下所示:

      BL16 = _mm_unpacklo_epi8(Blue8, Zero);
      BH16 = _mm_unpackhi_epi8(Blue8, Zero);
      GL16 = _mm_unpacklo_epi8(Green8, Zero);
      GH16 = _mm_unpackhi_epi8(Green8, Zero);
      RL16 = _mm_unpacklo_epi8(Red8, Zero);
      RH16 = _mm_unpackhi_epi8(Red8, Zero);
      MaxL16 = _mm_unpacklo_epi8(Max8, Zero);
      MaxH16 = _mm_unpackhi_epi8(Max8, Zero);
    

  这计算Avg就回至渠道成了:

     AvgL16 = _mm_srli_epi16(_mm_add_epi16(_mm_add_epi16(BL16, RL16), _mm_slli_epi16(GL16, 1)), 2);
     AvgH16 = _mm_srli_epi16(_mm_add_epi16(_mm_add_epi16(BH16, RH16), _mm_slli_epi16(GH16, 1)), 2);

  中间两独Green相加是为此运动或一直相加对进度没有啥影响的。

     
 接下来的优化则是本例的一个表征有了。我们来详细分析。

     
 我们清楚,SSE对于跳转是坏无友善之,他万分擅长序列化处理一个政工,虽然他供了广大比指令,但是多情况下复杂的跳转SSE还是无为力,对于本例,情况较特别,如果要运用SSE的于指令也是可直接实现之,实现之法门时,使用比较指令得到一个Mask,Mask中称于结实的值会为FFFFFFFF,不适合的为0,然后把这个Mask和后要计算的之一值进行And操作,由于和FFFFFFFF进行And操作不见面改操作数本身,和0进行And操作则变为0,在群状下,就是无论你符合条件与否,都进展末端的乘除,只是不符合条件的计不会见潜移默化结果,这种计算可能会见失效SSE优化的片段提速效果,这个将要具体情况具体分析了。

     
注意观察本例的代码,他的本心是使尽要命价值与有分量的值不等同,则展开末端的调动操作,否则不进行调节。可后面的调整操作着发出极特别价值减去该分量的操作,也就是象征一旦尽要命价值和拖欠分量相同,两者相减则也0,调整量此时为为0,并无影响结果,也就一定给无调节,因此,把这个标准判断失丢,并无见面影响结果。同时考虑到实在情况,最要命价值在很多气象呢只有会以及有一个重量相同,也就是说只出1/3的几率不履行跳转后底语句,在本例中,跳反后底代码执行复杂度并无高,去丢这些条件判断用增加并代码所吃的性质和减少3独判断的岁月已经以一个水平及了,因此,完全可去除这些判断语句,这样就非常适合于SSE实现了。

  接着分析,由于代码中发生((Max –
Blue) * AmtVal) >> 14,其中AmtVal
= (Max – Avg) * Adjustment,展开就为:  ((Max – Blue) * (Max – Avg) *
Adjustment)>>14;这三个数据相乘很非常程度达到会见超出short所能发表的界定,因此,我们还欲对端的16各数据开展扩展,扩展至32个,这样就多矣森指令,那么来无起免欲扩大的方法啊。经过一番思考,我提出了下述解决方案:

   
在超越快指数模糊算法的落实同优化(10000*10000每当100ms左右实现 一中和被,我当文章最后提到了极的一个发令:_mm_mulhi_epi16(a,b),他能一次性处理8只16个数据,其计算结果一定给对(a*b)>>16,但此很明a和b必须是short类型所能够达的克。

       注意我们的此表达式:

              ((Max – Blue) * (Max –
Avg) * Adjustment)>>14

   
   首先,我们用他恢弘为活动16员之结果,变为如下:

         ((Max – Blue) * 4 * (Max –
Avg) * Adjustment)>>16

     
Adjustment我们已经用他限定于了[-128,128]中间,而(Max –
Avg)理论及的无比老价值为255 –
85=170,(即RGB分量有一个是255,其他的还也0),最小值为0,因此,两者在个别范围外之大成未会见胜出short所能够表达的限定,而(Max-Blue)的顶可怜价值也255,最小值为0,在乘机以4为以short类型所能表达的界定外。所以,下一致步你们知道了呢?

     
 经过上述分析,下面就四履C代码可由于下述SSE函数实现:

    int AmtVal = (Max - Avg) * Adjustment;                                //    Get adjusted average
    if (Blue != Max)    Blue += (((Max - Blue) * AmtVal) >> 14);
    if (Green != Max)    Green += (((Max - Green) * AmtVal) >> 14);
    if (Red != Max)        Red += (((Max - Red) * AmtVal) >> 14);

  对应的SSE代码为:

    AmtVal = _mm_mullo_epi16(_mm_sub_epi16(MaxL16, AvgL16), Adjustment128);
    BL16 = _mm_adds_epi16(BL16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxL16, BL16), 2), AmtVal));
    GL16 = _mm_adds_epi16(GL16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxL16, GL16), 2), AmtVal));
    RL16 = _mm_adds_epi16(RL16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxL16, RL16), 2), AmtVal));

    AmtVal = _mm_mullo_epi16(_mm_sub_epi16(MaxH16, AvgH16), Adjustment128);
    BH16 = _mm_adds_epi16(BH16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxH16, BH16), 2), AmtVal));
    GH16 = _mm_adds_epi16(GH16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxH16, GH16), 2), AmtVal));
    RH16 = _mm_adds_epi16(RH16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxH16, RH16), 2), AmtVal));

  最后一步就是是将这些16各的数额更转移为8个之,注意原始代码中出Clamp操作,这个操作实际是个耗时的历程,而SSE天然的具有抗饱和的函数。

  Blue8 = _mm_packus_epi16(BL16, BH16);
  Green8 = _mm_packus_epi16(GL16, GH16);
  Red8 = _mm_packus_epi16(RL16, RH16);

 _mm_packus_epi16这个的用法和含义自己去MSDN搜索一下吧,实在是懒得解释了。

   最终优化速度:5ms。

   来单速度比较:

版本 VB6.0 C++,float优化版本 C++定点版 C++/SSE版
速度 400ms 70ms 35ms 5ms

 

     

  上面的VB6.0的耗时凡是原作者的代码编译后的实行进度,如果我自己去用VB6.0去优化他的语,有信念会成功70ms以内的。

  但好歹,SSE优化的速度提升是高大的。

结论:

       简单的剖析了当然饱和度算法的落实,分享了该SSE实现的过程,对于那些刚接触SSE,想做图像处理的冤家起得的扶助。

     
  源代码下载地址:http://files.cnblogs.com/files/Imageshop/Vibrance.rar

       
写的真好累,休息去了,觉得对君有效之请吃自身购买杯啤酒或咖啡吧。

语言 6

 

发表评论

电子邮件地址不会被公开。 必填项已用*标注

网站地图xml地图