CN109389625B - Three-dimensional image registration method based on multi-scale descriptor screening and mismatching - Google Patents

Abstract

The invention discloses a three-dimensional image registration method for screening out mismatching based on a multi-scale descriptor, which better describes the characteristics of corresponding key points and preliminarily obtains corresponding matching points by constructing a novel multi-scale descriptor; traversing a point cloud set with similar multi-scale descriptors in the point cloud to be registered by taking the point cloud to be registered as a characteristic, greatly improving the operation efficiency of rough point cloud registration, reducing the calculated amount of a computer and bringing great convenience to point cloud registration; the method can obtain a more accurate registration effect in a shorter time, has better robustness, and is suitable for the field of precision measurement with noise, complex structure and high registration requirement.

Description

Three-dimensional image registration method based on multi-scale descriptor screening and mismatching

Technical Field

The invention belongs to the field of image registration, and particularly relates to a three-dimensional image registration method for screening out mismatching based on a multi-scale descriptor.

Background

With the rapid development of electronic and computer technologies and the demand of people for high-quality products, foreign machine vision has entered the high-speed development period in the 90 s of the 20 th century and is widely applied to the field of industrial control. Machine vision is an essential means for industrial automation and intelligence, and is equivalent to the extension of human vision on a machine. The machine vision has the advantages of high automation, high efficiency, high precision and the like, and can play an important role in the field of automation in China. Compared with the traditional two-dimensional machine vision, the three-dimensional machine vision has higher precision and better detection effect, and can meet the more rigorous production process. However, the point cloud data acquired by the three-dimensional machine vision is large in amount, and no obvious topological relation exists between points, so that a misregistration pair is easy to occur in the registration process, the registration accuracy of the point cloud data is reduced, and the measurement and detection effects on the target object are influenced.

In the intelligent manufacturing processes of aerospace, automobile traffic, bridges, ships and the like, complex special-shaped curved surfaces are processed and detected, the curved surfaces are compositely bent and randomly stacked in space, the edge curvature is complex, and the two-dimensional machine vision cannot meet the high-precision and high-efficiency production and processing requirements. However, the accuracy of the existing point cloud registration algorithm is not high, the registration accuracy and efficiency are influenced by the existence of a large number of misregistration pairs, and the detection and reconstruction of the target object are also seriously influenced.

Currently, the most widely used solution to this problem is the RANSAC algorithm, which, although it works well in most cases, grows exponentially with the ratio η of outlier point pairs, and thus is very sensitive to the selection of sample points in the point cloud registration and less robust. In the registration process, the distribution of the density of the point cloud is irregular and the structure of the point cloud is not regular, so that the registration accuracy of the three-dimensional image is much lower than that of the two-dimensional image.

Disclosure of Invention

The invention provides a three-dimensional image registration method for screening out mismatching based on a multi-scale descriptor, which screens by using the space geometric relationship between registration point pairs, further can directly eliminate the mismatching points and improves the registration efficiency and the registration precision.

A three-dimensional image registration method based on multi-scale descriptor screening mismatching comprises the following steps:

step 1: three-dimensional point clouds X and Y of view-adjacent scenes are acquired, and a multi-scale descriptor (Q) is constructed for each point in the three-dimensional point clouds X and Y_i，N_i)；

Wherein L represents a neighborhood search radius mark of a three-dimensional point cloud X or Y, the values are 1, 2, …, L represents the total number of neighborhood search radii, and L is more than or equal to 2；T_lWhich represents a normalized vector of the vector,

λ_l1，λ_l2，λ_l3is a matrix C_lEigenvalues after SVD decomposition, lambda_l1≥λ_l2≥λ_l3，n_l1，n_l2，n_l3Are each lambda_l1，λ_l2，λ_l3The corresponding feature vector is used as a basis for determining the feature vector,

S_l＝{p_i|||p_i-p₀||≤r_l}，|S_li represents S_lMaximum value of p_iRepresenting a point p in a three-dimensional point cloud X or Y₀Neighborhood point of r_lRepresents a point p₀The l-th neighborhood search radius of (c),

r_Lfor a set maximum neighborhood search radius, n_lRespectively representing normal vectors corresponding to the search neighborhood radius l of the points in the three-dimensional point cloud, and taking the value as n_l3；

|S_lI is expressed as p₀Is the center of a sphere r_lThe maximum value of the distances from all points in the sphere to the sphere center;

