CN110533765B - 一种多次单应诱导的多平面结构物体的三维重建方法 - Google Patents

Abstract

本发明公开了一种多次单应诱导的多平面结构物体的三维重建方法，步骤包括：步骤1，对图像对的特征进行提取与匹配；步骤2，确定主平面；步骤3，对辅助平面进行估计；步骤4，对辅助平面进行校正；步骤5，确定模型平面区域，将多平面方程和其对应的内点集，由内点集在平面上构成凸集合，将该凸集合称为高置信区域；然后对所有图像进行线段检测，利用得到的线段和确定的平面间的交线来扩展高置信区域，以便确定每个平面的重建区域范围；最后对确定的每个平面纹理贴图，完成对低纹理场景的建筑物三维重建。本发明的方法，能够完成对低纹理场景的建筑物三维重建。

Description

一种多次单应诱导的多平面结构物体的三维重建方法

技术领域

本发明属于基于图像序列的建筑物三维重建技术领域，涉及一种多次单应诱导的多平面结构物体的三维重建方法。

背景技术

建筑物的三维重建是获取建筑物三维结构信息的重要手段，在城市规划、古建筑物保护、灾害监控通信设施建设和数字化城市等领域具有非常广泛的应用。如在工程估计上，对一些破损的建筑物，采用拍照重建，就可以快速地估计出修复建筑物所要的费用，提高了工程估计的高效性、灵活性和可靠性。

但是对于拍摄的建筑物场景中存在单一纹理和重复纹理的情况，现有的方法会因重建过程中所提取的特征点较少，致使重建的三维点云太过稀疏，无法很好地重构出完整的三维模型。

发明内容

本发明的目的是提供一种多次单应诱导的多平面结构物体的三维重建方法，解决了现有技术对于拍摄的建筑物场景，难以应对单一纹理和重复纹理现象，致使重建的三维点云太过稀疏，无法很好地重构出完整三维模型的问题。

本发明所采用的技术方案是，一种多次单应诱导的多平面结构物体的三维重建方法，按照以下步骤实施：

步骤1，对图像对的特征进行提取与匹配；

步骤2，确定主平面；

步骤3，对辅助平面进行估计；

步骤4，对辅助平面进行校正；

步骤5，确定模型平面区域，

通过步骤3和步骤4确定模型的多平面方程P_kc，kc＝1,2,...,N_c和其对应的内点集O_kc，kc＝1,2,...,N_c，由内点集在平面上构成凸集合，将该凸集合称为高置信区域；然后对所有图像进行线段检测，利用得到的线段和确定的平面间的交线来扩展高置信区域，以便确定每个平面的重建区域范围；最后对确定的每个平面纹理贴图，完成对低纹理场景的建筑物三维重建。

本发明的有益效果是，在利用点、线特征的基础上，结合点特征和线特征的优点以确定平面区域，使用更高层次的平面特征来描述建筑物进行重建，从多平面几何结构模型特征分析解决了重建场景中平面重建不完整的问题，尤其是针对环境中存在单一纹理和重复纹理的情况，本发明算法效果更加明显。

附图说明

图1是本发明实施例的建筑物重建流程框图；

图2是本发明方法实施例的特征点中多次单应诱导平面结果，其中，图2a是图像对极几何约束剔除误匹配，图2b是图像对第一次估计单应变换内点，图2c是图像对第二次估计单应变换内点；

图3是本发明方法实施例的高置信区域扩展平面；

图4是本发明方法实施例的盒状物模型重建结果，图4a和4b分别为重建模型的不同视角展示。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

参照图1，本发明的方法，按照以下步骤实施：

步骤1，对图像对的特征进行提取与匹配，

1.1)设输入的图像序列为{I_ts＝[Im_ts(x,y)]_m×n|ts＝1,2,...,T}，ts是指图像序数，T为图像序列中的图像数量，(x,y)代表当前图像像素横纵坐标值，m和n是图像的长宽尺度；依次取相邻的两幅图像，即I_t＝[Im_t(x,y)]_m×n和I_t+1＝[Im_t+1(x,y)]_m×n，t＝1,2,...,T-1，对该相邻两幅图像分别提取其SIFT特征(SIFT特征的提取方法是公开方法，可在相关的技术文献中检索到)，设I_t的SIFT特征点集为S_t＝{s_t(1),s_t(2),...,s_t(n_t)}，t＝1,2,...,T-1，这里的s_t(n_t)代表每个特征点；I_t+1的SIFT特征点集为S_t+1＝{s_t+1(1),s_t+1(2),...,s_t+1(n_t+1)}，t＝1,2,...,T-1，按照邻近距离约束匹配的算法得到该相邻两幅图像一一对应的特征点匹配集合，计算式如下：

