【AI PC 端算法优化】一步步优化自然饱和度算法

技术讨论 kira ⋅ 于 5个月前 ⋅ 856 阅读

文章来源: BBuf GiantPandaCV@微信公众号

上一节的RGB转灰度图算法我又做了两个相关优化,加入了多线程以及去掉了上次SSE计算中的一些重复计算,现在相对于传统实现已经可以获得4倍加速。同时我也在做一个AVX2的优化,所以不久后我将发布一个RGB转灰度图算法优化的升级版,尝试触摸这一个算法的优化极限,我会尽快做完实验发出来的。今天我先介绍一个有趣的自然饱和度算法,并讲解如何一步步进行优化。

上节内容:【AI PC 端算法优化】一步步优化 RGB 转灰度图算法

1. 原始实现

今天要介绍的自然饱和度算法是一个开源图像处理软件PhotoDemon(地址:https://github.com/tannerhelland/PhotoDemon)上的,原版是C#的,代码如下:

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

然后将其翻译为C++就获得了原始实现,代码如下:

//Adjustment如果为正值,会增加饱和度void VibranceAlgorithm_FLOAT(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Adjustment) { float VibranceAdjustment = -0.01 * Adjustment; 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 = max(max(Blue, Green), Red);   float AmtVal = (abs(Max - Avg) / 127.0f) * VibranceAdjustment;   if (Blue != Max) Blue += (Max - Blue) * AmtVal;   if (Green != Max) Green += (Max - Green) * AmtVal;   if (Red != Max) Red += (Max - Red) * AmtVal;   if (Red < 0) Red = 0;   else if (Red > 255) Red = 255;   if (Green < 0) Green = 0;   else if (Green > 255) Green = 255;   if (Blue < 0) Blue = 0;   else if (Blue > 255) Blue = 255;   LinePD[0] = Blue;   LinePD[1] = Green;   LinePD[2] = Red;   LinePS += 3;   LinePD += 3;  } }}

代码看起来非常简单,我们可以使用这个代码去对人像进行处理,效果如下:
file

原图
file

Adjustment=50,面色红润有精神
file
Adjustment=-50,面色苍白

接下来看一下这个算法原始实现的速度测试:
file

2. 自然饱和度算法优化第一版

首先,我们可以考虑去掉算法中的浮点运算,即是将float AmtVal = (abs(Max - Avg) / 127.0f) * VibranceAdjustment;这里的127.0f优化为乘法,怎么优化为乘法呢?这里是接近的,如果我们把它换成,那么我们可以用位运算来代替这个除法,实测将换成基本不影响算法的效果,所以这里直接采用了这种优化技巧。另外,Adjustment默认的范围为[-100,100],如果把它的范围线性扩大一些,比如扩大倍变成,这样在最后我们一次性移位,减少中间的损失。再然后,我们将这VibranceAdjustment里面的*0.01变成*0.01=1.28/128,然后把128放到下面的计算中并将VibranceAdjustment 重新设置为:int VibranceAdjustment = -1.28 * Adjustment;。最后还有一个点就是「这个算法中的绝对值运算完全可以去掉,因为平均值肯定是小于最大值的」。可能有点小晕哈,但是看代码很容易就理解了,下面给出优化后的代码:

void VibranceAlgorithm_INT(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Adjustment) { int VibranceAdjustment = -1.28 * Adjustment; 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) * VibranceAdjustment;   if (Blue != Max) Blue += (((Max - Blue) * AmtVal) >> 14);   if (Green != Max) Green += (((Max - Green) * AmtVal) >> 14);   if (Red != Max) Red += (((Max - Red) * AmtVal) >> 14);   if (Red < 0) Red = 0;   else if (Red > 255) Red = 255;   if (Green < 0) Green = 0;   else if (Green > 255) Green = 255;   if (Blue < 0) Blue = 0;   else if (Blue > 255) Blue = 255;   LinePD[0] = Blue;   LinePD[1] = Green;   LinePD[2] = Red;   LinePS += 3;   LinePD += 3;  } }}

