SSE图像算法优化系列十九：一种局部Gamma校订对比度加强算法及其SSE优化。

　　这是一篇2010年比较古老的文章了，是在QQ群里一位群友提到的，无聊下载看了下，其实也没有啥高深的理论，抽空实现了下，虽然不高大上，仍是花了点时间和心思优化了代码，既然这样，就顺便分享下优化的思路和经历。

　　文章的名字为：Contrast image correction method，因为本人博客的后台文件已经快超过博客园所允许的最大空间，这里就不直接上传文章了，你们能够直接点我提供的连接下载。

　　文章的核心就是对普通的伽马校订作改进和扩展，通常来讲，伽马校订具备如下的标准形式：

　　

　　其中I(i,j)为输入图像，O(i,j)为输出图像，γ为控制参数，当γ大于1时，图像总体变亮，当γ小于1大于0时，图像总体变暗，γ小于0算法无心义。　　

　　这个算法对于图像总体偏暗或总体偏亮时，经过调节参数γ能够得到较为满意的效果，可是若是图像中同时存在欠曝或过曝的区域，同一个参数就没法同时满意的效果了，所以，可引入一种γ随图像局部区域信息变化的算法来获取更为满意的效果，一种经常使用的形式以下：

　　

　　Moroney在其论文Local colour correction using nonlinear masking提出了以下公式：

　　

 　　其中的mask获取方式为：先对原图进行反色处理，而后进行必定半径的高斯模糊。

　　这样作的道理以下：若是mask的值大于128，说明那个点是个暗像素同时周边也是暗像素，所以γ值须要小于0以便将其增亮，mask值小于128，对应的说明当前点是个较亮的像素，且周边像素也较亮，mask值为128则不产生任何变化，同时，mask值离128越远，校订的量就越大，而且还有个特色就是纯白色和纯黑色不会有任何变化（这其实也是会产生问题的）。

　　以下图所示，直观的反应了不一样的mask值的映射结果。

  

　　简单写一段测试代码，看看这个的效果如何：

int IM_LocalExponentialCorrection(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride) { unsigned char *Mask = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char)); IM_Invert(Src, Mask, Width, Height, Stride); // Invert Intensity
    IM_ExpBlur(Mask, Mask, Width, Height, Stride, 20);        // Blur 
    for (int Y = 0; Y < Height; Y++) { unsigned char *LinePS = Src + Y * Stride; unsigned char *LinePD = Dest + Y * Stride; unsigned char *LinePM = Mask + Y * Stride; for (int X = 0; X < Width; X++) { LinePD[0] = IM_ClampToByte(255 * pow(LinePS[0] * IM_INV255, pow(2, (128 - LinePM[0]) / 128.0f)));        // Moroney论文的公式
            LinePD[1] = IM_ClampToByte(255 * pow(LinePS[1] * IM_INV255, pow(2, (128 - LinePM[1]) / 128.0f))); LinePD[2] = IM_ClampToByte(255 * pow(LinePS[2] * IM_INV255, pow(2, (128 - LinePM[2]) / 128.0f))); LinePS += 3;    LinePD += 3;    LinePM += 3; } } free(Mask); return IM_STATUS_OK; }

　　基本按照论文的公式写的代码，未作优化，测试两张图片看看。

   

　　　　　　　　　　原图1　　　　　　　　　　　　　　　　　　　　　　　　　        Moroney论文的结果

　　彷佛效果还不错。

　　做为一种改进，Contrast image correction method一文做者对上述公式进行了2个方面的调整，以下所示：

　　

　　第一，高斯模糊的mask使用双边滤波来代替，由于双边滤波的保边特性，这样能够减小处理后的halo瑕疵。这没啥好说的。

　　第二，常数2使用变量α代替，而且是和图像内容相关的，具体算式以下：

　　当图像的总体平均值小于128时，使用计算，当平均值大于128时，使用计算，论文做者给出了这样作的理由：对于低对比度的图像，应该须要较强烈的校订，所以α值应该偏大，而对于有较好对比度的图，α值应该偏向于1，从而产生不多的校订量。

　　对于第二条，实际上存在很大的问题，好比对于咱们上面进行测试的原图1，因为他上半部分为天空，下半部分比较暗，且基本各占通常，所以其平均值很是靠近128，所以计算出的α也很是接近1，这样若是按照改进后的算法进行处理，则基本上图像无什么变化，显然这是不符合实际的需求的，所以，我的认为做者这一改进是不合理的，还不如对全部的图像该值都取2，靠mask值来修正对比度。

