Objective Quality Assessment of Tone-Mapped Images

Abstract

将高动态范围(HDR)图像转换为低动态范围(LDR)图像的色调映射操作符(TMOs)为HDR图像在标准LDR显示上的可视化提供了实用的工具。不同的TMOs创建不同的色调映射图像，一个自然的问题是哪一个有最好的质量。没有适当的质量衡量，不同的tmo无法进行比较，进一步的改进也没有方向。主观评分可能是一种可靠的评价方法，但它昂贵且耗时，更重要的是，难以嵌入到优化框架中。本文将基于结构相似度的多尺度信号保真度测度和基于自然图像强度统计的自然度测度相结合，提出了一种客观的图像质量评价算法。使用独立的受试者评分图像数据库进行的验证表明，主观评分与提出的色调映射图像质量指数(TMQI)之间存在良好的相关性。此外，我们通过两个例子来展示TMQI的扩展应用:tmo的参数调整和多色调映射图像的自适应融合。

I. INTRODUCTION

近年来，人们对高动态范围(HDR)图像越来越感兴趣，HDR图像的强度水平范围可能在10,000到1[1]之间。这使得真实场景中亮度变化的准确表征得以实现，范围从直射阳光到微弱的星光[1]。随着成像和计算机图形技术的最新进展，HDR图像正变得越来越广泛。在实践中经常遇到的一个常见问题是，如何在标准显示设备上可视化HDR图像，这些设备被设计成显示低动态范围(LDR)图像。为了克服这个问题，越来越多的色调映射算子(TMOs)被开发出来，可以将HDR图像转换为LDR图像，例如[2]-[5]。由于动态范围的缩小，色调映射过程不可避免地会造成信息丢失。如果有多个TMO可用，我们会问哪个TMO忠实地保留了原始HDR图像中的结构信息，哪个TMO生成了看起来自然逼真的LDR图像。

过去的TMO评估大多依赖于人的主观评价。在[6]中，对6个tmo进行了关于相似性和偏好的知觉评价。在[7]报告了8个tmo的综述和主观比较。在主观实验中采用配对比较法对[8]中6个tmo进行HDR监测比较。在[9]中，14名受试者被要求根据基本图像属性和LDR图像的自然度对7个TMOs生成的2个建筑室内场景进行评分。在[10]中进行了更全面的主观评价，将14个TMOs生成的色调映射图像展示给2组10人的观察者，对LDR图像的整体质量、亮度、对比度、细节再现和颜色进行评分。在[11]中，受试者被要求从2个不同参数设置的TMOs中选择最佳LDRs来优化算法。主观测试的价值怎么估计也不过分。然而，它们有根本的局限性。首先，这是昂贵和耗时的。其次，很难整合到优化框架中，自动改进TMOs并调整其参数设置。此外，HDR图像中包含的重要图像结构可能在色调映射图像中缺失，但人类观察者可能不知道它们的存在。从这个意义上说，主观评价不应被视为色调映射图像质量的金标准。

典型的客观图像质量评估(IQA)方法假设参考图像和测试图像具有相同的动态范围[12]，因此不能直接应用于评价色调映射图像。目前对HDR图像的客观评价方法较少。HDR可见差异预测器(HDR- vdp)[1]，[13]是一种基于人类视觉系统(HVS)的保真度度量，旨在区分可见(阈上)和不可见(阈下)失真。度量反映了在检测概率方面的失真感知。由于HDR- vdp的设计目的是预测相同动态范围内两幅HDR图像之间差异的可见性，因此不适合将HDR图像与LDR图像进行比较。在[14]中提出了一种动态范围无关的方法，该方法改进了HDR-VDP，并产生了三种类型的质量图，分别表明可见特征的丢失、不可见特征的放大和对比极性的逆转。这些质量图显示了良好的相关性与主观分类的图像退化类型，包括模糊，锐化，对比度逆转，和无失真。但是，它并没有为整个图像提供单一的质量评分，因此无法对图像的整体质量进行主观评价。

当前工作的目的是开发一个客观的IQA模型的色调映射LDR图像，使用其对应的HDR图像作为参考。我们的工作受到了IQA文献中两个成功设计原则的启发。首先是结构相似度(SSIM)方法[15]及其多尺度衍生方法[16]，[17]，认为视觉的主要目的是从视觉场景中提取结构信息，因此结构保真度可以很好地预测感知质量。第二种是自然场景统计(NSS)方法，该方法认为视觉系统高度适应自然视觉环境，并使用与自然图像统计的背离作为感知质量[18]的度量。在这里，我们提出了一种方法，结合了多尺度结构保真度测量和统计自然度测量，导致色调映射图像质量指数(TMQI)。此外，我们还证明了TMQI可以用于TMOs的参数优化和多色调映射图像的自适应融合。

II. QUALITY ASSESSMENT METHOD

由于动态范围的缩小，TMOs无法保存HDR图像中的所有信息，这些图像的LDR版本的人类观察者可能没有意识到这一点。因此，结构保真度在评价色调映射图像[19]的质量中起着重要的作用。另一方面，结构保真度本身并不足以提供一个全面的质量评估。一个好的色调映射图像应该在结构保真度和统计自然度之间取得良好的折衷，这有时是相互竞争的因素。

