SSE图像算法优化系列二十：一种快速简单而又有效的低照度图像恢复算法。

　又有很久没有动笔了，主要是最近没研究什么东西，而且现在主流的趋势都是研究深度学习去了，但自己没这方面的需求，同时也就很少有动力再去看传统算法，今天一个人在家，还是抽空分享一个简单的算法吧。

　　前段日子在看水下图像处理方面的资料时，在github搜到一个链接，里面居然有好几篇文章附带的代码，除了水下图像的文章外，我看到了一篇《Adaptive Local Tone Mapping Based on Retinex for High Dynamic Range Images  》的文章，也还不算老，2013年的，随便看了下，内容比较简单，并且作者还分享了Matlab和Java的相关代码，因此，我也尝试这用C和SIMD做了优化，感觉对于低照度的图像处理起来效果还是很不错，因此，记录下优化和开发过程的一些收获。

　　上述github链接为：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