CN110458174A - A kind of unordered accurate extracting method of cloud key feature points - Google Patents

Abstract

A kind of unordered accurate extracting method of cloud key feature points, it includes have the following steps: 1) centered on each cloud point, calculating a ring k neighborhood union and construct triangle meshes；2) discrete curvature, discrete normal and the neighborhood point distance for calculating each cloud point are used to construct window function, calculate gradient matrix；3) angle point response characteristic of division point is calculated to gradient matrix.The present invention remains local detail feature using angle point, folding line point and planar point in change of gradient identification point cloud data, and recognition result is accurate, reliable；The 3 D laser scanning detection that can be widely applied to the components such as aviation, automobile, has very strong applicability.

Description

A kind of unordered accurate extracting method of cloud key feature points

Technical field

The present invention relates to computer graphics and three-dimensional point cloud feature extraction field, and in particular to a kind of unordered point Yun Guanjian The accurate extracting method of characteristic point.

Background technique

Airplane component has a large scale, weak rigidity and the characteristics of structure is complicated, not only needs to measure the sky of single component Between size, and the assembly features of complex component are also required to measure.The eyeball cloud amount of these structural members is big and unordered, quasi- Really extracting key feature points is the necessary means simplified data, improve Data Analysis Services efficiency.Feature extraction is metrical information Most basic in processing, most critical part, reverse-engineering, the precision of image procossing depend critically upon the precision of feature extraction.

Extracting field in cloud key feature both at home and abroad has had certain research achievement, mainly uses three kinds of modes: First is that the feature extraction based on normal vector, by combining clustering algorithm determining method to swear, angle is more than that the point of threshold value is characterized a little, this Kind method is suitable for the part of simple rule shape, but lower for the accuracy of identification of local fine feature.Second is that being based on curvature Feature extraction, identified according to curvature information mutation at characteristic point.Although this method has noise spot relatively strong Resistance, but extract characteristic point inside can have pseudo feature point, need additional filtering algorithm to be simplified.Third is that being based on certain Kind mathematical model or operator improve, and this method improves Feature point recognition precision by designing specific solution process, but The problem of being omitted there are still wrong identification and characteristic point.

Summary of the invention

The purpose of the present invention is the accuracy of identification of existing local fine feature is extracted for existing cloud key feature It is lower and there are pseudo feature point and characteristic points to omit, the problem of needing additional filtering algorithm to be simplified, provide a kind of nothing The accurate extracting method of sequence point cloud key feature points, it is therefore intended that avoid the non-characteristic point in 2 ring neighborhood of real features point by mistake The problem of identification, accurately identifies angle point, planar point and folding line point using change of gradient and guarantees local detail feature.

The technical scheme is that

Unordered cloud key feature of one kind accurately identifies and extracting method, it is characterized in that carrying out first to each cloud point adjacent Domain search, calculates neighborhood distance and angle sorts counterclockwise to neighborhood point；Then triangle gridding is constructed to a ring k neighborhood, The discrete normal direction for calculating vertex calculates discrete curvature using Taylor expansion；Finally, to neighborhood distance, discrete normal direction and discrete song Rate is weighted construction gradient matrix and calculates angle point response, passes through given threshold identification feature point.Specific step is as follows institute It states:

Step 1 calculates each cloud point P_iA ring k neighborhood P_i-k。

(1) neighborhood k value value range [10,20] are set, each data point is found apart from k+5 nearest neighborhood point, right Neighborhood point is ranked up counterclockwise, calculates neighbouring 2 neighborhood points and point cloud point P_iThe vector angle θ of composition；

(2) calculating the smallest neighborhood point of space length is starting point, calculates threshold angle θ₀=360/k, along counterclockwise Direction is less than threshold angle θ when the angle theta of 2 neighborhood points₀, the long biggish neighborhood point of mould is deleted, another point is labeled as qualified point Retained, qualification point is as new starting point.Circulation (2) terminates when neighborhood point is reduced to k.

Step 2, the neighborhood relationships obtained based on step 1 construct triangle meshes to calculate each cloud point P_iIt is discrete Normal vector calculates discrete curvature using Taylor expansion, specifically:

(1) as shown in Figure 1, calculating all neighborhood triangle side length a/b/c, the shape weight of any neighborhood triangle is calculated Δ_j, formula is as follows:

(2) normal vector N is calculated according to neighborhood triangular apex coordinate value_fj, according to the shape weight that (1) obtains, calculate point The discrete normal direction N of cloud point Pi_p。

