SSE圖像算法優化系列二十：一種快速簡單而又有效的低照度圖像恢復算法。

　又有好久沒有動筆了，主要是最近沒研究什麼東西，並且如今主流的趨勢都是研究深度學習去了，但本身沒這方面的需求，同時也就不多有動力再去看傳統算法，今天一我的在家，仍是抽空分享一個簡單的算法吧。

　　前段日子在看水下圖像處理方面的資料時，在github搜到一個連接，裏面竟然有好幾篇文章附帶的代碼，除了水下圖像的文章外，我看到了一篇《Adaptive Local Tone Mapping Based on Retinex for High Dynamic Range Images  》的文章，也還不算老，2013年的，隨便看了下，內容比較簡單，而且做者還分享了Matlab和Java的相關代碼，所以，我也嘗試這用C和SIMD作了優化，感受對於低照度的圖像處理起來效果仍是很不錯，所以，記錄下優化和開發過程的一些收穫。

　　上述github連接爲： https://github.com/IsaacChanghau/OptimizedImageEnhance

　　論文提出的主要內容時這對高動態圖像的顯示問題，他結合傳統的Retinex技術提出了全局自適應和局部自適應的HDR實現過程，我也實現了整個的代碼，但感受前面的全局適應方案特別對於低照度圖像有着很是明顯的調節做用，所以，我重點談下整個。

　　直接應用原文的英文算了：

　Global adaptation takes place like an early stage of the human visual system [4]. The human visual system senses rightness as an approximate logarithmic function according o the Weber-Fechner law [5]. To globally compress the ynamic range of a HDR scene, we use the following function n (4) presented in [5].

       

　　用中文解釋下上面的公式，也是本文最重要的一個公式。

　 是全自適應輸出的結果，咱們這裏就是須要獲得他，表示輸入圖像的luminance值（亮度值），表示輸入圖像亮度值對的最大值，表示輸入亮度對數的平均值，以下式所示：

　　

　　其中N表示像素的總數，而δ通常是個很小的值，其做用主要是爲了不對純黑色像素進行log計算時數值溢出，這個問題在圖像處理時很是常見。

　　在log域進行計算，這個HDR算法中基本是個定律了。 

　　直接應用原文的話，上述算式的主要做用是：

      The input world luminance values and the maximum luminance values are divided by the log-average luminance of he scene. This enables (4) to adapt to each scene. As the log-verage luminance converges to the high value, the function converges from the shape of the logarithm function to the near function. Thus, scenes of the low log-average luminance reboosted more than scenes with high values. As a result, the overall scene luminance values are adequately compressed in ccordance with the log-average luminance of the scene.

 　　特別注意的是 scenes of the low log-average luminance reboosted more than scenes with high values. 這句話，他的意思是說低照度的亮度部分比高照度的部分要能獲得更大程度的提高，因此對於低照度圖，上述公式能起到很好的加強做用。而算式中使用了全局的對數平均值，這就有了必定的自適應性。

　　我貼一段稍微修改了的做者共享的matlab代碼做爲本算法的參考代碼：

function outval = ALTM_Retinex(I)
II = im2double(I);
Ir=double(II(:,:,1)); Ig=double(II(:,:,2)); Ib=double(II(:,:,3));
% Global Adaptation
Lw = 0.299 * Ir + 0.587 * Ig + 0.114 * Ib;% input world luminance values
Lwmax = max(max(Lw));% the maximum luminance value
[m, n] = size(Lw);
Lwaver = exp(sum(sum(log(0.001 + Lw))) / (m * n));% log-average luminance
Lg = log(Lw / Lwaver + 1) / log(Lwmax / Lwaver + 1);
gain = Lg ./ Lw;
gain(find(Lw == 0)) = 0;
outval = cat(3, gain .* Ir, gain .* Ig, gain .* Ib);
figure;
imshow(outval)

　　很簡單的代碼，貼一個這個代碼的結果:

          

　 網上K的一個常常用來測試的美女，也不知道是那位性福男士的女友，左圖很是的暗淡無光澤，處理後的圖像飽和度和亮度都有較大的提高。

　 M代碼一般只是用來學習用的，並不具備工程價值，M代碼也基本不考慮速度和內存佔用，所以通常都要根據M代碼自行修改成C或者C++代碼纔有可能實用。

　我貼出部分我寫的C代碼來分析下提速的辦法。

