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