前言
在ImageShop的博客上看到这个算法，本着跟大佬一起学习SSE的目的，开始了这个博客。我的代码实现：https://github.com/BBuf/Image-processing-algorithm-Speed 。作者的源码失效了，这个代码是我参考着作者博客复现出的。
算法原理
实际上就是下面这个代码，具体原理不清楚，这篇博客的目的本来也就是学习如何SSE优化代码，所以也无伤大雅。
//Adjustment如果为正值，会增加饱和度
//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;
}
}
}
优化一
首先可以去掉算法里面的浮点运算，就是把float AmtVal = (abs(Max - Avg) / 127.0f) * VibranceAdjustment; 这里的/127.0优化为乘法，同时注意到VibranceAdjustment在程序中没有改变，所以可以把他们整合到循环外。最后修改为：float AmtVal = (abs(Max - Avg) / 127.0f) * VibranceAdjustment; ，之后考虑用一下/127似乎不太好处理，我们把/127改成/128，这基本不影响效果，同时Adjustment默认的范围为[-100,100]，把它也线性扩大一点，比如扩大1.28倍，扩大到[-128，128]，这样在最后我们一次性移位，减少中间的损失，所以使用了这些优化的代码如下：
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;
}
}
}
优化二-SSE
首先对于这样单像素点，且领域无关的算法，为了利用SSE提高程序运行速度，一个核心步骤就是把各个颜色分量分离为单独的连续的变量，这部分的原理解释请看：https://www.cnblogs.com/Imageshop/p/7234463.html 。然后在程序处理完后，我们又需要把单独连续的变量分解为RGB的形式。这2个过程的代码如下：
RGBRGB—>RRGGBB
Src1 = _mm_loadu_si128((__m128i *)(LinePS + 0)); //B1,G1,R1,B2,G2,R2,B3,G3,R3,B4,G4,R4,B5,G5,R5,B6
Src2 = _mm_loadu_si128((__m128i *)(LinePS + 16));//G6,R6,B7,G7,R7,B8,G8,R8,B9,G9,R9,B10,G10,R10,B11,G11
Src3 = _mm_loadu_si128((__m128i *)(LinePS + 32));//R11,B12,G12,R12,B13,G13,R13,B14,G14,R14,B15,G15,R15,B16,G16,R16
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)));
RRGGBB---->RGBRGB
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操作来完成的。
接下来来看其他语句的实现，由于字节数据类型的表达范围非常有限，除了少有的几个有限的操作能针对字节类型直接处理外，比如本例的丘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为例，其内部布局为：
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| B0 | 0 | B1 | 0 | B2 | 0 | B3 | 0 | B4 | 0 | B5 | 0 | B6 | 0 | B7 | 0 |
知道了这一个SIMD操作后，我们就需要来实现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);
接下来带有循环的判断和计算AmtVal是比较麻烦的，这里参考ImageShop的优化方式如下：
下面这一段转载自ImageShop的博客：
我们知道，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位，这样就多了很多指令，那么有没有不需要扩展的方式呢。经过一番思索，我提出了下述解决方案：
在超高速指数模糊算法的实现和优化（1000010000在100ms左右实现 一文中，我在文章最后提到了终极的一个指令：_mm_mulhi_epi16（a,b），他能一次性处理8个16位数据，其计算结果相当于对于(ab）>>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));
最后就是将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);
完整代码并带有速度测试函数请看我的github: https://github.com/BBuf/Image-processing-algorithm-Speed/blob/master/speed_vibrance_algorithm.cpp
速度测试
| 优化方式 | 图像分辨率 | 速度 |
|---|---|---|
| C语言实现+单线程 | 4032*3024 | 70.17ms |
| 浮点数改成整形运算+单线程 | 4032*3024 | 36.30ms |
| SSE优化+单线程 | 4032*3024 | 8.72ms |
京公网安备 11010502036488号