CN109389625A - A kind of three-dimensional image registration method screening out error hiding based on multiple dimensioned description - Google Patents
A kind of three-dimensional image registration method screening out error hiding based on multiple dimensioned description Download PDFInfo
- Publication number
- CN109389625A CN109389625A CN201811169910.7A CN201811169910A CN109389625A CN 109389625 A CN109389625 A CN 109389625A CN 201811169910 A CN201811169910 A CN 201811169910A CN 109389625 A CN109389625 A CN 109389625A
- Authority
- CN
- China
- Prior art keywords
- point
- registration
- indicate
- value
- circular areas
- 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)
- Image Analysis (AREA)
Abstract
The invention discloses a kind of three-dimensional image registration method for screening out error hiding based on multiple dimensioned description, this method preferably describes the feature of corresponding key point, and tentatively obtain Corresponding matching point by constructing a kind of novel multiple dimensioned description;And traverse as feature has the point of similar multiple dimensioned description to converge therewith in subject to registration cloud, greatly improve an operational efficiency for cloud rough registration, reduce the calculation amount of computer, bring great convenience for point cloud registering;This method can obtain more accurate registration effect in the shorter time, and robustness is more preferable, suitable for there are noise, structure is complicated, to being registrated demanding field of precision measurement.
Description
Technical field
The invention belongs to image registration field, in particular to a kind of three-dimensional figure that error hiding is screened out based on multiple dimensioned description
As method for registering.
Background technique
Demand with electronics, the rapid development of computer technology and people to high quality of products, external machine vision
Enter high-speed development period in the 1990s, and is widely used in industrial control field.Machine vision is to realize that industry is automatic
Change and intelligentized necessary means, are equivalent to extension of the human vision on machine.Machine vision has increasingly automated, efficient
The advantages that rate, high-precision, will also play great effect in China's automatic field.It is regarded compared to traditional two-dimentional machine
Feel, the precision of three-dimensional machine vision is higher, detection effect is more preferable, is able to satisfy harsher production technology.But three-dimensional machine
The point cloud data amount that vision obtains is big and between points without apparent topological relation, therefore we are during registration
It is easy to appear and mismatches quasi- right, this will be such that our registration accuracy reduces, and influence measurement and detection effect to object.
In aerospace, the intelligence manufactures technique such as mechanical transport, bridge ship, adding for Complex Different Shape curved surface can be all encountered
Work and detection, curved surface complex bend is simultaneously folded in space arbitrarily product, and edge is complicated, and two-dimentional machine vision is difficult to meet its height
Precision, efficient production and processing requirement.But existing point cloud registration algorithm accuracy rate is not high, largely mismatches quasi- to depositing
In the precision and efficiency for affecting registration, detection and reconstruction to object, which also result in, to be seriously affected.
It is directed to this problem at present, most widely used solution is RANSAC algorithm, although RANSAC algorithm is big
It works well in most cases, but its runing time can be exponentially increased with the ratio η of abnormal point pair, therefore the calculation
Method is very sensitive to the selection of sample point in point cloud registering, and robustness is poor.During registration, the distribution of point cloud density is
Irregular and its structure also lacks regularity, so that the accuracy rate of 3-D image registration is more much lower than two dimensional image.
Summary of the invention
The present invention provides a kind of three-dimensional image registration method that error hiding is screened out based on multiple dimensioned description, this method benefits
It is screened with the space geometry relationship between registration point pair, and then can directly reject and mismatch a little, improve registration efficiency and match
Quasi- precision.
A kind of three-dimensional image registration method screening out error hiding based on multiple dimensioned description, comprising the following steps:
Step 1: obtaining the three-dimensional point cloud X and Y of the adjacent scene in visual angle, and more for each point construction in three-dimensional point cloud X and Y
Scale describes son (Qi, Ni);
Wherein, l indicates the neighborhood search radius label at the midpoint three-dimensional point cloud X or Y, and value is followed successively by 1,2 ..., and L, L are indicated
Neighborhood search radius sum, L >=2;TlIndicate normalized vector,λl1, λl2, λl3For Matrix ClIt carries out
Characteristic value after SVD decomposition, λl1≥λl2≥λl3, nl1, nl2, nl3Respectively λl1, λl2, λl3Corresponding feature vector,Sl={ pi|||pi-p0||≤rl, | Sl| indicate SlMaximum value, piIndicate three
Certain point p in dimension point cloud X or Y0Neighborhood point, rlIndicate point p0First of neighborhood search radius,rLFor setting
Maximum neighborhood search radius, nlRespectively indicate the point normal vector corresponding with the search radius of neighbourhood l of point in three-dimensional point cloud, value
For nl3;
|Sl| it is with p0For the centre of sphere, rlFor the maximum value in the distance of all the points in the ball of the centre of sphere to the centre of sphere;
Step 2: multiple dimensioned description of each point in traversal three-dimensional point cloud X and Y will have identical multiple dimensioned retouch in X and Y
The point for stating son is recorded as crucial point set A respectivelyi0And Bi0;
Step 3: according to principal curvatures signature search key point set Ai0And Bi0In initial matching point set Ai1And Bi1;
Initial matching point set Ai1In point ai1And Bi1In bi1It is matching double points;
From crucial point set Ai0Choose arbitrary point ai, traverse Bi0Midpoint is selected and point aiThe point b to matchj, following to meet
The point of condition constructs initial matching point set Ai1And Bi1:
Wherein, k1(ai) and k2(ai) it is illustrated respectively in point aiPlace, utilizes point aiNeighborhood point obtain principal curvatures minimum
With principal curvatures maximum, k1(bj) and k2(bj) it is illustrated respectively in point bjPlace, utilizes point bjNeighborhood point obtain principal curvatures it is minimum
Value and principal curvatures maximum, δ1For the first error threshold value of setting;
Step 4: using similar triangles threshold value to initial matching point set Ai1And Bi1In matching double points screened, obtain
Registration point is obtained to collection Ai2And Bi2;
Step 5: A is calculated based on SVDi2And Bi2Between the first spin matrix R1With the first translation matrix T1, utilize
One spin matrix R1With the first translation matrix T1, successively calculate Ai2In each point ai2Change point a3iCoordinate, while calculating Ai2
In each point in Bi2In match point b2iWith the Euclidean distance d between correspondent transform point3i;
Ai3=R1*Ai2+T1
Wherein, (a3ix,a3iy,a3iz) and (bi2x,bi2y,bi2z) it is respectively a3iAnd bi2Coordinate on x, y, z axis;
Step 6: calculatingThe initial value of k is that 1, N is current registration point to collection Ai2And Bi2Middle registration point
Logarithm, and judge dk+1-dk< ξ3It is whether true, if so, then with current registration point to collection Ai2And Bi2As three-dimensional point cloud X and Y
Registration data otherwise enable k=k+1, return step 3 chooses each neighborhood of a point point again, calculates principal curvatures.
Further, rotation angle threshold θ is utilizedmRegistration is screened again to doing, by current registration point to collection Ai2And Bi2
In, it is unsatisfactory for the corresponding registration point of following formula and is rejected to as abnormal point or noise spot:
Δθi< θm
Wherein, Δ θiIndicate current registration point to collection Ai2And Bi2Rotation angle difference between middle i-th pair registration point pair;
Current registration point is to collection Ai2And Bi2Middle i-th pair registration point is to (xi,yi) between rotation angle difference solution
Journey is as follows:
1) from current registration point to collection Ai2And Bi2In, non-coplanar two pairs of registrations pair are randomly selected, it is non-coplanar with four
Point constitutes a spherical surface, finds out affiliated spherical equation and centre of sphere O by analytic geometry and establishes space coordinates using point O as origin;
2) using the Z axis for the space coordinates built as the second spin matrix R2Rotary shaft, and R2=A2B2, A2And B2Point
The second spin matrix R is not indicated2Split-matrix, B2By current registration point to collection Ai2And Bi2Middle a pair of registration point is to (xk,yk)
Meet B2xk=ykWhen, carry out calculating acquisition, A2Rotary shaft be Z axis, meanwhile, according to the registration distance threshold ξ of setting, calculate
It is registrated angle threshold εi,
3) rotation byThe border circular areas surrounded, whenThe left margin of the border circular areas surrounded with
yiWhen tangent, byThe rotation angle of the border circular areas surrounded is θ1, at this time byThe circle surrounded
The center of circle in domain is point M, by point M and yiDistance rsA quadratic equation with one unknown is constructed, using equation of the point M in great circle O,
Two equation groups of simultaneous obtain the x and y coordinates value of point M, find out θ1;
Indicate point xiBy B2Rigid body translation is formed with B2xiFor the centre of sphere, εiFor the border circular areas of radius;
Because of B2It is the rigid body translation spin matrix when having noise or abnormal point,It indicates: xiThrough
Cross B2Rigid body translation, but because there are noise or abnormal point, can only find out a satisfactory range, and this range
It is exactly with B2xiFor the center of circle, εiFor the border circular areas of radius;
The great circle O is to give directions xiIn shown space coordinates, using point O as the circle where origin;
4) continue to rotateThe border circular areas surrounded, byThe right margin of the border circular areas surrounded
With yiWhen tangent,The rotation angle of the border circular areas surrounded is θ2, at this time byThe circle surrounded
The center of circle in domain is point N, by point N and yiDistance rsA quadratic equation with one unknown is constructed, equation of the N in great circle O, connection are utilized
Vertical two equation groups obtain the x and y coordinates value of point N, find out θ2, obtain angle difference Δ θi;
Wherein, rsIt isThe radius of the border circular areas surrounded, г indicate the B in space coordinates O2xiTo point O
Distance, rs=г * tan εi。
Further, the principal curvatures calculating process of each point is as follows in the three-dimensional point cloud:
Step 3.1: in three-dimensional point cloud arbitrary point and the neighborhood of a point point be fitted to a local secondary curved surface G (x,
y);
G (x, y)=ax2+bxy+cy2+dx+ey
Step 3.2: derivation being carried out to local quadratic surface G (x, y), it is basic to obtain first fundamental constant E, F, G and second
Constant L, M, N;
E=Gxx, F=Gx*Gy, G=Gyy;
L=nGxx, M=nGxy, N=nGyy;
Wherein, Gx、GyRespectively G (x, y) asks x, y and once leads, Gxx、GyyIt asks secondary to x, y for G (x, y) to lead, GxyFor G
(x, y) successively to x and y derivation, n is normal vector of the local secondary curved surface G (x, y) at selected arbitrary point;
Step 3.3: calculate Gaussian curvature K and mean curvature H:
Step 3.4: calculating principal curvatures minimum k1With principal curvatures maximum k2:
Further, described to utilize similar triangles threshold value to initial matching point set Ai1And Bi1In matching double points carry out
The process of screening is as follows:
Step 4.1: in Ai1In randomly select 3 not conllinear point ai,aj,ak, in Bi1Find corresponding 3 point bu,bv,
bw;
Step 4.2: enabling (ai,bu) and (aj,bv) it is known matching double points, it is carried out according to Similar Principle of Triangle following
It calculates:
Δaij=| | ai-aj||;Δaik=| | ai-ak||;Δajk=| | aj-ak||;
Δbuv=| | bu-bv||;Δbuw=| | bu-bw||;Δbwv=| | bw-bv||;
Step 4.3: by Similar Principle of Triangle, calculating two triangle similarity Ω;
Wherein, Δ aij,Δaik,Δajk;Δbuv,Δbuw,ΔbwvFor each side length of triangle, δ2It is the second of setting
Error threshold;
Step 4.4: if ai,aj,ak, bu,bv,bwMeet Ω < δ2, then it is assumed that akAnd bwBelong to a pair of of registration point pair.
Further, the L value is 4, and r4=0.004m.
Further, the δ2Value is any one in 0.02,0.05,0.1.
Further, the three-dimensional point cloud of the adjacent scene in visual angle is obtained using 3D profile scanning sensor.
Beneficial effect
The present invention provides a kind of three-dimensional image registration methods that error hiding is screened out based on multiple dimensioned description, pass through construction
A kind of novel multiple dimensioned description, preferably describes the feature of corresponding key point, and tentatively obtain Corresponding matching point;And with this
Be characterized in subject to registration cloud of traversal has the point of similar multiple dimensioned description to converge therewith, greatly improves a cloud rough registration
Operational efficiency reduces the calculation amount of computer, brings great convenience for point cloud registering;This method can be when shorter
Between in obtain more accurate registration effect, and robustness is more preferable, suitable for there are noise, structure is complicated, to high with alignment request
Field of precision measurement.
After essence registration, the point cloud pair being registrated can inevitably exist compared to traditional method for registering and mismatch standard
Pair and noise, accuracy rate and precision of the objective in measurement and detection certainly will be will affect.Application now is at most
It is ICP iterative method to carry out optimal registration, but largely mismatches quasi- pair of presence, largely will increase algorithm and fall into part
A possibility that optimal or algorithm can not restrain seriously affects the reliability of registration.Therefore, the present invention is based on rotation angle by introducing
Degree threshold value mismatches quasi- right further to screen out, so that the quantity for mismatching quasi- pair and noise spot largely reduces, ensure that registration knot
The accuracy of fruit makes it have good robustness;On the other hand, quasi- pair and noise spot are mismatched due to having screened out a part,
The input quantity of a cloud is decreased, the calculation amount in registration process is reduced, improves the efficiency of registration.
Detailed description of the invention
Fig. 1 is the flow diagram of the method for the invention;
Fig. 2 is that of the present invention utilize rotates angle progress registration process schematic diagram, wherein figure (a) is space coordinate
It is schematic diagram, figure (b) is that threshold xi is converted into angle εiSchematic diagram, figure (c) be xiAnd xkRespectively in R and R2Signal under transformation
Figure, figure (d) are B2xiRotary course schematic diagram, figure (e) beThe border circular areas surrounded is in space coordinates
Relation schematic diagram.
Specific embodiment
Below in conjunction with drawings and examples, the present invention is described further.
As shown in Figure 1, a kind of three-dimensional image registration method that error hiding is screened out based on multiple dimensioned description, including following step
It is rapid:
Step 1: obtaining the three-dimensional point cloud X and Y of the adjacent scene in visual angle, and more for each point construction in three-dimensional point cloud X and Y
Scale describes son (Qi, Ni);
The three-dimensional point cloud of the adjacent scene in visual angle is obtained in this example using 3D profile scanning sensor;
Wherein, l indicates the neighborhood search radius label at the midpoint three-dimensional point cloud X or Y, and value is followed successively by 1,2 ..., and L, L are indicated
Neighborhood search radius sum, L >=2;TlIndicate normalized vector,λl1, λl2, λl3For Matrix ClIt carries out
Characteristic value after SVD decomposition, λl1≥λl2≥λl3, nl1, nl2, nl3Respectively λl1, λl2, λl3Corresponding feature vector,Sl={ pi|||pi-p0||≤rl, | Sl| indicate SlMaximum value, piIndicate three
Certain point p in dimension point cloud X or Y0Neighborhood point, rlIndicate point p0First of neighborhood search radius,rLFor setting
Maximum neighborhood search radius, nlRespectively indicate the point normal vector corresponding with the search radius of neighbourhood l of point in three-dimensional point cloud, value
For nl3;
In this example, the L value is 4, and r4=0.004m.
|Sl| it is with p0For the centre of sphere, rlFor the maximum value in the distance of all the points in the ball of the centre of sphere to the centre of sphere;
Step 2: multiple dimensioned description of each point in traversal three-dimensional point cloud X and Y will have identical multiple dimensioned retouch in X and Y
The point for stating son is recorded as crucial point set A respectivelyi0And Bi0;
Step 3: according to principal curvatures signature search key point set Ai0And Bi0In initial matching point set Ai1And Bi1;
Initial matching point set Ai1In point ai1And Bi1In bi1It is matching double points;
From crucial point set Ai0Choose arbitrary point ai, traverse Bi0Midpoint is selected and point aiThe point b to matchj, following to meet
The point of condition constructs initial matching point set Ai1And Bi1:
Wherein, k1(ai) and k2(ai) it is illustrated respectively in point aiPlace, utilizes point aiNeighborhood point obtain principal curvatures minimum
With principal curvatures maximum, k1(bj) and k2(bj) it is illustrated respectively in point bjPlace, utilizes point bjNeighborhood point obtain principal curvatures it is minimum
Value and principal curvatures maximum, δ1For the first error threshold value of setting;
First error threshold value δ1Size one best value of effect is acquired by experiment simulation;
The principal curvatures calculating process of each point is as follows in the three-dimensional point cloud:
Step 3.1: in three-dimensional point cloud arbitrary point and the neighborhood of a point point be fitted to a local secondary curved surface G (x,
y);
G (x, y)=ax2+bxy+cy2+dx+ey
Step 3.2: derivation being carried out to local quadratic surface G (x, y), it is basic to obtain first fundamental constant E, F, G and second
Constant L, M, N;
E=Gxx, F=Gx*Gy, G=Gyy;
L=nGxx, M=nGxy, N=nGyy;
Wherein, Gx、GyRespectively G (x, y) asks x, y and once leads, Gxx、GyyIt asks secondary to x, y for G (x, y) to lead, GxyFor G
(x, y) successively to x and y derivation, n is normal vector of the local secondary curved surface G (x, y) at selected arbitrary point;
Step 3.3: calculate Gaussian curvature K and mean curvature H:
Step 3.4: calculating principal curvatures minimum k1With principal curvatures maximum k2:
Step 4: using similar triangles threshold value to initial matching point set Ai1And Bi1In matching double points screened, obtain
Registration point is obtained to collection Ai2And Bi2;
It is described to utilize similar triangles threshold value to initial matching point set Ai1And Bi1In the process screened of matching double points
It is as follows:
Step 4.1: in Ai1In randomly select 3 not conllinear point ai,aj,ak, in Bi1Find corresponding 3 point bu,bv,
bw;
Step 4.2: enabling (ai,au) and (aj,av) it is known matching double points, it is carried out according to Similar Principle of Triangle following
It calculates:
Δaij=| | ai-aj||;Δaik=| | ai-ak||;Δajk=| | aj-ak||;
Δbuv=| | bu-bv||;Δbuw=| | bu-bw||;Δbwv=| | bw-bv||;
Step 4.3: by Similar Principle of Triangle, calculating two triangle similarity Ω;
Wherein, Δ aij,Δaik,Δajk;Δbuv,Δbuw,ΔbwvFor each side length of triangle, δ2It is the second of setting
Error threshold;
In this example, the δ2Value is any one in 0.02,0.05,0.1.
Step 4.4: if ai,aj,ak, bu,bv,bwMeet Ω < δ2, then it is assumed that akAnd bwBelong to a pair of of registration point pair.
Step 5: A is calculated based on SVDi2And Bi2Between the first spin matrix R1With the first translation matrix T1, utilize
One spin matrix R1With the first translation matrix T1, successively calculate Ai2In each point ai2Change point a3iCoordinate, while calculating Ai2
In each point in Bi2In match point b2iWith the Euclidean distance d between correspondent transform point3i;
Ai3=R1*Ai2+T1
Wherein, (a3ix,a3iy,a3iz) and (bi2x,bi2y,bi2z) it is respectively a3iAnd bi2Coordinate on x, y, z axis;
Step 6: calculatingThe initial value of k is that 1, N is current registration point to collection Ai2And Bi2Middle registration point
Logarithm, and judge dk+1-dk< ξ3It is whether true, if so, then with current registration point to collection Ai2And Bi2As three-dimensional point cloud X and Y
Registration data otherwise enable k=k+1, return step 3 chooses each neighborhood of a point point again, calculates principal curvatures.
In this example, ξ3Value is 0.01.
Utilize rotation angle threshold θmRegistration is screened again to doing, by current registration point to collection Ai2And Bi2In, it is unsatisfactory for
The corresponding registration point of following formula is rejected to as abnormal point or noise spot:
Δθi< θm
Wherein, Δ θiIndicate current registration point to collection Ai2And Bi2Rotation angle difference between middle i-th pair registration point pair;
Current registration point is to collection Ai2And Bi2Middle i-th pair registration point is to (xi,yi) between rotation angle difference solution
Journey is as follows:
1) from current registration point to collection Ai2And Bi2In, non-coplanar two pairs of registrations pair are randomly selected, it is non-coplanar with four
Point constitutes a spherical surface, finds out affiliated spherical equation and centre of sphere O by analytic geometry and establishes space coordinates using point O as origin,
As shown in Fig. 2 (a), and from Fig. 2 (a), it can be seen that with yiThe public area intersected is had with sphere O for the sphere of the centre of sphere
Domain, this region are considered yiMeet the region of registration condition;
As long as xiBy R1、T1Transformed point belong to yiFor the centre of sphere, ξ is in the sphere of radius, then it is assumed that xiAnd yi
It is a pair of of registration pair, i.e., the satisfaction to threshold xi is converted into angle εiSolution, as shown in Fig. 2 (b);
Consider be based on a screening conditions (rotation transformation), to meet rotation and translation transformation when angle threshold value,
The threshold value of rotation angle when must first meet based on rotation transformation, most homogeneous problem is converted are as follows: ∠ (R1xi,
yi)≤εi, wherein εiFor the angle threshold for judging whether it is registration point pair;
By Fig. 2 (b) it follows that
Based on above-mentioned analysis, it does such as down conversion:
1. R1It is decomposed into two matrix As, B, i.e. R1=A*B;(B is set to meet condition: ∠ (Bxi,yi)≤εi;Make A simultaneously
For BxiRotary shaft, most homogeneous problem is converted into the rotation angle asked and A is made to meet condition as far as possible;
2. when no noise and abnormal point, xiAnd yiIt can perfectly be stitched together very much.Define one
R in the ideal situation2=A2B2.By a pair of registration to xk,yk, can be in the hope of B2(B2xk=yk);Using the point of i ≠ k, with
B2xkFor rotary shaft (i.e. A in Fig. 2 (c)2Rotary shaft), find out it and make B2xiWith yiThe rotation angle θ that perfection is mixed2(because of point
B2xiAnd yiCoordinate in space coordinates is known, A2Rotary shaft can be regarded as the Z axis in space coordinates, by
The relationship of analytic geometry can directly find out rotation angle θ2, that is, acquire A2:
A2=exp (θ2[B2xk]X)
Wherein, [X]XIndicate the vector product of X, exp () indicates index mapping relationship.
3. but noise and abnormal point are inevitable under normal conditions, therefore utilize R2To estimate xkAfter rotation
One range limit;Make ∠ (Bxk,yk)≤εi, it is known that BxkIt is in Fig. 2 (c) in a hypographous border circular areas of tool.Together
Reason, for the point of i ≠ k, can obtain BxiAlso in figure (c) in a hypographous border circular areas of tool.
It is denoted as:
2) using the Z axis for the space coordinates built as the second spin matrix R2Rotary shaft, and R2=A2B2, A2And B2Point
The second spin matrix R is not indicated2Split-matrix, B2By current registration point to collection Ai2And Bi2Middle a pair of registration point is to (xk,yk)
Meet B2xk=ykWhen, carry out calculating acquisition, A2Rotary shaft be Z axis, meanwhile, according to the registration distance threshold ξ of setting, calculate
It is registrated angle threshold εi,
3) rotation byThe border circular areas surrounded, whenThe left margin of the border circular areas surrounded with
yiWhen tangent (shown in such as Fig. 2 (d)), byThe rotation angle of the border circular areas surrounded is θ1, at this time by
The center of circle of the border circular areas surrounded is point M, by point M and yiDistance rsA quadratic equation with one unknown is constructed, is existed using point M
Equation in great circle O, two equation groups of simultaneous obtain the x and y coordinates value of point M, find out θ1;
Indicate point xiBy B2Rigid body translation is formed with B2xiFor the centre of sphere, εiFor the border circular areas of radius, such as scheme
Shown in 2 (e);
Because of B2It is the rigid body translation spin matrix when having noise or abnormal point,It indicates: xiThrough
Cross B2Rigid body translation, but because there are noise or abnormal point, can only find out a satisfactory range, and this range
It is exactly with B2xiFor the centre of sphere, εiFor the border circular areas of radius;
The great circle O is to give directions xiIn the space coordinates, using point O as the circle where origin;
4) continue to rotateThe border circular areas surrounded, byThe right margin of the border circular areas surrounded
With yiWhen tangent,The rotation angle of the border circular areas surrounded is θ2, at this time byThe circle surrounded
The center of circle in region is point N, by point N and yiDistance rsA quadratic equation with one unknown is constructed, using equation of the N in great circle O,
Two equation groups of simultaneous obtain the x and y coordinates value of point N, find out θ2, obtain angle difference Δ θi;
Wherein, rsIt isThe radius of the border circular areas surrounded, г indicate the B in space coordinates O2xiTo point O
Distance, rs=г * tan εi。
Specific embodiment described herein is only an example for the spirit of the invention.The neck of technology belonging to the present invention
The technical staff in domain can make various modifications or additions to the described embodiments or replace by a similar method
In generation, however, it does not deviate from the spirit of the invention or beyond the scope of the appended claims.
Claims (7)
1. a kind of three-dimensional image registration method for screening out error hiding based on multiple dimensioned description, which is characterized in that including following step
It is rapid:
Step 1: obtaining the three-dimensional point cloud X and Y of the adjacent scene in visual angle, and multiple dimensioned for each point construction in three-dimensional point cloud X and Y
Son (Q is describedi, Ni);
Ni={ nl};
Wherein, 1 indicate that the neighborhood search radius at the midpoint three-dimensional point cloud X or Y marks, value is followed successively by 1,2 ..., and L, L indicate neighborhood
Search radius sum, L >=2;TlIndicate normalized vector,λl1, λl2, λl3For Matrix ClCarry out SVD points
Characteristic value after solution, λl1≥λl2≥λl3, nl1, nl2, nl3Respectively λl1, λl2, λl3Corresponding feature vector,Sl={ pi|||pi-p0||≤rl, | Sl| indicate SlMaximum value, piIndicate three
Certain point p in dimension point cloud X or Y0Neighborhood point, rlIndicate point p0The 1st neighborhood search radius,rLFor setting
Maximum neighborhood search radius, nlRespectively indicate the point normal vector corresponding with the search radius of neighbourhood 1 of point in three-dimensional point cloud, value
For nl3;
Step 2: multiple dimensioned description of each point in traversal three-dimensional point cloud X and Y will have identical multiple dimensioned description in X and Y
Point be recorded as crucial point set A respectivelyi0And Bi0;
Step 3: according to principal curvatures signature search key point set Ai0And Bi0In initial matching point set Ai1And Bi1;
From crucial point set Ai0Choose arbitrary point ai, traverse Bi0Midpoint is selected and point aiThe point b to matchj, to meet the following conditions
Point construct initial matching point set Ai1And Bi1:
Wherein, k1(ai) and k2(ai) it is illustrated respectively in point aiPlace, utilizes point aiNeighborhood point obtain principal curvatures minimum and master
Curvature maximum, k1(bj) and k2(bj) it is illustrated respectively in point bjPlace, utilizes point bjNeighborhood point obtain principal curvatures minimum and
Principal curvatures maximum, δ1For the first error threshold value of setting;
Step 4: using similar triangles threshold value to initial matching point set Ai1And Bi1In matching double points screened, be registrated
Point is to collection Ai2And Bi2;
Step 5: A is calculated based on SVDi2And Bi2Between the first spin matrix R1With the first translation matrix T1, utilize the first rotation
Torque battle array R1With the first translation matrix T1, successively calculate Ai2In each point ai2Change point a3iCoordinate, while calculating Ai2In it is every
A point is in Bi2In match point b2iWith the Euclidean distance d between correspondent transform point3i;
Ai3=R1*Ai2+T1
Wherein, (a3ix, a3iy, a3iz) and (bi2x, bi2y, bi2z) it is respectively a3iAnd bi2Coordinate on x, y, z axis;
Step 6: calculatingThe initial value of k is that 1, N is current registration point to collection Ai2And Bi2Middle registration point pair
Number, and judge dk+1-dk< ξ3It is whether true, ξ3Indicate the error threshold greater than zero, if so, then with current registration point to collection Ai2
And Bi2As the registration data of three-dimensional point cloud X and Y, otherwise, k=k+1 being enabled, return step 3 chooses each neighborhood of a point point again,
Calculate principal curvatures.
2. the method according to claim 1, wherein utilizing rotation angle threshold θmRegistration is screened again to doing,
By current registration point to collection Ai2And Bi2In, the corresponding registration point of following formula is unsatisfactory for as abnormal point or noise spot progress
It rejects:
Δθi< θm
Wherein, Δ θiIndicate current registration point to collection Ai2And Bi2Rotation angle difference between middle i-th pair registration point pair;
Current registration point is to collection Ai2And Bi2Middle i-th pair registration point is to (xi, yi) between rotation angle difference solution procedure such as
Under:
1) from current registration point to collection Ai2And Bi2In, two pairs of non-coplanar registrations pair are randomly selected, with four non-coplanar point structures
At a spherical surface, affiliated spherical equation and the centre of sphere 0 are found out by analytic geometry and establish space coordinates with point 0 for origin;
2) using the Z axis for the space coordinates built as the second spin matrix R2Rotary shaft, and R2=A2B2, A2And B2Table respectively
Show the second spin matrix R2Split-matrix, B2By current registration point to collection Ai2And Bi2Middle a pair of registration point is to (xk, yk) full
Sufficient B2xk=ykWhen, carry out calculating acquisition, A2Rotary shaft be Z axis, meanwhile, according to the registration distance threshold ξ of setting, calculating is matched
Quasi- angle threshold εi,
3) rotation byThe border circular areas surrounded, whenThe left margin and y of the border circular areas surroundediPhase
When cutting, byThe rotation angle of the border circular areas surrounded is θ1, at this time byThe border circular areas surrounded
The center of circle is point M, by point M and yiDistance rsA quadratic equation with one unknown is constructed, equation of the point M in great circle 0, simultaneous are utilized
Two equation groups obtain the x and y coordinates value of point M, find out θ1;
Indicate point xiBy B2Rigid body translation is formed with B2xiFor the centre of sphere, εiFor the border circular areas of radius;
The great circle 0 is to give directions xiIn the space coordinates, to put 0 circle where origin;
4) continue to rotateThe border circular areas surrounded, byThe right margin and y of the border circular areas surroundedi
When tangent,The rotation angle of the border circular areas surrounded is θ2, at this time byThe border circular areas surrounded
The center of circle is point N, by point N and yiDistance rsA quadratic equation with one unknown is constructed, equation of the N in great circle 0, simultaneous two are utilized
A equation group obtains the x and y coordinates value of point N, finds out θ2, obtain angle difference Δ θi;
Wherein, rsIt isThe radius of the border circular areas surrounded, r indicate the B in space coordinates 02xiTo point 0 away from
From rs=r*tan εi。
3. the method according to claim 1, wherein in the three-dimensional point cloud each point principal curvatures calculating process
It is as follows:
Step 3.1: in three-dimensional point cloud arbitrary point and the neighborhood of a point point be fitted to a local secondary curved surface G (x, y);
G (x, y)=ax2+bxy+cy2+dx+ey
Step 3.2: derivation being carried out to local quadratic surface G (x, y), obtains first fundamental constant E, F, G and the second fundamental constant
L,M,N;
E=Gxx, F=Gx*Gy, G=Gyy;
L=nGxx, M=nGxy, N=nGyy;
Wherein, Gx、GyRespectively G (x, y) asks x, y and once leads, Gxx、GyyIt asks secondary to x, y for G (x, y) to lead, GxyFor G (x, y)
Successively to x and y derivation, n is normal vector of the local secondary curved surface G (x, y) at selected arbitrary point;
Step 3.3: calculate Gaussian curvature K and mean curvature H:
Step 3.4: calculating principal curvatures minimum k1With principal curvatures maximum k2:
4. the method according to claim 1, wherein described utilize similar triangles threshold value to initial matching point set
Ai1And Bi1In the process screened of matching double points it is as follows:
Step 4.1: in Ai1In randomly select 3 not conllinear point ai, aj, ak, in Bi1Find corresponding 3 point bu, bv, bw;
Step 4.2: enabling (ai, bu) and (aj, bv) it is known matching double points, following calculate is carried out according to Similar Principle of Triangle:
Δaij=| | ai-aj||;Δaik=| | ai-ak||;Δajk=| | aj-ak||;
Δbuv=| | bu-bv||;Δbuw=| | bu-bw||;Δbwv=| | bw-bv||;
Step 4.3: by Similar Principle of Triangle, calculating two triangle similarity Ω;
Wherein, Δ aij, Δ aik, Δ ajk;Δbuv, Δ buw, Δ bwvFor each side length of triangle, δ2For the second error of setting
Threshold value;
Step 4.4: if ai, aj, ak, bu, bv, bwMeet Ω < δ2, then it is assumed that akAnd bwBelong to a pair of of registration point pair.
5. method according to claim 1-4, which is characterized in that the L value is 4, and r4=0.004m.
6. according to the method described in claim 4, it is characterized in that, the δ2Value is any one in 0.02,0.05,0.1
It is a.
7. the method according to claim 1, wherein obtaining the adjacent scene in visual angle using 3D profile scanning sensor
Three-dimensional point cloud.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811169910.7A CN109389625B (en) | 2018-10-08 | 2018-10-08 | Three-dimensional image registration method based on multi-scale descriptor screening and mismatching |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811169910.7A CN109389625B (en) | 2018-10-08 | 2018-10-08 | Three-dimensional image registration method based on multi-scale descriptor screening and mismatching |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109389625A true CN109389625A (en) | 2019-02-26 |
CN109389625B CN109389625B (en) | 2021-09-14 |
Family
ID=65419294
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811169910.7A Active CN109389625B (en) | 2018-10-08 | 2018-10-08 | Three-dimensional image registration method based on multi-scale descriptor screening and mismatching |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109389625B (en) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110276790A (en) * | 2019-06-28 | 2019-09-24 | 易思维(杭州)科技有限公司 | Point cloud registration method based on shape constraining |
CN110312070A (en) * | 2019-04-23 | 2019-10-08 | 维沃移动通信有限公司 | A kind of image processing method and terminal |
CN110335297A (en) * | 2019-06-21 | 2019-10-15 | 华中科技大学 | A kind of point cloud registration method based on feature extraction |
CN111611996A (en) * | 2020-04-22 | 2020-09-01 | 青岛联合创智科技有限公司 | Computing method of point cloud characteristic point descriptor |
CN117495932A (en) * | 2023-12-25 | 2024-02-02 | 国网山东省电力公司滨州供电公司 | Power equipment heterologous point cloud registration method and system |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104143210A (en) * | 2014-07-31 | 2014-11-12 | 哈尔滨工程大学 | Multi-scale normal feature point cloud registering method |
CN105118059A (en) * | 2015-08-19 | 2015-12-02 | 哈尔滨工程大学 | Multi-scale coordinate axis angle feature point cloud fast registration method |
US20160321838A1 (en) * | 2015-04-29 | 2016-11-03 | Stmicroelectronics S.R.L. | System for processing a three-dimensional (3d) image and related methods using an icp algorithm |
US9547838B2 (en) * | 2013-11-06 | 2017-01-17 | Oracle International Corporation | Automated generation of a three-dimensional space representation and planogram verification |
CN108022262A (en) * | 2017-11-16 | 2018-05-11 | 天津大学 | A kind of point cloud registration method based on neighborhood of a point center of gravity vector characteristics |
CN108537805A (en) * | 2018-04-16 | 2018-09-14 | 中北大学 | A kind of target identification method of feature based geometry income |
-
2018
- 2018-10-08 CN CN201811169910.7A patent/CN109389625B/en active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US9547838B2 (en) * | 2013-11-06 | 2017-01-17 | Oracle International Corporation | Automated generation of a three-dimensional space representation and planogram verification |
CN104143210A (en) * | 2014-07-31 | 2014-11-12 | 哈尔滨工程大学 | Multi-scale normal feature point cloud registering method |
US20160321838A1 (en) * | 2015-04-29 | 2016-11-03 | Stmicroelectronics S.R.L. | System for processing a three-dimensional (3d) image and related methods using an icp algorithm |
CN105118059A (en) * | 2015-08-19 | 2015-12-02 | 哈尔滨工程大学 | Multi-scale coordinate axis angle feature point cloud fast registration method |
CN108022262A (en) * | 2017-11-16 | 2018-05-11 | 天津大学 | A kind of point cloud registration method based on neighborhood of a point center of gravity vector characteristics |
CN108537805A (en) * | 2018-04-16 | 2018-09-14 | 中北大学 | A kind of target identification method of feature based geometry income |
Non-Patent Citations (2)
Title |
---|
JUN LU ET.AL: "Fast point cloud registration algorithm using multiscale angle features", 《JOURNAL OF ELECTRONIC IMAGING》 * |
张曼 等: "一种基于尺度空间的三维点云数据配准算法", 《***仿真学报》 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110312070A (en) * | 2019-04-23 | 2019-10-08 | 维沃移动通信有限公司 | A kind of image processing method and terminal |
CN110335297A (en) * | 2019-06-21 | 2019-10-15 | 华中科技大学 | A kind of point cloud registration method based on feature extraction |
CN110335297B (en) * | 2019-06-21 | 2021-10-08 | 华中科技大学 | Point cloud registration method based on feature extraction |
CN110276790A (en) * | 2019-06-28 | 2019-09-24 | 易思维(杭州)科技有限公司 | Point cloud registration method based on shape constraining |
CN111611996A (en) * | 2020-04-22 | 2020-09-01 | 青岛联合创智科技有限公司 | Computing method of point cloud characteristic point descriptor |
CN111611996B (en) * | 2020-04-22 | 2023-06-20 | 青岛联合创智科技有限公司 | Calculation method of point cloud characteristic point descriptors |
CN117495932A (en) * | 2023-12-25 | 2024-02-02 | 国网山东省电力公司滨州供电公司 | Power equipment heterologous point cloud registration method and system |
CN117495932B (en) * | 2023-12-25 | 2024-04-16 | 国网山东省电力公司滨州供电公司 | Power equipment heterologous point cloud registration method and system |
Also Published As
Publication number | Publication date |
---|---|
CN109389625B (en) | 2021-09-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109389625A (en) | A kind of three-dimensional image registration method screening out error hiding based on multiple dimensioned description | |
Yang et al. | Automated registration of dense terrestrial laser-scanning point clouds using curves | |
CN103400388B (en) | Method for eliminating Brisk key point error matching point pair by using RANSAC | |
CN109767463A (en) | A kind of three-dimensional point cloud autoegistration method | |
CN104809738B (en) | A kind of air bag overall size detection method based on binocular vision | |
CN104930985B (en) | Binocular vision 3 D topography measurement method based on space-time restriction | |
CN104167003A (en) | Method for fast registering remote-sensing image | |
CN110473239A (en) | A kind of high-precision point cloud registration method of 3 D laser scanning | |
CN104121902B (en) | Implementation method of indoor robot visual odometer based on Xtion camera | |
CN110136178B (en) | Three-dimensional laser point cloud registration method and device based on endpoint fitting | |
CN104786098B (en) | Geometric error six-position recognition method of multi-axis numerical control machine tool rotary table | |
CN110211129B (en) | Low-coverage point cloud registration algorithm based on region segmentation | |
CN109118528A (en) | Singular value decomposition image matching algorithm based on area dividing | |
CN111223133A (en) | Registration method of heterogeneous images | |
CN101315698A (en) | Characteristic matching method based on straight line characteristic image registration | |
CN108986149A (en) | A kind of point cloud Precision Registration based on adaptive threshold | |
CN108830888A (en) | Thick matching process based on improved multiple dimensioned covariance matrix Feature Descriptor | |
CN105787451A (en) | Fingerprint matching method based on multi-judgment point mode | |
CN107633507A (en) | LCD defect inspection methods based on contour detecting and characteristic matching | |
CN114648445B (en) | Multi-view high-resolution point cloud splicing method based on feature point extraction and fine registration optimization | |
CN116721144A (en) | Cone hole size measurement method based on point cloud slicing | |
CN114627250B (en) | Human body standing posture three-dimensional reconstruction and measurement method based on Kinect | |
CN104282001A (en) | Method for enhancing image feature two-value descriptor performance | |
CN106251350B (en) | A kind of RCS picture based on infinitesimal inertia and objective contour method for registering | |
CN116402792A (en) | Space hole site butt joint method based on three-dimensional point cloud |
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 |