A. Structural Fidelity

SSIM方法为测量图像之间的结构保真度提供了一种实用的设计理念和方法。原始的SSIM算法是局部应用的，包含三个比较分量:亮度、对比度和结构。由于TMOs的目的是改变局部强度和对比度，直接比较局部和对比度是不合适的。设x和y分别是从HDR和tone-mapped LDR图像中提取的两个局部图像patch。我们将局部结构保真度定义为

其中，σx、σy和σxy分别为HDR和LDR图像中两个对应patch的局部标准差和互相相关，C1和C2为正稳定常数。与SSIM定义[15]相比，缺少了亮度比较分量，(1)中的第二项结构比较分量完全相同。(1)中的第一项比较了信号强度，并基于两个直观考虑对SSIM定义的第一项进行了修改。首先，当HDR和LDR图像块的信号强度都显著(高于可见性阈值)或都不显著(低于可见性阈值)时，不应惩罚HDR和LDR图像块的信号强度差异。其次，算法应该惩罚信号强度在其中一个图像块显著而在另一个图像块不显著的情况。这与原始SSIM定义中的相应术语不同，在原始SSIM定义中，信号强度的任何变化都会受到惩罚。

为了区分信号强度的显著性和不显著性，我们将局部标准差σ通过非线性映射，从而得到(1)中所采用的σ'值。非线性映射应设计为显著性信号强度映射为1，不显著性信号强度映射为0。中间有一个平稳的过渡。因此，非线性映射与对比度的视觉敏感性有关，这是视觉心理物理学[21]文献中广泛研究的课题。实际上，HVS并没有一个固定的对比度检测阈值，但在观察对比度变化时通常遵循一个逐渐增加的概率。描述信号强度检测概率的心理测量函数已被用于从心理物理实验中获得的数据建模。一般来说，心理测量函数类似于s形[22]，[23]，感觉阈值通常定义在50%的检测概率水平上。一个普遍采用的心理测量函数被称为Galton 's ogive[21]，它采用累积正态分布函数的形式给出

其中p为检测概率密度，s为正弦刺激的振幅，τs为调制阈值，θs为控制检测概率变化斜率的正态分布的标准差。经研究发现，比例

大概是一个常数，被称为克罗泽定律[21][24]。k的典型值在2.3 ~ 4之间，k = 3使得误报概率非常小[21]。

视觉对比灵敏度通常用调制阈值的倒数τs来量化，后者是空间频率的函数，即对比灵敏度函数(CSF)[21]。[25]给出了一个与各种心理实验数据相吻合的CSF公式

其中f为空间频率。这个函数被归一化为峰值为1，因此只提供整个频谱的相对灵敏度。在实践中，它需要被缩放到λ常数，以适应心理数据。在我们的实现中，我们遵循Kelly的CSF测量[26]。将此与(4)相结合，得到

这个阈值是基于单纯正弦刺激下的对比敏感度测量计算出来的。为了将其转换为用信号的标准差测量的信号强度阈值，我们需要考虑到信号的幅度与对比度和平均信号强度同时缩放，正弦信号的幅度和标准差之间有一个√2的因子。因此，信号标准差上定义的阈值可以计算为

局部结构保真度测度Slocal应用于图像，使用滑动窗口，在图像空间中运行。这就形成了一幅反映结构保真度在空间上变化的地图。图像细节的可见性取决于图像的采样密度、图像与观察者之间的距离、显示的分辨率以及观察者视觉系统的感知能力。单尺度方法无法捕捉这种变化。遵循多尺度[16]和信息加权SSIM[17]的思路，我们采用多尺度方法，对图像进行迭代低通滤波和下采样，形成图像金字塔结构[27]，如图1所示。在每个比例尺上生成局部结构保真度图。图2显示了从两个不同的tmo创建的LDR图像在多个比例尺上计算的这类地图的两个例子。观察这些保真度地图并研究它们与感知图像保真度之间的关系是很有趣的。例如，在图(b)中，最亮的窗口区域的结构细节缺失，但在图(a)中更明显。再例如，在图(a)中，右上方的黑暗区域有不容易识别的详细结构，但在图(b)中更直观。所有这些都在结构保真度图中清晰地反映出来。

在每个比例尺上，地图通过平均得出一个分数:

其中xi和yi分别是被比较的HDR和LDR图像中的第i个patch, Nl是第l个尺度上的patch个数。在文献中，先进的池化策略如基于信息内容的池化[17]已经被证明可以提高IQA算法的性能。然而，在我们当前的实验中，这些先进的池化方法在提出的结构保真度度量中并没有带来显著的性能增益。采用[16]中的方法结合尺度级结构保真度得分计算整体结构保真度

其中L为总量表数，βl为分配给第L个量表的权重。

