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