《测绘学报》

朱艳¹, 颜青松¹,曲英杰¹,陈欣¹,邓非^1,2                            

1.武汉大学测绘学院, 湖北 武汉 430079;
2.国土资源部城市土地资源监测与仿真重点实验室, 深圳 518000

基金项目：国土资源部城市土地资源监测与仿真重点实验室开放基金（KF2018-03-025）

摘要：在三维重建中，网格优化通常用于解决密集点云构建的三角网格含有较多噪声且缺少细节的问题。现有的变分优化方法利用完整的影像数据对初始网格进行影像一致性优化，但在一定程度上忽视了影像信息的冗余以及视图的质量对网格优化的影响。对此，本文提出主视图选取与从视图选取策略，以提升网格优化的效率与质量。首先综合影像梯度幅值与轮廓检测，构建马尔科夫随机场，为每个三角面选取主视图；其次根据相应的观测条件为每个主视图选取从视图；然后计算主、从视图间归一化加权的影像一致性；最后利用梯度下降法最小化表面能量函数，实现网格优化。试验在定性和定量上证实了本文方法的有效性，表明本文方法能恢复更多精细细节，且优化的时间更短、精度更高。

关键词：视图选取    影像一致性    变分优化    三维重建    

近20年来，随着多视立体(multi-view stereo，MVS)基准评价集^[1-3]的提出，利用高分辨率多视图影像进行三维重建的研究取得了重大进展^[4]。相比于借助激光雷达等主动测距的方法，基于影像的被动式的多视立体方法具有重建成本低、速度快、非接触和无伤害等优势，适用于高精度大规模户外场景重建^[5-6]。车载移动测量技术^[7]与无人机航拍技术^[8]的发展使得更多科研人员将其应用于地面与航拍影像的三维重建，促进了智慧城市等的建设。对此，文献[9]提出了一套由带有相机标定参数的影像进行多视立体重建的完整流程，主要包括密集匹配、网格重建与网格优化3个步骤。其中，网格重建是指利用密集匹配得到的稠密点云进行构网操作。常用的方法有：Poisson表面重建^[10-11]等隐式拟合法，即通过Poisson方程建立指示函数与表面法向间的关联，进而提取零交叉等值面；Delaunay四面体图割的方法^[12]，即先对点云进行Delaunay四面体空间剖分，再基于可视性信息进行最小s-t图割，从而提取三角网格。Poisson表面重建等方法没有考虑点在影像上的可视性信息，且重建结果的平滑程度较高，不能很好地重建出复杂表面；Delaunay四面体图割的方法则加入了点的可视性，对噪声、异常值以及点云密度的变化更稳健，故文献[9]在其提出的多视立体重建流程中选用此方法进行网格重建。

网格重建得到的三角网格通常存在噪声且缺少细节，因此需要进行网格优化以提升网格质量。对此，国内外已有一些相关研究进展。早期，文献[13]用偏微分方程定义水平集，在变分框架下利用影像一致性进行网格优化。文献[14]将影像一致性的积分范围从三维网格转到二维影像，得到了更准确的优化结果。文献[15]进一步用梯度代替偏微分方程，提高了优化的稳健性。然而以上基于水平集的优化方法通常较为复杂、难以理解，故文献[16-17]利用基于网格的多视立体变分方法，通过相机与网格间的观测关系以及相应的影像数据，对网格的顶点位置进行调整。同样地，文献[9]通过定义表面能量函数对连续水平集进行离散化，将优化直接作用于网格顶点，在网格的准确性与完整性上均取得了较佳的效果。

诸多学者进一步对该网格优化方法进行了各方面的改进：文献[18]在表面能量函数中添加轮廓匹配项以强化边缘，并将影像一致性作为正则项的权重以实现自适应平滑；文献[19]利用光束法平差同时优化相机标定参数与表面网格；文献[20]使用各向异性的影像一致性测度与表面降噪方法，更好地保持尖锐特征；文献[21]对网格区域进行活动或非活动的二元标记，仅对标记为活动的网格区域进行优化，提高了优化效率；文献[22]结合影像一致性与语义一致性，同时进行网格几何优化和语义分割优化；文献[23]在优化前进行移动物体检测，避免相应影像区域影响场景中静态模型的优化；文献[24]利用深度图获取遮挡感知掩膜，得到有效像素，处理由深度的不连续性造成的遮挡问题。

