CN102799763A - Point cloud posture standardization-based method for extracting linear characteristic of point cloud - Google Patents

Abstract

The invention provides a point cloud posture standardization-based method for extracting the linear characteristic of point cloud, which comprises six steps, is used for extracting the linear characteristic in disordered and three-dimensional point cloud, can conveniently measure the relative posture of a target, and belongs to the technical field of the three-dimensional measurement and the machine vision. The method comprises the following steps of: firstly, building a KD-TREE structure of point cloud, so that the searching speed of the adjacent point set of the point cloud can be improved; secondly, building the adjacent point set of each point according to the density of the whole point cloud, obtaining the main direction of the point set, and building a Householder transformation matrix to adjust the posture of the point cloud; thirdly, carrying out surface fitting on the adjacent point set to obtain two main curvature of the point based on a curved surface equation, and selecting the main curvature with the higher absolute value as the curvature estimation of the point; and finally, obtaining the curvature estimation value of all point cloud, and taking the point which is larger than a given threshold value as a linear characteristic point, so that the linear characteristic can be extracted.

Description

A kind of based on standardized some cloud line of a cloud attitude feature extracting method

Technical field

The present invention relates to a kind of based on standardized some cloud line of a cloud attitude feature extracting method, it be based on that the impact point cloud carries out that targeted attitude is measured or the various visual angles coupling in extract a kind of method of some cloud line characteristic, the relative attitude that is applied to a cloud target is measured.Belong to three-dimensional measurement and machine vision technique field.

Technical background

The technology of obtaining of three dimensional point cloud is comparative maturity, and common method has based on binocular stereo vision to be obtained high-precision dot cloud information, obtain the some cloud information of object fast through laser scanning methods, and other three-dimensional point cloud obtains technology.Based on three dimensional point cloud, can realize processing dimension measurement, reverse-engineering, object pose measurement etc., in these technology, the gordian technique of bottom is feature extraction.The line characteristic has again that data volume is little, the expression structure characteristic is important, be easy to advantage and using values such as inventory analysis processing.

Feature extraction has several kinds of representational methods about the three-dimensional point cloud line, mainly contains following several types:

Carry out projection to three dimensions point cloud to specific plane, planar utilize the method for seeking limit or peak value to carry out feature point extraction.But the isoplanar as a result of plane projection is selected and the original point cloud structure is closely related, and the result of extraction is prone to receptor site cloud structure attitude and chooses the multiple influence in plane and stable inadequately.

As the unique point discriminant criterion, for the bigger unique point cloud of structural change, this angle can be bigger with point of proximity cloud vector angle for utilization point cloud vector, and this angle can be smaller for plane or smooth surface.Have stability with some cloud vector angle as the feature selecting foundation to spatial alternation, but bigger to a have calculated amount of cloud computing vector, and noise spot is to the obvious effect of result of calculation.

Use the wider feature extracting method that is based on curvature in addition, describe a smooth degree on cloud surface with curvature.The general line feature extraction algorithm based on curvature can be according to a cloud constructing curve, based on surface equation calculation level cloud Curvature Estimation value.But general algorithm selects still the space of improvement to be arranged in point set search, the adjustment of some cloud principal direction, Curvature Estimation mode.

Summary of the invention

Technical matters: the invention provides a kind of based on standardized some cloud line of a cloud attitude feature extracting method; It converges functions such as principal direction aligning selection principal curvatures through adding point; Improving a lot aspect characteristic mass and the algorithm stability based on curvature method extraction effect than existing, can be follow-up relative attitude measurement, field stitching unique point cloud preferably is provided.

Technical scheme: the targeted attitude based on cloud data is measured or field stitching; Be widely used in practice; Directly carry out analyzing and processing with scanning institute invocation point cloud, it is more to put cloud quantity on the one hand, speed and complexity that influence is handled; On the other hand because the noise spot influence that a large amount of visual fields difference causes makes final measuring accuracy be difficult to ensure.Therefore propose cloud data is carried out the thought of line feature extraction, used the measurement or the identification of the less more convenient target of cloud data.

