CN109345571B - Point cloud registration method based on extended Gaussian image - Google Patents

Point cloud registration method based on extended Gaussian image Download PDF

Info

Publication number
CN109345571B
CN109345571B CN201811190655.4A CN201811190655A CN109345571B CN 109345571 B CN109345571 B CN 109345571B CN 201811190655 A CN201811190655 A CN 201811190655A CN 109345571 B CN109345571 B CN 109345571B
Authority
CN
China
Prior art keywords
point cloud
gaussian
grid
normals
rotation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811190655.4A
Other languages
Chinese (zh)
Other versions
CN109345571A (en
Inventor
杜志江
高永卓
董为
徐威
李明洋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN201811190655.4A priority Critical patent/CN109345571B/en
Publication of CN109345571A publication Critical patent/CN109345571A/en
Application granted granted Critical
Publication of CN109345571B publication Critical patent/CN109345571B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds

Abstract

The invention discloses a point cloud registration method based on an extended Gaussian image, and relates to a point cloud registration method of an image. The invention aims to solve the problems that the cost of a rough registration algorithm SAC-IA algorithm commonly used in the field of point cloud registration is too high, a Gaussian sphere is equally divided into a polyhedron or longitude and latitude based on an EGI algorithm, the time is long, and when a rotation matrix is calculated based on the EGI algorithm, the optimal solution is directly searched in a feasible region, and the dual requirements of efficiency and precision are difficult to balance. The process is as follows: let the coordinate set of the normal line under each body coordinate system be NTAnd NS(ii) a Gridding and dividing the Gaussian expansion image; sorting the preliminarily matched grids in a descending order according to the number of normals in the grids; obtaining a corrected grid; calculating a rotation matrix; outputs iter, Corr and M. The invention is used in the cross field of computer graphics and machining technology.

Description

Point cloud registration method based on extended Gaussian image
Technical Field
The invention belongs to the cross field of computer graphics and machining technology, and particularly relates to a point cloud registration method of an image.
Background
In the process of processing a workpiece, casting is often used for manufacturing a blank, and then a material reduction processing mode such as turning and grinding is used for processing a part into a final shape. In the process, the condition of the casting needs to be detected, the distribution of the allowance is determined, and then the distribution of the allowance is further optimized, so that the next step of cutting and machining is facilitated. The laser scanner is more suitable for detecting the casting compared with manual work due to the characteristics of convenience in surface point cloud data rapidness, low time cost and higher precision. After the scanner acquires the point cloud data of the component, the point cloud is registered and compared with a model designed in CAD software by using a point cloud registration algorithm, and the method can efficiently and reliably analyze the machining allowance after the machining.
The field of point cloud registration makes many relevant researches on point cloud registration of various different scenes and different appearances. A rough registration algorithm SAC-IA commonly used in the field of point cloud registration minimizes the distance between a target point cloud and a model point cloud in an iterative mode, and registration is completed. The algorithm overhead is excessive. There are documents that a class of algorithms that use Extended Gaussian Images (EGI) to estimate a rotation matrix and use an ICP algorithm to complete registration equally divides a Gaussian sphere into polyhedrons or longitudes, the areas of each surface are the same, however, it is difficult to obtain a specific coordinate of each edge angle of a polyhedron, and when the number of normals is huge, it is necessary to judge and calculate whether each normal points to the surface one by one, which takes a long time. When the rotation matrix is calculated, the optimal solution is directly searched in a feasible domain, and the dual requirements of efficiency and precision are difficult to balance.
EGI is an extended Gaussian image;
disclosure of Invention
The invention aims to solve the problems that the cost of a rough registration algorithm SAC-IA algorithm commonly used in the field of point cloud registration is too high, a Gaussian sphere is equally divided into a polyhedron or longitude and latitude based on an EGI algorithm, the time is long, and when a rotation matrix is calculated based on the EGI algorithm, the optimal solution is directly searched in a feasible region, and the dual requirements of efficiency and precision are difficult to balance, and provides a point cloud registration method based on an extended Gaussian image.
A point cloud registration method based on an extended Gaussian image comprises the following specific processes:
step one, recording a coordinate set of normals in model point cloud and blank point cloud under respective body coordinate system as NTAnd NS(get normal position);
step two, respectively carrying out gridding division on Gaussian expansion images (EGI) generated by the model point cloud and the blank point cloud;
step three, based on step two to grid Pi T、
Figure BDA0001827422200000013
And a grid Pi s、
Figure BDA0001827422200000012
Carrying out preliminary matching and matching the preliminarily matched grids Pi T、
Figure BDA0001827422200000021
And a grid Pi s、
Figure BDA0001827422200000022
Sorting in descending order according to the number of normals in the grid, and recording as P1 T″、
Figure BDA00018274222000000216
And P1 s″、
Figure BDA00018274222000000217
Wherein, Pi SIs composed of
Figure BDA0001827422200000023
The grid with the largest number of medium normals; pi TIs composed of
Figure BDA00018274222000000221
The grid with the largest number of medium normals;
Figure BDA0001827422200000025
is composed of
Figure BDA0001827422200000026
The grid with the largest number of medium normals;
Figure BDA0001827422200000027
is composed of
Figure BDA0001827422200000028
The grid with the largest number of medium normals;
Figure BDA0001827422200000029
the ith base mesh of the gaussian expanded image generated for the blank point cloud,
Figure BDA00018274222000000210
the ith base mesh of the gaussian expanded image generated for the model point cloud,
Figure BDA00018274222000000211
a jth base mesh of a gaussian expanded image generated for the blank point cloud,
Figure BDA00018274222000000212
generating a jth basic grid of a Gaussian expansion image for the model point cloud;
i＝1,2,...,12；j＝i+1,i+2,...,12；i≠j；
step four, P is paired1 T″、
Figure BDA00018274222000000218
And P1 s″、
Figure BDA00018274222000000219
Performing correction to obtain corrected P1 T′、
Figure BDA00018274222000000213
And P1 s′、
Figure BDA00018274222000000220
Step five, according to the corrected P1 T′、
Figure BDA00018274222000000214
And P1 s′、
Figure BDA00018274222000000215
Computing gaussian extended image (EGI) and blank point cloud enabling model point cloud generationA rotation matrix of a gaussian extended image (EGI) registration;
sixthly, outputting iter, Corr and M based on the step rotation matrix;
m is a final output rotation matrix;
corr is the correlation between a Gaussian expansion image generated by the blank point cloud and a Gaussian expansion image generated by the model point cloud;
iter is the number of iterations.
The invention has the beneficial effects that:
the algorithm for point cloud registration in allowance analysis of machining occasions is provided, and the defects of the conventional EGI-based point cloud registration algorithm are improved to a certain extent in the parts of a gridding algorithm, corresponding point identification, rotation matrix calculation, an iteration strategy and the like. The calculation complexity is obviously lower than that of the traditional algorithm, and the algorithm overhead is reduced. Under the background of using laser scanning to realize precision machining allowance analysis, a dividing mode of a spherical surface is changed by adopting an HEALPIX algorithm, a KDtree is constructed to determine a normal included by each surface, dividing efficiency is improved, and time is shortened. And the patches are further subdivided to improve the precision, and a dynamic average value calculation mode is adopted to prevent excessively fine division fuzzy peak values. And iteration and variation are used, so that the precision is improved on one hand, and the dual requirements of local optimization, efficiency balance and precision are avoided on the other hand.
Equally dividing a Gaussian sphere into polyhedrons or changing the longitude and the latitude based on an EGI algorithm into a HEALPIX algorithm, and changing the division mode of the sphere; when the rotation matrix is calculated by the EGI algorithm, the optimal solution is directly searched in the feasible region, and a definite calculation formula is provided instead to calculate the rotation matrix and an iteration strategy. The algorithm overhead of the invention is reduced by about 80-85% compared with the SAC-IA algorithm. The accuracy of the present invention is measured in RMS (root mean square), which is about an order of magnitude higher.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a schematic diagram of the present invention dividing a Gaussian extended image into 12 grids by a 2HEALPIX method;
FIG. 3 is a schematic diagram of the present invention for segmenting a Gaussian extended image into 48 grids by a 2HEALPIX method;
FIG. 4 is a schematic diagram of the present invention dividing a Gaussian extended image into 192 grids by a 2HEALPIX method;
FIG. 5 is a schematic diagram of the present invention dividing a Gaussian extended image into 768 grids by a 2HEALPIX method;
FIG. 6 is a schematic view of a wind blade.
Detailed Description
The first embodiment is as follows: the embodiment is described with reference to fig. 1, and a specific process of the point cloud registration method based on the extended gaussian image in the embodiment is as follows:
step one, recording a coordinate set of normals in model point cloud and blank point cloud under respective body coordinate system as NTAnd NS(get normal position);
the model point cloud is a point cloud of an ideal mechanical part model;
the blank point cloud is the point cloud of a blank (the blank is processed in a factory and is obtained by scanning through a scanner);
step two, respectively carrying out gridding division on Gaussian expansion images (EGI) generated by the model point cloud and the blank point cloud;
step three, based on step two to grid Pi T、
Figure BDA0001827422200000031
And a grid Pi s、
Figure BDA0001827422200000032
Carrying out preliminary matching and matching the preliminarily matched grids Pi T、
Figure BDA0001827422200000033
And a grid Pi s、
Figure BDA0001827422200000034
Sorting in descending order according to the number of normals in the grid, and recording as P1 T″、
Figure BDA0001827422200000035
And P1 s″、
Figure BDA0001827422200000036
Wherein, Pi SIs composed of
Figure BDA0001827422200000037
The grid with the largest number of medium normals; pi TIs composed of
Figure BDA0001827422200000038
The grid with the largest number of medium normals;
Figure BDA0001827422200000039
is composed of
Figure BDA00018274222000000310
The grid with the largest number of medium normals;
Figure BDA00018274222000000311
is composed of
Figure BDA00018274222000000312
The grid with the largest number of medium normals;
Figure BDA00018274222000000313
the ith base mesh of the gaussian expanded image generated for the blank point cloud,
Figure BDA00018274222000000314
the ith base mesh of the gaussian expanded image generated for the model point cloud,
Figure BDA00018274222000000315
a jth base mesh of a gaussian expanded image generated for the blank point cloud,
Figure BDA00018274222000000316
generating a jth basic grid of a Gaussian expansion image for the model point cloud;
i＝1,2,...,12；j＝i+1,i+2,...,12；i≠j；
step four, P is paired1 T″、
Figure BDA00018274222000000317
And P1 s″、
Figure BDA00018274222000000318
Performing correction to obtain corrected P1 T′、
Figure BDA00018274222000000319
And P1 s′、
Figure BDA00018274222000000320
Step five, according to the corrected P1 T′、
Figure BDA00018274222000000321
P1 s′、
Figure BDA00018274222000000322
Calculating a rotation matrix enabling registration of a Gaussian expansion image (EGI) generated by the model point cloud and a Gaussian expansion image (EGI) generated by the blank point cloud;
sixthly, outputting iter, Corr and M based on the step rotation matrix;
m is a final output rotation matrix;
corr is the correlation between a Gaussian expansion image generated by the blank point cloud and a Gaussian expansion image generated by the model point cloud;
iter is the number of iterations.
The second embodiment is as follows: the first difference between the present embodiment and the specific embodiment is: and in the second step, the Gaussian expansion image (EGI) generated by the model point cloud and the blank point cloud is respectively subjected to gridding division, and the specific process is as follows:
HEALPIX (structural equivalent Area isotope labeling) was originally studied in the context of cosmic microwaves to support high resolution discretization of spheres. HEALPIX can be divided into several levels according to resolution, and the lowest level is 12 grids. HEALPIX is divided into spherical surfaces as shown in FIGS. 2, 3, 4, and 5.
HEALPIX has the following characteristics:
the grid blocks are constructed in a layered mode, and data are convenient to access.
The method for generating the grid block is simple and efficient in operation.
The spherical surface is divided into curved surface quadrangles with equal size, and the area of each grid is the same.
Dividing a Gaussian expansion image generated by the model point cloud and the blank point cloud into 12 basic grids, wherein HEALPIX (HEALPIX is an acronym for Hierarchical Equal Area ideal stereoscopic image of a space. As acquired in the name, this pixel processing a sub-division of a spatial surface in a white Area of a spatial surface) is divided into 4 grids for each grid of the previous level layer by layer in a recursive manner until m grids are reached, and m is 12 × 4nN is an integer, m is 12 × 4n＝3072；
Calculating the number of normals in each grid by adopting a kd-tree algorithm (the number of normals is determined by the number of points in the model, and the normal corresponding to each point is calculated once for each point), and storing the normals in a vector form (each element represents the number of the normals of the grid); simultaneously recording the coordinates of the normal lines under respective body coordinate systems;
indexes are established for each grid in a mode of 'ring' and 'nesting', so that data can be accessed and processed quickly in different modes.
According to a nested index mode of HEALPIX, a Gaussian expansion image H generated by blank point cloudSGaussian expansion image H generated by model point cloudTDivided into 12 parts (each part has 4) according to 12 basic gridsn(256) Number of normals of each grid), respectively
Figure BDA0001827422200000041
And
Figure BDA0001827422200000042
i＝1,2,...,12；
wherein
Figure BDA0001827422200000043
The 1 st base mesh of the gaussian expanded image generated for the blank point cloud,
Figure BDA0001827422200000044
the ith base mesh of the gaussian expanded image generated for the blank point cloud,
Figure BDA0001827422200000051
generating a 12 th basic grid of a Gaussian expansion image for the blank point cloud;
Figure BDA0001827422200000052
the 1 st base mesh of the gaussian expanded image generated for the model point cloud,
Figure BDA0001827422200000053
the ith base mesh of the gaussian expanded image generated for the model point cloud,
Figure BDA0001827422200000054
a 12 th base mesh of a gaussian extended image generated for the model point cloud.
Other steps and parameters are the same as those in the first embodiment.
The third concrete implementation mode: the present embodiment differs from the first or second embodiment in that: two pairs of grids P based on the steps in the third stepi T、
Figure BDA0001827422200000055
And a grid Pi s、
Figure BDA0001827422200000056
Performing preliminary matching, and performing the preliminary matchingGrid P ofi T、
Figure BDA0001827422200000057
And a grid Pi s、
Figure BDA0001827422200000058
Sorting in descending order according to the number of normals in the grid, and recording as P1 T″、
Figure BDA0001827422200000059
And P1 s″、
Figure BDA00018274222000000510
The specific process is as follows:
corresponding point identification
Finding the corresponding mesh can be a challenging task, especially when the normals of the objects are distributed to some extent evenly. Due to the distribution characteristics of normal vectors, a high density mesh may have a similar density mesh in the opposite direction, thereby causing errors in the rotation matrix calculation. In addition, a good algorithm should be insensitive to rotation and position changes of the grid. Here we use a simple and efficient method to find the grid in the target EGI that is related to the grid in the source EGI based on the number of normals and the length relationship between the selected bins. The corresponding point identification is carried out according to the following steps:
step three one, traverse
Figure BDA00018274222000000511
And
Figure BDA00018274222000000512
respectively find
Figure BDA00018274222000000513
And
Figure BDA00018274222000000514
grid formation with the greatest number of normals
Figure BDA00018274222000000515
And
Figure BDA00018274222000000516
i＝1,2,...,12；
wherein P isi SIs composed of
Figure BDA00018274222000000517
The grid with the largest number of medium normals; pi TIs composed of
Figure BDA00018274222000000538
The grid with the largest number of medium normals;
step three, two, pair
Figure BDA00018274222000000519
And
Figure BDA00018274222000000520
sorting in descending order according to the number of the normal lines respectively;
finding a point P starting from the center of the Gaussian expansion imagei s、
Figure BDA00018274222000000521
And Pi T、
Figure BDA00018274222000000522
Is marked as Vi S、
Figure BDA00018274222000000523
And Vi T、
Figure BDA00018274222000000524
Wherein
Figure BDA00018274222000000525
Is composed of
Figure BDA00018274222000000526
The grid with the largest number of medium normals;
Figure BDA00018274222000000540
is composed of
Figure BDA00018274222000000528
The grid with the largest number of medium normals;
by each Pi sJudgment of others
Figure BDA00018274222000000539
Whether or not to satisfy
Figure BDA00018274222000000530
If not, discarding
Figure BDA00018274222000000531
If so, retaining
Figure BDA00018274222000000532
j＝i+1,i+2,...,12；i≠j；
By each Pi TJudgment of others
Figure BDA00018274222000000541
Whether or not to satisfy
Figure BDA00018274222000000534
If not, discarding
Figure BDA00018274222000000535
If so, retaining
Figure BDA00018274222000000536
j＝i+1,i+2,...,12；i≠j；
Step three, calculating the distance between any two grids in the remaining grids
Figure BDA00018274222000000537
And
Figure BDA0001827422200000061
Figure BDA0001827422200000062
a jth base mesh of a gaussian expanded image generated for the blank point cloud,
Figure BDA0001827422200000063
generating a jth basic grid of a Gaussian expansion image for the model point cloud;
Figure BDA0001827422200000064
is PSThe distance between any two of the remaining meshes;
Figure BDA0001827422200000065
is PTThe distance between any two of the remaining meshes;
DSis PSA set of distances between any two of the remaining meshes;
DTis PTA set of distances between any two of the remaining meshes;
judgment of DSAnd DTWhether there is a similar value in
Figure BDA0001827422200000066
And
Figure BDA0001827422200000067
if it is
Figure BDA0001827422200000068
And
Figure BDA0001827422200000069
if the difference is within 5%, the grid P isi T、
Figure BDA00018274222000000610
And mesh Pi s、
Figure BDA00018274222000000620
Performing primary matching;
tests have shown that this method effectively excludes closely spaced dense points, while two dense points with higher normal density are identified in two EGIs by their geometric distance, which is rotation and noise insensitive.
Mesh P after preliminary matchingi T、
Figure BDA00018274222000000612
And a grid Pi s、
Figure BDA00018274222000000613
Sorting in descending order according to the number of normals in the grid, and recording as P1 T″、
Figure BDA00018274222000000614
And P1 s″、
Figure BDA00018274222000000615
Other steps and parameters are the same as those in the first or second embodiment.
The fourth concrete implementation mode: the difference between this embodiment mode and one of the first to third embodiment modes is: the fourth step is to P1 T″、
Figure BDA00018274222000000616
And P1 s″、
Figure BDA00018274222000000617
Performing correction to obtain corrected P1 T′、
Figure BDA00018274222000000618
And P1 s′、
Figure BDA00018274222000000619
The specific process is as follows:
coordinate optimization of selected points
Taking 3072 meshes generated by the HEAPIX as an example, even if the number of segmentations is large, the error is still large if the angle of each two adjacent mesh centers to the sphere is about 4 degrees. Note also that a dense grid may form a dense normal region with neighboring grid clusters. After finding the paired dense regions in both EGIs. Instead of using the geometric center point of the highest scoring mesh in the region, which is usually not equal, we can improve the accuracy by computing the center point of the entire region on the basis of the maximum density. It indicates that we need to further improve the precision, and it is mentioned that simply improving the density of the divided grids may reduce the number of normals contained in each grid of the normal dense region, which is difficult to identify and is not beneficial to the post-registration. The dynamic averaging method is used here, i.e. the size of each "window" is fixed, and the position of the window center is varied, visually as if the window were moved in very small steps within the found dense area, the window size is fixed, ensuring that there are always enough normals to fall within the window, and the window position is varied until a value is reached that maximizes the number of normals that fall within the window. The specific algorithm steps are as follows.
Sorting P in descending order according to normal quantity1 T″、
Figure BDA0001827422200000071
And P1 s″、
Figure BDA0001827422200000072
Each of the four grids performs the following steps:
step four, finding the separation P0The most recent n grids, forming a set { P }0,...,PnFind the set { P }0,...,PnNormal bit of each grid inSetting (knowing the normal position by knowing the normal quantity in front) to form a normal set NP；
P0Is P1 T″、
Figure BDA0001827422200000073
P1 s″、
Figure BDA0001827422200000074
One of (a);
step four, finding a point P which points to the center of the Gaussian expansion image as a starting point0Constructing a plane A perpendicular to V, and the plane A is used for gathering the normal line N by Gaussian expansion of the spherical center of the imagePProjecting to a plane A vertical to V to obtain a normal set NPTwo-dimensional coordinate set N corresponding to each normal in the imagei；
Step three, generating a set PA (a plurality of grids) of uniformly distributed points in the plane A, covering the normal lines in the plane A by the set PA of points, determining the number of the normal lines in a circle with a fixed radius and taking each point in the set PA as the center by using a KD-tree method, and forming an array HA；
Step four, finding the array HAPoint (x) with the largest number of normals0,y0) Will (x)0,y0) Projecting along the vector V back to the Gaussian expansion image sphere to obtain P0New coordinate value, finish to P0Correcting; (x)0,y0)∈Ni。
Other steps and parameters are the same as those in one of the first to third embodiments.
The fifth concrete implementation mode: the difference between this embodiment and one of the first to fourth embodiments is: the fixed radius is:
Figure BDA0001827422200000075
wherein, tau is positive integer (used for controlling the size of the circle); (x)i,yi)∈Ni(ii) a range is distance; x is the number ofiSet N as normalPTwo-dimensional abscissa, y, corresponding to the ith normal lineiSet N as normalPThe ith normal line.
Other steps and parameters are the same as in one of the first to fourth embodiments.
The sixth specific implementation mode: the difference between this embodiment and one of the first to fifth embodiments is: according to the corrected P in the step five1 T′、
Figure BDA0001827422200000076
P1 s′、
Figure BDA0001827422200000077
Calculating a rotation matrix enabling registration of a Gaussian expansion image (EGI) generated by the model point cloud and a Gaussian expansion image (EGI) generated by the blank point cloud; the specific process is as follows:
the higher fractional mesh is aligned first and then the two meshes are kept aligned, along an axis pointing to the previously aligned mesh. The two EGIs rotated until the two grids with lower scores align. The two rotation vectors are represented as a rotation matrix, and the alignment of the two point clouds is accomplished by applying this matrix to each point in the target point cloud. The specific formula is as follows:
respectively pointing to P from the center of the sphere of the Gaussian expansion image1 T′、
Figure BDA0001827422200000081
And P1 s′、
Figure BDA0001827422200000082
Is marked as V1 T、
Figure BDA0001827422200000083
And V1 s、
Figure BDA0001827422200000084
Defining a rotation vector V1And V2Angle of rotation theta1And theta2；
Calculating a rotation vector V1Angle of rotation theta1The formula is as follows:
Figure BDA0001827422200000085
Figure BDA0001827422200000086
Figure BDA0001827422200000087
the solving formula is as follows:
Figure BDA0001827422200000088
Figure BDA0001827422200000089
wherein
Figure BDA00018274222000000810
Is an intermediate variable;
initial rotation matrix R1By a rotation vector V1And a rotation angle theta1The formula is obtained as follows:
Figure BDA00018274222000000811
wherein I is an identity matrix, V1 HIs a rotation vector V1Transpose of (V)1xIs a rotation vector V1Component in the x-direction, V1yIs a rotation vector V1Component in the y-direction, V1zIs a rotation vector V1A component in the z direction;
calculating a rotation vector V2Angle of rotation theta2The formula is as follows:
V2＝V1 T
Figure BDA00018274222000000812
Figure BDA00018274222000000813
Figure BDA00018274222000000814
wherein
Figure BDA00018274222000000815
Is an intermediate variable;
initial rotation matrix R2By a rotation vector V2And a rotation angle theta2The formula is obtained as follows:
Figure BDA0001827422200000091
wherein
Figure BDA0001827422200000092
Is a rotation vector V2Transpose of (V)2xIs a rotation vector V2Component in the x-direction, V2yIs a rotation vector V2Component in the y-direction, V2zIs a rotation vector V2A component in the z direction;
according to R1And R2Obtaining a rotation matrix R:
R＝R1R2。
other steps and parameters are the same as those in one of the first to fifth embodiments.
The seventh embodiment: the difference between this embodiment and one of the first to sixth embodiments is: in the sixth step, based on the step rotation matrix, iter, Corr and M are output; the specific process is as follows:
step six, defining the iteration times iter as 0 and defining the hysteresis stall as 0; definition M ═ I;
wherein I is an identity matrix; m is a final output rotation matrix;
update Hs＝R×Hs
Calculate HTAnd after update HsThe formula is as follows:
Figure BDA0001827422200000093
the correlation function is the maximum correlation of two EGI images when the correlation function is 1.
Sixthly, updating M to be R multiplied by M, and if M is larger than or equal to the threshold value 0.98, outputting iter, Corr and M;
if M is less than the threshold value of 0.98, executing step six and three;
sixthly, recording coordinate sets of normals in the model point cloud and the blank point cloud under respective body coordinate systems as NTAnd Ns′(get normal position); step six or four is executed;
Ns′＝R×Ns
sixthly, re-executing the step two to the step six to obtain new HsAnd HTCalculating new HsAnd HTThe correlation function is as follows:
Figure BDA0001827422200000101
adding 1 to Stall;
executing the step six and five;
sixthly, comparing Corr and Corr':
if Corr × 1.05 is less than Corr', perform Stall decrease by 1; execute the step six
If Corr is less than or equal to Corr ', let Corr equal to Corr', update M ═ R × M, let N bes＝Mrand×Ns(ii) a Adding 1 to Stall; execute the step six
(Corr is greater than Corr' not considered);
the M israndA rotation matrix corresponding to the rotation angle theta;
theta is a randomly generated rotation angle around XYZ three axes of the point cloud body coordinate system;
sixthly, judging whether Stall is greater than 3 and Corr' is less than or equal to 0.98, if so, judging that M is equal to Mrand× M, order Ns＝Mrand×Ns(ii) a Sixthly, executing the step;
if not, executing step six and seven;
sixthly, if the stall is 0, iter + 1;
if Corr' is less than the threshold value of 0.98 and iter is less than or equal to 15, repeating the steps sixty-four to sixty-seven, and outputting iter, Corr and M;
if not, iter, Corr and M are output.
Other steps and parameters are the same as those in one of the first to sixth embodiments.
The specific implementation mode is eight: the present embodiment differs from one of the first to seventh embodiments in that: a rotation matrix M corresponding to the rotation angle thetarandThe solving process is as follows:
randomly and uniformly generating theta (theta)1,θ2,θ3)，θ1,θ2,θ3∈(-π,π)
Figure BDA0001827422200000102
In the formula, theta is a rotation angle around XYZ three axes of the point cloud body coordinate system generated randomly, theta1For rotation angle, theta, about the X-axis of the coordinate system of the point cloud body2For rotation angle, theta, about the Y-axis of the coordinate system of the point cloud body3Is the rotation angle around the Z axis of the point cloud body coordinate system.
Other steps and parameters are the same as those in one of the first to seventh embodiments.
The following examples were used to demonstrate the beneficial effects of the present invention:
the first embodiment is as follows:
simulating the object wind driven generator blade; FIG. 6 shows a blade of a wind turbine;
simulation flow:
point cloud of wind turbine bladesAs a model point cloud NT；
Randomly generated 50 rotation matrixes corresponding to rotation angles of XYZ axes act on wind turbine blade points
Cloud, generating a target point cloud NS。
Repeating the operation for 50 times according to an algorithm, and outputting the output result obtained each time to iter, Corr and; m
RMS and run time were calculated for each time separately. Compared with other algorithms, the following results were obtained.
Figure BDA0001827422200000111
The present invention is capable of other embodiments and its several details are capable of modifications in various obvious respects, all without departing from the spirit and scope of the present invention.