其中，th是阈值，取经验值，优选取值为th＝0.8；k_t是指两个特征向量之间的距离，公式(1)代表点s_t(i_t)的次邻近距离，公式(2)为其最邻近距离，p_t是次近邻与最近邻距离的比值，

如果公式(3)的结果p_t(i_t)＜th，则表明s_t(i_t)与s_t+1(k_t*)配对成功，将s_t+1(k*)从S_t+1＝{s_t+1(1),s_t+1(2),...,s_t+1(n_t+1)}中取出，继续S_t＝{s_t(1),s_t(2),...,s_t(n_t)}中下一个特征点的匹配；以此类推，对I_t的特征点集S_t中的所有特征点，按照公式(1)和公式(2)及公式(3)进行匹配，得到配对成功的特征点集，设为

1.2)采用RANSAC算法(RANSAC算法为现有公开技术，可在相关技术文献中检索到)，去除特征点集

中误匹配的特征点，得到正确匹配的特征点集，将其称为内点集合，记作

t＝1,2,...,T-1，其中，

和

是归一化坐标，

j＝1，2,...,m_t，t＝1,2,...,T-1；

1.3)采用三角测量算法(三角测量算法为现有公开技术，可在相关文献中检索到)计算出内点集合

对应在三维空间中的点集Ω_t＝{(x_t(j),y_t(j),z_t(j))，j＝1，2,...,m_t}，t＝1,2,...,T-1；

步骤2，确定主平面，

2.1)求主平面的单应变换矩阵H₀，

将步骤1中得到的T-1个内点集合

t＝1,2,...,T-1的所有对应点汇集在一起，构建总的内点集合为：

其中ks代表内点数，M是内点的总个数，即M＝m₁+m₂+...+m_t+...+m_T-1；对应在三维空间中的T-1个点集Ω_t＝{(x_t(j)，y_t(j)，z_t(j))，j＝1，2，...，m_t)}，t＝1,2,...,T-1所有对应点汇集在一起，构建总的内点集合为：

Ω＝{Ω₁,Ω₂,...,Ω_t,...,Ω_T-1}＝{(x(ks),y(ks),z(ks))，ks＝1，2,...,m₁,...,M}；

求解公式(4)中的主平面单应变换矩阵H₀：

其中,

将公式(4)展开，则有：

(5)

此处A₀是一个2M×9的矩阵，求A₀的SVD奇异值分解，对应于最小特征值的单位特征矢量便是解h₀，将求得的解h₀带入H₀，得到单应变换矩阵H₀；

2.2)拟合主平面，

将求解得到的单应变换矩阵H₀带入公式(4)，挑出内点集

中，满足公式(4)的点，设为

M₀是满足求解主平面的内点个数；对应的三维空间上的点集设为Ω₀＝{(x(k₀),y(k₀),z(k₀))，k₀＝1，2,...,M₀}，用Ω₀拟合空间上的主平面P₀，主平面P₀的表达式如下：

a₀·x(k₀)+b₀·y(k₀)+c₀·z(k₀)-d₀＝p(k₀)·P₀＝0 k₀＝1，2,...,M₀ (6)

其中，P₀＝[a₀,b₀,c₀,d₀]^T是平面方程的系数，a0,b₀,c₀,d₀均为参变量，都是为了确定平面指标，p(k₀)＝(x(k₀),y(k₀),z(k₀))，k₀＝1,2,...,m₀是Ω₀中的点；由于公式(6)是超定方程组，所以在此采用最小二乘法(最小二乘法为现有公开技术，可在相关学术论文中检索到)求解得到主平面方程P₀＝[a₀，b₀，c₀，d₀]^T，本公式中的T是转置符号；

步骤3，对辅助平面进行估计，

在步骤2中确定主平面方程P₀＝[a₀,b₀,c₀,d₀]^T之后，接下来就是对其他的辅助平面进行估计：

3.1)令确定平面数k_p＝1；

3.2)求剩余内点集：

其中

即从

中去除已经完成拟合平面的内点以及在拟合平面附近的内点，得到剩余内点集

统计

中的内点个数N_p，如果N_p＜ε，ε是经验值，此处优选值为5，转步骤4；否则，转步骤3.3)；

3.3)求辅助平面的单应变换矩阵

对

的内点按照如下公式(7)计算h_kp：

其中，N_p是剩余内点集

中的内点个数；

此处A_kp是一个2N_p×9的矩阵，求A_kp的SVD奇异值分解，对应于最小特征值的单位特征矢量便是解