int IM_ALTM_Retinex(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
	float Avg = 0;
	unsigned char *OldY = (unsigned char *)malloc(Height * Width * sizeof(unsigned char));
	unsigned char *NewY = (unsigned char *)malloc(Height * Width * sizeof(unsigned char));
	float *LogG = (float *)malloc(Height * Width * sizeof(float));
	float *Table = (float *)malloc(256 * sizeof(float));
	IM_GetLuminance(Src, OldY, Width, Height, Stride);						//	input world luminance values

	//	 *************************************************   Global Adaptation     ******************************************
	int MaxL = IM_GetMaxValue(Src, Width, Height, Width);					//	the maximum luminance value

	for (int Y = 0; Y < 256; Y++)				Table[Y] = log(0.000001F + Y * IM_INV255);
	for (int Y = 0; Y < Height * Width; Y++)	Avg += Table[OldY[Y]];		//	sum of log luminance value

	Avg = exp(Avg / (Height * Width));										//	log - average luminance

	for (int Y = 0; Y < 256; Y++)				Table[Y] = log(Y* IM_INV255 / Avg + 1) / log(MaxL * IM_INV255 / Avg + 1);
	for (int Y = 0; Y < Height * Width; Y++)	LogG[Y] = Table[OldY[Y]];	//	globally compress the dynamic range of a HDR scene we use the following function in(4) presented in[5].

	IM_Normalize(LogG, NewY, Width, Height);								//	after normalization, an output image is obtained from the processed luminance values and the input chrominance values

	IM_ModifyYComponent(Src, NewY, Dest, Width, Height, Stride);

	if (OldY != NULL)		free(OldY);
	if (LogG != NULL)		free(LogG);
	if (NewY != NULL)		free(NewY);

	return 0;
}

　　第一步，獲得input world luminance，使用IM_GetLuminance函數，具體的實現細節能夠參考SSE圖像算法優化系列十五：YUV/XYZ和RGB空間相互轉化的極速實現(此後老闆不用再擔憂算法轉到其餘空間通道的耗時了)一文自行實現。

　　第二步：獲得the maximum luminance value，使用IM_GetMaxValue函數，我分享下這個函數的SSE實現代碼。

//	求16個字節數據的最大值, 只能針對字節數據。
inline int _mm_hmax_epu8(__m128i a)
{
	__m128i b = _mm_subs_epu8(_mm_set1_epi8(255), a);
	__m128i L = _mm_unpacklo_epi8(b, _mm_setzero_si128());
	__m128i H = _mm_unpackhi_epi8(b, _mm_setzero_si128());
	return 255 - _mm_extract_epi16(_mm_min_epu16(_mm_minpos_epu16(L), _mm_minpos_epu16(H)), 0);
}

int IM_GetMaxValue(unsigned char *Src, int Length)
{
	int BlockSize = 16, Block = Length / BlockSize;
	__m128i MaxValue = _mm_setzero_si128();
	for (int Y = 0; Y < Block * BlockSize; Y += BlockSize)
	{
		MaxValue = _mm_max_epu8(MaxValue, _mm_loadu_si128((__m128i *)(Src + Y)));
	}
	int Max = _mm_hmax_epu8(MaxValue);
	for (int Y = Block * BlockSize; Y < Length; Y++)
	{
		if (Max < Src[Y])	Max = Src[Y];
	}
	return Max;
}

int IM_GetMaxValue(unsigned char *Src, int Width, int Height, int Stride)
{
	int Max = 0;
	for (int Y = 0; Y < Height; Y++)
	{
		Max = IM_Max(Max, IM_GetMaxValue(Src + Y * Stride, Width));
		if (Max == 255) break;
	}
	return Max;
}

 　　用到了函數的多態特性，其中求最大值的思路在 一種可實時處理 O(1)複雜度圖像去霧算法的實現一文中有較爲詳細的說明。

　　 單行像素裏求取最大值，能夠充分利用_mm_max_epu8這條SIMD指令，一次性進行16次比較，速度大爲提升，而最後從SIMD寄存器裏求出16個字節數據的最大值則充分利用了_mm_minpos_epu16這個SIMD指令，使人感到困惑的是爲何系統只提供_mm_minpos_epu16指令，而沒有_mm_minpos_epu8或者_mm_minpos_epu32這樣的，16有什麼特殊的場合用的最廣呢。

　　求對數值這對於8位圖像，最快速的方式也只有查表，而且查表能夠明確的說是沒有好的SIMD加速方法，除非是16個字節的小表（返回值也是字節的），能夠利用_mm_shuffle_epi8加速外，其餘的都不行，而這個的應用場合極爲少見。

　　查完表後計算出對數的平均值Avg，最後按照公式（4）獲得全局的自適應輸出值。

    是浮點數，咱們須要將其轉換爲圖像數據，這裏有個額外的選項，一種是在轉換前對浮點數仍是作點處理，不少M代碼裏都有這個過程，一般叫作SimplestColorBalance.m（我看到過不少次了），這個其實就有點相似於圖像處理的自動色階過程，把必定百分比的低亮度和高亮度值刪除掉，而後中間的值進行拉昇。以後的normalization就是正常的線性處理了，這個浮點處理過程也可使用SIMD指令加速。