(3) as shown in Fig. 2, calculating point cloud point P according to the normal vector that (2) obtain_iTangent plane F, to P_iWith P_i-kPlace Area of space carry out parametric surface S (u, v) fitting.Construct P_iAs the local coordinate system xyz of origin, surface coordinates system u/v It projects on the face F and corresponds to x/y axis.Calculate neighborhood point P_i-kProject to the subpoint Q of tangent plane F_j-k, calculate any subpoint Q_j-k With x-axis angle σ_j, Q_j-kWith y-axis angle η_j。

(4) Taylor expansion is carried out to curved surface S, calculation formula is as follows:

R_nIt is a higher-order shear deformation, this value is ignored in subsequent calculating.It is write above formula as matrix form, calculates each data point P_iDiscrete single order { the S of point_u,S_vAnd second order { S_uu,S_uv,S_vvDerivative, calculation formula is as follows:

Wherein, D_j=| P_i-k-P_i|, Δ u_j=| Q_j-k-P_i|cosσ_j, Δ v_j=| Q_j-k-P_i|cosη_j。

(5) the discrete first derivative obtained according to (4) calculates outer En Jiateng transformation matrix W, and formula is as follows:

In formula: E, F, G are the first kind fundamental quantity of curved surface, L, N, M₁、M₂For the second class fundamental quantity of curved surface.

(6) the characteristic value k of outer En Jiateng transformation matrix W is calculated₁,k₂, obtain point cloud point P_iDiscrete Gaussian curvature κ, formula It is as follows:

κ=k₁·k₂

Step 3 is weighted construction gradient matrix to neighborhood distance, discrete normal direction and discrete curvature and calculates characteristic value, It is specific as follows by given threshold identification feature point:

(1) neighborhood point spacing d, the discrete normal direction obtained using step 2 and discrete Gaussian curvature calculation window function are calculated W, formula are as follows:

Wherein, σ_k,σ_n, and σ_cIt is bandwidth parameter, σ_k=max (d), σ_n=180, σ_c=max (κ_pi-κ_pi-k)。

(2) window function obtained by (1) calculates P_iAlong any P_i-kThe change of gradient of point vector direction, formula is such as Under

(3) by P_iSingle order and Second-Order Discrete derivative substitute into above formula, calculate P_iIntegral gradient E in a ring neighborhood, formula It is as follows:

Above formula is rewritten, neighborhood point P is calculated_i-kDirection cosines cos α, cos β, cos at local coordinate system xyz γ obtains P_iCovariance matrix M (the P of point gradient_i), calculation formula is as follows:

(4) as shown in figure 3, calculating gradient matrix M (P_i) three characteristic values and be ranked up, obtain λ₁≥λ₂≥λ₃If Determine threshold value Ω, made the following judgment according to size relation:

1) angle point: λ₁≥λ₂≥λ₃≥Ω。

2) seamed edge point: λ₁≥λ₂, λ₁>=Ω, λ₃≠0；Or λ₂≥λ₁, λ₂>=Ω, λ₃≠0。

3) planar point: λ₁≤ Ω, λ₂≤ Ω, λ₃≤Ω。

The curved surface S (u, v) needs to carry out least square and solves acquisition parametric equation, and detailed process is as follows:

(1) surface equation general expression is denoted as:

Z=a₁x²+a₂xy+a₃y²+a₄x+a₅y+a₆

(2) by P_iPoint and any 5 neighborhoods point P_i-kIt brings into above formula respectively, obtains hexa-atomic linear function group and asked Solution obtains coefficient sets A={ a1 ... a6 }.Step (2) are repeated to carry out being no more than N=C_k ⁵Secondary permutation and combination obtains N group coefficient Collection.

(3) all neighborhood point P are calculated_i-kTo the sum of the Z-direction projector distance of every group of coefficient fitting surface, minimum value institute is selected Corresponding surface coefficients are as fitting result.

Simplify the deterministic process of characteristic value size relation, detailed process by calculating angle point receptance function R are as follows:

(1) to gradient matrix E (P_i) derivation is carried out, angle point response R is obtained, calculation formula is as follows:

(2) it is [0.01,0.1] according to the value range of Ω, judges the size relation of angle point analog value R and threshold value Ω such as Under:

1) work as R > Ω: P_iPoint angle point.

2) when R <-Ω: P_iPoint is folding line point.

3) as | R | < Ω: P_iPoint is planar point.

The beneficial effects of the present invention are:

One, the present invention carries out the extraction of characteristic point using the variation of gradient, as shown in figure 4, the other methods that compare are in office There is higher accuracy of identification, local detail Feature point recognition rate is 97% or more in portion's details.

Two, for scale, sagging and chamfering category feature point, (feature shows as a plurality of contour line to the present invention, and between distance Every the minimum accuracy of identification up to 0.5mm) in 0.15mm (1.5 times of eyeball cloud equalization point spacing) hereinafter, can be than calibrated The non-characteristic point between neighbouring two profiles is really rejected, without additional filtering algorithm.

Detailed description of the invention

Mono- ring neighborhood triangle meshes of Fig. 1 construct schematic diagram.

The building of mono- ring neighborhood local coordinate system of Fig. 2 and neighborhood point perspective view.

Fig. 3 gradient matrix extracts characteristic point schematic diagram.

Fig. 4 the method for the present invention figure compared with other feature extraction algorithm performances.

Fig. 5 the method for the present invention operational flowchart.

Specific embodiment

The present invention is further illustrated for Structure Figure and embodiment below.

As shown in Figs. 1-5.

The neighborhood triangle gridding topology that the present invention is applicable in please refers to shown in Fig. 1, the building process of local coordinate system and one Ring neighborhood point projecting direction referring to figure 2., calculates the characteristic value that gradient matrix obtains and please refers to figure to extract the process of characteristic point 3。

Unordered cloud key feature of one kind accurately identifies and extracting method, as shown in figure 5, it includes the following steps.

Step 1 calculates each cloud point P_iA ring k neighborhood P_i-k。

(1) neighborhood k value value range [10,20] are set, each data point is found apart from k+5 nearest neighborhood point, right Neighborhood point is ranked up counterclockwise, calculates neighbouring 2 neighborhood points and point cloud point P_iThe vector angle θ of composition；

(2) calculating the smallest neighborhood point of space length is starting point, calculates threshold angle θ₀=360/k, along counterclockwise Direction is less than threshold angle θ when the angle theta of 2 neighborhood points₀, the long biggish neighborhood point of mould is deleted, another point is labeled as qualified point Retained, qualification point is as new starting point.Circulation (2) terminates when neighborhood point is reduced to k.

Step 2, the neighborhood relationships obtained based on step 1 construct triangle meshes to calculate each cloud point P_iIt is discrete Normal vector calculates discrete curvature using Taylor expansion, specifically:

(1) as shown in Figure 1, calculating all neighborhood triangle side length a/b/c, the shape weight of any neighborhood triangle is calculated Δ_j, formula is as follows:

(2) normal vector N is calculated according to neighborhood triangular apex coordinate value_fj, according to the shape weight that (1) obtains, calculate point The discrete normal direction N of cloud point Pi_p。

(3) as shown in Fig. 2, calculating point cloud point P according to the normal vector that (2) obtain_iTangent plane F, to P_iWith P_i-kPlace Area of space carry out parametric surface S (u, v) fitting.

The curved surface S (u, v) can be solved by least square and be obtained parametric equation, and detailed process is as follows:

(1) surface equation general expression is denoted as:

Z=a₁x²+a₂xy+a₃y²+a₄x+a₅y+a₆

(2) by P_iPoint and any 5 neighborhoods point P_i-kIt brings into above formula respectively, obtains hexa-atomic linear function group and asked Solution obtains coefficient sets A={ a1 ... a6 }.Step (2) are repeated to carry out being no more than N=C_k ⁵Secondary permutation and combination obtains N group coefficient Collection.

(3) all neighborhood point P are calculated_i-kTo the sum of the Z-direction projector distance of every group of coefficient fitting surface, minimum value institute is selected Corresponding surface coefficients are as fitting result.

Construct P_iAs the local coordinate system xyz of origin, surface coordinates system u/v, which is projected on the face F, corresponds to x/y axis.It calculates Neighborhood point P_i-kProject to the subpoint Q of tangent plane F_j-k, calculate any subpoint Q_j-kWith x-axis angle σ_j, Q_j-kWith y-axis angle η_j。

(4) Taylor expansion is carried out to curved surface S, calculation formula is as follows:

R_nIt is a higher-order shear deformation, this value is ignored in subsequent calculating.It is write above formula as matrix form, calculates each data point P_iDiscrete single order { the S of point_u,S_vAnd second order { S_uu,S_uv,S_vvDerivative, calculation formula is as follows:

Wherein, D_j=| P_i-k-P_i|, Δ u_j=| Q_j-k-P_i|cosσ_j, Δ v_j=| Q_j-k-P_i|cosη_j。

(5) the discrete first derivative obtained according to (4) calculates outer En Jiateng transformation matrix W, and formula is as follows:

(6) the characteristic value k of outer En Jiateng transformation matrix W is calculated₁,k₂, obtain point cloud point P_iDiscrete Gaussian curvature κ, formula It is as follows:

κ=k₁·k₂

Step 3 is weighted construction gradient matrix to neighborhood distance, discrete normal direction and discrete curvature and calculates characteristic value, It is specific as follows by given threshold identification feature point:

(1) neighborhood point spacing d, the discrete normal direction obtained using step 2 and discrete Gaussian curvature calculation window function are calculated W, formula are as follows:

Wherein, σ_k,σ_n, and σ_cIt is bandwidth parameter, σ_k=max (d), σ_n=180, σ_c=max (κ_pi-κ_pi-k)。

(2) window function obtained by (1) calculates P_iAlong any P_i-kThe change of gradient of point vector direction, formula is such as Under

(3) by P_iSingle order and Second-Order Discrete derivative substitute into above formula, calculate P_iIntegral gradient E in a ring neighborhood, formula It is as follows:

Above formula is rewritten, neighborhood point P is calculated_i-kDirection cosines cos α, cos β, cos at local coordinate system xyz γ obtains P_iCovariance matrix M (the P of point gradient_i), calculation formula is as follows:

(4) as shown in figure 3, calculating gradient matrix M (P_i) three characteristic values and be ranked up, obtain λ₁≥λ₂≥λ₃If Determine threshold value Ω, made the following judgment according to size relation:

1) angle point: λ₁≥λ₂≥λ₃≥Ω。

2) seamed edge point: λ₁≥λ₂, λ₁>=Ω, λ₃≠0；Or λ₂≥λ₁, λ₂>=Ω, λ₃≠0。

3) planar point: λ₁≤ Ω, λ₂≤ Ω, λ₃≤Ω。

The deterministic process of characteristic value size relation, tool can also be simplified when specific implementation by calculating angle point receptance function R Body process are as follows:

(1) to gradient matrix E (P_i) derivation is carried out, angle point response R is obtained, calculation formula is as follows:

(2) it is [0.01,0.1] according to the value range of Ω, judges the size relation of angle point analog value R and threshold value Ω such as Under:

1) work as R > Ω: P_iPoint angle point.

2) when R <-Ω: P_iPoint is folding line point.

3) as | R | < Ω: P_iPoint is planar point.

Part that the present invention does not relate to is same as the prior art or can be realized by using the prior art.

Claims (4)

1. a kind of unordered cloud key feature accurately identifies and extracting method, it is characterized in that carrying out neighborhood to each cloud point first Search, while calculating neighborhood distance and angle and being sorted counterclockwise to neighborhood point；Then the triangulation network is constructed to a ring k neighborhood Lattice calculate the discrete normal direction on vertex, calculate discrete curvature using Taylor expansion；Finally, to neighborhood distance, discrete normal direction and discrete Curvature is weighted construction gradient matrix and calculates angle point response, passes through given threshold identification feature point.

2. according to the method described in claim 1, it is characterized in that it the following steps are included:

Step 1 calculates each cloud point P_iA ring k neighborhood P_i-k。

(1) neighborhood k value value range [10,20] are set, each data point is found apart from k+5 nearest neighborhood point, to neighborhood Point is ranked up counterclockwise, calculates neighbouring 2 neighborhood points and point cloud point P_iThe vector angle θ of composition；

(2) calculating the smallest neighborhood point of space length is starting point, calculates threshold angle θ₀=360/k, in the counterclockwise direction when The angle theta of 2 neighborhood points is less than threshold angle θ₀, the long biggish neighborhood point of mould is deleted, another point is protected labeled as qualified point It stays, qualification point is as new starting point；Circulation (2) terminates when neighborhood point is reduced to k；

Step 2, the neighborhood relationships obtained based on step 1 construct triangle meshes to calculate each cloud point P_iDiscrete normal direction Amount calculates discrete curvature using Taylor expansion, specifically:

(1) all neighborhood triangle side length a/b/c are calculated, the shape weight Δ of any neighborhood triangle is calculated_j, formula is as follows:

(2) the shape weight obtained according to (1) calculates point cloud point P_iDiscrete normal direction N_p；

In formula: N_fjFor the normal vector of neighborhood triangle, calculates the vector multiplication cross that any two vertex of neighborhood triangle is constituted and obtain.

(3) point cloud point P was calculated according to the normal vector that (2) obtain_iTangent plane F, to P_iWith P_i-kThe area of space at place carries out Parametric surface S (u, v) fitting, constructs P_iAs the local coordinate system xyz of origin, the direction selection of x-axis is that most short mould is long | P_i- P_i-k| vector direction.

(4) Taylor expansion is carried out to curved surface S, calculation formula is as follows:

U/v is the bent line coordinate direction of curved surface in formula, projects on the face F and corresponds to x/y axis.R_nIt is a higher-order shear deformation, subsequent meter This value is ignored in calculation.It is write above formula as matrix form, calculates each data point P_iDiscrete single order { the S of point_u,S_vAnd second order { S_uu, S_uv,S_vvDerivative, calculation formula is as follows:

Wherein, D_j=| P_i-k-P_i|, Δ u_j=| Q_j-k-P_i|cosσ_j, Δ v_j=| Q_j-k-P_i|cosη_j；Q in formula_j-kFor neighborhood point P_i-k Project to the subpoint of tangent plane F, σ_jFor any subpoint Q_j-kWith x-axis angle, η_jFor any subpoint Q_j-kWith y-axis angle.

(5) the discrete first derivative obtained according to (4) calculates outer En Jiateng transformation matrix W, and formula is as follows:

In formula: E, F, G are the first kind fundamental quantity of curved surface, L, N, M₁、M₂For the second class fundamental quantity of curved surface.

(6) the characteristic value k of outer En Jiateng transformation matrix W is calculated₁,k₂, obtain point cloud point P_iDiscrete Gaussian curvature κ, formula is such as Under:

κ=k₁·k₂

Step 3 is weighted construction gradient matrix to neighborhood distance, discrete normal direction and discrete curvature and calculates characteristic value, passes through Given threshold identification feature point, specific as follows:

(1) neighborhood point spacing d, the discrete normal direction obtained using step 2 and discrete Gaussian curvature calculation window function w are calculated, it is public Formula is as follows:

Wherein, σ_k,σ_n, and σ_cIt is bandwidth parameter, σ_k=max (d), σ_n=180, σ_c=max (κ_pi-κ_pi-k)。

(2) window function obtained by (1) calculates P_iAlong any P_i-kThe change of gradient of point vector direction, formula are as follows

(3) by P_iSingle order and Second-Order Discrete derivative substitute into above formula, calculate P_iIntegral gradient E in a ring neighborhood, formula is such as Under:

Above formula is rewritten, P is obtained_iCovariance matrix M (the P of point gradient_i), calculation formula is as follows:

Cos α in formula, cos β, cos γ respectively correspond neighborhood point P_i-kDirection cosines at local coordinate system x/y/z, Sx, Sy, Sz is the first-order partial derivative of curved surface S.

(4) gradient matrix M (P is calculated_i) three characteristic values and be ranked up, obtain λ₁≥λ₂≥λ₃, given threshold Ω, according to big Small relationship makes the following judgment:

1) angle point: λ₁≥λ₂≥λ₃≥Ω；

2) seamed edge point: λ₁≥λ₂, λ₁>=Ω, λ₃≠0；Or λ₂≥λ₁, λ₂>=Ω, λ₃≠0；

3) planar point: λ₁≤ Ω, λ₂≤ Ω, λ₃≤Ω。

3. according to the method described in claim 2, it is characterized in that the curved surface S (u, v) need to carry out least square solution obtain Obtain parametric equation:

(1) surface equation general expression is denoted as:

Z=a₁x²+a₂xy+a₃y²+a₄x+a₅y+a₆

(2) by P_iPoint and any 5 neighborhoods point P_i-kIt brings into above formula respectively, obtains hexa-atomic linear function group and solved, obtained To coefficient sets A={ a1 ... a6 }；Step (2) are repeated to carry out being no more than N=C_k ⁵Secondary permutation and combination obtains N group coefficient set；

(3) all neighborhood point P are calculated_i-kTo the sum of the Z-direction projector distance of every group of coefficient fitting surface, select corresponding to minimum value Surface coefficients as fitting result.

4. according to the method described in claim 2, it is characterized in that being closed by calculating angle point receptance function R to simplify characteristic value size The deterministic process of system:

(1) to gradient matrix E (P_i) derivation is carried out, angle point response R is obtained, calculation formula is as follows:

(2) it is [0.01,0.1] according to the value range of Ω, judges that angle point analog value R and the size relation of threshold value Ω are as follows:

1) work as R > Ω: P_iPoint angle point.

2) when R <-Ω: P_iPoint is folding line point.

3) as | R | < Ω: P_iPoint is planar point.