尽管上述方法从能量函数构成、网格分割、遮挡检测等各方面对网格优化的方法进行了完善，但仍然在一定程度上忽视了影像信息的冗余以及视图的质量对网格优化的影响。对此，本文提出在计算影像一致性的过程中，为网格每个顶点的邻接三角面进行视图的选取，从而使更优质的视图更多地参与网格优化。虽然密集匹配已提出在像素级别上进行视图选取以得到准确点云^[25-26]，但是网格优化方面的文献总是局限于选取公共的、有合理视差的三维点数目最多的视图对。如文献[9]直接采用稀疏点云估计中的视图对；文献[6, 19, 23]则利用摄影基线、主视方向间的夹角、影像重叠度等对视图对在全局范围内进行选取；文献[27]利用图简化方法对大规模、不规则采样的影像集进行压缩，得到影像子集。而本文提出的为每个三角面选取视图的变分优化方法，是从网格中每个面的角度出发，进行局部视图选取。

文献[9]提出的网格优化方法在多个标准数据集上得到的结果均具有较优的准确性与完整性，故本文将其作为基础方法。该方法借助带有相机标定参数的影像数据对网格进行优化：利用相机与表面间的观测关系构建表面能量函数，表示表面在多张影像上的投影的匹配程度，并将网格重建得到的三角网格作为函数自变量的初值，通过应用梯度下降法使能量最小化，从而实现网格优化，达到减少噪声、增强细节的目的。

表面能量函数包括影像一致性项与正则项

 (1)

式中，为模型表面；E_photo为影像一致性项；E_fair为正则项；λ为比例参数。此处正则项采用伞算子近似拉普拉斯算子，进而模拟离散化的薄板能量函数，对表面进行平滑约束^[28]。由于本文主要研究利用视图对的影像一致性对网格进行优化，故接下来重点阐述影像一致性的计算及如何最小化能量以优化网格。

影像一致性用于度量同一位置的表面在可视的不同影像上的投影间的相似性，为此，定义x∈S为表面上一点，N为表面在该点处朝外的法向量，g_i, j(x, N)描述(x, N)在视图i与视图j间的影像一致性(其中视图包括影像及对应的相机标定参数信息)，v_i, j^S(x)∈0, 1描述x在i和j上的可视性^[13]。则描述影像一致性的原始能量函数为

 (2)

式(2)在三维空间上积分困难，为了更方便、有效地计算影像一致性，引入重投影误差，将积分空间转换到二维影像上^[29]。如图 1所示，x_i为视图i对应的影像上的一个像素点，反投影到表面S上得到点x，再由点x投影到视图j对应的影像上得到像素点x_j，x_j与x_i的匹配程度即反映重投影误差。本文称i为主视图，j为主视图的从视图。由此，式(2)可改写为

 (3)

式中，I_i为视图i对应的影像；为由表面诱导的与I_i相对应的I_j(Π_i表示由表面投影到I_i，Π_i^-1表示由I_i反投影到表面)；Ω_i, j^S为重投影误差在I_i上的定义域；h(I_i, I_i, j^S)(x_i)描述(I_i, I_i, j^S)关于x_i的影像一致性。另定义M为重投影误差测度，则视图对(i, j)的重投影误差可表示为M(i, j)=∫h(i, j)(x_i)dx_i。考虑对噪声和光照的变化的稳健性，重投影误差测度M采用影像间的负的零均值归一化互相关(zero-mean normalized cross correlation，ZNCC)^[29]。

