CN103914874B - 一种无特征提取的紧致sfm三维重建方法 - Google Patents

Abstract

本发明公开了一种无特征提取的紧致SFM三维重建方法，包括：输入关于某场景的n幅图像，n≥2；建立与某个相机坐标系相一致的世界坐标系；以三维场景的深度和相机投影矩阵作为变量，构造类似光流估计的目标函数，采用由粗到细的金字塔方法，设计迭代算法对目标函数进行优化，输出表示场景三维信息的深度和代表相机相对位姿信息的相机投影矩阵；根据表示场景三维信息的深度，实现紧致的射影、相似或者欧几里德重建。本发明能够一步完成紧致SFM三维重建。由于通过一步优化实现紧致三维信息的估计，以目标函数值作为指标，能够得到最优解，至少是局部最优解，比现有方法有很大改进，已初步得到实验验证。

Description

一种无特征提取的紧致SFM三维重建方法

技术领域

本发明涉及图像三维重建领域，更具体地，涉及一种无特征提取的紧致SFM三维重建方法。

背景技术

基于计算机视觉的三维重建是指利用数码相机或者摄像机获取图像，构建算法以估计所拍摄场景或者目标的三维信息，实现表达三维客观世界的目的，其应用范围包括机器人导航、汽车自动或者辅助驾驶、虚拟现实、数字媒体创作、计算机动画、基于图像的绘制(image-based rendering)和文化遗产的保存等。

基于运动的建模(Structure from Motion，SFM)是目前常用的一种三维重建方法，即从两幅、多幅图像或者视频估计场景或者目标的三维信息。已有实现SFM三维重建的技术手段有以下特点：基于特征点的、稀疏的和分两步完成。已有SFM三维重建分两步完成：首先从图像检测并匹配具有尺度或者仿射等不变性(invariance)的特征点，包括Harris特征点、Kanade-Lukas-Tomasi(KLT)特征和Lowe尺度不变特征(scale invariant featuretransform，SIFT)，然后估计所检测特征量的三维信息和相机的姿态(包括位置和角度)。

已有的SFM三维重建算法分为两步完成，不能真正达到最优化效果。由于从图像中检测到特征点的二维坐标有误差，在其基础上即使采用优化算法重建其三维信息，也无法获得全局意义上的优化结果。由于特征点的匹配精度通常比较低，因此不可避免的造成低精度的三维重建。

三维重建效果是稀疏的(sparse)；由于只对所提取特征点估计其三维信息，不能实现紧致的(dense)三维重建，即不能估计出所有像素点的三维信息。对于30万像素的480*640图像，在保证一定正确匹配率的前提下，通常只能检测到200～300个甚至更少的特征点，相对于30万像素的图像来说，特征点是非常稀疏的，绝大部分像素都没有直接估计其三维信息。虽然可以进一步在特征点基础上，利用估计出的极线约束(epipolar constraint)等技术手段进一步估计其它点的三维信息，实现紧致或者半紧致(quasi dense)重建，但是由于所估计的特征点的三维信息和相机姿态存在一定的误差，影响后续其它点的三维估计效果。

发明内容

为了克服现有技术的不足，本发明提出一种无特征提取的紧致SFM三维重建方法。采用该SFM三维重建方法，不需要特征点检测及匹配，采用一步优化即可实现紧致三维重建。

为了实现上述目的，本发明的技术方案为：

一种无特征提取的紧致SFM三维重建方法，包括以下步骤：

S1.输入关于某场景的n幅图像，n≥2；

S2.建立与某个相机坐标系相一致的世界坐标系，设世界坐标系与第一相机的坐标系相一致，即世界坐标系的原点、x轴和y轴与第一相机的相机中心、第一相机成像平面的x轴和y轴重合，其z轴垂直指向第一相机的成像平面；

S3.以三维场景的深度和相机投影矩阵作为变量，所述三维场景的深度是指第1幅图像像素点对应的三维空间点具有的深度q；所述相机投影矩阵是指其它(n-1)幅图像的3×4矩阵P_i，2≤i≤n；

S4.构造类似光流估计的目标函数，所述目标函数是连续域上的变分目标函数或其离散形式的目标函数；

S5.采用由粗到细的金字塔方法，在连续域或者离散域上设计迭代算法对目标函数进行优化，输出表示场景三维信息的深度和代表相机相对位姿信息的相机投影矩阵；

S6.根据表示场景三维信息的深度，实现紧致的射影、相似或者欧几里德重建。

该方法能够一步完成SFM三维重建，由于通过一步优化实现三维信息的估计，以目标函数值作为指标，能够得到最优解，至少是局部最优解，比已有方法有很大改进，且已初步得到实验验证。

上述的相机是指某幅图像对应的相机，在本发明中，关于场景的第一幅图像对应的相机是第一相机，第一相机的坐标系与世界坐标系相一致；各幅图像均对应一个3×4相机投影矩阵。

在本发明中采用该类方式建立世界坐标系是为了计算方便，在实际中，可以任意建立世界坐标系，如果任意建立坐标系，则待估计的参数包括n个相机投影矩阵，刻画每个三维点需要三个坐标参数。在本发明虽然没有给出这种方案的技术细节，但是任意建立世界坐标系的方案与上述建立世界坐标系的方案基本原理一样。

为了实现射影三维重建，则首先要进行参数化设定，即在实现射影三维重建中，参数化具体为：在建立世界坐标系的同时，其第一相机的相机投影矩阵为[I₃ 0]∈R^3，4，其中I₃是一个3×3的单位阵，0是一个3×1的零向量；其它相机投影矩阵P_i∈R^3，4，2≤i≤n，作为待估计的未知参数；场景的三维结构由定义在第一副图像上的三维场景的深度决定：假设与第一幅图像像素(x,y)相对应的三维空间点的三维场景的深度为q_x,y，则该三维点的三维坐标为

(q_x,y×x,q_x,y×y,q_x,y) (1)

在射影三维重建中，相机投影矩阵P_i和三维场景的深度q_x,y作为待估计的未定参数，为了表达式的简练，在不造成误解的情况下，省略下标x,y。

实现连续域上射影三维重建的具体实现过程为：

构造的连续域上的目标函数具体为：

f(P₂,…,P_n,q)＝f_dats+f_{smooth_uv}+f_{smooth_depth} (2)

其中

对上述目标函数的说明如下：(a)为梯度算子， 为拉普拉斯算子；(b)目标函数分为三部分，数据项f_data，偏移平滑项f_{smooth_uv}和深度平滑项f_{smooth_depth}，其中α、β、τ₁和τ₂是非负权重；(c)图像有k个色彩分量C₁,…,C_k，代表第一幅图像在位置(x,y)的色彩I分量值，相应地，为第i幅图像在位置(uⁱ,vⁱ)的色彩I分量值；(d)鲁棒函数ρ的引入是为了克服深度发生剧变带来的影响，鲁棒函数ρ为Charbonnier函数其中∈是一个足够小的的正数，ε＜10^-6；或者为Lorentzian函数 σ为某个常数；当不引入鲁棒函数，则ρ(x)＝x；(e)uⁱ和vⁱ是定义在图像域上、以相机投影矩阵P_i和深度q为参数的函数：和代表与第一幅图像像素(x,y)相对应的三维点在第i幅图像的成像位置

其中P_i,j为第i个相机投影矩阵P_i的第j个行向量；为了表达式的简练，在不造成误解的情况下，在和中省略下标x,y；

在连续域上设计的迭代优化算法具体为：因为三维场景的深度是定义在第一幅图像域上的连续函数，在极值点必须满足欧拉-拉格朗日方程；同时，在极值点对相机投影矩阵的偏导数为0；在图像的离散格点上，联合欧拉-拉格朗日方程和对相机投影矩阵偏导数为0的两类方程，并采用增量方式表示形式，能够把求解相机投影矩阵和三维场景的深度增量的迭代过程转化为求解如下线性方程组

Hδθ+b＝0 (6)其中向量θ由相机投影矩阵P_i2≤i≤n和三维场景的深度q按次序构造而成；这样，每次迭代归结为求解

δθ＝-H^-1b (7)，

从而确定相应的增量δP_i和δq；根据所求解的增量更新参数P_i和q，P_i←δP_i+P_i，q←δq+q，直到收敛；

即算法1的具体过程为：

输入：n幅图像，初始化三维场景的深度q和相机投影矩阵P_i，2≤i≤n；

输出：相机投影矩阵P_i(2≤i≤n_i)、三维场景的深度q和场景的三维表示；

1、迭代

1)、由欧拉-拉格朗日方程和目标函数对相机投影矩阵的偏导数为0确定式子(7)中的H和b；

2)、由式子(7)计算增量δθ，并确定相应的增量δP_i和δq；

3)、更新参数P_i，2≤i≤n和q：P_i←δP_i+P_i，q←δq+q；直到收敛

2、根据收敛后的三维场景的深度q，由式子(1)计算场景的三维表示。

上述彩色图像，可以采用常见的RGB或者HSV表示，以RGB格式为例，图像有三个分量，分别为红色(R)、绿色(G)和蓝色(B)分量；色彩分量可以是不同格式的组合，比如采用R、G、B和H四个分量。上述鲁棒函数有很多种选择并不限于上述所列出的函数。

上述式(1)中构造的目标函数f(P,q)的依据在一定程度上与光流计算类似，即灰度不变假设和像素偏移量uⁱ-x和vⁱ-y的平滑假设，分别对应目标函数中的第一部分和第二部分即数据项和平滑项。式子(1)中的第三部分对应深度的平滑假设。

更进一步的，所述目标函数中的数据项和偏移平滑项能够采用其它类似的变化形式：

其中和

鲁棒函数的引入还能够以其它变化形式出现，式子(3)数据项的另外一种变化形式为：

与连续域上情形相似，构造离散形式的目标函数具体为：

f(P₂,…,P_n,q)＝f_data+f_{smooth_uv}+f_{smooth_depth} (11)

其中

离散目标函数(11)及其变化形式的迭代优化算法具体为：离散形式的目标函数(11)在本质上是一个非线性最小二乘问题，能够采用常规的Levenberg-Marquadt算法或高斯-牛顿算法，每次迭代过程归结为求解一个线性方程组(15)：

δθ＝-(H+μI)^-1b (15)

其中H是海森矩阵或者高斯-牛顿海森矩阵，b是梯度向量，μ是非负数，取决于采用Levenberg-Marquadt算法或高斯-牛顿算法，从而确定相应的增量δP_i和δq；根据增量更新参数P_i和q，P_i←δP_i+P_i，q←δq+q，直到收敛；

算法2具体实现过程：

输入：多幅图像，以及相机投影矩阵P_i和三维场景的深度q的初始化，2≤i≤n；

输出：相机投影矩阵P_i(2≤i≤n_i)、三维场景的深度q和场景的三维表示；

1、迭代

1)、计算式子(15)中的高斯-牛顿海森矩阵H和梯度向量b；

2)、由式子(15)计算增量δθ，并分别确定相应的增量δP_i和δq；

3)、更新参数P_i和q：P_i←δP_i+P_i，q←δq+q，2≤i≤n；直到收敛；

2、根据收敛后的三维场景的深度q，由式子(1)计算场景的三维表示。

更进一步的，由粗到细的金字塔方法步骤具体为：计算图像的n层金字塔表示；在最粗图像层，初始化n-1个相机投影矩阵为P_i＝[I₃ 0]，2≤i≤n，所有点的三维场景的深度初始化为1；由粗到细依次估计相机投影矩阵和三维场景的深度，并且对相机投影矩阵和三维场景的深度分别修正和插值，以此作为下一精细图像层迭代过程的初始值；

关于不同精度层之间三维场景的深度的插值，采用双线性插值、双三次插值方法或者其它类似插值方法实现；

关于不同精度层之间相机投影矩阵的修正，设相邻两级精度的图像在x和y方向的像素比为s₁和s₂，s₁,s₂＜1，在较低精度图像估计得到某个相机的相机投影矩阵为P^(k+1)，其中上标(k+1)代表图像金字塔结构的第k+1层，那么对应第k层图像的相机投影矩阵为

由粗到细金字塔方法的具体迭代算法如下：

即算法3的具体过程为：

输入：多幅图像；

输出：相机投影矩阵P_i(2≤i≤n_i)、三维场景的深度q和场景的三维表示；

1、计算图像的m层金字塔表示；

2、迭代：图像层k从第m层依次到第1层

(1)如果k≠m

以上一层估计的三维场景的深度q^(k+1)为基准，采用插值方法计算本层的三维场景的深度q^(k)，以其作为三维场景的深度的初始化；利用上一层图像估计的相机投影矩阵P_i ^{(k +1)}，2≤i≤n，根据式子(16)计算本层的相机投影矩阵P_i ^(k)，以其作为相机投影矩阵的初始化；

否则，在第m层图像

初始化：第i个相机的相机投影矩阵设置为P_i ^(m)＝[I₃ 0]，2≤i≤n，所有点的三维场景的深度信息设置为q^(m)＝1；

结束

(2)采用算法1或者算法2估计该层相机投影矩阵P_i ^(k)，2≤i≤n和三维场景的深度q^(k)；

结束迭代

3、输出相机投影矩阵和三维场景的深度：P_i←P_i ⁽¹⁾(2≤i≤n)，q←q⁽¹⁾。

4、根据三维场景的深度q，由式子(1)计算场景的三维表示。

更进一步的，实现相似三维重建或者欧几里德三维重建的具体过程为：

参数化具体为：相机投影矩阵由相机内部参数和相机外部参数描述：

P＝K[R t]

其中相机内部参数α_x、α_y、s、p_x和p_y包含在3×3矩阵相机外部参数由3×3的旋转矩阵R和3×1平移向量t确定，旋转矩阵R由三个角度参数确定，即分别绕x轴、y轴和z轴的旋转角度γ_x、γ_y和γ_z：

当相机内部参数和外部参数都是未知的，内部参数α_x、α_y、s、p_x、p_y、平移向量t、旋转角度γ_x、γ_y和γ_z和三维场景的深度q为待估计的未定参数，实现相似三维重建；

当相机内部参数是已知的，而外部参数是未知的，平移向量t、旋转角度γ_x、γ_y和γ_z和三维场景的深度q为待估计的未定参数，实现相似三维重建；

当相机的内部参数和外部参数都已知的条件下，实现欧几里德三维重建中，在这种情形中，三维场景的深度q为待估计的未定参数。

更进一步的，还能够推广到大基线情形，具体步骤为：在射影几何中，大基线情形是指相机之间的相对运动比较大，造成图像之间有显著的区别，在大基线情形中，具体来说，SFM三维重建分为三步：

第一步，从图像提取特征并匹配，提取的特征为：Harris特征、SIFT特征或KLT特征；

第二步，在所提取特征的基础上，估计特征点的三维信息和相机的相机投影矩阵；

第三步，在前面两步的基础上，利用算法3实现紧致SFM三维重建；

其中，以第二步估计得到的相机投影矩阵作为第三步的相机投影矩阵初始值，对第二步估计得到的三维场景的深度进行插值，作为第三步的三维场景的深度初始值。

与现有技术相比，本发明的有益效果为：本发明提出无特征提取的紧致SFM三维重建方案，与现有的SFM三维重建方法相比，能够一步完成紧致SFM三维重建。由于通过一步优化实现三维信息的估计，以目标函数值作为指标，能够得到最优解，至少是局部最优解，比已有方法有很大改进，且已初步得到实验验证。

附图说明

图1为本发明的三维重建实现流程图。

具体实施方式

下面对本发明做进一步的描述，但本发明的实施方式并不限于此。

如图1，S1.输入关于某场景的n幅图像，n≥2；

S2.建立与某个相机坐标系相一致的世界坐标系，设世界坐标系与第一相机的坐标系相一致，即世界坐标系的原点、x轴和y轴与第一相机的相机中心、第一相机成像平面的x轴和y轴重合，其z轴垂直指向第一相机的成像平面；

S3.以三维场景的深度和相机投影矩阵作为变量，所述三维场景的深度是指第1幅图像像素点对应的三维空间点具有的深度q；所述相机投影矩阵是指其它(n-1)幅图像的3×4矩阵P_i，2≤i≤n；

S4.构造类似光流估计的目标函数，所述目标函数是连续域上的变分目标函数或其离散形式的目标函数；

S5.采用由粗到细的金字塔方法，在连续域或者离散域上设计迭代算法对目标函数进行优化，输出表示场景三维信息的深度和代表相机相对位姿信息的相机投影矩阵；

S6.根据表示场景三维信息的深度，实现紧致的射影、相似或者欧几里德重建。

下面列举各种三维重建的实现方式

A.连续域中基于两幅灰度图像的射影三维重建

在本小节，详细介绍本发明的核心模型：在连续域上从两幅灰度图像实现射影(projective)三维重建。假设第一、二幅图像分别是I¹和I²(上标表示图像序号)，第一幅图像在位置(x,y)的灰度值为相应地，第二幅图像在位置(u,v) 的灰度值为虽然目前获取的图像大都是数字图像，即图像定义在离散格点(lattice)上，在本节提出的发明方案中假设图像是定义在连续域上，采用数值优化算法实现三维重建。

为了具体刻画三维信息，如下建立世界坐标系：其原点、x轴和y轴与第一个相机的相机中心(camera centre)、第一个相机成像平面的x轴和y轴重合，其z轴垂直指向第一个相机的成像平面。

根据射影几何原理，假设第一个相机的相机投影矩阵为[I₃ 0]∈R^3，4，其中I₃是一个3×3的单位阵，0是一个3×1的零向量；假设第二个相机的相机投影矩阵为一个3×4的矩阵P∈R^3，4。同时假设与第一幅图像像素(x,y)相对应的三维点具有深度q_x,y，即其z坐标为q_x,y；相应地，该点的三维坐标为

(q_x,y×x,q_x,y×y,q_x,y) (1)

这样，所拍摄场景或者目标的三维结构由定义在第一幅图像上的深度q_x,y描述；同时，第二个相机的相机投影矩阵P描述两个相机之间的相对运动信息。射影三维重建的目的即估计场景的深度信息q_x,y(在知道深度信息q_x,y的条件下，可以由式子(1)计算出相应的三维坐标)和相机投影矩阵P。

根据射影几何原理，与第一幅图像像素(x,y)相对应的三维点在第二幅图像中的成像位置(u_x,y,v_x,y)如下：

其中P_j为相机投影矩阵P的第j个行向量。为了表达式的简练，在不造成误解的情况下，省略下标x,y。从(2)可以看出，u和v是定义在图像域上、以相机投影矩阵P和深度q为参数的函数：u(P,q)和v(P,q)。

为了实现射影三维重建，构造如下变分目标函数：

其中为梯度算子，目标函数(3)中第一部分的依据是光流计算中的灰度不变假设(gray value constancy assumption)，即同一个三维空间点在不同图像具有相同的灰度值。如果只有灰度不变假设，该优化问题是一个病态问题。为此，在目标函数(3)中引入第二部分，假设两幅图像中的成像具有平滑的偏移量u-x和v-y，该平滑假设由目标函数(3)中的第二部分刻画。目标函数(3)的两部分分别叫做数据项和平滑项，与光流计算中数据项和平滑项相对应。

本发明采用迭代算法实现目标函数(3)的优化，其核心思想是变分法中的欧拉-拉格朗日方程(Euler-Lagrange equation)。为了符合变分法中欧拉-拉格朗日方程的标准形式，定义优化目标函数(3)中的积分量L如下(暂时忽略相机投影矩阵参数P，只考虑深度参数q)：

其中和

根据欧拉-拉格朗日方程，目标函数(3)取得极值必须满足：

具体地，欧拉-拉格朗日方程(5)为：

同时，根据导数在极值点等于0的性质可知，目标函数(3)满足

即

在迭代算法中，由前面迭代步骤已经得到P和q的当前估计，每次迭代过程的目的是求解P和q的增量δP和δq，并如下更新参数

为此，把式子(6)和(7)中的用代替，其中I_x和I_y是第二幅图像在位置(u,v)的偏导数；用q+δq的一阶差分近似q_x和q_y；用取代联立式子(6)和(7)的增量形式，可以在图像格点上建立n+12个以增量δP和δq为变量的线性方程，该线性方程组的变量数目也是n+12。其中n为图像像素的数目，即深度q有n个变量；相机投影矩阵有12个变量。

为了表达式的简练，把相机投影矩阵P的12个变量和深度q的n个变量构成一个维数为n+12的向量θ。通过整理可以把(6)和(7)的增量形式表示成如下线性方程组

Hδθ+b＝0 (9)

这样，每次迭代归结为求解

δθ＝-H^-1b (10)

综合前面的内容，给定两幅灰度图像，以及相机投影矩阵和深度的初始值，可以采用如下算法实现SFM射影三维重建。

即算法1的具体过程为：

输入：两幅灰度图像I¹和I²，以及相机投影矩阵P和深度q的初始化

输出：相机投影矩阵P和深度q

迭代

1、由式子(6)和式子(7)确定式子(9)中的H和b；

2、由式子(10)计算增量δθ，并确定相应的增量δP和δq；

3、由式子(8)更新参数P和q；

直到收敛

为了克服局部极值点的困难，采用从粗到细(from coarse to fine)的金字塔方法。即首先在较低分辨率的图像实现三维重建；然后对所估计深度进行插值，并且修正相机投影矩阵，以此作为下一层较高分辨率三维重建的初始解；直到最高分辨率。

关于深度的插值，可以采用双线性插值(bilinear interpolation)、双三次插值(bicubic interpolation)或者其它类似插值方法实现。

关于不同精度层之间相机投影矩阵的修正。假设相邻两层的图像在x和y方向的像素比为s₁和s₂(s₁,s₂＜1)，在较低精度层估计的第二个相机投影矩阵为P⁽ⁱ⁺¹⁾(其中上标i+1代表图像金字塔结构的第i+1层)，那么对应第i层图像的第二个相机投影矩阵为

即算法2的具体过程为：

输入：两幅灰度图像

输出：相机投影矩阵P和深度q

1、计算两幅图像的m层金字塔表示；

2、for i＝m：-1：1

if i≠m

以上一层估计的深度q⁽ⁱ⁺¹⁾为基准，采用插值方法计算在当前层的深

度q⁽ⁱ⁾，以其作为深度的初始化；

利用上一层估计的相机投影矩阵P⁽ⁱ⁺¹⁾，根据式子(11)计算当前层

的相机投影矩阵P⁽ⁱ⁾，以其作为第二个相机投影矩阵的初始化；

else

初始化：第二个相机的相机投影矩阵设置为P^(m)＝[I₃ 0]和所有点

的深度信息为q^(m)＝1；

end

采用算法1估计该层的相机投影矩阵P⁽ⁱ⁾和深度q⁽ⁱ⁾；

end

B.基于两幅灰度图像的射影三维重建的离散化形式

由于数字图像本身以离散的形式存在，下面对目标函数(3)直接给出其离散形式：

目标函数(12)本质上是一个非线性最小二乘问题，可以采用高斯-牛顿(Gauss-Newton)算法(或者其它类似算法，如Levenberg Marquadt，LM算法)实现迭代优化。为了简单起见，首先考虑以u和v作为变量的情形。对于目标函数(12)的第一部分数据项，采用高斯-牛顿近似方法，可以得到相关的高斯-牛顿海森(Hessian)矩阵H_data和梯度向量b_data：

对于目标函数(12)的第二部分平滑项，采用一阶差分代替其中的偏导数。具体来说，和同样地，和这里可以采用其它形式的差分以近似偏导数，如这些细微的变化不会改变算法的整体结构。

由于差分计算涉及到其邻居像素，对平滑项的推导不如数据项(13)简洁，需要对单个像素点进行分析。跟迭代算法1类似，由前一次迭代已经得到u和v的估计，每次迭代的目的是确定其增量δu和δv。这样，目标函数(12)中的平滑项可以表示成：

上述两项之和可以用高斯-牛顿海森矩阵H_x,y和梯度向量b_x,y表示：

其中δ_x,y＝[δu_x-1,y δu_x,y δu_x,y-1 δv_x-1,y δv_x,y δv_x,y-1]^T，

对所有像素建立类似(14)的高斯-牛顿近似，联立得到平滑项的高斯-牛顿海森矩阵H_smooth和梯度向量b_smooth

注意，式子(15)中求和项里的加法与通常意义上的加法有所不同，因为每个δ_x,y涉及到其邻居像素，如δ_x,y的定义可以看出，因此式子(15)中的求和运算需要考虑δ_x,y之间的对齐。

综合式子(13)和式子(15)，可以求得目标函数(12)的高斯-牛顿海森矩阵H和梯度向量b：

增量δu和δv可以通过求解如下方程得到：

-H^-1b (17)

在本发明提出的方案中，最终的参数是相机投影矩阵P和深度q。跟连续域的推导一样，为了简练的表达式，建立一个包括相机投影矩阵P和深度q的参数向量θ。由复合函数法则，增量δθ为：

δθ＝-(J^THJ)^-1J^Tb (18)

其中雅可比矩阵

即算法3的具体实现过程为：

输入：两幅灰度图像I¹和I²，以及相机投影矩阵P和深度q的初始化

输出：相机投影矩阵P和深度q

迭代

1、计算式子(18)中的H、J和b；

2、由(18)计算增量δθ，并分别给出相应的增量δP和δq；

3、由式子(8)更新参数P和q；

直到收敛

同样，与连续域的三维重建一样，可以采用由粗到精的金字塔方法实现离散形式的三维重建，其基本框架与算法2一样。不同之处在于在离散三维重建中，采用算法2实现各层的三维重建。为了避免冗余，由粗到细的金字塔三维重建离散算法这里予以省略。

C.基于两幅彩色图像的三维重建

对于基于彩色图像的三维重建，其原理和基于灰度图像的三维重建一样。彩色图像的表示方案有很多种，比如RGB(Red Green Blue)、HSV(Hue Saturation Value)、HSL(Hue Saturation Lightness)和HSI(Hue Saturation Intensity)。下面以在RGB彩色图像的离散三维重建为例，说明如何处理基于彩色图像的三维重建。构造离散目标函数如下(可以类似地构造连续域上的目标函数)：

目标函数(19)的优化算法与目标函数(12)完全相同，这里予以省略。另外，也同样可以实现连续域上彩色图像的三维重建，其实现算法与连续域灰度图像的三维重建类似。

D.基于多幅图像的三维重建

基于n(n>2)幅图像的三维重建，其基本算法与基于两幅图像的三维重建一样，构造一个类似(3)或者(12)的目标函数，包括数据项和平滑项。与基于两幅图像的三维重建一样，可以把世界坐标系设置为第一个相机的坐标系，这样第一个相机的投影矩阵为[I₃ 0]∈R^3，4，其它n-1相机投影矩阵和深度q为待估计参数。

基于多幅图像的三维重建，数据项和平滑项的构造有很多不同的选择。这里分别对数据项和平滑项给出两种不同的构造方案，但不限于此。数据项的第一种构造方案如下

其中和上标表示图像序列，(uⁱ,vⁱ)代表第i幅图像上与(x,y)相对应的坐标。第二种方案是

在第一种方案，式(20)中，与两幅图像(12)的数据项非常相似，即灰度不变假设，假设同一个点在所有图像中具有相同的灰度值。在第二种方案(21)中，灰度不变假设略有改变，只考虑相邻图像之间的灰度不变假设。第二种方案比较适合基于视频序列的三维重建。

在平滑项中，给出两种构造方案(22)和(23)

在(22)中，偏移量都是以第一幅图像为基准，而在(23)中考虑相邻图像之间的偏移量。

基于多幅图像的三维重建实现算法与两幅图像的三维重建类似，其具体算法予以省略。

E.在三维重建中引入梯度(gradient)不变假设

为了克服不同光照条件带来图像的变化，在优化目标函数中引入梯度(gradient)不变假设，即假设同一个三维点在不同图像中的梯度保持不变。以基于两幅灰度图像的三维重建为例，说明如何在三维重建中引入梯度不变假设。构造目标函数如下：

对于式子(24)的迭代优化算法与算法3一致，这里予以忽略。

F.在优化目标函数中引入鲁棒函数

由于优化目标函数(3)和(12)，以及在其它类似的构造中，数据项和平滑项都是以平方项的形式出现，这样的模型在场景或者目标的深度出现不平滑的情况时，效果比较差。为此，在优化目标函数中引入某个鲁棒函数ρ。以(12)为例，

鲁棒函数ρ的选取有很多，比如Charbonnier函数其中∈是一个很小的正常数；Lorentzian函数σ为某个常数。

G.在深度上引入平滑约束

除了可以在像素偏移量上引入平滑约束之外，另外一个方案是直接把平滑约束引入到三维目标的深度q上。以离散情形(12)为例，可以构造如下目标函数：

其中

另外一种方案是同时对像素偏移量和深度引入平滑约束，优化目标函数变为：

H.基于两幅灰度图像的相似三维重建和欧几里德三维重建

在上面的发明方案中，重建的结果是射影三维结构，不是常见的欧几里德(Euclidean)三维结构，这节提出相似(similarity)三维重建和欧几里德三维重建方案。

根据射影几何原理，相机的投影矩阵可以由相机内部参数和相机外部参数描述：

P＝K[R t] (26)

其中相机内部参数α_x、α_y、s、p_x和p_y包含在3×3矩阵相机外部参数由旋转矩阵R和平移向量t确定，其中R刻画世界坐标系统到相机坐标系统的旋转变换。

与前面射影重建相似，假设世界坐标系统与第一个相机的坐标系统相同，这样相机之间的运动信息完全由第二个相机的外部参数描述。为了简单起见，假设两个相机具有相同的内部参数。为了把第一相机投影矩阵化为标准型[I₃ 0]，相应地设第二相机形式为

P＝K[R t]K^-1 (27)

关于描述场景或者目标的深度参数q与射影三维重建的一样。

实现相似三维重建的方案与射影重建相似，即优化目标函数(3)或者(12)，不同之处在于第二相机投影矩阵P的形式(27)。这里，只给出相似三维重建的离散实现方案。与射影三维重建的连续方案一样，可以实现连续域的相似三维重建。

旋转矩阵R由三个角度参数确定，即分别绕x轴、y轴和z轴的旋转角度γ_x、γ_y和γ_z：

假设在射影三维重建中n+12维参数向量θ的前12个量为相机投影矩阵参数，其余的为n个深度参数。同时，假设在相似三维重建中，参数向量θ′的前11个量为3个旋转角度参数γ＝[γ_x,γ_y,γ_z]^T、3个平移向量参数t和5个相机内部参数α_x、α_y、s、p_x和p_y，构成一个11维的向量θ″，后面有n个深度参数。关于式子(27)，定义雅可比矩阵

其中I_n为n×n的单位阵。那么，在相似三维重建中，迭代过程求解参数增量为：

δθ′＝-(J″^TJ^THJJ″)^-1J″^TJ^Tb (28)

其中H、J和b是式子(18)的H、J和b。

当已知相机的部分参数的条件下，比如当相机是标定时，即已知相机内部参数，同样可以实现上述相似三维重建。当相机的内部参数和外部参数都已知的条件下，可以实现欧几里德三维重建。

I.在大基线(large baseline)情形下的SFM三维重建

在射影几何中，大基线(large baseline)是指相机之间的相对运动比较大，造成图像之间有显著的区别，究其原因，可能是因为相机之间的旋转角度或者平移过大，也可能是因为相机之间的焦距差别太大。在大基线情形中，可以综合现有的三维重建方法和本发明前面提出的方案，实现大基线的SFM三维重建。具体来说，SFM三维重建分为三步，第一步从图像提取特征并匹配，如Harris特征、SIFT特征或KLT特征；第二步，在所提取特征的基础上，估计特征点的三维信息和相机的相机投影矩阵；第三步，在前面两步的基础上，利用前面提出的方法实现紧致SFM三维重建。其中，以第二步估计得到的相机投影矩阵作为第三步的相机投影矩阵初始值，对第二步估计得到的三维场景的深度进行插值，作为第三步的三维场景的深度初始值。

专业术语的中英文对照

基于运动的建模 structure from motion(SFM)

紧致的 dense

相机投影矩阵 camera projection matrix

射影的 projective

相似 similarity

欧几里德 Euclidean

格点 lattice

变分 variational

欧拉-拉格朗日方程 Euler-Lagrange equation

从粗到细 from coarse to fine

金字塔 pyramid

高斯-牛顿算法 Gauss-Newton algorithm

高斯-牛顿海森矩阵 Gauss-Newton Hessian matrix

海森矩阵 Hessian matrix

内部参数 internal parameters

外部参数 external parameters

大基线 large baseline

以上所述的本发明的实施方式，并不构成对本发明保护范围的限定。任何在本发明的精神原则之内所作出的修改、等同替换和改进等，均应包含在本发明的权利要求保护范围之内。

Claims (8)

1.一种无特征提取的紧致SFM三维重建方法，其特征在于，包括以下步骤：

S1.输入关于某场景的n幅图像，n≥2；

S2.建立与某个相机坐标系相一致的世界坐标系，设世界坐标系与第一相机的坐标系相一致，即世界坐标系的原点、x轴和y轴与第一相机的相机中心、第一相机成像平面的x轴和y轴重合，其z轴垂直指向第一相机的成像平面；

S3.以三维场景的深度和相机投影矩阵作为变量，所述三维场景的深度是指第1幅图像像素点对应的三维空间点具有的深度q；所述相机投影矩阵是指其它(n-1)幅图像的3×4矩阵P_i，2≤i≤n；

S4.构造类似光流估计的目标函数，所述目标函数是连续域上的变分目标函数或其离散形式的目标函数；

S5.采用由粗到细的金字塔方法，在连续域或者离散域上设计迭代算法对目标函数进行优化，输出表示场景三维信息的深度和代表相机相对位姿信息的相机投影矩阵；

S6.根据表示场景三维信息的深度，实现紧致的射影、相似或者欧几里德重建。

2.根据权利要求1所述的无特征提取的紧致SFM三维重建方法，其特征在于，在实现射影三维重建中，参数化具体为：在建立世界坐标系的同时，其第一相机的相机投影矩阵为[I₃ 0]∈R^3，4，其中I₃是一个3×3的单位阵，0是一个3×1的零向量；其它相机投影矩阵P_i∈R^3，4，2≤i≤n，作为待估计的未知参数；场景的三维结构由定义在第一副图像上的三维场景的深度决定：假设与第一幅图像像素(x,y)相对应的三维空间点的三维场景的深度为q_x,y，则该三维空间点的三维坐标为

(q_x,y×x,q_x,y×y,q_x,y) (1)

在射影三维重建中，相机投影矩阵P_i和三维场景的深度q_x,y作为待估计的未定参数，为了表达式的简练，在不造成误解的情况下，省略下标x，y。

3.根据权利要求2所述的无特征提取的紧致SFM三维重建方法，其特征在于，实现连续域上射影三维重建的具体实现过程为：

构造的连续域上的目标函数具体为：

f(P₂,…,P_n,q)＝f_data+f_{smooth_uv}+f_{smooth_depth} (2)

其中

$f_{data}=\∫\∫dxdy\mathrm{\ρ}\left({\mathrm{\Σ}}_{i=2}^n{\mathrm{\Σ}}_{I=C_1}^{C_k}\right({\left(I_{x,y}^1-I_{u^i,v^i}^i\right)}^2+\mathrm{\β}{\left(\▿I_{x,y}^1-\▿I_{u^i,v^i}^i\right)}^2\left)\right)---\left(3\right)$

$f_{\mathrm{smooth}\_\mathrm{uv}}=\mathrm{\α}\∫\∫\mathrm{dxdy\ρ}\left({\mathrm{\Σ}}_{i=2}^n({|\▿(u^i-x)|}^2+{|\▿(u^i-y)|}^2)\right)---\left(4\right)$

$f_{smooth\_depth}={\mathrm{\τ}}_1\∫\∫dxdy\mathrm{\ρ}\left(\right|\▿q{|}^2+{\mathrm{\τ}}_2\left|\mathrm{\Δ}q{|}^2\right)---\left(5\right)$

对上述目标函数的说明如下：(a)为梯度算子， 为拉普拉斯算子；(b)目标函数分为三部分，数据项f_data，偏移平滑项f_{smooth_uv}和深度平滑项f_{smooth_depth}，其中α、β、τ₁和τ₂是非负权重；(c)图像有k个色彩分量C₁,…,C_k，代表第一幅图像在位置(x,y)的色彩I分量值，相应地，为第i幅图像在位置(uⁱ,vⁱ)的色彩I分量值；(d)鲁棒函数ρ的引入是为了克服深度发生剧变带来的影响，鲁棒函数ρ为Charbonnier函数其中∈是一个足够小的的正数，∈＜10^-6；或者为Lorentzian函数 σ为某个常数；当不引入鲁棒函数，则ρ(x)＝x；(e)uⁱ和vⁱ是定义在图像域上、以相机投影矩阵P_i和深度q为参数的函数：和代表与第一幅图像像素(x,y)相对应的三维空间点在第i幅图像的成像位置

$\left\{\begin{array}{c}{u}_{x,y}^i=\frac{P_{i,1}{\[q_{x,y}\×x,q_{x,y}\×y,q_{x,y},1\]}^T}{P_{i,3}{\[q_{x,y}\×x,q_{x,y}\×y,q_{x,y},1\]}^T}\\ v_{x,y}^i=\frac{P_{i,2}{\[q_{x,y}\×x,q_{x,y}\×y,q_{x,y},1\]}^T}{P_{i,3}{\[q_{x,y}\×x,q_{x,y}\×y,q_{x,y},1\]}^T}\end{array}\right.,$

其中P_i,j为第i个相机投影矩阵P_i的第j个行向量；为了表达式的简练，在不造成误解的情况下，在和中省略下标x,y；

在连续域上设计的迭代优化算法具体为：因为三维场景的深度是定义在第一幅图像上的连续函数，在极值点必须满足欧拉-拉格朗日方程；同时，在极值点对相机投影矩阵的偏导数为0；在图像的离散格点上，联合欧拉-拉格朗日方程和对相机投影矩阵偏导数为0的两类方程，并采用增量方式表示形式，能够把求解相机投影矩阵和三维场景深度增量的迭代过程转化为求解如下线性方程组

Hδθ+b＝0 (6)

其中向量θ由相机投影矩阵P_i，2≤i≤n和三维场景的深度q按次序构造而成；这样，每次迭代归结为求解

δθ＝-H^-1b (7)，

从而确定相应的增量δP_i和δq；根据所求解的增量更新相机投影矩阵P_i和q，P_i←δP_i+P_i，q←δq+q，直到收敛；

即定义算法1的具体过程为：

输入：n幅图像，初始化三维场景的深度q和相机投影矩阵P_i，2≤i≤n；

输出：相机投影矩阵P_i(2≤i≤n_i)、三维场景的深度q和场景的三维表示；

1、迭代

1)、由欧拉-拉格朗日方程和目标函数对相机投影矩阵的偏导数为0确定式子(7)中的H和b；

2)、由式子(7)计算增量δθ，并确定相应的增量δP_i和δq；

3)、更新参数P_i，2≤i≤n和q：P_i←δP_i+P_i，q←δq+q；

直到收敛

2、根据收敛后的三维场景的深度q，由式子(1)计算场景的三维表示。

4.根据权利要求3所述的无特征提取的紧致SFM三维重建方法，其特征在于，目标函数(2)中的数据项和偏移平滑项能够采用其它类似的变化形式：

$f_{data}=\∫\∫dxdy\mathrm{\ρ}\left({\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{I=C_1}^{C_k}\right({\left({\stackrel{\‾}{I}}_{x,y}-I_{u^i,v^i}^i\right)}^2+\mathrm{\β}{\left(\▿{\stackrel{\‾}{I}}_{x,y}-\▿I_{u^i,v^i}^i\right)}^2\left)\right)---\left(8\right)$

$f_{smooth\_uv}=\mathrm{\α}\∫\∫dxdy\mathrm{\ρ}\left({\mathrm{\Σ}}_{i=2}^n\right(|\▿(u^i-u^{i-1}){|}^2+|\▿(v^i-v^{i-1}\left){|}^2\right))---\left(9\right)$

其中和

鲁棒函数的引入还能够以其它变化形式出现，式子(3)的数据项的另外一种变化形式为：

$f_{data}=\∫\∫dxdy\{\mathrm{\ρ}({\mathrm{\Σ}}_{i=2}^n{\mathrm{\Σ}}_{I=C_1}^{C_k}{\left(I_{x,y}^1-I_{u^i,v^i}^i\right)}^2+\mathrm{\β}\mathrm{\ρ}({\Σ}_{i=2}^n{\Σ}_{I=C_1}^{C_k}{\left(\▿I_{x,y}^1-\▿I_{u^i,v^i}^i\right)}^2\left)\right\}---\left(10\right).$

5.根据权利要求2所述的无特征提取的紧致SFM三维重建方法，其特征在于，构造离散形式的目标函数具体为：

f(P₂,…,P_n,q)＝f_data+f_{smooth_uv}+f_{smooth_depth} (11)

其中

$f_{data}={\mathrm{\Σ}}_{x,y}\mathrm{\ρ}\left({\mathrm{\Σ}}_{i=2}^n{\mathrm{\Σ}}_{I=C_1}^{C_k}\right({\left(I_{x,y}^1-I_{u^i,v^i}^i\right)}^2+\mathrm{\β}{\left(\▿I_{x,y}^1-\▿I_{u^i,v^i}^i\right)}^2\left)\right)---\left(12\right)$

$f_{smooth\_uv}={\mathrm{\α\Σ}}_{x,y}\mathrm{\ρ}\left({\mathrm{\Σ}}_{i=2}^n\right(|\▿(u^i-x){|}^2+|\▿(v^i-y){|}^2))---\left(13\right)$

$f_{smooth\_depth}=x_1{\mathrm{\Σ}}_{x,y}\mathrm{\ρ}\left(\right|\▿q{|}^2+{\mathrm{\τ}}_2\left|\mathrm{\Δ}q{|}^2\right)---\left(14\right)$

离散目标函数(11)及其变化形式的迭代优化算法具体为：离散形式的目标函数(11)在本质上是一个非线性最小二乘问题，能够采用常规的Levenberg-Marquadt算法或高斯-牛顿算法，每次迭代过程归结为求解一个线性方程组(15)：

δθ＝-(H+μI)^-1b (15)，

其中H是海森矩阵或者高斯-牛顿海森矩阵，b是梯度向量，μ是非负数，取决于采用Levenberg-Marquadt算法或高斯-牛顿算法，从而确定相应的增量δP_i和δq；根据增量更新参数P_i和q，P_i←δP_i+P_i，q←δq+q，直到收敛；

定义算法2具体实现过程：

输入：n幅图像，以及相机投影矩阵P_i和三维场景的深度q的初始化，2≤i≤n；输出：相机投影矩阵P_i(2≤i≤n_i)、三维场景的深度q和场景的三维表示；

1、迭代

1)、计算式子(15)中的高斯-牛顿海森矩阵H和梯度向量b；

2)、由式子(15)计算增量δθ，并分别确定相应的增量δP_i和δq；

3)、更新参数P_i和q：P_i←δP_i+P_i，q←δq+q，2≤i≤n；

直到收敛；

2、根据收敛后的三维场景的深度q，由式子(1)计算场景的三维表示。

6.根据权利要求3或5所述的无特征提取的紧致SFM三维重建方法，其特征在于，由粗到细的金字塔方法步骤具体为：计算图像的n层金字塔表示；在最粗图像层，初始化n-1个相机投影矩阵为P_i＝[I₃ 0]，2≤i≤n，所有点的三维场景的深度初始化为1；由粗到细依次估计相机投影矩阵和三维场景的深度，并且对相机投影矩阵和三维场景的深度分别修正和插值，以此作为下一精细图像层迭代过程的初始值；

关于不同精度层之间三维场景的深度的插值，采用双线性插值或双三次插值方法实现；

关于不同精度层之间相机投影矩阵的修正，设相邻两级精度的图像在x和y方向的像素比为s₁和s₂，s₁,s₂＜1，在较低精度图像层估计得到某个相机的相机投影矩阵为P^(k+1)，其中上标(k+1)代表图像金字塔结构的第k+1层，那么对应第k层图像的相机投影矩阵为

$P^{\left(k\right)}=\left[\begin{array}{ccc}{s}_1& 0& 0\\ 0& s_2& 0\\ 0& 0& 1\end{array}\right]P^{(k+1)}\left[\begin{array}{cccc}1/s_1& 0& 0& 0\\ 0& 1/s_2& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]---\left(16\right)$

由粗到细金字塔方法的具体迭代算法如下：

即定义算法3的具体过程为：

输入：n幅图像；

输出：相机投影矩阵P_i(2≤i≤n_i)、三维场景的深度q和场景的三维表示；

1、计算图像的m层金字塔表示；

2、迭代：图像层k从第m层依次到第1层

(1)如果k≠m

以上一层估计的三维场景的深度q^(k+1)为基准，采用插值方法计算在本层的三维场景的深度q^(k)，以其作为三维场景的深度的初始化；利用上一层图像估计的相机投影矩阵P_i ^(k+1)，2≤i≤n，根据式子(16)计算本层的相机投影矩阵P_i ^(k)，以其作为相机投影矩阵的初始化；

否则，在第m层图像

初始化：第i个相机的相机投影矩阵设置为P_i ^(m)＝[I₃ 0]，2≤i≤n，所有点的三维场景的深度信息设置为q^(m)＝1；

结束

(2)采用算法1或者算法2估计该层相机投影矩阵P_i ^(k)，2≤i≤n和三维场景的深度q^(k)；

结束迭代

3、输出相机投影矩阵和三维场景的深度：P_i←P_i ⁽¹⁾(2≤i≤n)，q←q⁽¹⁾；

4、根据三维场景的深度q，由式子(1)计算场景的三维表示。

7.根据权利要求6所述的无特征提取的紧致SFM三维重建方法，其特征在于，实现相似三维重建或者欧几里德三维重建的具体过程为：

参数化具体为：相机投影矩阵由相机内部参数和相机外部参数描述：

P＝K[R t]

其中相机内部参数α_x、α_y、s、p_x和p_y包含在3×3矩阵相机外部参数由3×3的旋转矩阵R和3×1平移向量t确定，旋转矩阵R由三个角度参数确定，即分别绕x轴、y轴和z轴的旋转角度γ_x、γ_y和γ_z：

$R=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left({\mathrm{\γ}}_x\right)& -\sin \left({\mathrm{\γ}}_x\right)\\ 0& \sin \left({\mathrm{\γ}}_x\right)& \cos \left({\mathrm{\γ}}_x\right)\end{array}\right]\left[\begin{array}{ccc}\cos \left({\mathrm{\γ}}_y\right)& 0& \sin \left({\mathrm{\γ}}_y\right)\\ 0& 1& 0\\ -\sin \left({\mathrm{\γ}}_y\right)& 0& \cos \left({\mathrm{\γ}}_y\right)\end{array}\right]\left[\begin{array}{ccc}\cos \left({\mathrm{\γ}}_z\right)& -\sin \left({\mathrm{\γ}}_z\right)& 0\\ \sin \left({\mathrm{\γ}}_z\right)& \cos \left({\mathrm{\γ}}_z\right)& 0\\ 0& 0& 1\end{array}\right]$

当相机内部参数和外部参数都是未知的，内部参数α_x、α_y、s、p_x、p_y、平移向量t、旋转角度γ_x、γ_y和γ_z和三维场景的深度q为待估计的未定参数，实现相似三维重建；

当相机内部参数是已知的，而外部参数是未知的，平移向量t、旋转角度γ_x、γ_y和γ_z和三维场景的深度q为待估计的未定参数，实现相似三维重建；

当相机的内部参数和外部参数都已知的条件下，实现欧几里德三维重建中，在这种情形中，三维场景的深度q为待估计的未定参数。

8.根据权利要求6所述的无特征提取的紧致SFM三维重建方法，其特征在于，所述无特征提取的紧致SFM三维重建方法还能够推广到大基线情形，具体步骤为：在射影几何中，大基线情形是指相机之间的相对运动比较大，造成图像之间有显著的区别，在大基线情形中，具体来说，SFM三维重建分为三步：

第一步，从图像提取特征并匹配，提取的特征为：Harris特征、SIFT特征或KLT特征；

第二步，在所提取特征的基础上，估计特征点的三维信息和相机投影矩阵；

第三步，在前面两步的基础上，利用算法3实现紧致SFM三维重建；

其中，以第二步估计得到的相机投影矩阵作为第三步的相机投影矩阵初始值，对第二步估计得到的三维场景的深度进行插值，作为第三步的三维场景的深度初始值。