CN102799763B - A kind of based on a cloud attitude standardized some cloud line feature extraction method - Google Patents

Abstract

A kind of based on a cloud attitude standardized some cloud line feature extraction method, it has six large steps, is applicable to the extraction of line features in unordered three-dimensional point cloud, facilitates the measurement of realize target relative attitude, belong to three-dimensional measurement and technical field of machine vision.First the method builds the KD-TREE structure of a cloud, to improve the search speed that a cloud closes on point set.What then build each point according to entirety point cloud density closes on point set, obtains the principal direction of this point set and builds a Householder transformation matrix adjustment point cloud attitude.Then to closing on the face matching of point set march, and then obtaining two principal curvaturess of this point based on surface equation, selecting principal curvatures absolute value the greater as this Curvature Estimation.Finally, obtain the Curvature Estimation value of all some clouds, be greater than the point of given threshold value as line feature point, realize line feature extraction.

Description

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

Technical field

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

Technical background

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

There is several representational method about three-dimensional point cloud line feature extraction, mainly contain following a few class:

Three dimensions point cloud is projected to specific plane, planar utilizes the method finding limit or peak value to carry out feature point extraction.But the result isoplanar of plane projection is selected and original point cloud structure is closely related, the result easy receptor site cloud structure attitude of extraction and choose the multiple impact of plane and stable not.

Utilize some cloud vector with point of proximity cloud vector angle as unique point discriminant criterion, for the unique point cloud that structure changes greatly, this angle can be comparatively large, and can be smaller for this angle plane or smooth surface.With a cloud vector angle as feature selecting according to having stability to spatial alternation, but to have a calculated amount of cloud computing vector comparatively large, and noise spot is to the obvious effect of result of calculation.

In addition applying wider is feature extracting method based on curvature, describes a planarization on cloud surface by curvature.The general line feature extraction algorithm based on curvature can according to a cloud constructing curve, based on surface equation calculation level cloud Curvature Estimation value.But general algorithm still can be improved space on point set search, the adjustment of some cloud principal direction, Curvature Estimation way selection.

Summary of the invention

Technical matters: the invention provides a kind of based on a cloud attitude standardized some cloud line feature extraction method, it converges the functions such as principal direction aligning selection principal curvatures by adding point, improve a lot in characteristic mass and algorithm stability based on Curvature Methods extraction effect than existing, can be that follow-up relative attitude is measured, field stitching provides good unique point cloud.

Technical scheme: the targeted attitude based on cloud data is measured or field stitching, be widely used in practice, direct scanning institute invocation point cloud carries out analyzing and processing, point cloud quantity is more on the one hand, the speed of impact process and complexity, on the other hand due to the noise spot impact that a large amount of visual fields difference causes, final measuring accuracy is made to be difficult to ensure.Therefore propose the thought of cloud data being carried out to line feature extraction, use measurement or the identification of the more convenient target of less cloud data.

According to the corresponding relation of a cloud surface data mechanism characteristics and curvature, the present invention proposes a kind of based on a cloud attitude standardized some cloud line feature extraction method, be applicable to the extraction of line features in unordered three-dimensional point cloud, first the method builds the KD-TRE E structure of a cloud, to improve the search speed that a cloud closes on point set.What then build each point according to entirety point cloud density closes on point set, obtains the principal direction of this point set and builds a Householder transformation matrix adjustment point cloud attitude.Then to closing on the face matching of point set march, and then obtaining two principal curvaturess of this point based on surface equation, selecting principal curvatures absolute value the greater as this Curvature Estimation.Finally, obtain the Curvature Estimation value of all some clouds, be greater than the point of given threshold value as line feature point, realize line feature extraction.

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

Step one: what build certain point closes on point set: after denoising operation is involved after filtration to a cloud, KD-TREE algorithm is used to build the tree construction of original point cloud, original point cloud is sub-divided into zones of different by the coordinate distribution according to a cloud, because segmentation process is based on coordinate information, directly can realize the search of closest approach according to zone-address information, significantly to improve search speed.What rapid build went out specified point closes on point set.

Step 2: calculate and close on point set principal direction: use PCA principal component analysis (PCA), build covariance matrix according to point set coordinate, its minimal eigenvalue characteristic of correspondence vector is principal direction.Set up an office p _ithe point set that closes on be (r is some cloud number in point set).I.e. basis the coordinate of point, calculates the principal direction of 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 being obtained point set by 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]$ Result is the matrix of 3 × 3, obtains eigenvalue λ to its feature decomposition ₁, λ ₂, λ ₃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, builds Householder matrix, adjusts, the principal direction of point set is become (0,0,1) to point set.First vectorial normalization is obtained from Householder construction method matrix, make z=(0,0,1) ^t, then transformation matrix R=I-2bb ^t.Point cloud after adjustment is

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

$\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],$ The determination of target equation coefficient is finally converted into separates linear equation (A ^ta) X=A ^tl.Wherein X=[a, b, c, d, e, f] ^t, L=[z ₁, z ₂, z _r] ^t.Directly solution linear equation can obtain the quadric coefficient of wanted matching.

Step 5: some curvature estimation: set point is projected along Z-direction to quadric surface.Calculate two, subpoint place principal curvatures according to surface equation, select absolute value the greater as Curvature Estimation value.For specified point p _i(x _i, y _i, z _i), generally can not, on fit Plane, now need to replace, by p with the projection of this point in plane _ialong Z-direction to curved surface projection, then subpoint is (x _i, y _i, ax _i ²+ by _i ²+ cx _iy _i+ dx _i+ ey _i+ f).For the special parameter curved surface meeting z=z (x, 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 principal curve value meets equation:

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

Separate k ₁, k ₂two principal curvaturess of this point on curved surface.Product k ₁k ₂be called Gaussian curvature, generally represent with K, average 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 Curvature Estimation index respectively.Traversal obtains the Curvature Estimation value of each some cloud according to the method described above, and all extract the maximum curvature point of 5% as unique point, paired observation extraction effect, known principal curvatures absolute value higher value effect is best.Because for certain a bit, have ready conditions from understanding when it has a direction change very violent and become unique point, and for all larger point of another different directions, numerically its mean curvature Gaussian curvature of possibility can be larger, but from the viewpoint of feature extraction, should select can the index of reaction structure change more accurately.

Step 6: unique point is screened: repeat the Curvature Estimation value that above-mentioned 5 steps obtain each point in a cloud respectively, setting appropriate threshold, Curvature Estimation value is greater than the point of threshold value as line feature point, realizes line feature extraction.

Beneficial effect: based in the object three-dimensional splicing of cloud data and attitude measurement, uses the feature extracted from original point cloud to carry out processing completing more accurately and quickly and process and reduce the different impact brought in visual angle.The present invention gives a kind of based on a cloud attitude standardized some cloud line feature extraction method, its advantage is:

1 carries out the search of neighborhood with KD-TREE searching method, substantially increases search efficiency.

2 calculate the attitude by principal direction adjustment point cloud before curvature, reduce the impact brought too greatly due to error of fitting.And this conclusion have passed through theoretical analysis and experimental verification.

The scheme of 3 Curvature Estimation has selected the larger principal curvatures of the most obvious absolute value of border contrast effect, makes feature extraction efficiency higher.

Accompanying drawing illustrates:

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

Embodiment: below according to the embodiment that the method step introduction shown in Fig. 1 is concrete.

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

Step one: what build certain point closes on point set: KD-TREE algorithm is that binary tree principle is in an application in cloud space, the distribution of original point cloud is sub-divided into different regions by the coordinate distribution according to a cloud, due in cutting procedure again based on coordinate information, directly can realize the search of closest approach according to packet zone address information, so that n some time to become log from n ₂n, can save the plenty of time.That specifies certain region of search can search out specified point quickly according to KD-TREE algorithm closes on point set.

Step 2: calculate and close on point set principal direction: set up an office p _ithe point set that closes on be (r is some cloud number in point set).I.e. basis the coordinate of point, calculates the principal direction of 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 being obtained point set by 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]$ Result is the matrix of 3 × 3, obtains eigenvalue λ to its feature decomposition ₁, λ ₂, λ ₃with characteristic of correspondence vector α ₁, α ₂, α ₃if, λ ₁=min (λ ₁, λ ₂, λ ₃), then principal direction is α ₁.

Step 3: point set attitude standardization: in order at utmost reduce due to an error that cloud principal direction deviation is brought, selecting to carry out adjustment to original point cloud is that its principal direction overlaps with Z axis positive dirction.This conversion and become find space conversion matrices make vector transfer to (0,0,1).First vectorial normalization is obtained from Householder construction method matrix, make z=(0,0,1) ^t, then transformation matrix R=I-2bb ^t.Point cloud after adjustment is

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

$\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],$ The determination of target equation coefficient is finally converted into separates linear equation (A ^ta) X=A ^tl.Wherein X=[a, b, c, d, e, f] ^t, L=[z ₁, z ₂, z _r] ^t.Directly solution linear equation can obtain the quadric coefficient of wanted matching.

Step 5: some curvature estimation: obtained certain so far and put the Quadratic Surface Equation after closing on the adjustment of point set principal direction, by differential geometric knowledge, the curvature of arbitrfary point on curved surface can have been calculated according to surface equation.For specified point p _i(x _i, y _i, z _i), generally can not, on fit Plane, now need to replace, by p with the projection of this point in plane _ialong Z-direction to curved surface projection, then subpoint is (x _i, y _i, ax _i ²+ by _i ²+ cx _iy _i+ dx _i+ ey _i+ f).For the special parameter curved surface meeting z=z (x, 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 principal curve value meets equation:

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

Separate k ₁, k ₂two principal curvaturess of this point on curved surface.Product k ₁k ₂be called Gaussian curvature, generally represent with K, average 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 Curvature Estimation index respectively.Traversal obtains the Curvature Estimation value of each some cloud according to the method described above, and all extract the maximum curvature point of 5% as unique point, paired observation extraction effect, known principal curvatures absolute value higher value effect is best.Because for certain a bit, have ready conditions from understanding when it has a direction change very violent and become unique point, and for all larger point of another different directions, numerically its mean curvature Gaussian curvature of possibility can be larger, but from the viewpoint of feature extraction, should select can the index of reaction structure change more accurately.

Step 6: unique point is screened: repeat the Curvature Estimation value that above-mentioned 5 steps obtain each point in a cloud respectively, setting appropriate threshold, Curvature Estimation value is greater than the point of threshold value as line feature point, realizes line feature extraction.

Claims (1)

1., based on a cloud attitude standardized some cloud line feature extraction method, it is characterized in that: the method concrete steps are as follows:

Step one: what build certain point closes on point set: after denoising operation is involved after filtration to a cloud, KD-TREE algorithm is used to build the tree construction of original point cloud, original point cloud is sub-divided into zones of different by the coordinate distribution according to a cloud, because segmentation process is based on coordinate information, the search of closest approach is directly realized according to zone-address information, significantly to improve search speed, what rapid build went out specified point closes on point set;

Step 2: calculate and close on point set principal direction: use PCA principal component analysis (PCA), build covariance matrix according to point set coordinate, its minimal eigenvalue characteristic of correspondence vector is principal direction; Set up an office p _ithe point set that closes on be r is some cloud number in point set, i.e. basis the coordinate of point, calculates the principal direction of 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 being obtained point set by 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^i\end{array}\right]\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]$ Result is the matrix of 3 × 3, obtains eigenvalue λ to its feature decomposition ₁, λ ₂, λ ₃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, builds Householder matrix, adjusts, the principal direction of point set is become (0,0,1) to point set; First vectorial normalization is obtained by Householder construction method matrix, make z=(0,0,1) ^t, then transformation matrix R=I-2bb ^t, the some cloud after adjustment is $Q_i^r={(R\·{\left(P_i^r\right)}^T)}^T;$

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

$\left\{\begin{array}{c}{z}_1={{\mathrm{ax}}_1}^2+{{\mathrm{by}}_1}^2+{\mathrm{cx}}_1{y}_1+{\mathrm{dx}}_1+{\mathrm{ey}}_1+f\\ z_2={{\mathrm{ax}}_2}^2+{{\mathrm{by}}_2}^2+{\mathrm{cx}}_2{y}_2+{\mathrm{dx}}_2+{\mathrm{ey}}_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],$ The determination of target equation coefficient is finally converted into separates linear equation (A ^ta) X=A ^tl; Wherein X=[a, b, c, d, e, f] ^t, L=[z ₁, z ₂, z _r] ^t, Directly solution linear equation obtains the quadric coefficient of wanted matching;

Step 5: some curvature estimation: projected along Z-direction to quadric surface by set point, calculates two, subpoint place principal curvatures according to surface equation, selects absolute value the greater as Curvature Estimation value; For specified point p _i(x _i, y _i, z _i), can not, on fit Plane, now need to replace, by p with the projection of this point in plane _ialong Z-direction to curved surface projection, then subpoint is (x _i, y _i, ax _i ²+ by _i ²+ cx _iy _i+ dx _i+ ey _i+ f); For the special parameter curved surface meeting z=z (x, 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 principal curve value meets equation:

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

Separate k ₁, k ₂two principal curvaturess of this point on curved surface; Product k ₁k ₂be called Gaussian curvature, represent with K, average 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 Curvature Estimation index respectively; Traversal obtains the Curvature Estimation value of each some cloud according to the method described above, and all extract the maximum curvature point of 5% as unique point, paired observation extraction effect, then use principal curvatures absolute value higher value effect best; Because for certain a bit, have ready conditions from understanding when it has a direction change very violent and become unique point, and for all larger point of another different directions, numerically its mean curvature Gaussian curvature can be large, but from the viewpoint of feature extraction, selecting can the index of reaction structure change more accurately;

Step 6: unique point is screened: repeat the Curvature Estimation value that above-mentioned 5 steps obtain each point in a cloud respectively, setting appropriate threshold, Curvature Estimation value is greater than the point of threshold value as line feature point, realizes line feature extraction.