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

Abstract

The invention belongs to crossing field of computer vision and intelligent robots and discloses a variational mechanism-based large-area indoor scene reconstruction method. The method comprises the following steps: step 1, acquiring calibration parameters of a camera, and building an aberration correcting model; step 2, building a camera position and gesture depiction and camera projection model; step 3, utilizing an SFM-based monocular SFM (Space Frequency Modulation) algorithm to realize camera position and gesture estimation; step 4, building a variational mechanism-based depth map estimation model, and performing solving on the model; and step 5, building a key frame selection mechanism to realize three-dimensional scene renewal. According to the invention, an RGB (Red Green Blue) camera is adopted to acquire environmental data, and a variational mechanism-based depth map generation method is proposed through utilizing a high-precision monocular positioning algorithm, so that quick large-area indoor three-dimensional scene reconstruction is realized, and problems of three-dimensional reconstruction algorithm cost and real-time performance are effectively solved.

Description

Indoor environment three-dimensional rebuilding method based on variation mechanism

Technical field

The invention belongs to the crossing domain of computer vision and intelligent robot, relate to a kind of indoor environment three-dimensional reconstruction technology, relate in particular to a kind of method for reconstructing of the indoor scene on a large scale based on variation mechanism.

Technical background

(modeling of environment 3 D stereo progressively becomes this area research focus, causes numerous scholars' concern for Simultaneous Localization And Mapping, SLAM) deepening continuously of research with map building along with locating simultaneously.G.Klein equals at first to propose to locate simultaneously in augmented reality (AR) field in 2007 that (Parallel Tracking and Mapping, concept PTAM) is to solve environment Real-time modeling set problem with map building.PTAM follows the tracks of video camera with the map generation and is divided into two separate threads, when utilizing the FastCorner method to upgrade detected characteristics point, (Bundle Adjustment BA), constantly realizes the renewal of camera pose and three-dimensional feature point map to adopt optimum part and overall light beam method of adjustment.This method has been set up the environment three-dimensional map based on sparse some cloud, but this map lacks the three-dimensional description directly perceived to environment.People such as Pollefeys have realized the three-dimensional reconstruction of large-scale outdoor scene by Multi-sensor Fusion.But there is the high complexity that calculates in this method and to shortcomings such as noise sensitivities.In the progress that some trial property have also been arranged aspect real-time follow-up and the dense environmental model reconstruct, still only be confined to the reconstruct of some simple objects, and can only under the particular constraints condition, can obtain higher precision at present.People such as Richard A.Newcombe, utilization is based on SFM(Structure from Moving) the SLAM algorithm obtain space sparse features point cloud, adopt multiple dimensioned radially basic interpolation, use Implicit Surface Polygonization method in graph image, structure three dimensions initialization grid map, and in conjunction with scene flows constraint and high precision TV-L1 optical flow algorithm renewal grid vertex coordinate, to reach the purpose of approaching real scene.This algorithm can obtain high-precision environmental model, but because its algorithm complex is higher, under two graphic hardware processors (GPU) acceleration situation, handles the time that a two field picture still need spend several seconds.

Summary of the invention

At the above-mentioned problems in the prior art, the invention provides a kind of quick three-dimensional reconstructing method based on variation mechanism, to be implemented in the three-dimensional modeling under the indoor complex environment.This method has reduced required deal with data amount when guaranteeing environmental information, can realize quick indoor 3 D scene rebuilding on a large scale.Solve three-dimensional reconstruction algorithm cost and real-time problem effectively, improved the reconstruction precision.

The technical solution used in the present invention is as follows:

Utilize the PTAM algorithm to estimate means as the camera pose, and choose suitable image sequence structure at the key frame place based on the depth map estimated energy function of variation pattern, use primal dual algorithm to optimize above-mentioned energy function, be implemented in obtaining of current key frame place environment depth map.Because this algorithm utilizes contiguous frames information structuring energy function, and effectively utilized the relevance between the certain viewing angles coordinate system, and the translating camera perspective projection relation, make data item contain and look the imaging constraint more, reduced the computation complexity that algorithm model is found the solution.Under the unified calculation framework, the present invention utilizes graphics accelerator hardware to realize the parallel optimization of algorithm, has effectively improved the algorithm real-time.

A kind of method of the indoor environment three-dimensional reconstruction based on 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.

In computer vision is used, by the geometric model of camera imaging, effectively set up the mapping relations between the pixel and space three-dimensional point in the image.The geometric parameter that constitutes camera model must just can obtain with calculating by experiment, and the process of finding the solution above-mentioned parameter just is referred to as camera calibration.The demarcation of camera parameter in the present invention is unusual the key link, and the precision of calibrating parameters directly influences the accuracy of net result three-dimensional map.

The detailed process of camera calibration is:

(1) prints a chessboard template.The present invention adopts an A4 paper, chessboard be spaced apart 0.25cm.

(2) from a plurality of angle shot chessboards.During shooting, 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) ask for the inner parameter of camera, method is as follows:

