CN111815757A

CN111815757A - 基于图像序列的大型构件三维重建方法

Info

Publication number: CN111815757A
Application number: CN202010496824.8A
Authority: CN
Inventors: 曹衍龙; 董献瑞; 王敬; 梁立威
Original assignee: Shandong Industrial Technology Research Institute of ZJU
Current assignee: Shandong Industrial Technology Research Institute of ZJU
Priority date: 2019-06-29
Filing date: 2020-06-03
Publication date: 2020-10-23
Anticipated expiration: 2040-06-03
Also published as: CN111815757B

Abstract

基于图像序列的大型构件三维重建方法，包括以下步骤：S1、使搭载有相机的无人机绕目标构件飞行，获得待重构的图像序列；S2、采用SIFT算法和SURF算法联合提取图像特征点；S3、基于SIFT角点与SURF角点得到的稀疏特征点，通过计算本质矩阵和基础矩阵估计相机运动，注册三维空间点获得三维场景的稀疏点云；S4、判断优化后的稀疏点云是否存在对称重复结构；S5、以稀疏点云作为种子点与参考图像输入，基于多视图的稠密三维点构建方法进行稠密重建，获得低分辨率的深度图。本发明的优点在于：给出了一种基于图像序列的三维点恢复及矫正方法，实现了从图像序列到空间稀疏三维点的构建。

Description

基于图像序列的大型构件三维重建方法

技术领域

本发明涉及基于视觉图像对大型构件进行三维重建的方法。

背景技术

基于图像的三维重建技术整体过程为：先从不同视角拍摄物体图像获取图像输入，而后提取图像集的特征点，通过图像间的特征点匹配关系建立相互关系，再利用多视图几何原理计算获取三维点坐标并估计相机位姿，再通过光束平差法进行三维信息优化得到相应的稀疏三维结果输出；其中特征点提取匹配过程中会受到光照、遮挡或场景中某结构重复出现等因素的影响，得到错误的结果，该结果将严重影响相机投影矩阵以及空间点坐标的解算结果，故而如何正确地提取并匹配特征点以及正确地获取三维几何关系是计算机视觉领域的一大难点。

通过上述流程，即运动恢复结构的方法，虽能有效地恢复相机位姿，注册空间三维点，但由于该方法所恢复的物体三维点云十分稀疏，难以直观辨认，亦难以直接用来重构场景的三维模型，故而需要结合场景稠密重建技术和表面重建技术等，来获取场景表面三维几何网格和纹理等方面的完整表达，达到进一步的使用。

在文物保护与数字展示中，可通过对文物进行图像或图像序列的拍摄来获取相应的颜色、纹理及几何信息等来对文物进行三维模型重构，后续便可用于破损文物的修复工作，或用于数字模型的3D展示及教育讲解等；在军事中也可利用图像的三维重建技术来对城市场景、地形地貌场景等进行重构生成，来用于军事作战人员的军事演习，以提高作战能力等；在其他方面，亦可于制造上应用于物体的缺陷检测定位等。

在基于图像的三维重建方法研究中，依据摄像机数目的不同可分为：单目视觉法、双目视觉法和多目视觉法。

单目视觉主要通过固定单一视点或多视点拍摄获取图像集，利用明暗度法、光度立体视觉方法、纹理法、轮廓法、调焦法和运动法等来获得三维重建结果；其中，运动法即运动恢复结构(Structure From Motion，SFM)，可根据图像间的匹配关系和几何约束来注册空间三维场景点；而该方法结果的好坏极大程度地依赖于特征提取与匹配的准确性；且其在输入图像集合数量较大时，由于运算量大，致使计算效率较低；同时由于单目在利用运动恢复结构方法进行场景点重构的过程中会对尺度归一处理，使其重建所得稀疏场景会缺失真实尺度信息。

双目视觉装置是将两个相机固定在两个不同视点，并对所获图像进行匹配处理取得图像像素偏移量，从而将视差结果转变为空间深度信息；该方法能够得到较好的重建效果，但运算量大，重建效果取决于基线距离。

同样，多目视觉是通过多个摄像机间的图像几何约束关系来获取重建结果的；对于单个建筑物重建，想要获取整体重建信息，需要有较大的基线，但大基线会减弱重建效果，故而双目或多目视觉方法不适用于进行大规模三维场景重建的应用。

对于基于图像的三维重建研究而言，其历史可追溯到上世纪六十年代；麻省理工学院的Roberts在他的博士论文中完成了对三维场景重建的分析，并通过计算机编程实现了从输入的二维图像中获取到真实物体的三维立体结构，从而正式拉开了基于图像的三维重建相关领域的研究帷幕。

之后各学者便在其研究基础上提出了其他描述数字图像重建物体三维模型过程的基础理论，其中以Marr的理论体系影响最为深远，更成为视觉研究领域的主流思想；但该视觉理论在通过多幅图像获取物体精确的三维几何信息上存在一定的困难，且无法定量地求取景物中物体的性质；后期，随着学者们将主动视觉的研究方法引入，使得很多Marr视觉理论框架下的病态问题转为良态。

随着前期三维重建基础理论的研究，许多方面都开始趋于成熟；而在此基础上，卡内基梅隆大学的Tomasi和Kanade等人于1992 年开发完成了第一个基于图像的三维重构系统；其主要假定摄像机为正交投影模型，而后利用仿射分解来解析物体三维信息并估计相机内外参数；但由于其假设相机为正交投影，仅适用于物体在较远处的情况，又因为其利用的特征点数量较少，从而使得整体模型的生成质量较低；后续K.U.Leuven大学的Pollefeys等人便基于图像重建技术实现了物体三维表面自动生成系统，并在欧洲的考古学、文物保护等领域得到了成功的应用；其主要结合相机自定标技术，将围绕待建物体拍摄所得的图像集作为输入，通过图像间对应点的稠密匹配，来重建表面模型以及相应的相机参数。

美国华盛顿大学Snavely等人开发的Photo Tourism三维重建系统及其随后放出的开源系统Bundler对后续的研究起到了很大的促进作用；Photo Tourism System可直接利用无序图像集来计算获取图像的相机参数并注册稀疏三维结构点，其图形界面更能使用户交互地移动三维空间，直观感知图像与三维场景的变换；而其开源的Bundler程序运行效率虽然相对于传统的软件有相应的提高，但由于该方法采用增量式方法处理图像，使得计算复杂度较高，在图像数量增加的情况下，其运行效率会显著降低。

2011年Wu开发出来的基于图像的三维重建系统Visual SFM，也采用了增量式方法，以普通图像输入，能得到度量重建的结果； Wu进一步在2013年为提高重建精度，提出了一个新的光束法平差策略，实现了重建速度和重建精度之间的良好平衡；通过不断地重新三角化原先未能完成三角化的特征匹配对，来保持重建的高准确度,并通过实验验证了其算法性能。

上述方法虽能取得较好的重建结果，但在整体重建技术中仍存在很多问题；例如增量式方法中图像对特征点匹配阶段匹配时间效率较低，需对特征点描述子采取聚类操作等来加速匹配效果；场景中存在的重复性结构会引发误匹配导致重建结果错乱等，其需从几何结构、时间序列或背景内容等方面上着手来获取正确的重建结果；以及整体的稠密重建过程效率较低等。

发明内容

本发明的目的在于提供一种通过无人机搭载相机对大型构件拍摄获得图像序列，基于图像序列对大型构件进行高效、完整的三维重建的方法；

基于图像序列的大型构件三维重建方法，包括以下步骤：

S1、使搭载有相机的无人机绕目标构件飞行，获得待重构的图像序列；

S2、采用SIFT算法和SURF算法联合提取图像特征点；结合关联度预匹配与层级哈希检索方式来进行特征匹配，关联度预匹配剔除匹配度低于预设特征点阈值的两幅图像的特征匹配；

S3、基于SIFT角点与SURF角点得到的稀疏特征点，通过计算本质矩阵和基础矩阵估计相机运动，注册三维空间点获得三维场景的稀疏点云；对稀疏点云进行捆绑调整优化；

S4、判断优化后的稀疏点云是否存在对称重复结构，若存在，则利用图像间的背景差异点对重建结果进行矫正，获得矫正后的稀疏点云；若无对称重复结构，则获得稀疏点云；

S5、以稀疏点云作为种子点与参考图像输入，基于多视图的稠密三维点构建方法进行稠密重建，获得低分辨率的深度图；以低分辨率的深度图作为输入，再基于多视图的稠密三维点构建算法进行稠密重建，获得高分辨率的稠密点云；基于稠密点云进行表面模型的构建，获得三维重建结果。

优选的，步骤5中将低分辨率的深度图结果作为高分辨率图像的初值进行深度优化，来获取最终的高分辨率的稠密点云，其主要步骤如下：

5.1、选取合适的图像层级L对稀疏种子点点集P进行扩散生长，获取点云p′，并将其投影到图像中构造像素集M。

5.2、取高分辨率图像层级L＝L-1，并对构造的像素集M进行置信度等插值拟合和邻近扩散，从而构造像素集M′。

5.3、对图像层级L进行视图选取，对不同的patch获取其置信度值T。

