【五仁科技】高动态范围（HDR）成像——debevec算法

首先简述问题和任务：正常情况下拍的一张照片无法同时将明亮的部分和阴暗的部分都显示清楚，我们需要将曝光时间长的照片（用于看清暗部）和曝光时间短的照片（用于看清亮部）进行融合，得到一张不论明亮都非常清晰的照片。

组合在不同曝光设置下获取的同一场景的不同图像的过程称为高动态范围 (HDR)成像。

前置知识

（1）我们平时所看到的图片文件，比如jpg、png那些，像素值是0~255的整数，我们称为 LDR（低动态范围）图片；而文件后缀名为 .hdr 的图片的像素值是小数，不是整数，在正常情况下是不能直接打开的，需要通过某个算法将小数的像素值转化为整数，才能在计算机中正常打开，这个算法称为色调映射算法。而某些能直接打开.hdr文件的软件，如PS（Photoshop），其实是有一个自己内置的色调映射算法。
（2）辐射强度的概念：单位时间内物体单位表面积辐射出某特定波长射线的能量，称为单色波长的辐射强度。常用单位为：W/m2，瓦特/平方米。以下将辐射强度称为辐照度。
（3）响应曲线的概念：当我们拍照的时候，现实场景中的光照会转化照片中的像素值，这可看作是一种映射关系，即函数。我们用字母表示物理量，辐照度为 E，曝光时间为 Δt，曝光为 X，那么照片的像素值（整数）对曝光X的响应是一个函数，将现实场景映射到照片中，这个函数就称为相机的响应函数（CRF），也可以叫做响应曲线。在曝光时间 Δt 不变的前提下，当某个光源变成两倍亮时（辐照度变为两倍），照片中的像素值往往不会变成原来的两倍，所以说响应曲线是一个非线性函数。

任务分解

我们的目的是将多张不同曝光时间的照片合成一张照片，合成之后的那张照片既能看清楚很亮的区域，也能看清楚很暗的区域。

我们的步骤一共有三步，请记住：
① 生成相机响应曲线
② 将多张照片融合，得到一个后缀为 .hdr 的图片文件，像素值为小数
③ 色调映射，将上一步得到的 hdr 图片转化成能查看的 jpg、png 图片文件

在计算机 opencv 库中，上述的三个步骤可以对应三个函数。
需要注意的是，由于计算机的函数可以重载，所以同名函数可以传入不同的参数。

为什么要说这个呢，因为上述的第 ② 步可以不依赖于第 ① 步。正常情况下，确实应该先求出响应曲线，然后第 ② 步利用这个响应曲线来融合照片。但是，如果没有求出第 ① 步的响应曲线，也不妨碍第 ② 步的执行，前面提到，真实的响应曲线一般是非线性的，换而言之，我们可以假设一个最理想但不真实的响应曲线，他是线性的，然后用这个线性响应曲线来融合图片。

简单的说，就是第 ② 步的融合函数有两个重载，一个需要传入响应曲线后进行融合，这个响应曲线是第 ① 步生成的；另一个重载不需要传入响应曲线，它直接使用一个固定的线性响应曲线，当然，这种固定的线性响应曲线一般不符合真实场景，所以融合出的图片也容易不符合预期。

Debevec算法

下面我们来介绍一下 Debevec 算法：
这个算法包含两个部分，对应第 ① 步和第 ② 步，即生成响应曲线和生成 hdr 图片。

如何求相机响应曲线。

已有的材料是几张有亮有暗的照片，已知是这几张照片的各个像素的像素值，以及这几张照片的曝光时间 Δt。
各符号的意义：曝光 X，辐照度 E，曝光时间 Δt，像素值 Z
前提假设：X=E·Δt，请务必记住这个假设，它是Debevec算法成功的前提，以后会介绍当这个假设不成立时，Debevec 算法会如何失败。
X=E·Δt 可以理解为，当辐照度 E 固定时，曝光度 X 随着曝光时间 Δt 线性增大，从人体直观的感受上来看，就是曝光时间越久的照片越亮。事实上，严格满足 X=E·Δt 几乎是不可能的，但是，越接近这个式子，Debevec 算法的效果就越好。

