CN106780459A - A kind of three dimensional point cloud autoegistration method - Google Patents

Abstract

The invention discloses a kind of three dimensional point cloud autoegistration method, comprise the following steps：Two panels point cloud subject to registration sample and obtains characteristic point, invariable rotary feature of these characteristic points is calculated respectively, invariable rotary feature to characteristic point in two panels point cloud carries out matching search, obtains the initial corresponding relation between characteristic point；Then Mismatching point present in first match point collection is judged and is removed using stochastical sampling unification algorism that the feature point correspondence for being optimized is calculated between two panels point cloud and is generally rigid transformation relation, realized rough registration；And a kind of detection algorithm of rigid transformation uniformity is provided, and restrictive detection is carried out to rough registration result using the local coordinate system transformation relation between matching characteristic point, complete the checking of rough registration result correctness；And using the rigid transformation relation between ICP algorithm optimization cloud data, the final automatic accuracy registration for realizing point cloud.

Description

A kind of three dimensional point cloud autoegistration method

Technical field

The invention belongs to three-dimensional measurement field, more particularly, to a kind of three dimensional point cloud autoegistration method.

Background technology

Threedimensional model all has great importance in fields such as industrial detection, historical relic's protection, biomedicines.As three-dimensional is surveyed The development of amount technology, the treatment technology to point cloud model has turned into the focus of Recent study.An important step in points cloud processing Rapid is exactly that the same coordinate system, i.e. point cloud registering are arrived in the cloud data unification for obtaining different points of view in Same Scene.

Point cloud registration method is broadly divided into following a few classes：(1) manual registration：The corresponding points between two panels point cloud are chosen manually Split is realized, the method needs the effect of manual intervention, split to often rely on the experience of operator；(2) based on external auxiliary Method for registering：Using turntable or on measured object, the mode of sticking sign point realizes auxiliary registration, and such method will to environment Ask higher, the precision of registration is easily influenceed by external auxiliary；(3) method for registering based on itself shape characteristic：Using measured object The shape characteristic of body itself realizes the matching between point cloud and then realizes the autoregistration of point cloud, without using extra equipment or Person's index point, disclosure satisfy that the demand of different measurement occasions, but the stability and precision of algorithm and the point cloud registering skill for having auxiliary Art is compared to also to be strengthened.

For the method for registering based on itself shape characteristic, existing research unit carries out correlative study, and achieves certain Achievement；Rusu by the use of characteristic point and surrounding neighbors point normal vector angle as feature construction a kind of quick point feature Nogata Figure (Fast Point Feature Histograms, FPFH) is used for the lookup and matching of characteristic point, completes initial registration and obtains After comparatively ideal initial position, registration accuracy is improved using accuracy registration algorithm；Chen and Besl propose iteration closest approach Algorithm (Iterative closest point, ICP) carries out the accuracy registration of two amplitude point cloud data, and the method passes through iteration The closest approach of corresponding points concentration is found as corresponding points, rigid transformation matrix is continued to optimize, accuracy registration is finally obtained.But the party Initial position requirement of the method to two panels point cloud is higher, it is necessary to initial relative position is substantially close between putting cloud.

Existing point cloud registration method generally by the way of rough registration and accuracy registration are combined, wherein, rough registration one As find corresponding points, the initial positional relationship that estimation obtains between two amplitude point clouds by calculating the feature of point；Essence is to use with criterion ICP and its innovatory algorithm are further optimized to rough registration result, realize the accuracy registration of cloud data；But existing utilization The feature of point carry out the method for point cloud registering for the disturbing factors such as a cloud noise, external isolated point and point cloud variable density compared with It is sensitivity, substantial amounts of Mismatching point easily occurs, reduces the precision and stability of registration.

The content of the invention

It is an object of the invention to be directed to existing technical problem, there is provided a kind of three dimensional point cloud autoegistration method, Registration accuracy is high and with very strong stability, completes the automatic decision of registration result correctness.

