基于各向异性全变分的GBR参数估计

Estimation of GBR Parameters Based on Anisotropic Total Variation

摘要: GBR参数估计是求解未标定光度立体技术的一般方法。本文针对GBR参数估计中模糊深度不连续的问题，引入具有保持深度不连续的各向异性扩散全变分做正则建立新的变分模型，并运用具有超线性收敛的拟牛顿法对模型进行求解。实验表明各向异性扩散全变分方法能有效降低法线误差，提高深度图不连续，且拟牛顿法求解该模型所用时间远低于EM方法和DM方法。

Abstract: GBR parameter estimation is a general method for solving uncalibrated photometric stereoscopic technique. In order to solve the problem of fuzzy depth discontinuity in GBR parameter estimation, a new variational model is established by introducing the anisotropic diffusion total variational model with depth discontinuity, and the quasi-Newton method with superlinear convergence is used to solve the model. Experiments show that the anisotropic diffusion total variational method can effectively reduce the normal error and improve the discontinuity of depth map, and the time of quasi-Newton method to solve the model is much lower than that of EM method and DM method.

1. 引言

三维重建是指由二维图像得到三维场景的非接触式的重建过程，有很高的实用价值，如应用于医学影像、文物数字化和虚拟现实等 [1] 。光度立体三维重建是利用同一视角具有不同光度信息的图像进行物体表面恢复的方法，光照信息未知的光度立体问题称为未标定光度立体问题。若图像通过正交投影获得，Hayakawa [2] 表明可以在这一问题中应用奇异值分解法(Singular value decomposition, SVD)，再增加表面可积性约束条件，则模糊性减小为通用浅浮雕(Generalized bas-relief, GBR)模糊性 [3] 。在未知光源条件下，使用正交投影会产生GBR问题，使用透视投影不会产生相同问题 [4] ，因此透视投影模型可以解决GBR模糊性问题，但是该模型比远光源的正交投影模型复杂得多，且计算困难。因此，GBR参数估计是未标定光度立体问题的一般解决方向。

Alldrin等人 [5] 通过最小化反照率分布的熵值求解GBR参数。Favaro和Papadhimitri [6] 利用局部漫反射最大值解决GBR歧义。Yvain Quéau [7] 通过法向量场的全变分(Total Variation, TV)最小化来估计GBR参数，提高了重建精度。

GBR参数估计是解决未标定光度立体模糊性问题的关键方向，其仍是一个有待研究的问题。TV最小化广泛应用于图像去噪领域 [8] ，已取得很大进展，其特点是能很好的保持边缘特性，所以可以应用到三维重建领域。本研究借鉴Yvain Quéau等人的思想，在GBR模糊性的基础上引入深度函数u的TV最小化，利用深度正则化估计GBR参数，从而提出了基于各向异性全变分(Anisotropic Total Variation, ATV)的三维重建模型，提高了重建精度。

2. GBR变换

在正交投影下的未标定光度立体问题中，光源信息未知，此时需要同时估计法向量和光源信息。这是一个不适定问题，这一不适定问题可以减小为GBR模糊性。

GBR模糊性是指当从某一特定角度观察物体时，不能对物体进行准确的深度估计。当使用正交投影观察一个具有朗伯反射率的物体时，其深度的确定会存在歧义。对于物体和通过GBR变换后的物体，选择不同的光照角度，这两种不同的物体在某些视角下可以得到完全相同的图片。

记I为图像灰度矩阵，矩阵中每一列为在不同光照信息下得到的图像灰度值， $n\left(x,y\right)={\left(n_x,n_y,n_z\right)}^{\text{T}}$ 是像素点 $\left(x,y\right)$ 的表面法向量， $\rho$ 是将像素点的反照率按列排列的向量生成的对角矩阵，N是像素点的归一化法向量排成的矩阵， $L^i$ 是第i个光源向量， $i=1,2,\cdots ,m$ ，其中 $m\ge 3$ 是光源数量，L是这些光源排列成的光照矩阵，令 $M=\rho N$ ，则I可表示为

$I=ML.$ (1)