前面说过，相机响应曲线是曝光 X 到像素值 Z 的函数，我们设为 f，即 $Z = f(X) = f(E·Δt)$，我们对应到具体的某一张照片上，即 $Z_{ij} = f(E_i·Δt_j)$，其中 $Z_{ij}$ 为第 j 张照片的第 i 个像素，$E_i$ 为照片中第 i 个像素所对应的位置在现实世界中的辐照度，$Δt_j$ 为第 j 张照片的曝光时间。

这里可能会有疑问，一张照片不是有长和宽吗，这明明需要 2 个坐标才能表示一个像素的位置，为什么下标 i 可以表示照片某个像素的位置呢？虽然图片是二维像素矩阵，但实际上可以压成一维，比如一张分辨率 2×5 的图片，实际上也可以表示为一个 10 维的向量，这样我们就能用一个下标来表示任意位置的像素了。第 6 个像素即对应于分辨率 2×5 图片中第 2 行第 1 列的像素。

公式推演

下面我们对 $Z_{ij} = f(E_i·Δt_j)$ 这条式子作变换，首先是初始公式：

$$Z_{ij} = f(E_i·Δt_j) \tag{1}$$

因为我们假设 f 是单调的（环境越亮，照片就越亮），则它是可逆的，我们可以把 (1) 重写为：

$$f^{-1}(Z_{ij}) = E_i·Δt_j$$

两边同时取对数，得：

$$ \ln{f^{-1}(Z_{ij})} = \ln{E_i} + \ln{Δt_j}$$

我们定义复合函数 $g(x) = \ln{f^{-1}(x)}$，得：

$$ g(Z_{ij}) = \ln{E_i} + \ln{Δt_j} \tag{2}$$

我们由（1）式得到了（2）式，这就是我们的核心公式，我们将以（2）式为基础来建立模型。

针对（2）式出现的字母，回想我们已知和未知的数据：各个点的像素值 $Z_{ij}$ 和各张图片的曝光时间 $Δt_j$已知，各个位置的辐照度 $E_i$ 和函数 g 未知。在这里再次明确我们的任务是求响应函数 f，注意 f 的值域是离散的，只有 0~255 的整数，相对应的，函数 g 的定义域也是只有 0~255 的整数，所以我们只需要求出函数 g 的 256 个点，即 $\{~g(0),g(1),\cdots,g(255)~\}$，就可以大致绘出响应曲线（小学生也懂的描点绘图法）。

于是，我们把每个像素点对应的 $Z_{ij}$ 值和 $Δt_j$ 值代入（2）式，能得到多条方程，他们组成一个方程组，其中未知数为 $E_i$ 和 $\{~g(0),g(1),\cdots,g(255)~\}$，它们也是我们所要求出的方程组的解。显然，只要方程的数量等于未知数的数量，这个方程组就有唯一解。

幸运的是，方程的数量往往是过多的，比如我们要合成 2 张 1920×1080 的照片，那么每张照片就有 2073600 像素，i 的范围是 1~2073600，j 的范围是 1~2，相当于用 4147200 个方程解 2073856（2073600+256）个未知数，这种方程个数大于未知数个数的方程组叫做超定方程组，一般是不存在精确解的，但是超定方程组一定存在最小二乘解，这也是机器学习中在拟合时常用的技巧。如果还记得线性代数中奇异值分解（SVD）的知识，就能知道可以用 SVD 求伪逆的方法来求解最小二乘法问题。

由最小二乘的知识可以知道，用上述方法求出的解使得下列式子的 O 值最小：

$$O = \sum_{i=1}^{N}\sum_{j=1}^{P}[g(Z_{ij}) - \ln{E_i} - \ln{Δt_j}]^2 \tag{3}$$

其中 P 是图片数，N 是一张图片的像素数。在机器学习中我们可以把 O 值理解为损失函数。

但是仅此而已还体现不出算法的精妙性，所以Debevec算法对上式做了以下变形：

$$ O = \sum_{i=1}^{N}\sum_{j=1}^{P}[g(Z_{ij}) - \ln{E_i} - \ln{Δt_j}]^2+\lambda\sum_{z=Z_{min}+1}^{Z_{max}-1}g''(z)^2 \tag{4}\\ g''(z) = g(z-1) - 2g(z) + g(z+1) $$