To achieve the above object, technical solution of the invention is that the point cloud based on invariable rotary Feature Descriptor is automatic Method for registering, first two panels point cloud subject to registration sample obtains characteristic point, and the rotation of these characteristic points is calculated respectively not Become feature, the invariable rotary feature to characteristic point in two panels point cloud carries out matching search, obtain initial right between characteristic point Should be related to；Then first match point is concentrated using stochastical sampling unification algorism (Random Sample Consensus, RANSAC) The Mismatching point of presence is judged and is removed that the feature point correspondence for being optimized is calculated between two panels point cloud Transformation relation is generally rigid, rough registration is realized；And a kind of detection algorithm of rigid transformation uniformity is provided, using matching characteristic point Between local coordinate system transformation relation restrictive detection is carried out to rough registration result, carry out the checking of rough registration result correctness； It is final to realize that point the automatic of cloud is accurately matched somebody with somebody finally using the rigid transformation relation between iteration closest approach algorithm optimization cloud data It is accurate.

Specifically include following steps：

Step 1：Source point cloud P subject to registration and impact point cloud Q that reading is collected；

Step 2：The density of source point cloud P and impact point cloud Q is calculated respectively, and is randomly selected from the source point cloud P some Individual point composition source feature point set S₁, several points are randomly selected from impact point cloud Q and constitute target signature point set S₂；

Step 3：According to the density of source point cloud and impact point cloud, source feature point set S is calculated respectively₁With target signature point set S₂ In each characteristic point local rotation translation invariant geomagnetic coordinates；

Step 4：Local rotation translation invariant geomagnetic coordinates according to characteristic point calculates the high dimensional feature description of each characteristic point； Feature Points Matching is carried out to impact point cloud Q and source point cloud P, just matching point set C is obtained；

Step 5：The mistake that just match point is concentrated obtained in the step 4 is removed using stochastical sampling unification algorism Match somebody with somebody, calculated using the variation in rigidity algorithm for estimating based on singular value decomposition method and obtain spin matrix R and translation matrix T, obtain institute State the rough registration result between impact point cloud and the source point cloud；

Step 6：The rough registration knot obtained using the registering arithmetic code error detection algorithm determination step 5 based on rigid transformation uniformity Whether fruit is correct, is as a result correctly then transferred to step 7, the as a result incorrect result for then returning to registration failure；

Step 7：Spin matrix R and translation matrix T that above-mentioned estimation is obtained are optimized using iteration closest approach algorithm iteration, it is real The automatic accuracy registration of the existing impact point cloud and source point cloud.

In general, by the contemplated above technical scheme of the present invention compared with prior art, can obtain down and show Beneficial effect：

(1) the three dimensional point cloud autoegistration method that the present invention is provided, by the selected characteristic point set from cloud data, For each characteristic point builds local rotation translation invariant geomagnetic coordinates, and the description of the high dimensional feature on this feature point is calculated, The replacement to Traditional Man index point is realized, be can be used in the middle of the application scenarios such as cultural relic digitalization scanning, can subtracted Few artificial workload for pasting, removing index point；

(2) the three dimensional point cloud autoegistration method that the present invention is provided, it is special using resulting feature point set and higher-dimension Description is levied, can effectively realize that the corresponding points between different cloud datas are searched and matched, realize putting cloud under arbitrary initial attitude The autoregistration of data；

(3) the three dimensional point cloud autoegistration method that the present invention is provided, by will be using part invariable rotary coordinate system The local spin matrix result being calculated and the global rotation translation matrix obtained using RANSAC algorithm solution Whether result is compared, and judges there is rotation translation uniformity between it, can with so that effectively whether judge rough registration result Success, realizes accuracy registration.

Brief description of the drawings

Fig. 1 is the overall flow schematic diagram of the three dimensional point cloud autoegistration method provided according to the present invention.

Specific embodiment

In order to make the purpose , technical scheme and advantage of the present invention be clearer, it is right below in conjunction with drawings and Examples The present invention is further elaborated.It should be appreciated that specific embodiment described herein is only used to explain the present invention, not For limiting the present invention.As long as additionally, technical characteristic involved in invention described below each implementation method that Conflict is not constituted between this can be just mutually combined.

The three dimensional point cloud autoegistration method that the embodiment of the present invention is provided, its flow is as shown in figure 1, including as follows Step：

Step 1：Impact point cloud and source point cloud that reading is collected by 3-D measuring apparatus；

