HDR：Recovering High Dynamic Range Radiance Maps from Photographs_综合

Recovering High Dynamic Range Radiance Maps from Photographs

本文以胶片为例，提出了通过多帧非线性RGB图像重建线性的高动态范围图像(High Dynamic Range Image)，之后在这篇文章的基础上延伸出来众多标定数码相机的相机响应曲线(Camera Response Function)来重建线性的高动态范围图像方法，本文算法算是计算相机响应曲线的开山之作。

文章目录

Recovering High Dynamic Range Radiance Maps from Photographs
- 几个概念
- - 传感器的非线性响应(Non-Linear Photosite Response)
  - 相机响应曲线(Camera Response Function)
  - 胶片的曝光响应曲线(Huter and Driffield Curves)
- Film Response Recovery
- Constructing the High Dynamic Range Radiance Map

几个概念

传感器的非线性响应(Non-Linear Photosite Response)

图像传感器的光电转换是线性过程，一般而言，图像传感器中势肼收集的光生电荷的数量与经过A/D转换后的数字信号之间也呈线性关系，但传感器的动态范围(Dynamic Range)和噪声(Noise)会造成光生电荷与数字信号间的非线性

动态范围对传感器响应的影响

文章CMOS图像传感器中的噪声来源分析中分析了图像传感器本身的满肼容量(Full-Well Capacity)决定了图像传感器动态范围(Dynamic Range)的上限，当光生电荷量超过满肼容量时，会造成过曝(Blooming Effect)：
$y^=min(y,dynamicrange)\hat{y}=min(y,dynamic \ range)$

满肼容量对势肼中光生电荷数量的影响

HDR：Recovering High Dynamic Range Radiance Maps from Photographs

常见sensor的满肼容量示意图

噪声对图像传感器响应的影响

文章噪声噪声分析中分析了图像传感器中噪声的呈现泊松-高斯分布，噪声的跳变会带来传感器响应的非线性，尤其是在接近0或动态范围上限的区域的非线性，具体内容可以参见Foi的论文《Clipped noisy images: Heteroskedastic modeling and practical denoising》，简单来说就是当势肼中的光生电荷数量接近0或动态范围上限时，由于存在均值为0的噪声，使得接近0的点由于负向噪声的影响，被Clipping至0，接近动态范围上限的点由于正向噪声的影响，被Clipping至动态范围上限，从而造成传感器响应的非线性：
$y^=max(min(y+noise,dynamicrange),0)\hat{y}=max(min(y+noise,dynamic \ range),0)$
其中 $y$ 为无噪声影响的响应值， $y^\hat{y}$ 为受噪声影响之后的响应值，具体非线性响应曲线如下：

横坐标为无噪声影响的响应值，纵坐标为受噪声影响的响应值

和满肼容量造成传感器响应的非线性不同，由噪声造成的传感器响应的非线性可以通过对噪声的建模进行逆变换，以消除造成带来的非线性，但整个逆变换过程相对比较复杂，后续ISP(Image Singal Processing)一般不考虑消除这一部分的非线性来提升后续图像处理的质量。

相机响应曲线(Camera Response Function)

参见文章ISP Pipeline对ISP当中图像处理模块的介绍可知，通过图像传感器采集下来的原始图像信号(Image Raw Data)，在后续的图像处理中非进过多次非线性处理，比如Gamma校正、局部对比度增强等，这些非线性处理，造成最终输出的标准RGB图像(Standard RGB)与原始Raw数据间的非线性，即相机响应曲线即是图像传感器采集下来的原始Raw数据与经过一系列非线性处理得到的sRGB图像间全部亮度变换的总和。

HDR：Recovering High Dynamic Range Radiance Maps from Photographs

高通 MSM8998 ISP Pipeline

胶片的曝光响应曲线(Huter and Driffield Curves)

