正文

SSE图像算法优化系列八：自然饱和度(Vibrance)算法的模拟实现及其SSE优化（附源码，可作为SSE图像入门，Vibrance算法也可用于简单的肤色调整）。

IP云
IP云V管理员/2 阅读
08/03
小程序：扫一扫查出行
【扫一扫了解最新限行尾号】
复制小程序

Vibrance这个单词搜索翻译一般振动，抖动或者是响亮、活力，但是官方的词汇里还从来未出现过自然饱和度这个词，也不知道当时的Adobe中文翻译人员怎么会这样处理。但是我们看看PS对这个功能的解释：

       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.

确实是和饱和度有关的，这样理解中文的翻译反而倒是合理，那么只能怪Adobe的开发者为什么给这个功能起个名字叫Vibrance了。

闲话不多说了，其实自然饱和度也是最近几个版本的PS才出现的功能，在调节有些图片的时候会有不错的效果，也可以作为简单的肤色调整的一个算法，比如下面这位小姑娘，用自然饱和度即可以让她失血过多，也可以让他肤色红晕。

          

原图                                                                                                  面色苍白                                                                                   肤色红晕一点

那么这个算法的内在是如何实现的呢，我没有仔细的去研究他，但是在开源软件PhotoDemon-master（开源地址：https://github.com/tannerhelland/PhotoDemon，visual basic 6.0的作品，我的最爱）提供了一个有点相似的功能，我们贴出他对改效果的部分注释：

'***************************************************************************
'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官方文档的描述有类似之处。

我们在贴出他的核心代码：

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

很简单的算法，先求出每个像素RGB分量的最大值和平均值，然后求两者之差，之后根据输入调节量求出调整量。

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 = ; Y < Height; Y++)
{
unsigned char * LinePS = Src + Y * Stride;
unsigned char * LinePD = Dest + Y * Stride;
for (int X = ; X < Width; X++)
{
int Blue = LinePS[], Green = LinePS[], Red = LinePS[];
int Avg = (Blue + Green + Green + Red) >> ;
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[] = IM_ClampToByte(Blue);
LinePD[] = IM_ClampToByte(Green);
LinePD[] = IM_ClampToByte(Red);
LinePS += ;
LinePD += ;
}
}

这个的结果和PS的是比较接近的，最起码趋势是非常接近的，但是细节还是不一样，不过可以断定的是方向是对的，如果你一定要复制PS的结果，我建议你花点时间改变其中的一些常数或者计算方式看看。应该能有收获，国内已经有人摸索出来了。

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

第一步优化，去除掉不必要计算和除法，很明显，这一句是本段代码中耗时较为显著的部分

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

/127.0f可以优化为乘法，同时注意VibranceAdjustment在内部不变，可以把他们整合到循环的最外层，即改为：

float VibranceAdjustment = -0.01 * Adjustment / 127.0f;

再注意abs里的参数， Max - Avg，这有必要取绝对值吗，最大值难道会比平均值小，浪费时间，最后改为：

float AmtVal = (Max - Avg) * VibranceAdjustment;

  这是浮点版本的简单优化，如果不勾选编译器的SSE优化，直接使用FPU，对于一副3000*2000的24位图像耗时在I5的一台机器上运行用时大概70毫秒，但这不是重点。

我们来考虑某些近似和定点优化。

第一我们把/127改为/128，这基本不影响效果，同时Adjustment默认的范围为[-100,100]，把它也线性扩大一点，比如扩大1.28倍，扩大到[-128，128]，这样在最后我们一次性移位，减少中间的损失，大概的代码如下：

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 <= ) || (Height <= )) return IM_STATUS_INVALIDPARAMETER;
if (Channel != ) return IM_STATUS_INVALIDPARAMETER; Adjustment = -IM_ClampI(Adjustment, -, ) * 1.28; // Reverse the vibrance input; this way, positive values make the image more vibrant. Negative values make it less vibrant.
for (int Y = ; Y < Height; Y++)
{
unsigned char *LinePS = Src + Y * Stride;
unsigned char *LinePD = Dest + Y * Stride;
for (int X = ; X < Width; X++)
{
int Blue, Green, Red, Max;
Blue = LinePS[]; Green = LinePS[]; Red = LinePS[];
int Avg = (Blue + Green + Green + Red) >> ;
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) >> );
if (Green != Max) Green += (((Max - Green) * AmtVal) >> );
if (Red != Max) Red += (((Max - Red) * AmtVal) >> );
LinePD[] = IM_ClampToByte(Blue);
LinePD[] = IM_ClampToByte(Green);
LinePD[] = IM_ClampToByte(Red);
LinePS += ;
LinePD += ;
}
}
return IM_STATUS_OK;
}