Step 2：The several points randomly selected respectively in source point cloud and impact point cloud constitute two groups of feature point sets, specifically For：

According to default characteristic point oversampling ratio h₁And h₂, stochastical sampling is obtained from source point cloud P and impact point cloud Q respectively Corresponding feature point set

Wherein, n₁Refer to that the feature of the feature point set of source point cloud P is counted out, n₂Refer to the spy of the feature point set of impact point cloud Q Levy and count out, n₁=h₁·N₁, n₂=h₂·N₂。

Step 3：The feature point set obtained to sampling, calculates the local rotation translation invariant coordinates of wherein each characteristic point System；Specifically include following sub-step：

Step 3.1：For each point p in source point cloud_iWith each point q in impact point cloud Q_i, using nearest neighbor search side Method calculates point p_iNearest neighbor point p_i', and point q_iNearest neighbor point q_i′

Obtain point cloud p_iPoint cloud density d_p, and point cloud q_iPoint cloud density d_q；

Wherein, | | p_i-p_i' | | it is to give directions p_iWith p_i' between Euclidean distance, | | q_i-q_i' | | point q between referring to_iWith point q_i′ Between Euclidean distance；I=1,2,3 ..., N, N are the number at cloud data midpoint；N is included in source point cloud P₁It is individual, to source point For cloud P, N=N₁；Q includes N in impact point cloud₂It is individual, for impact point cloud Q, N=N₂；

Step 3.2：Any selected characteristic point concentrates a point p_i, construction point p_iCovariance matrix：

Wherein, w_i=1/ | | p_j-p_i||,Radius r_locIt is p_iNeighborhood calculate half Footpath；In embodiment, r is taken_loc=15d_p；

Step 3.3：The characteristic value and characteristic vector of covariance matrix are solved using following formula；

COV(p_i) V=EV

Wherein, E be by3 × 3 diagonal matrixs for constituting

It is characteristic valueCorresponding characteristic vector, l=1,2,3；

Step 3.4：WillIt is set to x⁺,y⁺,z⁺,Opposite vector be set to x^-,y^-,z^-；

The direction of z-axis is determined by following formula in local coordinate system：

Step 3.5：X-axis direction in local coordinate system is determined using the method for step 3.4；Y-axis direction in local coordinate system Determined by z × x, thus set up with p_iIt is the local rotation translation invariant geomagnetic coordinates F of origin_i={ x, y, x },；Using with it is above-mentioned The same method in step 3.2~3.4, sets up with point q_iIt is the local rotation translation invariant geomagnetic coordinates of origin.

Step 4：According to the local rotation translation invariant geomagnetic coordinates of the characteristic point being calculated, the higher-dimension for calculating characteristic point is special Description is levied, Feature Points Matching is carried out to two amplitude point clouds, obtain just matching point set；Specifically include following sub-step：

Step 4.1：To any feature point p in source point cloud feature point set_i, its high dimensional feature is described as high dimension vector form In the higher dimensional space that the feature geometries corresponding to impact point cloud feature point set are constituted search withClosest feature And secondary near feature

In the present invention, the search of arest neighbors feature is by quick approximate KNN searching method (Fast Approximate Nearest Neighbor Search) carry out speed-up computation；

Step 4.2：Calculate featureNearest feature is arrived respectivelySecondary near featureEuclidean distanceWithMeter CalculateWithRatio, and judge size of the ratio with default distance than threshold tau；

Specifically, judged to obtain according to following formulaWithBetween whether have correct corresponding relation

IfThen showWithThe match is successful, and otherwise it fails to match；IfWithThe match is successful, then its Corresponding source point cloud is also correctly corresponding with the characteristic point in impact point cloud respectively, and the point is to constituting a matching double points；

Other characteristic points to source point cloud are also carried out above-mentioned treatment, and just matching point set C is made up of all of matching double points.

Step 5：Using stochastical sampling unification algorism (Random Sample Consensus, RANSAC) removal just matching The error hiding that point is concentrated, spin matrix R and translation matrix T is calculated using singular value decomposition method (SVD), obtains source point cloud and target Rough registration relation between point cloud, specifically includes following sub-step：

Step 5.1：Set the number of times S of stochastical sampling_numWith initial residual error E_rr, adopted at random from just matching point set C Sample, chooses 3 pairs of match points as initial point；

If the match point of certain characteristic point chosen has multiple, wherein conduct is randomly selected from corresponding point set should Search the corresponding points of point；

Step 5.2：Spin matrix R and translation matrix T are solved using svd algorithm, it is specific as follows：

Step 5.21：If it is P={ p to obtain matching point set by step 5.1₁,p₂,p₃, P '={ p '₁,p′₂, p '₃, utilize Following formula calculates the barycenter of point set；

Wherein, p_kWith p '_kIt is any three-dimensional coordinate to match point；

Step 5.22：Point set P and P ' is translated relative to respective barycenter using following formula, obtains new point set Q={ q₁,q₂, q₃, and Q '={ q '₁,q′₂,q′₃}；

q_k=p_k-g,q′_k=p '_k- g ', (k=1,2,3)；

Step 5.23：3 × 3 matrix M is calculated using following formula：

Step 5.24：Singular value decomposition M=U Λ V are carried out to Metzler matrix^T；

Wherein, subscript T is the transposition of matrix, and U, V are 3 × 3 unitary matrice, and Λ is 3 × 3 diagonal matrix, defines 3 × 3 Diagonal matrix A is：

Step 5.25：The translation matrix T of 3 × 3 spin matrix R and 3 × 1 is calculated using following formula：

R=UAV^T, T=g '-Rg；

Step 5.3：All of matching double points are concentrated to first match point, the range error after rotation translation is calculated according to following formula d_err；

d_err=| | p_τ-(R·q_τ+T)||₂；

If range error d_errLess than given initial residual error E_rr, then judge that the match point is interior point；According to above-mentioned side The number m ' that method is obtained all interior points and put in counting；

In embodiment, residual error is distance threshold；(p_τ,q_τ) be just match point concentrate a matching double points, τ= 1,2 ... r, r are the number of the just matching double points that match point is concentrated；

Step 5.4：If interior number m ' is adopted more than interior given number m, the point concentrated using the interior point data Solve spin matrix R, translation matrix T again with svd algorithm；

Residual error E ' is calculated according to following formula_rr；

If residual error is less than initial residual error, the rotation translation matrix that will be now calculated is used as best estimate Objective model parameter, updates rotation translation matrix R, T and initial residual error E '_rr；Otherwise go to step 5.5；If interior point is individual Number then passes directly to step 5.5 less than the number of given interior point；

Step 5.5：Repeat stochastical sampling S_numSecondary, repeat step 5.2~5.4 obtains the interior point corresponding to group sampling Number, the interior points to all sampling are ranked up, and choose the most sampling resultses of interior points as optimal sampling；

Using the interior point data collection C obtained under the optimal sampling, spin matrix R, translation matrix T are solved according to step 5.2, As optimal rotation translation matrix R_ranAnd T_ran；Wherein C₁In comprising S to matching double points.

Step 6：Whether the detection algorithm verification rough registration relation using rigid transformation uniformity is correct；

Step 6.1：Matching double points (the p in just matching point set C is calculated using following formula_loc,q_loc) rotation between its match point Turn translation matrix (R_loc,T_loc)：

Wherein, R_locIt is the local spin matrix between two match points, T_locIt is the part translation between two match points Matrix；

F_ploc、F_qlocIt is respectively in point p_loc,q_locThe local coordinate system of foundation, loc=1,2 ... S, S are in matching point set C The number of matching double points；

Step 6.2：The spin matrix R for obtaining will be solved_locBe converted to Eulerian angles to represent, i.e.,

R_loc→(α_loc,β_loc,γ_loc)；The spin matrix R that will be obtained in step 5.6_ranIt is converted into Eulerian angles to represent, i.e., R_ran→(α_ran,β_ran,γ_ran)；

Step 6.3：Angle difference d between Eulerian angles is calculated according to following formula_a：

Wherein, Δ (η₁,η₂)²=(η₁-η₂)²,η₁,η₂It is Eulerian angles；

Step 6.4：Calculated according to following formulaWith T_ran=(t_x,t_y,t_z)^TThe distance between Difference d_t：

Step 6.5：According to d_aAnd d_tWhether given threshold value σ is respectively less than_aAnd σ_t, judge the rotation that two methods are solved Whether translation matrix is consistent；

In embodiment, σ_aTake 0.5236, σ_tTake 15D_den, D_denIt is a cloud density, i.e., spin matrix difference is less than 30 °, translation When matrix distance is less than 15 times of point cloud density, then show what rotation translation matrix that this step calculated and step 5 were obtained Rotation translation matrix is consistent；

Step 6.6：Step 6.1 is also performed to other matching double points in the interior point data collection C that is obtained under optimal sampling to arrive Step 6.5, obtains the local rotation translation matrix and above-mentioned optimal rotation translation matrix (R between all matching double points_ran,T_ran) Between comformity relation；

Final statistics obtains all matching double points for meeting comformity relation in the interior point data collection C obtained under optimal sampling Number s；

Step 6.7：Calculating meets the number s and the number for matching matching double points in point set C of the matching double points of comformity relation The uniformity ratio λ=s/S of mesh S；

If λ >=threshold tau_λ, then show to be solved using RANSAC the spin matrix and the local invariable rotary coordinate of utilization for obtaining The spin matrix that system is calculated is consistent, then be judged to registration success；Otherwise, it is determined that registration failure；In embodiment, threshold value τ_λ=0.7.

Step 7：Using the rigid transformation relation between improved ICP algorithm optimization point cloud, the automatic accurate of point cloud is realized Registration；Rigid transformation relation includes rotation translation matrix；Comprise the following steps that：

Step 7.1：Setpoint distance threshold value ω as iteration ends condition；Wherein, ω>0；Distance threshold ω is according to source point The point cloud density d of cloud P_p；In embodiment, ω=5d is taken_p；

Step 7.2：If randomly selected in source point cloud doing as point to be matched；

Step 7.3：The corresponding points of point to be matched in searching source point cloud in impact point cloud with backwards projection method；

Step 7.4：Face distance metric is arrived as the object function solved needed for ICP algorithm, continuous iteration using based on point Calculate source point cloud to the rigid transformation relation of impact point cloud；

Step 7.5：When target function value is less than distance threshold ω, stop iteration；And will now solve the rigidity for obtaining Transformation relation completes the automatic accuracy registration of point cloud as final result.

As it will be easily appreciated by one skilled in the art that the foregoing is only presently preferred embodiments of the present invention, it is not used to The limitation present invention, all any modification, equivalent and improvement made within the spirit and principles in the present invention etc., all should include Within protection scope of the present invention.

Claims (7)

1. a kind of three dimensional point cloud autoegistration method, it is characterised in that comprise the following steps：

Step 1：Source point cloud P subject to registration and impact point cloud Q that reading is collected；

Step 2：The density of the source point cloud P and impact point cloud Q is calculated respectively, and is randomly selected from the source point cloud P some Individual point composition source feature point set S₁, several points are randomly selected from impact point cloud Q and constitute target signature point set S₂；

Step 3：According to described density, source feature point set S is calculated respectively₁With target signature point set S₂In each characteristic point office Portion's rotation translation invariant geomagnetic coordinates；

Step 4：The high dimensional feature description of each characteristic point is calculated according to the local rotation translation invariant geomagnetic coordinates；To impact point Cloud Q and source point cloud P carries out Feature Points Matching, obtains just matching point set C；

Step 5：The error hiding of C is concentrated using the stochastical sampling unification algorism removal just match point, using singular value decomposition method Calculate and obtain spin matrix R and translation matrix T, obtain the rough registration result of the impact point cloud and the source point cloud；

Step 6：Judge whether the rough registration result meets rigid transformation condition for consistence, if so, then entering step 7；If it is not, Then judge registration failure；

Step 7：The spin matrix R and translation matrix T is optimized using iteration closest approach algorithm iteration, obtain impact point cloud Q with The accuracy registration result of source point cloud P.

2. three dimensional point cloud autoegistration method as claimed in claim 1, it is characterised in that the step 2 includes as follows Sub-step：

Step 2.1：The nearest neighbor point of each point in source point cloud P and impact point cloud Q is calculated using nearest neighbor search method；

And calculate the point cloud density of each point cloud

Wherein, | | calculate point-nearest neighbor point | | refer to calculate the Euclidean distance between point and nearest neighbor point；N is included in source point cloud P₁ It is individual, for source point cloud P, N=N₁；Q includes N in impact point cloud₂It is individual, for impact point cloud Q, N=N₂；

Step 2.2：Setting characteristic point oversampling ratio h₁And h₂, stochastical sampling obtains source respectively from source point cloud P and impact point cloud Q The feature point set of point cloud PWith the feature point set of impact point cloud Q

Wherein, i is the numbering of characteristic point in feature point set；n₁Refer to that the feature of the feature point set of source point cloud P is counted out, n₂Refer to The feature of the feature point set of impact point cloud Q is counted out, n₁=h₁·N₁, n₂=h₂·N₂。

3. three dimensional point cloud autoegistration method as claimed in claim 2, it is characterised in that the step 3 is included such as Lower sub-step：

Step 3.1：Arbitrarily choose the feature point set S of the source point cloud P₁In a point p_i, according to p_iPoint cloud density obtain p_i's Neighborhood calculates radius, using following formula construction point p_iCovariance matrix：

$COV\left(p_i\right)=\underset{\left|\right|p_j-p_i\left|\right|<r_{loc}}{\&Sigma;}{w}_i(p_j-p_i){(p_j-p_i)}^T/\underset{\left|\right|p_j-p_i\left|\right|<r_{loc}}{\&Sigma;}{w}_i$

Wherein, w_i=1/ | | p_j-p_i| |,Radius r_locIt is p_iNeighborhood calculate radius；

To the feature point set S of impact point cloud Q₂In a point q_i, it is also carried out that acquisition q is processed as above_iCovariance matrix COV (q_i)；

Step 3.2：Solve the characteristic value and characteristic vector of covariance matrix；

COV(p_i) V=EV

Wherein, E be by3 × 3 diagonal matrixs for constituting It is characteristic valueCorresponding feature Vector, l=1,2,3,

Step 3.3：According to the characteristic vectorThe x of the local rotation translation invariant geomagnetic coordinates of structure, y, z coordinate axle, Set up with point p_iIt is the local rotation translation invariant geomagnetic coordinates F of origin_i={ x, y, x }；

Step 3.4：According to the feature point set S of impact point cloud Q₂Midpoint q_iCovariance matrix, using with step 3.2~3.3 Method, set up with point q_iIt is the local rotation translation invariant geomagnetic coordinates of origin.

4. three dimensional point cloud autoegistration method as claimed in claim 3, it is characterised in that the step 4 is included such as Lower sub-step：

Step 4.1：To any feature point p in the source point cloud feature point set_i, its high dimensional feature is described as high dimension vector form In the higher dimensional space that the feature geometries corresponding to impact point cloud feature point set are constituted search withClosest featureAnd Secondary near feature

Step 4.2：Calculate featureNearest feature is arrived respectivelySecondary near featureEuclidean distanceWith

Judged according to following formulaWithBetween whether have correct corresponding relation

IfThen showWithThe match is successful, and otherwise it fails to match；IfWithThe match is successful, then it is right respectively The source point cloud answered is also correctly corresponding with the characteristic point in impact point cloud, and the point is to constituting a matching double points；By all of matching Point is to constituting just matching point set.

5. three dimensional point cloud autoegistration method as claimed in claim 4, it is characterised in that the step 5 is included such as Lower sub-step：

Step 5.1：Set the number of times S of sampling_numWith distance threshold μ；Every time during sampling, random choosing is concentrated from the just match point T is taken to match point as initial point；

Step 5.2：Spin matrix R and translation matrix T is solved using singular value decomposition method；

Step 5.3：All of match point is concentrated to first match point, the range error d after rotation translation is calculated_err；If range error d_errLess than the distance threshold μ, then judge that the match point is interior point, the number m ' for obtaining all interior points and being put in counting；

Wherein, d_err=| | p_i-(R·q_i+T)||₂, (p_i,q_i) be just match point concentrate a pair of match points, i=1,2 ... r, r It is the quantity of the matching double points that first match point is concentrated；

Step 5.4：Repeat stochastical sampling S_numIt is secondary, to sampling results each time, this sampling institute is obtained according to step 5.2~5.3 Corresponding interior points；

Interior points to all sampling are ranked up, and choose the most sampling resultses of interior points as optimal sampling；Using optimal The interior point data collection C obtained under sampling₁, the method according to step 5.2 obtains spin matrix R, translation matrix T, used as optimal rotation Turn translation transformation matrix, be denoted as R_ranAnd T_ran。

6. three dimensional point cloud autoegistration method as claimed in claim 5, it is characterised in that the step 5.2 is included such as Lower sub-step：

Step 5.21：The barycenter of just matching point set is calculated according to following formula；

$p_0=\frac{1}{t}{\mathrm{\&Sigma;}}_{k=1}^t{p}_k,q_0=\frac{1}{t}{\mathrm{\&Sigma;}}_{k=1}^t{q}_k;$

Wherein, just matching point set is P={ p₁,p₂,…,p_t, P '={ q₁,q₂,…,q_t}；

Step 5.22：Just matching point set P and P ' is translated relative to the barycenter obtained in the step 5.21, new point set is obtained P '={ p '₁,p′₂,…,p′_tAnd Q '={ q '₁,q′₂,…,q′_t}；

p′_k=p_k-p₀,q′_k=q_k-q₀

Step 5.23：Calculate the matrix M of t × t：

$M={\mathrm{\&Sigma;}}_{k=1}^t{p}_k^{\&prime;}{q}_k^{\&prime;T};$

Step 5.24：Singular value decomposition M=U Λ V are carried out to Metzler matrix^T, wherein subscript T is the transposition of matrix, and U, V are 3 × 3 tenth of the twelve Earthly Branches Matrix, Λ is 3 × 3 diagonal matrix, defines 3 × 3 diagonal matrix A：

$A=\left\{\begin{array}{cc}{I}^3\det \left(U\right)& \det \left(V\right)\&GreaterEqual;0\\ diag(1,1,-1)& \det \left(U\right)\det \left(V\right)<0\end{array}\right.$

Step 5.25：Calculate the spin matrix R of the t × t and translation matrix T of t × 1：

R=UAV^T, T=q₀-R·p₀。

7. three dimensional point cloud autoegistration method as claimed in claim 5, it is characterised in that the step 6 is included such as Lower sub-step：

Step 6.1：To the interior point data collection C obtained under the optimal sampling₁In each matching double points be handled as follows：

To the interior point data collection C₁In a matching double points (p_loc,q_loc), solve the local spin moment between two match points Battle array R_locWith local translation matrix T_loc：

$R_{loc}={\left(F_{p_{loc}}\right)}^T{F}_{q_{loc}}$

T_loc=p_loc-q_locR_loc

WhereinIt is respectively in point p_loc、q_locThe local coordinate system of foundation, loc=1,2 ... S, S are the interior points According to collection C₁The number of middle matching double points；

Step 6.2：Calculate the spin matrix R of this step acquisition_locThe spin matrix R obtained with step 5.4_ranBetween angle away from From：

$d_a=arccos\left(\frac{trace\left(R_{loc}{R}_{ran}^T\right)-1}{2}\right);$

Step 6.3：Calculate the translation matrix T of this step acquisition_locThe translation matrix T obtained with step 5.4_ranBetween translation away from From：

$d_t=\frac{\left|\right|T_{loc}-T_{ran}\left|\right|}{d_p};$

Step 6.4：If above-mentioned d_a、d_tRespectively less than given threshold value σ_a、σ_t, then the rotation translation that step 5 is solved with step 6 is shown Matrix is consistent；

Step 6.5：Internal point data collection C₁In other matching double points carry out the treatment of step 6.1~step 6.4, owned Rotation translation matrix (the R that local rotation translation matrix between matching double points is obtained with step 5_ran,T_ran) between uniformity Relation；

Point data collection C in statistics₁In all matching double points for meeting comformity relation number s；

Step 6.6：Calculate s and interior point data collection C₁The uniformity ratio λ=s/S of the number S of middle matching double points；

If λ >=threshold tau_λ, then profit in the local spin matrix that is calculated using local invariable rotary coordinate system and step 5 is shown It is consistent that the global rotation translation matrix for obtaining is solved with RANSAC algorithm, is judged to registration success；Otherwise, sentence It is set to registration failure.