在我们的结构保真度模型的实现中有几个参数。首先，在计算Slocal时，我们设置C1 = 0.01和C2 = 10，我们发现结构保真度模型的整体性能对这些参数不敏感，在一个数量级内。其次，为了创建每个比例尺上的保真度地图，我们采用与SSIM算法[15]中相同的设置，即使用大小为11×11，标准差为1.5的高斯滑动窗口。第三，与[16]一样，我们假设观察距离为32圈/度，在不混叠的情况下可以表示高达16圈/度分辨率的信号，因此，当将(4)中的CSF应用于最精细尺度测量时，我们使用16圈/度作为空间频率参数。应用于随后更精细尺度的空间频率参数分别为8、4、2、1个周期/度。第四，取(6)中的平均强度值为LDR图像动态范围的均值，即μ = 128。第五，结合跨量表测量，我们设定L = 5， {βl}={0.0448, 0.2856, 0.3001, 0.2363, 0.1333}，遵循[16]中报道的心理物理实验结果。最后，为了评估彩色图像的质量，我们首先将其从RGB颜色空间转换到Yxy空间，然后只在Y分量上应用提出的结构保真度度量。

B. Statistical Naturalness

一个高质量的色调映射LDR图像不仅应该忠实地保持HDR图像的结构保真度，而且看起来也很自然。然而，自然性是一个难以定量定义的主观量。自然图像的统计对图像处理应用和理解生物视觉[28]都有重要意义。在[29]中对色调映射图像的自然度进行了一项有趣的研究，该研究提供了图像自然度与不同图像属性(如亮度、对比度、颜色再现、可见性和细节再现)之间的相关性的有用信息。结果表明，在所有被测属性中，亮度和对比度与感知自然度的相关性较大。这促使我们基于这两个属性构建统计自然性模型。这种选择在定义统计图像自然度的一般概念时可能过于简化(并且可能无法推广到其他使用自然度概念的图像处理应用程序)，但它在我们的模型的简单性和捕捉与我们试图解决的色调映射评估问题(亮度映射是所有色调映射操作中不可避免的问题)相关的最重要的自然性成分的能力之间提供了一个理想的折衷方案。它也是对第II-A节中描述的结构保真度测量的最好补充，其中亮度建模和评估是缺失的。

我们的统计自然度模型是建立在对来自[30]、[31]的约3000张8位/像素灰度图像进行统计的基础上的，这些图像代表了许多不同类型的自然场景。图3显示了这些图像的均值和标准差的直方图，它们是反映图像整体强度和对比度的有用度量。我们发现这些直方图可以很好地拟合使用高斯和贝塔概率密度函数

最近的研究表明，亮度和对比度在自然图像统计和生物计算方面都是很大程度上独立的量。因此，它们的联合概率密度函数将是两者的乘积。因此，我们将统计自然度度量定义为

C. Quality Assessment Model

Section II-A中引入的结构保真度度量S和Section II-B中描述的统计自然度度量N表征了色调映射图像质量的不同方面。它们可以单独使用，也可以联合使用作为向量值测度。然而，在许多实际应用程序中，用户更喜欢用一个分数来表示图像的整体质量。因此，这些参数应该以某种方式组合。在IQA的文献中，已经有早期的工作将图像统计和结构和对比[33]的测量结合起来，尽管是在不同的背景下。在这里，我们定义了一个三参数函数来对联合度量进行标量化，从而得到色调映射图像质量指数(TMQI)

其中0≤a≤1调整两组分的相对重要性，α和β分别决定它们的灵敏度。因为S和N的上界都是1，所以整体质量度量也是1的上界。

(14)中的参数有待确定。在我们的实现中，它们经过了优化，以最适合[34]作者提供的主观评估数据。在他们的实验中，受试者被要求同时观看由两个不同TMOs应用在同一张HDR图像上生成的两张LDR图像，然后选择整体质量更好的那一张。两项研究已经完成，涉及两组受试者。第一项研究是在浙江大学进行的，59名天真的志愿者被邀请做配对比较任务，并填写偏好矩阵。第二项研究是使用亚马逊土耳其机器人(Amazon Mechanical Turk)进行的，这是一项主观评价的在线服务。每对比较分配给150名匿名受试者。数据库包含6个数据集，每个数据集包含5个知名TMOs生成的图像，由Drago et. al. [4]， Durand & Dorsey [35]， Fattal et. al. [5]， Reinhard et. al. [2]， Mertens et. al.[36]引入。然后使用偏好矩阵计算每个文件夹中的主观排名分数。

利用主观数据找到(14)中的最佳参数本质上是一个回归问题。与传统回归问题的主要区别是，这里我们只提供了图像之间的相对排序数据，而不是与单个图像相关的质量分数。我们开发了一种学习方法，其中参数是学习从迭代方法。在每次迭代中，从一个随机选择的数据集中随机选择一对图像。如果模型生成的客观得分与主观排序的顺序相同，则模型参数不发生变化;否则，每个参数将向修正模型误差的方向更新一小步。迭代继续进行，直到收敛。在实验中，我们发现这种迭代学习过程具有很好的收敛性。此外，为了保证方法的稳健性，我们进行了留一交叉验证过程，将数据库(6个数据集)分为5个训练集和1个测试集，重复相同的过程6次，每一次都有不同的训练集和测试集划分。虽然每次都得到一组不同的参数，但它们彼此相当接近，并且对所有训练集和测试集产生相同的排序顺序。最后选取a = 0.8012， α = 0.3046， β = 0.7088作为最终的模型参数。

