SSE图像算法优化系列八：自然饱和度(Vibrance)算法的模拟实现及其SSE优化（附源码，可作为SSE图像入门，Vibrance算法也可用于简单的肤色调整）。
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 + 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分量的最大值和平均值，然后求两者之差，之后根据输入调节量求出调整量。
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;
}
}
这个的结果和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 <= 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; }
这样优化后，同样大小的图像算法用时35毫秒，效果和浮点版本的基本没啥区别。
最后我们重点来讲讲SSE版本的优化。
对于这种单像素点、和领域无关的图像算法，为了能利用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 |
我们需要把它们变为：
| 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 |
处理完后我们又要把他们恢复到原始的BGR布局。
为了实现这个功能，我参考了采石工大侠的有关代码，分享如下：
我们先贴下代码：
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)));
首先，一次性加载48个图像数据到内存，正好放置在三个__m128i变量中，同时另外一个很好的事情就是48正好能被3整除，也就是说我们完整的加载了16个24位像素，这样就不会出现断层，只意味着下面48个像素可以和现在的48个像素使用同样的方法进行处理。
如上代码，则Src1中保存着:
| 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 |
Src2中保存着:
| 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 |
Src3中的数据则为：
| 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 |
为了达到我们的目的，我们就要利用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中的数据为：
| 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 | 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，想做图像处理的朋友有一定的帮助。
写的真的好累，休息去了，觉得对你有用的请给我买杯啤酒或者咖啡吧。
