CN109345571A - A kind of point cloud registration method based on extension Gaussian image - Google Patents
A kind of point cloud registration method based on extension Gaussian image Download PDFInfo
- Publication number
- CN109345571A CN109345571A CN201811190655.4A CN201811190655A CN109345571A CN 109345571 A CN109345571 A CN 109345571A CN 201811190655 A CN201811190655 A CN 201811190655A CN 109345571 A CN109345571 A CN 109345571A
- Authority
- CN
- China
- Prior art keywords
- grid
- point cloud
- gaussian
- image
- normal
- 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.)
- Granted
Links
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
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Processing Or Creating Images (AREA)
- Image Analysis (AREA)
Abstract
A kind of point cloud registration method based on extension Gaussian image, the present invention relates to the point cloud registration methods of image.The purpose of the present invention is to solve the common rough registration algorithm SAC-IA algorithm expense in existing point cloud registering field is excessive, Gaussian sphere is divided into polyhedron or longitude and latitude based on EGI algorithm, spend the time longer, and optimal solution is directly searched in feasible zone when based on EGI algorithm calculating spin matrix, it is difficult to the problem of the double requirements of balance efficiency and precision.Process are as follows: coordinate set of the normal under respective body coordinate system is denoted as NTAnd NS；Gridding point is carried out to Gaussian spread image；By the grid after preliminary matches according to grid inter normal quantity descending sort；Obtain revised grid；Calculate spin matrix；Export iter, Corr and M.The present invention is used for the crossing domain of computer graphics and mechanical manufacturing technology.
Description
Technical field
The invention belongs to the crossing domains of computer graphics and mechanical manufacturing technology, and in particular to the point cloud registering of image
Method.
Background technique
Blank manufacture first usually is carried out using casting in workpiece process, subtracts material using turning, grinding etc. later and processes
Part is processed into final shape by mode.It needs to detect casting situation during this, determines the distribution of surplus, later
The distribution of surplus is advanced optimized, next step machining is convenient for.Laser scanner because its convenient for surface point cloud data it is rapid, when
Between at low cost, the higher characteristic of precision, compared with the detection being manually more suitable for casting.The point cloud number of component is got in scanner
According to later, using point cloud registration algorithm, the model designed in the cloud and CAD software is registrated and is compared, this method
It is capable of the surplus of efficient, reliable analysis machining later.
Point cloud registering field is directed to various different scenes and the point cloud registering of different shapes has done many correlative studys.Point
Cloud is registrated the common rough registration algorithm SAC-IA in field, by way of iteration, minimizes between target point cloud and model point cloud
Distance, complete registration.The algorithm expense is excessive.There is document to mention to estimate using Extended Gaussian Images (EGI)
A kind of algorithm that ICP algorithm completes registration is used in combination in meter spin matrix, Gaussian sphere is divided into polyhedron or longitude and latitude, each
Face area is identical, however the specific coordinate for obtaining each polyhedral corner angle is more difficult, while when normal quantity is more huge
When big, need one by one face judge to calculate whether each normal is directed toward the face, spend the time longer.It calculates direct when spin matrix
Optimal solution is searched in feasible zone, it is difficult to the double requirements of balance efficiency and precision.
EGI is extension Gaussian image；
Summary of the invention
The purpose of the present invention is to solve the common rough registration algorithm SAC-IA algorithm expenses in existing point cloud registering field
It is excessive, Gaussian sphere is divided by polyhedron or longitude and latitude based on EGI algorithm, spends the time longer, and calculate based on EGI algorithm
Optimal solution is directly searched for when spin matrix in feasible zone, it is difficult to the problem of the double requirements of balance efficiency and precision, and propose
A kind of point cloud registration method based on extension Gaussian image.
A kind of point cloud registration method detailed process based on extension Gaussian image are as follows:
Step 1: coordinate set of the normal in model point cloud and blank point cloud under respective body coordinate system is denoted as NT
And NS(obtaining normal position)；
Step 2: the Gaussian spread image (EGI) generated respectively to model point cloud and blank point cloud carries out gridding minute；
Step 3: based on step 2 to grid Pi T、With grid Pi s、Preliminary matches are carried out, by the net after preliminary matches
Lattice Pi T、With grid Pi s、According to grid inter normal quantity descending sort, it is denoted as P1 T″、And P1 s″、
Wherein, Pi SForThe most grid of middle normal quantity；Pi TForThe most grid of middle normal quantity；
ForThe most grid of middle normal quantity；ForThe most grid of middle normal quantity；
For blank point cloud generate Gaussian spread image i-th of basic grid,The height generated for model point cloud
I-th of basic grid of this expanded images,
For blank point cloud generate Gaussian spread image j-th of basic grid,The height generated for model point cloud
J-th of basic grid of this expanded images；
I=1,2 ..., 12；J=i+1, i+2 ..., 12；i≠j；
Step 4: to P1 T″、And P1 s″、It is modified, obtains revised P1 T′、And P1 s′、
Step 5: according to revised P1 T′、And P1 s′、Calculate the Gaussian spread image that model point cloud can be made to generate
(EGI) and blank point cloud generate Gaussian spread image (EGI) registration spin matrix；
Step 6: being based on step spin matrix, iter, Corr and M are exported；
M is the spin matrix of final output；
Corr is the correlation for the Gaussian spread image that the Gaussian spread image that blank point cloud generates and model point cloud generate；
Iter is the number of iterations.
The invention has the benefit that
The algorithm for point cloud registering in machining occasion surplus analysis is proposed, is known in GRIDDING WITH WEIGHTED AVERAGE, corresponding points
Not, spin matrix calculates and the parts such as iterative strategy, the shortcomings that for point cloud registration algorithm based on EGI in the past, has carried out one
Determine the improvement of degree.It is significantly lower than previous algorithm on computation complexity, reduces algorithm expense.Use laser scanning reality
Under the background of existing Precision Machining surplus analysis, by using HEALPix algorithm, changes the division mode of spherical surface, construct KDtree
It determines the normal that each face includes, improves division efficiency, spend the time shorter.It segments dough sheet further to improve precision, adopts simultaneously
With the calculation of " moving average ", prevent meticulous division from obscuring peak value.Using iteration plus variation, essence is on the one hand improved
Degree, on the other hand avoids suboptimization, the double requirements of balance efficiency and precision.
Gaussian sphere is divided into polyhedron based on EGI algorithm or longitude and latitude is changed to change spherical surface using HEALPix algorithm
Division mode；Optimal solution is searched for when EGI algorithm calculates spin matrix directly in feasible zone to be changed to propose clear calculation formula meter
Calculate spin matrix and iterative strategy.Inventive algorithm expense about reduces the calculating time of 80-85% compared with SAC-IA algorithm.This hair
Bright precision is with RMS (root mean square) metering, about raising an order of magnitude.
Detailed description of the invention
Fig. 1 is flow chart of the present invention；
Fig. 2 is that the present invention passes through 2HEALPix method to Gaussian spread image segmentation into 12 grid schematic diagrames；
Fig. 3 is that the present invention passes through 2HEALPix method to Gaussian spread image segmentation into 48 grid schematic diagrames；
Fig. 4 is that the present invention passes through 2HEALPix method to Gaussian spread image segmentation into 192 grid schematic diagrames；
Fig. 5 is that the present invention passes through 2HEALPix method to Gaussian spread image segmentation into 768 grid schematic diagrames；
Fig. 6 is wind electricity blade schematic diagram.
Specific embodiment
Specific embodiment 1: embodiment is described with reference to Fig. 1, one kind of present embodiment is based on extension Gaussian image
Point cloud registration method detailed process are as follows:
Step 1: coordinate set of the normal in model point cloud and blank point cloud under respective body coordinate system is denoted as NT
And NS(obtaining normal position)；
The model point cloud is the point cloud of ideal machine part model；
Blank point cloud is the point cloud (blank that factory process goes out, what scanner scanning obtained) of blank；
Step 2: the Gaussian spread image (EGI) generated respectively to model point cloud and blank point cloud carries out gridding minute；
Step 3: based on step 2 to grid Pi T、With grid Pi s、Preliminary matches are carried out, after preliminary matches
Grid Pi T、With grid Pi s、According to grid inter normal quantity descending sort, it is denoted as P1 T″、And P1 s″、
Wherein, Pi SForThe most grid of middle normal quantity；Pi TForThe most grid of middle normal quantity；
ForThe most grid of middle normal quantity；ForThe most grid of middle normal quantity；
For blank point cloud generate Gaussian spread image i-th of basic grid,The Gauss generated for model point cloud
I-th of basic grid of expanded images,
For blank point cloud generate Gaussian spread image j-th of basic grid,The height generated for model point cloud
J-th of basic grid of this expanded images；
I=1,2 ..., 12；J=i+1, i+2 ..., 12；i≠j；
Step 4: to P1 T″、And P1 s″、It is modified, obtains revised P1 T′、And P1 s′、
Step 5: according to revised P1 T′、P1 s′、Calculate the Gaussian spread image that model point cloud can be made to generate
(EGI) and blank point cloud generate Gaussian spread image (EGI) registration spin matrix；
Step 6: being based on step spin matrix, iter, Corr and M are exported；
M is the spin matrix of final output；
Corr is the correlation for the Gaussian spread image that the Gaussian spread image that blank point cloud generates and model point cloud generate；
Iter is the number of iterations.
Specific embodiment 2: the present embodiment is different from the first embodiment in that: it is right respectively in the step 2
The Gaussian spread image (EGI) that model point cloud and blank point cloud generate carries out gridding minute, detailed process are as follows:
HEALPix (Hierarchical Equal Area isoLatitude Pixelization) is initially applied to
Research in terms of cosmic microwave background, to support the high-resolution discretization of spherical surface.HEALPix can be divided into number according to resolution ratio
Grade, lowermost level are 12 grids.HEALPix divides spherical surface as shown in Fig. 2,3,4,5.
HEALPix has the feature that
Grid block is layered construction, and data are easily accessible.
Generation grid block method is simple, and operation is efficient.
Spherical surface such as is divided at the curved surface quadrangle of sizes, and the area of each grid is identical.
The Gaussian spread image that model point cloud and blank point cloud generate is divided into 12 basic grids, HEALPix respectively
(HEALPix is an acronym for Hierarchical Equal Area isoLatitudePixelization of
a sphere.As suggested in the name,this pixelization produces a subdivision of
A spherical surface in which each pixel covers the same surface area.) by passing
The layer-by-layer each grid dividing to upper level of the mode returned is 4, until reaching m grid, m=12 × 4n, n value is whole
Number；M=12 × 4n=3072；
Calculating the normal quantity in each grid using kd-tree (k-d tree) algorithm, (normal quantity is model midpoint
What quantity determined, its primary corresponding normal is calculated to each point), being stored as vector, (each element represents the gridding method
The quantity of line) form；Coordinate of the normal under respective body coordinate system is recorded simultaneously；
Index is established for each grid with by the way of " nested " using " ring-type ", convenient for quickly accessing by different modes
Data and data are handled.
According to the nested-attribute index mode of HEALPix, the Gaussian spread image H that blank point cloud is generatedSIt is raw with model point cloud
At Gaussian spread image HTBeing divided into 12 parts according to 12 basic grids, (each part has 4n(256) normal of a grid
Quantity), respectivelyWithI=1,2 ..., 12；
WhereinFor blank point cloud generate Gaussian spread image the 1st basic grid,For the generation of blank point cloud
Gaussian spread image i-th of basic grid,For the 12nd facilities network of the Gaussian spread image that blank point cloud generates
Lattice；For model point cloud generate Gaussian spread image the 1st basic grid,The Gaussian spread generated for model point cloud
I-th of basic grid of image,For the 12nd basic grid of the Gaussian spread image that model point cloud generates.
Other steps and parameter are same as the specific embodiment one.
Specific embodiment 3: the present embodiment is different from the first and the second embodiment in that: base in the step 3
In step 2 to grid Pi T、With grid Pi s、Preliminary matches are carried out, by the grid P after preliminary matchesi T、And grid
Pi s、According to grid inter normal quantity descending sort, it is denoted as P1 T″、And P1 s″、Detailed process are as follows:
Corresponding points identification
Finding corresponding grid may be a challenging job, especially when the normal of object is to a certain degree
On when being uniformly distributed.Due to the characteristic distributions of normal vector, a highdensity grid may have similar density in the reverse direction
Grid, thus spin matrix calculating in generate error.In addition, good algorithm should change not the rotation and position of grid
It is sensitive.Herein, we are closed using a kind of simple and effective method according to the length between the quantity of normal and selected lattice
System, finds the grid in target EGI relevant to grid in the EGI of source.Corresponding points identification, carries out according to following step:
Step 3 one, traversalWithIt finds respectivelyWithMiddle method
The most grid of line number amount is constitutedWithI=1,2 ..., 12；
Wherein Pi SForThe most grid of middle normal quantity；Pi TForThe most grid of middle normal quantity；
It is step 3 two, rightWithRespectively according to normal quantity
Carry out descending sort；
It finds and is directed toward P by starting point of the Gaussian spread image centre of spherei s、And Pi T、Vector be denoted as Vi S、And Vi T、
WhereinForThe most grid of middle normal quantity；ForThe most grid of middle normal quantity；
With each Pi sJudge otherWhether meetIf conditions are not met, giving upSuch as
Fruit meets, and retainsJ=i+1, i+2 ..., 12；i≠j；
With each Pi TJudge otherWhether meetIf conditions are not met, giving upSuch as
Fruit meets, and retainsJ=i+1, i+2 ..., 12；i≠j；
Step 3 three the distance between calculates any two grid in remaining grid
With
For blank point cloud generate Gaussian spread image j-th of basic grid,The height generated for model point cloud
J-th of basic grid of this expanded images；
For PSIn the distance between any two grid in remaining grid；
For PTIn the distance between any two grid in remaining grid；
DSFor PSIn the distance between any two grid in remaining grid set；
DTFor PTIn the distance between any two grid in remaining grid set；
Judge DSAnd DTIn whether have similar valueWithIfWithDifference is within 5%, then grid Pi T、With
Grid Pi s、Preliminary matches；
By test it can be found that this method can effectively exclude the point off density of hypotelorism, while
Higher two points off density of normal direction density are identified by their geometric distance in two EGI, this method is to rotation
Turn and insensitive for noise.
By the grid P after preliminary matchesi T、With grid Pi s、According to grid inter normal quantity descending sort, it is denoted as
P1 T″、And P1 s″、
Other steps and parameter are the same as one or two specific embodiments.
Specific embodiment 4: unlike one of present embodiment and specific embodiment one to three: the step 4
In to P1 T″、And P1 s″、It is modified, obtains revised P1 T′、And P1 s′、Detailed process are as follows:
The coordinate optimizing of Chosen Point
By taking 3072 grids that HEALPix is generated as an example, even if large number of, the grid element center that every two is closed on of segmentation
It is about 4 degree to spherical angle, error is still larger.Notice that the high grid of a density may be with adjacent grid cluster simultaneously
Form the intensive region of normal.After finding out close quarters pairs of in two EGI.We can be by with maximum close
Degree is the central point of criterion calculation whole region to improve precision, rather than uses the geometry of the grid of highest scoring in the region
Central point, the two central points are usually unequal.Then show we it is necessary to further increase precision, in mention simple improve and draw
The mesh-density divided may be reduced the normal quantity that each grid of normal close quarters includes, it is difficult to be identified, be unfavorable for
Later period registration.Here the method for using consecutive mean, that is, fix the size of each " window ", and changes the position of window center
It sets, intuitively sees and moved such as window in discovery close quarters with very small step-length, the size of window is fixed, and is guaranteed
Always there are enough normals to fall in window, and the normal fallen in window can be made until reaching one by changing the window's position
Until quantity maximizes.Specific algorithm steps are as follows.
It will be by the P after the descending sort of normal quantity1 T″、And P1 s″、Each in four grids is performed both by following step
It is rapid:
Step 4 one is found from P0N nearest grid, composition set { P0,...,Pn, find set { P0,...,Pn}
In each grid normal position (knowing that normal quantity is known that normal position in front), constitute normal collection NP；
P0For P1 T″、P1 s″、In one；
Step 4 two finds and is directed toward P by starting point of the Gaussian spread image centre of sphere0Vector V, construct perpendicular to V plane
A, and plane A passes through the Gaussian spread image centre of sphere, by normal collection NPProjection obtains normal collection N to the plane A perpendicular to VPIn it is every
The corresponding two-dimensional coordinate set N of one normali；
Step 4 three, the set PA (many grids) that equally distributed point is generated in plane A, the set PA covering of point
Normal in plane A is determined using KD-tree method using point each in set PA as the method line number in center radii fixus circle
Amount forms array HA；
Step 4 four finds array HAMost point (the x of middle normal quantity0,y0), by (x0,y0) along vector V it is projected back in Gauss
Expanded images spherical surface, obtains P0New coordinate value is completed to P0Amendment；(x0,y0)∈Ni。
Other steps and parameter are identical as one of specific embodiment one to three.
Specific embodiment 5: unlike one of present embodiment and specific embodiment one to four: the fixation half
Diameter are as follows:
Wherein, τ value is positive integer (for controlling round size)；(xi,yi)∈Ni；Range is distance；xiFor normal
Collect NPIn the corresponding two-dimentional abscissa of i-th of normal, yiFor normal collection NPIn the corresponding two-dimentional ordinate of i-th of normal.
Other steps and parameter are identical as one of specific embodiment one to four.
Specific embodiment 6: unlike one of present embodiment and specific embodiment one to five: the step 5
It is middle according to revised P1 T′、P1 s′、Calculate the Gaussian spread image (EGI) and blank point that model point cloud can be made to generate
The spin matrix for Gaussian spread image (EGI) registration that cloud generates；Detailed process are as follows:
It is aligned higher score grid first, then keeps the two Grid Aligns, the grid being previously aligned along direction
Axis.Two EGI are rotated, until lower two Grid Aligns of score.Two rotating vectors are expressed as a spin moment
Battle array, by the alignment that each point of this matrix application in target point cloud is completed to two clouds.Specific formula is as follows:
P will be respectively directed to using the Gaussian spread image centre of sphere as starting point1 T′、And P1 s′、Vector be denoted as V1 T、With
V1 s、Define rotating vector V1And V2, rotation angle θ1And θ2；
Calculate rotating vector V1, rotation angle θ1Formula is as follows:
Solution formula are as follows:
WhereinFor intermediate variable；
Initial rotation vector R1By rotating vector V1And rotation angle θ1It obtains, formula are as follows:
Wherein I is unit matrix, V1 HFor rotating vector V1Transposition, V1xFor rotating vector V1Component in the x direction, V1y
For rotating vector V1Component in y-direction, V1zFor rotating vector V1Component in a z-direction；
Calculate rotating vector V2, rotation angle θ2Formula is as follows:
V2=V1 T
WhereinFor intermediate variable；
Initial rotation vector R2By rotating vector V2And rotation angle θ2It obtains, formula are as follows:
WhereinFor rotating vector V2Transposition, V2xFor rotating vector V2Component in the x direction, V2yFor rotating vector V2
Component in y-direction, V2zFor rotating vector V2Component in a z-direction；
According to R1And R2Obtain spin matrix R:
R=R1R2。
Other steps and parameter are identical as one of specific embodiment one to five.
Specific embodiment 7: unlike one of present embodiment and specific embodiment one to six: the step 6
In be based on step spin matrix, export iter, Corr and M；Detailed process are as follows:
Step 6 one defines the number of iterations iter=0, defines sluggishness stall=0；Define M=I；
Wherein, I is unit matrix；M is the spin matrix of final output；
Update Hs=R × Hs
Calculate HTWith H after updatesCorrelation, formula is as follows:
When degree of correlation function is the degree of correlation maximum that 1 is two EGI images.
Step 6 two updates M=R × M, if M is more than or equal to threshold value 0.98, exports iter, Corr and M；
If M is less than threshold value 0.98, step 6 three is executed；
Coordinate set of the normal in model point cloud and blank point cloud under respective body coordinate system is denoted as by step 6 three
NTAnd Ns′(obtaining normal position)；Execute step 6 four；
Ns′=R × Ns
Step 6 four re-execute the steps two to step 6 one, obtains new HsAnd HT, calculate new HsAnd HTCorrelation
Property, degree of correlation function are for example as follows:
Stall is enabled to add 1；
Execute step 6 five；
Step 6 five compares Corr and Corr ':
If Corr × 1.05 are less than Corr ', execute Stall and subtract 1；Execute step 6 six
If Corr is less than or equal to Corr ', Corr is allowed to be equal to Corr ', updates M=R × M；Enable Ns=Mrand×Ns；It enables
Stall adds 1；Execute step 6 six
(Corr is greater than Corr ' and does not consider)；
The MrandFor the corresponding spin matrix of rotation angle θ；
θ is to be generated at random around a rotation angle for tri- axis of cloud body coordinate system XYZ；
Step 6 six judges whether Stall is greater than 3, and Corr ' is less than or equal to 0.98, if it is satisfied, M=Mrand×M；
Enable Ns=Mrand×Ns；Execute step 6 seven；
If conditions are not met, executing step 6 seven；
Step 6 seven, another stall=0, iter+1；
If meeting Corr ' less than threshold value 0.98, and iter is less than or equal to 15, repeats step 6 four to step 6
Seven, export iter, Corr and M；
If conditions are not met, output iter, Corr and M.
Other steps and parameter are identical as one of specific embodiment one to six.
Specific embodiment 8: unlike one of present embodiment and specific embodiment one to seven: the rotation angle θ
Corresponding spin matrix MrandSolution procedure are as follows:
It is random uniformly to generate θ=(θ1,θ2,θ3), θ1,θ2,θ3∈(-π,π)
In formula, θ is to be generated at random around a rotation angle for tri- axis of cloud body coordinate system XYZ, θ1For around a cloud ontology coordinate
It is the rotation angle of X-axis, θ2For around a rotation angle for cloud body coordinate system Y-axis, θ3For around a cloud body coordinate system Z axis rotation angle.
Other steps and parameter are identical as one of specific embodiment one to seven.
Beneficial effects of the present invention are verified using following embodiment:
Embodiment one:
Simulation object blade of wind-driven generator；If Fig. 6 is blade of wind-driven generator；
Simulation flow:
Using blade of wind-driven generator point cloud as model point cloud NT；
50 be randomly generated act on blade of wind-driven generator point to the corresponding spin matrix of the rotation angle of XYZ axis
Cloud generates target point cloud NS。
According to algorithm, rerun 50 times, by the output result obtained every time export iter, Corr and；M
Calculate separately RMS each time and runing time.With other algorithm comparisons, following result is obtained.
The present invention can also have other various embodiments, without deviating from the spirit and substance of the present invention, this field
Technical staff makes various corresponding changes and modifications in accordance with the present invention, but these corresponding changes and modifications all should belong to
The protection scope of the appended claims of the present invention.
Claims (8)
1. a kind of point cloud registration method based on extension Gaussian image, it is characterised in that: the method detailed process are as follows:
Step 1: coordinate set of the normal in model point cloud and blank point cloud under respective body coordinate system is denoted as NTAnd NS；
Step 2: the Gaussian spread image generated respectively to model point cloud and blank point cloud carries out gridding minute；
Step 3: based on step 2 to grid Pi T、With grid Pi s、Preliminary matches are carried out, by the grid after preliminary matches
Pi T、With grid Pi s、According to grid inter normal quantity descending sort, it is denoted as P1 T″、And P1 s″、
Wherein, Pi SForThe most grid of middle normal quantity；Pi TForThe most grid of middle normal quantity；
ForThe most grid of middle normal quantity；ForThe most grid of middle normal quantity；
For blank point cloud generate Gaussian spread image i-th of basic grid,The Gaussian spread generated for model point cloud
I-th of basic grid of image；
For blank point cloud generate Gaussian spread image j-th of basic grid,Expand for the Gauss that model point cloud generates
Open up j-th of basic grid of image；
I=1,2 ..., 12；J=i+1, i+2 ..., 12；i≠j；
Step 4: to P1 T″、And P1 s″、It is modified, obtains revised P1 T′、And P1 s′、
Step 5: according to revised P1 T′、P1 s′、Calculate the Gaussian spread image and blank that model point cloud can be made to generate
The spin matrix for the Gaussian spread image registration that point cloud generates；
Step 6: being based on spin matrix, iter, Corr and M are exported；
M is the spin matrix of final output；
Corr is the correlation for the Gaussian spread image that the Gaussian spread image that blank point cloud generates and model point cloud generate；
Iter is the number of iterations.
2. a kind of point cloud registration method based on extension Gaussian image according to claim 1, it is characterised in that: the step
The Gaussian spread image generated respectively to model point cloud and blank point cloud in two carries out gridding minute, detailed process are as follows:
The Gaussian spread image that model point cloud and blank point cloud generate is divided into 12 basic grids respectively, HEALPix is by passing
The layer-by-layer each grid dividing to upper level of the mode returned is 4, until reaching m grid, m=12 × 4n, n value is whole
Number；
Normal quantity in each grid is calculated using kd-tree algorithm, is stored as vector form；Normal is recorded each simultaneously
Coordinate under from body coordinate system；
According to the nested-attribute index mode of HEALPix, the Gaussian spread image H that blank point cloud is generatedSThe height generated with model point cloud
This expanded images HT12 parts are divided into according to 12 basic grids, respectivelyWith
WhereinFor blank point cloud generate Gaussian spread image the 1st basic grid,The height generated for blank point cloud
I-th of basic grid of this expanded images,For the 12nd basic grid of the Gaussian spread image that blank point cloud generates；H1 T
For model point cloud generate Gaussian spread image the 1st basic grid,The Gaussian spread image generated for model point cloud
I-th of basic grid,For the 12nd basic grid of the Gaussian spread image that model point cloud generates.
3. a kind of point cloud registration method based on extension Gaussian image according to claim 1 or claim 2, it is characterised in that: described
Based on step 2 to grid P in step 3i T、With grid Pi s、Preliminary matches are carried out, by the grid P after preliminary matchesi T、With grid Pi s、According to grid inter normal quantity descending sort, it is denoted as P1 T″、And P1 s″、Detailed process are as follows:
Step 3 one, traversalWithIt finds respectivelyWithMiddle method line number
Most grids is measured to constituteWith
Wherein Pi SForThe most grid of middle normal quantity；Pi TForThe most grid of middle normal quantity；
It is step 3 two, rightWithIt is carried out respectively according to normal quantity
Descending sort；
It finds and is directed toward P by starting point of the Gaussian spread image centre of spherei s、And Pi T、Vector be denoted as Vi S、And Vi T、
WhereinForThe most grid of middle normal quantity；ForThe most grid of middle normal quantity；
With each Pi sJudge otherWhether meetIf conditions are not met, giving upIf full
Foot retains
With each Pi TJudge otherWhether meetIf conditions are not met, giving upIf full
Foot retains
Step 3 three the distance between calculates any two grid in remaining gridWith
For blank point cloud generate Gaussian spread image j-th of basic grid,The Gaussian spread generated for model point cloud
J-th of basic grid of image；
For PSIn the distance between any two grid in remaining grid；
For PTIn the distance between any two grid in remaining grid；
DSFor PSIn the distance between any two grid in remaining grid set；
DTFor PTIn the distance between any two grid in remaining grid set；
Judge DSAnd DTIn whether have similar valueWithIfWithDifference is within 5%, then grid Pi T、With grid
Pi s、Preliminary matches；
By the grid P after preliminary matchesi T、With grid Pi s、According to grid inter normal quantity descending sort, it is denoted as P1 T″、
And P1 s″、
4. a kind of point cloud registration method based on extension Gaussian image according to claim 3, it is characterised in that: the step
To P in four1 T″、And P1 s″、It is modified, obtains revised P1 T′、And P1 s′、Detailed process are as follows:
It will be by the P after the descending sort of normal quantity1 T″、And P1 s″、Each in four grids is performed both by following steps:
Step 4 one is found from P0N nearest grid, composition set { P0,...,Pn, find set { P0,...,PnIn it is every
The normal position of one grid constitutes normal collection NP；
P0For P1 T″、P1 s″、In one；
Step 4 two finds and is directed toward P by starting point of the Gaussian spread image centre of sphere0Vector V, construct the plane A perpendicular to V, and flat
Face A is by the Gaussian spread image centre of sphere, by normal collection NPProjection obtains normal collection N to the plane A perpendicular to VPIn each method
The corresponding two-dimensional coordinate set N of linei；
Step 4 three, the set PA that equally distributed point is generated in plane A, the normal in the set PA overlay planes A of point, make
Determine that point each using in set PA as the normal quantity in center radii fixus circle, forms array H with KD-tree methodA；
Step 4 four finds array HAMost point (the x of middle normal quantity0,y0), by (x0,y0) along vector V it is projected back in Gaussian spread
Image spherical surface, obtains P0New coordinate value is completed to P0Amendment；(x0,y0)∈Ni。
5. a kind of point cloud registration method based on extension Gaussian image according to claim 4, it is characterised in that: the fixation
Radius are as follows:
Wherein, τ value is positive integer；(xi,yi)∈Ni；Range is distance；xiFor normal collection NPIn i-th of normal corresponding two
Tie up abscissa, yiFor normal collection NPIn the corresponding two-dimentional ordinate of i-th of normal.
6. a kind of point cloud registration method based on extension Gaussian image according to claim 5, it is characterised in that: the step
According to revised P in five1 T′、P1 s′、It calculates the Gaussian spread image that model point cloud can be made to generate and blank point cloud is raw
At Gaussian spread image registration spin matrix；Detailed process are as follows:
P will be respectively directed to using the Gaussian spread image centre of sphere as starting point1 T′、And P1 s′、Vector be denoted as V1 T、And V1 s、Define rotating vector V1And V2, rotation angle θ1And θ2；
Calculate rotating vector V1, rotation angle θ1Formula is as follows:
Solution formula are as follows:
WhereinFor intermediate variable；
Initial rotation vector R1By rotating vector V1And rotation angle θ1It obtains, formula are as follows:
Wherein, I is unit matrix, V1 HFor rotating vector V1Transposition, V1xFor rotating vector V1Component in the x direction, V1yFor
Rotating vector V1Component in y-direction, V1zFor rotating vector V1Component in a z-direction；
Calculate rotating vector V2, rotation angle θ2Formula is as follows:
V2=V1 T
WhereinFor intermediate variable；
Initial rotation vector R2By rotating vector V2And rotation angle θ2It obtains, formula are as follows:
WhereinFor rotating vector V2Transposition, V2xFor rotating vector V2Component in the x direction, V2yFor rotating vector V2In y
Component on direction, V2zFor rotating vector V2Component in a z-direction；
According to R1And R2Obtain spin matrix R:
R=R1R2。
7. a kind of point cloud registration method based on extension Gaussian image according to claim 6, it is characterised in that: the step
It is based on step spin matrix in six, exports iter, Corr and M；Detailed process are as follows:
Step 6 one defines the number of iterations iter=0, defines sluggishness stall=0；Define M=I；
Wherein, I is unit matrix；M is the spin matrix of final output；
Update Hs=R × Hs
Calculate HTWith H after updatesCorrelation, formula is as follows:
Step 6 two updates M=R × M, if M is more than or equal to threshold value 0.98, exports iter, Corr and M；
If M is less than threshold value 0.98, step 6 three is executed；
Coordinate set of the normal in model point cloud and blank point cloud under respective body coordinate system is denoted as N by step 6 threeTWith
Ns′；Execute step 6 four；
Ns′=R × Ns
Step 6 four re-execute the steps two to step 6 one, obtains new HsAnd HT, calculate new HsAnd HTCorrelation, phase
Pass degree function is for example as follows:
Stall is enabled to add 1；
Execute step 6 five；
Step 6 five compares Corr and Corr ':
If Corr × 1.05 are less than Corr ', execute Stall and subtract 1；Execute step 6 six
If Corr is less than or equal to Corr ', Corr is allowed to be equal to Corr ', updates M=R × M；Enable Ns=Mrand×Ns；Stall is enabled to add
1；Execute step 6 six；
The MrandFor the corresponding spin matrix of rotation angle θ；
θ is to be generated at random around a rotation angle for tri- axis of cloud body coordinate system XYZ；
Step 6 six judges whether Stall is greater than 3, and Corr ' is less than or equal to 0.98, if it is satisfied, M=Mrand×M；Enable Ns=
Mrand×Ns；Execute step 6 seven；
If conditions are not met, executing step 6 seven；
Step 6 seven, another stall=0, iter+1；
If meeting Corr ' less than threshold value 0.98, and iter is less than or equal to 15, step 6 four is repeated to step 6 seven, it is defeated
Iter, Corr and M out；
If conditions are not met, output iter, Corr and M.
8. a kind of point cloud registration method based on extension Gaussian image according to claim 7, it is characterised in that: the rotation
The corresponding spin matrix M of angle θrandSolution procedure are as follows:
It is random uniformly to generate θ=(θ1,θ2,θ3), θ1,θ2,θ3∈(-π,π)
In formula, θ is to be generated at random around a rotation angle for tri- axis of cloud body coordinate system XYZ, θ1For around a cloud body coordinate system X-axis
Rotation angle, θ2For around a rotation angle for cloud body coordinate system Y-axis, θ3For around a cloud body coordinate system Z axis rotation angle.
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 true CN109345571A (en) | 2019-02-15 |
| CN109345571B 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) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110136182A (en) * | 2019-05-28 | 2019-08-16 | 北京百度网讯科技有限公司 | Method for registering, device, equipment and the medium of 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 |
| CN113408688A (en) * | 2021-06-29 | 2021-09-17 | 哈尔滨工业大学 | Unknown environment-oriented multi-radioactive source online searching method |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20140334685A1 (en) * | 2013-05-08 | 2014-11-13 | Caterpillar Inc. | Motion estimation system utilizing point cloud registration |
| 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 |
-
2018
- 2018-10-12 CN CN201811190655.4A patent/CN109345571B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20140334685A1 (en) * | 2013-05-08 | 2014-11-13 | Caterpillar Inc. | Motion estimation system utilizing point cloud registration |
| 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)
| Title |
|---|
| CALABRETTA M R，AND ETC: "Mapping on the HEALPix grid", 《MONTHLY NOTICES OF THE ROYAL ASTRONOMICAL SOCIETY》 * |
| GORSKI K M，AND ETC: "HEALPix: A FRAMEWORK FOR HIGH-RESOLUTION DISCRETIZATION AND FAST ANALYSIS OF DATA DISTRIBUTED ON THE SPHERE", 《 ASTROPHYSICAL JOURNAL》 * |
| 张学昌等: "基于扩展高斯球的点云数据与CAD模型配准", 《机械工程学报》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110136182A (en) * | 2019-05-28 | 2019-08-16 | 北京百度网讯科技有限公司 | Method for registering, device, equipment and the medium of laser point cloud and 2D image |
| 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 |
| CN113408688A (en) * | 2021-06-29 | 2021-09-17 | 哈尔滨工业大学 | Unknown environment-oriented multi-radioactive source online searching method |
| CN113408688B (en) * | 2021-06-29 | 2022-06-07 | 哈尔滨工业大学 | Unknown environment-oriented multi-radioactive source online searching method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109345571B (en) | 2020-10-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104361632B (en) | A kind of triangle gridding filling-up hole method based on Hermite RBFs | |
| CN104616349B (en) | Scattered point cloud data based on local surface changed factor simplifies processing method | |
| CN105740798B (en) | A kind of point cloud object scene recognition methods based on structural analysis | |
| CN110516388A (en) | Surface tessellation point cloud model ring cutting knife rail generating method based on reconciliation mapping | |
| CN100485662C (en) | Characteristic analytical method for product point clouds surface based on dynamic data access model | |
| CN109345571A (en) | A kind of point cloud registration method based on extension Gaussian image | |
| CN107093205A (en) | A kind of three dimensions building window detection method for reconstructing based on unmanned plane image | |
| CN101751695A (en) | Estimating method of main curvature and main direction of point cloud data | |
| CN110838115A (en) | Ancient cultural relic three-dimensional model change detection method by contour line extraction and four-dimensional surface fitting | |
| CN113327276B (en) | Mobile measurement-oriented general mass point cloud data registration method | |
| CN113034554B (en) | Whale optimized broken warrior body fragment registration method based on chaos reverse learning | |
| CN110827302A (en) | Point cloud target extraction method and device based on depth map convolutional network | |
| CN108230452A (en) | A kind of model filling-up hole method based on textures synthesis | |
| CN109299301A (en) | A kind of method for searching three-dimension model based on distribution of shapes and curvature | |
| CN117274339A (en) | Point cloud registration method based on improved ISS-3DSC characteristics combined with ICP | |
| Zhang et al. | Efficient and accurate Hausdorff distance computation based on diffusion search | |
| Cheng et al. | A novel point cloud simplification method using local conditional information | |
| Zhao | Point cloud denoising algorithm with geometric feature preserving | |
| Deng et al. | General multidimensional cloud model and its application on spatial clustering in Zhanjiang, Guangdong | |
| CN112991402B (en) | Wen Wudian cloud registration method and system based on improved differential evolution algorithm | |
| CN110207618A (en) | The surface line data extraction method of three-dimensional scanning measurement data | |
| Malinverni et al. | Automatic land use/land cover classification system with rules based both on objects attributes and landscape indicators | |
| CN110414379A (en) | In conjunction with the building extraction algorithm of elevation map Gabor textural characteristics and LiDAR point cloud feature | |
| Zhang et al. | Detection and filling of pseudo-hole in complex curved surface objects | |
| Cheng et al. | Segment-tree stereo matching algorithm based on improved matching costs |
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 |