将求得的h_kp带入

由此得到单应变换矩阵

则有：

3.4)拟合辅助平面

将求解得到的单应矩阵

带入如下公式(9)：

将内点集

中满足公式(9)的点挑出来，构成辅助平面的内点集，设为

此处M_p表示满足公式(9)的内点个数，对应的三维空间上的点集设为

用

拟合空间上的主平面

表达式如下：

其中，

是平面方程的系数，p(k_p)＝(x(k_p),y(k_p),z(k_p))，k_p＝1,2,...,M_p是

中的点，M_p是满足平面方程的点总个数；

由于公式(10)是超定方程组，在此采用最小二乘法(最小二乘法为现有公开技术，可在相关学术论文中检索到)求解得到主平面方程为

3.5)令k_p＝k_p+1，转步骤3.2)重新循环；

步骤4，对辅助平面进行校正，

基于前述从步骤2得到的一个主平面，以及从步骤3得到的N_c-1个辅助平面，考虑到噪声对平面方程的影响，需要对其进行校正，具体过程是：

4.1)计算主平面与辅助平面之间的夹角，

i_c＝0,1,...,N_c-1，j_c＝i_c+1,i_c+2,...,N_c

其中，

代表主平面，

代表辅助平面；

4.2)判断辅助平面是否需要校正，

4.2.1)如果θ(i_c,j_c)接近垂直，则有：

|θ(i_c,j_c)-90°|＜δ (12)

其中，δ是偏差量，取经验值，根据数据受干扰程度的大小选择，在此优选值为δ＝1.5°；则转步骤4.3)进行校正；

4.2.2)如果θ(i_c,j_c)接近平行，则有：

|θ(i_c,j_c)-180°|＜δ (13)

其中，δ是偏差量，取经验值，根据数据受干扰程度的大小选择，在此优选值为δ＝1.5°；

4.2.3)如果两个平面既不是垂直或者也不是平行关系，即θ(i_c,j_c)即不满足公式(12)也不满足公式(13)，则不进行校正，直接转入步骤5；

4.3)对互相垂直平面进行校正，

为描述方便起见，不失一般性，设第一个近似垂直主平面的辅助平面为P₁，经过校正之后为

与主平面P₀垂直,计算式如下：

其中，

是校正后辅助平面参数，确定平面参数有

是主平面参数；

是第一个辅助平面对应的三维空间点集合Ω₁中的点.

为校正后第一个辅助平面后的方程组系数，求解该最优化问题即可得到；如果后面再没有与主平面垂直的辅助平面，则转入步骤5；(此处的第一辅助平面即是前面步骤3和步骤4描述的辅助平面，他们都是用P₁，P₂依次表示的)，

如果还存在与主平面P₀和第一辅助平面P₁均垂直的平面P₂，则将其称为第二辅助平面P₂，求解其校正平面方程

的优化模型如下：

其中，

是第二个辅助平面对应的三维空间点集合Ω₂中的点；

为校正后第二个辅助平面的方程组系数，求解该最优化问题即可得到；然后转入步骤5；

4.4)对互相平行的平面进行校正，

为描述方便起见，不失一般性，以主平面为例，设待校正平面P_i与主平面P₀平行，则待校正平面P_i的法向量和主平面P₀的法向量相同，即

所以只需要确定

即可，平行于主平面的辅助平面(这里需要校正的辅助平面都可以称作待校正平面)校正模型表达式如下：

其中，

是第一个辅助平面对应的三维空间点集合Ω_ι中的点，

为校正后第一个辅助平面后的方程组系数，求解该最优化问题即可得到；

将所有相互平行的平面均按照上述的方式校正之后，转入步骤5；

实施例中，图2a是图像对极几何约束剔除误匹配，图2b是图像对第一次估计单应变换内点，图2c是图像对第二次估计单应变换内点。

步骤5，确定模型平面区域，

参照图3，通过步骤3和步骤4确定模型的多平面方程P_kc，kc＝1,2,...,N_c和其对应的内点集O_kc，kc＝1,2,...,N_c，由内点集在平面上构成凸集合，将该凸集合称为高置信区域；然后对所有图像进行线段检测，利用得到的线段和确定的平面间的交线来扩展高置信区域，以便确定每个平面的重建区域范围；最后对确定的每个平面纹理贴图，

至此，完成对低纹理场景的建筑物三维重建。

传统的重建方法对于环境纹理单一的环境不能实现重建。本发明算法通过更高级的面匹配特征，重建出建筑物模型，实施例见图4a和4b，明显可见，其中重建出的建筑物平面的完整度和准确性高。

Claims (5)

1.一种多次单应诱导的多平面结构物体的三维重建方法，其特征在于，按照以下步骤实施：

步骤1，对图像对的特征进行提取与匹配；

步骤2，确定主平面；

步骤3，对辅助平面进行估计；

步骤4，对辅助平面进行校正；

步骤5，确定模型平面区域，

通过步骤3和步骤4确定模型的多平面方程P_kc，kc＝1,2,...,N_c和其对应的内点集O_kc，kc＝1,2,...,N_c，由内点集在平面上构成凸集合，将该凸集合称为高置信区域；然后对所有图像进行线段检测，利用得到的线段和确定的平面间的交线来扩展高置信区域，以便确定每个平面的重建区域范围；最后对确定的每个平面纹理贴图，完成对低纹理场景的建筑物三维重建。

2.根据权利要求1所述的多次单应诱导的多平面结构物体的三维重建方法，其特征在于：所述的步骤1中，具体过程是，

1.1)设输入的图像序列为{I_ts＝[Im_ts(x,y)]_m×n|ts＝1,2,...,T}，ts是指图像序数，T为图像序列中的图像数量，(x,y)代表当前图像像素横纵坐标值，m和n是图像的长宽尺度；依次取相邻的两幅图像，即I_t＝[Im_t(x,y)]_m×n和I_t+1＝[Im_t+1(x,y)]_m×n，t＝1,2,...,T-1，对该相邻两幅图像分别提取其SIFT特征，设I_t的SIFT特征点集为S_t＝{s_t(1),s_t(2),...,s_t(n_t)}，t＝1,2,...,T-1，这里的s_t(n_t)代表每个特征点；I_t+1的SIFT特征点集为S_t+1＝{s_t+1(1),s_t+1(2),...,s_t+1(n_t+1)}，t＝1,2,...,T-1，按照邻近距离约束匹配的算法得到该相邻两幅图像一一对应的特征点匹配集合，计算式如下：

其中，th是阈值，k_t是指两个特征向量之间的距离，公式(1)代表点s_t(i_t)的次邻近距离，公式(2)为其最邻近距离，p_t是次近邻与最近邻距离的比值，

如果公式(3)的结果p_t(i_t)＜th，则表明s_t(i_t)与s_t+1(k_t*)配对成功，将s_t+1(k*)从S_t+1＝{s_t+1(1),s_t+1(2),...,s_t+1(n_t+1)}中取出，继续S_t＝{s_t(1),s_t(2),...,s_t(n_t)}中下一个特征点的匹配；以此类推，对I_t的特征点集S_t中的所有特征点，按照公式(1)和公式(2)及公式(3)进行匹配，得到配对成功的特征点集，设为

1.2)采用RANSAC算法，去除特征点集

中误匹配的特征点，得到正确匹配的特征点集，将其称为内点集合，记作

其中，

和

是归一化坐标，

1.3)采用三角测量算法计算出内点集合

对应在三维空间中的点集Ω_t＝{(x_t(j),y_t(j),z_t(j))，j＝1，2,...,m_t}，t＝1,2,...,T-1。

3.根据权利要求2所述的多次单应诱导的多平面结构物体的三维重建方法，其特征在于：所述的步骤2中，具体过程是，

2.1)求主平面的单应变换矩阵H₀，

将步骤1中得到的T-1个内点集合

的所有对应点汇集在一起，构建总的内点集合为：

其中ks代表内点数，M是内点的总个数，即M＝m₁+m₂+...+m_t+...+m_T-1；

对应在三维空间中的T-1个点集Ω_t＝{(x_t(j)，y_t(j)，z_t(j))，j＝1，2，...，m_t}，t＝1,2,...,T-1所有对应点汇集在一起，构建总的内点集合为：Ω＝{Ω₁,Ω₂,...,Ω_t,...,Ω_T-1}＝{(x(ks),y(ks),z(ks))，ks＝1，2,...,m₁,...,M}；

求解公式(4)中的主平面单应变换矩阵H₀：

其中,

将公式(4)展开，则有：

此处A₀是一个2M×9的矩阵，求A₀的SVD奇异值分解，对应于最小特征值的单位特征矢量便是解h₀，将求得的解h₀带入H₀，得到单应变换矩阵H₀；

2.2)拟合主平面，

将求解得到的单应变换矩阵H₀带入公式(4)，挑出内点集

中，满足公式(4)的点，设为

M₀是满足求解主平面的内点个数；对应的三维空间上的点集设为Ω₀＝{(x(k₀),y(k₀),z(k₀))，k₀＝1，2,...,M₀}，用Ω₀拟合空间上的主平面P₀，主平面P₀的表达式如下：