将I应用SVD方法进行主成分分析，仅取前三个最大的奇异值及其对应的特征向量，这样可以保留矩阵的主特征，则图像矩阵I可表示为

$I=UD{V}^{\text{T}},$ (2)

其中， $U\in R^{R\times 3}$ ， $D\in R^{3\times 3}$ ， $V\in R^{m\times 3}$ ，R为像素点总数量。记P和Q是满足 $P{Q}^{\text{T}}=D$ 的两个3 × 3矩阵，则存在无穷多个这样的 $\left(P,Q\right)$ 矩阵，此时得到法向量值和实际值相差的一个3 × 3参数的可逆线性变换A，即：

$I=\left(U{P}^{\text{T}}\right)\left(Q{V}^{\text{T}}\right)=\hat{M}\hat{L}={\hat{M}}^1\left(A{A}^{-1}\right){\hat{L}}^1=\left({\hat{M}}^{1}A\right)\left(A^{-1}{\hat{L}}^1\right),$ (3)

其中 $\left({\hat{M}}^1,{\hat{L}}^1\right)$ 是 $\hat{M}$ 和 $\hat{L}$ 的一组解。假设表面连续，即在此基础上增加表面可积性约束，则有 $A=G$ ， $A^{-1}=G^{-1}$ ，即可以将此线性变换模糊性减少为GBR模糊性。GBR矩阵及其逆矩阵为：

$A=G=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ \mu & \vartheta & \lambda \end{array}\right),A^{-1}=G^{-1}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ -\frac{\mu }{\lambda }& -\frac{\vartheta }{\lambda }& \frac{1}{\lambda }\end{array}\right).$ (4)

3. 基于各向异性TV的GBR参数估计

3.1. 基于各向异性TV的GBR参数估计模型

基于L1范数的各向异性全变分定义为：

$ATV\left(u\right)={\displaystyle {\int }_{\Omega }\left(\left|{\nabla }_{x}u\right|+\left|{\nabla }_{y}u\right|\right)\text{d}X}.$ (5)

Yvain Quéau分别观察了 $\mu ,\vartheta$ 和 $\lambda$ 对 $TV\left(u\right)$ 的影响，当观察其中一个参数的影响时，固定其他两个参数为实际值，发现 $\mu$ 和 $\vartheta$ 越接近实际值时， $TV\left(u\right)$ 的值越小。ATV的特性是可以更好地保留物体的边缘结构信息，因此本文考虑通过深度场的ATV最小化估计 $\mu$ 和 $\vartheta$ 。

GBR变换影响深度及其梯度的变化，假设深度函数 $u(x,y)$ 分段两次可微，则几乎处处满足可积性约束 $u_{xy}=u_{yx}$ 。Yuille和Snow提出了在最小二乘的意义上M场满足可积性约束且M和L满足朗伯体定律的一个封闭解 [9] ，称为 $\left(M^0,L^0\right)$ ，取其作为初始解，用 $u^0\left(x,y\right)$ 表示这一方法估计的法向量场 $N^0$ 积分得到的深度场，其中， $M^0=\left(M_1^0,M_2^0,M_3^0\right)$ ， $N^0=\left(N_1^0,N_2^0,N_3^0\right)$ 。

在GBR变换的作用下， $M^0$ 经变换后得到的M场为：

$\bar{M}=\left(M_1^0+\mu M_3^0,M_2^0+\vartheta M_3^0,\lambda M_3^0\right).$ (6)

对应的深度 $\bar{u}$ 的梯度场为：

$\nabla \bar{u}=-\frac{1}{\lambda }{\left(\frac{M_1^0}{M_3^0}+\mu ,\frac{M_2^0}{M_3^0}+\vartheta \right)}^{\text{T}}.$ (7)

令 $U^0={\left(x,y,u^0\left(x,y\right)\right)}^{\text{T}}$ ， $g\left(\mu ,\vartheta ,\lambda \right)=\frac{1}{\lambda }\left(-\mu ,-\vartheta ,1\right)$ ， $g\left(\mu ,\vartheta ,\lambda \right)$ 对应 $G^{-1}$ 矩阵的第三行元素，则通过简单的微积分得到GBR变换下的深度函数：