下面看一下速度测试:
file
可以看到稍加优化,速度快了近2倍,还是比较可观的。

3. 自然饱和度算法优化第二版

在上面算法的基础上如果使用多线程(OpenMP)来优化的话那么会获得多少加速呢?我们来试试,源码如下:

void VibranceAlgorithm_INT_OpenMP(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Adjustment) { int VibranceAdjustment = -1.28 * Adjustment; for (int Y = 0; Y < Height; Y++) {  unsigned char *LinePS = Src + Y * Stride;  unsigned char *LinePD = Dest + Y * Stride;#pragma omp parallel for num_threads(4)  for (int X = 0; X < Width; X++) {   int Blue, Green, Red, Max;   Blue = LinePS[X*3 + 0], Green = LinePS[X*3 + 1], Red = LinePS[X*3 + 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) * VibranceAdjustment;   if (Blue != Max) Blue += (((Max - Blue) * AmtVal) >> 14);   if (Green != Max) Green += (((Max - Green) * AmtVal) >> 14);   if (Red != Max) Red += (((Max - Red) * AmtVal) >> 14);   if (Red < 0) Red = 0;   else if (Red > 255) Red = 255;   if (Green < 0) Green = 0;   else if (Green > 255) Green = 255;   if (Blue < 0) Blue = 0;   else if (Blue > 255) Blue = 255;   LinePD[X*3 + 0] = Blue;   LinePD[X*3 + 1] = Green;   LinePD[X*3 + 2] = Red;  } }}

我们来看一下加速效果:
file
可以看到使用OpenMP开启4线程,可以将我们的算法又优化接近2倍,仍然是可观的。接下来我们开始今天的主角,使用SSE指令集对这段代码进行优化。

4. 自然饱和度算法优化第三版

注意,在这个例子中,我们一次性加载48个图像数据到内存中,刚好可以放在3个__m128i变量中,同时看了我第一篇优化的人应该知道48正好被3整除,也就是说我们完整的加载了16个24位像素,这不会出现上一篇文章中的断层现象,使得下面48个像素可以和现在的48个像素使用同样的方法进行处理。上篇文章传送门:【AI PC端算法优化】一,一步步优化RGB转灰度图算法

首先,对于这样单像素点且邻域无关的算法,为了利用SSE提高运行速度,一个核心步骤就是把各个颜色分量分离为单独的连续的变量。然后在计算完之后,我们又需要把单独连续的变量重新分解成BGR(「注意OpenCV默认读图方式是BGR」)的形式,这两部分的代码实现如下:

  • BGRBGR->BBGGRR

    Src1 = _mm_loadu_si128((m128i )(LinePS + 0)); //B1,G1,R1,B2,G2,R2,B3,G3,R3,B4,G4,R4,B5,G5,R5,B6Src2 = _mm_loadu_si128((__m128i )(LinePS + 16));//G6,R6,B7,G7,R7,B8,G8,R8,B9,G9,R9,B10,G10,R10,B11,G11Src3 = _mm_loadu_si128((m128i *)(LinePS + 32));//R11,B12,G12,R12,B13,G13,R13,B14,G14,R14,B15,G15,R15,B16,G16,R16Blue8 = _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)));

  • BBGGRR->BGRBGR

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

这两个过程就是巧妙的利用__mm_shuffle_epi8指令完成的,而那个_mm_or_si128指令就是实现了将这次操作时所有的B或者G或者R放在了一个SSE向量中,对照后面的注释就很好理解了,也可以自己手推一下这个过程。如果想看详细步骤可以参考ImageShop大佬的博客,链接如下:https://www.cnblogs.com/Imageshop/p/7234463.html

接下来我们看看其它代码的实现,由于uchar数据类型的表示范围非常有限,除了少数几个操作能针对字节类型直接处理外(比如这段代码中的求RGB的Max值,就可以直接用下面的SIMD指令实现)

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

其他的一些操作无法在这样的范围(uchar)内进行了,所以我们需要将数据扩展到short或者int/float范围内,在SSE中完成这种操作是有直接命令的,例如byte扩展到short,则可以用_mm_unpacklo_epi8_mm_unpackhi_epi8配合Zero来实现:

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

其中_mm_unpacklo_epi8是将两个__m128i的低8位交错布置形成一个新的128位数据,如果其中一个参数为0,则就是把另外一个参数的低8个字节无损的扩展为16位了,以上述BL16为例,其内部布局为:
file

接下来我们需要来实现int Avg = (Blue + Green + Green + Red) >> 2;这行代码了,可以看到里的计算是无法在uchar范围内完成的,中间的Blue + Green + Green + Red在大部分情况下都会超出255并且绝对小于,因此我们需要扩展数据到16位(short),按照上面介绍的指令集对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);

上面的分析都非常常规,接下来要到本文的硬核部分了,首先SSE对于跳转处理是很不好做的,它擅长的是序列化的处理一件事情,虽然SSE提供了比较指令,但很多复杂的跳转情况,SSE仍然无能为力。而我们这段代码中跳转情况不算复杂,我们可以使用SSE中的比较指令得到个MaskMask中符合比较结果的值会为F,不符合的为0,然后把这个Mask和后面需要计算的某个值进行And操作,由于和F进行And操作不会改变操作数本身,和0进行And操作则变为0。因此,这种操作方式是无论你符合条件与否,都要进行后面的计算,只是不符合条件的计算并不会影响结果,这种多余的计算可能会低效SSE优化的部分提速效果,这个就要具体情况具体分析了。

然后参考ImageShop博主的思路,仔细观察我们的代码可以发现,这里跳转语句的本意是如果最大值和某个分量的值不相同,则进行后面的调整操作,否则就不进行调整。而后面的调整操作中有最大值减去这个分量的逻辑,也就是说如果满足条件两者相减则为,调整量此时也为,并不会对结果产生影响。「基于此,我们可以直接把这个条件判断去掉,这并不会影响结果。」 同时我们能节省一些SSE指令,并且也更加适合SSE的处理。

继续分析,由于代码中有((Max - Blue) * AmtVal) >> 14这行逻辑,其中AmtVal = (Max - Avg) * Adjustment,展开即为:

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

这三个数据相乘很大程度上会超出short所能表达的范围,因此,我们还需要对上面的位数据进行扩展,扩展到位,这样就多了很多指令,那么有没有不需要扩展的方式呢?ImageShop博主提出了一种方式,这里搬运一下:

在SSE指令集中有一个指令叫做_mm_mulhi_epi16,我们看一下这个指令能干什么?
file

_mm_mulhi_epi16 指令

这个指令剋一次性处理8个16位的数据,其计算结果相当于对于),但和必须是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类型所能表达的范围内。所以SSE代码就呼之欲出了。

  • 原始代码段。

    int AmtVal = (Max - Avg)  Adjustment;                                //    Get adjusted averageif (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));

最后一步就是将获得的B8,G8,R8别转换为不连续存储的形式即是BGR的格式,然后再存储即可。

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)));_mm_storeu_si128((__m128i *)(LinePD + 0), Dest1);_mm_storeu_si128((__m128i *)(LinePD + 16), Dest2);_mm_storeu_si128((__m128i *)(LinePD + 32), Dest3);

完整代码实现请看:https://github.com/BBuf/Image-processing-algorithm-Speed/blob/master/speed_rgb2gray_sse.cpp

我们来看一下速度测试:
file

5. 结论

======

这篇文章介绍了如何一步步优化一个自然饱和度算法,从原始算法的115.36ms优化到了13.04ms,「加速比达到了9.09倍」,还是比较可观的。

6. 参考

======

成为第一个点赞的人吧 :bowtie:
回复数量: 0
暂无回复~
您需要登陆以后才能留下评论!