III. VALIDATION OF QUALITY ASSESSMENT METHOD

验证过程是通过比较我们的客观质量评估结果与主观数据进行的。采用两个评价指标，分别给出如下。

1. 斯皮尔曼秩相关系数(SRCC)定义为

式中di为第I个图像在主观评价和客观评价中的排名差。SRCC是一种基于非参数秩序的相关度量，独立于主观和客观评分之间的任何单调非线性映射。

2) Kendall’s rank-order correlation coefficient (KRCC)是另一个非参数秩相关度量，计算为

式中，Nc和Nd分别为数据集中一致(秩序一致)和不一致(秩序不一致)对的个数。

拟议的TMQI是正在测试的唯一客观质量测量。据我们所知，几乎没有其他方法被提出来比较不同动态范围的图像。唯一的例外是[14]中提出的方法，它创建概率映射来区分可见(阈值以上)和不可见(阈值以下)降级。在分类图像畸变类型时，概率图被证明是有用的，但并不意味着要汇集在一起来产生色调映射图像的整体质量分数。因此，与所提出的方法进行直接比较是不可能的。

在我们的验证过程中进行了三个实验，每个实验都使用不同的学科排名数据库。第一个数据库来自[34]，在第II-C节讨论的参数训练步骤中也使用了[34]。我们在第II-C节中描述的留一交叉验证方法为六个测试数据集中的每一个创建SRCC和KRCC值，其中对于每个数据集，参数是使用其他五个数据集训练的。表1分别显示了主观排名和我们模型预测的KRCC和SRCC值的均值和标准差。

在第二个实验中，我们使用了[10]，[37]中引入的数据库，其中我们使用了10个朴素受试者从同一HDR图像中创建的14张色调映射图像的整体质量排名数据。图像主观排名与我们的结构保真度、统计自然度和整体质量分数之间的KRCC和SRCC值见表二，在这里我们观察到，结构保真度测量单独可以提供主观排名的合理预测。统计自然度测量本身并不能很好地预测整体质量排名，但它补充了结构保真度测量。当这两种度量相结合时，可以更好地预测整体图像质量。值得注意的是，这里的测试数据在训练过程中没有使用，但得到的KRCC和SRCC值与使用第一个数据库进行训练时得到的测试值是相当的。这说明在第II-C节中描述的训练方法具有良好的泛化能力。

第三个实验是使用我们自己开发的数据库进行的。为20名受试者提供了15组色调映射图像，每组图像包括8张由8个TMOs从同一HDR图像生成的图像。由Reinhard et al. [2]， Drago et al. [4]， Durand & Dorsey [35]， Mantiuk et al.[38]和Pattanaik et al.[39]开发的五个tmo的结果使用公开可用的软件Qtpfsgui[40]计算。此外，使用Adobe Photoshop中的内置TMOs创建了另外三张图像，分别是“曝光和伽马”、“均衡直方图”和“局部适应”。所有8个tmo中使用的参数都被设置为它们的默认值，没有进行优化。选用参考HDR图像代表不同的室内和室外场景，均可在线获得[10]，[41]-[43]。在主观测试中，20名观察者被要求对每个图像集中的8幅图像进行从最好到最差的排序。然后对每个图像的主观排名进行平均，得到其在集合中的平均排名分数。

为了评估TMQI方法，我们计算了每个图像集的平均排名分数和客观质量度量之间的KRCC和SRCC值。结果见表三。为了提供一个评价TMQI表现的锚点，我们将其与普通受试者的行为进行了比较。为此，我们首先计算平均排名分数和每个个体受试者对每个图像集给出的排名分数之间的KRCC和SRCC值。然后我们计算这些KRCC和SRCC值对受试者的平均值和标准差，如表三所示。在最后一行给出了所有15幅图像集合的平均KRCC和SRCC值。可以看出，对于所有的图像集，TMQI的KRCC和SRCC值与所有受试者的平均值KRCC和SRCC值都在±1的标准差范围内。这表明TMQI的行为与一般受试者非常相似。

由于TMQI算法不涉及任何昂贵的搜索或迭代过程，因此具有较高的计算效率。我们的未经优化的MATLAB实现在Intel四核2.67 GHz计算机上平均需要0.75秒和2.7秒来评估大小分别为512×512和1024×1024的图像。图4显示了20次HDR-LDR比较的运行时与图像像素数之间的散点图。结果表明，TMQI算法的计算复杂度与图像中的像素数近似成线性关系。相对较低的计算成本使其很容易适应实际应用，涉及迭代优化过程。

IV. APPLICATIONS OF QUALITY ASSESSMENT METHOD

客观的IQA测量方法的应用范围超出了对图像的评价和算法的比较。更广泛的应用扩展到开发新的图像处理算法优化的新型IQA措施。在本节中，我们使用两个示例来演示TMQI的潜力。

A. Parameter Tuning in TMO Algorithm