$\bar{u}=g\left(\mu ,\vartheta ,\lambda \right)U^0.$ (8)

GBR参数的估计可以描述为(9)式：

$\begin{array}{l}\hat{u}=\arg \min ATV\left(\bar{u}\right),\\ {\Vert \bar{M}\bar{L}-I\Vert }^2={\Vert \rho L\frac{\left(-{\partial }_x\bar{u},-{\partial }_y\bar{u},1\right)}{\sqrt{\nabla {\bar{u}}^2+1}}\bar{L}-I\Vert }^2=\min .\end{array}$ (9)

由于 ${\Vert M^0{L}^0-I\Vert }^2={\Vert \left(M^{0}G\left(\hat{\mu },\hat{\vartheta },\hat{\lambda }\right)\right)\left(G^{-1}\left(\hat{\mu },\hat{\vartheta },\hat{\lambda }\right)L^0\right)-I\Vert }^2={\Vert \bar{M}\bar{L}-I\Vert }^2$ ，因此，已知 ${\Vert {\rho }^0\frac{\left(-{\partial }_x{u}^0,-{\partial }_y{u}^0,1\right)}{\sqrt{\nabla u^0{}^2+1}}{L}^0-I\Vert }^2=\min$ ，则有 ${\Vert \bar{\rho }\frac{\left(-{\partial }_x\bar{u},-{\partial }_y\bar{u},1\right)}{\sqrt{\nabla {\bar{u}}^2+1}}\bar{L}-I\Vert }^2=\min$ 。

由于 $\lambda$ 的符号与物体的凸性直接相关，而凹凸模糊性不能在没有物体附加信息的情况下求解，因此通常假设 $\lambda >0$ 。由此，(9)的解可由式(10)和(11)给出：

$\hat{u}=g\left(\hat{\mu },\hat{\vartheta },\hat{\lambda }\right)U^0,$ (10)

其中，

$\left(\hat{\mu },\hat{\vartheta },\hat{\lambda }\right)=\underset{\begin{array}{c}\mu ,\vartheta \\ \lambda >0\end{array}}{\arg }\min ATV\left(g\left(\mu ,\vartheta ,\lambda \right)U^0\right).$ (11)

为了准确估计GBR参数，关注的是式(11)。由(7)式可知， $\mu ,\vartheta$ 的估计独立于 $\lambda$ 的估计，同时可知，最小化 $ATV\left(u\right)$ 导致 $\lambda$ 趋于无穷大，因此，仅利用深度函数的ATV最小化估计 $\mu$ 和 $\vartheta$ 的值是可行的，此时，先取 $\lambda$ 为固定值1。通过下式估计 $\mu$ 和 $\vartheta$ ：

$\left(\hat{\mu },\hat{\vartheta }\right)=\underset{\mu ,\vartheta }{\arg }\min ATV\left(g\left(\mu ,\vartheta ,1\right)U^0\right).$ (12)

利用反照率熵最小化的方法估计GBR变换中的参数 $\lambda$ ，对 $\lambda$ 可能值进行由粗到精的估计，将式(12)和熵最小化模型 [5] 相结合，则有以下模型：

$\begin{array}{l}\left(\hat{\mu },\hat{\vartheta }\right)=\underset{\mu ,\vartheta }{\arg }\min ATV\left(g\left(\mu ,\vartheta ,1\right)U^0\right),\\ \hat{\lambda }=\underset{\lambda }{\arg }\min f\left(\Vert M^{0}G\left(\hat{\mu },\hat{\vartheta },\lambda \right)\Vert \right),\end{array}$ (13)

其中， $f\left(\Vert M^{0}G\left(\hat{\mu },\hat{\vartheta },\lambda \right)\Vert \right)$ 表示变换后反照率分布的熵值。

3.2. 拟牛顿法求解模型