这样优化后，同样大小的图像算法用时35毫秒，效果和浮点版本的基本没啥区别。

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

对于这种单像素点、和领域无关的图像算法，为了能利用SSE提高程序速度，一个核心的步骤就是把各颜色分量分离为单独的连续的变量，对于24位图像，我们知道图像在内存中的布局为：

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748
B1G1R1B2G2R2B3G3R3B4G4R4B5G5R5B6G6R6B7G7R7B8G8R8B9G9R9B10G10R10B11G11R11B12G12R12B13G13R13B14G14R14B15G15R15B16G16R16

我们需要把它们变为：

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748
B1B2B3B4B4B5B7B8B9B10B11B12B13B14B15B16G1G2G3G4G5G6G7G8G9G10G11G12G13G14G15G16R1R2R3R4R5R6R7R8R9R10R11R12R13R14R15R16

处理完后我们又要把他们恢复到原始的BGR布局。

为了实现这个功能，我参考了采石工大侠的有关代码，分享如下：

我们先贴下代码：

Src1 = _mm_loadu_si128((__m128i *)(LinePS + ));
Src2 = _mm_loadu_si128((__m128i *)(LinePS + ));
Src3 = _mm_loadu_si128((__m128i *)(LinePS + )); Blue8 = _mm_shuffle_epi8(Src1, _mm_setr_epi8(, , , , , , -, -, -, -, -, -, -, -, -, -));
Blue8 = _mm_or_si128(Blue8, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-, -, -, -, -, -, , , , , , -, -, -, -, -)));
Blue8 = _mm_or_si128(Blue8, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-, -, -, -, -, -, -, -, -, -, -, , , , , ))); Green8 = _mm_shuffle_epi8(Src1, _mm_setr_epi8(, , , , , -, -, -, -, -, -, -, -, -, -, -));
Green8 = _mm_or_si128(Green8, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-, -, -, -, -, , , , , , , -, -, -, -, -)));
Green8 = _mm_or_si128(Green8, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-, -, -, -, -, -, -, -, -, -, -, , , , , ))); Red8 = _mm_shuffle_epi8(Src1, _mm_setr_epi8(, , , , , -, -, -, -, -, -, -, -, -, -, -));
Red8 = _mm_or_si128(Red8, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-, -, -, -, -, , , , , , -, -, -, -, -, -)));
Red8 = _mm_or_si128(Red8, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-, -, -, -, -, -, -, -, -, -, , , , , , )));

首先，一次性加载48个图像数据到内存，正好放置在三个__m128i变量中，同时另外一个很好的事情就是48正好能被3整除，也就是说我们完整的加载了16个24位像素，这样就不会出现断层，只意味着下面48个像素可以和现在的48个像素使用同样的方法进行处理。

如上代码，则Src1中保存着:

12345678910111213141516
B1G1R1B2G2R2B3G3R3B4G4R4B5G5R5B6

Src2中保存着:

12345678910111213141516
G6R6B7G7R7B8G8R8B9G9R9B10G10R10B11G11

Src3中的数据则为：

12345678910111213141516
R11B12G12R12B13G13R13B14G14R14B15G15R15B16G16R16

为了达到我们的目的，我们就要利用SSE中强大的shuffle指令了，如果能够把shuffle指令运用的出神入化，可以获取很多很有意思的效果，有如鸠摩智的小无相功一样，可以催动拈花指发、袈裟服魔攻等等，成就世间能和我鸠摩智打成平成的没有几个人一样的丰功伟绩。哈哈，说远了。

简单的理解shuffle指令，就是将__m128i变量内的各个数据按照指定的顺序进行重新布置，当然这个布置不一定要完全利用原有的数据，也可以重复某些数据，或者某些位置无数据，比如在执行下面这条指令

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中的数据为：

