這是一篇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下。