5.4、当所得置信度值大于阈值τ，则根据patch深度信息进行深度值的偏移调整；当小于阈值τ，则需将其加入优化队列Q进行进一步优化。

5.5、进行光度一致性约束，对整体优化队列进行优化获取L层的参考视图深度图结果。

优选的，基于多视图的稠密三维点构建方法为：将稀疏三维点作为种子点连同每个视角图像作为输入，而后通过全局与局部视图选择、空间区域生长以及多视图匹配优化处理来获得各视角深度图的最终输出。

优选的，基于双目立体匹配，对三维重建结果进行尺度恢复，执行的操作为：利用双目相机拍摄目标物体的局部，而后对物体进行局部重建，以获得局部的真实尺寸信息；再相对应地量取整体重建所得模型中对应的局部尺寸长度，从而获得尺寸比例来对整体重建模型进行缩放，恢复整体模型的真实尺度，后续便可直接对整体模型进行尺寸量取来获得真实的尺寸长度。

本发明的优点在于：

1、给出了一种基于图像序列的三维点恢复及矫正方法，实现了从图像序列到空间稀疏三维点的构建。利用联合的特征点提取算法，增加特征点数量，为稀疏点构建提供更多信息；改进了特征点匹配方式加快特征点匹配速度；并利用图像的显著背景特征解决了传统增量式重建方法在对称重复性结构物体重建中易出现错误情况的问题，从而得到正确的重建结果。

2、提出了一种基于多尺度图像的稠密三维点云构建方法和等值面提取的改进方法。该方法先利用了低分辨率图像进行深度图恢复，而后将其结果作为高分辨图像的初始值进行重建，其在获得良好重建结果的同时能够提高重建效率，且能提高后续的尺寸量取的准确性。

3、由于通过图像序列重建的三维模型缺失尺度信息，故而本文利用了双目相机低成本且便捷的特点，通过立体匹配算法获取局部真实尺度的点云，将其用于恢复重建模型的真实尺寸。实验表明该方法流程具有良好的精度结果。

附图说明

图1是基于图像序列的三维重构流程图

图2是SIFT特征点提取流程图

图3是高斯金字塔与高斯差分金字塔图像图

图4是高斯金字塔图像显示图

图5是极值点检测示意图

图6是y方向上的模板图示图

图7是混合模板图示

图8是SIFT算法与SURF算法尺度空间构造差异图

图9是SURF特征点主方向计算示意图

图10是SURF特征点描述子构造示意图

图11是SIFT角点提取图示图

图12是SRUF角点提取图示图

图13是角点提取结果图

图14是特征点匹配结果图

图15是描述子归一化与二值化图

图16是SIFT哈希映射示意图

图17是特征匹配耗时图

图18是特征匹配tracks数量结果图

图19是特征匹配匹配对数情况图

图20是匹配结果图

图21是图像间的连接图

图22是两图间的匹配示意图

图23是三图间的匹配关系图

图24是多视图的特征点跟踪的示意图

图25是对极几何示意图

图26是三角测量示意图

图27是三角测量不确定性

图28是过渡相机坐标系建立图

图29是过渡世界坐标系建立图

图30是半平面表示图

图31是半平面旋转示意图

图32是稀疏点重建结果图

图33是对称或重复性结构图

图34是错误重建结果图

图35是连接图分割为蓝绿两部分图

图36是背景点产生搭接图

图37是不同组连接图示意图

图38是未矫正重建结果图

图39是矫正后重建结果图

图40是稠密重建算法流程图

图41是尺度度量示意图

图42是patch模型图

图43是空间区域扩散生长示意图

图44是稠密重建结果示意图

图45是不同尺度图像重建结果图

图46是结果显示图

图47是八叉树示意图

图48是八叉树空间划分流程图

图49是符号距离场中的零平面示意图

图50是基础函数示意图

图51是权重函数示意图

图52是网格划分二义性图

图53是八叉树邻近面等值线划分不一致图

图54是体素示意图

图55是消除细化结果图

图56是退化网格去除图

图57是网格化结果图

图58放置已知标志物图

图59是高性能IMU或GPS模块图

图60是数据项的计算过程示意图

图61是16个搜索方向示意图

图62是双目立体匹配RGB图、视差图与三维立体图

图63是尺度缩放示意图

图64是双目点云与重建模型长度量取图

图65是尺寸比例缩放后模型的尺寸长度对比图

具体实施方式

结合附图，进一步说明本发明。

基于图像序列的重建方法并不依赖于额外设备来获取位置、方向或者几何结构等信息，而是运用计算机视觉与几何技术等通过图像本身得到。本发明主要利用了增量式运动恢复结构方法来估计相机位姿并注册空间稀疏三维点。由于重建流程强依赖于图像间特征点匹配的准确性，当图像间存在相似重复结构时，会引入错误匹配对，从而使重建结果错误。对于这种存在重复性结构的物体，本发明利用了图像间的独立观测点来进行重复性结构的检测和重建结果的矫正。

基于图像序列的大型构件三维重建方法，包括以下步骤：

S5、以稀疏点云作为种子点与参考图像输入，基于多视图的稠密三维点构建方法进行稠密重建，获得低分辨率的深度图；以低分辨率的深度图作为输入，再基于多视图的稠密三维点构建算法进行稠密重建，获得高分辨率的稠密点云；基于稠密点云进行表面模型的构建，获得三维重建结果；

图像特征点提取

特征角点检测在计算机领域中具有广泛的应用，其提取匹配亦是三维重建流程中关键的一步，其正确性将直接影响重建精度。在离线三维场景重建中，应用较多且效果较好的为SIFT算法和SURF算法两种。此处联合这两种算法进行特征点提取匹配。

常用的特征点提取算法

SIFT算法

SIFT算法是一种对亮度变化、旋转变换、平移变化、仿射变换等都具有良好稳定性的特征点提取算法。其主要计算流程如图2所示。

1.尺度空间极值检测

尺度空间极值检测是在图像中寻找对尺度变化保持不变的潜在极值点。其先利用式(3-1)通过高斯平滑图像和降采样两种方式，搭建高斯金字塔，如图3、4所示；

L(x，y，σ)＝G(x，y，σ)*I(x，y) (3-1)

式中：*为卷积运算符号，

为高斯卷积核，I(x，y)为灰度图像，σ为尺度空间因子， L(x，y，σ)为尺度空间函数。

假设高斯金字塔共O组(降采样获得)，S层(高斯平滑获得)，则可得如式(3-2)：

式中：o为组别的索引值,s为组内层数的索引。特征点的尺度坐标σ根据组别与层数的索引共同获得。为得稳定的特征点，再利用式(3-3)建立高斯差分尺度空间(DOG)。

D(x，y，σ)＝(G(x，y，kσ)-G(x，y，σ))*I(x，y)＝L(x，y，kσ)-L(x，y，σ) (3-3)

式中：k＝1/s。即对各组相邻层高斯图像相减来获取高斯差分尺度空间。其示意图如图3所示。

随后的极值点检测便是在所构建的高斯差分金字塔图像上进行。如图5所示，将某一像素点与其同层和上下层相邻位置的26个像素点进行比较，若该点像素值为极值，则列为候选关键点，记录其尺度和位置信息。

2.精确特征点位置确定

为将所得的离散极值点精确化，需通过二次函数进行精确定位，同时剔除低对比度的候选点和不稳定的边缘响应点，以提高提取的稳定性和抗噪能力。

1)去除低对比度点

将尺度空间函数D(x，y，σ)在局部极值点(x₀，y₀，σ₀)处泰勒二阶展开，如式(3-4)所示：

式中：X＝(x，y，σ)为偏移量。对该式求导并令其为0，则可得到精确的极值位置偏移量

如式(3-5)所示：

当偏移量大于0.5，像素点更接近于其他的像素位置，需用插值方法进行替代。将该点带入式(3-4)可得：

设置阈值m＝0.04，若

则保留该点，则而反之。

2)消除边缘效应点

由于DoG函数所得的边缘极值点对噪声敏感，较不稳定，故而需对不稳定的边缘点进行剔除。根据该函数峰值沿边缘方向主曲率大，垂直于边缘的方向曲率小的特性，可利用Hessian矩阵计算主曲率来进行剔除。其中计算各特征点的hessian矩阵如式(3-7)所示：

特征点在两个方向上的主曲率，与hessian矩阵特征值成正比。设α为最大的特征值，β为仅此于α的特征值，令 α＝γβ,则有如式(3-8)：

两者比值表示两主曲率的差异程度，比值越大，即主曲率相差越大，此处取 γ＝10，计算

结果，若大于等于

则该点为边缘点，进行剔除。通过上述对特征点的精确定位和不稳定点的剔除可使该算法拥有抗噪性。

3.特征点主方向分配

确定特征点位置后，需利用其邻域内像素点信息为其分配主方向。对于(x，y)处的梯度值和梯度方向计算如式(3-9)：

式中：L为与特征点的尺度最接近的高斯图像，采用梯度方向直方图统计邻域内各像素的梯度方向。其范围为0到360度，每10 度一个柱体，共36个，将直方图中的峰值作为主方向，若仍有大于峰值80％的局部峰值，则可定位辅助方向，从而使得SIFT特征点有旋转不变性。