12345678910111213141516
B1B2B3B4B5B600000 0 0 00  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))得到的临时数据为：

12345678910111213141516
000000B7B8B9B10B1100000

如果把这个临时结果和之前的Blue8进行或操作甚至直接进行加操作，新的Blue8变量则为：

12345678910111213141516
B1B2B3B4B5B6B7B8B9B10B1100000

最后这一句和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临时的得到的变量为：

12345678910111213141516
00000000000B12B13B14B15B16

再次和之前的Blue8结果进行或操作得到最终的结果：

12345678910111213141516
B1B2B3B4B5B6B7B8B9B10B11B12B13B14B15B16

对于Green和Red分量，处理的方法和步骤是一样的，只是由于位置不同，每次进行shuffle操作的常数有所不同，但原理完全一致。

如果理解了由BGRBGRBGR ---》变为了BBBGGGRRR这样的模式的原理后，那么由BBBGGGRRR-->变为BGRBGRBGR的道理就非常浅显了，这里不赘述，直接贴出代码：

Dest1 = _mm_shuffle_epi8(Blue8, _mm_setr_epi8(, -, -, , -, -, , -, -, , -, -, , -, -, ));
Dest1 = _mm_or_si128(Dest1, _mm_shuffle_epi8(Green8, _mm_setr_epi8(-, , -, -, , -, -, , -, -, , -, -, , -, -)));
Dest1 = _mm_or_si128(Dest1, _mm_shuffle_epi8(Red8, _mm_setr_epi8(-, -, , -, -, , -, -, , -, -, , -, -, , -))); Dest2 = _mm_shuffle_epi8(Blue8, _mm_setr_epi8(-, -, , -, -, , -, -, , -, -, , -, -, , -));
Dest2 = _mm_or_si128(Dest2, _mm_shuffle_epi8(Green8, _mm_setr_epi8(, -, -, , -, -, , -, -, , -, -, , -, -, )));
Dest2 = _mm_or_si128(Dest2, _mm_shuffle_epi8(Red8, _mm_setr_epi8(-, , -, -, , -, -, , -, -, , -, -, , -, -))); Dest3 = _mm_shuffle_epi8(Blue8, _mm_setr_epi8(-, , -, -, , -, -, , -, -, , -, -, , -, -));
Dest3 = _mm_or_si128(Dest3, _mm_shuffle_epi8(Green8, _mm_setr_epi8(-, -, , -, -, , -, -, , -, -, , -, -, , -)));
Dest3 = _mm_or_si128(Dest3, _mm_shuffle_epi8(Red8, _mm_setr_epi8(, -, -, , -, -, , -, -, , -, -, , -, -, )));

核心还是这些常数的选取。

以上是处理的第一步，看上去这个代码很多，实际上他们的执行时非常快的，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为例，其内部布局为：

12345678910111213141516
B10B20B30B30B40B50B60B70

如果我们需要进行在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) >> );
if (Green != Max) Green += (((Max - Green) * AmtVal) >> );
if (Red != Max) Red += (((Max - Red) * AmtVal) >> );

对应的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), ), AmtVal));
GL16 = _mm_adds_epi16(GL16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxL16, GL16), ), AmtVal));
RL16 = _mm_adds_epi16(RL16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxL16, RL16), ), 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), ), AmtVal));
GH16 = _mm_adds_epi16(GH16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxH16, GH16), ), AmtVal));
RH16 = _mm_adds_epi16(RH16, _mm_mulhi_epi16(_mm_slli_epi16(_mm_sub_epi16(MaxH16, RH16), ), 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.0C++,float优化版本C++定点版C++/SSE版
速度400ms70ms35ms5ms

上面的VB6.0的耗时是原作者的代码编译后的执行速度，如果我自己去用VB6.0去优化他的话，有信心能做到70ms以内的。

但无论如何，SSE优化的速度提升是巨大的。

结论：

       简单的分析了自然饱和度算法的实现，分享了其SSE实现的过程，对于那些刚刚接触SSE，想做图像处理的朋友有一定的帮助。

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

写的真的好累，休息去了，觉得对你有用的请给我买杯啤酒或者咖啡吧。