step 2: traversing the multi-scale descriptors of each point in the three-dimensional point cloud X and Y, and respectively recording the points with the same multi-scale descriptors in the X and the Y as a key point set A_i0And B_i0；

And step 3: searching key point set A according to principal curvature characteristics_i0And B_i0Initial matching point set A in (1)_i1And B_i1；

Initial matching point set A_i1Point a in_i1And B_i1B in (1)_i1Is a matching point pair;

from the set of keypoints A_i0Selecting an arbitrary point a_iGo through B_i0Midpoint, select and point a_iMatched point b_jAn initial matching point set A is constructed with points satisfying the following conditions_i1And B_i1：

Wherein k is₁(a_i) And k₂(a_i) Respectively, at the point a_iAt point a_iThe principal curvature minimum and the principal curvature maximum, k, obtained from the neighborhood points of₁(b_j) And k₂(b_j) Respectively, at point b_jTo utilize point b_jThe principal curvature minimum and the principal curvature maximum, δ, obtained from the neighborhood points of₁Is a set first error threshold;

and 4, step 4: initial matching point set A using similar triangle threshold_i1And B_i1Screening the matched point pairs to obtain a registration point pair set A_i2And B_i2；

And 5: computing A based on SVD_i2And B_i2First rotation matrix R in between₁And a first translation matrix T₁Using a first rotation matrix R₁And a first translation matrix T₁Sequentially calculating A_i2Each point a in_i2Change point a of_3iWhile calculating A_i2Each point in B_i2Matching point b in (1)_2iEuclidean distance d from corresponding transformation point_3i；

A_i3＝R₁*A_i2+T₁

Wherein (a)_3ix,a_3iy,a_3iz) And (b)_i2x,b_i2y,b_i2z) Are respectively a_3iAnd b_i2Coordinates on the x, y, z axes;

step 6: computing

k is 1 as the initial value, and N is the current registration point pair set A_i2And B_i2Centering the registration point pair, and judging d_k+1-d_k<ξ₃If yes, using the current registration point pair set A_i2And B_i2And (3) as registration data of the three-dimensional point clouds X and Y, if not, making k equal to k +1, returning to the step 3, reselecting the neighborhood points of each point, and calculating the principal curvature.

Further, using a rotation angle threshold θ_mScreening the registration pairs again, and setting the current registration point pair A_i2And B_i2And (3) rejecting the registration point corresponding to the point which does not meet the following formula as an abnormal point or a noise point:

Δθ_i<θ_m

wherein, Delta theta_iRepresenting the current registration point pair set a_i2And B_i2The rotation angle difference between the ith pair of registration points;

current registration point pair set a_i2And B_i2The ith pair of registration point pairs (x)_i,y_i) The solving process of the rotation angle difference between the two is as follows:

1) set A from current registration point pair_i2And B_i2Two pairs of non-coplanar registration pairs are randomly selected, four non-coplanar points form a spherical surface, the spherical surface equation and the spherical center O are solved by analytic geometry, and a space coordinate system is established by taking the point O as an origin;

2) taking the Z axis of the established space coordinate system as a second rotation matrix R₂And R is a rotational axis of₂＝A₂B₂，A₂And B₂Respectively represent second rotation matrices R₂Decomposition matrix of, B₂Set A by current registration point pair_i2And B_i2One pair of registration point pairs (x)_k,y_k) Satisfies B₂x_k＝y_kThen, calculate to obtainTo obtain a₂The rotation axis of the system is Z axis, and simultaneously, according to the set registration distance threshold value xi, the registration angle threshold value epsilon is calculated_i，

3) Rotate by

A circular area enclosed when

Left boundary and y of the enclosed circular area_iWhen tangent, by

The rotation angle of the enclosed circular area is theta₁At this time, the product is made by

The center of the circle of the circular area is the point M which is formed by the points M and y_iA distance r of_sConstructing a one-dimensional quadratic equation, using the equation of the point M in the great circle O, simultaneously establishing two equation sets to obtain the x and y coordinate values of the point M, and solving theta₁；

Represents point x_iThrough B₂Rigid body transformation to form B₂x_iIs the center of a sphere, epsilon_iA circular area of radius;

because B₂Is a rigid body transformation rotation matrix in the presence of noise or outliers,