　　那么对于彩色图像，咱们有两种方法，一种是直接对RGB各份量处理，如上面的代码所示，另一种就是把他转换到YCBCR或者LAB或者YUV等空间，而后只处理亮度通道，最后在转换到RGB空间，那么本文对个人有用的帮助就是提供了一个恢复色彩饱和度的方法。通常来讲在对Y份量作处理后，再转换到RGB空间，图像会出现饱和度必定程度丢失的现象，看上去图像彷佛色彩不足。以下图中间图所示，所以，论文提出了下面的修正公式：

　　

　　经测试，这样处理后的图色彩仍是很鲜艳的，和直接三通道分开处理的差很少（直接三通道分开处理有可能会致使严重偏色，而只处理Y则不会）。

　　    

　　　　　　　　　　　原图　　　　　　　　　　　　　  　直接处理Y通道再转换到RGB空间                                                改进后的效果

　　咱们贴出按照上述思路改进后的代码：

int IM_LocalExponentialCorrection(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
    unsigned char *OldY = NULL, *Mask = NULL, *Table = NULL;
    OldY = (unsigned char *)malloc(Height * Width * sizeof(unsigned char));
    Mask = (unsigned char *)malloc(Height * Width * sizeof(unsigned char));
    IM_GetLuminance(Src, OldY, Width, Height, Stride);            //    获得Y通道的数据
    IM_GuidedFilter(OldY, OldY, Mask, Width, Height, Width, IM_Max(IM_Max(Width, Height) * 0.01, 5), 25, 0.01f);    //    经过Y通道数据处理获得255-Mask值
    unsigned char *NewY = Mask;
    for (int Y = 0; Y < Height * Width; Y++)
    {
        NewY[Y] = IM_ClampToByte(255 * pow(OldY[Y] * IM_INV255, pow(2, (128 - (255 - Mask[Y])) / 128.0f)));
    }

    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        unsigned char *LinePD = Dest + Y * Stride;
        unsigned char *LinePO = OldY + Y * Width;
        unsigned char *LinePN = NewY + Y * Width;
        for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3, LinePO++, LinePN++)
        {
            int Old = LinePO[0], New = LinePN[0];
            if (Old == 0)
            {
                LinePD[0] = 0;    LinePD[1] = 0;    LinePD[2] = 0;
            }
            else
            {
                LinePD[0] = IM_ClampToByte((New * (LinePS[0] + Old) / Old + LinePS[0] - Old) >> 1);
                LinePD[1] = IM_ClampToByte((New * (LinePS[1] + Old) / Old + LinePS[1] - Old) >> 1);
                LinePD[2] = IM_ClampToByte((New * (LinePS[2] + Old) / Old + LinePS[2] - Old) >> 1);
            }
        }
    }
    free(OldY);
    free(Mask);
    return IM_STATUS_OK;
}

　　代码并不复杂，基本就是按照公式一步一步编写的，其中IM_GetLuminance和IM_GuidedFilter为已经使用SSE优化后的算法，对于本文一直使用的测试图675*800大小的图，测试时间大概再40ms，而上述两个SSE的代码耗时才5ms不到，所以，能够进一步优化。

　　第一个须要优化的固然就是那个NewY[Y]的计算过程了，里面的pow函数是很是耗时的，仔细观察算式里只有两个变量，切他们都是[0,255]范围内的，所以创建一个256*256的查找表就能够了，以下所示：

    Table = (unsigned char *)malloc(256 * 256 * sizeof(unsigned char));
    for (int Y = 0; Y < 256; Y++)
    {
        float Gamma = pow(2, (128 - (255 - Y)) / 128.0f);
        for (int X = 0; X < 256; X++)
        {
            Table[Y * 256 + X] = IM_ClampToByte(255 * pow(X * IM_INV255, Gamma));
        }
    }
    
    for (int Y = 0; Y < Height * Width; Y++)
    {
        NewY[Y] = Table[Mask[Y] * 256 + OldY[Y]];
    }
　　 free(Table);

　　速度一会儿跳到了15ms，因为是查表，基本上无SSE优化的发挥地方。

　　接着再看最后的饱和度校订部分的算法，核心代码即：

    LinePD[0] = IM_ClampToByte((New * (LinePS[0] + Old) / Old + LinePS[0] - Old) >> 1);
    LinePD[1] = IM_ClampToByte((New * (LinePS[1] + Old) / Old + LinePS[1] - Old) >> 1);
    LinePD[2] = IM_ClampToByte((New * (LinePS[2] + Old) / Old + LinePS[2] - Old) >> 1);