Claims (7)

1. A point cloud registration method based on an extended Gaussian image is characterized in that: the method comprises the following specific processes:
step one, recording a coordinate set of normals in model point cloud and blank point cloud under respective body coordinate system as NTAnd NS；
Step two, respectively carrying out gridding division on the Gaussian expansion image generated by the model point cloud and the blank point cloud, wherein the specific process comprises the following steps:
respectively dividing the Gaussian expansion image generated by the model point cloud and the blank point cloud into 12 basic grids, dividing each grid of the upper level into 4 grids by the HEALPIX layer by layer in a recursion mode until m grids are reached, wherein m is 12 × 4nN is an integer;
calculating the number of normals in each grid by adopting a kd-tree algorithm, and storing the normals in a vector form; simultaneously recording the coordinates of the normal lines under respective body coordinate systems;
according to a nested index mode of HEALPIX, a Gaussian expansion image H generated by blank point cloudSGaussian expansion image H generated by model point cloudTIs divided into 12 parts according to 12 basic grids
Figure FDA0002505886740000011
And
Figure FDA0002505886740000012
wherein
Figure FDA0002505886740000013
The 1 st base mesh of the gaussian expanded image generated for the blank point cloud,
Figure FDA0002505886740000014
the ith base mesh of the gaussian expanded image generated for the blank point cloud,
Figure FDA0002505886740000015
generating a 12 th basic grid of a Gaussian expansion image for the blank point cloud;
Figure FDA0002505886740000016
the 1 st base mesh of the gaussian expanded image generated for the model point cloud,
Figure FDA0002505886740000017
the ith base mesh of the gaussian expanded image generated for the model point cloud,
Figure FDA0002505886740000018
a 12 th base mesh of a gaussian extended image generated for the model point cloud;
step three, based on step two to grid Pi T、
Figure FDA0002505886740000019
And a grid Pi s、
Figure FDA00025058867400000110
Carrying out preliminary matching and matching the preliminarily matched grids Pi T、
Figure FDA00025058867400000111
And a grid Pi s、
Figure FDA00025058867400000112
Sorting in descending order according to the number of normals in the grid, and recording as P1 T″、
Figure FDA00025058867400000123
And P1 s″、
Figure FDA00025058867400000124
Wherein, Pi SIs composed of
Figure FDA00025058867400000113
The grid with the largest number of medium normals; pi TIs composed of
Figure FDA00025058867400000114
The grid with the largest number of medium normals;
Figure FDA00025058867400000115
is composed of
Figure FDA00025058867400000116
The grid with the largest number of medium normals;
Figure FDA00025058867400000117
is composed of
Figure FDA00025058867400000118
The grid with the largest number of medium normals;
Figure FDA00025058867400000119
the ith base mesh of the gaussian expanded image generated for the blank point cloud,
Figure FDA00025058867400000120
generating an ith base mesh of a Gaussian expansion image for the model point cloud;
Figure FDA00025058867400000121
a jth base mesh of a gaussian expanded image generated for the blank point cloud,
Figure FDA00025058867400000122
generating a jth basic grid of a Gaussian expansion image for the model point cloud;
i＝1,2,...,12；j＝i+1,i+2,...,12；i≠j；
step four, P is paired1 T″、
Figure FDA0002505886740000021
And P1 s″、
Figure FDA0002505886740000022
Performing correction to obtain corrected P1 T′、
Figure FDA0002505886740000023
And P1 s′、
Figure FDA0002505886740000024
Step five, according to the corrected P1 T′、
Figure FDA0002505886740000025
P1 s′、
Figure FDA0002505886740000026
Calculating a rotation matrix which can enable the Gaussian expansion image generated by the model point cloud and the Gaussian expansion image generated by the blank point cloud to be registered;
sixthly, outputting iter, Corr and M based on the rotation matrix;
m is a final output rotation matrix;
corr is the correlation between a Gaussian expansion image generated by the blank point cloud and a Gaussian expansion image generated by the model point cloud;
iter is the number of iterations.
2. The point cloud registration method based on the extended Gaussian image as claimed in claim 1, wherein: two pairs of grids P based on the steps in the third stepi T、
Figure FDA0002505886740000027
And a grid Pi s、
Figure FDA0002505886740000028
Carrying out preliminary matching and matching the preliminarily matched grids Pi T、
Figure FDA0002505886740000029
And a grid Pi s、
Figure FDA00025058867400000210
Sorting in descending order according to the number of normals in the grid, and recording as P1 T″、
Figure FDA00025058867400000211
And P1 s″、
Figure FDA00025058867400000212
The specific process is as follows:
step three one, traverse
Figure FDA00025058867400000213
And
Figure FDA00025058867400000214
respectively find
Figure FDA00025058867400000215
And
Figure FDA00025058867400000216
grid formation with the greatest number of normals
Figure FDA00025058867400000217
And
Figure FDA00025058867400000218
wherein P isi SIs composed of
Figure FDA00025058867400000219
The grid with the largest number of medium normals; pi TIs composed of
Figure FDA00025058867400000220
The grid with the largest number of medium normals;
step three, two, pair
Figure FDA00025058867400000221
And
Figure FDA00025058867400000222
sorting in descending order according to the number of the normal lines respectively;
finding a point P starting from the center of the Gaussian expansion imagei s、
Figure FDA00025058867400000223
And Pi T、
Figure FDA00025058867400000224
Is marked as Vi S、
Figure FDA00025058867400000225
And Vi T、
Figure FDA00025058867400000226
Wherein
Figure FDA00025058867400000227
Is composed of
Figure FDA00025058867400000228
The grid with the largest number of medium normals;
Figure FDA00025058867400000229
is composed of
Figure FDA00025058867400000230
The grid with the largest number of medium normals;
by each Pi sJudgment of others
Figure FDA00025058867400000231
Whether or not to satisfy
Figure FDA00025058867400000232
If not, discarding
Figure FDA00025058867400000233
If so, retaining
Figure FDA00025058867400000245
By each Pi TJudgment of others
Figure FDA00025058867400000235
Whether or not to satisfy
Figure FDA00025058867400000236
If not, discarding
Figure FDA00025058867400000237
If so, retaining
Figure FDA00025058867400000244
Step three, calculating the distance between any two grids in the remaining grids
Figure FDA00025058867400000239
And
Figure FDA00025058867400000240
Figure FDA00025058867400000241
a jth base mesh of a gaussian expanded image generated for the blank point cloud,
Figure FDA00025058867400000242
generating a jth basic grid of a Gaussian expansion image for the model point cloud;
Figure FDA0002505886740000031
is PSThe distance between any two of the remaining meshes;
Figure FDA0002505886740000032
is PTThe distance between any two of the remaining meshes;
DSis PSAny two of the remaining gridsA set of distances between the grids;
DTis PTA set of distances between any two of the remaining meshes;
judgment of DSAnd DTWhether there is a similar value in
Figure FDA0002505886740000033
And
Figure FDA0002505886740000034
if it is
Figure FDA0002505886740000035
And
Figure FDA0002505886740000036
if the difference is within 5%, the grid P isi T、
Figure FDA0002505886740000037
And mesh Pi s、
Figure FDA0002505886740000038
Performing primary matching;
mesh P after preliminary matchingi T、
Figure FDA0002505886740000039
And a grid Pi s、
Figure FDA00025058867400000310
Sorting in descending order according to the number of normals in the grid, and recording as P1 T″、
Figure FDA00025058867400000311
And P1 s″、
Figure FDA00025058867400000312
3. The point cloud registration method based on the extended Gaussian image as claimed in claim 2, wherein: the fourth step is to P1 T″、
Figure FDA00025058867400000313
And P1 s″、
Figure FDA00025058867400000314
Performing correction to obtain corrected P1 T′、
Figure FDA00025058867400000315
And P1 s′、
Figure FDA00025058867400000316
The specific process is as follows:
sorting P in descending order according to normal quantity1 T″、
Figure FDA00025058867400000317
And P1 s″、
Figure FDA00025058867400000318
Each of the four grids performs the following steps:
step four, finding the separation P0The most recent n grids, forming a set { P }0,...,PnFind the set { P }0,...,PnThe normal position of each grid in the set constitutes a normal set NP；
P0Is P1 T″、
Figure FDA00025058867400000319
P1 s″、
Figure FDA00025058867400000320
One of (a);
step four, finding a point P which points to the center of the Gaussian expansion image as a starting point0Constructing a plane A perpendicular to V, and the plane A is used for gathering the normal line N by Gaussian expansion of the spherical center of the imagePProjecting to a plane A vertical to V to obtain a normal set NPTwo-dimensional coordinate set N corresponding to each normal in the imagei；
Step three, generating a set PA of uniformly distributed points in the plane A, covering the normal in the plane A by the set PA of points, determining the number of the normal in a circle with a fixed radius by using each point in the set PA as a center by using a KD-tree method, and forming an array HA；
Step four, finding the array HAPoint (x) with the largest number of normals0,y0) Will (x)0,y0) Projecting along the vector V back to the Gaussian expansion image sphere to obtain P0New coordinate value, finish to P0Correcting; (x)0,y0)∈Ni。
4. The point cloud registration method based on the extended Gaussian image as claimed in claim 3, wherein: the fixed radius is:
Figure FDA0002505886740000041
wherein tau is a positive integer; (x)i,yi)∈Ni(ii) a range is distance; x is the number ofiSet N as normalPTwo-dimensional abscissa, y, corresponding to the ith normal lineiSet N as normalPThe ith normal line.
5. The point cloud registration method based on the extended Gaussian image as claimed in claim 4, wherein: according to the corrected P in the step five1 T′、
Figure FDA0002505886740000042
P1 s′、
Figure FDA0002505886740000043
Calculating a rotation matrix which can enable the Gaussian expansion image generated by the model point cloud and the Gaussian expansion image generated by the blank point cloud to be registered; the specific process is as follows:
respectively pointing to P from the center of the sphere of the Gaussian expansion image1 T′、
Figure FDA0002505886740000044
And P1 s′、
Figure FDA0002505886740000045
Is marked as V1 T、
Figure FDA0002505886740000046
And V1 s、
Figure FDA0002505886740000047
Defining a rotation vector V1And V2Angle of rotation theta1And theta2；
Calculating a rotation vector V1Angle of rotation theta1The formula is as follows:
Figure FDA0002505886740000048
Figure FDA0002505886740000049
Figure FDA00025058867400000410
the solving formula is as follows:
Figure FDA00025058867400000411
Figure FDA00025058867400000412
wherein
Figure FDA00025058867400000413
Is an intermediate variable;
initial rotation matrix R1By a rotation vector V1And a rotation angle theta1The formula is obtained as follows:
Figure FDA00025058867400000414
wherein I is an identity matrix and V1 HIs a rotation vector V1Transpose of (V)1xIs a rotation vector V1Component in the x-direction, V1yIs a rotation vector V1Component in the y-direction, V1zIs a rotation vector V1A component in the z direction;
calculating a rotation vector V2Angle of rotation theta2The formula is as follows:
V2＝V1 T
Figure FDA0002505886740000051
Figure FDA0002505886740000052
Figure FDA0002505886740000053
wherein
Figure FDA0002505886740000054
Is an intermediate variable;
initial rotation matrix R2By a rotation vector V2And a rotation angle theta2The formula is obtained as follows:
Figure FDA0002505886740000055
wherein
Figure FDA0002505886740000056
Is a rotation vector V2Transpose of (V)2xIs a rotation vector V2Component in the x-direction, V2yIs a rotation vector V2Component in the y-direction, V2zIs a rotation vector V2A component in the z direction;
according to R1And R2Obtaining a rotation matrix R:
R＝R1R2。
6. the point cloud registration method based on the extended Gaussian image as claimed in claim 5, wherein: in the sixth step, based on the step rotation matrix, iter, Corr and M are output; the specific process is as follows:
step six, defining the iteration times iter as 0 and defining the hysteresis stall as 0; definition M ═ I;
wherein I is an identity matrix; m is a final output rotation matrix;
update Hs＝R×Hs
Calculate HTAnd after update HsThe formula is as follows:
Figure FDA0002505886740000057
sixthly, updating M to be R multiplied by M, and if M is larger than or equal to the threshold value 0.98, outputting iter, Corr and M;
if M is less than the threshold value of 0.98, executing step six and three;
sixthly, recording coordinate sets of normals in the model point cloud and the blank point cloud under respective body coordinate systems as NTAnd Ns′(ii) a Step six or four is executed;
Ns′＝R×Ns
sixthly, re-executing the step two to the step six to obtain new HsAnd HTCalculating new HsAnd HTThe correlation function is as follows:
Figure FDA0002505886740000061
adding 1 to Stall;
executing the step six and five;
sixthly, comparing Corr and Corr':
if Corr × 1.05 is less than Corr', perform Stall decrease by 1; execute the step six
If Corr is less than or equal to Corr ', let Corr equal to Corr', update M ═ R × M, let N bes＝Mrand×Ns(ii) a Adding 1 to Stall; step six is executed;
the M israndA rotation matrix corresponding to the rotation angle theta;
theta is a randomly generated rotation angle around XYZ three axes of the point cloud body coordinate system;
sixthly, judging whether Stall is greater than 3 and Corr' is less than or equal to 0.98, if so, judging that M is equal to Mrand× M, order Ns＝Mrand×Ns(ii) a Sixthly, executing the step;
if not, executing step six and seven;
sixthly, if the stall is 0, iter + 1;
if Corr' is less than the threshold value of 0.98 and iter is less than or equal to 15, repeating the steps sixty-four to sixty-seven, and outputting iter, Corr and M;
if not, iter, Corr and M are output.
7. The point cloud registration method based on the extended Gaussian image as claimed in claim 6, wherein: a rotation matrix M corresponding to the rotation angle thetarandThe solving process is as follows:
randomly and uniformly generating theta (theta)1,θ2,θ3)，θ1,θ2,θ3∈(-π,π)
Figure DEST_PATH_FDA0001827422190000071
In the formula, theta is a rotation angle around XYZ three axes of the point cloud body coordinate system generated randomly, theta1For rotation angle, theta, about the X-axis of the coordinate system of the point cloud body2For rotation angle, theta, about the Y-axis of the coordinate system of the point cloud body3Is the rotation angle around the Z axis of the point cloud body coordinate system.
CN201811190655.4A 2018-10-12 2018-10-12 Point cloud registration method based on extended Gaussian image Active CN109345571B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811190655.4A CN109345571B (en) 2018-10-12 2018-10-12 Point cloud registration method based on extended Gaussian image

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811190655.4A CN109345571B (en) 2018-10-12 2018-10-12 Point cloud registration method based on extended Gaussian image