//	將浮點數按照最大值和最小值線性的插值到0和255之間的像素值。
int IM_Normalize(float *Src, unsigned char*Dest, int Width, int Height)
{
	if ((Src == NULL) || (Dest == NULL))					return IM_STATUS_NULLREFRENCE;
	if ((Width <= 0) || (Height <= 0))						return IM_STATUS_INVALIDPARAMETER;
	
	float Min = 0, Max = 0;
	
	IM_GetMinMaxValue(Src, Width * Height, Min, Max);

	if (Max == Min)
	{
		memset(Dest, 128, Width * Height);
	}
	else
	{
		float Inv = 255.0f / (Max - Min);							//	不建議用整數的乘以255，由於可能會溢出，反正最後的除法要調用浮點版本的，這裏就用浮點乘不是更好嗎
		int BlockSize = 8, Block = (Height * Width) / BlockSize;
		for (int X = 0; X < Block * BlockSize; X += BlockSize)		//	正規化
		{
			__m128 SrcV1 = _mm_load_ps(Src + X);
			__m128 SrcV2 = _mm_load_ps(Src + X + 4);
			__m128 Diff1 = _mm_sub_ps(SrcV1, _mm_set1_ps(Min));
			__m128 Diff2 = _mm_sub_ps(SrcV2, _mm_set1_ps(Min));
			__m128i Value1 = _mm_cvtps_epi32(_mm_mul_ps(Diff1, _mm_set1_ps(Inv)));
			__m128i Value2 = _mm_cvtps_epi32(_mm_mul_ps(Diff2, _mm_set1_ps(Inv)));
			_mm_storel_epi64((__m128i *)(Dest + X), _mm_packus_epi16(_mm_packus_epi32(Value1, Value2), _mm_setzero_si128()));
		}
		for (int X = Block * BlockSize; X < Height * Width; X++)
		{
			Dest[X] = (int)((Src[X] - Min) * Inv + 0.5f);
		}
	}
	return IM_STATUS_OK;
}

　　就是簡單的SIMD指令運用，沒啥好說的。

　　因爲以前處理獲得值仍是luminance值，若是是灰度圖，處理到這裏就能夠結束了，可是若是是彩色RGB圖，咱們有不少方案來得到最終的RGB份量值，這裏我提三種方式。

　　第一種，如上述C++代碼所示，使用IM_ModifyYComponent這樣的方式，細節上爲把原圖轉換到YUV或者HSV這中帶亮度的顏色空間中，而後用新獲得的luminance值代替Y通道或V通道，而後在轉換會RGB空間。

　　第二種方法，如上述Matalb代碼所示，用新的luminance值和原始luminance值的比值做爲三通到的加強係數，這樣三通道能夠獲得一樣程度的加強。

　　第三種方法就是在SSE圖像算法優化系列十九：一種局部Gamma校訂對比度加強算法及其SSE優化一文中提到的 算式。

　　第一種方法容易出現結果圖色彩偏淡，第二種每一個份量易出現過飽和，第三種可能要稍微好一點，建議使用第三種。

　　但總的來講差別可能都不會太大。

　    貼一些改過程的圖片看看效果，你們也能夠本身用下面的圖片作測試，看看結果如何。

  

  

   

 

   

 

  

   

　　倒數第二張是一幅灰度圖，可見其加強效果是很是明顯的，而最後一張是一副本來就不錯的彩色照片，通過算法處理後總體變化不大，稍微更加鮮豔一點，這是有好處的，說明改算法對本來就不錯的圖，處理後質量不會有明顯的差別，這正是咱們所但願的。

　　倒數第三張圖的效果和matlab處理的有這必定的卻別，特別是用方框框注的那一塊區域，能夠看到M代碼出現了明顯的紅色色塊，這主要是因爲M代碼使用的各份量乘以同一個係數致使的溢出形成的，而使用第三種方法澤能有效地避免該問題。

　　那麼在原文中，做者還提出局部的自適應處理，主要也是基於報邊濾波器和Log域進行了處理，你們能夠直接運行做者提供的matlab代碼試一試，可是那個代碼彷佛對原文作了很多的添加，特備是係數計算方面，不知道他的依據是什麼。

　　論文中還說到，全局處理容易出現halo現象，可是在咱們的測試圖中均未出現，其實這主要是因爲我這些圖都是LDR圖，像素取值範圍有限，即便對比度低，也不會出現HDR數據中數值之間可能出現的巨大差別，所以，對於高動態圖像，後續的局部處理可能尤其必要，下半年個人部分時間可能會在這方面作點研究。

　　速度測試：1080P的圖像處理時間約爲20ms。

　　算法比較簡單，有興趣的朋友自行編程C代碼，應該不難實現，測試DEMO見下述附件的Enhance-->全局低照度加強菜單。

　　Demo下載地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar