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下。