CN101751697A - 一种基于统计模型的三维场景重建方法 - Google Patents

Abstract

本发明涉及一种基于统计模型的三维场景重建方法，技术特征在于步骤为：使用Harris角点检测算法提取每幅图像中的角点，然后生成三维点集合X与摄像机参数集合M；利用马尔可夫链蒙特卡罗MCMC方法估计图像角点与三维点之间的匹配概率，利用图像角点与三维点之间的匹配概率对图像角点进行加权平均得到虚拟测量点矩阵V；采用可处理遮挡的射影因式分解算法对虚拟测量点进行射影重构，加入确定性退火算法迭代求解全局最优射影重构结果，利用基于绝对对偶二次曲面的摄像机自标定算法将射影重构提升为度量重构。将原本一次性计算的过程转化为迭代求精的过程。即使匹配关系未知或初始匹配结果较差，通过迭代反馈的方法，最终仍可精确得到三维重建结果。

Description

一种基于统计模型的三维场景重建方法

技术领域

本发明涉及一种基于统计模型的三维场景重建方法，属于计算机视觉中基于图像的三维场景重建领域，特别是涉及图像中特征点之间匹配关系的建立与优化问题。

背景技术

目前在三维重建领域内，Noah Snavely，Steven M.Seiz和Richard Szeliski在文献Modeling the World from Internet Photo Collections(IJCV，vol.80，pp.189-210，2008)中提出了一种利用互联网上大规模图片进行三维重建的完整算法。该算法对每幅图像提取局部不变特征，通过kd树在两两图像间建立匹配并采用归一化8点算法计算基本矩阵，然后利用特征点跟踪得到所有图像间的特征匹配关系，摄像机内参数从图像EXIF标签读取，最后利用bundle adjustment(BA)算法迭代恢复出三维点结构与摄像机参数。但是如果图像之间视角变化较大，则上述算法匹配结果的精度无法满足要求，产生的大量错误匹配会直接影响三维重建的结果。Frank Dellaert在文献Monte Carlo EM forData-Association and its Applications in Computer Vision(doctoral dissertation，tech.report CMU-CS-01-153，Carnegie Mellon University，September，2001)中引入蒙特卡罗EM算法解决计算机视觉领域中的匹配问题，通过概率方法建立起图像特征点与三维点之间的软匹配并进行迭代优化。但是其马尔可夫链蒙特卡罗(MCMC)算法采样效率不高，且仅使用局部优化算法BA进行模型求解，致使迭代容易陷入局部极值。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于统计模型的三维场景重建方法，可以在图像间匹配关系未知或不准确的条件下有效恢复精确匹配关系与三维场景模型。

本发明的思想在于：

技术方案

一种基于统计模型的三维场景重建方法，其特征在于步骤如下：

步骤1使用Harris角点检测算法提取每幅图像中的角点：针对无序数字图像集合，采用Harris角点检测算法提取每幅图像的角点作为图像特征点，以角点集合U＝{u_ik|k∈1…n_j，i∈1…m}表示，其中k为图像中的角点序号，i为图像序号，n_i为第i幅图像中的角点数量，m为图像数量，u_ik为第i幅图像中第k个角点的二维坐标；

步骤2生成三维点集合X与摄像机参数集合M：采用满足高斯分布N(0，1)的随机函数生成N个三维点，得到三维点集合X＝{x_i|i∈1…n}，x_i为第i个三维点的坐标，N的值等于所有图像中最大的角点数量；

所述的摄像机参数集合M，M＝{m_i|i∈1…m}中每幅图像的摄像机参数初始值m_i取为

步骤3利用马尔可夫链蒙特卡罗MCMC方法估计图像角点与三维点之间的匹配概率：利用基本投影公式计算得到三维点在所有图像的重投影点集合H＝{h_ij|j∈1…n，i∈1…m}，h_ij为第j个三维点在i幅图像上的投影点；计算图像角点u_ik和三维点x_j之间的权重：

$w(u_{\mathrm{ik}},x_j)=\log \left(\sqrt{2\mathrm{\π}}\mathrm{\σ}\right)+\frac{{\left|\right|u_{\mathrm{ik}}-h_{\mathrm{ij}}\left|\right|}^2}{{\mathrm{\σ}}^2}$

其中：σ为干扰图像的高斯噪声的方差，初始取为σ⁰；

然后利用MCMC统计模拟方法对图像角点u_ij与三维点x_j之间的匹配概率f_ijk进行概率估计：对于任意一幅图像i的图像角点与三维点匹配向量j_i＝{j_ik|k∈1…n_i}，匹配向量满足约束：

图像i的采样过程如下：

步骤a：为图像i中每个角点随机选择一个三维点，建立角点与三维点间的匹配边，得到匹配向量j_i ⁰，约束条件为：不同角点不可对应同一三维点；

步骤b：从图像i的角点集合u_i中随机选择一点u_ik，计算u_ik与任意三维点之间的转移概率：

$q^{j_i}(u_{\mathrm{ik}},x_j)=\left\{\begin{array}{ccc}\frac{\exp (-w(u_{\mathrm{ik}},x_j))}{{\mathrm{\Σ}}_{j\≠j_{\mathrm{ik}}}\exp (-w(u_{\mathrm{ik}},x_j))}& \mathrm{if}& j{\≠j}_{\mathrm{ik}}\\ 0& \mathrm{if}& j=j_{\mathrm{ik}}\end{array}\right.$

以概率

从三维点集合X中选择一点x_j，以点x_j与u_ik建立匹配边；

步骤c：选择x_j的原有匹配边(x_j，u_ik′)对应的角点u_ik′作为下一次循环的图像角点；

步骤d：重复步骤b直至某次迭代中选择的图像角点在之前已被选择，得到由图像角点和三维点形成的闭环；

步骤e：删除在建立闭环过程中不属于闭环中的临时边，得到新的闭环C；

步骤f：删除闭环C中在步骤a的原有匹配关系，得到图像特征点与三维点之间新的匹配j′_i；

重复3000～5000次步骤a～f，保存每次采样生成的匹配向量j_i ^r，并计算匹配概率：

其中δ为克罗内克函数，若j_ik＝j，δ(j_ik，j)＝1；否则δ(j_ik，j)＝0；

步骤4：利用图像角点与三维点之间的匹配概率对图像角点进行加权平均得到虚拟测量点矩阵V

其中：v_ij为V中第i行第j个元素；

步骤5采用可处理遮挡的射影因式分解算法对虚拟测量点进行射影重构：具体步骤如下：

步骤a：利用归一化8点算法计算任意两幅图像间的基本矩阵F_ij与极点e_ij；

步骤b：计算每个虚拟测量点v_ip的射影深度信息λ_ip

在得到所有射影深度信息后，利用v_ip＝λ_ipv_ip计算新的虚拟测量点；

步骤c：对虚拟测量矩阵V中的缺失部分进行线性拟合：首先对虚拟测量矩阵V进行T次采样，每次随机抽取4列，得到一个列四元组；对每个列四元组求取其生成子空间，得到B_t，t∈1…T；对B_t进行SVD分解：B_t＝S_tV_tD_t ^T，B_t ^⊥取为S_t的最后4列，其中⊥表示线性空间的补集；将每个B_t ^⊥按列组合，对得到的联合矩阵使用SVD分解：

B取为S的最后四列；将V中存在缺失元素的列表示为B的列的线性组合，即V_j＝k₁b₁+k₂b₂+k₃b₃+k₄b₄，利用列中已知元素计算相应系数，进而计算列中未知元素，最终得到完全矩阵V^*；

步骤d：对V^*进行SVD分解求取三维点坐标与摄像机参数：对V^*进行SVD分解，V^*＝USV^T，取S的前4列得到矩阵摄像机参数三维点坐标X＝V^T；利用BA算法对三维点坐标和摄像机参数进行优化，以得到的优化结果替换掉原有三维点坐标与摄像机参数；

步骤6加入确定性退火算法迭代求解全局最优射影重构结果：重复步骤3～5进行循环迭代，σ在迭代过程中以等比例策略，将σ参数取为原σ的0.85～0.95，在每次迭代中，比较当前σ与给定阈值σ_min大小，若σ≤σ_min则迭代结束，进行后续步骤。否则置σ^t+1＝λσ^t，进入下一次迭代过程；

步骤7利用基于绝对对偶二次曲面的摄像机自标定算法将射影重构提升为度量重构：首先计算绝对对偶二次曲面Q_∞ ^*，然后将Q_∞ ^*分解为

得到单应矩阵

$H=\left[\begin{array}{cc}K& 0\\ -p^{T}K& 1\end{array}\right]$

其中：p为无穷远平面的前三个坐标，将射影重构变换到度量重构，得到欧式空间中三维场景重建的三维结构M×H与摄像机参数H^-1X。

有益效果

本发明提出的基于统计模型的三维场景重建方法的优越性在于：

目前的三维重建算法在初始步骤中就对各幅图像的特征点建立确定的匹配关系，并且其在后续步骤中不再变化。当图像之间的摄像机视角变化很小时，如视频序列，这种方法可以取得较好结果。然而对于视角变化较大(30°以上)的图像，匹配结果误差较大，存在大量错误匹配。而特征点匹配又是三维场景重建中非常关键的一步，匹配质量对重建结果有着非常大的影响，如果存在较多错误匹配，将无法得到正确的三维重建结果。

基于统计模型的三维场景重建技术不要求初始就建立高精度特征匹配关系，甚至不要求建立匹配关系。如果匹配关系完全未知，则随机给定初始三维点坐标与摄像机参数；若已有初始匹配估计，则利用传统SfM算法估计相应参数。我们将原本一次性计算的过程转化为迭代求精的过程。即使匹配关系未知或初始匹配结果较差，通过迭代反馈的方法，最终仍可精确得到三维重建结果。

传统的SfM算法直接求解最大似然估计，在匹配关系未知时，需考虑所有可能的匹配关系。若图像数量为m，三维点个数为n，则总的匹配数目为n！^m，其随m、n的增大呈现爆炸性增长，精确求解最大似然估计变得不可行。而这里通过引入匹配向量进行迭代求解，解决了上述问题。

由于匹配存在互斥，即不存在某点对应多个点的关系，无法得到图像特征点与三维点之间匹配概率的准确解析式，这里利用MCMC统计方法得到其近似解，只要采样数目足够大，近似解可无限接近于真实值。然而传统的MCMC采样方法存在收敛速度慢且效率低下的缺点，我们通过避免选择原有匹配边，使得采样效率得到明显提升，在相同精度要求下，原始算法需要5000步左右采样步数才能达到稳定状态，而智能交换环仅需2000步左右即可。

传统方法利用BA算法求解三维点坐标与摄像机参数，BA算法是一种局部迭代算法，容易陷入局部最优，需要有良好的初值估计。这里采用可处理遮挡的射影因式分解算法，具有全局最优特性。通过对遮挡点进行填充，弥补了传统因式分解算法需要满足所有图像点可见的缺陷。

从整体引入了确定性退火策略进行全局优化。在迭代初始时MCMC算法能够在采样空间中更加均匀采样，不容易陷入局部极小点。随着迭代的深入，逐步降低温度，减小噪声方差σ，最终使得三维重建结果收敛到全局极小点。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

1、使用Harris角点检测算法提取图像的角点：

对西安钟楼拍摄48幅图像，并利用1-48之间的整数对这48幅图像进行编码标记，对每幅图像采用harris角点检测算法提取角点作为图像特征点，以角点集合U＝{u_ik|k∈1…n_i，i∈1…48}表示，其中k为图像中的角点序号，i为图像序号，n_i为第i幅图像中的角点数量，u_ik为第i幅图像中第k个角点的二维坐标。

2、生成三维点集合X与摄像机参数集合M：

比较各幅图像中的角点数量，得到最大角点数量N，这里N＝213。使用满足高斯分布N(0，1)的随机函数生成213个三维点，得到三维点集合X，X＝{x_i|i∈1…213}。48幅图像集合中每幅图像的摄像机参数m_i取为

得到摄像机参数集合M，M＝{m_i|i∈1…48}。

3、利用马尔可夫链蒙特卡罗(MCMC)方法估计图像特征点与三维点之间的匹配概率：

设置初始噪声方差σ＝20，利用基本投影公式计算213个三维点集合X在48幅图像上的重投影点集合H，H＝{h_ij|i∈1…48，j∈1…213}。利用式错误！未找到引用源。计算每幅图像中的角点与任意三维点之间的匹配权重。依次对每幅图像的角点与三维点之间的匹配关系进行MCMC采样，一次采样步骤如下：

步骤a：在213个三维点中为图像i中每个角点随机选择一个三维点，建立角点与三维点间的匹配边，得到匹配向量j_i ⁰，约束条件为：不同角点不可对应同一三维点；

步骤b：从图像i的角点集合u_i中随机选择一点u_ik，计算u_ik与任意三维点之间的转移概率：

$q^{j_i}(u_{\mathrm{ik}},x_j)=\left\{\begin{array}{ccc}\frac{\exp (-w(u_{\mathrm{ik}},x_j))}{{\mathrm{\Σ}}_{j\≠j_{\mathrm{ik}}}\exp (-w(u_{\mathrm{ik}},x_j))}& \mathrm{if}& j{\≠j}_{\mathrm{ik}}\\ 0& \mathrm{if}& j=j_{\mathrm{ik}}\end{array}\right.$

以概率

从213个三维点中选择一点x_j，以点x_j与u_ik建立匹配边；

步骤c：选择x_j的原有匹配边(x_j，u_ik′)对应的角点u_ik′作为下一次循环的图像角点；

步骤d：重复步骤b直至某次迭代中选择的图像角点在之前已被选择，得到由图像角点和三维点形成的闭环；

步骤e：删除在建立闭环过程中不属于闭环中的临时边，得到新的闭环C；

步骤f：删除闭环C中在步骤a的原有匹配关系，得到图像特征点与三维点之间新的匹配j′_i；

对于每幅图像，首先进行100次预采样，使得图像角点与三维点间的匹配状态达到稳定，在此过程中不保存生成的匹配向量。然后继续采样2000次，记录匹配向量，利用

计算图像角点与三维点之间的匹配概率f_ijk。

4、利用匹配概率对图像特征点进行加权平均得到虚拟测量点矩阵V：

根据计算虚拟测量点矩阵。

5、采用可处理遮挡的射影因式分解算法对虚拟测量点进行射影重构

步骤a：计算48幅图像集合中每两幅图像之间的基本矩阵F_ij与极点e_ij

步骤b：计算每个虚拟测量点的射影深度λ_ip，

利用v_ip＝λ_ipv_ip更新每个虚拟测量点；

步骤c：对虚拟测量点矩阵V进行400次采样。每次采样随机抽取矩阵的4列，得到一个列四元组。对每个列四元组求取其对应生成子空间B_t，t∈1…400。对每个B_t进行SVD分解：B_t＝S_tV_tD_t ^T，B_t⊥取为S_t的最后4列。对联合矩阵使用SVD分解：

B取为S的最后四列。将V中存在缺失元素的列表示为B的列的线性组合，即V_j＝k₁b₁+k₂b₂+k₃b₃+k₄b₄，利用列中已知元素计算相应系数，进而对列中未知元素进行填充，最终得到完全填充矩阵V^*；

步骤d：对V^*进行SVD分解，V^*＝USV^T，取S的前4列得到矩阵

摄像机参数

三维点坐标X＝V^T。利用BA算法对三维点坐标和摄像机参数进行优化，以得到的优化结果替换掉原有三维点坐标与摄像机参数；

6、加入确定性退火算法迭代求解全局最优射影重构结果：

对步骤3～5进行直接迭代的一个突出缺点就是只能保证收敛到局部极小值，这里引入退火策略。将噪声参数σ看作随温度T变化的退火因子，当T＝∞时，σ＝∞；当T＝0时，T＝1。初始时，σ取值较大(一般设为20-30个像素)，然后对步骤3～5进行循环迭代，σ在迭代过程中逐渐减小。参数减小包含多种策略，如线性策略、等比例策略等。线性策略是每次减小固定值，即σ^t+1＝σ^t-c；而等比例策略是将参数取为原来的一定比例(一般设为0.85～0.95)，即σ^t+1＝λσ^t。这里采用等比例退火策略。在每次迭代中，比较当前σ与给定阈值σ_min大小，若σ≤σ_min则迭代结束，进行后续步骤。否则置σ^t+1＝λσ^t，进入下一次迭代过程。

本实施例中重复步骤3～5，每次查看噪声方差σ的值，如果σ＞1，则令σ＝0.93σ，转步骤3进行下一次迭代；否则结束循环。

7、利用基于绝对对偶二次曲面的摄像机自标定算法将射影重构提升为度量重构：

计算绝对对偶二次曲面Q_∞ ^*，通过将Q_∞ ^*分解为

得到单应矩阵H，最终将射影重构变换到度量重构，得到欧式空间中三维场景重建的三维点坐标MH与摄像机参数H^-1X。

Claims (1)

1.一种基于统计模型的三维场景重建方法，其特征在于步骤如下：

步骤1使用Harris角点检测算法提取每幅图像中的角点：针对无序数字图像集合，采用Harris角点检测算法提取每幅图像的角点作为图像特征点，以角点集合U＝{u_ik|k∈1...n_i，i∈1...m}表示，其中k为图像中的角点序号，i为图像序号，n_i为第i幅图像中的角点数量，m为图像数量，u_ik为第i幅图像中第k个角点的二维坐标；

步骤2生成三维点集合X与摄像机参数集合M：采用满足高斯分布N(0，1)的随机函数生成N个三维点，得到三维点集合X＝{x_i|i∈1...n}，x_i为第i个三维点的坐标，N的值等于所有图像中最大的角点数量；

所述的摄像机参数集合M，M＝{m_i|i∈1...m}中每幅图像的摄像机参数初始值m_i取为 $\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right];$

步骤3利用马尔可夫链蒙特卡罗MCMC方法估计图像角点与三维点之间的匹配概率：利用基本投影公式计算得到三维点在所有图像的重投影点集合H＝{h_ij|j∈1...n，i∈1...m}，h_ij为第j个三维点在i幅图像上的投影点；计算图像角点u_ik和三维点x_j之间的权重：

$w(u_{\mathrm{ik}},x_j)=\log \left(\sqrt{2\mathrm{\π}}\mathrm{\σ}\right)+\frac{{\left|\right|u_{\mathrm{ik}}-h_{\mathrm{ij}}\left|\right|}^2}{{\mathrm{\σ}}^2}$

其中：σ为干扰图像的高斯噪声的方差，初始取为σ⁰；

然后利用MCMC统计模拟方法对图像角点u_ik与三维点x_j之间的匹配概率f_ijk进行概率估计：对于任意一幅图像i的图像角点与三维点匹配向量j_i＝{j_ik|k∈1...n_i}，匹配向量满足约束：k₁≠k₂

j_ik1≠j_ik2，图像i的采样过程如下：

步骤a：为图像i中每个角点随机选择一个三维点，建立角点与三维点间的匹配边，得到匹配向量j_i ⁰，约束条件为：不同角点不可对应同一三维点；

步骤b：从图像i的角点集合u_i中随机选择一点u_ik，计算u_ik与任意三维点之间的转移概率：

$q^{j_i}(u_{\mathrm{ik}},x_j)=\left\{\begin{array}{ccc}\frac{\exp (-w(u_{\mathrm{ik}},x_j))}{{\mathrm{\Σ}}_{j\≠j_{\mathrm{ik}}}\exp (-w(u_{\mathrm{ik}},x_j))}& \mathrm{if}& j\≠j_{\mathrm{ik}}\\ 0& \mathrm{if}& j=j_{\mathrm{ik}}\end{array}\right.$

以概率

从三维点集合X中选择一点x_j，以点x_j与u_ik建立匹配边；

步骤c：选择x_j的原有匹配边(x_j，u_ik′)对应的角点u_ik′作为下一次循环的图像角点；

步骤d：重复步骤b直至某次迭代中选择的图像角点在之前已被选择，得到由图像角点和三维点形成的闭环；

步骤e：删除在建立闭环过程中不属于闭环中的临时边，得到新的闭环C；

步骤f：删除闭环C中在步骤a的原有匹配关系，得到图像特征点与三维点之间新的匹配j_i′；

重复3000～5000次步骤a～f，保存每次采样生成的匹配向量j_i ^r，并计算匹配概率：

$f_{\mathrm{ijk}}\≈\frac{1}{R}\underset{r=1}{\overset{R}{\mathrm{\Σ}}}\mathrm{\δ}(j_{\mathrm{ik}}^r,j),$ 其中δ为克罗内克函数，若j_ik＝j，δ(j_ik，j)＝1；否则δ(j_ik，j)＝0；

步骤4：利用图像角点与三维点之间的匹配概率对图像角点进行加权平均得到虚拟测量点矩阵V

$v_{\mathrm{ij}}=\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{f}_{\mathrm{ijk}}{u}_{\mathrm{ik}},$ 其中：v_ij为V中第i行第j个元素；

步骤5采用可处理遮挡的射影因式分解算法对虚拟测量点进行射影重构：具体步骤如下：

步骤a：利用归一化8点算法计算任意两幅图像间的基本矩阵F_ij与极点e_ij；

步骤b：计算每个虚拟测量点v_ip的射影深度信息λ_ip

${\mathrm{\λ}}_{\mathrm{ip}}=\frac{\left(F_{\mathrm{ij}}^T{v}_{\mathrm{jp}}\right)\·(e_{\mathrm{ij}}\×v_{\mathrm{ip}})}{{\left|\right|e_{\mathrm{ij}}\×v_{\mathrm{jp}}\left|\right|}^2}{\mathrm{\λ}}_{\mathrm{jp}},$ 在得到所有射影深度信息后，利用v_ip＝λ_ipv_ip计算新的虚拟测量点；

步骤c：对虚拟测量矩阵V中的缺失部分进行线性拟合：首先对虚拟测量矩阵V进行T次采样，每次随机抽取4列，得到一个列四元组；对每个列四元组求取其生成子空间，得到B_t，t∈1...T；对B_t进行SVD分解：B_t＝S_tV_tD_t ^T，B_t ^⊥取为S_t的最后4列，其中⊥表示线性空间的补集；将每个B_t ^⊥按列组合，对得到的联合矩阵使用SVD分解：

B取为S的最后四列；将V中存在缺失元素的列表示为B的列的线性组合，即V_j＝k₁b₁+k₂b₂+k₃b₃+k₄b₄，利用列中已知元素计算相应系数，进而计算列中未知元素，最终得到完全矩阵V^*；

步骤d：对V^*进行SVD分解求取三维点坐标与摄像机参数：对V^*进行SVD分解，V^*＝USV^T，取S的前4列得到矩阵

摄像机参数 $M=U\hat{S},$ 三维点坐标X＝V^T；利用BA算法对三维点坐标和摄像机参数进行优化，以得到的优化结果替换掉原有三维点坐标与摄像机参数；

步骤6加入确定性退火算法迭代求解全局最优射影重构结果：重复步骤3～5进行循环迭代，σ在迭代过程中以等比例策略，将σ参数取为原σ的0.85～0.95，在每次迭代中，比较当前σ与给定阈值σ_min大小，若σ≤σ_min则迭代结束，进行后续步骤。否则置σ^t+1＝λσ^t，进入下一次迭代过程；

步骤7利用基于绝对对偶二次曲面的摄像机自标定算法将射影重构提升为度量重构：首先计算绝对对偶二次曲面Q_∞ ^*，然后将Q_∞ ^*分解为 $H\tilde{I}{H}^T(\tilde{I}=\left[\begin{array}{cc}{I}_{3\×3}& 0\\ 0^T& 0\end{array}\right]),$ 得到单应矩阵

$H=\left[\begin{array}{cc}K& 0\\ {-p}^{T}K& 1\end{array}\right]$

其中：p为无穷远平面的前三个坐标，将射影重构变换到度量重构，得到欧式空间中三维场景重建的三维结构M×H与摄像机参数H^-1X。