it is shown that: x is the number of_iThrough B₂Rigid body transformation, but because of noise or abnormal points, only one satisfactory range can be found, and the range is B₂x_iAs a center of a circle, epsilon_iIs a circle with a radiusAn area;

the great circle O is pointing x_iIn the space coordinate system, a circle with the point O as an origin is located;

4) continue to rotate

A circular area enclosed by

Right boundary and y of the enclosed circular area_iWhen the two-way pipe is tangent to each other,

the rotation angle of the enclosed circular area is theta₂At this time, the product is made by

The circle center of the circular area is point N, and the point N and the point y are_iA distance r of_sConstructing a one-dimensional quadratic equation, using the equation of N in the great circle O, simultaneously establishing two equation sets to obtain the x and y coordinate values of the point N, and solving theta₂To obtain an angle difference Delta theta_i；

Wherein r is_sIs that

Radius of a circular area enclosed, r, being B in a space coordinate system O₂x_iDistance to point O, r_s＝г*tanε_i。

Further, the calculation process of the principal curvature of each point in the three-dimensional point cloud is as follows:

step 3.1: fitting any point in the three-dimensional point cloud and a neighborhood point of the point to form a local quadric surface G (x, y);

G(x,y)＝ax²+bxy+cy²+dx+ey

step 3.2: performing derivation on the local quadric G (x, y) to obtain a first basic constant E, F, G and a second basic constant L, M, N;

E＝G_xx，F＝G_x*G_y，G＝G_yy；

L＝n·G_xx，M＝n·G_xy，N＝n·G_yy；

wherein G is_x、G_yEach G (x, y) is derived once for x and y, G_xx、G_yySecond derivative of x, y for G (x, y), G_xySequentially differentiating x and y for G (x, y), wherein n is a normal vector of the local quadric surface G (x, y) at any selected point;

step 3.3: calculating gaussian curvature K and mean curvature H:

step 3.4: calculating a principal curvature minimum k₁And a maximum value k of principal curvature₂：

Further, the initial matching point set A is thresholded by using similar triangles_i1And B_i1The process of screening the matching point pairs in (1) is as follows:

step 4.1: in A_i1In the method, 3 points a which are not collinear are randomly selected_i,a_j,a_kIn B_i1Find 3 points b corresponding to it_u,b_v,b_w；

Step 4.2: order (a)_i,b_u) And (a)_j,b_v) Is a known pair of matching points, the following calculations are made according to the triangle similarity principle:

△a_ij＝||a_i-a_j||；△a_ik＝||a_i-a_k||；△a_jk＝||a_j-a_k||；

△b_uv＝||b_u-b_v||；△b_uw＝||b_u-b_w||；△b_wv＝||b_w-b_v||；

step 4.3: calculating the similarity omega of the two triangles according to a triangle similarity principle;

wherein, Δ a_ij,△a_ik,△a_jk；△b_uv，△b_uw,△b_wvFor each side length, delta, of a triangle₂Is a set second error threshold;

step 4.4: if a is_i,a_j,a_k，b_u,b_v,b_wSatisfy omega<δ₂Then, consider a_kAnd b_wBelonging to a pair of registration point pairs.

Further, L is 4 and r₄＝0.004m。

Further, the delta₂The value is any one of 0.02, 0.05 and 0.1.

Further, a 3D contour scanning sensor is adopted to obtain three-dimensional point cloud of a scene adjacent to a visual angle.

Advantageous effects

The invention provides a three-dimensional image registration method for screening out mismatching based on a multi-scale descriptor, which better describes the characteristics of corresponding key points and preliminarily obtains corresponding matching points by constructing a novel multi-scale descriptor; traversing a point cloud set with similar multi-scale descriptors in the point cloud to be registered by taking the point cloud to be registered as a characteristic, greatly improving the operation efficiency of rough point cloud registration, reducing the calculated amount of a computer and bringing great convenience to point cloud registration; the method can obtain a more accurate registration effect in a shorter time, has better robustness, and is suitable for the field of precision measurement with noise, complex structure and high registration requirement.