也称Characteristic Film Curve，是胶片的曝光量与胶片沉积银颗粒的密度的关系。首先简单介绍胶片的构成，如下图所示，胶片由保护膜、感光层(主要是卤化银和明胶)、片基(支撑作用，常用玻璃)、防光晕层(一部分透过感光层的光线会被片基或相机内部反射回胶片，如果被感光层再次吸收会在光源周围的位置产生光晕，防光晕层就是为了吸收这些反射光线)。

HDR：Recovering High Dynamic Range Radiance Maps from Photographs

胶片结构示意图

当光线通过感光层，胶片上的辐照度为 $E$ ，在曝光时间 $Δt\Delta{t}$ 内累计曝光量为 $X=EΔtX=E\Delta t$ ， $X$ 与最终胶片上沉积的银粒子的质量或者密度间的关系即为胶片的曝光响应曲线，可以参考Kodak Pro 800的曝光响应曲线，下图是不同胶片的曝光响应曲线，其中Fog是由自然界的背景辐射造成，Base是由胶片的自然衰减造成。另外，也常使用胶片的显影率曲线来作为胶片的曝光响应曲线。

HDR：Recovering High Dynamic Range Radiance Maps from Photographs

不同胶片的曝光响应曲线

HDR：Recovering High Dynamic Range Radiance Maps from Photographs

胶片的显影率

可以观察到在胶片的动态范围当中大部分响应是接近线性的，因此同时假设胶片机与数码相机采集到的原始图像数据与曝光量之间近似是线性关系。

Film Response Recovery

HDR：Recovering High Dynamic Range Radiance Maps from Photographs

图像采集过程的常规处理流程