4.特征点描述子生成

获取特征点位置、尺度和方向信息后，为使其具有一定的亮度、仿射不变性等，需要建立描述子。此处利用邻域内各像素的梯度信息确定一个4*4*8＝128维特征向量。根据之前得到的主方向旋转图像，以该像素点为中心在尺度图像中选取一个16*16的窗口，划分为4*4个子区域，每个子区域包含4*4个窗口，而后统计各个子区域中8个方向上的梯度值绘制直方图，将其按照顺序排列所得的 128维特征向量即sift特征描述子。

通过上述步骤便可完成SIFT角点的提取。

SURF算法

相对于SIFT算法，SURF(speed up robust features)算法先利用了hessian矩阵确定候选点，而后进行非极大值抑制，降低了计算复杂度，其整体算法如下：

1.构建海森矩阵与高斯金字塔尺度空间

假设函数f(x，y)，则hessian矩阵H为该函数的偏导。对于图像上某个像素点的hessian矩阵可定义为如式(3-10)：

故而hessian矩阵的判别式为如式(3-11)：

判别式的值即hessian矩阵的特征值，可根据其结果正负来判断极值点是否存在。在SURF算法中，用图像像素值I(x，y)作为函数值f(x，y)，用二阶标准高斯函数进行滤波，通过特定核的卷积计算二阶偏导数来计算hessian矩阵的三个矩阵元素

L_xx、L_xy、L_yy，得hessian矩阵如式(3-12)：

为使得特征点具有尺度不变性，利用式(3-13)对图像进行高斯滤波。

式中：L(x，t)表示图像的状态函数，通过高斯核G(t)与图像函数I(x，t)在点x的卷积实现。g(t)为高斯函数， t为高斯方差。通过该方法可为图像中的像素点计算hessian矩阵的决定值来判别是否为特征点。其中，为平衡准确值与近似值间的误差引入随尺度变化的权值，则可用如式(3-14)表示hessian矩阵判别。

det(H_approx)＝D_xxD_yy-(0.9D_xy)² (3-14)

由上可知，hessian矩阵的获取需要进行高斯平滑滤波和二阶偏导数的求取。这两步可用一个模板卷积来进行处理，y方向上可用如图6所示的模板进行处理。

其中：左图为高斯平滑处理只在y方向上求取二阶偏导数的模板，右图为处理后的结果。同理，对于x和y方向上的混合模板可用如图7所示的模板进行处理：

上述方式可得一张近似hessian行列式图。在高斯金字塔的构建过程中，如图8所示，SIFT算法改变图像尺寸，且运算中反复使用高斯函数对子层进行平滑处理的方式来构建；而SURF算法则使原始图像保持不变，只改变滤波器大小，从而允许尺度空间多层图像同时被处理，不需要对图像进行二次采样，提高了算法效率。

2.利用非极大值抑制确定特征点的位置

该阶段与sift算法相同，将上述得到的像素点和它三维邻域内的26个像素点进行大小比较，保留极值点。而后利用线性插值得到亚像素级的特征点，同时去除小于一定阈值的特征点，以保证一定的抗噪性。

3.特征点的主方向确定

为保证旋转不变特性，SURF对特征点邻域内的harr小波特征进行统计。以特征点为中心，计算半径为6s(s为特征点所在的尺度值)的邻域内，统计60°扇形内所有特征点的水平和垂直harr(尺寸为4s)小波响应总和。并给以响应值高斯权重系数，使得越靠近特征点，响应贡献越大。而后将60°范围内的响应值相加形成新矢量，将圆区域内最长矢量的方向为特征点主方向。其示意图如图9所示：

4.构造SURF特征点描述子

SURF特征点描述子为64维向量。先在特征点邻域内取边长为20s的正方形框，方向为特征点主方向，再将其划分成4*4个子区域，每个子区域将存在25个像素点，对这25个像素点进行水平与垂直方向4个haar小波特征结果计算，分别为：水平方向之和∑dx, 水平方向绝对值之和∑|dx|，垂直方向之和∑dy，以及垂直方向绝对值之和∑|dy|。该过程如图10所示：

由于各子区域有4个值，则各正方形框有16*4＝64个值，故而各SURF特征点描述子可用64维向量表示，为SIFT描述子的一半，从而减少了计算时间，提升了效率。

两种算法对比与实验

在本节中通过SIFT与SURF算法对求是牌坊进行特征点提取，所得的图示结果如图11、12所示。

通过两图对比可以发现，SIFT角点在图像中的显示更加均匀，而SURF角点在图像中会大量集中在物体边缘处显示；在图像的特征点提取数量方面，如图13所示，SURF角点提取的数量明显大于SIFT角点；而在匹配点对的结果数量上看，如图14所示，SIFT角点的结果又稍稍比SURF角点好，由于本文所用的重建方法为离线重建，故而为能有更好地重建结果，将两种特征点提取结果联合考虑，来提供更丰富的图像信息。

图像特征点匹配

通过3.2中的方法得到特征点后，需对图像对进行特征点匹配，以获取图像间关系。而在匹配过程中，衡量特征点相似性的度量主要为欧式距离、汉明距离和余弦距离三种，而度量方式的选择一般取决于特征点描述子的类型。

改进的特征点匹配

由于SIFT与SURF特征点描述子不是通过二进制表达的，需利用欧氏距离度量两个特征点的相似性。而计算高维特征向量间的欧式距离需计算平方根，较为耗时，为提高匹配速度，本文采用层级哈希检索方式。同时，由于无序图像序列中，需对所有图像进行两两匹配，当两视图间的关联度不大，如拍摄角度变化、拍摄区域不同等情况下，两者间的匹配点较少时便无需进行全部特征点的匹配，便可利用预匹配来加快匹配速度。本文主要结合了预匹配与层级哈希检索方式来进行特征匹配，在一定程度上提升了匹配速度。

预匹配

预匹配过程主要判断两图像对的关联度。其主要步骤如下：

1、对图像对图I₁、图I₂提取角点_，分别得到M₁和M₂个特征点；

2、当匹配次数M₁×M₂小于阈值M，则直接进行两两匹配；否则进行预判断；

3、图I₁对图I₂匹配：取出SIFT角点M_SIFI个和SURF角点M_SURF个于图I₂中寻找对应匹配点，若匹配个数N_SIFT1和N_SURF1大于阈值n₁和n₂，则进行下一步匹配；否则跳过该图像对的匹配。

4、图I₂对图I₁匹配与3同理，先取出SIFT角点M_SIFI个，SURF角点M_SURF个于图I₁中寻找对应匹配点，若匹配个数N_SIFT2和N_SURF2大于阈值n₁和n₂，则需要对图I₁与图I₂进行完全匹配；否则跳过该图像对匹配。

通过该过程的预匹配，可在一定程度上加快图像对的匹配过程。

基于层级哈希检索的匹配

该方法主要利用了局部敏感哈希(LSH)的特性，将角点描述子转换为二进制表示，并对描述子进行聚类操作，来加快匹配速度。

其主要步骤如下：

1、图像特征描述子归一化

为对图像集I中所有图像进行统一关联，且在后续将描述子转换为二进制，此处对所有特征点描述子去均值归一化，将所有描述子的取值范围变为[-0.5,0.5]。由于本文利用了SIFT和SURF两种算法，故此处需将两者分别进行特征描述子的处理。

2、描述子二值化

先通过std：：normal_distribution用均值为0，方差为1的正态分布函数生成n×n维矩阵，而后将其与步骤1中归一化后的n维描述子向量[D₀，D₁，…，D_n]相乘，取得n维的向量值[C₀，C₁，…，C_n]，而后n位的二进制码中每一位d_i可由式(3-16)计算得到：

此处对于流程中n值，若为SIFT角点则取值128，若为SURF角点取值为64。以SIFT描述子为例，其二值化过程如图15所示。

3、哈希映射

为减少候选匹配的角点数量，仅与相似度高的特征点匹配，此处通过哈希映射将相似性高的特征点映射到一个桶中，从而使得待匹配点仅与落在同一个桶中的候选匹配点进行匹配。

1)利用std：：normal_distribution函数构建l×k×n维矩阵。其中：l表示构建的l张哈希表，k表示映射后的二进制编码位数为k位，即k个哈希函数。

2)矩阵相乘：将n维的二进制描述子与每张哈希表对应的k×n矩阵相乘得k维的向量值[C₀，C₁，…，C_k]。

3)二值化处理：可由式(3-15)，对向量值二值化，获得k位二进制编码。

4)桶分配：k位二进制编码对应2^k个桶，故而对于每张哈希表，一个特征点都将分配到2^k个桶中一个。一个特征点在 l张哈希表中，将被分配到l个桶中。最终只需取出与待匹配点落在同一桶中的候选匹配点进行匹配便可，从而减少匹配数量，加快匹配速度。以SIFT角点为例，其示意图如图16所示。

4、最近邻哈希排序