许多TMOs包含一个或多个参数，其最佳值通常与图像有关。在没有人为干扰的情况下，选择这些参数往往是一项具有挑战性的任务，这可能会导致截然不同的结果。客观的质量度量为自动选择这些参数提供了一个有用的工具。这里我们以[4]中提出的TMO为例，它在不同位置使用不同基底的对数函数自适应地改变动态范围。给出了算法

其中Lw和Lwmax分别为场景的世界亮度和最大亮度，Ld和Ldmax分别为显示亮度和显示的最大亮度，b为调优参数。色调映射图像的感知质量随b的变化显著，但文献中，通过对多幅图像[4]，[40]的实证实验，通常将b值固定在0.8左右。

在图5(a)和图5(b)中，我们分别绘制了“Desk”和“Bristol Bridge”图像的TMQI随b的函数变化情况(由于b = 1超出了算法的建议值范围，因此不进行超过b = 1的计算)。从图中可以看出，质量分数作为b的函数时，表现出了很大的不同。根据图的结果，我们分别选取了b = 0.8和b = 1作为这两幅图像的最佳值。这些结果证实了最优的b值接近以往研究中选取的经验值(0.8左右)，但不同的图像会有所不同。所选的三个b值对应的色调映射LDR图像分别如图6和图7所示。仔细检查这些图像表明，最好的b值导致了保存结构细节和产生自然图像之间的良好平衡。

B. Adaptive Fusion of Tone-Mapped Images

当在不同的HDR图像上使用不同的TMO进行实验时，我们经常发现很难挑选出一个对所有HDR图像产生最佳效果的TMO。此外，在单个HDR图像中，当考虑图像中的不同区域时，最佳TMO也可能会有所不同。为了充分利用多TMOs的优势，可以采用图像融合技术对多幅色调映射图像进行融合，而客观质量度量在这一过程中起着重要作用。

给定由不同TMOs创建的多色调映射图像，我们首先应用拉普拉斯金字塔变换，将这些图像分解成不同的尺度。在金字塔域，这导致了相同尺度和相同空间位置的多个系数，每个系数对应不同的TMO。图8的前两行给出了例子，展示了四尺度的拉普拉斯金字塔分解，其中精细尺度系数(尺度1-3)代表图像细节，而最粗尺度系数(尺度4)保持了空间上的局部平均强度。然后利用融合策略将多个系数在每个尺度上的每个位置合并成一个系数，然后利用拉普拉斯金字塔逆变换重建融合图像。典型的融合方案是局部选取最显著的图像特征[44]。最广泛采用的方法包括取系数的平均值或选择绝对值最大的系数之一。

在此，我们提出了一种不同的融合方案。总体思路是利用TMQI作为融合过程中的权重因子。设Sj和c j分别是局部结构保真度和从被融合的第j个色调映射图像计算的拉普拉斯金字塔变换系数。熔融系数计算为

这适用于除最粗糙的量表外的所有量表，我们使用统计自然度测度作为权重因子:

其中Nj为第j个色调映射图像的统计自然度得分。

本文提出的拉普拉斯金字塔域融合方法如图8的下一行所示，融合后的图像保留了(f)中最亮区域(顶部的亮区)的细节，同时在较暗的区域保持较高的对比度，如(a)。图9是一个自然场景的例子，其中一个色调映射的图像(a)更好地保存了结构细节，而另一个(b)给出了更自然的整体外观(但失去了结构信息，特别是在最亮的区域)。(c)、(d)、(e)分别给出了三种不同的图像融合算法生成的三幅融合图像。该方法生成的图像在结构保留性和统计自然性之间取得了最佳的平衡，并在TMQI中获得了最佳的质量分数。

为了进一步验证所提出的融合方案，我们进行了一个额外的主观实验，邀请10名受试者对5组色调映射图像进行排序，每组包含8幅图像。其中7张图像是使用第三节第三个实验中使用的TMOs生成的。在这7个TMOs中选择其中的2个，使用所提出的融合方法生成第8幅图像。表四比较了源图像和相应融合图像的平均主观评分，评分越低，质量越好。可以看出，融合后的图像几乎总是显著高于被融合的两个源图像。

V. CONCLUSION

我们开发了一个客观的模型来评估色调映射图像的质量，通过结合多尺度结构保真度测量和统计自然度测量。该方法不仅提供了图像的整体质量评分，而且还创建了反映结构保真度在尺度和空间上变化的多尺度质量地图。我们的实验表明，TMQI与图像质量的主观评价有一定的相关性。此外，我们还演示了TMQI在色调映射算法的自动参数调整和融合多个色调映射图像方面的有用性。

作为该研究课题的首批尝试之一，我们的方法有几个局限性，可能在未来得到解决或改进。首先，TMQI仅用于评估灰度图像，但大多数自然场景的HDR图像都是彩色的。评估色调映射彩色图像的一种简单方法是将TMQI单独应用于每个颜色通道，然后将它们组合起来。颜色保真度和颜色自然度的测量可以改进质量的测量。

其次，在当前的结构保真度图的池化方法中采用了简单平均的方法。结合视觉注意模型的先进的池化方法可以提高质量预测的性能。

