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位，按上述办法，对Blue8Green8Red8Max8进行扩展，如下所示：

    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

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