深度函数 $u\left(x,y\right)$ 的ATV最小化可以减少未标定光度立体的模糊性。在这一过程中，深度函数的ATV最小化只依赖于GBR参数，由此产生的问题是凸的，可以使用拟牛顿的凸优化方法求解式(13)，目标函数为：

$E\left(\mu ,\vartheta \right)=ATV\left(\bar{u}\right)={\displaystyle {\int }_{\Omega }\left(\left|{\nabla }_x\bar{u}\right|+\left|{\nabla }_y\bar{u}\right|\right)\text{d}X}.$ (14)

$E\left(\mu ,\vartheta \right)$ 分别对 $\mu$ 和 $\vartheta$ 求偏导，得到偏微分方程：

$\begin{array}{l}{\partial }_{\mu }E={\displaystyle {\int }_{\Omega }\frac{M_1^0+\mu M_3^0}{\left|M_1^0+\mu M_3^0\right|}\text{d}X},\\ {\partial }_{\vartheta }E={\displaystyle {\int }_{\Omega }\frac{M_2^0+\vartheta M_3^0}{\left|M_2^0+\vartheta M_3^0\right|}\text{d}X}.\end{array}$ (15)

已知目标函数及其偏微分方程，可用拟牛顿法快速估计 $\mu$ 和 $\vartheta$ ，为加快估计速度，选择Tikhonov正则化 [10] 的解 ${\hat{\mu }}_0$ 和 ${\hat{\vartheta }}_0$ 作为一组初始解，其中， ${\hat{\mu }}_0$ 和 ${\hat{\vartheta }}_0$ 分别代表 $-\frac{M_1^0}{M_3^0}$ 和 $-\frac{M_2^0}{M_3^0}$ 的平均值。

以下为具体算法： $\left[\mu ,\vartheta \right]\leftarrow \left(M^0,L^0\right)$

1：给定初始值 $x_0=\left[{\mu }_0,{\vartheta }_0\right]\leftarrow \left[-mean\left(M_1^0./M_3^0\right),-mean\left(M_2^0./M_3^0\right)\right]$ 和容许误差 $\epsilon >0$ ， $k=0$ ， $H_0=I$ ,计算梯度 $g_k$ ；

2：若 $\Vert g_k\Vert <\epsilon$ ，则迭代停止；

3：计算迭代方向 $d_k=-H_k{g}_k$ ；

4：从点 $x_0$ 出发，沿着 $d_k$ 做一维搜索求步长 ${\alpha }_k$ ，令 $x_{k+1}=x_k+{\alpha }_k{d}_k$ ；

5：计算梯度 $g_{k+1}$ ，令 $s_k=x_{k+1}-x_k$ ， $y_k=g_{k+1}-g_k$ ；

6：若 $s_k^{\text{T}}{y}_k>0$ ，则令 $H_{k+1}=\left(I-\frac{\Delta x\Delta g^{\text{T}}}{\Delta g^{\text{T}}\Delta x}\right)H_k\left(I-\frac{\Delta g\Delta x^{\text{T}}}{\Delta g^{\text{T}}\Delta x}\right)+\frac{\Delta x\Delta x^{\text{T}}}{\Delta g^{\text{T}}\Delta x}$ ；否则，令 $H_{k+1}=H_K$ ；

7：令 $k=k+1$ ，转步骤2。

4. 实验结果与分析

漫反射极大值(Diffuse Maxima, DM)方法求解GBR参数是把最亮的像素看作朝向光源的像素，这是已知的有效方法；熵最小化(Entropy Minimization, EM)方法通过最小化反照率的熵来估计GBR参数，考虑的是全局属性，与本文方法有相同之处。因此本文所提算法的结果将与EM、DM方法进行对比。

为定量、定性评估构造方法的有效性和鲁棒性，选取4个图像数据集进行了实验，其中，Owl数据集由12张分辨率为512 × 340的图像组成；Frog、Pig和Scholar数据集均由20张灰度图像组成，分辨率分别为690 × 590、760 × 550和1070 × 790。这些物体几乎是朗伯体物体。下面给出每一个数据集在不同光照条件下的三幅图像，如图1所示。

Figure 1. Images of objects under different lightsource directions