如上述，除欠曝和过曝区域，胶片曝光后产生的银与曝光量 $X$ 间的关系为： $D∝ln(X)=ln(EΔt)D\propto ln (X)=ln(E \Delta t)$ ，经过一列处理后得到最终的数字图像 $Z$ 和曝光量间不再是简单的对数关系，可以用一种更复杂的非线性函数 $f$ 来表示，即：
$Zij=f(X)=f(EiΔtj)(1)Z_{ij} = f(X)=f(E_{i}\Delta t_{j})\tag{1}$
即有：
$f?1(Zij)=EiΔtj(2)f^{-1}(Z_{ij})=E_{i}\Delta t_{j} \tag{2}$
假设在整个曝光过程中，场景的辐照度 $E$ 恒定，则 $E=f?1(Z)ΔtE=\frac{f^{-1}(Z)}{\Delta t}$ ，对上式(2)取对数：
$lnf?1(Zij)=lnEi+lnΔtjg=:lnf?1?g(Zij)=lnEi+lnΔtj(3)\begin{aligned} &lnf^{-1}(Z_{ij})=ln E_{i}+ln\Delta t_{j}\\ &g=:lnf^{-1} \longrightarrow g(Z_{ij})=lnE_{i}+ln\Delta t_{j}\tag{3} \end{aligned}$
其中 $i$ 对应像素位置索引， $j$ 对应曝光时间索引，上式中 $Zij、ΔtjZ_{ij}、\Delta t_{j}$ 已知，一般 $f$ 是单调递增的，即 $g$ 也是单调递增的，值得注意的是，由于 $Z$ 的取值有限，一般为0~255的整数，因此恢复函数 $g$ 仅需要恢复 $g (Z)$ 的有限个数量的值及 $N$ 个 $lnE_{i}$ 的值，为使方程超定，方便求解最小二乘解，需要满足 $N(P-1)＞(Z_{max}-Z_{min})$ ,其中 $P$ 为图像张数，即有损失函数：
$O=∑i=1N∑j=1P[g(Zij)?ln?Ei?ln?Δtj]2(4)O=\sum_{i=1}^{N} \sum_{j=1}^{P}\left[g\left(Z_{i j}\right)-\ln E_{i}-\ln \Delta t_{j}\right]^{2}\tag{4}$
最后Debevec再次假设响应曲线 $g$ 应该是一条光滑的曲线，添加了函数 $g$ 的二阶导作为正则化项，以约束 $g$ 保持平滑：
$O=∑i=1N∑j=1P[g(Zij)?ln?Ei?ln?Δtj]2+λ∑Z=Zmin?+1Zmax??1g′′(z)2(5)O=\sum_{i=1}^{N} \sum_{j=1}^{P}\left[g\left(Z_{i j}\right)-\ln E_{i}-\ln \Delta t_{j}\right]^{2}+\lambda \sum_{Z=Z_{\min }+1}^{Z_{\max }-1} g^{\prime \prime}(z)^{2}\tag{5}$
第二项用于约束函数 $g$ 的平滑性，再次观察前面给出的响应曲线，通常在 $Z_{min}、Z_{max}$ 附近有较高的二阶导，中间部分接近线性，通过引入权重 $w$ 来更好拟合中部位置的曲线，作者选择的权重函数如下：
$z>12(Zmin?+Zmax?)(6)w(z)=\left\{\begin{array}{ll} z-Z_{\min } & \text { for } z \leq \frac{1}{2}\left(Z_{\min }+Z_{\max }\right) \\ Z_{\max }-z & \text { for } z>\frac{1}{2}\left(Z_{\min }+Z_{\max }\right) \end{array}\right.\tag{6}$
最终需要优化求解的目标函数为：
$O=∑i=1N∑j=1P{w(Zij)[g(Zij)?ln?Ei?ln?Δtj]}2+λ∑Z=Zmin?+1Zmax??1[w(z)g′′(z)]2(7)O=\sum_{i=1}^{N} \sum_{j=1}^{P}\left\{w\left(Z_{i j}\right)\left[g\left(Z_{i j}\right)-\ln E_{i}-\ln \Delta t_{j}\right]\right\}^{2}+\lambda \sum_{Z=Z_{\min }+1}^{Z_{\max }-1}\left[w(z) g^{\prime \prime}(z)\right]^{2}\tag{7}$
HDR：Recovering High Dynamic Range Radiance Maps from Photographs

低动态范围的图像序列，曝光时间在1/1000~30秒

HDR：Recovering High Dynamic Range Radiance Maps from Photographs

由此方法恢复出的DCS460响应曲线

Constructing the High Dynamic Range Radiance Map

求解出相机响应曲线之后，根据式(3)便可求得像素 $i$ 对应的辐照度 $E_{i}$ ：
$lnEi=g(Zij)?lnΔtjlnE_{i}=g(Z_{ij})-ln\Delta t_{j}$
为了增加所求辐照度 $E_{i}$ 的鲁棒性，用 $P$ 帧图像的加权平均，作为像素点 $i$ 对应场景的最终辐照度 $E_{i}$ ：
$ln?Ei=∑j=1Pw(Zij)(g(Zij)?ln?Δtj)∑j=1Pw(Zij)(8)\ln E_{i}=\frac{\sum_{j=1}^{P} w\left(Z_{i j}\right)\left(g\left(Z_{i j}\right)-\ln \Delta t_{j}\right)}{\sum_{j=1}^{P} w\left(Z_{i j}\right)}\tag{8}$
其中 $w$ 为与像素值 $Z$ 相关的权重函数，可以选择式(6)所示函数，也可以选择其他权重函数，最终的目的是赋予靠近中灰度值的像素点更高的权重，靠近 $Z_{max}、Z_{min}$ 的像素点更低的权重，由式(8)可从不同曝光的图像序列中重建线性的高动态范围图像。最后用直方图均衡化对重建得到的高动态范围图像进行动态范围压缩，以符合人眼视觉对辐照度 $E$ 的反应，也便于后续编码、存储等。

HDR：Recovering High Dynamic Range Radiance Maps from Photographs
上图中：a:曝光合适的低动态范围图像；b:本方法重建后的高动态范围图像并现在缩放至低动态范围来显示；c本方法重建后的高动态范围图像并当中低动态范围图像直接显示；d:低动态范围图像序列的灰度图组成的伪彩色图；e:本方法重建后的高动态范围图像用直方图均衡化做动态范围压缩后的显示；f:其他方法生成的高动态范围图像