CN103247075A

CN103247075A - Variational mechanism-based indoor scene three-dimensional reconstruction method

Info

Publication number: CN103247075A
Application number: CN201310173608XA
Authority: CN
Inventors: 贾松敏; 王可; 李雨晨; 李秀智
Original assignee: Beijing University of Technology
Current assignee: Beijing University of Technology
Priority date: 2013-05-13
Filing date: 2013-05-13
Publication date: 2013-08-14
Anticipated expiration: 2033-05-13
Also published as: CN103247075B

Abstract

The invention belongs to the intersecting field of computer vision and intelligent robot, and discloses a reconstruction method of a large-scale indoor scene based on a variational mechanism, including: step 1, obtaining calibration parameters of a camera, and establishing a distortion correction model; step 2, establishing Camera pose description and camera projection model; step 3, use the SFM-based monocular SLAM algorithm to realize camera pose estimation; step 4, establish a depth map estimation model based on a variational mechanism, and solve the model; step 5, establish a key The frame selection mechanism realizes the update of the 3D scene. The present invention uses RGB cameras to acquire environmental data, and proposes a depth map generation method based on a variational mechanism for the use of high-precision monocular positioning algorithms, which realizes large-scale rapid indoor 3D scene reconstruction and effectively solves the 3D reconstruction algorithm Cost and timeliness issues.

Description

3D Reconstruction Method of Indoor Environment Based on Variational Mechanism

技术领域technical field

本发明属于计算机视觉与智能机器人的交叉领域，涉及一种室内环境三维重建技术，尤其涉及一种基于变分机制的大范围室内场景的重建方法。The invention belongs to the intersecting field of computer vision and intelligent robots, relates to a three-dimensional reconstruction technology of an indoor environment, and in particular relates to a reconstruction method of a large-scale indoor scene based on a variational mechanism.

技术背景technical background

随着同时定位与地图创建（Simultaneous Localization And Mapping，SLAM）研究的不断深入，环境三维立体化建模已逐步成为该领域研究热点，引起众多学者的关注。G.Klein等于2007年在增强现实（AR）领域首先提出同时定位与地图创建（Parallel Tracking and Mapping,PTAM）的概念，以解决环境实时建模问题。PTAM将摄像机跟踪与地图生成划分为两个独立线程，利用FastCorner方法更新检测特征点的同时，采用最优的局部与全局光束平差法(Bundle Adjustment,BA)，不断实现相机位姿与三维特征点地图的更新。该方法基于稀疏点云建立了环境三维地图，但该地图缺乏对环境的直观三维描述。Pollefeys等人通过多传感器融合实现了大型室外场景的三维重建。但该方法存在计算的高复杂性以及对噪音敏感等缺点。目前在实时跟踪和稠密环境模型重构方面也有了一些尝试性的进步，但是仅仅局限于一些简单物体的重构，并且只能在特定约束条件下可以获得较高的精度。Richard A.Newcombe等人，利用基于SFM（Structure from Moving）的SLAM算法获取空间稀疏特征点云，采用多尺度径向基插值，运用图形图像学中隐式曲面多边形化方法，构造三维空间初始化网格地图，并结合场景流约束与高精度TV-L1光流算法更新网格顶点坐标，以达到逼近真实场景的目的。该算法能获取高精度的环境模型，但由于其算法复杂度较高，在两个图形硬件处理器（GPU）加速情况下，处理一帧图像仍需花费几秒钟的时间。With the continuous deepening of Simultaneous Localization And Mapping (SLAM) research, three-dimensional environmental modeling has gradually become a research hotspot in this field, attracting the attention of many scholars. In 2007, G. Klein et al. first proposed the concept of simultaneous positioning and map creation (Parallel Tracking and Mapping, PTAM) in the field of augmented reality (AR) to solve the problem of real-time modeling of the environment. PTAM divides camera tracking and map generation into two independent threads, uses the FastCorner method to update the detection feature points, and uses the optimal local and global bundle adjustment (Bundle Adjustment, BA) to continuously realize the camera pose and 3D features Point map updates. This method builds a 3D map of the environment based on sparse point clouds, but the map lacks an intuitive 3D description of the environment. Pollefeys et al. achieved 3D reconstruction of large outdoor scenes through multi-sensor fusion. However, this method has disadvantages such as high computational complexity and sensitivity to noise. At present, there have been some tentative progress in real-time tracking and reconstruction of dense environment models, but they are only limited to the reconstruction of some simple objects, and can only obtain high accuracy under certain constraints. Richard A. Newcombe et al. used the SFM (Structure from Moving) SLAM algorithm to obtain spatially sparse feature point clouds, adopted multi-scale radial basis interpolation, and used the implicit surface polygonization method in graphics and imaging to construct a three-dimensional space initialization network. Grid map, combined with scene flow constraints and high-precision TV-L1 optical flow algorithm to update grid vertex coordinates, in order to achieve the purpose of approaching the real scene. This algorithm can obtain a high-precision environment model, but due to its high algorithm complexity, it still takes a few seconds to process one frame of image under the acceleration of two graphics hardware processors (GPU).

发明内容Contents of the invention

针对现有技术中存在的上述问题，本发明提供了一种基于变分机制的快速三维重建方法，以实现在室内复杂环境下的三维建模。该方法保证环境信息的同时降低了所需处理数据量，能实现大范围的快速室内三维场景重建。有效地解决了三维重建算法成本与实时性问题，提高了重建精度。Aiming at the above-mentioned problems in the prior art, the present invention provides a fast 3D reconstruction method based on a variational mechanism to realize 3D modeling in complex indoor environments. This method reduces the amount of data to be processed while ensuring the environmental information, and can realize rapid indoor three-dimensional scene reconstruction in a large range. It effectively solves the cost and real-time problems of the 3D reconstruction algorithm, and improves the reconstruction accuracy.

本发明采用的技术方案如下：The technical scheme that the present invention adopts is as follows:

利用PTAM算法作为相机位姿估计手段，并在关键帧处选取适当图像序列构造基于变分模式的深度图估计能量函数，运用原始对偶算法优化上述能量函数，实现在当前关键帧处环境深度图的获取。由于该算法利用邻近帧信息构造能量函数，且有效利用了特定视角坐标系统间的关联性，以及摄像机透视投影变换关系，使得数据项蕴含了多视成像约束，降低了算法模型求解的计算复杂度。在统一计算架构下，本发明利用图形加速硬件实现了算法的并行优化，有效提高了算法实时性。The PTAM algorithm is used as the camera pose estimation method, and an appropriate image sequence is selected at the key frame to construct a depth map estimation energy function based on the variational mode, and the primal dual algorithm is used to optimize the above energy function to realize the environment depth map at the current key frame. Obtain. Since the algorithm uses adjacent frame information to construct energy functions, and effectively utilizes the correlation between coordinate systems of specific viewing angles and the camera perspective projection transformation relationship, the data items contain multi-view imaging constraints, which reduces the computational complexity of the algorithm model solution . Under the unified computing framework, the invention realizes the parallel optimization of the algorithm by using the graphics acceleration hardware, and effectively improves the real-time performance of the algorithm.

一种基于变分机制的室内环境三维重建的方法，其特征在于包括以下步骤：A method for three-dimensional reconstruction of an indoor environment based on a variational mechanism, characterized in that it comprises the following steps:

步骤一，获取相机的标定参数，并建立畸变矫正模型。Step 1, obtain the calibration parameters of the camera, and establish a distortion correction model.

在计算机视觉应用中，通过相机成像的几何模型，有效建立图像中像素点与空间三维点之间的映射关系。构成相机模型的几何参数须通过实验与计算才能得到，求解上述参数的过程就称之为相机标定。在本发明中相机参数的标定是非常关键的环节，标定参数的精度直接影响最终结果三维地图的准确性。In computer vision applications, through the geometric model of camera imaging, the mapping relationship between pixels in the image and three-dimensional points in space is effectively established. The geometric parameters that constitute the camera model must be obtained through experiments and calculations. The process of solving the above parameters is called camera calibration. In the present invention, the calibration of the camera parameters is a very critical link, and the accuracy of the calibration parameters directly affects the accuracy of the final three-dimensional map.

相机标定的具体过程为：The specific process of camera calibration is as follows:

（1）打印一张棋盘模板。本发明采用一张A4纸，棋盘的间隔为0.25cm。(1) Print a chessboard template. The present invention adopts a piece of A4 paper, and the interval of chessboard is 0.25cm.

（2）从多个角度拍摄棋盘。拍摄时，应尽量让棋盘占满屏幕，并保证棋盘的每一个角都在屏幕中，一共拍摄6张模板图片。(2) Photograph the chessboard from multiple angles. When shooting, try to make the chessboard fill the screen as much as possible, and ensure that every corner of the chessboard is on the screen, and shoot a total of 6 template pictures.

（3）检测出图像中的特征点，即棋盘的每一个黑色交叉点。(3) Detect the feature points in the image, that is, every black intersection of the chessboard.

（4）求取相机的内部参数，方法如下：(4) Obtain the internal parameters of the camera, the method is as follows:

RGB相机标定参数主要为相机内参。相机的内参矩阵K为：RGB camera calibration parameters are mainly camera internal parameters. The internal parameter matrix K of the camera is:

$K K = = [\begin{matrix} {f f}_{u u} & 00 & {u u}_{00} \\ 00 & {f f}_{v v} & {v v}_{00} \\ 00 & 00 & 11 \end{matrix}]$

式中，u、v为相机平面坐标轴，(u₀,v₀)是相机像平面中心坐标，(f_u,f_v)是相机的焦距。In the formula, u and v are the coordinate axes of the camera plane, (u ₀ , v ₀ ) is the center coordinate of the camera image plane, and (f _u , f _v ) is the focal length of the camera.

根据标定参数，RGB图像中点与三维空间点的映射关系如下：RGB图像中点p＝(u,v)在相机坐标系下的坐标P_3D＝(x,y,z)表示为：According to the calibration parameters, the mapping relationship between the point in the RGB image and the point in the three-dimensional space is as follows: the coordinate P _3D = (x, y, z) of the point p in the RGB image = (u, v) in the camera coordinate system is expressed as:

$\{\begin{matrix} x x = = ((u u - - {u u}_{00})) * * z z / / {f f}_{u u} \\ y the y = = ((v v - - {v v}_{00})) * * z z / / {f f}_{v v} \\ z z = = d d \end{matrix}$

式中，d表示深度图像中点p的深度值。where d represents the depth value of point p in the depth image.

本发明中相机坐标系如图2所示，向下为y轴正方向，向前为z轴正方向，向右为x正方向。将相机的起始点位置设定为世界坐标系原点，世界坐标系的X、Y、Z方向与相机的定义相同。In the present invention, the camera coordinate system is shown in FIG. 2 , downward is the positive direction of the y-axis, forward is the positive direction of the z-axis, and rightward is the positive direction of the x-axis. Set the starting point of the camera as the origin of the world coordinate system. The X, Y, and Z directions of the world coordinate system are the same as those defined by the camera.

FOV（Field of Viewer）相机矫正模型为：The FOV (Field of Viewer) camera correction model is:

${u u}_{d d} = = [\begin{matrix} {u u}_{00} \\ {v v}_{00} \end{matrix}] + + [\begin{matrix} {f f}_{u u} & 00 \\ 00 & {f f}_{v v} \end{matrix}] \frac{{r r}_{d d}}{{r r}_{u u}} {x x}_{u u}$

${r r}_{d d} = = \frac{11}{ω ω} arctan arctan (({22 r r}_{u u} tan the tan \frac{ω ω}{22}))$

${r r}_{u u} = = \frac{tan the tan (({r r}_{d d} ω ω))}{22 tan the tan \frac{ω ω}{22}}$

式中，x_u为z=1面的像素坐标，u_d为原始图像中像素坐标，ω为FOV相机畸变系数。In the formula, x _u is the pixel coordinates of z=1 plane, u _d is the pixel coordinates in the original image, and ω is the FOV camera distortion coefficient.

步骤二，建立相机位姿描述与相机投影模型。Step 2, establish the camera pose description and camera projection model.

在已建立起的世界坐标系下，相机位姿可以表示为如下矩阵：In the established world coordinate system, the camera pose can be expressed as the following matrix:

${T T}_{cw cw} = = [\begin{matrix} {R R}_{cw cw} & {t t}_{cw cw} \\ 00 & 11 \end{matrix}]$

式中，“cw”表示从世界坐标系到当前相机坐标系，T_cw∈SE(3)，SE(3)为刚体的旋转平移变换空间。T_cw可由如下六元组μ＝(μ₁,μ₂,μ₃,μ₄,μ₅,μ₆)表示，即：In the formula, "cw" means from the world coordinate system to the current camera coordinate system, T _cw ∈ SE(3), SE(3) is the rotation-translation transformation space of the rigid body. T _cw can be represented by the following six-tuple μ=(μ ₁ , μ ₂ , μ ₃ , μ ₄ , μ ₅ , μ ₆ ), namely:

${T T}_{cw cw} = = exp exp ((\overset{^^}{μ μ}))$

$\overset{^^}{u u} = = [\begin{matrix} 00 & {μ μ}_{66} & - - {μ μ}_{55} & {μ μ}_{11} \\ {μ μ}_{66} & 00 & {μ μ}_{44} & {μ μ}_{22} \\ {μ μ}_{55} & - - {μ μ}_{44} & 00 & {μ μ}_{33} \\ 00 & 00 & 00 & 00 \end{matrix}]$

式中，μ₁,μ₂,μ₃分别为Kinect在全局坐标系下的平移量，μ₄,μ₅,μ₆表示局部坐标系下坐标轴的旋转量。In the formula, μ ₁ , μ ₂ , μ ₃ are the translation amount of Kinect in the global coordinate system respectively, and μ ₄ , μ ₅ , μ ₆ represent the rotation amount of the coordinate axis in the local coordinate system.

相机的位姿T_cw建立了当前坐标系下空间点云坐标p_c到世界坐标p_w的变换关系，即：The pose T _cw of the camera establishes the transformation relationship from the spatial point cloud coordinates p _c to the world coordinates p _w in the current coordinate system, namely:

p_c＝T_cwp_w p _c =T _cw p _w

在当前标系下，三维空间点云到z=1平面上投影定义为：Under the current coordinate system, the projection of the 3D space point cloud onto the z=1 plane is defined as:

π(p)＝(xz,yz)^T π(p)=(xz,yz) ^T

式中，p∈R³为三维空间点，x，y，z为该点的坐标值。根据当前坐标点深度值d，利用逆向投影法确定当前空间三维点坐标p，其坐标关系可表示为：In the formula, p∈R ³ is a point in three-dimensional space, and x, y, z are the coordinate values of the point. According to the depth value d of the current coordinate point, use the reverse projection method to determine the coordinate p of the three-dimensional point in the current space, and its coordinate relationship can be expressed as:

π^-1(u,d)＝dK^-1uπ ^-1 (u,d) = dK ^-1 u

步骤三，利用基于SFM的单目SLAM算法实现相机位姿估计。Step 3, use the SFM-based monocular SLAM algorithm to realize camera pose estimation.

目前，单目视觉SLAM算法主要包括基于滤波与SFM(Structure from Moving)的SLAM算法。本发明采用PTAM算法实现对相机的定位。该算法是一种基于SFM的单目视觉SLAM方法，通过将系统划分为相机跟踪与地图创建两个独立的线程。在相机跟踪线程，系统利用相机获取当前环境纹理信息，并构建四层高斯图像金字塔，运用FAST-10角点检测算法提取当前图像中特征信息，采用块匹配的方式建立角点特征间的数据关联。在此基础上，根据当前投影误差，建立位姿估计模型实现相机的精确定位，并结合特征匹配信息与三角测量算法生成当前三维点云地图。相机位姿估计的具体过程为：At present, monocular vision SLAM algorithms mainly include SLAM algorithms based on filtering and SFM (Structure from Moving). The invention adopts the PTAM algorithm to realize the positioning of the camera. This algorithm is a SFM-based monocular vision SLAM method, which creates two independent threads by dividing the system into camera tracking and map. In the camera tracking thread, the system uses the camera to obtain the current environmental texture information, and constructs a four-layer Gaussian image pyramid, uses the FAST-10 corner detection algorithm to extract feature information in the current image, and uses block matching to establish data association between corner features . On this basis, according to the current projection error, a pose estimation model is established to realize the precise positioning of the camera, and the current 3D point cloud map is generated by combining feature matching information and triangulation algorithm. The specific process of camera pose estimation is as follows:

（1）稀疏地图的初始化(1) Initialization of sparse map

PTAM算法利用标准立体相机算法模型建立当前环境初始化地图，并在此基础上结合新增加关键帧不断更新三维地图。在地图的初始化过程中，通过人为选择两个独立关键帧，利用图像中FAST角点匹配关系，采用基于随机采样一致性(Random Sample Consensus,RANSAC)的五点法实现上述关键帧间重要矩阵F估计，并计算当前特征点处的三维坐标，同时，结合RANSAC算法选取适当空间点建立当前一致性平面，以确定全局世界坐标系，实现地图的初始化。The PTAM algorithm uses the standard stereo camera algorithm model to establish the current environment initialization map, and on this basis, it combines the newly added key frames to continuously update the three-dimensional map. In the initialization process of the map, by artificially selecting two independent key frames, using the FAST corner point matching relationship in the image, and using the five-point method based on Random Sample Consensus (RANSAC) to realize the important matrix F between key frames. Estimate and calculate the three-dimensional coordinates of the current feature points, and at the same time, combine the RANSAC algorithm to select appropriate spatial points to establish the current consistent plane to determine the global world coordinate system and realize the initialization of the map.

（2）相机位姿估计(2) Camera pose estimation

系统利用相机获取当前环境纹理信息，并构建四层高斯图像金字塔，运用FAST-10角点检测算法提取当前图像中特征信息，采用块匹配的方式建立角点特征间数据关联。在此基础上，根据当前投影误差，建立位姿估计模型，其数学描述如下：The system uses the camera to obtain the current environmental texture information, and builds a four-layer Gaussian image pyramid, uses the FAST-10 corner detection algorithm to extract feature information in the current image, and uses block matching to establish data association between corner features. On this basis, according to the current projection error, a pose estimation model is established, and its mathematical description is as follows:

$ξ ξ = = \underset{ξ ξ}{arg arg min min} ΣObj Σ Obj ((\frac{| | {e e}_{j j} | |}{{σ σ}_{j j}},, {σ σ}_{T T}))$

${e e}_{j j} = = (\begin{matrix} {u u}_{i i} \\ {v v}_{i i} \end{matrix}) - - Kπ Kπ ((exp exp ((\overset{^^}{ξ ξ})) p p))$

式中，e_j是投影误差，∑Obj(·,σ_T)为Tukey双权目标函数，σ_T为特征点的匹配标准差的无偏估计值，ξ为当前位姿6元组表示，

为由ξ组成的反对称矩阵。In the formula, e _j is the projection error, ∑Obj(·,σ _T ) is the Tukey dual-weight objective function, σ _T is the unbiased estimate of the matching standard deviation of the feature point, ξ is the current pose 6-tuple representation,

is an antisymmetric matrix composed of ξ.

根据上述位姿估计模型，选取位于图像金字塔顶层的50个特征匹配点，实现对相机的初始化位姿估计。更进一步，该算法结合相机初始位姿，采用极线收索的方式，建立图像金字塔中角点特征亚像素精度匹配关系，并将上述匹配对带入位姿估计模型，实现相机的精确重定位。According to the above pose estimation model, 50 feature matching points located at the top of the image pyramid are selected to realize the initial pose estimation of the camera. Furthermore, the algorithm combines the initial pose of the camera and adopts the method of epipolar retrieval to establish the sub-pixel precision matching relationship of the corner features in the image pyramid, and brings the above matching pairs into the pose estimation model to realize the precise repositioning of the camera .

（3）相机位姿优化(3) Camera pose optimization

系统经初始化后，地图创建线程将等待新的关键帧进入。若相机与当前关键帧间图像帧数超出阈值条件，且相机跟踪效果最佳时，将自动执行添加关键帧过程。此时，系统将会对新增加关键帧中所有FAST角点进行Shi-Tomas评估，以获取当前具有显著特征的角点信息，并选取与之最近的关键帧利用极线收索与块匹配方法建立特征点映射关系，结合位姿估计模型实现相机精确重定位，同时将匹配点投影到空间，生成当前全局环境三维地图。After the system is initialized, the map creation thread will wait for new keyframes to come in. If the number of image frames between the camera and the current key frame exceeds the threshold condition, and the camera tracking effect is the best, the process of adding key frames will be automatically performed. At this time, the system will perform Shi-Tomas evaluation on all FAST corner points in the newly added keyframes to obtain the current corner point information with salient features, and select the nearest keyframe to use epipolar search and block matching methods Establish the mapping relationship of feature points, combine the pose estimation model to achieve precise camera repositioning, and project the matching points into space to generate a 3D map of the current global environment.

为了实现全局地图的维护，在地图创建线程等待新关键帧进入的过程中，系统将利用局部与全局的Levenberg-Marquardt集束调整算法实现当前地图的一致性优化。该集束调整算法的数学描述为：In order to maintain the global map, the system will use the local and global Levenberg-Marquardt bundle adjustment algorithm to achieve the consistency optimization of the current map while the map creation thread is waiting for new keyframes to enter. The mathematical description of the cluster adjustment algorithm is:

${{{{= = {ξ ξ}_{11} . . . . . . {ξ ξ}_{N N}}},, {{{p p}_{11} . . . . . . {p p}_{M m}}}}} = = \underset{{{{{μ μ}},, {{p p}}}}}{arg arg min min} {Σ Σ}_{i i = = 11}^{N N} \underset{j j &Element; &Element; {s the s}_{i i}}{Σ Σ} Obj Obj ((\frac{| | {e e}_{ji the ji} | |}{{σ σ}_{ji the ji}},, {σ σ}_{T T}))$

式中，σ_ji为在第i个关键帧中，FAST特征点的匹配标准差的无偏估计，ξ_i表示第i个关键帧位姿的6元组表示，p_i为全局地图中的点。where σ _ji is the unbiased estimate of the matching standard deviation of the FAST feature points in the ith keyframe, _ξi is the 6-tuple representation of the pose of the ith keyframe, p _i is the point in the global map .

步骤四，建立基于变分机制的深度图估计模型，并求解该模型。Step 4: Establish a variational mechanism-based depth map estimation model and solve the model.

在PTAM精确位姿估计前提下，本发明基于多视重建方法，利用变分机制建立深度求解模型。该方法基于光照不变性与深度图平滑性假设，建立L1型数据惩罚项与变分规则项，该模型由在光照不变性假设的前提下建立数据惩罚项，并利用数据惩罚项保证当前深度图的平滑性，其数学模型如下：On the premise of PTAM accurate pose estimation, the present invention is based on a multi-view reconstruction method, and uses a variational mechanism to establish a depth solution model. Based on the assumption of illumination invariance and depth map smoothness, this method establishes L1 type data penalty items and variation rule items. The model establishes data penalty items under the premise of illumination invariance assumption, and uses data penalty items to ensure that the current depth map The smoothness of , its mathematical model is as follows:

E_d＝∫_Ω(E_data+λE_reg)dxE _d ＝∫ _Ω (E _data +λE _reg )dx

式中，λ为数据惩罚项E_data与变分规则项E_reg间的权重系数，

为深度图取值范围。In the formula, λ is the weight coefficient between the data penalty item E _data and the variation rule item E _reg ,

Value range for the depth map.

通过选取当前关键帧为深度图估计算法的参考帧I_r，利用其相邻像序列I＝{I₁,I₂,...,I_n}，结合投影模型建立数据惩罚项E_data，其数学描述为：By selecting the current key frame as the reference frame I _r of the depth map estimation algorithm, using its adjacent image sequence I={I ₁ ,I ₂ ,..., _In }, combined with the projection model to establish the data penalty item E _data , its The mathematical description is:

${E E.}_{data data} = = \frac{11}{| | I I ((r r)) | |} \underset{{I I}_{i i} &Element; &Element; I I}{Σ Σ} | | {I I}_{r r} ((x x)) - - {I I}_{i i} (({x x}^{' '})) | |$

式中，|I(r)|为当前图像序列中与参考帧具有重合信息的图像帧数量，x′为在当深度d下参考帧x处在I_i处的投影坐标，即：In the formula, |I(r)| is the number of image frames in the current image sequence that have overlapping information with the reference frame, and x' is the projection coordinate of the reference frame x at I _i at the current depth d, that is:

${x x}^{' '} = = {π π}^{- - 11} (({KT KT}_{r r}^{i i} π π ((x x,, d d))))$

在深度图平滑性假设前提下，为了确保在图像中边界处的不连续性，引入加权Huber算子构建变分规则项，其数学描述为：Under the assumption of smoothness of the depth map, in order to ensure the discontinuity at the boundary in the image, a weighted Huber operator is introduced to construct a variational rule item, whose mathematical description is:

E_reg＝g(u)||▽d(u)||_α E _reg ＝g(u)||▽d(u)|| _α

式中，▽d为深度图的梯度，g(u)为像素梯度权重系数，且g(u)＝exp(-a||▽I_r(u)||)In the formula, ▽d is the gradient of the depth map, g(u) is the weight coefficient of the pixel gradient, and g(u)=exp(-a||▽I _r (u)||)

Huber算子||x||_α的数学描述为：The mathematical description of the Huber operator ||x|| _α is:

${| | | | x x | | | |}_{α α} = = \{\begin{matrix} \frac{{| | | | x x | | | |}^{22}}{22 α α},, | | | | x x | | | | \leq \leq α α \\ | | | | x x | | | | - - \frac{α α}{22},, others others \end{matrix}$

式中，α为常量。In the formula, α is a constant.

根据Legendre-Fenchel变换，能量函数可表示为： $g {| | &dtri; d | |}_{α} = < g &dtri; d, q > - δ (q) - \frac{α}{2} {| | q | |}^{2}$ According to the Legendre-Fenchel transformation, the energy function can be expressed as: $g {| | &dtri; d | |}_{α} = < g &dtri; d, q > - δ (q) - \frac{α}{2} {| | q | |}^{2}$

式中， $δ (q) = \{\begin{matrix} \frac{α}{2} & α < | | q | | \leq 1 \\ \infty & others \end{matrix}$ In the formula, $δ (q) = \{\begin{matrix} \frac{α}{2} & α < | | q | | \leq 1 \\ \infty & others \end{matrix}$

上述Huber算子的引入为三维重建过程提供了光滑性保证，同时也为确保深度图中存在非连续边界，提高了三维地图创建质量。The introduction of the above-mentioned Huber operator provides a smoothness guarantee for the 3D reconstruction process, and also ensures the presence of discontinuous boundaries in the depth map, improving the quality of 3D map creation.

针对上述数学模型求解复杂度高、计算量大的问题，引入辅助变量建立凸优化模型，采用交替下降法实现对上述模型的优化，其具体过程如下：Aiming at solving the problems of high complexity and large amount of calculation in the above mathematical model, an auxiliary variable is introduced to establish a convex optimization model, and the optimization of the above model is realized by using the alternating descent method. The specific process is as follows:

(1)固定h，求解：(1) Fix h, solve:

$\underset{q q}{arg arg max max} {{\underset{d d}{arg arg min min} {E E.}_{d d,, q q}}}$

${E E.}_{d d,, q q} = = {&Integral; &Integral;}_{Ω Ω} ((< < g g &dtri; &dtri; d d,, q q > > + + \frac{11}{22 θ θ} {((d d - - h h))}^{22} - - δ δ ((q q)) - - \frac{α α}{22} {| | | | q q | | | |}^{22})) dx dx$

式中，θ为二次项常系数，g为变分规则项中梯度权重系数。In the formula, θ is the constant coefficient of the quadratic term, and g is the gradient weight coefficient in the variation rule term.

根据拉格朗日极值法，上述能量函数达到极值的条件为：According to the Lagrangian extreme value method, the condition for the above energy function to reach the extreme value is:

$\frac{{&PartialD; &PartialD; E E.}_{d d,, q q}}{&PartialD; &PartialD; q q} = = g g &dtri; &dtri; d d - - αq αq = = 00$

$\frac{{&PartialD; &PartialD; E E.}_{d d,, q q}}{&PartialD; &PartialD; d d} = = g g div div q q + + \frac{11}{θ θ} ((d d - - h h)) = = 00$

式中，divq为q的散度。In the formula, divq is the divergence of q.

结合偏导数离散化描述，上述极值条件可表示为：Combined with the discretization description of partial derivatives, the above extreme value condition can be expressed as:

$\frac{{q q}^{n no + + 11} - - {q q}^{n no}}{{ϵ ϵ}_{q q}} = = g g &dtri; &dtri; d d - - {αq αq}^{n no + + 11}$

$\frac{{d d}^{n no + + 11} - - {d d}^{n no}}{{ϵ ϵ}_{d d}} = = g g div div p p + + \frac{11}{θ θ} (({d d}^{n no + + 11} - - h h))$

此时可采用原始对偶算法实现能量函数的迭代优化，即：At this time, the original dual algorithm can be used to realize the iterative optimization of the energy function, namely:

${p p}^{n no + + 11} = = \frac{(({p p}^{n no} + + {ϵ ϵ}_{q q} g g &dtri; &dtri; {d d}^{n no})) / / ((11 + + {ϵ ϵ}_{q q} α α))}{max max ((11,, (({p p}^{n no} + + {ϵ ϵ}_{q q} g g &dtri; &dtri; {d d}^{n no})) / / ((11 + + {ϵ ϵ}_{q q} α α))))}$

${d d}^{n no + + 11} = = \frac{{d d}^{n no} + + {ϵ ϵ}_{d d} (({g g div div q q}^{n no + + 11} + + {h h}^{n no} / / θ θ))}{((11 + + {ϵ ϵ}_{d d} / / θ θ))}$

式中，ε_q、ε_d为常数，分别表示最大化与最小化梯度描述系数。In the formula, ε _q and ε _d are constants, representing the maximum and minimum gradient description coefficients respectively.

(2)固定d，求解：(2) Fix d, solve:

$\underset{h h}{arg arg min min} {E E.}_{h h}$

${E E.}_{h h} = = {&Integral; &Integral;}_{Ω Ω} ((\frac{θ θ}{22} {((d d - - h h))}^{22} + + \frac{λ λ}{| | I I ((r r)) | |} {Σ Σ}_{i i = = 00}^{n no} | | {I I}_{i i} ((x x)) - - {I I}_{ref ref} ((x x,, h h)) | |)) dx dx$

在上述能量函数求解过程中，为了有效减少算法的复杂度，同时保证重建过程中的部分细节信息。本发明将深度取值范围[d_min,d_max]划分为S个采样平面，采用穷举的方式获取当前能量函数的最优解。其中步长的选择为：In the process of solving the above energy function, in order to effectively reduce the complexity of the algorithm, and at the same time ensure part of the detailed information in the reconstruction process. The present invention divides the depth value range [d _min , d _max ] into S sampling planes, and obtains the optimal solution of the current energy function in an exhaustive manner. The choice of step size is:

${d d}_{inc inc}^{k k} = = \frac{{Sd SD}_{min min} {d d}_{max max}}{((S S - - k k)) {d d}_{min min} + + {d d}_{max max}}$

式中，为k与k-1采样平面间隔。In the formula, Sampling plane intervals for k and k-1.

步骤五，建立关键帧选取机制，实现三维场景的更新。Step five, establishing a key frame selection mechanism to realize updating of the 3D scene.

考虑系统冗余信息的消除，为了提高重建结果的清晰度以及实时性，减少系统在计算负担，本发明只在关键帧处实现对三维场景的估计，并更新和维护所生成的三维场景。当新增一帧KeyFrame数据后，根据式

将当前新增KeyFrame数据转换到世界坐标系中，完成场景数据的更新。Considering the elimination of redundant information in the system, in order to improve the clarity and real-time performance of the reconstruction results and reduce the computational burden of the system, the present invention only realizes the estimation of the 3D scene at key frames, and updates and maintains the generated 3D scene. After adding a frame of KeyFrame data, according to the formula

Convert the current newly added KeyFrame data to the world coordinate system to complete the update of the scene data.

利用深度模型中数据惩罚项，建立当前帧与关键帧间的信息重合程度评估函数，即：Using the data penalty item in the depth model, the evaluation function of the information coincidence degree between the current frame and the key frame is established, namely:

$N N = = \underset{x x &Element; &Element; {R R}^{22}}{Σ Σ} c c ((x x))$

$c c ((x x)) = = \{\begin{matrix} 11,, & | | {I I}_{r r} ((x x)) - - {I I}_{i i} (({x x}^{' '})) | | < < ζ ζ \\ 00,, & others others \end{matrix}$

式中，ζ为常数。In the formula, ζ is a constant.

若此时N小于图像大小的0.7时，即确定当前帧为新关键帧。If N is smaller than 0.7 of the image size at this time, it is determined that the current frame is a new key frame.

本发明的有益效果是：本发明采用RGB相机获取环境数据。针对利用高精度单目定位算法，提出一种基于变分机制的深度图生成方法，实现了大范围的快速室内三维场景重建，有效地解决了三维重建算法成本与实时性问题。The beneficial effect of the present invention is that: the present invention adopts RGB camera to obtain environmental data. For the use of high-precision monocular positioning algorithm, a depth map generation method based on variational mechanism is proposed, which realizes a large-scale rapid indoor 3D scene reconstruction, and effectively solves the cost and real-time problems of 3D reconstruction algorithms.

附图说明Description of drawings

图1为基于变分模型的室内三维场景重建方法流程图；Fig. 1 is the flowchart of the indoor 3D scene reconstruction method based on the variational model;

图2为相机坐标系示意图；Figure 2 is a schematic diagram of the camera coordinate system;

图3为本发明应用实例的三维重建实验结果。Fig. 3 is the experimental result of three-dimensional reconstruction of the application example of the present invention.

具体实施方式Detailed ways

图1是基于变分模型的室内三维场景重建方法流程图，包括以下步骤：Figure 1 is a flowchart of a method for indoor 3D scene reconstruction based on a variational model, including the following steps:

步骤二，建立相机位姿描述与相机投影模型。Step 2, establish camera pose description and camera projection model.

下面给出本发明的一个应用实例。An application example of the present invention is given below.

本实例采用的RGB相机为Point Grey Flea2，图像辨率为640×480，最高帧频为30fps，水平视场角为65°，焦距大约为3.5mm。所使用的PC机配备有GTS450GPU和i5四核CPU。The RGB camera used in this example is Point Gray Flea2, the image resolution is 640×480, the maximum frame rate is 30fps, the horizontal field of view is 65°, and the focal length is about 3.5mm. The PC used is equipped with GTS450GPU and i5 quad-core CPU.

在实验过程中，通过彩色相机获取环境深度信息，结合相机位姿估计算法实现对自身精确定位。当进入关键帧后，选择关键帧周围20帧图像作为本文深度估计算法的输入。在深度估计算法执行过程中，令d⁰＝h⁰且q⁰＝0，计算

以获取当前深度图的初始化输入，并迭代优化E_d,q与E_h直到收敛。同时，该在算法迭代过程中不断减小θ值，增加二次函数在算法执行过程中的权重，有效提高了算法收敛速度。最终实验结果如图3所示，实验表明该方法能有效实现环境的稠密三维重建，并进一步验证了该方法的可行性。During the experiment, the depth information of the environment is obtained through the color camera, and the camera pose estimation algorithm is used to realize the precise positioning of itself. After entering the key frame, select 20 frames of images around the key frame as the input of the depth estimation algorithm in this paper. During the execution of the depth estimation algorithm, let d ⁰ =h ⁰ and q ⁰ =0, calculate

To obtain the initialization input of the current depth map, and iteratively optimize E _{d, q} and E _h until convergence. At the same time, the value of θ is continuously reduced during the algorithm iteration process, and the weight of the quadratic function in the algorithm execution process is increased, which effectively improves the algorithm convergence speed. The final experimental results are shown in Figure 3. The experiment shows that the method can effectively realize the dense 3D reconstruction of the environment, and further verifies the feasibility of the method.

Claims

1. method based on the indoor environment three-dimensional reconstruction of variation mechanism is characterized in that may further comprise the steps:

Step 1 is obtained the calibrating parameters of camera, and sets up the distortion correction model;

The detailed process of camera calibration is:

(1) prints a chessboard template;

(2) from a plurality of angle shot chessboards, should allow chessboard take screen, and each angle that guarantees chessboard be taken 6 template picture altogether all in screen as far as possible;

(3) detect unique point in the image, i.e. each black point of crossing of chessboard;

(4) inner parameter of asking for, method is as follows:

RGB camera calibration parameter is mainly the camera confidential reference items, and the confidential reference items matrix K of camera is:

K = [\begin{matrix} f_{u} & 0 & u_{0} \\ 0 & f_{v} & v_{0} \\ 0 & 0 & 1 \end{matrix}]

In the formula, u, v are the camera plane coordinate axis, (u ₀, v ₀) be that camera is as planar central coordinate, (f _u, f _v) be the focal length of camera;

According to calibrating parameters, the mapping relations of RGB image mid point and three dimensions point are as follows: RGB image mid point p=(u, v) the coordinate P under camera coordinates system _3D=(x, y z) are expressed as:

\{\begin{matrix} x = (u - u_{0}) * z / f_{u} \\ y = (v - v_{0}) * z / f_{v} \\ z = d \end{matrix}

In the formula, d represents the depth value of depth image mid point p;

Camera coordinates system is y axle positive dirction downwards, is forward z axle positive dirction, is to the right the x positive dirction; The initial point position of camera is set at the world coordinate system initial point, and the X of world coordinate system, Y, Z direction are identical with the definition of camera;

FOV camera correction model is:

u_{d} = [\begin{matrix} u_{0} \\ v_{0} \end{matrix}] + [\begin{matrix} f_{u} & 0 \\ 0 & f_{v} \end{matrix}] \frac{r_{d}}{r_{u}} x_{u}

r_{d} = \frac{1}{ω} \arctan ({2 r}_{u} \tan \frac{ω}{2})

r_{u} = \frac{\tan (r_{d} ω)}{2 \tan \frac{ω}{2}}

In the formula, x _uBe the pixel coordinate of z=1 face, u _dBe pixel coordinate in the original image, ω is FOV camera distortion factor;

Step 2 is set up the camera pose and is described and the camera projection model, and direction is as follows:

Under the world coordinate system of having set up, the camera pose can be expressed as matrix:

T_{cw} = [\begin{matrix} R_{cw} & t_{cw} \\ 0 & 1 \end{matrix}]

In the formula, cw represents that being tied to current camera coordinates from world coordinates is T _Cw∈ SE (3), SE (3) are the rotation translation transformation space of rigid body; T _CwCan be by following hexa-atomic group of μ=(μ ₁, μ ₂, μ ₃, μ ₄, μ ₅, μ ₆) expression, that is:

T_{cw} = \exp (\hat{μ})

\hat{μ} = [\begin{matrix} 0 & μ_{6} & - μ_{5} & μ_{1} \\ μ_{6} & 0 & μ_{4} & μ_{2} \\ μ_{5} & - μ_{4} & 0 & μ_{3} \\ 0 & 0 & 0 & 0 \end{matrix}]

In the formula, μ ₁, μ ₂, μ ₃Be respectively the translational movement of Kinect under global coordinate system, μ ₄, μ ₅, μ ₆The rotation amount of coordinate axis under the expression local coordinate system;

The pose T of camera _CwSet up spatial point cloud coordinate p under the current coordinate system _cTo world coordinates p _wTransformation relation, that is:

p _c＝T _cwp _w

Under current mark system, the projection to the z=1 plane of three dimensions point cloud is defined as:

π(p)＝(xz,yz) ^T

In the formula, p ∈ R ³Be the three dimensions point, x, y, z are the coordinate figure of this point; According to current coordinate points depth value d, utilize reverse sciagraphy to determine current space three-dimensional point coordinate p, its coordinate relation can be expressed as:

π ^-1(u,d)＝dK ^-1u

Step 3 is utilized based on the monocular SLAM algorithm of SFM and is realized the estimation of camera pose;

Step 4 is set up the depth map estimation model based on variation mechanism, and is found the solution this model;

Step 5 is set up the key frame selection mechanism, realizes the renewal of three-dimensional scenic, and method is as follows:

Realize the estimation to three-dimensional scenic at the key frame place, and upgrade and safeguard the three-dimensional scenic that generates; After newly-increased frame KeyFrame data, according to formula

Current newly-increased KeyFrame data-switching in world coordinate system, is finished the renewal of contextual data;

Utilize data penalty term in the depth model, set up present frame and overlap the degree valuation functions with information between key frame, that is:

N = \underset{x &Element; R^{2}}{Σ} c (x)

c (x) = [\begin{matrix} 1, & | I_{r} (x) - I_{i} (x^{'}) | < ζ \\ 0, & others \end{matrix}]

In the formula, ζ is constant;

If N was less than 0.7 o'clock of the image size at this moment, determine that namely present frame is new key frame.

2. the method for a kind of indoor environment three-dimensional reconstruction based on variation mechanism according to claim 1 is characterized in that, the step 3 utilization realizes that based on the monocular SLAM algorithm of SFM camera pose estimation approach is further comprising the steps of:

(1) initialization of sparse map

The PTAM algorithm utilizes standard stereoscopic camera algorithm model to set up current environment initialization map, and brings in constant renewal in three-dimensional map in conjunction with increasing key frame newly on this basis; In the initialization procedure of map, by two independent key frames of artificial selection, utilize FAST corners Matching relation in the image, employing realizes the estimation of the important matrix F between above-mentioned key frame based on the conforming five-spot of stochastic sampling, and calculate the three-dimensional coordinate at current unique point place, simultaneously, set up current consistance plane in conjunction with the suitable spatial point of RANSAC algorithm picks, to determine overall world coordinate system, realize the initialization of map;

(2) the camera pose is estimated

System utilizes camera to obtain the current environment texture information, and makes up this image pyramid of four floor heights, uses the FAST-10 Corner Detection Algorithm to extract characteristic information in the present image, and the mode of employing piece coupling is set up the data association between the angle point feature; On this basis, according to current projection error, set up the pose estimation model, its mathematical description is as follows:

ξ = \underset{ξ}{\arg \min} ΣObj (\frac{| e_{j} |}{σ_{j}}, σ_{T})

e_{j} = (\begin{matrix} u_{i} \\ v_{i} \end{matrix}) - Kπ (\exp (ξ) p)

In the formula, e _jBe projection error, Σ Obj (, σ _T) be the two power of Tukey objective function function, σ _TBe the unbiased estimator of the match-on criterion difference of unique point, ξ is current pose 6 element group representations,

Be the antisymmetric matrix of being formed by ξ;

According to above-mentioned pose estimation model, choose 50 characteristic matching points that are positioned at the image pyramid top layer, realize the initialization pose of camera is estimated; Further, the initial pose of this algorithm combining camera adopts polar curve to receive the mode of rope, sets up angle point feature sub-pixel precision matching relationship in the image pyramid, and with above-mentioned coupling to bringing the pose estimation model into, realize the accurate reorientation of camera;

(3) the camera pose is optimized

System is after initialization, and the key frame that the map building thread waits is new enters; If number of image frames exceeds threshold condition between camera and current key frame, and the camera tracking effect will automatically perform the key frame process of adding when best; At this moment, system will carry out the Shi-Tomas assessment to all FAST angle points in the key frame that increases newly, to obtain current angle point characteristic information with notable feature, and choose nearest with it key frame and utilize polar curve receipts rope and block matching method to set up the unique point mapping relations, realize the accurate reorientation of camera in conjunction with the pose estimation model, simultaneously match point is projected to the space, generate current global context three-dimensional map;

In order to realize the maintenance of global map, in the process that the new key frame of map building thread waits enters, the local Levenberg-Marquardt boundling adjustment algorithm with the overall situation of system's utilization realizes the global coherency optimization of current map; The mathematical description of this boundling adjustment algorithm is:

{{ξ_{2} . . ξ_{N}}, {p_{1} . . p_{M}}} = \underset{{{μ}, {p}}}{\arg \min} Σ_{i = 1}^{N} \underset{j &Element; S_{i}}{Σ} Obj (\frac{| e_{ji} |}{σ_{ji}}, σ_{T})

In the formula, σ _JiFor in i key frame, the nothing of the match-on criterion difference of FAST unique point is estimated ξ partially _i6 element group representations of representing i key frame pose, p _iBe the point in the global map.

3. the method for a kind of indoor environment three-dimensional reconstruction based on variation mechanism according to claim 1 is characterized in that, step 4 is set up and find the solution based on the method for the depth map estimation model of variation mechanism as follows:

Based on the depth map estimation model of variation mechanism, under the prerequisite of illumination unchangeability hypothesis, to set up the data penalty term, and utilize the data penalty term to guarantee the flatness of current depth map, its mathematical model is as follows:

E _d＝∫ _Ω(E _data+λE _reg)dx

In the formula, λ is data penalty term E _DataWith variation regularization term E _RegBetween weight coefficient,

Be the depth map span;

By choosing the reference frame I that current key frame is the depth map algorithm for estimating _r, utilize its adjacent picture sequence I={I ₁, I ₂..., I _n, set up data penalty term E in conjunction with projection model _Data, its mathematical description is:

E_{data} = \frac{1}{| I (r) |} \underset{I_{i} &Element; I}{Σ} | I_{r} (x) - I_{i} (x^{'}) |

In the formula, | I (r) | for have the image frames numbers of the information of coincidence in the present image sequence with reference frame, x ' is for being in I at reference frame x under depth d _iThe projection coordinate at place, that is:

x^{'} = π^{- 1} ({KT}_{r}^{i} π (x, d))

Under depth map smoothness assumption prerequisite, in order to ensure the uncontinuity of boundary in image, to introduce Weighted H uber operator and make up the variation regularization term, its mathematical description is:

E _reg＝g(u)||▽d(u)|| _α

In the formula, ▽ d is the gradient of depth map, and g (u) is the pixel gradient weight coefficient, g (u)=exp (a|| ▽ I _r(u) ||)

The Huber operator || x|| _αMathematical description be:

{| | x | |}_{α} = \{\begin{matrix} \frac{{| | x | |}^{2}}{2 α}, | | x | | \leq α \\ | | x | | - \frac{α}{2}, others \end{matrix}

In the formula, α is constant;

According to the Legendre-Fenchel conversion, energy function is transformed to:

g {| | &dtri; d | |}_{α} = < g &dtri; d, q > - δ (q) - \frac{α}{2} {| | q | |}^{2}

In the formula,

δ (q) = \{\begin{matrix} \frac{α}{2} & α < | | q | | \leq 1 \\ \infty & other \end{matrix}

In view of above-mentioned mathematical model is found the solution the complexity height, calculated amount is big, introduce auxiliary variable and set up protruding optimization model, adopting alternately, descent method realizes that to above-mentioned Model Optimization detailed process is as follows:

(1) fixing h, find the solution:

\underset{q}{\arg \max} {\underset{d}{\arg \min} E_{d, q}}

E_{d, q} = {&Integral;}_{Ω} (< g &dtri; d, q > + \frac{1}{2 θ} {(d - h)}^{2} - δ (q) - \frac{α}{2} {| | q | |}^{2}) dx

In the formula, g is gradient weight coefficient in the variation regularization term, and θ is the quadratic term constant coefficient;

According to Lagrangian extremum method, the condition that above-mentioned energy function reaches extreme value is:

\frac{&PartialD; E_{d, q}}{&PartialD; q} = g &dtri; d - αq = 0

\frac{{&PartialD; E}_{d, q}}{&PartialD; d} = g div q + \frac{1}{θ} (d - h) = 0

In the formula, divq is the divergence of q;

Describe in conjunction with the partial derivative discretize, above-mentioned extremum conditions can be expressed as:

\frac{q^{n + 1} - q^{n}}{ϵ_{q}} = g &dtri; d - {αq}^{n + 1}

\frac{d^{n + 1} - d^{n}}{ϵ_{d}} = g div p + \frac{1}{θ} (d^{n + 1} - h)

Adopt primal dual algorithm to realize the iteration optimization of energy function, that is:

p^{n + 1} = \frac{(p^{n} + ϵ_{q} g &dtri; d^{n}) / (1 + ϵ_{q} α)}{\max (1, (p^{n} + ϵ_{q} g &dtri; d^{n}) / (1 + ϵ_{q} α))}

d^{n + 1} = \frac{d^{n} + ϵ_{d} ({g div q}^{n + 1} + h^{n} / θ)}{(1 + ϵ_{d} / θ)}

In the formula, ε _q, ε _dBe constant, expression maximizes and minimizes gradient and describes coefficient respectively;

(2) fixing d, find the solution:

\underset{h}{\arg \min E_{h}}

E_{h} = {&Integral;}_{Ω} (\frac{θ}{2} {(d - h)}^{2} + \frac{λ}{| I (r) |} Σ_{i = 0}^{n} | I_{i} (x) - I_{ref} (x, h) |) dx

In above-mentioned energy function solution procedure, in order effectively to reduce the complexity of algorithm, guarantee the part detailed information in the process of reconstruction simultaneously, with degree of depth span [d _Min, d _Max] be divided into S sample plane, adopt exhaustive mode to obtain the optimum solution of current energy function; Being chosen as of step-length wherein:

d_{inc}^{k} = \frac{{Sd}_{\min} d_{\max}}{(S - k) d_{\min} + d_{\max}}

In the formula,

Be k and k-1 sample plane interval.