为加快匹配搜索速度，此处开辟长度为2^k的vector，在汉明空间中遍历计算候选点二进制描述子与待匹配点二进制描述子的距离，按照距离的大小存放到对应的vector空间中，再根据设定的阈值K，按距离从小到大依次取前K个邻近点；通过两者欧式距离的计算来确定最终描述子间的距离，取其中的距离最小值a与次小值b的两个邻近特征点，若两者比值c＝a/b大于一定的阈值 C则匹配成功，储存该匹配对，反之跳过。

实验结果

通过预匹配和哈希匹配在一定程度上可加速匹配，此处将其与普通的Brute匹配进行比较；实验中M设为1000000，M₁设为300，M₂设为300，n₁设为3，n₂设为3，l设为6，k设为8，对于n值，SIFT描述子为128维， surf描述子为64维。其相应的实验结果如图17、18、19所示，实验比较了基于哈希检索的方法与Brute匹配方法分别在有预匹配与无预匹配下的结果，从图中可知，哈希检索方法能在保证匹配所得tracks差别不大的情况下极大地提升匹配速度；而加入预匹配过程也在一定程度上加快了整体的匹配速度。如表3.1所示，为图像20张时的具体数据结果。

表3.1特征点匹配情况分析

图像特征匹配的校正

匹配中，由于存在环境光线亮度变化、拍摄角度变化和位置大跨度等以及相应计算方法的局限性，上述所得匹配点对仍有错误结果存在。而错误匹配将影响重建结果，故需消除尽可能消除误匹配点。此处主要利用以下约束与方法进行误匹配的消除：

1、一致性约束。场景中点可见且无遮挡时，图像对中的两图像特征点一一对应，即图1的特征点A对应图2中的B点，那么图 2中的B点也应对应图1中的A点。若匹配时两者并不对应，则去除该匹配对。

2、差异性约束。匹配时，点对的最近邻与次近邻的度量距离之比即c＝a/b需大于相应的阈值M。

3、随机采样一致性方法：通过连续迭代最小点集来估计目标函数的初始值，并以初始值为依据将数据分为内点与外点，在达到最大采样次数后取内点数最多的一次模型参数作为最终的结果进行数据分类。本文中通过随机采样两视图中的8组匹配对估计基础矩阵F，而后用sampson distance度量数据是否为外点，统计内点数量，取内点数最多的模型参数作为最终结果。

通过上述方法可得如表3.2所示的结果，图20显示了匹配对1中的情况，其中图20a为未消除误匹配的显示，图20b为消除后的图像显示，从两图中可见其中多对不正常的匹配对被剔除。故而通过以上约束与方法能较好地去除误匹配结果。

表3.2匹配结果情况

空间三维点注册与相机位姿估计

获取图像匹配对集后，可构造图像间的连接关系，如图21所示，从而通过计算本质矩阵、基础矩阵等估计相机运动，并注册空间三维点，但此处由于SIFT与SURF角点所得为稀疏特征点，故而注册所得为场景的稀疏点云。

图像特征点跟踪

图像特征点匹配节中所得匹配仅为两两图像的匹配，如图22所示；而对于多幅图像亦存在共同特征点，以三幅图像为例，如图 23所示，图I₁中的a点，图I₂中的b点与图I₃中的c点对应于空间中的同一点，三者可构成连接关系。对图像特征点跟踪即找到多视图间中的共同匹配点，以确定多个视图间的对应关系，如图24所示。

对于特征点的跟踪可利用其配对的传递性将所有共同特征点串联成一条条tracks。如图24所示，图I₁中的点

图 I₂中的点

两者匹配记为

同时图I₃中的点

与图I₂中的点

匹配，得

图I₄中的点

与图I₃中的点

匹配，得

则利用其传递性可得点

这4个特征点互为匹配点对，即表示空间中同一点。以此类推，通过特征点的互匹配关系串成各自的track，对所有特征点遍历后，可得所有点对的串联情况，从而获得多视图间的连接关系，实现整体特征点的跟踪。

对极几何与三角测量

对极几何

当已知两视图中的若干匹配点对，便可用对极几何约束估计出两帧之间的相机运动。

对极几何，用于表述两个相机帧的几何约束，即两相机图像中其对应点的视线延长线会交于空间中的一点，如图25所示。图中

O₁与O₂为两相机的光心点；M是空间中的三维点；m₁和m₂分别表示M投影在两个相机图像上的点，两者呈配对关系；e₁和e₂为极点，是光心连线与图像平面的交点；O₁与O₂的连线为基线；由M、O₁与O₂三点构成的平面为极平面；l₁和l₂为相机图像平面与极平面的交线；m₁的匹配点必定位于与之对应的外极线l₂上；同理m₂的匹配点也必然位于l₁上。

假设空间点M＝(x，y，z)^T，K₁，K₂为相机内参，R为相机旋转向量，t为相机平移向量，则点m₁和m₂的像素位置关系可表示为如式(3-16)：

m₁＝Rm₂+t (3-16)

等式两边同乘

得式(3-17)：

式中：t^表示反对称矩阵。整理可得：

令E＝t^R，称E为本质矩阵，上式可写成：

又因为

将两式代入式子可得：

整理可得：

令

称F为基础矩阵。式(3-19)和(3-21)便为对极几何约束。

而求解基础矩阵并结合相机内参可获取相机外参R、t。基础矩阵F的秩为2，7个自由度，常用解法为8点法、基于RANSAC的估计法等。本文主要结合8点法与RANSAC方法进行求解。

除基础矩阵与本质矩阵外，单应性矩阵H也描述了两平面间的映射关系，应用于相机发生纯旋转等情况，其亦可通过匹配点对求解。但相机纯旋转为退化状态，无平移则无法准确注册三维空间点。在初始化阶段需去除该情况。

三角测量

对极几何可估计两帧间的相机运动，对于空间点的注册需利用三角测量法。如图26所示：

三角测量，通过相机光心与图像点所成两视线交点来确定空间点位置。即图I₁中的点p₁与图I₂中的点p₂对应，直线O₁p₁与O₂p₂理论上交于空间中的P点，为两特征点在空间中映射点。但由于误差存在，两条不相交，故而可利用式(3-22)通过最小二乘法求解最优位置坐标。

x＝TX (3-22)

式中：x＝(u，v，1)为图像点的齐次坐标，X＝(x，y，z)为三维点坐标，P为3×4的相机外参矩阵。可将式(3-22) 分解为：

uP³X-P¹X＝0

vP³X-p²X＝0 (3-23)

其中：P₁、P₂、P₃分别为P的第1、2、3行，两组相机可得4个方程，从而线性求得X。为使X的求解结果更精准，此处利用最小二乘法求解。

利用三角测量法可得三维点的空间位置，但也存在不确定性。如图27所示：

当平移较小，由于像素噪声使像素偏移δx，视线角变化δθ，则测量所得深度会变化δd；但平移较大时，深度变化会较小。故而，同样分辨率下，平移越大测量结果越精确。但增大平移会引发遮挡，因此需在可匹配情况下尽可能使两图像间的平移距离最大化。

初始对选择

获取整体tracks信息后，便可利用对极几何和三角测量方法，进行初始重构。由于后续的重建恢复都基于初始重构进行，故而初始相机对的选择很大程度上影响最终重建结果。在初始相机对的选择上，需从以下几方面考虑：1)两视图间的匹配点对足够多：匹配点对数n大于阈值N；此处N设置为50。2)两视图所得单应性矩阵对应的内点数量少：满足单应性矩阵表明内点多位于同一平面或相机为纯旋转，属于退化情况对重建精度存在影响。此处设置内点数占比阈值为0.6；3)两视图间的平移足够大：平移足够大可使结果精确，其平移是否满足通过两视角与空间三维点所成的夹角确定。设置角度阈值为5°。4)可注册足够多三维空间点：注册点数为匹配对数的一半以上则可作为初始图像对。若以上条件都可满足，则为较好的初始对。若有不满足，则可设置相应的分数机制选取。

通过初始对的选择、对极几何的求解和三角测量的使用便可获得初始重构结果。

3D点-2D点运动估计

初始重构后，需不断添入新视角，从而逐步获得全部情况。在新视角添入的选择上，需选取与已有的视图间存在较多共同特征点对的。与初始重构中使用对极几何求解不同，新视角的运动估计时，由于已知的空间点位置与其在新视角上的投影位置，构成3D-2D的投影关系，则可通过PNP估计新视角运动。此处主要引用keip方法来进行求解。其主要步骤如下：

1.建立过渡相机坐标系τ

如图28所示，记原相机坐标系为υ，并假设空间中点P₁，P₂，P₃和相机内参已知，来获取矢量

通过相机坐标系υ中的

可建立过渡相机坐标系

其中

这样便可构建变换矩阵

特征点矢量亦可通过式(3-24)转换到过渡帧τ下。

2.建立过渡世界坐标系η

如图29所示，构建过渡世界坐标系η。该世界坐标系

其中

通过变换矩阵

可将原世界坐标系中的点转换到η坐标系下，

3.建立半平面并使用相应的参数求解两过渡坐标系变换

定义半平面∏，该平面包含点P₁、P₂、C，和单位向量

和

如图30所示：

其中：P₁，P₂，C点构成三角形，P₁，P₂点的间距d₁₂可知，