a₀·x(k₀)+b₀·y(k₀)+c₀·z(k₀)-d₀＝p(k₀)·P₀＝0 k₀＝1，2,...,M₀ (6)

其中，P₀＝[a₀,b₀,c₀,d₀]^T是平面方程的系数，a₀,b₀,c₀,d₀均为参变量，都是为了确定平面指标，p(k₀)＝(x(k₀),y(k₀),z(k₀))，k₀＝1,2,...,m₀是Ω₀中的点；

由于公式(6)是超定方程组，所以在此采用最小二乘法求解得到主平面方程P₀＝[a₀,b₀,c₀,d₀]^T，本公式中的T是转置符号。

4.根据权利要求3所述的多次单应诱导的多平面结构物体的三维重建方法，其特征在于：所述的步骤3中，具体过程是，

3.1)令确定平面数k_p＝1；

3.2)求剩余内点集：

其中

即从

中去除已经完成拟合平面的内点以及在拟合平面附近的内点，得到剩余内点集

统计

中的内点个数N_p，如果N_p＜ε，ε是经验值，此处优选值为5，转步骤4；否则，转步骤3.3)；

3.3)求辅助平面的单应变换矩阵

对

的内点按照如下公式(7)计算h_kp：

其中，N_p是剩余内点集

中的内点个数；

此处A_kp是一个2N_p×9的矩阵，求A_kp的SVD奇异值分解，对应于最小特征值的单位特征矢量便是解

将求得的h_kp带入

由此得到单应变换矩阵

则有：

3.4)拟合辅助平面

将求解得到的单应矩阵

带入如下公式(9)：

将内点集

中满足公式(9)的点挑出来，构成辅助平面的内点集，设为

此处M_p表示满足公式(9)的内点个数，对应的三维空间上的点集设为

用

拟合空间上的主平面

表达式如下：

其中，

是平面方程的系数，p(k_p)＝(x(k_p),y(k_p),z(k_p))，k_p＝1,2,...,M_p是

中的点，M_p是满足平面方程的点总个数；

由于公式(10)是超定方程组，在此采用最小二乘法求解得到主平面方程为

3.5)令k_p＝k_p+1，转步骤3.2)重新循环。

5.根据权利要求4所述的多次单应诱导的多平面结构物体的三维重建方法，其特征在于：所述的步骤4中，

基于前述从步骤2得到的一个主平面，以及从步骤3得到的N_c-1个辅助平面，进行校正的具体过程是：

4.1)计算主平面与辅助平面之间的夹角，

其中，

代表主平面，

代表辅助平面；

4.2)判断辅助平面是否需要校正，

4.2.1)如果θ(i_c,j_c)接近垂直，则有：

|θ(i_c,j_c)-90°|＜δ (12)

其中，δ是偏差量，转步骤4.3)进行校正；

4.2.2)如果θ(i_c,j_c)接近平行，则有：

|θ(i_c,j_c)-180°|＜δ (13)

其中，δ是偏差量；

4.2.3)如果两个平面既不是垂直或者也不是平行关系，即θ(i_c,j_c)即不满足公式(12)也不满足公式(13)，则不进行校正，直接转入步骤5；

4.3)对互相垂直平面进行校正，

为描述方便起见，不失一般性，设第一个近似垂直主平面的辅助平面为P₁，经过校正之后为

与主平面P₀垂直,计算式如下：

其中，

是校正后辅助平面参数，确定平面参数有

[a₀ b₀ c₀ d₀]＝p₀是主平面参数；

是第一个辅助平面对应的三维空间点集合Ω₁中的点.

为校正后第一个辅助平面后的方程组系数，求解该最优化问题即可得到；如果后面再没有与主平面垂直的辅助平面，则转入步骤5；

如果还存在与主平面P₀和第一辅助平面P₁均垂直的平面P₂，则将其称为第二辅助平面P₂，求解其校正平面方程

的优化模型如下：

其中，

是第二个辅助平面对应的三维空间点集合Ω₂中的点；

为校正后第二个辅助平面的方程组系数，求解该最优化问题即可得到；然后转入步骤5；

4.4)对互相平行的平面进行校正，

设待校正平面P_i与主平面P₀平行，则待校正平面P_i的法向量和主平面P₀的法向量相同，即

所以只需要确定

即可，平行于主平面的辅助平面校正模型表达式如下：

其中，

是第一个辅助平面对应的三维空间点集合Ω_ι中的点，

为校正后第一个辅助平面后的方程组系数，求解该最优化问题即可得到；

将所有相互平行的平面均按照上述的方式校正之后，转入步骤5。

Images