为了对能量函数进行离散化，将表面S表示为三角网格，有n个顶点, i∈[1, n]。对S上任意一点x，有x=∑_iϕ_i(x)v_i，且∑_iϕ_i(x)=1(当v_i为包含有x的三角面的顶点时，ϕ_i(x)为x关于v_i的重心坐标，否则ϕ_i(x)=0)。省略v_i与ϕ_i(x)的下标、Ω_i, j^S的上标，对式(3)关于网格顶点求导^[30]，可得

 (4)

式中，f_i, j(x_i)=∂₂M(x_i)DI_j(x_j)DΠ_j(x)d_i；d_i为由i对应的相机中心指向x的向量；∂₂M为M关于第2张影像(即从视图对应的影像)的导数，D表示对应函数的雅可比矩阵。

综上可知，基础方法对每一个视图对(i, j)，由主视图i对应的影像I_i逐像素x_i反投影到表面点x，借助重投影误差计算(i, j)关于x_i的影像一致性的梯度，并将该梯度按重心坐标分配到x所在三角面的3个顶点上。换言之，获取的每个顶点的最终梯度，即为可视的所有视图对关于该顶点对应的投影点的影像一致性的梯度和。再应用梯度下降法，按梯度反向迭代移动网格顶点，使顶点移动到更准确的位置，同时能量迭代下降至最小化，最终得到影像一致性最优的优化网格。

根据上述基础方法，在影像一致性的计算过程中，将需要经由网格进行重投影操作的视图定义为主视图，并将主视图重投影的对象视图定义为从视图，主视图与对应的从视图构成视图对。由式(4)可知，在基础方法的优化过程中，顶点的所有可视视图对均同等地参与其影像一致性的计算。考虑到视图的质量对影像一致性计算的影响，本文方法提出视图选取的策略，即从三角面的角度出发，新增主视图选取与从视图选取两个步骤。

如图 2所示，本文方法的主要输入为初始网格与带有相机标定参数的校准影像，首先在所有视图中为网格的每个三角面f选取最佳视图作为主视图i，即主视图选取(图中不同颜色代表选取的不同主视图)；然后根据观测条件为主视图i的所有从视图j设置相应权重，以此排序进而选取从视图，同时为被选取的i与j构成的视图对(i, j)关于三角面f的3个顶点的影像一致性计算进行归一化加权，即从视图选取；最后，按照基础方法求取能量函数在网格顶点上的梯度，并利用梯度下降法进行迭代求解，由初始网格得到优化网格。

据此，在本文方法中，基础方法中的式(4)可改写为

 (5)

式中，p_i(x)∈{0, 1}表示主视图i是否为点x所在的三角面f的最佳视图，若是，则值为1，否则值为0；ω_i, j(x)表示视图对(i, j)关于点x所在的三角面f的3个顶点的影像一致性的归一化权重。接下来分别详细介绍主视图选取与从视图选取。

2.1  主视图选取

构建关于三角面的最佳视图的马尔科夫随机场，为网格的每个三角面f_i分配标签l_i，以此为每个三角面选取最佳视图，并将其作为用于计算影像一致性的视图对的主视图。为此，定义能量函数

 (6)

式中，l为所有面的最佳视图标签集；f_i和f_j为有公共三角边的邻接面；μ为比例参数，用于平衡数据项E_data与拓扑平滑约束项E_smooth。数据项用于度量利用影像一致性进行优化时，三角面在其可视视图上的质量；拓扑平滑约束项用于使邻接面的最佳视图标签尽可能保持一致。

创建无向图，将三角面作为图的节点，邻接三角面的公共边作为图的边，利用基于图割的α扩张算法^[31]最小化能量函数(6)，得到最优标签集，从而得到每个三角面的最佳视图，实现主视图选取。

基于波特模型，拓扑平滑约束项定义为

 (7)

式中，[·]为艾弗森括号，若l_i≠l_j，则值为1，否则为0。

综合影像梯度幅值与轮廓检测，将轮廓因子C作为惩罚项与梯度幅值M相乘，由此可定义三角面f在其可视视图i上的质量

 (8)

由面f在其所有可视视图中的最大质量max(Q_f)与最小质量min(Q_f)得到归一化质量

 (9)