间的角度记作β，

由于β参数在使用中都以cotβ表示，故直接求取cotβ如式(3-25)：

式中：b表示cotβ。而后定义参数α∈[0，π]用以表示∠P₂P₁C，可得：

此处相机中心C在半平面∏中表示为：

故而τ在半平面∏中的三个基向量可表示为：

再定义自由参数θ表示半平面∏绕

的旋转角度，如图31所示：

其旋转矩阵如式：

其中，若

则θ∈[0，π]，则而反之。故相机中心C在η坐标系下表示为：

而从η变换到τ的变换矩阵如式(3-30)。由于

可得式(3-31)

假设

则两坐标系表示为如式(3-32)(3-33)(3-34)：

把(3-32)带入(3-34)中可得四阶多项式：

a₄·cos⁴θ+a₃·cos³θ+a₂cos²θ+a₁cosθ+a₀＝0 (3-35)

其中：

通过上述四阶多项式求解可得cotα值。而相机中心相对于世界参考坐标系的坐标和方向，可表示为如式(3-37)：

C＝P₁+N^T·C^η

R＝N^T·Q^T·T (3-37)

由此便可估得新添视角位姿，而后注册处空间三维点。

三维点优化

空间点注册与视角不断添入过程中，远处点在视角间的夹角小出，存在的误差大，需予以滤除；对于其余的相机位姿以及稀疏空间点，为能够获得更好的结果，需要对数据进行全局或局部的非线性优化，该过程称为捆绑调整优化。其主要通过求解相机矩阵T_k和三维点M_i使观测所得的图像点m_ki与其重投影图像点P_k(M_k)的均方距离最小化。对于m个相机和n个点而言可得如式(3-38)：

而后便可通过梯度下降法、高斯牛顿法和levenberg-Marquardt法等方法来进行非线性优化求解。本文调用了ceres库中的LM法求解。

以浙江大学求是牌坊为例，对其采集图像后通过上述步骤便可实现整体的结构重建，其重建结果如图32所示。从重建结果可见，在纹理较丰富的牌匾等位置处重建出的三维点较多，而在纹理较少平面出重建所得的三维点较少。

对称重复性结构的检测与矫正

图像匹配基本准确的情况时，上述重建流程能获取较好的结构重建结果，但当场景中存在重复性或对称性结构时，会对重建结果产生影响。如图33所示：图33a与图33b显示了求是大讲堂牌坊在不同角度拍摄的图像，由于两者主体图像基本相似，引发错误匹配使得注册结果错误，利用一般的重建流程所得结果如图34所示。

其图像整体采集本为两侧采集，但由于牌坊的正反两面的图像几近一致，按一般重建流程将使所有相机位姿和空间点仅恢复在一个侧面，产生错误。故而需要通过其他方法对其进行矫正。本文主要结合Heinly方法，利用图像间的背景差异点对重建结果进行矫正。如图33a、b所示，两图像间主体虽然相似，但其背景仍旧存在区别，如图33a左图中牌坊的右侧是树木，而在右图中牌坊的右侧则是建筑，故而可利用背景间的差异来检测图像是否一致，对整体进行矫正。相应的矫正流程主要包括图像连接图分割、重复性结构冲突测与组别点云融合三部分。

图像连接图分割

匹配过程中，图像间构成两两相关的连接图。如图21所示。但因相似性结构的存在使匹配错误，导致连接图中部分连接有误，故需对其进行分割，生成不同的子图进行检测；连接图分割时，若对每一条连接线都进行分割，则会产生大量组别，如若存在m个图像，则最多会有2^m-1种分割情况；为化简分割情况，此处利用最小生成树的方式，提取骨架进行分割。如图35所示。最小生成树生成过程中，以每幅图像为节点，以两图间三维点的观测情况为边的权重，若存在两图正确匹配，则会存在大量的共同点，故此处利用两者的非共同观测点作为边的权重进行处理，如式(3-39)所示：

式中e_ij表示图i和图j间的权重结果，O_i表示图i中可视三维点的集合。若两者匹配正确，则权重较小且趋近于0即为A边；若两者匹配错误，则权重较大；若匹配中有重复性结构存在，则权重介于两者之间，记为B边。理想情况下，可通过最小生成树由A边将正确结构的图集连为一组，而用B边连接存在重复性结构的组别。在断除B边后便得到相似性结构的图分割结果，而后对相似结构的图集进行进一步的重复性结构检测，以获取细分，取得较好的结果，从而减少分割情况。其分割如图35所示。

重复性结构冲突检测

对于B边的检测，需通过评估两子图间的冲突来判断，以确定是否分割，从而进一步减少分割的情况。冲突检测过程主要利用子图集中的独立点来进行判断；

对于场景空间点而言，可分为独立点和共同点；如子图集i和j，其独立点和共同观测点的划分如式(3-40)。

式中：C为两图集的共同观测部分，U₁和U₂分别为图集i和图集j的独立观测部分，独立点包括显著的背景不同点和重复结构中细节的不同点等。为提高鲁棒性，减少噪声影响，对图集中可观测到的三维点而言，需保证一个三维点能被共同观测到的数量γ要超过一定的阈值N，从而一定程度上保证检测的可靠性。

冲突检测的过程便是将各自图像中的独立点投影到图像对中的另一图像上；在图像对的选择上，为缓解遮挡或大角度变换等情况的影响，仅考虑视角变化较小的图像对。此处通过式(3-41)利用两视图光心位置C_i和C_j与其共同观测点平均中心

位置来计算视图间的夹角α。当α大于相应的阈值θ时，则计算该图像对，则而反之。此处设置阈值为15°。通过视图夹角的计算和阈值设定，亦可在检测出冲突的基础上避免了对分割中所有图像对进行检测的可能，在一定程度上减少了计算量。

满足视角要求后，可进行冲突检测。对于图像对i和j，将视图i中的独立点U₁投影到视图j中，若U₁中的点与视图j中的点有所重叠，则两图像在重建的过程中，由于重复性结构的存在导致相机运动估计结果错误。如图36所示，左图圈中的点通过投影变换将会与右图圈中的点产生搭接，从而检测出冲突，其中在检测投影点是否存在搭接上，主要利用了SLICO 方法来对图像进行整体划分，用于判断点是否在同一区域。但由于共同观测点附近的部分独立观测点，可能是由于遮挡或角度等原因产生，故而在冲突判断上需剔除与共同观测点C相近的独立观测点。

而在冲突度量上，利用图像对间投影所产生的搭接点的最小数量值t作为冲突的度量结果，其表达如式(3-42)：

N₁＝near(U₁，proj(U₂))，N₂＝near(U₂，proj(U₁))

t＝min(N₁，N₂)

其中：若某个点与大量点存在冲突，则将其设置为额外特殊点并剔除；τ为冲突阈值。对于每种分割情况而言，可通过计算该分割下两子图中所有匹配对情况的平均冲突值

将其作为该分割情况的冲突结果t；并将其与阈值τ比较来判断是否冲突；若冲突则断开连接分为两部分再进行进一步的分割检测，直到所有的冲突结果都小于τ。如图37所示，为通过最小生成树分割得到不同的组别的部分结果显示。

组别点云融合

获得各组无冲突结果的组别后需对其进行各组别的融合。融合流程中主要利用最初两两图像对匹配的点对结果；即在特征点匹配阶段获得的内点点对，由于特征点跟踪和视角添加环节被恢复到了其他位置而使得三角化时生成了额外三维点，让两者失去连接关系。故而需重新利用这些点对间的连接关系，通过两者生成的点云，用sim3求解相似变换矩阵来对其进行融合，从而获得最终的重建结果。

在相似变换矩阵求解中，主要通过两子图间各自的独立点，利用RANSAC方法来获取矩阵参数；而在内点判断环节，需将两子图集间包含的所有三维点都代入计算(此处包含共同观测点和独立观测点)，来获取内点个数，作为所获矩阵参数正确性的度量。

在对两两组别融合更新后，为进一步确定融合的准确性，需再次计算变换后两子图间的冲突值t，若变换后两子图集间的冲突值小于阈值，则融合正确。否者，需再次估计变换矩阵结果。获得正确空间三维点后，后续若再有新视角添入便可利用图像间的独立观测点来进行运动估计，从而获取正确的估计结果。

对比图38与图39c圈中结果，可见在一般流程下其所获结构重建结果只能恢复出一个面，也并没有恢复出另一个面的边缘等信息，所以其显示结果并不完整。而从图39中可见通过上述方法可将重建结果进行矫正，其较好地矫正了相机的运动估计结果和三维点的空间位置，恢复了正确的重建结果，从而可进一步进行后续的整体重建与应用。

基于图像序列的物体稠密表面重建及尺度恢复

三维重建的最终目的是实现场景或物体的整体重构，使其可辨认，可进行进一步应用。而通过增量式运动恢复结构的方法获取和矫正所得的三维点为稀疏三维点，结果不完整，且可视效果差，需将所得三维稀疏点扩展成稠密点，以获得对物体表面的更准确的描述和表达。

基于多尺度图像的稠密三维点构建

基于多视图的稠密三维点构建

