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));