After the accurate registration, a registered point cloud pair is obtained, and compared with the traditional registration method, the method has the defects that a misregistration pair and noise exist inevitably, and the accuracy and precision of the three-dimensional target during measurement and detection are influenced definitely. At present, the optimal registration is performed by using an ICP iterative method most, but the existence of a large number of misregistration pairs greatly increases the possibility that the algorithm falls into local optimization or the algorithm cannot converge, and seriously affects the reliability of registration. Therefore, the misregistration pairs are further screened out by introducing the threshold value based on the rotation angle, so that the number of the misregistration pairs and noise points is greatly reduced, the accuracy of a registration result is ensured, and the registration result has good robustness; on the other hand, as a part of misregistration pairs and noise points are screened, the input quantity of point cloud is reduced, the calculated quantity in the registration process is reduced, and the registration efficiency is improved.

Drawings

FIG. 1 is a schematic flow diagram of the process of the present invention;

FIG. 2 is a schematic diagram of the registration process using rotation angles according to the present invention, in which diagram (a) is a schematic diagram of a spatial coordinate system, and diagram (b) is a schematic diagram of a threshold ξ converted into an angle ε_iIn (c) is x_iAnd x_kAt each of R and R₂Schematic view under transformation, where (d) is B₂x_iThe rotation process of (a) is as shown in (e)

And the relation of the enclosed circular area in a space coordinate system is shown schematically.

Detailed Description

The invention will be further described with reference to the following figures and examples.

As shown in fig. 1, a three-dimensional image registration method based on multi-scale descriptor screening for mismatching includes the following steps:

step 1: three-dimensional point clouds X and Y of view-adjacent scenes are acquired, and a multi-scale descriptor (Q) is constructed for each point in the three-dimensional point clouds X and Y_i，N_i)；

In the embodiment, a 3D contour scanning sensor is adopted to obtain three-dimensional point cloud of a scene with an adjacent view angle;

wherein L represents a neighborhood search radius mark of a middle point of the three-dimensional point cloud X or Y, the values are 1, 2, …, L represents the total number of neighborhood search radii, and L is more than or equal to 2; t is_lWhich represents a normalized vector of the vector,

λ_l1，λ_l2，λ_l3is a matrix C_lEigenvalues after SVD decomposition, lambda_l1≥λ_l2≥λ_l3，n_l1，n_l2，n_l3Are each lambda_l1，λ_l2，λ_l3The corresponding feature vector is used as a basis for determining the feature vector,

S_l＝{p_i|||p_i-p₀||≤r_l}，|S_li represents S_lMaximum value of p_iRepresenting a point p in a three-dimensional point cloud X or Y₀Neighborhood point of r_lRepresents a point p₀The l-th neighborhood search radius of (c),

r_Lfor a set maximum neighborhood search radius, n_lRespectively representing points in a three-dimensional point cloudThe normal vector corresponding to the search neighborhood radius l of the point takes the value of n_l3；

In this example, L is 4 and r₄＝0.004m。

|S_lI is expressed as p₀Is the center of a sphere r_lThe maximum value of the distances from all points in the sphere to the sphere center;

step 2: traversing the multi-scale descriptors of each point in the three-dimensional point cloud X and Y, and respectively recording the points with the same multi-scale descriptors in the X and the Y as a key point set A_i0And B_i0；

And step 3: searching key point set A according to principal curvature characteristics_i0And B_i0Initial matching point set A in (1)_i1And B_i1；

Initial matching point set A_i1Point a in_i1And B_i1B in (1)_i1Is a matching point pair;

from the set of keypoints A_i0Selecting an arbitrary point a_iGo through B_i0Midpoint, select and point a_iMatched point b_jAn initial matching point set A is constructed with points satisfying the following conditions_i1And B_i1：

Wherein k is₁(a_i) And k₂(a_i) Respectively, at the point a_iAt point a_iThe principal curvature minimum and the principal curvature maximum, k, obtained from the neighborhood points of₁(b_j) And k₂(b_j) Respectively, at point b_jTo utilize point b_jThe principal curvature minimum and the principal curvature maximum, δ, obtained from the neighborhood points of₁Is a set first error threshold;

first error threshold δ₁The size of the obtained value is obtained by experimental simulation to obtain a best effect value;

the calculation process of the main curvature of each point in the three-dimensional point cloud is as follows:

step 3.1: fitting any point in the three-dimensional point cloud and a neighborhood point of the point to form a local quadric surface G (x, y);

G(x,y)＝ax²+bxy+cy²+dx+ey

step 3.2: performing derivation on the local quadric G (x, y) to obtain a first basic constant E, F, G and a second basic constant L, M, N;

E＝G_xx，F＝G_x*G_y，G＝G_yy；