图1. 不同光源方向下的物体图像

为定量分析本文方法的实验结果，将标定光度立体三维重建的结果作为真实值，将实验结果与真实值相比较，引入法线的平均角度误差(mean angular error，MAE)作为一个评价指标，定义如下：

$\text{MAE}=\frac{1}{R}{\cos }^{-1}\left(N_r\cdot {\hat{N}}_r\right),$ (16)

其中R是像素点总数量，r是图片上的像素点。使用法线角度误差小于5°的像素点百分比和运行时间作为另两个评价指标。

为了客观比较，计算了EM、DM和本文方法估计的法线的MAE、误差小于5°的像素点百分比和CPU运行时间，如表1所示。

Table 1. Reconstruction error and calculation time of each method

表1. 各方法的重建误差和计算时间

表1列出了三种方法估计的法线的误差数据和CPU运行时间。可以看到，Owl数据集和Pig数据集通过本文方法得到的法线的平均角度误差低于EM方法和DM方法的误差，本文方法的结果取得了更高的精度；Owl数据集中，DM方法得到的误差小于5°的像素点百分比为0.96%，本文方法得到的误差小于5°的像素点百分比为21.98%，是前者的近23倍；Scholar数据集中，EM方法的运行时间为147.30s，本文方法的运行时间为23.53s，是前者的近六分之一。因此，无论从法线的MAE还是误差小于5°的像素点百分比上看，本文方法与其他两种方法相比都有较好的估计结果，且计算时间更短。ATV具有比熵更好的微分性质，故本文构造的方法的运行时间大大短于EM方法。EM方法利用全局优化来求解三个GBR参数，DM方法通过确定图像局部亮度最大值区域求解GBR参数，误差均较小，但当图像分辨率高时，都会导致时间较长，而本文提出的方法可以缩短运算时间。

以下展示Owl数据集的真实反照率图和在三种方法下生成的反照率图，见图2，由图可知本文方法生成的反照率图更接近反照率实际值，与其他两种方法相比误差更小。

(a) 真实反照率图 (b) EM方法 (c) DM方法 (c) 本文方法

Figure 2. The real albedo map and the albedo map generated by the three methods

图2. 真实反照率图以及三种方法下生成的反照率图

将本文方法与EM、DM方法进行定性比较，结果如图3所示。图3展示了数据集中四个物体的真实法向图、三种方法得到的法向图以及生成法向图和真实法向图之间的误差图。法线角度误差低的地方表

(a) 真实法向图

(b) EM方法生成的法向图及其误差图

(c) DM方法生成的法向图及其误差图

(d) 本文方法生成的法向图及其误差图

Figure 3. The real normal graph, the normal graph generated by three methods and the error graph

图3. 真实法向图、三种方法生成的法向图及误差图

现为蓝色，误差高的地方表现为红色。通过定性比较可知，本文方法估计的法向量误差较小，能有效提高表面重建的精度。

在图4中，展示了Pig数据集的三维重建结果在三个不同视角下观察到的图片，其中深度是由本文方法估计的法线积分后得到的。

(a) EM方法 (b) DM方法(c) 本文方法

Figure 4. 3D reconstruction results of each method are in the front, right front and right top view of the picture

图4. 各方法的三维重建结果在正面、右前、正上方视角下的图片

从视觉效果上看，相较于EM、DM方法，本文方法的三维重建实验结果较好地恢复了物体实际的三维形状，深度的凸出效果明显，更好地恢复了物体的表面细节，重建精度高，表现出了良好的综合性能。

5. 结论

本文提出了一种基于各向异性全变分的三维重建方法，该方法在GBR模糊性的基础上，利用ATV最小化估计GBR参数 $\mu$ 和 $\vartheta$ ，通过实验证明，提出的方法取得了较高的重建精度，当图像分辨率更高时，在保证重建精度的同时有效缩短了计算时间，提高了计算效率。这表明ATV在解决未标定光度立体三维重建问题中是适用的，其可以有效降低GBR模糊性，提高计算效率。

基金项目

吉首大学科研项目(Jdy23029)。