该方法所使用的基于多视图的稠密重建算法流程如图40所示：

由增量式运动恢复结构的方式可获得稀疏三维场景点以及相应相机运动估计结果，方法将稀疏三维点作为种子点连同每个视角图像作为输入，而后通过全局与局部视图选择、空间区域生长以及多视图匹配优化处理来获得各视角深度图的最终输出。

视图选择

过程中，将待恢复深度的视角图像作为参考视图，记为R。而图集中将存在许多图像与该参考视图相互关联匹配，若所有相关图像都带入优化阶段将会使优化效率较低，且由于某些较差视角的存在会降低整体优化精度，故而需对视图进行选择。主要分为全局视图的选择与局部视图的选择：

3)全局视图选择

该过程主要从以下三方面考虑，从而保证稳定的匹配优化效果：

1)场景重叠度：通过比较待选视图V与参考视图R之间的共同特征点数量来选择共同观测点多、重叠度大的视图。

2)视差角：根据三角测量原理，足够的视差角可获得精度更高的重建结果。

3)尺度近似程度：视角中远近尺度不同亦会影响重建精度，因此需要选择物体尺度近似的视角作为候选。

为衡量视角是否作为全局候选视图，设置如式(4-1)的分数函数。

其中：N表示所选全局视图；F_X表示视图中可视三维点；ω_N(f)用以度量视差角情况，其表达式如式(4-2)：

其中：ω_α(f，V_i，V_j)＝min((α/α_max)²，1)；a表示两视角与三维点所成的视差角；阈值α_max设置为10°，用以剔除过小或过大的视角。

式(4-1)中：ω_s(f)用以度量两视图间的尺度关系。如图41所示，设像素点p(x，y)在相机坐标系下的位置为 p_c(X_c，Y_c，Z_c)，则一个像素点在空间中所占直径大小为Z_c/f；以S_R(f)和S_V(f)分别作为参考视图和待选视图的直径大小，用两者比值r＝S_R(f)/SV(f)作为权重函数通过下式(4-3)来衡量结果；

利用式(4-1)的分数函数可选择分数最高的视图集合N_G，作为全局视图。而在全局视角的处理上，由于不同视角的尺度不一致，会对立体匹配结果有所影响，故而需通过合适的滤波和降采样将尺度进行变换。此处主要选择全局视图中的最小分辨率作为衡量，将参考视图R下采样到与最小分辨率图像相近的尺度；若参考视图R为最小尺度_，则将高尺度的图像下采样到参考视图 R相近的分辨率大小进行匹配。整体尺度的计算上，需对待选视图V与参考视图R间共同特征点通过式(4-4)进行尺度统计确定最终的尺度。

式中：V_min表示

若scale_R(V_min)小于阈值t，则将参考视图R下采样到 scale_R(V_min)＝t的尺度。对scale_R(V)大于1.2的视图V都上采样到与参考视图分辨率近似的尺度。

4)局部视图选择

全局视图的选择，获得了与参考视图R相关性较高的匹配视图，缩小整体匹配范围；但对一个patch进行扩散时，为进一步得到更高的相关性匹配，提高立体匹配效率，用于深度图的估计和优化，需在全局视图中进一步筛选以获取局部视图。即选取其中的A张视图用于深度图的计算优化。选择中，主要针对待估计patch像素点的初始深度和法向选择光度一致和视差角较大的视图，直到数量满足A张或全局视图被挑选完毕。且对不同特征点patch进行扩散时，需选择不同的视角；同样，在深度估计优化过程中也需要调整视角，故需对A视图集合进行不断地更新迭代。

其中，以归一化互相关值(NCC值)度量视图R与视图V间的光度一致性关系。若NCC值大于一定的阈值，则可作为局部视图。

与此同时，再通过式(4-5)计算视图V与参考视图R的极线夹角来度量两视图间是否有较好的视差角。

l_R(V)＝g_R(V)·∏_V′∈Aω_e(V，V′V) (4-5)

式中：ω_e(V，V′)＝min(γ/γ_max，1),γ表示两视图的极线夹角，此处设置γ_max为10°。

全局视角的选择是针对某个参考视图R而言；而局部视图的选择，需要在patch扩散和深度估计优化过程中不断进行。

空间区域扩散生长

空间区域扩散生长，如图42所示，其将上章所得的注册点投影到参考视图R中，若三维点在视图R的可视锥内，则将其作为种子点，并以该点为中心设置n×n大小的patch。基于连续性与表面平滑性假设，相邻像素间的深度与法向值近似，可将种子点的深度与法向值作为邻近像素的初值。对于该patch，其中心种子点的空间位置可表示为如式(4-6)：

其中：O_R表示参考视图R的光心，

表示像素点的光线方向。假设patch中深度在水平与竖直方向上的深度线性变化值分别为h_s(s，t)和h_t(s，t)，如图43所示，则一个三维点位置在patch中的变化可表达为如式(4-7)所示：

式中：

在表面平滑且n足够小时，

故而邻近像素初值设置为

的值；同时对其中的像素点进行几何投影可获得局部视图中对应视角下的亚像素级光度值，可将其用于光度值一致性的优化；其相应的光度模型如式(4-8)：

I_R(s+i，t+j)＝c_k(s，t)·I_k(P_k(X_R(s+i，t+j))) (4-8)

式中：

c_k(s，t)表示视图R中的像素投影到其他邻近视图k中时所简化的反射效应系数；其中，每个patch为n²个像素，3通道彩色值，m张邻近视图，共可获得3n²m个方程，可解3m+3 个未知量：每张邻近视图h、h_s、h_t的3m个值以及每个彩色通道的ck₃个值。若忽略式(4-8)中的(s，t)并进行线性展开可得如式(4-9)。

该式在给定h、h_s、h_t初值后便可利用最小二乘法求解dh，dh_s，dh_t值，来不断更新h,h_s，ht值估计深度。在优化中，由于NCC值计算的复杂性，利用平方差的和值(SSD值)为光度一致性的度量来进行最小值优化，将NCC值作为最终的置信度结果。

扩散中，若邻近像素点仍未被处理，则将该像素的邻近像素点加入到优化队列中；若已存在初值，则选取置信度高的值作为最终值进行替换；直到扩散队列中所有像素点都优化完毕。但在优化过程中当重新计算到一个像素点亦或有些patch在优化过程中由于遮挡、光照变化等原因在某些局部视图中不可见等，则需对局部视图进行重新选择。为加速优化迭代过程加速收敛：需先迭代5次，计算每次迭代后的NCC值，若该值低于一定的阈值N_NCC则剔除该视图；若两次迭代后的NCC变化小于阈值τ_NCC，则优化收敛停止迭代；如此直到20次迭代结束；其中，每5次迭代中，仅优化深度信息；5次迭代结束或局部视图更换时，需同时优化深度、法向值和颜色反射系数；若14次迭代后，NCC变化值仍高于阈值τ_NCC，则重新选择视图，直到所有视图选择完毕。优化收敛后，将其NCC值作为其置信度保留。

由于该优化过程为非线性，良好的初值可避免陷入局部最优值，以获取较好的结果。但当物体为非平滑表面时，该方法效果较局限。

改进的稠密重建方法

基于上述的扩散生长方案可获得每个视角对应的深度图以及其转换后的稠密点云结果，如图44所示。由于该点云图是从各个视角的图像中获取所得，利用高分辨率图像可获得高精度的重建结果，但与此同时，随着分辨率的增大会使得其计算量和储存量上升；为在重建效果与效率间权衡，通过调整相应的算法流程，在保持一定重建完整性和细节的同时在一定程度上提高效率。

如图45和表4.1所示，其分别表示了不同尺度大小的图像通过稠密重建后所获得的结果。左侧为各视角的深度图结果，右侧为对应的稠密点云重建结果。其相应的图像分辨率大小分别为1920×1080，960×540；从图表中可见：对于高分辨图像 45a，虽能得到精确度高且纹理丰富的三维信息，得到的扩散像素点数量更多，但其耗费时间更大，且由于丰富的纹理细节存在导致重建的完整性较差，且其噪声显示的波动较明显。而对于低分辨率图，如图45b所示，其整体重建完整性高于高分辨率图像结果，且重建速度更快，但由于分辨率低使得纹理较为模糊，且其重建最终结果中一个像素所代表的尺度信息变大，从而会使测量的结果精度变低，波动性更大。

表4.1不同尺度图像重建结果

由于低分辨率图像重建所得结果较为完整，此处将低分辨率结果与高分辨率结果进行联合优化，将低分辨率图像构建的深度图结果作为高分辨率图像的初值进行深度优化，来获取最终的高分辨率深度图结果，其主要步骤如下：

1、选取合适的图像层级L对稀疏种子点点集P进行扩散生长，获取点云p′，并将其投影到图像中构造像素集M。

2、取高分辨率图像层级L＝L-1，并对构造的像素集M进行置信度等插值拟合和邻近扩散，从而构造像素集M′。

3、对图像层级L按式(4-1)～(4-5)中的方法进行视图选取，对不同的patch获取其置信度值T。

4、当所得置信度值大于阈值τ，则根据patch深度信息进行深度值的偏移调整；当小于阈值τ，则需将其加入优化队列Q进行进一步优化。