L＝n·G_xx，M＝n·G_xy，N＝n·G_yy；

wherein G is_x、G_yEach G (x, y) is derived once for x and y, G_xx、G_yySecond derivative of x, y for G (x, y), G_xySequentially differentiating x and y for G (x, y), wherein n is a normal vector of the local quadric surface G (x, y) at any selected point;

step 3.3: calculating gaussian curvature K and mean curvature H:

step 3.4: calculating a principal curvature minimum k₁And a maximum value k of principal curvature₂：

And 4, step 4: initial matching point set A using similar triangle threshold_i1And B_i1Screening the matched point pairs to obtain a registration point pair set A_i2And B_i2；

The initial matching point set A is subjected to similar triangle threshold value pair_i1And B_i1The process of screening the matching point pairs in (1) is as follows:

step 4.1:in A_i1In the method, 3 points a which are not collinear are randomly selected_i,a_j,a_kIn B_i1Find 3 points b corresponding to it_u,b_v,b_w；

Step 4.2: order (a)_i,b_u) And (a)_j,b_v) Is a known pair of matching points, the following calculations are made according to the triangle similarity principle:

△a_ij＝||a_i-a_j||；△a_ik＝||a_i-a_k||；△a_jk＝||a_j-a_k||；

△b_uv＝||b_u-b_v||；△b_uw＝||b_u-b_w||；△b_wv＝||b_w-b_v||；

step 4.3: calculating the similarity omega of the two triangles according to a triangle similarity principle;

wherein, Δ a_ij,△a_ik,△a_jk；△b_uv，△b_uw,△b_wvFor each side length, delta, of a triangle₂Is a set second error threshold;

in this example, the δ₂The value is any one of 0.02, 0.05 and 0.1.

Step 4.4: if a is_i，a_j，a_k，b_u，b_v,b_wSatisfy omega<δ₂Then, consider a_kAnd b_wBelonging to a pair of registration point pairs.

And 5: computing A based on SVD_i2And B_i2First rotation matrix R in between₁And a first translation matrix T₁Using a first rotation matrix R₁And a first translation matrix T₁Sequentially calculating A_i2Each point a in_i2Change point a of_3iWhile calculating A_i2Each point in B_i2Matching point b in (1)_2iEuclidean distance d from corresponding transformation point_3i；

A_i3＝R₁*A_i2+T₁

Wherein (a)_3ix，a_3iy，a_3iz) And (b)_i2x,b_i2y,b_i2z) Are respectively a_3iAnd b_i2Coordinates on the x, y, z axes;

step 6: computing

k is 1 as the initial value, and N is the current registration point pair set A_i2And B_i2Centering the registration point pair, and judging d_k+1-d_k<ξ₃If yes, using the current registration point pair set A_i2And B_i2And (3) as registration data of the three-dimensional point clouds X and Y, if not, making k equal to k +1, returning to the step 3, reselecting the neighborhood points of each point, and calculating the principal curvature.

In this example, xi₃The value is 0.01.

Using a rotation angle threshold theta_mScreening the registration pairs again, and setting the current registration point pair A_i2And B_i2And (3) rejecting the registration point corresponding to the point which does not meet the following formula as an abnormal point or a noise point:

Δθ_i<θ_m

wherein, Delta theta_iRepresenting the current registration point pair set a_i2And B_i2The rotation angle difference between the ith pair of registration points;

current registration point pair set a_i2And B_i2The ith pair of registration point pairs (x)_i,y_i) The solving process of the rotation angle difference between the two is as follows:

1) set A from current registration point pair_i2And B_i2In the method, two pairs of non-coplanar registration pairs are randomly selected, four non-coplanar points form a spherical surface, the spherical surface equation and the spherical center O are solved by analytic geometry, the point O is used as an origin, and a space coordinate system is established, as shown in fig. 2(a), and in fig. 2(a), the y is used_iThe sphere at the center of the sphere has a common area of intersection with the sphere O, which is considered to be y_iRegions that meet registration criteria;

as long as x_iThrough R₁、T₁Is attributed to the point of_iIs the center of the sphere and xi is the radius, then x is considered_iAnd y_iIs a pair of registration pairs, i.e. translating the satisfaction of the threshold xi into a pair angle epsilon_iAs shown in fig. 2 (b);

considering that based on a screening condition (rotation transformation), if the threshold of the included angle is satisfied when rotation and translation transformation are performed, the threshold of a rotation angle based on rotation transformation is necessarily satisfied first, and the problem of maximum consistency is converted into: angle (R)₁x_i,y_i)≤ε_iWherein, epsilon_iJudging whether the point pair is an angle threshold value of the registration point pair;

from fig. 2(b), it can be derived that:

based on the above analysis, the following transformations are made:

r₁Decomposed into two matrices A, B, i.e. R₁A ═ a × B; (so that B satisfies the condition:angle (Bx)_i,y_i)≤ε_i(ii) a While making A as Bx_iThe rotation axis of (a) converts the maximum consistency problem into a rotation angle that makes a satisfy the condition as much as possible;

② in the absence of noise and outliers, x_iAnd y_iCan be spliced together perfectly. Defining an ideal R₂＝A₂B₂. By a pair of registration pairs x_k,y_kThen, B can be obtained₂(B₂x_k＝y_k)；

Using points i ≠ k, denoted B₂x_kIs a rotating shaft (i.e. A in FIG. 2 (c))₂Rotational axis of) which is found so that B₂x_iAnd y_iRotation angle theta of perfect fit₂(because of the point B₂x_iAnd y_iThe coordinates in the spatial coordinate system are known, A₂The rotation axis of (2) can be regarded as the Z axis in the space coordinate system, and the rotation angle theta can be directly obtained by analyzing the geometrical relationship₂That is to obtain A₂：

Wherein, [ X ]]_XRepresents the vector product of X, exp () represents the exponential mapping relationship.

③ but in general noise and outliers are unavoidable, so R is used₂To estimate x_kA range limit after rotation; to make < (Bx)_k,y_k)≤ε_iKnowing Bx_kIs in a shaded circular area in fig. 2 (c). Similarly, for points where i ≠ k, we can obtain Bx_iAlso in a shaded circular area in figure (c).

Is recorded as:

2) taking the Z axis of the established space coordinate system as a second rotation matrix R₂And R is a rotational axis of₂＝A₂B₂，A₂And B₂Respectively represent second rotation matrices R₂Decomposition matrix of, B₂Set A by current registration point pair_i2And B_i2One pair of registration point pairs (x)_k,y_k) Satisfies B₂x_k＝y_kThen, a calculation is performed to obtain₂The rotation axis of the system is Z axis, and simultaneously, according to the set registration distance threshold value xi, the registration angle threshold value epsilon is calculated_i，

3) Rotate by

A circular area enclosed when

Left boundary and y of the enclosed circular area_iWhen tangent (as shown in FIG. 2 (d)), is composed of

The rotation angle of the enclosed circular area is theta₁At this time, the product is made by

The center of the circle of the circular area is the point M which is formed by the points M and y_iA distance r of_sConstructing a one-dimensional quadratic equation, using the equation of the point M in the great circle O, simultaneously establishing two equation sets to obtain the x and y coordinate values of the point M, and solving theta₁；

Represents point x_iThrough B₂Rigid body transformation to form B₂x_iIs the center of a sphere, epsilon_iA circular area of radius, as shown in fig. 2 (e);

because B₂Is a rigid body transformation rotation matrix in the presence of noise or outliers,

it is shown that: x is the number of_iThrough B₂Rigid body transformation, but because of noise or abnormal points, only one satisfactory range can be found, and the range is B₂x_iIs the center of a sphere, epsilon_iA circular area of radius;

the great circle O is pointing x_iIn the space coordinate system, a circle with a point O as an origin is located;

4) continue to rotate

A circular area enclosed by

Right boundary and y of the enclosed circular area_iWhen the two-way pipe is tangent to each other,

the rotation angle of the enclosed circular area is theta₂At this time, the product is made by

The circle center of the circular area is point N, and the point N and the point y are_iA distance r of_sConstructing a one-dimensional quadratic equation, using the equation of N in the great circle O, simultaneously establishing two equation sets to obtain the x and y coordinate values of the point N, and solving theta₂To obtain an angle difference Delta theta_i；

Wherein r is_sIs that

Radius of a circular area enclosed, r, being B in a space coordinate system O₂x_iDistance to point O, r_s＝г*tanε_i。

The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments or alternatives may be employed by those skilled in the art without departing from the spirit or ambit of the invention as defined in the appended claims.

Claims (6)

1. A three-dimensional image registration method for screening out mismatching based on multi-scale descriptors is characterized by comprising the following steps:

step 1: three-dimensional point clouds X and Y of view-adjacent scenes are acquired, and a multi-scale descriptor (Q) is constructed for each point in the three-dimensional point clouds X and Y_i，N_i)；

Wherein L represents a neighborhood search radius mark of a middle point of the three-dimensional point cloud X or Y, the values are 1, 2, …, L represents the total number of neighborhood search radii, and L is more than or equal to 2; t is_lWhich represents a normalized vector of the vector,

λ_l1，λ_l2，λ_l3is a matrix C_lEigenvalues after SVD decomposition, lambda_l1≥λ_l2≥λ_l3，n_l1，n_l2，n_l3Are each lambda_l1，λ_l2，λ_l3The corresponding feature vector is used as a basis for determining the feature vector,

S_l＝{p_i|||p_i-p₀||≤r_l}，|S_li represents S_lMaximum value of p_iRepresenting a point p in a three-dimensional point cloud X or Y₀Neighborhood point of r_lRepresents a point p₀The l-th neighborhood search radius of (c),

r_Lfor a set maximum neighborhood search radius, n_lRespectively representing normal vectors corresponding to the search neighborhood radius l of the points in the three-dimensional point cloud, and taking the value as n_l3；

Step 2: traversing the multi-scale descriptors of each point in the three-dimensional point cloud X and Y, and respectively recording the points with the same multi-scale descriptors in the X and the Y as a key point set A_i0And B_i0；

And step 3: searching key point set A according to principal curvature characteristics_i0And B_i0Initial matching point set A in (1)_i1And B_i1；

From the set of keypoints A_i0Selecting an arbitrary point a_iGo through B_i0Midpoint, select and point a_iMatched point b_jAn initial matching point set A is constructed with points satisfying the following conditions_i1And B_i1：

Wherein k is₁(a_i) And k₂(a_i) Respectively, at the point a_iAt point a_iThe principal curvature minimum and the principal curvature maximum, k, obtained from the neighborhood points of₁(b_j) And k₂(b_j) Respectively, at point b_jTo utilize point b_jThe principal curvature minimum and the principal curvature maximum, δ, obtained from the neighborhood points of₁Is a set first error threshold;

and 4, step 4: initial matching point set A using similar triangle threshold_i1And B_i1Screening the matched point pairs to obtain a registration point pair set A_i2And B_i2；

And 5: computing A based on SVD_i2And B_i2First rotation matrix R in between₁And a first translation matrix T₁Using a first rotation matrix R₁And a first translation matrix T₁Sequentially calculating A_i2Each point a in_i2Change point a of_3iWhile calculating A_i2Each point in B_i2Matching point b in (1)_2iEuclidean distance d from corresponding transformation point_3i；

A_i3＝R₁*A_i2+T₁

Wherein (a)_3ix,a_3iy,a_3iz) And (b)_i2x,b_i2y,b_i2z) Are respectively a_3iAnd b_i2Coordinates on the x, y, z axes;

step 6: computing

k is 1 as the initial value, and N is the current registration point pair set A_i2And B_i2Centering the registration point pair, and judging d_k+1-d_k<ξ₃Whether it is true, xi₃Indicating an error threshold greater than zero, and if true, set a with the current registration point pair_i2And B_i2As registration data of the three-dimensional point clouds X and Y, if not, making k equal to k +1, returning to the step 3, reselecting neighborhood points of each point, and calculating the principal curvature;

using a rotation angle threshold theta_mScreening the registration pairs again, and setting the current registration point pair A_i2And B_i2And (3) rejecting the registration point corresponding to the point which does not meet the following formula as an abnormal point or a noise point:

Δθ_i<θ_m

wherein, Delta theta_iRepresenting the current registration point pair set a_i2And B_i2The rotation angle difference between the ith pair of registration points;

current registration point pair set a_i2And B_i2The ith pair of registration point pairs (x)_i,y_i) The solving process of the rotation angle difference between the two is as follows:

1) set A from current registration point pair_i2And B_i2Two pairs of non-coplanar registration pairs are randomly selected, four non-coplanar points form a spherical surface, the spherical surface equation and the spherical center O are solved by analytic geometry, and a space coordinate system is established by taking the point O as an origin;

2) taking the Z axis of the established space coordinate system as a second rotation matrix R₂And R is a rotational axis of₂＝A₂B₂，A₂And B₂Respectively represent second rotation matrices R₂Decomposition matrix of, B₂Set A by current registration point pair_i2And B_i2One pair of registration point pairs (x)_k,y_k) Satisfies B₂x_k＝y_kThen, a calculation is performed to obtain₂The rotation axis of the system is Z axis, and simultaneously, according to the set registration distance threshold value xi, the registration angle threshold value epsilon is calculated_i，

3) Rotate by

A circular area enclosed when

Left boundary and y of the enclosed circular area_iWhen tangent, by

Circular shape formed by encirclingThe rotation angle of the region is theta₁At this time, the product is made by

The center of the circle of the circular area is the point M which is formed by the points M and y_iA distance r of_sConstructing a one-dimensional quadratic equation, using the equation of the point M in the great circle O, simultaneously establishing two equation sets to obtain the x and y coordinate values of the point M, and solving theta₁；

Represents point x_iThrough B₂Rigid body transformation to form B₂x_iIs the center of a sphere, epsilon_iA circular area of radius;

the great circle O is pointing x_iIn the space coordinate system, a circle with a point O as an origin is located;

4) continue to rotate

A circular area enclosed by

Right boundary and y of the enclosed circular area_iWhen the two-way pipe is tangent to each other,

the rotation angle of the enclosed circular area is theta₂At this time, the product is made by

The circle center of the circular area is point N, and the point N and the point y are_iA distance r of_sConstructing a one-dimensional quadratic equation, using the equation of N in the great circle O, simultaneously establishing two equation sets to obtain the x and y coordinate values of the point N, and solving theta₂To obtain an angle difference Delta theta_i；

Wherein r is_sIs that

Radius of a circular area enclosed, r, being B in a space coordinate system O₂x_iDistance to point O, r_s＝г*tanε_i。

2. The method of claim 1, wherein the principal curvature of each point in the three-dimensional point cloud is calculated as follows:

step 3.1: fitting any point in the three-dimensional point cloud and a neighborhood point of the point to form a local quadric surface G (x, y);

G(x,y)＝ax²+bxy+cy²+dx+ey

step 3.2: performing derivation on the local quadric G (x, y) to obtain a first basic constant E, F, G and a second basic constant L, M, N;

E＝G_xx，F＝G_x*G_y，G＝G_yy；

L＝n·G_xx，M＝n·G_xy，N＝n·G_yy；

wherein G is_x、G_yEach G (x, y) is derived once for x and y, G_xx、G_yySecond derivative of x, y for G (x, y), G_xySequentially differentiating x and y for G (x, y), wherein n is a normal vector of the local quadric surface G (x, y) at any selected point;

step 3.3: calculating gaussian curvature K and mean curvature H:

step 3.4: calculating a principal curvature minimum k₁And a maximum value k of principal curvature₂：

3. The method of claim 1, wherein the initial matching point set A is thresholded using similar triangles_i1And B_i1The process of screening the matching point pairs in (1) is as follows:

step 4.1: in A_i1In the method, 3 points a which are not collinear are randomly selected_i,a_j,a_kIn B_i1Find 3 points b corresponding to it_u,b_v,b_w；

Step 4.2: order (a)_i,b_u) And (a)_j,b_v) Is a known pair of matching points, the following calculations are made according to the triangle similarity principle:

△a_ij＝||a_i-a_j||；△a_ik＝||a_i-a_k||；△a_jk＝||a_j-a_k||；

△b_uv＝||b_u-b_v||；△b_uw＝||b_u-b_w||；△b_wv＝||b_w-b_v||；

step 4.3: calculating the similarity omega of the two triangles according to a triangle similarity principle;

wherein, Δ a_ij,△a_ik,△a_jk；△b_uv,△b_uw,△b_wvFor each side length, delta, of a triangle₂Is a set second error threshold;

step 4.4: if a is_i,a_j，a_k，b_u，b_v,b_wSatisfy omega<δ₂Then, consider a_kAnd b_wBelonging to a pair of registration point pairs.

4. The method of any one of claims 1-3, wherein L is 4 and r is₄＝0.004m。

5. The method of claim 3, wherein δ₂The value is any one of 0.02, 0.05 and 0.1.

6. The method of claim 1, wherein a 3D contour scanning sensor is used to acquire a three-dimensional point cloud of view-adjacent scenes.

Images