【16位RAW图像处理二】：一种自适应对数映射的高对比度图像显示技术及其速度优化。
办公室今天停电，幸好本本还有电，同事们好多都去打麻将去了，话说麻将这东西玩起来也还是有味的，不过我感觉我是输了不舒服，赢了替输的人不舒服，所以干脆拜别麻坛四五年了，在办公室一个人整理下好久前的一片论文的思想，和万千世界里有缘人共同分享下资源了。
论文的名字是<Adaptive Logarithmic Mapping For Displaying High Contrast Scenes>，相关的PDF文档可以在百度上下载到，翻译成中文的意思是一种显示高对比度场景的自适应对数映射算法，也是一篇很古老的算法文章的，看了下好像是2003年的，在Opencv 3.0中已经提供了这个算法的实现了，但是其实现的细节我觉得写的真恶心（我觉得Opencv所有的算法写的都恶心，饶了一堆弯才到算法的中点，感觉大的工程都是这样的），其算法位于OpenCV3.0\opencv\sources\modules\photo\src\tonemap文件中，我在实现这个算法时时参考了另外一个非常有名的开源软件：luminance hdr， 这个软件最新的版本集成了11中用于HDR显示的算法。过完年我要再次好好的看看这个软件的代码了。
话说回本文的重点，由于对论文的理解不是很深刻，部分内容仅以翻译为主。
一、应用背景
简单的说，就是我们认为显示给人眼看的亮度值Ld和场景亮度值Lw之间存在如下的关系：
其中Lmax为场景的最大亮度，这个映射关系式能够保证无论真实亮度范围有多大，最亮的部分总能映射为1（白色），而其他的亮度也能平缓的变化。虽然这个算式对一些图能获得较为满意的结果，但是我们也发现有些图的亮度压缩太过了，一些高对比度的内容也丢失了。
二、自适应对数映射
那么论文提出的色调映射方案遵循了下面几条规则：（1）不管原始数据如何分散，他必须都能输出连续的结果。（2）它应该具有自适应性和可扩展性，他必须能显示出场景的物理本质同时不得引入对比度反转和光晕。总的亮度也必须忠实于实际的内容。（3）算法也需要对用户友好，也就是说在大部分情况下需要能自动实现，少数情况需要调节一些比较直观的参数。
1、把场景亮度映射到图像亮度
输出图像的整体亮度主要是由场景的亮度特性决定的，所以找到一个从场景亮度到输出图像亮度之间的初始缩放因子很重要，这就类似于在摄影中曝光设置决定了所拍摄的图片的最终效果一样。现在的镜头都提供的不同的自动曝光选项，比如center- weighted, center-spot以及 matrix-metering等。同样，本文提出两种不同的方法适合于不同的用途。对于静态图片或者当用户不直接和场景交互，我们基于所有像素的亮度信息计算出场景的对数平均值作为初始缩放因子，对于需要交互的场合，缩放因子不是固定值，而是使用亮度的对数信息的高斯模糊后的结果，通常高斯模糊核能覆盖场景的15%范围即可以，当然也可以调整这个范围。
2、对比度调整
本文提出的核心的最具特色的色调映射函数就是根据每个像素的信息来自适应的调整函数中的对数基（从2到10）。这从本质上提供了很好的对比度和细节信息，同时能对高亮度值进行很大程度上的压缩。原则上，一个更宽或者更窄的对数基也可以使用，但是实际上，我们没有任何理由去使用它。当对数基小于2时，其值迅速增加，导致曝光调整很困难。另一方面，当对数基大于10时，亮度压缩的量很小，导致这个图片丢失了太多的对比度。同时，我们也观察到了高对数基时的一些颜色偏移现象。
如下左图所示，左图时基于2对数基的，有图是基于10对数基的（全图），很明显，两幅图的对比度和亮度区别很明显，但是他们都没有给出非常令人满意的结果。
为了实现不同像素不同的对数基以及像素的连续性的要求，我们参考了 Perlin和Hoffert的偏置对数函数，该函数是纹理分析的的标准工具并且在计算机视觉上广为应用。偏置函数定义在单位区间[0,1]之间的power函数，有一个参数b,直接决定了输出值的大小，具体形式如下：
其曲线如上右图所示。
2.3 算法细节
一般数据是基于RGB空间的，我们首先将数据转换到XYZ空间，其中的Y分量反应了每个像素的亮度值。我们首先基于Y分量计算出对数的最大亮度Lwmax和平均亮度Lwa。然后把公式（3）带入公式(1)并做适当的阔啊站，得到最终的计算显示亮度的公式如下所示：
相比于论文，上述公式后半部分是我自己添加上去的，主要是为了解释方便。
我们首先看下算式中的
，很明显，他的取值范围是2到10之间，这和论文前面的描述是一致的。然后底部的Log10也保证了整个算式的区间范围。
然后我们知道，对数计算式有如下特性：
因此把公式（4）的后半部分展开就到了论文的结果。
式中有多了个参数Ldmax,这个值表示显示设备的最大显示能力，对于普通的CRT显示器，我们直接取值为100。
论文后面还有关于Gamma校正的内容，那些都是辅助的了，实际上没啥意思，论文核心的就是上述公式。
三：具体实现即细节注意：
具体实现代码可以完美的参考luminance hdr，关键是要注意一些数据的范围要映射到0和1之间才能处理，特别是论文有些地方其实没有讲的特别清晰，比如在论文里有这样一句话：The XYZ luminance component Y of each pixel (Lw for world luminance) and the maximum luminance of the scene Lwmax are divided by the world adaptation luminance Lwa and eventually multiplied by an exposure factor set by the user。这里其实没有明确的说Lwa是什么（整篇论文都没有说）其实就是上面讲的对数平均值，还有最后直接用公式计算得到的Ld一般情况下是很小的，如何处理让其显示也是值得讲究的。我这里贴出对公式（4）计算的核心代码：
void ComputeScaleTable(float ScaleTable[], float Bias, float Saturation, float MaxLum, float AvLum) { float Divider, BiasP, NormalY, Interpol, NewLum; MaxLum = MaxLum / AvLum; // normalize maximum luminance by average luminance（divided by the world adaptation luminance Lwa） Divider = log10f(MaxLum + 1.0f); // 论文公式（4）左半部分的除数 BiasP = logf(Bias) / log(0.5); // 公式（3）中的上指标 for (int Index = 0; Index < 256; Index++) // Normal tone mapping of every pixel { // The XYZ luminance component Y of each pixel (Lw for world luminance) and the maximum luminance of the scene // Lwmax are divided by the world adaptation luminance Lwa and eventually multiplied by an exposure factor set by the user. NormalY = (Index / 255.0) / AvLum; // divided by the world adaptation luminance Lwa Interpol = logf(2.0f + powf(NormalY / MaxLum, BiasP) * 8.0f); // 论文公式（4）右半部分的除数 NewLum = (logf(NormalY + 1.0f) / Interpol) / Divider; // 论文公式（4） ScaleTable[Index] = powf(NewLum / (Index / 255.0 + 1e-4), Saturation); } }
以上代码是针对8位图像的，上面的/255.0就是归一化到[0,1]范围的作用。但是最后一行的NewLum / (Index / 255.0 + 1e-4)你们能理解是什么意思吗?
最后一行代码的Saturation的作用见<Gradient domain high dynamic range compression>一文，当大于1时，图像越饱和也越亮，小于1是图像变暗。
虽然论文描述的算法本意是用到HDR这种高动态范围的图像的，但是实际上我目前也只实现了8位的LDR的代码，但是对于8位的图像，特别是偏暗的图像还是有很好的增强的效果的，对于正常的图像，一般也不会出现特别不好的效果。
对于常态的图片，一般也能起到一定的视觉增强效果：
关于速度优化方面，如果是针对8位图像，则中间的很多浮点计算可以用查表代替，而XYZ和RGB空间转化，我前面一篇博客已经提到可以用SSE快速实现。处理1080P的图大概需要20ms。
也曾尝试不转到XYZ空间，直接提取出亮度信息，然后直接在RGB空间处理，似乎效果也还可以，但是有可能会出现较多的偏色。
明年有时间把这个算法扩展到16位图像上取看看，有什么效果。
8位测试工程：https://files.cnblogs.com/files/Imageshop/TonemapDrago.rar