5、根据式(4-9)进行光度一致性约束，对整体优化队列进行优化获取level L层的参考视图深度图结果。

虽然该算法在理论上可从最低分辨率向上提升到最高分辨率，但由于步骤2中的插值拟合和后期的一致性优化过程中由于低分辨率纹理模糊严重等原因，在向上扩散时易产生偏移较大点，从而影响整体的模型结构。但对于一些平面区域较多且面积较大的牌坊而言，可以适当从较低分辨率向上提升至较高分辨率，其能在重建完整性和效率上有较好的结果。而在牌坊场景的恢复中主要以1层级作为初始层级进行，并取阈值τ为0.8。如图46所示，便是以低分辨率结果为初值恢复高分辨图像模型的整体结果示意图，如表4.2所示为相应的结果显示，从结果图中可见，对于混合结果(图46c)，与原先的的高分辨率结果(图46a)比较，其混合结果在整体的完整性和重建的计算效率上有相应的提高；而与低分辨率结果(图46b)比较，其在细节显示和整体精度上也会有所提升，这样对于后续的尺寸量取，也能得到较好的量取结果。

表4.2不同尺度图像重建结果

由于本文的重建流程是对每个视图都先构建出深度图，后续还需根据两相机外参间的关系对其深度图进行融合，即将不同相机坐标系下的点云进行坐标变换，而后融合成稠密的三维点云。

稠密点云网格化处理

为更好地表述三维模型信息，需对稠密点云进行网格化处理。此处基于八叉树原理对三维空间进行划分，将上述所得的三维稠密点云进行部分融合，而后利用隐式曲面重构法对点云数据进行网格处理和表面提取，以获得最终的三维表面模型。

八叉树构建

八叉树是一种用于描述三维空间的树状数据结构，在三维空间数据处理中应用广泛。该树中，每个节点表示一个正方体的体积元素；每个节点有八个子节点，父节点的体积便为八个子节点体积元素的叠加；每个子节点仍可不断往下细分八个子节点，从而将整个模型空间划分为多个子空间，并表示为树的根节点与其八个子节点的连接，该树状存储结构便为八叉树，如图47所示。

该八叉树结构可应用于点云快速最近邻查找、点云压缩或碰撞检测等，其主要实现原理为：1、设定最大递归深度，即设定最小立方体体素大小。2、获取场景的最大尺寸，建立根节点体素。3、依序将点云存入能被包含且没有子节点的节点体素中。4、若未达最大递归深度，则细分叶子节点体素，将该节点所含点云元素分别存储到八个子立方体中。5、若子立方体所分得的点云元素数量不为零且与父立方体一致，则停止细分该子立方体，以避免无限切割的情形。6、重复步骤3，直到最大递归深度。

通过如上步骤便可利用八叉树来对三维点云进行管理。本节中为避免隐式曲面重构中由于元素距离表面太远而出现混淆现象，对每个体素都根据存入的点云元素尺度设置了相应的体积大小。对于每个点云采样点i都存在对应尺度值s_i，则可将该采样点的空间半径定为3s_i，再利用式(4-10)对空间进行划分。

式中：l表示采样点在八叉树中的层深，S_l表示l层级下的节点所具有的边长，其主要划分流程如图48所示。

读入采样点i，若为第一个采样点，则确立空间中心点p_i，并设置s_i的边长。若非第一个采样点，则根据式(4-10) 判断是否在现有树内部；若在内部：则将该点插入到八叉树内，且当尺度s_i＜S_l时需继续对八叉树进行细化，直到该采样点插入；若在外部：则需对八叉树进行扩展；即根据采样点空间位置与原八叉树根节点中心位置不断调整新根节点中心位置及其边长值，直到根节点包含所有采样点。再不断读入采样点，直至所有采样点都被存入八叉树结构中，则输出最终结果。

符号距离函数值获取

获得八叉树储存结构后，需对其进行网格化和表面提取处理以重构物体表面。与泊松重建方式相同，本文也使用了隐式曲面重构法。为获取如图49所示的零平面，构建带有符号距离隐函数，如式(4-11)：

式中：f_i为基础函数和ω_i为权重函数；其都与采样点i的空间位置p_i、法向量n_i和尺度s_i等参数相关。利用式(4-12)可将树中立方体八个角点的坐标点位置进行转换到采样点的局部坐标系下，让采样点处于原点位置，使其存在一致的正向x轴方向。

x_i＝R_i·(x-p_i) (4-12)

式(4-11)中基础函数f_i的表达式如式(4-13)：

式中：σ＝s_i，即采样点的尺度。其图像示意图如图50所示，仅在x轴方向上存在正负变换，而在y、z轴上都不存在正负变换。故而可将物体表面内部设为负值，外部为正值，以使物体表面为零平面。

而式(4-11)中权重函数ω_i的表达式如式(4-14)：

该函数图像示意图如图51所示，其可使符号距离值近零平面时更平滑，且外部值的平滑效果大于内部值。

通过构建如式(4-11)的隐函数，便可估计八叉树中体素角点的符号距离值。其主要流程如下：

1、计算角点坐标值：对每个八叉树的叶子节点所代表的体素计算其八个角点的坐标位置。

2、筛选采样点：由式(4-13)与(4-14)可知，隐函数值估计时，若与采样点距离过远，则其对计算结果无太大影响，故而忽略过远点的采样点的计算，其判别标准如式(4-17)，以体素角点与采样点所在的体素中心点距离以及边长的关系判断。

3、局部坐标系转换：为获取统一结果，将体素角点转换到采样点的局部坐标系下，进行隐函数值计算。

4、隐函数值计算：利用式(4-11)计算角点的隐函数值。若同时计算法向量，则需对式(4-11)求导以获取方向向量。

5、还原变换：在获得隐函数结果后，对于其法向量，则需将其转换回角点所在坐标系。

通过以上步骤便可获取相应的符号距离值结果。

移动立方体的等值面提取原理与改进

在计算得到各体素角点的隐函数值后，便需进行零等值面提取，来进行网格化处理。一般规则网格的网格化处理，可利用普通的移动立方体(Marching Cube)算法便可获得较好的结果。而八叉树划分的空间中，由于八叉树中各叶子节点所处深度层级不一，利用Marching Cube进行等值面提取时会产生二义性，如图52所示。

为消除该二义性，需不断检测该叶子结点的兄弟节点中是否存在细化的子节点；若存在细化子节点，则需获取细化面或细化边中的结果来获取最终的等值面。其主要表现如下：

一般Marching Cube算法中，由于八叉树中节点尺度不同导致相邻面的等值线划分粗细不一致。如图53所示，立方体M₁与立方体M₂为八叉树中两个互为兄弟节点的体素，而小立方体m₁为父体素M₂中的一个子体素，其中面A₄A₅A₆A₈为M₁和M₂两体素的相邻面，故而M₂中子体素m₁的面a₁a₂a₃a₅亦与M₁的面A₄A₅A₆A₈相邻，但两者在八叉树中属于不同的尺度。对于体素M₁而言，其面A₄A₅A₆A₈的等值线划分结果为p₁p₂，而在体素M₂中由于细化体素m₁的存在使得其对等值线进行了更细的划分，导致两个相邻面的等值线出现歧义。

为消除以上歧义现象，在进行等值面提取时需要进行如下步骤的判断：如图54所示，计算体素Ⅰ的等值线时，需判断其兄弟体素Ⅱ的相邻面是否存在细化；若存在细化，则需获取体素Ⅱ中细化后的等值结果如图53右图所示；而后将细化结果替换至体素Ⅰ中；若不存在细化，则可直接进行等值面提取。同理，在等值点提取过程中，也需要考虑细化情况，对于一条边需考虑共边的三个体素，对于一个面则需考虑共面的一个体素，从而使得最终结果图如图55所示。

由上述网格划分及等值面提取后，需再消除如图56中所示的退化三角网格，使网格结果更平滑简略。对于牌坊重建后的局部网格化结果如图57所示，在网格划分时进行邻近面的判断可使整体平面更平滑；而没有判断邻近面时，会让网格中有一些较大的空洞等，如图圈中所示。而通过图像序列的输入，三维空间稀疏点的注册，相机位姿的估计，稠密点云的恢复以及点云的网格化等步骤处理后，便可获得整体的重建结果。

基于双目立体匹配的尺度恢复

就单目相机重建而言，由于图像在三角测量阶段对尺度进行了统一，使得所得三维模型失去尺度信息；为使重建模型能准确反映实际物体尺寸，需重建模型的尺寸信息，即重建模型尺寸与实际物体尺寸间的缩放比。对于该缩放比，其可通过如下方法获取：

1、利用已知尺寸的标志物。如图58所示，可将一个或多个已知尺寸信息的标志物放置于待重建物体周围，使得该标志物能与目标物体同时重建；而后量取重建后标志物模型尺寸长度，即失去尺度信息后的长度尺寸；则整体重建模型与实际物体的三维尺寸缩放比可由该标志物长度尺寸与其实际尺寸比值获取。而后便可利用该缩放比对重建模型进行缩放，恢复目标物体模型的实际尺寸，便可通过模型的尺寸量取获得实际的物体尺寸长度。