Corresponding relation according to a cloud surface data mechanism characteristics and curvature; The present invention proposes a kind of based on standardized some cloud line of a cloud attitude feature extracting method; Be applicable to unordered three-dimensional point cloud center line Feature Extraction; This method at first makes up the KD-TRE E structure of a cloud, to improve the search speed that the some cloud closes on point set.Make up the point set that closes on of each point then according to integral body point cloud density, obtain the principal direction of this point set and make up Householder transformation matrix adjustment point cloud attitude.Then to closing on the match of point set march face, and then obtain two principal curvaturess of this point, select principal curvatures absolute value the greater as this Curvature Estimation based on surface equation.At last, obtain the Curvature Estimation value of whole somes clouds, as line feature point, realize the line feature extraction greater than the point of given threshold value.

The present invention is a kind of based on standardized some cloud line of a cloud attitude feature extracting method, and these method concrete steps are:

Step 1: the point set that closes on that makes up certain point: after a cloud process filtering and denoising operation; Use the KD-TREE algorithm to make up the tree construction of original point cloud; Coordinate distribution according to a cloud is sub-divided into zones of different with the original point cloud; Because the segmentation process is based on coordinate information, can directly realize the search of closest approach, significantly to improve search speed according to regional address information.Construct the point set that closes on of specified point fast.

Step 2: calculate and close on point set principal direction: use the PCA PCA, make up covariance matrix according to the point set coordinate, its minimal eigenvalue characteristic of correspondence vector is principal direction.P sets up an office _iThe point set that closes on do

(r is a some cloud number in the point set).The coordinate of promptly putting according to

, the principal direction

of calculating point set

If $P_i^r=\left[\begin{array}{ccc}{x}_1^i& y_1^i& z_1^i\\ x_2^i& y_2^i& z_2^i\\ .& .& .\\ .& .& .\\ .& .& .\\ x_r^i& y_r^i& z_r^i\end{array}\right],$ The covariance matrix that is got point set by the PCA algorithm is:

${\left(P_i^r\right)}^T\·P_i^r={\left[\begin{array}{ccc}{x}_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ x_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ .& .& .\\ .& .& .\\ .& .& .\\ x_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^r\end{array}\right]}^T\left[\begin{array}{ccc}{x}_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ x_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ .& .& .\\ .& .& .\\ .& .& .\\ x_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\end{array}\right]$ The result is 3 * 3 matrix, and its feature decomposition is got eigenvalue ₁, λ ₂, λ ₃With characteristic of correspondence vector α ₁, α ₂, α ₃, if λ ₁=min (λ ₁, λ ₂, λ ₃), then principal direction is α ₁

Step 3: point set attitude standardization: according to the principal direction of point set, make up the Householder matrix, point set is adjusted, the principal direction of point set is become (0,0,1).At first vectorial normalization is got

Can know by Householder construction method matrix, make z=(0,0,1) ^T,

Transformation matrix R=I-2bb then ^TAdjusted some cloud is

Step 4: point set surface fitting: adopt least square method, use surface equation z=ax ²+ by ²+ cxy+dx+ey+f match point set obtains the surface equation of point set.The hypothetical target surface equation is z=ax ²+ by ²+ cxy+dx+ey+f then has:

$\left\{\begin{array}{c}{z}_1={{\mathrm{ax}}_1}^2+{{\mathrm{by}}_1}^2+{\mathrm{cx}}_1{y}_1+{\mathrm{dx}}_1+e{y}_1+f\\ z_2={{\mathrm{ax}}_2}^2+{{\mathrm{by}}_2}^2+{\mathrm{cx}}_2{y}_2+{\mathrm{dx}}_2+e{y}_2+f\\ .\\ .\\ .\\ z_r={{\mathrm{ax}}_r}^2+{{\mathrm{by}}_r}^2+{\mathrm{cx}}_r{y}_r+{\mathrm{dx}}_r+{\mathrm{ey}}_r+f\end{array}\right.$

Order $A=\left[\begin{array}{cccccc}{x}_1& y_1& x_1{y}_1& x_1& y_1& 1\\ x_2& y_2& x_2{y}_2& x_2& y_2& 1\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ x_r& y_r& x_r{y}_r& x_r& y_r& 1\end{array}\right],$ Definite finally being converted into of target equation coefficient separated linear equation (A ^TA) X=A ^TL.X=[a, b, c, d, e, f] wherein ^T, L=[z ₁, z ₂, z _r] ^TDirectly separate linear equation can obtain the quadric coefficient of the match of wanting.

Step 5: some curvature is calculated: with set point along Z-direction to the quadric surface projection.Calculate two principal curvaturess in subpoint place according to surface equation, select absolute value the greater as the Curvature Estimation value.For specified point p _i(x _i, y _i, z _i), generally can be on fit Plane, need replace this moment with the projection of this point on the plane, with p _iTo curved surface projection, then subpoint is (x along Z-direction _i, y _i, ax _i ²+ by _i ²+ cx _iy _i+ dx _i+ ey _i+ f).For satisfying z=z (x, special parameter curved surface y), order

$p=\frac{\∂z}{\∂x},$ $q=\frac{\∂z}{\∂y},$

$r=\frac{{\∂}^{2}z}{{\∂x}^2},$ $s=\frac{{\∂}^{2}z}{\∂x\∂y},$ $t=\frac{{\∂}^{2}z}{{\∂y}^2}$

If have:

E=1+p ²,F=pq,G=1+q ²

$L=\frac{r}{\sqrt{1+p^2+q^2}},$ $M=\frac{s}{\sqrt{1+p^2+q^2}},$

$N=\frac{t}{\sqrt{1+p^2+q^2}}$

Then the principal curvatures value satisfies equation:

(EG-F ²)k ²-(LG-2MF+NE)k+(LN-M ²)=0

Separate k ₁, k ₂Be two principal curvaturess of this point on the curved surface.Product k ₁k ₂Be called Gaussian curvature, generally represent, average with K

Be called mean curvature, represent with H.Use mean curvature, Gaussian curvature, maximum principal curvatures, minimum principal curvatures, principal curvatures absolute value higher value as the Curvature Estimation index respectively.Traversal is obtained the Curvature Estimation value of each some cloud according to the method described above, all extracts 5% maximum curvature point as unique point, and the paired observation extraction effect can be known with principal curvatures absolute value higher value effect best.Because a bit for certain;, it becomes unique point when having a direction to change very acutely from understanding to have ready conditions; And for all bigger point of another different directions; Numerically its mean curvature Gaussian curvature of possibility can be bigger, but consider from the aspect of feature extraction, should select the index of reaction structure variation more accurately.

Step 6: unique point screening: repeat the Curvature Estimation value that above-mentioned 5 steps are obtained each point in the cloud respectively, set appropriate threshold, the Curvature Estimation value as line feature point, realizes the line feature extraction greater than the point of threshold value.

Beneficial effect: in object three-dimensional splicing and attitude measurement, use the characteristic of from the original point cloud, extracting to handle and to accomplish processing more accurately and quickly and reduce the influence that the visual angle difference is brought based on cloud data.The present invention has provided a kind of based on standardized some cloud line of a cloud attitude feature extracting method, and its advantage is:

1 usefulness KD-TREE searching method carries out the search of neighborhood, has improved search efficiency greatly.

Through the attitude of principal direction adjustment point cloud, reduce before 2 computing curvature because the influence that error of fitting is brought too greatly.And this conclusion has been passed through theoretical analysis and experimental verification.

The scheme of 3 Curvature Estimation has been selected the most significantly big principal curvatures of absolute value of border contrast effect for use, makes feature extraction efficient higher.

Description of drawings:

Fig. 1 is the process flow diagram of feature extracting method of the present invention

Embodiment: following according to the concrete embodiment of method step introduction shown in Figure 1.

The present invention is a kind of based on standardized some cloud line of a cloud attitude feature extracting method, and these method concrete steps are:

Step 1: the point set that closes on that makes up certain point: the KD-TREE algorithm is that the binary tree principle is in an application in cloud space; Coordinate distribution according to a cloud is sub-divided into different zones with the distribution of original point cloud; Owing to be based on coordinate information again in the cutting procedure; Can directly realize the search of closest approach according to the packet zone address information, be to go up the example time to become log from n with n point ₂N can practice thrift the plenty of time.Specify certain region of search can search out the point set that closes on of specified point quickly according to the KD-TREE algorithm.

Step 2: calculate and close on point set principal direction: p sets up an office _iThe point set that closes on do

(r is a some cloud number in the point set).The coordinate of promptly putting according to

, the principal direction of calculating point set

If $P_i^r=\left[\begin{array}{ccc}{x}_1^i& y_1^i& z_1^i\\ x_2^i& y_2^i& z_2^i\\ .& .& .\\ .& .& .\\ .& .& .\\ x_r^i& y_r^i& z_r^i\end{array}\right],$ The covariance matrix that is got point set by the PCA algorithm is:

${\left(P_i^r\right)}^T\·P_i^r={\left[\begin{array}{ccc}{x}_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ x_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ .& .& .\\ .& .& .\\ .& .& .\\ x_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^r\end{array}\right]}^T\left[\begin{array}{ccc}{x}_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ x_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ .& .& .\\ .& .& .\\ .& .& .\\ x_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\end{array}\right]$ The result is 3 * 3 matrix, and its feature decomposition is got eigenvalue ₁, λ ₂, λ ₃With characteristic of correspondence vector α ₁, α ₂, α ₃, if λ ₁=min (λ ₁, λ ₂, λ ₃), then principal direction is α ₁

Step 3: point set attitude standardization: owing to put the error that cloud principal direction deviation is brought, selecting the original point cloud is adjusted is that its principal direction overlaps with Z axle positive dirction in order at utmost to reduce.This conversion and become and seek the spatial alternation matrix and make vector transfer (0 to; 0,1).At first vectorial normalization is got Can know by Householder construction method matrix, make z=(0,0,1) ^T,

Transformation matrix R=I-2bb then ^TAdjusted some cloud is

Step 4: point set surface fitting: the hypothetical target surface equation is z=ax ²+ by ²+ cxy+dx+ey+f then has:

$\left\{\begin{array}{c}{z}_1={{\mathrm{ax}}_1}^2+{{\mathrm{by}}_1}^2+{\mathrm{cx}}_1{y}_1+{\mathrm{dx}}_1+e{y}_1+f\\ z_2={{\mathrm{ax}}_2}^2+{{\mathrm{by}}_2}^2+{\mathrm{cx}}_2{y}_2+{\mathrm{dx}}_2+e{y}_2+f\\ .\\ .\\ .\\ z_r={{\mathrm{ax}}_r}^2+{{\mathrm{by}}_r}^2+{\mathrm{cx}}_r{y}_r+{\mathrm{dx}}_r+{\mathrm{ey}}_r+f\end{array}\right.$

Order $A=\left[\begin{array}{cccccc}{x}_1& y_1& x_1{y}_1& x_1& y_1& 1\\ x_2& y_2& x_2{y}_2& x_2& y_2& 1\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ x_r& y_r& x_r{y}_r& x_r& y_r& 1\end{array}\right],$ Definite finally being converted into of target equation coefficient separated linear equation (A ^TA) X=A ^TL.X=[a, b, c, d, e, f] wherein ^T, L=[z ₁, z ₂, z _r] ^TDirectly separate linear equation can obtain the quadric coefficient of the match of wanting.

Step 5: some curvature is calculated: having obtained certain point so far closes on the adjusted quadric surface equation of point set principal direction, by differential geometric knowledge, can calculate the curvature of arbitrfary point on the curved surface according to surface equation.For specified point p _i(x _i, y _i, z _i), generally can be on fit Plane, need replace this moment with the projection of this point on the plane, with p _iTo curved surface projection, then subpoint is (x along Z-direction _i, y _i, ax _i ²+ by _i ²+ cx _iy _i+ dx _i+ ey _i+ f).For satisfying z=z (x, special parameter curved surface y), order

$p=\frac{\∂z}{\∂x},$ $q=\frac{\∂z}{\∂y},$

$r=\frac{{\∂}^{2}z}{\∂x^2},$ $s=\frac{{\∂}^{2}z}{\∂x\∂y},$ $t=\frac{{\∂}^{2}z}{\∂y^2}$

If have:

E=1+p ²,F=pq,G=1+q ²

$L=\frac{r}{\sqrt{1+p^2+q^2}},$ $M=\frac{s}{\sqrt{1+p^2+q^2}},$

$N=\frac{t}{\sqrt{1+p^2+q^2}}$

Then the principal curvatures value satisfies equation:

(EG-F ²)k ²-(LG-2MF+NE)k+(LN-M ²)=0

Separate k ₁, k ₂Be two principal curvaturess of this point on the curved surface.Product k ₁k ₂Be called Gaussian curvature, generally represent, average with K Be called mean curvature, represent with H.Use mean curvature, Gaussian curvature, maximum principal curvatures, minimum principal curvatures, principal curvatures absolute value higher value as the Curvature Estimation index respectively.Traversal is obtained the Curvature Estimation value of each some cloud according to the method described above, all extracts 5% maximum curvature point as unique point, and the paired observation extraction effect can be known with principal curvatures absolute value higher value effect best.Because a bit for certain;, it becomes unique point when having a direction to change very acutely from understanding to have ready conditions; And for all bigger point of another different directions; Numerically its mean curvature Gaussian curvature of possibility can be bigger, but consider from the aspect of feature extraction, should select the index of reaction structure variation more accurately.

Step 6: unique point screening: repeat the Curvature Estimation value that above-mentioned 5 steps are obtained each point in the cloud respectively, set appropriate threshold, the Curvature Estimation value as line feature point, realizes the line feature extraction greater than the point of threshold value.

Claims (1)

1. one kind based on standardized some cloud line of a cloud attitude feature extracting method, and it is characterized in that: these method concrete steps are following:

Step 1: the point set that closes on that makes up certain point: after a cloud process filtering and denoising operation; Use the KD-TREE algorithm to make up the tree construction of original point cloud; Coordinate distribution based on a cloud is sub-divided into zones of different with the original point cloud; Because the segmentation process is based on coordinate information; Directly realize the search of closest approach based on regional address information; Significantly to improve search speed, construct the point set that closes on of specified point fast;

Step 2: calculate and close on point set principal direction: use the PCA PCA, make up covariance matrix according to the point set coordinate, its minimal eigenvalue characteristic of correspondence vector is principal direction; P sets up an office _iThe point set that closes on do

R is some cloud number, i.e. a basis in the point set

The coordinate of point, the principal direction of calculating point set

If $P_i^r=\left[\begin{array}{ccc}{x}_1^i& y_1^i& z_1^i\\ x_2^i& y_2^i& z_2^i\\ .& .& .\\ .& .& .\\ .& .& .\\ x_r^i& y_r^i& z_r^i\end{array}\right],$ The covariance matrix that is got point set by the PCA algorithm is:

${\left(P_i^r\right)}^T\·P_i^r={\left[\begin{array}{ccc}{x}_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ x_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ .& .& .\\ .& .& .\\ .& .& .\\ x_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^r\end{array}\right]}^T\left[\begin{array}{ccc}{x}_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_1^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ x_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_2^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\\ .& .& .\\ .& .& .\\ .& .& .\\ x_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{x}_n^i& y_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{y}_n^i& z_r^i-\frac{1}{r}\underset{n=1}{\overset{r}{\mathrm{\Σ}}}{z}_n^i\end{array}\right]$ The result is 3 * 3 matrix, and its feature decomposition is got eigenvalue ₁, λ ₂, λ ₃With characteristic of correspondence vector α ₁, α ₂, α ₃, if λ ₁=min (λ ₁, λ ₂, λ ₃), then principal direction is α ₁

Step 3: point set attitude standardization: according to the principal direction of point set, make up the Householder matrix, point set is adjusted, the principal direction of point set is become (0,0,1); At first vectorial normalization is got

By Householder construction method matrix, make z=(0,0,1) ^T,

Transformation matrix R=I-2bb then ^T, adjusted some cloud does

Step 4: point set surface fitting: adopt least square method, use surface equation z=ax ²+ by ²+ cxy+dx+ey+f match point set obtains the surface equation of point set; The hypothetical target surface equation is z=ax ²+ by ²+ cxy+dx+ey+f then has:

$\left\{\begin{array}{c}{z}_1={{\mathrm{ax}}_1}^2+{{\mathrm{by}}_1}^2+{\mathrm{cx}}_1{y}_1+{\mathrm{dx}}_1+e{y}_1+f\\ z_2={{\mathrm{ax}}_2}^2+{{\mathrm{by}}_2}^2+{\mathrm{cx}}_2{y}_2+{\mathrm{dx}}_2+e{y}_2+f\\ .\\ .\\ .\\ z_r={{\mathrm{ax}}_r}^2+{{\mathrm{by}}_r}^2+{\mathrm{cx}}_r{y}_r+{\mathrm{dx}}_r+{\mathrm{ey}}_r+f\end{array}\right.$

Order $A=\left[\begin{array}{cccccc}{x}_1& y_1& x_1{y}_1& x_1& y_1& 1\\ x_2& y_2& x_2{y}_2& x_2& y_2& 1\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ x_r& y_r& x_r{y}_r& x_r& y_r& 1\end{array}\right],$ Definite finally being converted into of target equation coefficient separated linear equation (A ^TA) X=A ^TL; X=[a, b, c, d, e, f] wherein ^T, L=[z ₁, z ₂, z _r] ^T, directly separate linear equation obtain the quadric coefficient of the match of wanting;

Step 5: some curvature is calculated: with set point along Z-direction to the quadric surface projection, calculate two principal curvaturess in subpoint place according to surface equation, select absolute value the greater as the Curvature Estimation value; For specified point p _i(x _i, y _i, z _i), generally can be on fit Plane, need replace this moment with the projection of this point on the plane, with p _iTo curved surface projection, then subpoint is (x along Z-direction _i, y _i, ax _i ²+ by _i ²+ cx _iy _i+ dx _i+ ey _i+ f); For satisfying z=z (x, special parameter curved surface y), order

$p=\frac{\∂z}{\∂x},$ $q=\frac{\∂z}{\∂y},$

$r=\frac{{\∂}^{2}z}{{\∂x}^2},$ $s=\frac{{\∂}^{2}z}{\∂x\∂y},$ $t=\frac{{\∂}^{2}z}{\∂y^2}$

If have:

E=1+p ²,F=pq,G=1+q ²

$L=\frac{r}{\sqrt{1+p^2+q^2}},$ $M=\frac{s}{\sqrt{1+p^2+q^2}},$

$N=\frac{t}{\sqrt{1+p^2+q^2}}$

Then the principal curvatures value satisfies equation:

(EG-F ²)k ²-(LG-2MF+NE)k+(LN-M ²)=0

Separate k ₁, k ₂Be two principal curvaturess of this point on the curved surface; Product k ₁k ₂Be called Gaussian curvature, generally represent, average with K

Be called mean curvature, represent with H; Use mean curvature, Gaussian curvature, maximum principal curvatures, minimum principal curvatures, principal curvatures absolute value higher value as the Curvature Estimation index respectively; Traversal is obtained the Curvature Estimation value of each some cloud according to the method described above, all extracts 5% maximum curvature point as unique point, and the paired observation extraction effect then uses principal curvatures absolute value higher value effect best; Because a bit for certain;, it becomes unique point when having a direction to change very acutely from understanding to have ready conditions; And for all bigger point of another different directions; Numerically its mean curvature Gaussian curvature of possibility can be bigger, but consider from the aspect of feature extraction, should select the index of reaction structure variation more accurately;

Step 6: unique point screening: repeat the Curvature Estimation value that above-mentioned 5 steps are obtained each point in the cloud respectively, set appropriate threshold, the Curvature Estimation value as line feature point, realizes the line feature extraction greater than the point of threshold value.

Images