Publications (2)

Publication Number Publication Date
CN109345571A CN109345571A (en) 2019-02-15
CN109345571B true CN109345571B (en) 2020-10-02

Family

ID=65309481

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811190655.4A Active CN109345571B (en) 2018-10-12 2018-10-12 Point cloud registration method based on extended Gaussian image

Country Status (1)

Country Link
CN (1) CN109345571B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110136182B (en) * 2019-05-28 2021-06-04 北京百度网讯科技有限公司 Registration method, device, equipment and medium for laser point cloud and 2D image
CN113177974A (en) * 2021-05-19 2021-07-27 上海商汤临港智能科技有限公司 Point cloud registration method and device, electronic equipment and storage medium
CN113408688B (en) * 2021-06-29 2022-06-07 哈尔滨工业大学 Unknown environment-oriented multi-radioactive source online searching method

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106204557A (en) * 2016-06-30 2016-12-07 扬州大学 A kind of extracting method of the non-complete data symmetrical feature estimated with M based on extension Gaussian sphere
CN108319567A (en) * 2018-02-05 2018-07-24 北京航空航天大学 A kind of spatial target posture estimation uncertainty calculation method based on Gaussian process

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9355462B2 (en) * 2013-05-08 2016-05-31 Caterpillar Inc. Motion estimation system utilizing point cloud registration

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106204557A (en) * 2016-06-30 2016-12-07 扬州大学 A kind of extracting method of the non-complete data symmetrical feature estimated with M based on extension Gaussian sphere
CN108319567A (en) * 2018-02-05 2018-07-24 北京航空航天大学 A kind of spatial target posture estimation uncertainty calculation method based on Gaussian process

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
HEALPix: A FRAMEWORK FOR HIGH-RESOLUTION DISCRETIZATION AND FAST ANALYSIS OF DATA DISTRIBUTED ON THE SPHERE;Gorski K M，and etc;《 Astrophysical Journal》;20051231;第1-12页 *
Mapping on the HEALPix grid;Calabretta M R，and etc;《Monthly Notices of the Royal Astronomical Society》;20071231;第1-7页 *
基于扩展高斯球的点云数据与CAD模型配准;张学昌等;《机械工程学报》;20070630;第43卷(第6期);第142-148页 *