RGB camera calibration parameter is mainly the camera confidential reference items.The confidential reference items matrix K of camera is:

$K=\left[\begin{array}{ccc}{f}_u& 0& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right]$

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:

$\left\{\begin{array}{c}x=(u-u_0)*z/f_u\\ y=(v-v_0)*z/f_v\\ z=d\end{array}\right.$

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

Camera coordinates system is downwards y axle positive dirction as shown in Figure 2 among the present invention, 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(Field of Viewer) the camera correction model is:

$u_d=\left[\begin{array}{c}{u}_0\\ v_0\end{array}\right]+\left[\begin{array}{cc}{f}_u& 0\\ 0& f_v\end{array}\right]\frac{r_d}{r_u}{x}_u$

$r_d=\frac{1}{\mathrm{\ω}}\arctan \left({2r}_u\tan \frac{\mathrm{\ω}}{2}\right)$

$r_u=\frac{\tan \left(r_d\mathrm{\ω}\right)}{2\tan \frac{\mathrm{\ω}}{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.

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

$T_{\mathrm{cw}}=\left[\begin{array}{cc}{R}_{\mathrm{cw}}& t_{\mathrm{cw}}\\ 0& 1\end{array}\right]$

In the formula, " cw " expression is tied to current camera coordinates from world coordinates and 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_{\mathrm{cw}}=\exp \left(\widehat{\mathrm{\μ}}\right)$

$\hat{u}=\left[\begin{array}{cccc}0& {\mathrm{\μ}}_6& -{\mathrm{\μ}}_5& {\mathrm{\μ}}_1\\ {\mathrm{\μ}}_6& 0& {\mathrm{\μ}}_4& {\mathrm{\μ}}_2\\ {\mathrm{\μ}}_5& -{\mathrm{\μ}}_4& 0& {\mathrm{\μ}}_3\\ 0& 0& 0& 0\end{array}\right]$

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.

At present, monocular vision SLAM algorithm mainly comprises the SLAM algorithm based on filtering and SFM (Structure from Moving).The present invention adopts the realization of PTAM algorithm to the location of camera.This algorithm is a kind of monocular vision SLAM method based on SFM, by being that camera is followed the tracks of and two of map buildings thread independently with system divides.In the camera track thread, system utilizes camera to obtain the current environment texture information, and make up this image pyramid of four floor heights, and use the FAST-10 Corner Detection Algorithm to extract characteristic information in the present image, 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 accurate location that the pose estimation model is realized camera, and generate current three-dimensional point cloud map in conjunction with characteristic matching information and triangulation algorithm.The detailed process that the camera pose is estimated is:

(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 based on the stochastic sampling consistance (Random Sample Consensus, five-spot RANSAC) realizes that important matrix F is estimated between above-mentioned key frame, and calculates 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, adopts the mode of piece coupling to set up 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:

$\mathrm{\ξ}=\underset{\mathrm{\ξ}}{\arg \min }\mathrm{\ΣObj}(\frac{\left|e_j\right|}{{\mathrm{\σ}}_j},{\mathrm{\σ}}_T)$

$e_j=\left(\begin{array}{c}{u}_i\\ v_i\end{array}\right)-\mathrm{K\π}\left(\exp \left(\widehat{\mathrm{\ξ}}\right)p\right)$

In the formula, e _jBe projection error, ∑ Obj (, σ _T) be the two power of Tukey objective 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 map building thread will wait for that new key frame 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 all FAST angle points carry out the Shi-Tomas assessment in the key frame to increasing newly, to obtain current angle point 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, system will utilize local Levenberg-Marquardt boundling adjustment algorithm with the overall situation to realize the consistance optimization of current map.The mathematical description of this boundling adjustment algorithm is:

$\left\{\right\{={\mathrm{\&xi;}}_1...{\mathrm{\&xi;}}_N\},\{{\mathrm{<figure-callout\; id="1"\; label="p"\; filenames="HDA00003178676100021.png"\; state="\{\{state\}\}">p</figure-callout>}}_1...p_M\left\}\right\}=\underset{\left\{\right\{\mathrm{\&mu;}\},\{p\left\}\right\}}{\arg \min }\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j\&Element;s_i}{\mathrm{\&Sigma;}}\mathrm{Obj}(\frac{\left|e_{\mathrm{ji}}\right|}{{\mathrm{\&sigma;}}_{\mathrm{ji}}},{\mathrm{\&sigma;}}_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.

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

Estimate to the present invention is based on many apparent weights construction method under the prerequisite at the accurate pose of PTAM, utilize variation mechanism to set up degree of depth solving model.This method is based on illumination unchangeability and depth map smoothness assumption, set up L1 type data penalty term and variation regularization term, this model is by setting up the data penalty term under the prerequisite of illumination unchangeability hypothesis, and utilizes the data penalty term to guarantee the flatness of current depth map, and 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_{\mathrm{data}}=\frac{1}{\left|I\left(r\right)\right|}\underset{I_i\&Element;I}{\mathrm{\&Sigma;}}|I_r\left(x\right)-I_i\left(x^{\&prime;}\right)|$

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^{\&prime;}={\mathrm{\&pi;}}^{-1}\left({\mathrm{KT}}_r^i\mathrm{\&pi;}(x,d)\right)$

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, and g (u)=exp (a|| ▽ I _r(u) ||)

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

${\left|\right|x\left|\right|}_{\mathrm{\&alpha;}}=\left\{\begin{array}{c}\frac{{\left|\right|x\left|\right|}^2}{2\mathrm{\&alpha;}},\left|\right|x\left|\right|\&le;\mathrm{\&alpha;}\\ \left|\right|x\left|\right|-\frac{\mathrm{\&alpha;}}{2},\mathrm{others}\end{array}\right.$

In the formula, α is constant.

According to the Legendre-Fenchel conversion, energy function can be expressed as: $g{\left|\right|\&dtri;d\left|\right|}_{\mathrm{\&alpha;}}=<g\&dtri;d,q>-\mathrm{\&delta;}\left(q\right)-\frac{\mathrm{\&alpha;}}{2}{\left|\right|q\left|\right|}^2$

In the formula, $\mathrm{\&delta;}\left(q\right)=\left\{\begin{array}{cc}\frac{\mathrm{\&alpha;}}{2}& \mathrm{\&alpha;}<\left|\right|q\left|\right|\&le;1\\ \&infin;& \mathrm{others}\end{array}\right.$

The three-dimensional reconstruction process that is introduced as of above-mentioned Huber operator provides the slickness assurance, also has discontinuous border in the depth map for guaranteeing simultaneously, has improved three-dimensional map and has created quality.

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

(1) fixing h, find the solution:

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

$E_{d,q}={\&Integral;}_{\mathrm{\&Omega;}}(<g\&dtri;d,q>+\frac{1}{2\mathrm{\&theta;}}{(d-h)}^2-\mathrm{\&delta;}\left(q\right)-\frac{\mathrm{\&alpha;}}{2}{\left|\right|q\left|\right|}^2)\mathrm{dx}$

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

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-\mathrm{\&alpha;q}=0$

$\frac{{\&PartialD;E}_{d,q}}{\&PartialD;d}=g\mathrm{div}q+\frac{1}{\mathrm{\&theta;}}(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}{{\mathrm{\&epsiv;}}_q}=g\&dtri;d-{\mathrm{\&alpha;q}}^{n+1}$

$\frac{d^{n+1}-d^n}{{\mathrm{\&epsiv;}}_d}=g\mathrm{div}p+\frac{1}{\mathrm{\&theta;}}(d^{n+1}-h)$

Can adopt this moment primal dual algorithm to realize the iteration optimization of energy function, that is:

$p^{n+1}=\frac{(p^n+{\mathrm{\&epsiv;}}_{q}g\&dtri;d^n)/(1+{\mathrm{\&epsiv;}}_q\mathrm{\&alpha;})}{\max (1,(p^n+{\mathrm{\&epsiv;}}_{q}g\&dtri;d^n)/(1+{\mathrm{\&epsiv;}}_q\mathrm{\&alpha;}))}$

$d^{n+1}=\frac{d^n+{\mathrm{\&epsiv;}}_d({g\mathrm{div}q}^{n+1}+h^n/\mathrm{\&theta;})}{(1+{\mathrm{\&epsiv;}}_d/\mathrm{\&theta;})}$

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;}_{\mathrm{\&Omega;}}(\frac{\mathrm{\&theta;}}{2}{(d-h)}^2+\frac{\mathrm{\&lambda;}}{\left|I\left(r\right)\right|}\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}|I_i\left(x\right)-I_{\mathrm{ref}}(x,h)\left|\right)\mathrm{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.The present invention is 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_{\mathrm{inc}}^k=\frac{{\mathrm{Sd}}_{\min }{d}_{\max }}{(S-k)d_{\min }+d_{\max }}$

In the formula, Be k and k-1 sample plane interval.

Step 5 is set up the key frame selection mechanism, realizes the renewal of three-dimensional scenic.

The elimination of taking into account system redundant information, for sharpness and the real-time that improves reconstructed results, the minimizing system is in computation burden, and the present invention only realizes the estimation to three-dimensional scenic at the key frame place, and upgrades 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}{\mathrm{\&Sigma;}}c\left(x\right)$

$c\left(x\right)=\left\{\begin{array}{cc}1,& |I_r\left(x\right)-I_i\left(x^{\&prime;}\right)|<\mathrm{\&zeta;}\\ 0,& \mathrm{others}\end{array}\right.$

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.

The invention has the beneficial effects as follows: the present invention adopts the RGB camera to obtain environmental data.At utilizing high precision monocular location algorithm, a kind of degree of depth drawing generating method based on variation mechanism is proposed, realized large-scale quick indoor 3 D scene rebuilding, solved three-dimensional reconstruction algorithm cost and real-time problem effectively.

Description of drawings

Fig. 1 is the indoor method for reconstructing three-dimensional scene process flow diagram based on the variation model;

Fig. 2 is synoptic diagram for camera coordinates;

Fig. 3 is the three-dimensional reconstruction experimental result of application example of the present invention.

Embodiment

Fig. 1 is based on the indoor method for reconstructing three-dimensional scene process flow diagram of variation model, may further comprise the steps:

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

Step 2 is set up the camera pose and is described and the camera projection model.

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.

Provide an application example of the present invention below.

The RGB camera that this example adopts is Point Grey Flea2, and image distinguishes that rate is 640 * 480, and the highest frame frequency is 30fps, and the horizontal field of view angle is 65 °, and focal length is approximately 3.5mm.Employed PC is equipped with GTS450GPU and i5 four nuclear CPU.

In experimentation, obtain the environment depth information by color camera, combining camera pose algorithm for estimating is realized self accurate location.After entering key frame, 20 two field pictures are as the input of this paper depth estimation algorithm around the selection key frame.In the depth estimation algorithm implementation, make d ⁰=h ⁰And q ⁰=0, calculate

Import with the initialization of obtaining current depth map, and iteration optimization E _{D, q}With E _hUp to convergence.Simultaneously, should in the algorithm iteration process, constantly reduce the θ value, increase the weight of quadratic function in the algorithm implementation, effectively improve algorithm the convergence speed.Final experimental result as shown in Figure 3, experiment shows that this method can effectively realize the dense three-dimensional reconstruction of environment, the step of going forward side by side has been demonstrate,proved the feasibility of this method.

Claims (3)

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=\left[\begin{array}{ccc}{f}_u& 0& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right]$

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:

$\left\{\begin{array}{c}x=(u-u_0)*z/f_u\\ y=(v-v_0)*z/f_v\\ z=d\end{array}\right.$

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=\left[\begin{array}{c}{u}_0\\ v_0\end{array}\right]+\left[\begin{array}{cc}{f}_u& 0\\ 0& f_v\end{array}\right]\frac{r_d}{r_u}{x}_u$

$r_d=\frac{1}{\mathrm{\&omega;}}\arctan \left({2r}_u\tan \frac{\mathrm{\&omega;}}{2}\right)$

$r_u=\frac{\tan \left(r_d\mathrm{\&omega;}\right)}{2\tan \frac{\mathrm{\&omega;}}{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_{\mathrm{cw}}=\left[\begin{array}{cc}{R}_{\mathrm{cw}}& t_{\mathrm{cw}}\\ 0& 1\end{array}\right]$

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_{\mathrm{cw}}=\exp \left(\widehat{\mathrm{\&mu;}}\right)$

$\widehat{\mathrm{\&mu;}}=\left[\begin{array}{cccc}0& {\mathrm{\&mu;}}_6& -{\mathrm{\&mu;}}_5& {\mathrm{\&mu;}}_1\\ {\mathrm{\&mu;}}_6& 0& {\mathrm{\&mu;}}_4& {\mathrm{\&mu;}}_2\\ {\mathrm{\&mu;}}_5& -{\mathrm{\&mu;}}_4& 0& {\mathrm{\&mu;}}_3\\ 0& 0& 0& 0\end{array}\right]$

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}{\mathrm{\&Sigma;}}c\left(x\right)$

$c\left(x\right)=\left[\begin{array}{cc}1,& |I_r\left(x\right)-I_i\left(x^{\&prime;}\right)|<\mathrm{\&zeta;}\\ 0,& \mathrm{others}\end{array}\right]$

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:

$\mathrm{\&xi;}=\underset{\mathrm{\&xi;}}{\arg \min }\mathrm{\&Sigma;Obj}(\frac{\left|e_j\right|}{{\mathrm{\&sigma;}}_j},{\mathrm{\&sigma;}}_T)$

$e_j=\left(\begin{array}{c}{u}_i\\ v_i\end{array}\right)-\mathrm{K\&pi;}\left(\exp \left(\mathrm{\&xi;}\right)p\right)$

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:

$\left\{\right\{{\mathrm{\&xi;}}_2..{\mathrm{\&xi;}}_N\},\{p_1..p_M\left\}\right\}=\underset{\left\{\right\{\mathrm{\&mu;}\},\{p\left\}\right\}}{\arg \min }\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j\&Element;S_i}{\mathrm{\&Sigma;}}\mathrm{Obj}(\frac{\left|e_{\mathrm{ji}}\right|}{{\mathrm{\&sigma;}}_{\mathrm{ji}}},{\mathrm{\&sigma;}}_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_{\mathrm{data}}=\frac{1}{\left|I\left(r\right)\right|}\underset{I_i\&Element;I}{\mathrm{\&Sigma;}}|I_r\left(x\right)-I_i\left(x^{\&prime;}\right)|$

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^{\&prime;}={\mathrm{\&pi;}}^{-1}\left({\mathrm{KT}}_r^i\mathrm{\&pi;}(x,d)\right)$

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:

${\left|\right|x\left|\right|}_{\mathrm{\&alpha;}}=\left\{\begin{array}{c}\frac{{\left|\right|x\left|\right|}^2}{2\mathrm{\&alpha;}},\left|\right|x\left|\right|\&le;\mathrm{\&alpha;}\\ \left|\right|x\left|\right|-\frac{\mathrm{\&alpha;}}{2},\mathrm{others}\end{array}\right.$

In the formula, α is constant;

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

$g{\left|\right|\&dtri;d\left|\right|}_{\mathrm{\&alpha;}}=<g\&dtri;d,q>-\mathrm{\&delta;}\left(q\right)-\frac{\mathrm{\&alpha;}}{2}{\left|\right|q\left|\right|}^2$

In the formula, $\mathrm{\&delta;}\left(q\right)=\left\{\begin{array}{cc}\frac{\mathrm{\&alpha;}}{2}& \mathrm{\&alpha;}<\left|\right|q\left|\right|\&le;1\\ \&infin;& \mathrm{other}\end{array}\right.$

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 }\left\{\underset{d}{\arg \min }{E}_{d,q}\right\}$

$E_{d,q}={\&Integral;}_{\mathrm{\&Omega;}}(<g\&dtri;d,q>+\frac{1}{2\mathrm{\&theta;}}{(d-h)}^2-\mathrm{\&delta;}\left(q\right)-\frac{\mathrm{\&alpha;}}{2}{\left|\right|q\left|\right|}^2)\mathrm{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-\mathrm{\&alpha;q}=0$

$\frac{{\&PartialD;E}_{d,q}}{\&PartialD;d}=g\mathrm{div}q+\frac{1}{\mathrm{\&theta;}}(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}{{\mathrm{\&epsiv;}}_q}=g\&dtri;d-{\mathrm{\&alpha;q}}^{n+1}$

$\frac{d^{n+1}-d^n}{{\mathrm{\&epsiv;}}_d}=g\mathrm{div}p+\frac{1}{\mathrm{\&theta;}}(d^{n+1}-h)$

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

$p^{n+1}=\frac{(p^n+{\mathrm{\&epsiv;}}_{q}g\&dtri;d^n)/(1+{\mathrm{\&epsiv;}}_q\mathrm{\&alpha;})}{\max (1,(p^n+{\mathrm{\&epsiv;}}_{q}g\&dtri;d^n)/(1+{\mathrm{\&epsiv;}}_q\mathrm{\&alpha;}))}$

$d^{n+1}=\frac{d^n+{\mathrm{\&epsiv;}}_d({g\mathrm{div}q}^{n+1}+h^n/\mathrm{\&theta;})}{(1+{\mathrm{\&epsiv;}}_d/\mathrm{\&theta;})}$

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;}_{\mathrm{\&Omega;}}(\frac{\mathrm{\&theta;}}{2}{(d-h)}^2+\frac{\mathrm{\&lambda;}}{\left|I\left(r\right)\right|}\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}|I_i\left(x\right)-I_{\mathrm{ref}}(x,h)|)\mathrm{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_{\mathrm{inc}}^k=\frac{{\mathrm{Sd}}_{\min }{d}_{\max }}{(S-k)d_{\min }+d_{\max }}$

In the formula,

Be k and k-1 sample plane interval.

Images