则数据项定义为

 (10)

接下来，分别对数据项中的影像梯度幅值与轮廓检测进行具体说明。

2.1.1  影像梯度幅值

为了反映三角面在其可视视图上的质量，考虑影像分辨率以及相机与被观测表面间的邻近性，尽量为三角面选取影像分辨率较高的、相机距被观测表面较近的视图，为此计算三角面在影像上的投影所占的像素的梯度幅值之和^[32]。利用索贝尔算子近似获取影像梯度：索贝尔算子包含横向S_x与纵向S_y，如图 3所示，将之与视图i对应的影像进行平面卷积，分别得到影像在横向与纵向上的梯度，即▽(i)_x=S_x*i与▽(i)_y=S_y*i。故像素点p的梯度幅值为

 (11)

则三角面f在其可视视图i上的影像梯度幅值之和表示为

 (12)

式中，P_f, i为f在i的影像上的投影所占的像素集合。若三角面的投影区域不足一个像素，则按其投影重心对梯度幅值进行采样，再乘以其投影面积。

2.1.2  轮廓检测

注意到当表面转折处在影像上的投影表现为轮廓时，该投影区域的梯度将对应为影像梯度的局部最大值^[33]，对于此种情形，仅依赖上述影像梯度幅值来为三角面选取主视图是不完全可靠的。如图 4所示，三角面f位于表面转折处，视图i₁和视图i₂均为f的可视视图。由于f在i₂影像上的投影表现为位于轮廓附近，故该处影像梯度较大，此时若仅根据影像梯度幅值判断f在i₁和i₂上的质量，可能会得到i₂优于i₁的错误结论。因此，本文提出，为三角面在其可视视图上作轮廓检测。

为进行轮廓检测，定义v₁、v₂和v₃为S的3个顶点，v_f=(v₁+v₂+v₃)/3为f的中心，N_f为f朝外的法向量，N_v1、N_v2和N_v3分别为表面在v₁、v₂和v₃处朝外的法向量，c_i为视图i对应的相机中心，v_view=v_f-c_i为由三角面的中心到相机中心的向量。则N_f与v_view间的夹角α为

 (13)

式中，〈, 〉表示向量内积。N_v1、N_v2和N_v3两两间的平均夹角β为

 (14)

当α大于某一阈值α_c且β大于某一阈值β_c时，认为三角面f在视图i的影像上的投影表现为位于轮廓附近。由此定义轮廓因子C_f, i为

 (15)

式中，η为比例参数。至此，综合轮廓因子与梯度幅值，即可得到式(8)。

2.2  从视图选取

主视图选取后，为进一步提升优化效率及质量，在主视图的所有从视图中，考虑从视图联合主视图对三角面的观测条件，据此为两者构成的视图对关于该三角面的影像一致性设置权重，即式(5)中的ω_i, j(x)，并依据该权重进行从视图的选取。

为选取从视图，本文主要考虑主、从视图间的视差以及相机到三角面的距离差。如图 5所示，三角面f的主视图为i，与该主视图构成视图对的从视图为j，用主、从视图视线间的夹角表示主、从视图间的视差，即θ_i, j=∠(c_iv_f, c_jv_f)。相机到三角面用相应焦距f_i、f_j归一化后的距离分别为l_i=||c_iv_f||/f_i、l_j=||c_jv_f||/f_j，用||c_iv_f||归一化后的距离差为Δl_i, j=|l_i-l_j|/||c_iv_f||。则分别由视差、距离差定义权重w_p、w_l

 (16)

 (17)

式中，为参考视差；σ_p、σ_l为标准差。由此，定义从视图j相对于主视图i的权重w_i, j为

 (18)

式中，γ为比例参数。因此，当三角面的主、从视图间的视差越接近参考视差、相机到三角面的距离越一致时，该视图对该三角面的影像一致性的贡献越大，使得优化过程更为合理、有效。

假设三角面f的主视图有n个从视图，对该n个从视图按其权重进行降序排列，若n大于参考最大从视图数N，则选取前N个从视图。对选取的从视图的权重归一化后得到式(5)中的ω_i, j(x)

 (19)