注意到，相比于（3）式，（4）式后面加多了一项求和的东西，这一项有什么作用呢。

原来这一项是平滑项，之前的 O 值的表达式过于理想，没有考虑到现实环境中的噪声，噪声指的是环境中一些奇奇怪怪的因素所产生的干扰。如果不加平滑项，那么求出来的响应曲线可能不“丝滑”，可能连续的几个 g 值上下反复跳动，这显然是我们不想要见到的（我们希望它单调）。所以加入了平滑项，它令相邻的三个g值不会相差太远，（4）式中把进行 $g''(z)^2$ 求和，这样做的话，用 SVD 解出来的结果会自动地尽量符合 $g(Z_{ij}) = \ln{E_i} + \ln{Δt_j}$ ，同时使响应曲线尽量平滑，（4）式子中 $g''(z)^2$ 求和前面的系数 λ 是一个人为设置的超参数，可以理解为平滑度常数，需要根据环境中的噪声进行人为手动选择。~~算法总得掺点玄学，没点人工调参都不好意思发论文。~~

下面继续变。我们来考虑一下，是否照片中的所有像素都同等重要？什么意思呢，考虑照片中有一部分是非常黑或者非常白，它能给予我们的信息其实是非常有限的，当我们观赏一张照片时，我们往往是盯着那些既不太亮又不太暗的部分，那些区域能给予我们更多的信息。从另一角度来说，因为像素值范围是有限的，当它取边界值时，往往不能反映真实世界中的光照，于是我们可以认为这些太黑或太白的像素是不那么重要的。这里需要区分绘画和照相，绘画中，使用纯黑或纯白来增强表现力是常用技巧，而不同于照相。那么就有一个很简单的想法，我们给每个像素加权，让它体现不同的重要性即可。

所以，我们对 0~255 这 256 个像素值赋以不同的权重，越靠中间（127）的权重越高，越靠两边（0 或 255）的权重越低。那么 Debevec 算法的加权函数如下：

$$ w(z)= \begin{cases} z - Z_{min} & for \quad z\le\frac{1}{2}(Z_{min} + Z_{max}) \\ Z_{max} - z & for \quad z\gt\frac{1}{2}(Z_{min} + Z_{max}) \end{cases} $$

上式中，$z$ 是某一个像素值，$w(z)$ 是这个像素值的权重，$Z_{min}$ 和 $Z_{max}$ 是照片中所有像素的最小像素值和最大像素值，直接取 0 和 255 也可以。将权重函数加入到之前的最小二乘 O 值表达式就可以得到最终的式子了：

$$ O = \sum_{i=1}^{N}\sum_{j=1}^{P}\{w(Z_{ij})[g(Z_{ij}) - \ln{E_i} - \ln{Δt_j}]\}^2+\lambda\sum_{z=Z_{min}+1}^{Z_{max}-1}[w(z)g''(z)]^2 \tag{5} $$

我们要设法使得 O 值最小。下面将介绍如何用线性代数的知识来实现，方法是解线性方程组 $A\vec{x}=\vec{b}$ 。重申一遍，未知数为 $E_i$ 和 $\{~g(0),g(1),\cdots,g(255)~\}$，它们也是我们所要求出的方程组的解。因此，我们要想办法构造矩阵 A 和向量 $\vec{b}$ 。

构造线性方程组

显然，未知数要放在 $\vec{x}$ 向量里面，所以 $\vec{x}=[g(0)~g(1) \cdots g(255)~\ln{E_1}~\ln{E_2} \cdots \ln{E_n}]$，而方程是 $$w(Z_{ij})g(Z_{ij}) - w(Z_{ij})\ln{E_i} = w(Z_{ij})\ln{Δt_j}\tag{a}$$，该方程对应（5）式中的 $\sum_{i=1}^{N}\sum_{j=1}^{P}\{w(Z_{ij})[g(Z_{ij}) - \ln{E_i} - \ln{Δt_j}]\}^2$ 部分，也就是加权后的（2）式。未知数和它的系数已经放在了等式的左边，如果学过线性代数的小伙伴应该知道，我们很轻松就能使用一个像素构造矩阵 A 的一行。

假设当我们使用一张曝光时间为 0.2 秒的照片的第一个像素来构造方程时，假设它的像素值是6，则代入后的方程就是：

$$w(6)g(6) - w(6)\ln{E_i} = w(6)\ln{0.2}$$

将其具体展开就是：

$$0·g(0)+\cdots+w(6)g(6)+\cdots+0·g(255) - w(6)\ln{E_1} - 0·\ln{E_2}- \cdots -0·\ln{E_2} = w(6)\ln{0.2}$$

上面的线性方程放入增广矩阵 A|b 的一行就是：

$$[0\cdots0~w(6)~0\cdots0~-\!\!w(6)~0\cdots0~|~w(6)\!\ln{0.2}]$$

也就是说，每张照片的每个像素都可以构造矩阵 A|b 的一行，每一行只需要修改 3 个位置的值。以上面的例子来说，就是把 $w(6)$ 填入下标为 6 的位置（对应$g(6)$），把 $-w(6)$ 填入下标为 256 的位置（对应$\ln{E_1}$），并把 $w(6)\ln{0.2}$ 填入最后一列（对应向量$\vec{b}$）。其余的像素都仿照这样的方式填入矩阵 A|b 即可。

（a）方程对应了（5）式的前半部分，我们还需要构造一些方程对应（5）式的后半部分，也就是平滑项。

$$\lambda·w(z)g(z-1) - 2\lambda·w(z)g(z) + \lambda·w(z)g(z+1)=0\tag{b}$$

（b）方程对应了（5）式中的 $\sum_{z=Z_{min}+1}^{Z_{max}-1}[\lambda w(z)g''(z)]^2$ 部分，你可能注意到我把 $\lambda$ 乘进去了，也就是把 $\lambda$ 变成了 $\lambda^2$，这其实没所谓，因为它只是一个人为决定的超参数，平不平方也就是影响我们手动调参的过程而已，只要你喜欢，也可以把（b）方程中的 $\lambda$ 开平方。

至于（b）方程如何插入矩阵 A，其实跟前面的（a）方程是一样的，我们只需要在下标分别为 z-1、z、z+1 的地方填上系数 $\lambda·w(z)$、$-2\lambda·w(z)$、$\lambda·w(z)$即可。用 python 代码来构造 A|b 的话差不多就是下面这样（还不是可运行代码）：

k = 0
for i in range(p_num):
    for j in range(img_num):
        zij = Z[j][(int)((int)(pixel_num / p_num) * i / col)][(int)(pixel_num / p_num) * i % col][channel]
        wij = w(zij + 1, Zmin, Zmax)
        A[k][zij] = wij
        A[k][n + i] = -wij
        B[k][0] = wij * exposure_times[j]
        k = k + 1

A[k][(int)((Zmin + Zmax) / 2)] = 1
k = k + 1

for i in range(Zmin + 1, Zmax):
    A[k][i - 1] = l * w(i + 1, Zmin, Zmax)
    A[k][i] = -2 * l * w(i + 1, Zmin, Zmax)
    A[k][i + 1] = l * w(i + 1, Zmin, Zmax)
    k = k + 1

至此，我们就构造好了 $A\vec{x}=\vec{b}$ 中的矩阵 A 和向量 $\vec{b}$，而向量 $\vec{x}$ 是我们所求。一般情况下，矩阵 A 的行数远远大于其列数，也就是我们之前所说的超定方程组。使用线性代数的知识可以轻松解出 $\vec{x}$，并保证 $\vec{x}$ 是满足（5）式的最小二乘解，我们只需要取出 $\vec{x}$ 前 256 个数，即 $g(0)$ 到 $g(255)$，就能得到相机响应曲线了。

如何利用相机响应曲线将多张图片融合成一张hdr图片

我们上一步求出了函数 g 的 256 个值 g(0)~g(255)，他满足之前的（2）式：

$$ g(Z_{ij}) = \ln{E_i} + \ln{Δt_j} \tag{2}$$