第三，目前的统计自然度测量仅基于强度统计。在自然图像统计方面已有丰富的文献[28]，可以引入先进的统计模型(反映自然图像在空间、尺度和方向上的结构规律)来提高统计自然度的度量。

第四，利用TMQI作为新的优化目标，可以对许多现有的TMOs进行重新设计，以获得更好的图像质量。新的TMOs也可以通过利用所提议的质量评估方法的构建来开发。

最后，对目前的方法仅使用自然图像进行了应用和测试。HDR图像和TMOs的应用范围超出了自然图像。例如，现代医学成像设备经常捕获需要在可视化之前进行色调映射的HDR医学图像。TMQI和优化方法可以适用于这些扩展的应用程序。

代码

#!/usr/bin/python3
# -*- coding: utf-8 -*-

from __future__ import print_function, division, absolute_import
import numpy as np
from scipy.signal import convolve, gaussian
from scipy.ndimage.filters import generic_filter
from scipy.stats import norm, beta
from contracts import contract
from skimage.util import view_as_blocks
import sys

__version__ = "0.10.0"
__title__ = "TMQIr"
__summary__ = "TMQI revised"
__uri__ = "https://github.com/dvolgyes/TMQI"

# for the Python reimplementation, original authors in 'upstream'
__author__ = "David Völgyes"
__email__ = "david.volgyes@ieee.org"

__derived__ = True   # Meaning: reimplementation with deviation
__upstream_license__ = "BSD-like"  # see the website for exact details
__upstream_uri__ = "https://ece.uwaterloo.ca/~z70wang/research/tmqi/"
__upstream_doi__ = "10.1109/TIP.2012.2221725"
__upstream_ref__ = ('H. Yeganeh and Z. Wang,' +
                    '"Objective Quality Assessment of Tone Mapped Images,"' +
                    'IEEE Transactions on Image Processing,' +
                    'vol. 22, no. 2, pp. 657-667, Feb. 2013.')


class Metric(object):

    def __init__(self, *args, **kwargs):
        self.name = "Undefined"
        self.descriptors = list()
        self.lists = tuple()
        self.maps = tuple()
        self.no_reference = False
        self.full_reference = False
        self.luminance = False
        self.RGB = False
        self.cache = dict()

    @contract(RGB='array[NxMx3](float)')
    def _RGBtoY(self, RGB):
        M = np.asarray([[0.2126, 0.7152, 0.0722], ])
        Y = np.dot(RGB.reshape(-1, 3), M.T)
        return Y.reshape(RGB.shape[0:2])


def img_read(link, gray=False, shape=None, dtype=None, keep=False):
    if os.path.exists(link):
        if dtype is None:
            img = imread(link)
            if gray and len(img.shape) > 2:
                img = skimage.color.rgb2hsv(img)[..., 2]
        else:
            W, H = shape
            img = np.fromfile(link, dtype=dtype)
            if gray:
                img = img.reshape(H, W)
            else:
                img = img.reshape(H, W, -1)
    else:
        tempfile = wget.download(link, bar=None)
        img = img_read(tempfile, gray, shape, dtype)
        if not keep:
            os.remove(tempfile)
    return img.astype(np.float)


class TMQI(Metric):

    def __init__(self, *args, **kwargs):
        super().__init__()
        self.name = "TMQI"
        self.descriptors = ("Q", "S", "N")
        self.lists = ("s_local",)
        self.maps = ("s_local",)
        self.no_reference = False
        self.full_reference = True
        self.luminance = True
        self.RGB = False
        self.original = True

        if len(args) + len(kwargs) > 0:
            self.__call__(*args, **kwargs)

    @contract(hdrImage='array[NxMx3](float)|array[NxM](float),N>10,M>10',
              ldrImage='array[NxMx3](float)|array[NxM](float)')
    def __call__(self, hdrImage, ldrImage, window=None):
        # images must have same dimenions
        assert hdrImage.shape == ldrImage.shape

        if len(hdrImage.shape) == 3 and len(ldrImage.shape) == 3:
            # Processing RGB images
            L_hdr = self._RGBtoY(hdrImage)
            L_ldr = self._RGBtoY(ldrImage)
            return self._TMQI_gray(L_hdr, L_ldr, window)

        # input is already grayscale
        return self._TMQI_gray(hdrImage, ldrImage, window)

    @contract(hdrImage='array[NxM](float)',
              ldrImage='array[NxM](float)',
              window='None|array[UxV],U<N,V<M,U>=2,V>=2')
    def _TMQI_gray(self, hdrImage, ldrImage, window=None):
        a = 0.8012
        Alpha = 0.3046
        Beta = 0.7088
        lvl = 5  # levels
        weight = [0.0448, 0.2856, 0.3001, 0.2363, 0.1333]

        M, N = hdrImage.shape

        if window is None:
            gauss = gaussian(11, 1.5)
            window = np.outer(gauss, gauss)

        # unnecessary, it is just for the sake of parallels with the matlab code
        L_hdr = hdrImage
        L_ldr = ldrImage

        # Naturalness should be calculated before rescaling
        N = self._StatisticalNaturalness(ldrImage)

        # The images should have the same dynamic ranges, e.g. [0,255]

        factor = float(2**32 - 1.)

        if self.original:
            L_hdr = factor * (L_hdr - L_hdr.min()) / (L_hdr.max() - L_hdr.min())
        else:
            # but we really should scale them similarly...
            L_hdr = factor * (L_hdr - L_hdr.min()) / (L_hdr.max() - L_hdr.min())
            L_ldr = factor * (L_ldr - L_ldr.min()) / (L_ldr.max() - L_ldr.min())

        S, s_local, s_maps = self._StructuralFidelity(L_hdr, L_ldr, lvl, weight, window)
        Q = a * (S ** Alpha) + (1. - a) * (N ** Beta)
        return Q, S, N, s_local, s_maps,

    @contract(L_hdr='array[NxM](float),N>0,M>0',
              L_ldr='array[NxM](float),N>0,M>0')
    def _StructuralFidelity(self, L_hdr, L_ldr, level, weight, window):

        f = 32
        s_local = []
        s_maps = []
        kernel = np.ones((2, 2)) / 4.0

        for _ in range(level):
            f = f / 2
            sl, sm = self._Slocal(L_hdr, L_ldr, window, f)

            s_local.append(sl)
            s_maps.append(sm)

            # averaging
            filtered_im1 = convolve(L_hdr, kernel, mode='valid')
            filtered_im2 = convolve(L_ldr, kernel, mode='valid')

            # downsampling
            L_hdr = filtered_im1[::2, ::2]
            L_ldr = filtered_im2[::2, ::2]

        S = np.prod(np.power(s_local, weight))
        return S, s_local, s_maps

    @staticmethod
    @contract(img1='array[NxM](float),N>0,M>0',
              img2='array[NxM](float),N>0,M>0',
              sf='float,>0')
    def _Slocal(img1, img2, window, sf, C1=0.01, C2=10.):

        window = window / window.sum()

        mu1 = convolve(window, img1, 'valid')
        mu2 = convolve(window, img2, 'valid')

        mu1_sq = mu1 * mu1
        mu2_sq = mu2 * mu2
        mu1_mu2 = mu1 * mu2

        sigma1_sq = convolve(img1 * img1, window, 'valid') - mu1_sq
        sigma2_sq = convolve(img2 * img2, window, 'valid') - mu2_sq

        sigma1 = np.sqrt(np.maximum(sigma1_sq, 0))
        sigma2 = np.sqrt(np.maximum(sigma2_sq, 0))

        sigma12 = convolve(img1 * img2, window, 'valid') - mu1_mu2

        CSF = 100.0 * 2.6 * (0.0192 + 0.114 * sf) * np.exp(- (0.114 * sf) ** 1.1)
        u_hdr = 128 / (1.4 * CSF)
        sig_hdr = u_hdr / 3.

        sigma1p = norm.cdf(sigma1, loc=u_hdr, scale=sig_hdr)

        u_ldr = u_hdr
        sig_ldr = u_ldr / 3.

        sigma2p = norm.cdf(sigma2, loc=u_ldr, scale=sig_ldr)

        s_map = ((2 * sigma1p * sigma2p + C1) / (sigma1p**2 + sigma2p**2 + C1)
                 * ((sigma12 + C2) / (sigma1 * sigma2 + C2)))
        s = np.mean(s_map)
        return s, s_map

    @contract(L_ldr='array[NxM](float),N>0,M>0', win='int,>0')
    def _StatisticalNaturalness(self, L_ldr, win=11):
        phat1 = 4.4
        phat2 = 10.1
        muhat = 115.94
        sigmahat = 27.99
        u = np.mean(L_ldr)

        # moving window standard deviation using reflected image
        if self.original:
            W, H = L_ldr.shape
            w_extra = (11 - W % 11)
            h_extra = (11 - H % 11)
            # zero padding to simulate matlab's behaviour
            if w_extra > 0 or h_extra > 0:
                test = np.pad(L_ldr, pad_width=((0, w_extra), (0, h_extra)), mode='constant')
            else:
                test = L_ldr
            # block view with fixed block size, like in the original article
            view = view_as_blocks(test, block_shape=(11, 11))
            sig = np.mean(np.std(view, axis=(-1, -2)))
        else:
            # deviation: moving window with reflected borders
            sig = np.mean(generic_filter(L_ldr, np.std, size=win))

        beta_mode = (phat1 - 1.) / (phat1 + phat2 - 2.)
        C_0 = beta.pdf(beta_mode, phat1, phat2)
        C = beta.pdf(sig / 64.29, phat1, phat2)
        pc = C / C_0
        B = norm.pdf(u, muhat, sigmahat)
        B_0 = norm.pdf(muhat, muhat, sigmahat)
        pb = B / B_0
        N = pb * pc
        return N


class TMQIr(TMQI):

    def __init__(self, *args, **kwargs):
        super().__init__()
        self.name = "TMQIrev"
        self.descriptors = ("Q", "S", "N")
        self.lists = ("s_local",)
        self.maps = ("s_local",)
        self.no_reference = False
        self.full_reference = True
        self.luminance = True
        self.RGB = False
        self.original = False


if __name__ == "__main__":
    if len(sys.argv) == 1:
        import doctest
        doctest.testmod()

    if len(sys.argv) > 1:  # there are command line parameters
        # these imports are unnecessary if the code is used as a library
        from optparse import OptionParser
        from scipy.misc import imsave
        from imageio import imread
        import os.path
        import wget
        import skimage.color

        usage = ("usage: %prog [options] HDR_image LDR_image\n" +
                 "The images could be files or a http(s)/ftp link.")
        parser = OptionParser(usage=usage)

        parser.add_option("-t", "--type",
                          type="string",
                          dest="maptype",
                          help="s_map file type (default: float32)",
                          default="float32")

        parser.add_option("-m", "--smap_file",
                          type="string",
                          dest="smap",
                          help="s_map file name prefix. (default: s_map_)",
                          default="s_map_")

        parser.add_option("-p", "--precision",
                          type="int",
                          dest="precision",
                          help="precision (number of decimals) (default: 4)",
                          default=4)

        parser.add_option("-W", "--width",
                          type="int",
                          dest="width",
                          help="image width (mandatory for RAW files)"
                          " (default: None)",
                          default=None)

        parser.add_option("-H", "--height",
                          type="int",
                          dest="height",
                          help="image height (mandatory for RAW files)"
                          " (default: None)",
                          default=None)

        parser.add_option("-i", "--input_type",
                          type="string",
                          dest="input",
                          help="type of the input images: float32/float64"
                          " for RAW images\n"
                          "None for regular images opening with scipy"
                          " (e.g. png) (default: None)",
                          default=None)

        parser.add_option("-g", "--gray",
                          dest="gray",
                          action="store_true",
                          help="gray input (ligthness/brightness)"
                          "  (default: RGB)",
                          default=False)

        parser.add_option("-Q", "--report-Q",
                          dest="report_q",
                          action="store_true",
                          help="report quality index",
                          default=True)

        parser.add_option("-S", "--report-S",
                          dest="report_s",
                          action="store_true",
                          help="report structural similarity",
                          default=False)

        parser.add_option("-L", "--report-SL",
                          dest="report_sl",
                          action="store_true",
                          help="report maps (S_locals)",
                          default=False)

        parser.add_option("-N", "--report-N",
                          dest="report_n",
                          action="store_true",
                          help="report naturalness",
                          default=False)

        parser.add_option("-M", "--report-MAPS",
                          dest="report_maps",
                          action="store_true",
                          help="report maps",
                          default=False)

        parser.add_option("-q", "--no-report-Q",
                          dest="report_q",
                          action="store_false",
                          help="do not report quality index")

        parser.add_option("-s", "--no-report-S",
                          dest="report_s",
                          action="store_false",
                          help="do not report structural similarity")

        parser.add_option("-l", "--no-report-SL",
                          dest="report_sl",
                          action="store_false",
                          help="do not report maps (S_locals)")

        parser.add_option("-n", "--no-report-N",
                          dest="report_n",
                          action="store_false",
                          help="do not report naturalness")

        parser.add_option("--quiet",
                          dest="quiet",
                          action="store_true",
                          help="suppress variable names in the report")

        parser.add_option("--verbose",
                          dest="quiet",
                          action="store_false",
                          help="use variable names in the report (default)",
                          default=False)

        parser.add_option("--keep",
                          dest="keep",
                          action="store_true",
                          help="keep downloaded files (default: False)",
                          default=False)

        parser.add_option("-r","--revised",
                          dest="revised",
                          action="store_true",
                          help="Enable revised TMQI. (default: Original)")

        (options, args) = parser.parse_args()

        if len(args) != 2:
            print("Exactly two input files are needed: HDR and LDR.")
            sys.exit(0)

        if options.input is not None:
            W, H = options.width, options.height
            shape = (W, H)
            dtype = np.dtype(options.input)
        else:
            shape = None
            dtype = None

        hdr = img_read(args[0], options.gray, shape, dtype, options.keep)
        ldr = img_read(args[1], options.gray, shape, dtype, options.keep)

        if options.revised:
            Q, S, N, s_local, s_maps = TMQIr()(hdr, ldr)
        else:
            Q, S, N, s_local, s_maps = TMQI()(hdr, ldr)

        prec = options.precision
        Q, S, N = np.round(Q, prec), np.round(S, prec), np.round(N, prec)
        s_local_str = " ".join(map(str, np.round(s_local, prec)))

        result = ""

        if options.report_q:
            if not options.quiet:
                result += "Q: "
            result += "%s " % Q

        if options.report_s:
            if not options.quiet:
                result += "S: "
            result += "%s " % S

        if options.report_n:
            if not options.quiet:
                result += "N: "
            result += "%s " % N

        if options.report_sl:
            if not options.quiet:
                result += "S_locals: "
            result += "%s " % s_local_str
        print(result.strip())

        if options.report_maps:
            for idx, sm in enumerate(s_maps):
                filename = "%s%i.%s" % (options.smap, idx + 1, options.maptype)

                try:
                    out = sm.astype(options.maptype)
                    out.tofile(filename)
                except TypeError:
                    imsave(filename, sm)