式中，k表示选取的从视图数，当n>N时，k=N，否则k=n。至此，将选取的主、从视图构成视图对，实现了对影像一致性计算的归一化加权。

OpenMVS开源库实现了多视立体重建的完整流程^[9]，故本文利用OpenMVS，在CPU i7-6700，3.40 GHz；RAM 31.96 GB的环境下进行试验。试验数据采用Strecha提供的fountain-P11、Herz-Jesu-P25与castle-P30 3套数据集^[2]，分别含有11、25与30张影像，影像分辨率均为3072×2048。与文献[9]相同，本文初始网格的构建方法采用引言部分提到的Delaunay四面体图割的方法^[12]。对于不同数据集，本文涉及的各参数统一设置如下：比例参数λ=0.2，μ=0.1，η=0.2，γ=1；角度阈值α_c=75°，β_c=15°；参考视差=50°；标准差σ_p=45°，σ_l=0.25；参考最大从视图数N=8。

3.1  定性分析

3.1.1  方法有效性分析

为论证本文方法的各组成部分在提升网格优化质量方面的有效性，对fountain-P11的优化细节逐步进行对比分析：图 6和图 7(a)-(f)依次为fountain-P11的真实网格、初始网格、基础方法优化网格、只考虑影像梯度幅值选取主视图的优化网格、综合影像梯度幅值与轮廓检测选取主视图的优化网格、主视图选取加从视图选取的优化网格(即本文方法)。图 6和图 7(c)相比于图 7(b)，利用影像一致性对初始网格进行优化，恢复了更多细节；(d)相比于(c)，凹陷与凸出等纹理细节得到了增强；而(e)由于考虑了轮廓检测，更好地保持了边缘轮廓；(f)则通过对视图对的影像一致性计算进行归一化加权，进一步丰富了细节同时减少了噪声。值得注意的是图 6(d)由于仅考虑影像梯度幅值为部分三角面选取了不恰当的主视图，(e)则通过加入轮廓检测，有效地解决了该问题。

为进一步说明主视图选取在网格优化中所发挥的作用，选取Herz-Jesu-P25 3个不同部位的优化细节进行对比分析：图 8(a)为25张输入影像中的两张，(b)为主视图选取结果，(c)和(d)分别为基础方法与本文方法的优化网格。(b)中不同颜色代表为三角面选取的不同主视图，黑色代表对应的三角面无主视图，即位于地面等处的三角面在其所有可视视图上的质量均未达到最佳视图选取条件，从而忽略对这些三角面的影像一致性优化，在一定程度上实现了自适应分辨率优化，提高了优化效率。分析(c)和(d)中3个部位的优化细节图：位于模型正面及中部的三角面在多个视图上均呈现为较优观测，故基础方法与本文方法得到的优化结果在细节上的差异不是特别突出，如门廊上的人形雕塑；而靠近模型边界或位于模型侧面的三角面仅在少数几个甚至一个视图上呈现为较优观测，在其他可视视图上则质量较差，此时本文方法相比于基础方法在细节恢复上效果有显著提升，如门廊立柱及右上角的精细纹饰。说明通过综合影像梯度幅值与轮廓检测为每个三角面选取主视图，减少了冗余视图的影像一致性的计算，不仅提高了优化效率，还在一定程度上避免了质量较差的视图对优化造成干扰。

3.1.2  与其他方法对比

利用castle-P30将基础方法、本文方法与其他主流方法进行结果对比：表 1为不同方法得到的网格顶点数与三角面数；图 9(a)为30张输入影像中的一张，图 9(b)和(c)分别为基础方法与本文方法的优化网格，图 9(d)-(g)依次为Poisson、MVE、Photoscan和ContextCapture(原Smart3D)方法得到的结果(各参数均按超高精度进行设置)。对比图中的窗户及墙壁上的纹饰等处，本文方法相比于其他方法细节更清晰、丰富。可见本文方法能更有效地利用影像信息，恢复更多精细细节。