移个项，即

$$\ln{E_i} = g(Z_{ij}) - \ln{Δt_j} $$

解开 $\ln$，即

$$E_i = exp(g(Z_{ij}) - \ln{Δt_j}) $$

将每张照片的每个像素代入上式，我们可以求出某一个像素位置的辐照度 $E_i$，这是一个小数，我们将辐照度 $E_i$ 作为最终生成的 hdr 图片的像素值即可。具体的方法其实也是加权平均：

$$\ln{E_i}=\frac{\sum_{j=1}^{P}w(Z_{ij})(g(Z_{ij})-\ln{\Delta t_j})}{\sum_{j=1}^{P}w(Z_{ij})}$$

其中 P 是照片的数量。也就是说，针对某一个位置的像素点，将不同曝光图片的这个位置的辐照度进行加权平均，求出最终 hdr 图片的这个位置的像素值。

你也许有疑问，我们之前求 $\vec{x}$ 时不是把每个点的 $E_i$ 都求出来了吗，为什么现在又重新求一遍呢？原因之一是：我们往往不会将所有的像素点用于构造 $A\vec{x}=\vec{b}$ 这个线性方程组，毕竟一般情况下只要矩阵 A 的行数大于列数就存在最小二乘解，而求解的时间复杂度是很高的。如今随随便便一张照片都不止 100 万像素，几张图片的像素总量可能上千万，假如所有像素都用于构造 $A\vec{x}=\vec{b}$，解这个方程组可能得花上一年半载。实际操作时，我们往往是只选其中若干个位置的像素来构造矩阵，这样就减少了 $\vec{x}$ 的维度，同时也降低了矩阵 A 的行数和列数，降到能在短时间内求出解的程度，所以我们还要多花一部来求每个点的 $E_i$ 。那些选取的像素点可能是随机选的，也可能是有规律的（比如每间隔多少就选一个），甚至可能是人为挑选的，该算法的原论文就是人手动选点（搁这学术造假呢），选点的好坏也会影响相机响应曲线的好坏，从而影响最终图片的效果，这也是玄学之一。

到这里，我们已经介绍完了Debevec算法的两个最精髓的部分——求响应曲线和融合图片，请容许我再说一遍合成图片的步骤：
① 求响应曲线
② 融合出 HDR 图片
③ 色调映射，将 HDR 图片的像素值转化成 0~255 的整数

显然，我们现在剩下最后一个步骤：色调映射。遗憾的是，Debevec算法并不包括色调映射有关的内容，我们演示时，将直接调用色调映射函数的接口，而暂时不关注色调映射算法具体的实现。而且，色调映射实际上是一个不错的研究方向，有许多人在这领域中研究，所以也存在多种色调映射算法。

以下代码就是一个完整的高动态范围成像示例：

import numpy as np
import cv2

np.set_printoptions(suppress=True)

# 读图
img_fn = ["img_0.033.jpg", "img_0.25.jpg", "img_2.5.jpg", "img_15.jpg"]
img_list = [cv2.imread(fn) for fn in img_fn]
exposure_times = np.array([0.033, 0.25, 2.5, 15.0], dtype=np.float32)

# 对齐
alignMTB = cv2.createAlignMTB()
alignMTB.process(img_list, img_list)

# 计算响应曲线
calibrateDebevec = cv2.createCalibrateDebevec()
responseDebevec = calibrateDebevec.process(img_list, exposure_times)
print(type(responseDebevec),responseDebevec.shape)
print(responseDebevec)

# 融合
mergeDebevec = cv2.createMergeDebevec()
hdrDebevec = mergeDebevec.process(img_list, exposure_times, responseDebevec)
cv2.imwrite("hdrDebevec.hdr", hdrDebevec)

# 色调映射
tonemapMantiuk = cv2.createTonemapMantiuk(2.2, 0.85, 1.2)
ldrMantiuk = tonemapMantiuk.process(hdrDebevec)
ldrMantiuk = 3 * ldrMantiuk
cv2.imwrite("ldr-Mantiuk.jpg", ldrMantiuk * 255)

代码执行结果如下，把多张明暗不一的照片组合成一张照片：