　　注意到这里是以24位图像为例的，其实24位图像在进行SSE优化时有的时候比32位麻烦不少，由于32位一个像素4个字节，一个SSE变量正好能容纳4个像素，而24位一个像素3个字节，不少时候要在编程时把他补充一个alpha，而后处理玩后在把这个alpha去掉。

　　对于本例，注意到还有特殊性，在处理一个像素时还涉及到对应的Y份量的读取，因此有增长了复杂性。

　　咱们在看上下上面的公式，因为SSE没有整数除法指令，一般状况下要进行整除必须借助浮点版本的除法，所以必须有这种数据类型的转换，另外，咱们考虑把括号里的加法展开下，能够获得公式变为以下：

 LinePD[0] = IM_ClampToByte((New * LinePS[0] / Old + LinePS[0] + New - Old) >> 1);

　　这样展开从C的角度来讲不会产生什么大的性能差别，可是对于SSE编程却有好处，注意到New和LinePS[0] 的最大只都不会超过255，所以二者相乘也在ushort所能表达的范围内，可是若是带上原来的(LinePS[0] + Old) 则会超出ushort范围，对于没有超出USHORT类型的乘法，咱们能够借助_mm_mullo_epi16一次性实现8个数据的乘法，而后在根据须要把他们扩展位32位。

　　具体的优化细节还有不少值得探讨的，因为以前的不少系列文章里基本已经讲到部分优化技巧，所以本文仅仅贴出最后这一块的优化代码，具体细节有兴趣的朋友能够自行去研究：

　　　　　__m128i SrcV = _mm_loadu_epi96((__m128i *)LinePS);
        __m128i OldV = _mm_cvtsi32_si128(*(int *)LinePO);
        __m128i NewV = _mm_cvtsi32_si128(*(int *)LinePN);

        __m128i SrcV08 = _mm_unpacklo_epi8(SrcV, Zero);
        __m128i OldV08 = _mm_shuffle_epi8(OldV, _mm_setr_epi8(0, -1, 0, -1, 0, -1, 1, -1, 1, -1, 1, -1, 2, -1, 2, -1));
        __m128i NewV08 = _mm_shuffle_epi8(NewV, _mm_setr_epi8(0, -1, 0, -1, 0, -1, 1, -1, 1, -1, 1, -1, 2, -1, 2, -1));
        __m128i Temp08 = _mm_sub_epi16(_mm_add_epi16(SrcV08, NewV08), OldV08);
        __m128i Mul08 = _mm_mullo_epi16(SrcV08, NewV08);
        __m128i Value04 = _mm_div_epi32(_mm_unpacklo_epi16(Mul08, Zero), _mm_unpacklo_epi16(OldV08, Zero));
        __m128i Value48 = _mm_div_epi32(_mm_unpackhi_epi16(Mul08, Zero), _mm_unpackhi_epi16(OldV08, Zero));
        __m128i Value08 = _mm_srli_epi16(_mm_add_epi16(_mm_packus_epi32(Value04, Value48), Temp08), 1);

        __m128i SrcV12 = _mm_unpackhi_epi8(SrcV, Zero);
        __m128i OldV12 = _mm_shuffle_epi8(OldV, _mm_setr_epi8(2, -1, 3, -1, 3, -1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1));
        __m128i NewV12 = _mm_shuffle_epi8(NewV, _mm_setr_epi8(2, -1, 3, -1, 3, -1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1));
        __m128i Temp12 = _mm_sub_epi16(_mm_add_epi16(SrcV12, NewV12), OldV12);
        __m128i Mul12 = _mm_mullo_epi16(SrcV12, NewV12);
        __m128i Value12 = _mm_div_epi32(_mm_unpacklo_epi16(Mul12, Zero), _mm_unpacklo_epi16(OldV12, Zero));
        __m128i Value16 = _mm_srli_epi16(_mm_add_epi16(_mm_packus_epi32(Value12, Zero), Temp12), 1);
        _mm_storeu_epi96((__m128i*)LinePD, _mm_packus_epi16(Value08, Value16));

　　这里充分运用的shuffle指令来实现各类需求。

　　优化后速度能够提高到7ms左右。

　   本文最后的运行效果可下载测试：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

　　位于菜单Enhance --> LocalExponentialCorrection下。