Also Published As

Publication number Publication date
CN109345571A (en) 2019-02-15

Similar Documents

Publication Publication Date Title
CN110335297B (en) Point cloud registration method based on feature extraction
CN111299815B (en) Visual detection and laser cutting trajectory planning method for low-gray rubber pad
CN109559338B (en) Three-dimensional point cloud registration method based on weighted principal component analysis method and M estimation
CN109345571B (en) Point cloud registration method based on extended Gaussian image
Ji et al. A novel simplification method for 3D geometric point cloud based on the importance of point
CN107610061B (en) Edge-preserving point cloud hole repairing method based on two-dimensional projection
CN110807781B (en) Point cloud simplifying method for retaining details and boundary characteristics
CN106373118A (en) A complex curved surface part point cloud reduction method capable of effectively keeping boundary and local features
CN103793712A (en) Image recognition method and system based on edge geometric features
CN111145228A (en) Heterogeneous image registration method based on local contour point and shape feature fusion
CN113628263A (en) Point cloud registration method based on local curvature and neighbor characteristics thereof
CN110838115A (en) Ancient cultural relic three-dimensional model change detection method by contour line extraction and four-dimensional surface fitting
CN113034554B (en) Whale optimized broken warrior body fragment registration method based on chaos reverse learning
CN114663373A (en) Point cloud registration method and device for detecting surface quality of part
CN103778619A (en) Image matching method based on Zernike matrix
CN111710023A (en) Three-dimensional point cloud data feature point extraction method and application
Sun et al. Oriented point sampling for plane detection in unorganized point clouds
CN113706588B (en) Annular forging point cloud registration method based on improved four-point quick robust matching algorithm
CN114332172A (en) Improved laser point cloud registration method based on covariance matrix
CN108010114B (en) Geometric shape recognition method and feature recognition method for basic primitive point cloud curved surface
CN110310322A (en) Method for detecting assembly surface of 10-micron-level high-precision device
Liu et al. A fast weighted registration method of 3d point cloud based on curvature feature
CN106023314A (en) B spline master curve fitting method based on rotary axis direction mapping
CN110580497A (en) Spatial scene matching method based on rotation invariance
CN110533781B (en) Automatic labeling method for multi-class three-dimensional model components

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant