A Hybrid ℓ1-ℓ0 Layer Decomposition Model for Tone Mapping—CVPR论文阅读
A Hybrid ℓ1-ℓ0 Layer Decomposition Model for Tone Mapping
目录
A Hybrid ℓ1-ℓ0 Layer Decomposition Model for Tone Mapping
摘要
色调映射旨在从保留了视觉信息的高动态范围图像中复制标准动态范围图像。最新的色调映射算法主要将图像分解为基础层和细节层，并进行相应处理。由于缺乏对这两层的适当先验，这些方法可能具有光晕伪像和过度增强的问题。在本文中，我们提出了混合ℓ1-ℓ0分解模型来解决这些问题。具体来说，在基础层上施加一个ℓ1稀疏项以模拟其分段平滑性。 稀疏项作为结构先验被施加到细节层，这导致分段恒定效应。我们还根据我们的层分解模型提出了一种多尺度色调映射方案。实验表明，我们的色调映射算法在几乎没有光晕伪影的情况下达到了视觉上令人信服的结果，在主观和客观评估中均优于最新的色调映射算法。
1、介绍
现实世界中的场景可能会跨越亮度动态范围，该范围大大超出了大多数成像设备的响应范围[4]。 由于近十年来高动态范围（HDR）技术的飞速发展，可以通过包围曝光融合技术将场景的完整信息记录在辐射图中[2，7]。 然而，大多数显示设备具有有限的动态范围，并且不能如实地再现辐射度图中的信息。 因此，需要一种有效的色调映射算法来将HDR辐射度映射图转换为标准动态范围（SDR）图像而又不牺牲主要的视觉信息。
在过去的二十年中，文献中已经提出了大量的色调映射方法。尽管设计方法各不相同，但这些色调映射方法中的很大一部分还是基于层分解[8,14,23,29]。具体来说，图像被分解为基础层和细节层，然后分别进行处理。细节层得以保留或增强[8，14]，具有大空间平滑度和高范围变化的基础层被压缩。尽管大多数基于层分解的色调映射算法都可以在某种程度上提高辐射图的视觉可解释性，但是它们在获得自然和视觉上令人愉悦的结果方面仍然有局限性。一个典型的问题是小尺度纹理细节的过度增强。这是因为现有作品通常会忽略细节层的空间属性，这对色调映射图像有重大影响。此外，由于缺乏基本层的边缘保留特性，因此在某些色调映射算法中，光晕伪影也是一个问题[14]。为了获得辐射图的自然且无伪影的复制，必须将某些适当的先验条件合并到层分解框架中。
鉴于在HDR亮度图中记录了大量信息的事实，应该为视觉感知分配哪一部分信息具有较高的优先级，这是色调映射的重要问题。 在心理学中，发现人类视觉对边缘更为敏感[1，13]。 这种视觉机制有助于捕获场景的主要语义信息。 在本征分解[3，6]的研究中，通常假设反射层（类似于细节层的概念）中的边缘稀疏，这也表明了结构信息在图像中的重要性。 鉴于上述观察，色调映射运算符应首先解决结构重现问题。 由于图层分解框架中细节层的空间特性在很大程度上影响色调映射图像的视觉外观，因此我们考虑在细节层上先施加结构稀疏性。
尽管色调映射研究中很少报道在细节层使用空间先验，但Retinex分解[12，25]中长期采用ℓ1稀疏先验来对反射层的结构稀疏进行建模。 尽管ℓ1项保留了图像中的边缘，但其分段平滑性会导致较弱的结构先验。 另一方面，ℓ0稀疏项显示出很大的分段平化特性[34]，对于结构先验而言似乎是更好的选择。
在本文中，我们提出了一个混合的ℓ1-ℓ0层分解模型用于色调映射。 具体来说，对细节层施加ℓ0梯度稀疏项以对结构先验建模。 这样，细节层将主要包含结构信息，该信息将得到增强。 同时，为了减少光晕伪像，对基础层施加了ℓ1的梯度稀疏项以保留边缘。 基于我们的分解模型，开发了一种多尺度色调映射方案。 由于在层分解中使用了适当的先验，因此在主观和客观评估方面，我们的色调映射器均优于最新的算法。
本文的组织如下。 第2节回顾了一些相关的工作。 第3节介绍了建议的层分解模型。 第4节总结了我们的多尺度色调映射算法。第5节和第6节分别是实验和结论。
2、相关工作
我们的工作主要涉及色调映射、基于Retinex的层分解和边缘感知滤波。
色调映射。现有的色调映射算法可以分为全局方法和局部方法。全局色调映射方法可复制具有单个压缩曲线的SDR图像[28、32、33]。相反，局部色调映射方法以空间变化的方式执行此任务，并且在细节方面更好。局部方法通常基于层分解，其中基础层首先由保留边缘的滤波器估计，而细节层是原始图像与基础层之间的残差。不同的本地色调映射算法主要在滤波器设计技术上有所不同。在早期阶段，采用了基于内核的过滤器。 Reinhard等。建议使用具有空间自适应比例参数的基于高斯的滤波器[29]。杜兰德等。采用双边滤波器来估计基础层[8]。尽管此方法可以在某种程度上避免出现光晕伪影，但它会通过增强小尺度细节来过度增强图像。 Li等。提出了一种用于色调映射的多尺度小波方案[18]。梅兰（Meylan）等人。提出了一种基于Retinex的自适应滤波器，用于色调映射[23]。文献[14]中提出了一种用于色调映射的加权导引滤波器，由于小尺度细节的过度增加，它也存在过度增强的问题。还提出了基于全局优化的滤波器用于色调映射。 Farbman等。提出了加权最小二乘（WLS）滤波器[10]。该过滤器具有出色的平滑效果和强大的边缘保留特性。其他色调映射算法包括全局线性窗口方法[30]和基于PCA的方法[17]。
摘要现有的基于分层分解的色调映射方法在基础层上增加了保留边缘的优先权，而对细节层的关注较少。与此相反，我们的分解框架在细节层上强加了结构优先，以提高结果的视觉质量。
基于Retinex的分解。尽管Retinex分解最初源自视觉稳定性研究[16]，但它从单个图像估计光照和反射。 Retinex分解通常被公式化为具有不同先验性和反射性的变分模型。在开创性的工作中[15]，Kimmel等人。提出了一个基于ℓ2的Retinex分解模型，用于对比度增强，其中照明和反射被假定为全局平滑。 Ng等。假设反射层是分段平滑度，并用总变化项替代了ℓ2范数[25]。梁等。假设照明是分段平滑的，并提出了一种基于非线性扩散的照明估计方法[19]。该方法保留了照明层中的边缘并抑制了结果中的光晕伪影。最近，傅等人。提出在反射层上按亮度倒数加权的权重项，以对反射的分段常数假设建模[12]。
Edge-preserving平滑。边缘感知平滑是图像处理中的一项基本技术。最早的保边滤波器是考虑图像[9]局部距离变化的双边滤波器。Min等人提出了一种基于加权最小二乘[24]的快速全局平滑算法。其他代表过滤器是ℓ[34]基于滤波器和基于加权-ℓ过滤器[3]。
3、层分解方法
我们首先提出一个混合ℓ1 -ℓ0层分解模型并给出了求解器进行求解。然后，我们将这种分解方法扩展到一个多尺度框架，其中不同的组成部分的图像可以操纵的色调映射。
3.1 混合ℓ1 -ℓ0层分解模型
为了设计合适的层分解框架，我们建议在细节层上加上结构先验，在基础层上加上保边先验。分别用S、B和S−B表示原始图像、基础层和细节层。提出的分层优化模型如下:
其中p是像素索引，N是图像中的像素数。 第一项（Sp -Bp）2使基础层接近原始图像。 基础层的空间特性公式为ℓ1梯度稀疏项|∂iBp|，i = x，y，其中∂i是沿x或y方向的偏导数。 细节层的空间特性被表示为具有指示函数F（x）的spa0梯度稀疏项：
我们的层分解模型的优点在于ℓ1和ℓ0正则化的混合用法。 一方面，由于ℓ1稀疏度项的异常值排除性质[20]，基层的大梯度得以保留。 因此，基层是分段光滑的。 另一方面，已经证明ℓ0稀疏项会产生fl化效应[26，34]。 我们的模型应用ℓ0项来强制细节层的小纹理渐变为零，同时保持主要结构渐变不变。 如图1（b）所示，这种安排可产生分段恒定效应并成功地对结构先验建模。
细节层的另一个可能的选择是ℓ1梯度稀疏度优先，这已在Retinex研究中报道[12，25]。在[12]中，对反射/细节层施加fl1项以获得分段恒定效应。但是，ℓ1项有两个缺点。首先，其分段平滑度[21]的性质不足以产生分段恒定结果，如图1（c）所示。其次，在相同的参数设置下，ℓ1项不能强烈地规范细节层，这可能导致色调映射图像的过度增强，如图1（e）所示。为了显示ℓ1项和ℓ0项之间的差异，图1（f）中显示了从其合成的细节层提取的一维轮廓信号。信号的位置在图1（a）中用黄线表示。我们可以看到，ℓ0项使小的琐碎变化变小并保留了视觉上重要的边缘，而the1项对此无效。结果，使用ℓ0项可以避免过度增强的问题，并增加图像的视觉可解释性，如图1（d）所示。
3.2 模型求解
目标函数(1)凸是因为ℓ0标准正规化。采用交替方向乘法器(ADMM)框架[5]求解该优化模型。由于空间有限，我们只对每个子问题的求解做了简要的介绍。详细说明请参考补充材料。
为了清晰起见，我们首先将目标函数(1)以矩阵向量形式重写为:
其中s，b∈RN分别是（1）中S，B的级联向量形式，1∈R2N是全1的向量。 ▽表示两个梯度算子矩阵的级联▽= [▽⊤x，▽⊤y]⊤∈R2N×N。 F（▽（s-b））执行逐元素非零指示并输出二进制向量。 现在引入两个辅助变量c1，c2∈R2N分别替换▽b，▽（sb）。 我们模型得出的增强拉格朗日函数是：
其中yi,i = 1,2为拉格朗日对偶变量。在k次迭代中，通过最小化几个原始子问题和最大化对偶问题来优化函数(4)。
(1) Solving bk+1:
如果我们把向量ck 1分成两条等长ck - 1, 1, ck 1, 2, ck 2分割成ck 2、1和ck 2, 2,即1分割成yk 1, 1,即1、2、2 yk分割成2 yk, 1,即2,2,目标函数关于bk + 1是一个二次规划问题,可以有效地解决通过FFT变换(指我们辅料)。
(2) Solving C1k+1 :
关于C1k+ 1的目标函数可以通过软缩法求解：
(3) Solving C2k+1:
根据对[34]的分析，关于C2k + 1的目标函数可以逐项求解。 这等于求解N个独立的标量函数。 用下标j表示向量的第j个条目。
3.3 扩展到多尺度分解
通过将混合ℓ1-ℓ0分解模型（1）应用于亮度图，我们可以生成分段常量细节层和分段平滑基础层。 尽管此单尺度方案为色调映射提供了标准框架，但将分解应用于基础层反复导致多尺度分解，从而可以进一步改善色调映射结果。 以此方式，可以不同地操纵由不同比例图层表示的图像的不同属性，从而获得更加灵活有效的色调再现。 通过利用效率和有效性，我们采用了两级分解方案进行色调映射，如图2所示。它将产生一个1级细节层D1，一个2级细节层D2和一个2级基础 B2层。
如3.1节所述，D1的空间性质对色调映射图像有很大影响。我们应用该ℓ1 -ℓ0模型(1)第一个规模分解:
模型ℓ1,ℓ0(·)是(1)优化模型。第一级分解后,在细节层结构信息仍然D1和主要结构信息转移到基地层B1。
这种简化是基于这样一个事实，即我们的目标是在scale-2细节层D2中保留图像的纹理信息。因此,ℓ基于结构之前,不适用在这个分解的规模。因此，层D2存储了大部分的纹理信息，层B2包含了局部平均亮度。
综上所述，我们的两尺度分解方案产生了三层D1、D2和B2，满足:
图3显示了色调映射的结果我们的模型规模1和2尺度(我们的色调映射算法的细节将在第四部分讨论)。可以看出,尽管地结合的结果是可以接受的,双尺度结果保留更好的图像的中频分量和达到更自然的外观。
加速。对第二次尺度分解(9)的准确性没有严格要求。因此，我们采用加速方案。首先，我们线性向下采样B1层4倍。然后利用(9)中的分解模型得到B2的低分辨率图像，再对原始分辨率进行线性上采样。由于采样方案使图像中的边界区域略有模糊，最终我们以原始B1作为引导图像对B2进行快速联合双边滤波，得到清晰的边界信息[27]。
4. Tone Mapping
基于提出的层分解的结果，开发了一种色调映射算法，其主要步骤包括颜色转换，多尺度分解，细节层增强，基础层压缩和层重组。 尽管此框架在色调映射研究中很常见，但我们的方法在两个方面有所不同。 首先，我们的图层分解模型适用于图像的空间属性。 如第3.3节所述，我们的多尺度分解将结构信息，纹理信息和局部平均亮度分别部署到不同的层中，而现有的多尺度模型仅执行渐进平滑[10，14]。 其次，在我们的多尺度操作方法中，我们执行了层选择非线性处理，而其他工作仅执行了线性强度缩放[10]。
由于图像的动态范围主要嵌入在亮度域中，因此我们的核心算法仅处理亮度通道并保留色度分量。 具体来说，将输入的RGB辐射度图转换为HSV空间，并且仅将V通道进行色调映射。 在逆变换阶段，将饱和通道乘以0.6，以防止过饱和。
我们在辐射度图的亮度通道上的色调映射算法如图2所示。辐射度图的通道Vh首先转换为对数域，并归一化为（0，1）的范围。 此步骤模拟了人类视觉对亮度的响应，并初步减小了动态范围。 然后应用使用（8）和（9）的两尺度分解模型，得出三层D1，D2和B2。 由于基本层B2可被视为图像的局部亮度水平，因此我们通过伽玛函数对其进行压缩：
其中L是最大的亮度级别(在我们的例子中，由于归一化，L = 1)。对于首尺度细节层D1，我们使用非线性拉伸函数对其进行增强:
该函数与参数的旋向对定心为0的信号有拉伸作用。更小的内径可产生更大的拉伸度，反之亦然。由于分解模型(1)对D1施加了结构先验，因此通过拉伸函数来增强原始图像的结构残差。这样的安排会产生一个视觉上更吸引人的图像。然后，通过重构一个亮度SDR图像。
最后，将0.5%和99.5%强度下的V值分别映射为0和1。超出此范围的值将被剪切。
5. 实验和分析
本节提出了几个实验来验证我们的性能混合ℓ1 -ℓ0层分解模型(1)和拟议的色调映射算法。HDR数据库包含40幅辐射地图，这些地图来自不同的来源1 2用于评估。这40幅图片涵盖了室内和室外场景，有不同类型的物体，包括植物、汽车、天空和建筑。
5.1 参数选择
主要参数影响我们ℓ1 -ℓ0λ1分解模型(1),和λ2控制平滑程度上基本层和细节后,分别。接下来，我们使用平均局部熵(MLE)3来客观地测量两层的平滑度。MLE越大，平滑度越低，即。，图像中有更多纹理，反之亦然。
图4显示了在固定了模型模型1后，模型模型2对细节层的影响。从图中可以看出，不同的模态2值会导致不同程度的平坦/平滑效果D1。当bioc 2过大(0.008)时，部分结构完全扁平，MLE较低(0.97)。相比之下，当贴图2太小(0.0008)时，D1中出现了一些小的纹理梯度，MLE较大(2.33)，结构先验被较少表示。我们对我们的数据库进行了广泛的实验，发现当将icsp2设置为0.01 icsp1时，分解效果始终令人满意。图5为当风害2固定为0.01风害1时，参数风害1的影响。可以看出，主要控制D1的信号大小，而对的分段平滑程度控制较小。
其他需要确定的参数有(9)中的thomas3、(11)中的和(12)中的。在最终的基层B2中，控制层的平滑程度。我们发现，除了一些极端的设置，音素3不影响太多的色调映射图像。因此，贴图3被固定为0.1。主要控制第一细节层D1的拉伸程度。为避免过度放大，我们将其设置为中等值0.8。最后，将 γ设为2.2，作为Retinex分解研究的常用方法[12,15,25]。
5.2 分解层
为了验证我们的音调映射算法的多尺度分解性能，我们与Gu的多尺度音调映射器[14]进行了比较。在Gu的模型中，将梯度函数加权的局部引导滤波器反复应用于原始图像，得到2尺度(3层)的分解。注意，尽管Gu的模型声称有3个比例(4层)，最后一个比例的基础层是一个恒定的图像。因此，有效的尺度数是2。Gu的模型在基本层上加强了边缘保持特性，而没有在细节层上强加任何优先权。
在图6中，比较了Gu模型和我们模型的多尺度分解结果。一维辅助分析如图7所示，其中从每种方法的分解层中提取了一个一维轮廓信号（位置在图7（a）中用白线表示）。从图7（b）可以看出，顾氏模型在不考虑细节层的空间特性的情况下进行了渐进平滑。因此，如图6（d）所示，第一细节层（图7（b）中的红色曲线）充满了小波动，并且色调映射的图像被过度增强。此外，由于局部过滤的特性，顾的模型没有严格保留边缘。因此，色调映射结果具有光晕伪像（请参见图6（d）中的放大图）。相反，由于结构先验，我们的方法在第二层D2中分配了小规模的变化，并强制第一层D1呈分段恒定，如图7（c）所示。同时，我们的方法也是边缘保留的。如图6（h）所示，它不仅避免了光晕伪影，而且还获得了令人瞩目的视觉效果。
5.3 色调映射比较
我们比较我们的色调映射器与最先进的色调映射器[10,11,14,22,30,31]收集的数据库。这些色调映射器包括基于WLS滤波器的[10]方法、全局线性窗方法(GLW)[30]、视觉适应方法(VAD)[11]、向后兼容方法(BWC)[22]、引导滤波方法(GF)[14]、梯度重建方法(GR)[31]。更多的对比结果可以在补充文件中找到。GF是由我们实现的，因为源代码不可用。BWC是由pfstool4实现的。其他的是由作者的源代码实现的。
主观评价。 图8、9显示了两个图像上色调映射结果的比较。 我们可以看到，我们的方法在细节增强和自然保留之间实现了良好的平衡。 相反，其他色调映射器遭受不同类型的失真。 WLS失去局部对比度，而GLW遭受亮度失真。 VAD存在色偏问题，并且BWC过度柔化了图像。 GR和GF具有过度增强问题和光晕伪像。 在图10中，我们的色调映射器与Photomatix5的默认色调映射器进行了比较。 我们可以看到，两种方法都能获得令人满意的结果，而由于突出了结构信息，我们的方法获得了更高的视觉解释性。
客观评价。除了主观评估，我们使用色调映射图像质量指数（TMQI）[35]对色调映射器执行客观评估。 TMQI首先评估色调映射图像的结构保真度和自然度。然后，通过幂函数对这两个量度进行调整并取平均值，以得到最终分数，范围为0到1。TMQI的值越大，表示色调映射图像的质量越好，反之亦然。表1说明了在我们的数据库中使用40张HDR图像执行的每个色调映射器的平均TMQI得分。我们可以看到，我们的方法不仅获得了最高的TMQI分数（0.8851），而且还获得了最高的自然度指标（0.5547）。这些出色的标记客观地表明了我们的算法所获得的高视觉质量。另一方面，我们的色调映射器没有达到很高的保真度分数。这是因为保真度度量在不同比例的局部窗口中计算标准偏差。但是，我们的算法会规范化小规模细节，以避免过度增强，从而降低了保真度得分。
效率。所提出的色调映射器具有中等的计算复杂度。基于admm的求解器中最复杂的部分是FFT操作，其代价为O(Nlog(N))。表2比较了5个色调绘制器在1333×2000大小的图像上的运行时间(图8(a))。测试环境为一台PC, CPU为i7 6850k, RAM为16G。可以看出，与其他方法相比，我们的色调映射器运行时间适中。
6、结论
本文提出了一种新的混合ℓ1-ℓ0层分解模型，以解决色调映射的过度增强和光晕伪影问题。 此分解模型有效地在细节层之前强制执行结构，并在基础层之前强制执行边缘保留。 采用ADMM算法有效地求解了分解模型。 基于ℓ1-ℓ0层分解输出，提出了一种多尺度色调映射算法。 它在基础层执行动态范围缩小，在细节层执行结构增强。 由于两个先验的正确使用，我们的多尺度色调映射算法不仅避免了光晕伪影，而且还比现有作品获得了更具视觉吸引力的色调映射结果。
matlab代码
用到的函数：
function y = R_func(x,g,sigma,a,b)
index_low = abs(x-g) <= sigma;
temp = x(index_low);
temp_low = g + fd( abs(temp - g)./sigma, a) .* sigma .* sign(temp - g);
% temp_low = g + betainc( abs(temp - g)./sigma,a*2,a*4) .* sigma .* sign(temp - g);
index_high = abs(x-g) > sigma;
temp = x(index_high);
temp_high = g + (fe( abs(temp - g) - sigma, b) + sigma) .* sign(temp - g);
y = zeros(size(x));
y(index_low) = temp_low;
y(index_high) = temp_high;
function y = fd(x,a)
y = x.^a;
end
function y = fe(x,b)
y = b*x;
end
endfunction y = nor(x)
y = (x-min(x(:)))/(max(x(:))-min(x(:)));
endfunction [D1,D2,B2] = Layer_decomp(img,lambda1,lambda2,lambda3)
%% first scale decomposition
[hei,wid,~] = size(img);
[D1,B1] = L1L0(img,lambda1,lambda2);
%% second scale decomposition
scale = 4;
B1_d = imresize(B1,[round(hei/scale),round(wid/scale)],'bilinear');
[~,B2_d] = L1(B1_d,lambda3);
B2_r = imresize(B2_d,[hei,wid],'bilinear');
B2 = bilateralFilter(B2_r,nor(B1),0,1,min(wid,hei)/100,0.05);
D2 = B1 - B2;
end%% This function performs the hybrid L1-L0 decomposition
function [D,B]= L1L0(S,lambda1,lambda2)
iter = 15;
[hei,wid] = size(S);
fx = [1, -1];
fy = [1; -1];
otfFx = psf2otf(fx,[hei,wid]);
otfFy = psf2otf(fy,[hei,wid]);
DxDy = abs(otfFy).^2 + abs(otfFx).^2;
B = S;
C = zeros(hei,wid*2);
E = zeros(hei,wid*2);
L1 = zeros(hei,wid*2);
L2 = zeros(hei,wid*2);
ro1 = 1;
ro2 = 1;
DiffS = [-imfilter(S,fx,'circular'),-imfilter(S,fy,'circular')];
for i = 1:iter
CL = C + L1./ro1;
EL = DiffS - E - L2./ro2;
%% for B
C1L1 = CL(:,1:wid);
C2L2 = CL(:,1+wid:end);
E1L3 = EL(:,1:wid);
E2L4 = EL(:,1+wid:end);
Nomi = fft2(S) + ro1.*conj(otfFx).*fft2(C1L1) + ro1.*conj(otfFy).*fft2(C2L2) ...
+ ro2.*conj(otfFx).*fft2(E1L3) + ro2.*conj(otfFy).*fft2(E2L4);
Denomi = 1 + (ro1 + ro2) .* DxDy;
B_new = real(ifft2(Nomi./Denomi));
DiffB = [-imfilter(B_new,fx,'circular'),-imfilter(B_new,fy,'circular')];
%% for C1, C2, shrinkage
BL = DiffB - L1./ro1;
C_new = sign(BL) .* max(abs(BL) - lambda1./ro1 ,0);
%% for E1, E2
BL = DiffS - DiffB - L2./ro2;
E_new = BL;
temp = BL.^2;
t = temp < 2.*lambda2./ro2;
E_new(t) = 0;
%% for Li,i=1,2,3,4
L1_new = L1 + ro1 * (C_new - DiffB);
L2_new = L2 + ro2 * (E_new - DiffS + DiffB);
%% for ro
ro1 = ro1 *4;
ro2 = ro2 *4;
%% update
B = B_new;
C = C_new;
E = E_new;
L1 = L1_new;
L2 = L2_new;
end
D = S - B;
end%% This function performs the TV like decomposition
function [D,B] = L1(S,lambda)
iter = 10;
[hei,wid] = size(S);
fx = [1, -1];
fy = [1; -1];
otfFx = psf2otf(fx,[hei,wid]);
otfFy = psf2otf(fy,[hei,wid]);
DxDy = abs(otfFy).^2 + abs(otfFx).^2;
B = S;
C = zeros(hei,wid);
D = zeros(hei,wid);
T1 = zeros(hei,wid);
T2 = zeros(hei,wid);
ro = 1;
for i = 1:iter
%% for B
Nomi = fft2(S) + ro * conj(otfFx) .* fft2(C + T1./ro) + ro * conj(otfFy) .* fft2(D + T2./ro);
Denomi = 1 + ro * DxDy;
B_new = real(ifft2(Nomi./Denomi));
GradxB = -imfilter(B_new,fx,'circular');
GradyB = -imfilter(B_new,fy,'circular');
%% for C, D
BB1 = GradxB - T1./ro;
BB2 = GradyB - T2./ro;
C_new = sign(BB1) .* max(abs(BB1) - lambda.*1./ro ,0);
D_new = sign(BB2) .* max(abs(BB2) - lambda.*1./ro ,0);
%% for T1, T2
T1_new = T1 + ro * (C_new - GradxB);
T2_new = T2 + ro * (D_new - GradyB);
%% for ro
ro = ro * 2;
%% update
B = B_new;
C = C_new;
D = D_new;
T1 = T1_new;
T2 = T2_new;
end
D = S - B;
end
function y = compress(x,gamma,W)
if nargin<3
W = 1;
end
y = W * ((x./W) .^ (1./gamma));
endfunction y = clampp(x,a,b)
[d1,d2] = size(x);
low = round(a*d1*d2);
high = round(b*d1*d2);
so = sort(x(:));
low = so(low);
high = so(high);
x(x>high) = high;
x(x<low) = low;
y = x;
endfunction output = bilateralFilter( data, edge, edgeMin, edgeMax, sigmaSpatial, sigmaRange, ...
samplingSpatial, samplingRange )
if( ndims( data ) > 2 ),
error( 'data must be a greyscale image with size [ height, width ]' );
end
if( ~isa( data, 'double' ) ),
error( 'data must be of class "double"' );
end
if ~exist( 'edge', 'var' ),
edge = data;
elseif isempty( edge ),
edge = data;
end
if( ndims( edge ) > 2 ),
error( 'edge must be a greyscale image with size [ height, width ]' );
end
if( ~isa( edge, 'double' ) ),
error( 'edge must be of class "double"' );
end
inputHeight = size( data, 1 );
inputWidth = size( data, 2 );
if ~exist( 'edgeMin', 'var' ),
edgeMin = min( edge( : ) );
warning( 'edgeMin not set! Defaulting to: %f\n', edgeMin );
end
if ~exist( 'edgeMax', 'var' ),
edgeMax = max( edge( : ) );
warning( 'edgeMax not set! Defaulting to: %f\n', edgeMax );
end
edgeDelta = edgeMax - edgeMin;
if ~exist( 'sigmaSpatial', 'var' ),
sigmaSpatial = min( inputWidth, inputHeight ) / 16;
fprintf( 'Using default sigmaSpatial of: %f\n', sigmaSpatial );
end
if ~exist( 'sigmaRange', 'var' ),
sigmaRange = 0.1 * edgeDelta;
fprintf( 'Using default sigmaRange of: %f\n', sigmaRange );
end
if ~exist( 'samplingSpatial', 'var' ),
samplingSpatial = sigmaSpatial;
end
if ~exist( 'samplingRange', 'var' ),
samplingRange = sigmaRange;
end
if size( data ) ~= size( edge ),
error( 'data and edge must be of the same size' );
end
% parameters
derivedSigmaSpatial = sigmaSpatial / samplingSpatial;
derivedSigmaRange = sigmaRange / samplingRange;
paddingXY = floor( 2 * derivedSigmaSpatial ) + 1;
paddingZ = floor( 2 * derivedSigmaRange ) + 1;
% allocate 3D grid
downsampledWidth = floor( ( inputWidth - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
downsampledHeight = floor( ( inputHeight - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
downsampledDepth = floor( edgeDelta / samplingRange ) + 1 + 2 * paddingZ;
gridData = zeros( downsampledHeight, downsampledWidth, downsampledDepth );
gridWeights = zeros( downsampledHeight, downsampledWidth, downsampledDepth );
% compute downsampled indices
[ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 );
% ii =
% 0 0 0 0 0
% 1 1 1 1 1
% 2 2 2 2 2
% jj =
% 0 1 2 3 4
% 0 1 2 3 4
% 0 1 2 3 4
% so when iterating over ii( k ), jj( k )
% get: ( 0, 0 ), ( 1, 0 ), ( 2, 0 ), ... (down columns first)
di = round( ii / samplingSpatial ) + paddingXY + 1;
dj = round( jj / samplingSpatial ) + paddingXY + 1;
dz = round( ( edge - edgeMin ) / samplingRange ) + paddingZ + 1;
% perform scatter (there's probably a faster way than this)
% normally would do downsampledWeights( di, dj, dk ) = 1, but we have to
% perform a summation to do box downsampling
for k = 1 : numel( dz ),
dataZ = data( k ); % traverses the image column wise, same as di( k )
if ~isnan( dataZ ),
dik = di( k );
djk = dj( k );
dzk = dz( k );
gridData( dik, djk, dzk ) = gridData( dik, djk, dzk ) + dataZ;
gridWeights( dik, djk, dzk ) = gridWeights( dik, djk, dzk ) + 1;
end
end
% make gaussian kernel
kernelWidth = 2 * derivedSigmaSpatial + 1;
kernelHeight = kernelWidth;
kernelDepth = 2 * derivedSigmaRange + 1;
halfKernelWidth = floor( kernelWidth / 2 );
halfKernelHeight = floor( kernelHeight / 2 );
halfKernelDepth = floor( kernelDepth / 2 );
[gridX, gridY, gridZ] = meshgrid( 0 : kernelWidth - 1, 0 : kernelHeight - 1, 0 : kernelDepth - 1 );
gridX = gridX - halfKernelWidth;
gridY = gridY - halfKernelHeight;
gridZ = gridZ - halfKernelDepth;
gridRSquared = ( gridX .* gridX + gridY .* gridY ) / ( derivedSigmaSpatial * derivedSigmaSpatial ) + ( gridZ .* gridZ ) / ( derivedSigmaRange * derivedSigmaRange );
kernel = exp( -0.5 * gridRSquared );
% convolve
blurredGridData = convn( gridData, kernel, 'same' );
blurredGridWeights = convn( gridWeights, kernel, 'same' );
% divide
blurredGridWeights( blurredGridWeights == 0 ) = -2; % avoid divide by 0, won't read there anyway
normalizedBlurredGrid = blurredGridData ./ blurredGridWeights;
normalizedBlurredGrid( blurredGridWeights < -1 ) = 0; % put 0s where it's undefined
% for debugging
% blurredGridWeights( blurredGridWeights < -1 ) = 0; % put zeros back
% upsample
[ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 ); % meshgrid does x, then y, so output arguments need to be reversed
% no rounding
di = ( ii / samplingSpatial ) + paddingXY + 1;
dj = ( jj / samplingSpatial ) + paddingXY + 1;
dz = ( edge - edgeMin ) / samplingRange + paddingZ + 1;
% interpn takes rows, then cols, etc
% i.e. size(v,1), then size(v,2), ...
output = interpn( normalizedBlurredGrid, di, dj, dz );上面是所需的函数，下面是demo
clear;clc;
%% Parameters
lambda1 = 0.3;
lambda2 = lambda1*0.01;
lambda3 = 0.1;
%% read files
hdr = hdrimread('Street_Art.hdr');
[hei,wid,channel] = size(hdr);
tic;
%% transformation
hdr_h = rgb2hsv(hdr);
hdr_l = hdr_h(:,:,3);
hdr_l = log(hdr_l+0.0001);
hdr_l = nor(hdr_l);
%% decomposition
[D1,D2,B2] = Layer_decomp(hdr_l,lambda1,lambda2,lambda3);
%% Scaling
sigma_D1 = max(D1(:));
D1s = R_func(D1,0,sigma_D1,0.8,1);
% sigma_D2 = max(D2(:));
% D2s = R_func(D2,0,sigma_D2,0.9,1);
B2_n= compress(B2,2.2,1);
hdr_lnn = 0.8*B2_n + D2 + 1.2*D1s;
%% postprocessing
hdr_lnn = nor(clampp(hdr_lnn,0.005,0.995));
out_rgb = hsv2rgb((cat(3,hdr_h(:,:,1),hdr_h(:,:,2)*0.6,hdr_lnn)));
toc;
figure,imshow(out_rgb)
运行结果：
左边原图，右边TMO后图像
代码来自GitHub，我只是互联网的搬运工！
下一篇： 微信小程序 实现列表刷新的实例详解