表 1 castle-P30不同方法的顶点数与面数Tab. 1 The number of vertices and faces by different methods on castle-P30

表选项

3.2  定量分析

3.2.1  优化时间对比

表 2为fountain-P11、Herz-Jesu-P25与castle-P30在相同环境下分别使用基础方法与本文方法进行优化的执行时间，本文方法相比基础方法优化时间缩短了约40%，且输入影像数目越多，时间缩短比例越高。可见通过对主视图与从视图的选取，本文方法减少了不必要甚至不可靠的视图对的数量，从而简化了影像一致性的计算，加快了能量函数梯度下降的收敛过程，故本文方法可提高优化效率。

表 2 执行时间对比Tab. 2 Comparison of running time

表选项

3.2.2  优化精度对比

为定量评价网格优化的精度，利用fountain-P11的激光雷达扫描数据作为真值(ground truth)，比较从每个相机中心出发、经影像每个像素的视线分别到待测优化网格与真实网格的距离(即像素的待测深度值d与真实深度值d_GT)，具体表现为将深度值差值Δ_d=|d-d_GT|与标准差σ(本文设为0.001 5)进行比较，且令maxσ=30σ。据此绘制误差分布直方图，如图 10所示，表示深度值差值在一定误差范围内的像素数占影像总像素数的比例，第i组即满足(i-1)·3σ≤Δ_d < i·3σ，i=11时则表示Δ_d≥maxσ。在误差图上表现为：当Δ_d < maxσ时，对应像素的灰度值为255-255·Δ_d/maxσ，即优化结果越准确，越接近白色；当Δ_d≥maxσ时，对应像素为红色；当d_GT=0时，表示真值缺失，对应像素为绿色。图 11(a)-(d)第一行依次为初始网格、基础方法优化网格、本文方法优化网格及Poisson网格，第2行依次为第1行的网格所对应的误差图。

为评价优化精度，重点分析图 10中i=1时，即深度值差值在3σ内的像素占比：使用基础方法优化后的网格相比初始网格，比例提升了约6%；而本文方法在基础方法之上，比例再度提升了约5%，且本文方法相比Poisson(八叉树深度设为10，其他参数默认)，比例提升了约12%。对应于图 11中(c)与(a)、(b)和(d)相比，在墙面、雕饰等处更接近白色。可见通过对主视图与从视图的选取，本文方法降低了质量较差的视图对的影像一致性对网格优化的干扰，故本文方法可提高优化精度。另外，图 10中i=11时，即深度值差值大于或等于30σ的像素占比可用于表示优化完整性，该比例越小，则优化完整性越好，可见本文方法的完整性略低于基础方法。结合图 11进行分析，在地面部分，(c)的红色区域略多于(b)，因为本文方法在较大程度上依赖视图的质量，对于地面部分，观测角度倾向为斜视，缺少优质的下视视图对这些三角面进行影像一致性优化，故本文方法在优化完整性上还有待提高。

本文在多视立体变分优化的基础方法之上，提出了影像一致性优化的视图选取策略。首先综合影像梯度幅值与轮廓检测，为网格中的每个三角面选取主视图。其次，根据从视图的观测条件为其设置权重，进而选取从视图。再对选取的主、从视图所构成的视图对的影像一致性计算进行归一化加权，得到能量函数在网格顶点上的梯度。最后，利用梯度下降法进行迭代求解，由初始网格得到优化网格。试验结果及其分析表明本文方法通过视图选取提高了优化效率、强化了网格细节。同时也存在一定的局限性，如优化结果较依赖视图质量、优化的完整性有待提高。下一步将继续研究如何在视图选取与优化完整性上更好地取得平衡，以及尝试结合几何基元提取^[34]与影像一致性优化，进一步提升优化的视觉效果。

第一作者简介：朱艳(1996-), 女, 硕士生, 研究方向为摄影测量与计算机视觉。E-mail:yzhu@whu.edu.cn

通信作者：邓非, E-mail: fdeng@sgg.whu.edu.cn