一种基于论文Radiometric Self Calibration的相机响应函数计算方法

本文是Jesse Chen的原创文章。 在研究HDR技术时，会有一个比较成熟的HDR分支领域—相机响应函数的计算。作者在通读这个领域的经典论文时，注意到了经典论文Radiometric Self Calibration中给出的计算方法存在的计算不单调和不合理问题。针对这些问题，本博文的作者提出了新的计算方案，并编写程序，进行的仿真。

一种基于论文Radiometric Self Calibration的相机响应函数计算方法 论文Radiometric Self Calibration 的回顾Radiometric Self Calibration 论文算法流程Radiometric Self Calibration 论文算法改进

---

论文Radiometric Self Calibration 的回顾

我们都知道可通过一个响应函数 M=g(I) 描述成像系统最终产生的亮度观测值 M 和（经缩放的）环境光辐射I之间的关系。该响应函数一般都是单调递增的，因此我们可以用函数 g 的逆函数f=g−1将所有的测量值M映射到环境光辐射值 I ，I=f(M)。

论文Radiometric Self Calibration描述了如何从一组同一场景的不同曝光图像中找到函数 f 。

值得注意的是，一个图像系统和另一个图像系统的响应曲线会有明显的差异，但是响应曲线不具有奇异的形式。首要原因是，图像传感器的输出是单调或半单调的。除非在图像系统中剩余模块中有不同寻常的映射，随着环境光辐射（曝光）的增加测量到的亮度要么也增加或者维持不变。因此，论文作者认为几乎所有的响应函数都可以用高阶多项式建模。

I=f(M)=∑n=0NcnMn(1)

因此只要确定了上述多项式的最高阶数 N 以及对应的多项式系数cn就可以确定响应函数。

假定同一场景，分别用不同的曝光 eq 和 eq+1 得到了两张照片，设曝光比 Rq,q+1 为

Rq,q+1=eqeq+1(2)

因此图像系统的响应函数和曝光比之间满足：

f(Mp,q)f(Mp,q+1)=Rq,q+1(3)

对图像进行排序以满足 eq<eq+1 ，因此 0<Rq,q+1<1 。将多项式模式带入上式，我们得到

∑Nn=0cnMnp,q∑Nn=0cnMnp,q+1=Rq,q+1(4)

假定曝光比 Rq,q+1 是已知的，那么响应函数可以通过下面的最小二乘误差计算得到。

\varepsilon=\sum_{q=1}^{Q-1}\sum_{p=1}^{P}\left[\sum_{n=0}^{N}c_{n}M_{p,q}^{n}-R_{q,q+1}\sum_{n=0}^{N}c_{n}M_{p,q+1}^{n}\right]^{2}\label{a}\tag{5}

上式中， Q 是计算中使用的总图像数，p表示不同的像素位置， P 表示像素位置的总数。假设，作者对像素值M进行归一化处理，那么 0≤M≤1 。规定环境光辐射的缩放比例，使得满足 f(1)=Imax ，作者得到附加约束：

cN=Imax−∑n=0N−1cn(6)

对 ε 求 cn 的偏导数，并让这些偏导数等于0。可以得到求解 cn 的线性方程组。

∂ε∂cn=0(7)

论文作者认为对于大多数的图像系统，难以获得准确的曝光比。论文作者提出通过迭的方法去估计多次成像的曝光比。当前的曝光比估计值 R(k−1)q,q+1 被用于计算下一组多项式系数 c(k)n 。这组多项式系数用来进一步更新曝光比：

R_{q,q+1}^{(k)}=\sum_{p=1}^{P}\frac{\sum_{n=0}^{N}c_{n}^{(k)}M_{p,q}^n}{\sum_{n=0}^{N}c_{n}^{(k)}M_{p,q+1}^n}\label{b}\tag{8}

初始的曝光比 R(0)q,q+1 可以直接由照片的EXIF信息得到。算法的收敛条件是：

∣∣f(k)(M)−f(k−1)(M)∣∣<ϵ,∀M(9)

其中 ϵ 是一个指定的常数。

论文作者建议多项式中 N 最大取10，在计算的过程中要取使得平方误差最小的N。

Radiometric Self Calibration 论文算法流程

从Radiometric Self Calibration的描述中，我们不难得到如下的计算步骤：