2、利用差分GPS设备或高性能IMU设备。如图59所示，利用差分GPS设备记录每个图像拍摄点的高精度实际位置信息，再与计算所得的各个图像拍摄点的位置信息相联合，通过如式(4-18)的相似矩阵变换，便可获得缩放比例，对重建模型进行尺度恢复。

B＝s*R*A+T (4-18)

以上两种方法均能获取重建模型与实际物体间的尺寸缩放比例，对重建模型进行整体比例缩放，获取具有实际尺寸信息的三维模型。但方法1需人为地先放置已知尺寸标志物，若待重建物在高处等，则由于场地因素等限制难以直接放置；而方法二虽能获取最终的缩放比例，但差分GPS或高性能IMU设备价格高昂。故而为获得单目重建模型与实际物体的缩放比例，此处对比结合了双目相机的简便性、价格低廉、对于局部尺寸有较好精度的特性，以利用双目相机来进行尺度恢复。

基于半全局的立体匹配方法的视差估计

双目相机的深度图获取主要利用了立体匹配算法，而在传统立体匹配算法中，全局匹配算法花费时间较长精度较高完整性较好，局部匹配算法的效率快但精度较低完整性较差；综合考虑，此处主要选用精度较好效率和完整性也较高的半全局匹配算法(SGM)，其主要流程为代价计算、代价聚合、视差计算、视差细化4步。

匹配代价计算

SGM算法的匹配代价主要通过互信息计算。而互信息结果主要表示两图像间的相关性，该值越大，则两图像间越相关。由于互信息计算匹配代价对环境光不敏感，故而与基于灰度值计算的方法相比其更具鲁棒性；互信息，通过熵来定义。熵表示了系统的混乱程度，其值表现了图像中包含的信息量多少，可通过如式(4-19)计算：

式中：P_I表示图像灰度值的概率分布结果。则两图间的联合熵可通过如式(4-20)计算：

式中：

表示两图像的灰度值联合概率分布结果，其如式(4-21)：

式中：若函数T[·]中的等式成立时，则函数值为1，否则为0。通过熵的定义可用如式(4-22)表示互信息：

式中：

与

分别表示左右图的熵，

表示两图的联合熵。

由于式(4-22)需将视差图作为先验，故需利用近似计算方法来获取联合熵，其式如(4-23)：

式中：

为数据项，其计算如式(4-24)：

图60所示为

的计算过程：其先根据初始的视差结果图，对右图进行平移操作，使左右两图中的匹配点尽可能位于相同位置；而后由式(4-21)计算两图的联合概率分布；并用7×7的patch对图像进行高斯平滑滤波；最后取负对数得

同样，单张图的熵

可由如式(4-25)计算获得：

式中：

便为数据项，由式(4-26)计算所得。

由式(4-21)已求得

故可由式(4-27)获取单张图的灰度概率分布。

将式(4-23)与式(4-25)，代入式(4-22)可得如式(4-28)：

式中：

表示为如式(4-29)：

故而，基于互信息的像素点P，可由式(4-30)计算其匹配代价值。

式中：I_b为_左图；I_m为右图；视差由d表示；

表示对右图的修正；q表示左图中的像素P在右图中的可能匹配点，其通过式(4-31)定义。

q＝e_bm(P，d) (4-31)

式中：e_bm为极线约束的表达。

其中，初始视差图为随机生成所得，后续可通过不断地迭代优化，来获得较高精度的视差结果图。且可利用分级求解来提高计算效率。

代价聚合

由于逐像素进行匹配代价易受光照、噪声等因素影响，易造成误匹配。故而在视差计算时，需根据当前像素的邻近像素视差结果来构建代价函数，以增加平滑一致性约束，即构建如式(4-31)的全局能量函数E(D)来综合考量匹配代价。

式中：第一项表示像素P在视差为D时的匹配代价和；后两项表示P点邻域N_P内的像素q进行的视差平滑性约束项：当P与q两者的视差差值较小，且小于等于1时，该项惩罚常数为P₁；当两者视差差值大于1时，惩罚系数为P₂。根据该平滑性约束，可对图像进行平滑的同时较好的保留其边缘。

视差计算

通过对上述全局能量函数E(D)的最小化便可获得最优视差结果。即若想要像素点P的匹配代价结果最小，则需使点 P邻域N_P内的q点匹配代价最小化；同样，若想q点的代价结果最小，则需q邻域内的n点代价最小。这样该问题便可转换为NP完全问题，但由于在二维图像上不能直接地利用动态规划方法进行求解；若直接按行列求解计算，其随能避免算法复杂度过大问题，但忽略了行间约束，会使视差估计结果较差。故而在计算P点的匹配代价聚合结果时，可利用P点邻域在16个方向上搜索得到的不同代价结果构建一维动态规划问题，来优化求解最终的视差结果值，如图61所示。

则像素点P在视差为D时的聚合匹配代价值计算如式(4-32)：

S(P，d)＝∑_rL_r(P，d) (4-32)

式中：r表示搜索方向；L_r(P，d)表示r方向上的聚合匹配代价结果，其计算如式(4-33)；最终像素点p的聚合代价结果由各方向上的代价和表示。

L_r(P-r，d+1)，minL_r(p-r，i)+P₂) (4-33)

式中：第一项为像素点P在视差为d时的匹配代价，第二项为p-r处的像素点在不同视差值下聚合代价最小值。

视差细化

由上述完成视差计算后，需对视差图进行亚像素插值、邻域一致性检测等后处理，使得获取结果更精确。

算法结果

如图62所示，为对牌坊进行立体匹配的结果。图62a为用于双目立体匹配的RGB图；图62b为通过SGM算法所得的视差图；图 62c是将视差图结合双目基线和相应的相机内参通过几何变换所得的局部牌坊三维点云图。由图中可见，其对于正对的牌匾有较好的重建结果，可用于尺度的比例恢复。

基于视差图的三维重建及尺度恢复

通过基于无序图像序列的三维重建方法可获取整体的三维物体模型，而基于双目立体匹配方法所得仅为局部物体重建结果。现可将两者结合，即通过局部物体的实际三维尺寸与重建所得的整体三维模型尺寸进行比较获取该整体模型恢复至实际尺寸所需的缩放比例S。

如图63所示，左图中的立方体Ⅰ为右图中立方体Ⅱ的子块；其中，Ⅰ中的A，B两点对应立方体Ⅱ中的A′，B′ 两点。但Ⅰ中AB的长度测量结果为M，对应于双目所得的物体实际距离；Ⅱ中的A′B′长度为m，对应于重建所得的缺失尺度的三维模型；则M＝S×m，即两者相差一个尺度比例

当获得S尺度后，便可对立方体Ⅱ进行尺度缩放，而后便可量取边CD、CA′等尺寸的实际长度。

其相应的实验结果如图64所示。其分别显示了双目点云的尺寸测量结果和通过图像序列重建所得三维模型的长度量取结果，对于牌匾的长度分别为270.8mm和0.7988；则其缩放比例S＝270.8/0.7988＝339.01。利用该比例S对整体模型进行缩放处理，而后便可获取整体中的其他尺寸长度。如图65所示，为缩放后的三维模型量取其他尺寸长度的结果。牌匾中绿色框量取所得的长度为190.4mm，而其实际长度为194mm，两者基本相似结果误差为1.8％。故而该方法对于较大构件而言，可较为准确地获得缩放比例，从而对物体进行整体缩放，而后用于量取其他部分或整体的尺寸长度。

Claims

1.基于图像序列的大型构件三维重建方法，包括以下步骤：

2.如权利要求1所述的方法，其特征在于：步骤5中将低分辨率的深度图结果作为高分辨率图像的初值进行深度优化，来获取最终的高分辨率的稠密点云，其主要步骤如下：

5.1、选取合适的图像层级L对稀疏种子点点集P进行扩散生长，获取点云p’，并将其投影到图像中构造像素集M；

5.2、取高分辨率图像层级L＝L-1，并对构造的像素集M进行置信度等插值拟合和邻近扩散，从而构造像素集M′；

5.3、对图像层级L进行视图选取，对不同的patch获取其置信度值T；

5.4、当所得置信度值大于阈值τ，则根据patch深度信息进行深度值的偏移调整；当小于阈值τ，则需将其加入优化队列Q进行进一步优化；

3.如权利要求1所述的方法，其特征在于：基于多视图的稠密三维点构建方法为：将稀疏三维点作为种子点连同每个视角图像作为输入，而后通过全局与局部视图选择、空间区域生长以及多视图匹配优化处理来获得各视角深度图的最终输出。

4.如权利要求1所述的方法，其特征在于：基于双目立体匹配，对三维重建结果进行尺度恢复，执行的操作为：利用双目相机拍摄目标物体的局部，而后对物体进行局部重建，以获得局部的真实尺寸信息；再相对应地量取整体重建所得模型中对应的局部尺寸长度，从而获得尺寸比例来对整体重建模型进行缩放，恢复整体模型的真实尺度，后续便可直接对整体模型进行尺寸量取来获得真实的尺寸长度。