Step 1: 对曝光值不同的图像集( Mq,q=1⋯Q )进行必要的预处理。 Step 2: 根据EXIF信息计算曝光值( eq,q=1⋯Q )，从而得到初始的曝光比 R(0)q,q+1=eqeq+1，q=1⋯Q−1 。 Step 3: 设定N=2以及多项式距离 dthershold ( dthershold>0 )。 Step 4: 在当前的 N 下，初始化最大的内层迭代次数MaxIter。初始化两个多项式 f=0,fprev=0 ，以及迭代计数值 Iter 。 Step 5: 计算方程 Ax=b 中的矩阵 A 和向量b。根据 ∂ε∂ci=0 ，可以得到

∑q=1Q−1∑p=1P{(Mip,q−R(k−1)q,q+1Mip,q+1)⋅[∑n=0N−1cn(Mnp,q−R(k−1)q,q+1Mnp,q+1)+(Imax−∑n=0N−1cn)(MNp,q−R(k−1)q,q+1MNp,q+1)]}=0(10) 从上式不难得到 A 和b的计算方法。 Step 6: 解方程 Ax=b 得到多项式系数，从而获得本次迭代计算得到的多项式函数 f 。迭代计数Iter=Iter 1。 Step 7: 计算 v=[f(0),f(1255),f(2255),⋯,f(254255),f(1)]T vprev=[fprev(0),fprev(1255),fprev(2255),⋯,fprev(254255),fprev(1)]T d=∥∥v−vprev∥∥∞ if Iter<MaxIter && d>dthreshold ，

使用公式 (8) 更新曝光比

fprev=f

goto Step 5

else

计算当前 N 下的cost值 N=N 1if ( N≤8 ) goto Step 4 else goto Step 8.

end Step 8: 比较不同多项式阶数 N 下的cost值，选择cost最小的多项式作为最终的计算结果。

这里博文作者附上上述算法的核心代码段：

A = zeros(N, N); b = zeros(N, 1); Mmax = 1.0; for q=1:(Q - 1) Rq = R(q); M_q = M(:, q , c); M_q_p = M(:, q 1, c); indx = find(M_q > 0.0 & M_q_p > 0.0); e = M_q(indx).^N - (M_q_p(indx).^N)*Rq; for i=0:N-1 di = M_q(indx).^i - (M_q_p(indx).^i)*Rq; di = di - e; %init b b(i 1) = b(i 1) - sum((di.*e)*Mmax); %init A for j=0:N-1 A(i 1,j 1) = A(i 1,j 1) sum(di.*(M_q(indx).^j - (M_q_p(indx).^j)*Rq - e)); end end end coeff = A \ b; coeff_n = Mmax - sum(coeff);

Radiometric Self Calibration 论文算法改进

博文作者自己编程实现了Radiometric Self Calibration 论文中所述的算法，但是在编程过程中发现该算法计算的中间结果和最终计算结果发现两大问题： - 某些Radiometric值小于0，显然这是不合理的 - 计算出来的曲线不满足单调性 如下面的结果就是一个不合理的例子。

问题出在什么地方呢？ 博文作者认为原始论文没有考虑到求解最优值得时候，需要对优化进行约束。因此，博文作者基于Mitsunaga Nayar论文，对论文中的算法进行了小修改。

采用带约束的优化方法进行计算。采用这种方案后，无需针对不同的多项式系数分别进行迭代。博文作者用到的带约束最优化为：

ε(c)=∑q=1Q−1∑p=1P[∑n=0NcnMnp,q−Rq,q 1∑n=0NcnMnp,q 1]2

f(M)=∑n=0NcnMn objective:mincε(c) s.t.⎧⎩⎨⎪⎪⎪⎪f′(M)f′(0)f(0)≥0,0<M<1≥0≥0

需要注意的是以上约束条件中 f′(0) 是在 0 处的右导数。

博文作者重新编程，并取N=8重新进行仿真。得到如下的仿真结果：

可以对照Debevec的论文中的算法，对同一组数据进行计算得到下面的结果。

不难发现经博文作者修改后的算法和Debevec算法一样，都能满足 - Radiometric值均大于0 - 计算出来的曲线满足单调性

Mitsunaga Nayar对HDR的细分领域有着突出贡献，博文作者对Radiometric Self Calibration算法的修改也是基于Mitsunaga Nayar的基本框架。

---