CN116091562A - Building point cloud automatic registration method based on two-dimensional projection line segments - Google Patents

Building point cloud automatic registration method based on two-dimensional projection line segments Download PDF

Info

Publication number
CN116091562A
CN116091562A CN202211256855.1A CN202211256855A CN116091562A CN 116091562 A CN116091562 A CN 116091562A CN 202211256855 A CN202211256855 A CN 202211256855A CN 116091562 A CN116091562 A CN 116091562A
Authority
CN
China
Prior art keywords
point
points
intersection
homonymous
qualified
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.)
Pending
Application number
CN202211256855.1A
Other languages
Chinese (zh)
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.)
Nanyang Normal University
Original Assignee
Nanyang Normal University
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 Nanyang Normal University filed Critical Nanyang Normal University
Priority to CN202211256855.1A priority Critical patent/CN116091562A/en
Publication of CN116091562A publication Critical patent/CN116091562A/en
Pending legal-status Critical Current

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
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G06T7/344—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/06—Topological mapping of higher dimensional structures onto lower dimensional surfaces
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • 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)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention provides a method for automatically registering building point clouds based on two-dimensional projection line segments, which can realize the rapid registration of large-scale building point clouds and comprises the following steps: cloud of target points P t And point cloud to be registered P s Respectively preprocessing to obtain two-dimensional point cloud
Figure DDA0003889903400000011
And
Figure DDA0003889903400000012
extracting two-dimensional point cloud line segments to form a line segment set, and arranging the line segments in descending order according to the length; acquiring qualified intersection points, attached virtual points and included angles in a line segment set; comparing the angles to determine qualified intersection points and virtual points attached to the qualified intersection points and homonymous point combinations of the qualified intersection points and the homonymous point combinations, and solvingConversion parameters corresponding to each homonymous point combination; obtain two types of counter C combined by homonymous points 1 And C 2 The sum of the cumulative values of the number of all the qualified interior points and the homonymous intersection point pair corresponding to each qualified interior point; from C 1 And C 2 And finding out the homonymous intersection point pair corresponding to the counter meeting the maximum accumulated value of the number of the inner points or exceeding a certain threshold value requirement as the intersection point pair meeting the registration condition for registration.

Description

Building point cloud automatic registration method based on two-dimensional projection line segments
Technical Field
The invention belongs to the technical field of three-dimensional laser scanning, and particularly relates to a building point cloud automatic registration method based on two-dimensional projection line segments.
Background
Although the ground laser scanner can acquire fine three-dimensional information of the surface of an object, the point cloud data acquired by a single view is insufficient to cover the whole scene, so that the point cloud data acquired by a plurality of views are necessarily registered under the same coordinate system, and the process comprises two steps of coarse registration and fine registration. Fine registration typically employs an Iterative Closest Point (ICP) algorithm or variants thereof. There are a number of point cloud coarse registration methods that calculate transformation parameters by extracting different geometric features (points, lines, planes and specific objects). Generally, the point-based method is more universal, because the method is suitable for various scenes, key points are firstly extracted from two point clouds to be matched, local shape descriptors of the key points are calculated, and the similarity of the descriptors is compared to establish the corresponding relation between the key points of different point clouds, and although the method can obtain a better registration effect, the calculation complexity is higher. The line or plane-based method has stronger robustness to noise and point density changes, and for a scene with an artificial object, the extraction of line/plane features for point cloud registration is a good choice because of a large number of line and plane features in the scene. However, most of the existing methods based on the line or the plane are extracting three-dimensional structure line/three-dimensional plane, which means that these methods process the point cloud data in the three-dimensional space, and extracting the three-dimensional structure line/three-dimensional plane is quite time-consuming, so that the accuracy and the speed of three-dimensional point cloud registration are limited, and the application of the method is seriously affected.
Disclosure of Invention
The invention aims to solve the problems of low efficiency, poor precision, easy sinking into local optimum, non-convergence and the like in large-scale building point cloud registration, and provides a method for automatically registering building point clouds based on two-dimensional projection line segments, which realizes rapid registration of the large-scale building point clouds.
In order to achieve the above object, the present invention adopts the following scheme:
the invention provides a method for automatically registering building point clouds based on two-dimensional projection line segments, which is characterized by comprising the following steps: step 1, setting the target point cloud as P t The point cloud to be registered is P s For P t and Ps All are carried out in sequence: the pretreatment of spatial uniform downsampling, filtering the main horizontal plane in the Z direction and projecting to the X-Y plane to obtain a two-dimensional point cloud
Figure BDA0003889903380000021
and />
Figure BDA0003889903380000022
Step 2, setting the farthest distance from the point to the straight line, the radius of neighbor searching during clustering, the minimum number of points contained in the qualified line segment and the shortest length of the qualified line segment as extraction requirements, and respectively taking the points and the radius of the neighbor searching, the minimum number of points contained in the qualified line segment and the shortest length of the qualified line segment as extraction requirements
Figure BDA0003889903380000023
and />
Figure BDA0003889903380000024
Extracting two-dimensional point cloud line segments meeting the extraction requirements to form a two-dimensional point cloud line segment set L t and Ls For L t and Ls The point cloud line segments in (a) are arranged in descending order according to the length;
step 3, acquiring qualified intersection points, attached virtual points and included angles of the two-dimensional point cloud line segment sets: for L t Firstly, calculating the longest point cloud line segment
Figure BDA0003889903380000025
Straight line and L t Any point cloud segment left->
Figure BDA0003889903380000026
Angle formed by the straight line>
Figure BDA0003889903380000027
If it is
Figure BDA0003889903380000028
When the intersection point is qualified, the intersection point of the two straight lines and the attached two virtual points are calculated and compared: let the intersection point of two intersecting straight lines be p inters (c 1 ,c 2 ) Assuming that a virtual point on a straight line having a length r from the intersection point is p (x, y), the following relationship exists: />
Figure BDA0003889903380000029
Substituting the first polynomial in equation (7) into the second polynomial, the value of x can be found as:
Figure BDA00038899033800000210
wherein ,
Figure BDA0003889903380000031
each straight line is calculated to obtain two virtual points, and then the two virtual points are calculated to p for each straight line ,ent Distance d of (2) 1 and d2 Distance p on each straight line cent The nearest virtual point is the virtual point attached to the qualified intersection point, p cent The mass center of the two intersecting point cloud line segments;
once the two line segments are aligned, they are then processedThe mark avoids repeated comparison, and the like, any two line segments are compared to obtain L t All qualified intersection points in the list form an intersection point set I t Each qualified intersection point is attached with two virtual points, and an included angle set theta formed by the included angle of each qualified intersection point and the attached two virtual points is calculated t The method comprises the steps of carrying out a first treatment on the surface of the By using the L t The same method obtains L s Intersection point set I of all qualified intersection points in the list s Two virtual points attached to each qualified intersection point and a corresponding included angle set theta s ；
Step 4. Comparing theta t Any angle of
Figure BDA0003889903380000032
And theta s Any angle->
Figure BDA0003889903380000033
Difference (I) of->
Figure BDA0003889903380000034
1≤i≤n 1 ，n 1 Is I t The number of subsets of (1) j n 2 ，n 2 Is P s Number of pass-through intersections if |Δθ|<1 DEG, the corresponding pair of intersection points and the virtual points attached to the intersection points are considered to be the same name points, and the corresponding same name point relationship is the first +.>
Figure BDA0003889903380000035
And->
Figure BDA0003889903380000036
And->
Figure BDA0003889903380000037
And (3) with
Figure BDA0003889903380000038
Or second +.>
Figure BDA0003889903380000039
And->
Figure BDA00038899033800000310
And->
Figure BDA00038899033800000311
And->
Figure BDA00038899033800000312
Two combination forms; the conversion parameters of different combinations have great difference, so that the two homonymous point combination relations need to be considered respectively to solve the conversion parameters corresponding to each homonymous point combination: a rotation matrix R, a translation vector T and a conversion matrix A between two coordinate systems;
step 5, setting a counter C 1 Corresponding to the first form of
Figure BDA00038899033800000313
And->
Figure BDA00038899033800000314
And->
Figure BDA00038899033800000315
And->
Figure BDA00038899033800000316
Is combined with the homonymy point, and p is converted by using a conversion matrix A of the homonymy point combination s Unification to p t In the coordinate system, next, for I s Each intersection point of (1) t Find out that it is nearest and belongs to I t Crossing of->
Figure BDA00038899033800000317
And calculates the distance d between these two points, once d < the length threshold ε 7 Two points are taken as a pair of homonymous intersection points, < ->
Figure BDA00038899033800000318
As a qualified inner point and to count the counter C 1 The number of inner points of (1) is increased by 1, and I is traversed according to the previous process s Each of (3)After each intersection point, a counter C is obtained 1 The sum of all the accumulated values of the number of the inner points and the homonymous intersection point pair corresponding to each qualified inner point; in this way, a counter C is obtained for the second form of homonymous point combination 2 The sum of the cumulative values of the number of all the qualified interior points and the homonymous intersection point pair corresponding to each qualified interior point; from C 1 and C2 Finding out the corresponding homonymous intersection point pair of the counter meeting the maximum accumulated value of the number of the inner points or exceeding a certain threshold value requirement as the intersection point pair meeting the registration condition for registration, and paying attention to the accumulated value according to the fact that the value of the accumulated value is C 1 and C2 All elements being compared, not separately, and more than one group being combined, in general, with eligible homonyms, and C 1 and C2 The difference value between the element with the largest accumulated value and the following elements is very small; if a plurality of groups of homonymous intersection point pairs meeting registration conditions exist, the singular value decomposition is used for solving the conversion parameters to carry out integral registration, and the obtained rotation matrix and translation vector are used for carrying out P s And P t Aligned on the X-Y plane, and then calculates the translation delta Z, P of the two point clouds along the Z-axis direction s And P t And after the alignment of the X-Y plane and the translation in the Z direction, a final registration result is obtained. />
Further, the method for automatically registering the building point cloud based on the two-dimensional projection line segments, provided by the invention, can also have the following characteristics: in step 4, the conversion parameter solving method is as follows: for the first form of homonymous point combination
Figure BDA0003889903380000041
And (3) with
Figure BDA0003889903380000042
And->
Figure BDA0003889903380000043
And->
Figure BDA0003889903380000044
First, a rotation matrix R is calculated:
Figure BDA0003889903380000045
in the formula ,
Figure BDA0003889903380000046
i＝1、2；
the translation vector T is:
Figure BDA0003889903380000047
the transformation matrix a between the two coordinate systems is:
Figure BDA0003889903380000048
solving the second form homonymous point combination by adopting the same method
Figure BDA0003889903380000049
And->
Figure BDA00038899033800000410
And->
Figure BDA00038899033800000411
And->
Figure BDA00038899033800000412
Is used for the conversion parameters of (a).
Further, the method for automatically registering the building point cloud based on the two-dimensional projection line segments, provided by the invention, can also have the following characteristics: in step 4, at least three pairs of intersection points which are possible to be homonymous points and virtual points attached to the intersection points are obtained; each pair of intersection points and virtual points attached to the intersection points are obtained and correspond to two combination forms, and corresponding two groups of conversion parameters are obtained through solving; in step 5, two types of counter C with the same name point combination are obtained for the two groups of conversion parameters corresponding to each pair of intersection points obtained in step 4 1 and C2 The method comprises the steps of carrying out a first treatment on the surface of the The sub-counter corresponding to the kth group of homonymous point combinations in the first form homonymous point combination is marked as an element
Figure BDA0003889903380000051
1≤k≤n 3 ，n 3 The number of all homonymous point combinations contained in the homonymous point combinations in the first form; the second form employs the same method; then from all C 1 and C2 And finding out homonymous intersection point pairs corresponding to elements with the largest cumulative value of the number of the inner points or exceeding a certain threshold value as homonymous intersection point pairs meeting registration conditions for registration.
Further, the method for automatically registering the building point cloud based on the two-dimensional projection line segments, provided by the invention, can also have the following characteristics: in step 5, if there are multiple sets of homonymous intersection pairs meeting the registration condition, to solve the coordinate transformation in the presence of multiple sets of homonymous point combinations, we set
Figure BDA0003889903380000052
and />
Figure BDA0003889903380000053
Respectively P t and Ps The u-th pair of homonymous point sets, the rotation matrix is R, the translation vector is T, and then:
Figure BDA0003889903380000054
order the
Figure BDA0003889903380000055
and />
Figure BDA0003889903380000056
Centroid of +.>
Figure BDA0003889903380000057
and />
Figure BDA0003889903380000058
Will->
Figure BDA0003889903380000059
and />
Figure BDA00038899033800000510
Centralizing, including:
Figure BDA00038899033800000511
in consideration of the influence of errors, the conversion parameters are solved in a mode of reducing error functions, and the total error equation in coordinate conversion is as follows:
Figure BDA00038899033800000512
where n is the logarithm of the homonymy point.
The total error equation after centering is:
Figure BDA00038899033800000513
/>
when the value of delta is minimum, namely corresponding to the rotation matrix R and the translation vector T, the delta is appropriately transformed:
Figure BDA00038899033800000514
therefore, finding the minimum value of Δ is equivalent to finding the maximum value of Δ' which is:
Figure BDA0003889903380000061
wherein ,
Figure BDA0003889903380000062
singular value decomposition is performed on H, with:
[U ∑ V]＝SVD(H) (5-7)
the rotation matrix R can be expressed as:
R＝VU T (5-8)
after calculating the value of the rotation matrix R, the translation vector T can be found from equation (12) as:
Figure BDA0003889903380000063
the transformation matrix a between the two coordinate systems is:
Figure BDA0003889903380000064
using the rotation parameter obtained to determine r s And r t Aligned on the X-Y plane, and then, the translation amount delta Z of the two point clouds along the Z-axis direction is calculated.
Further, the method for automatically registering the building point cloud based on the two-dimensional projection line segments, provided by the invention, can also have the following characteristics: the translational amount deltaz is obtained by the following steps: counting the average elevation h of the flat ground of the two point clouds 1 and h2 Will h 1 And h 2 The difference is taken as DeltaZ; or using the difference in elevation measured by the scanner's own GPS as deltaz.
Further, the method for automatically registering the building point cloud based on the two-dimensional projection line segments, provided by the invention, can also have the following characteristics: in step 5, when the plurality of groups of homonymous intersection point pairs meeting the registration condition are utilized for integral registration, homonymous intersection point pairs are uniformly selected in a region, and then the conversion parameters are obtained for the selected homonymous point combinations by singular value decomposition.
Further, the method for automatically registering the building point cloud based on the two-dimensional projection line segments, provided by the invention, can also have the following characteristics: in step 5, registration is performed based on the homonymic intersection point pair corresponding to the counter with the largest number of internal points, and if the number of the counter with the largest number of internal points is multiple, one counter is selected randomly.
Further, the method for automatically registering the building point cloud based on the two-dimensional projection line segments, provided by the invention, can also have the following characteristics: in step 5 ε 7 Take 0.01m.
In addition, in the method for automatically registering building point clouds based on two-dimensional projection line segments, the following method can be adopted for preprocessing in step 1 to obtain
Figure BDA0003889903380000071
and />
Figure BDA0003889903380000072
Let the target point cloud be P t The point cloud to be registered is P s For P respectively t and Ps Performing spatially uniform downsampling, and setting a length threshold epsilon 1 As the minimum sampling interval, the distance between any two points in the point cloud after downsampling is not less than epsilon 1 Thereby obtaining the down-sampled point cloud +.>
Figure BDA0003889903380000073
and />
Figure BDA0003889903380000074
Next, the +.A. are counted separately according to a certain step size>
Figure BDA0003889903380000075
and />
Figure BDA0003889903380000076
The distribution histogram of the points in the Z direction is automatically removed from the main horizontal plane in the Z direction according to the distribution rule of the histogram, and the +.>
Figure BDA0003889903380000077
and />
Figure BDA0003889903380000078
The distribution histogram of the points in the Z direction is that the points of the ground and the ceiling have two main peaks in the Z direction, and according to the remarkable characteristic, the main horizontal plane is removed by utilizing the Z direction coordinate value, if the point cloud of a certain layer is obviously larger than a plurality of times (such as 2 times, etc.) of the average value, the layer is considered to have the main horizontal plane and is removed, thereby obtaining the filtered point cloud>
Figure BDA0003889903380000079
and />
Figure BDA00038899033800000710
For->
Figure BDA00038899033800000711
and />
Figure BDA00038899033800000712
Projection is carried out, and a projection equation is as follows:
ax+by+cz+d＝0 (1)
where a=b=d=0 and c=1, that is,
Figure BDA00038899033800000713
and />
Figure BDA00038899033800000714
Is projected to the X-Y plane, so as to obtain the projected two-dimensional point cloud +.>
Figure BDA00038899033800000715
and />
Figure BDA00038899033800000716
In addition, in the method for automatically registering building point clouds based on two-dimensional projection line segments, the step 2 can extract L in the following way t and Ls : respectively extracting by using random sampling consistency algorithm (RANSAC)
Figure BDA00038899033800000717
and />
Figure BDA00038899033800000718
Setting the iteration number as epsilon when executing RANSAC algorithm when extracting the two-dimensional point cloud line segments 2 Setting a length threshold epsilon 3 The length threshold epsilon is set as the furthest distance from the point to the straight line 4 As the radius of the neighbor search during clustering, a numerical threshold epsilon is set ； Setting a length threshold epsilon as the minimum point number contained in a qualified line segment 6 The shortest length of one qualified line segment is L, and the extracted two-dimensional point cloud line segment sets are respectively t and Ls For L t and Ls The point cloud line segments in (a) are arranged in descending order according to the length. To->
Figure BDA0003889903380000081
For example, random from->
Figure BDA0003889903380000082
Two points a (x a ,y a ,z a) and B(xb ,y b ,z b ) Then from A (x a ,y a ,z a) and B(xb ,y b ,z b ) The determined straight line L can be expressed as:
Figure BDA0003889903380000083
assume a plane psi c Through the spatial point P (x c ,y c ,z c ) And by
Figure BDA0003889903380000084
As its normal vector, let plane ψ c The intersection point with the straight line AB is C (x d ,y d ,z d ) Then point C (x d ,y d ,z d ) Can be calculated as follows:
Figure BDA0003889903380000085
based on the fact CP ζab, t can be expressed as:
Figure BDA0003889903380000086
then point P (x p ,y p ,z p ) To point C (x c ,y c ,z c ) Distance d of (2) pc The method comprises the following steps:
Figure BDA0003889903380000087
setting the maximum iteration number as epsilon 2 Setting a threshold epsilon 3 (e.g., 0.01 m) as the furthest distance from the point to the line, if d pc ≤ε 3 Regarding the corresponding point as a straight line inner point, and regarding the point cloud to be processed
Figure BDA0003889903380000088
One point p of (a) i A k-nearest neighbor search is performed to find the corresponding neighbor point (denoted +.>
Figure BDA0003889903380000089
) D is then aver The following can be calculated:
Figure BDA00038899033800000810
wherein K is
Figure BDA00038899033800000811
I is more than or equal to 1 and less than or equal to K. For example, n may be set to 2./>
Every time a collinear point P is extracted col After that, P is established col And establishes an empty cluster set E and an empty queue QWill then P col Is stored in Q; for each point Q in Q i Conducting a k-nearest neighbor search and storing the search results in Q i k In calculating Q i k Each point in (a) to Q i Distance d of (2) ik Setting a length threshold epsilon 4 As the furthest distance between two adjacent points, if d i k ≤ε 4 The corresponding point will be equal to Q i Together stored in cluster set E; calculating the distance between clusters in E according to equation (6), and combining the clusters at a distance of not more than ε 4 Until the distance between any two clusters is greater than epsilon 4 The method comprises the steps of carrying out a first treatment on the surface of the If the clustering result is not null, storing the result in P seg In calculating P seg Is defined as each subset P of seg i Length L of (2) seg i Sum of points N seg i Setting a threshold epsilon ； As the minimum point number per unit length, a threshold epsilon is set 6 As the shortest length of a pass line segment, if N seg i /L seg i ≥ε ； And L is seg i ≥ε 6 ，P seg i Is considered a qualified line segment, and then P is defined as seg The point corresponding to the qualified line segment in the list is from
Figure BDA0003889903380000091
Removing; will->
Figure BDA0003889903380000092
The rest points in (a) are used as original data, the process is repeated to enter the next loop, the iteration is ended when the clustering is empty for 5 times continuously, on the other hand, if the rest points are only originally +.>
Figure BDA0003889903380000093
The iterative process is also forced to end, and all qualified line segments are stored in the two-dimensional point cloud line segment set L t In the same way +.>
Figure BDA0003889903380000094
Two-dimensional point cloud segment set L in (1) s For L t and Ls The two-dimensional point cloud line segments are arranged in descending order according to the length.
Effects and effects of the invention
The method for automatically registering the building point clouds based on the two-dimensional projection line segments automatically and rapidly carries out target-free registration on the building point clouds data of different visual angles by utilizing the two-dimensional projection line segment characteristics, aims at the problems of low efficiency, poor precision, easy sinking into local optimum, non-convergence and the like commonly existing in the existing method, has particularly obvious advantages when aiming at registering the large-scale building point clouds, and meets practical requirements in the aspects of precision, applicability, robustness, automation level and the like, thereby serving the subsequent applications such as indoor and outdoor integrated modeling, building Information Model (BIM), indoor rapid navigation and positioning and the like.
Drawings
FIG. 1 is a flow chart of a method of building point cloud automatic registration based on two-dimensional projected line segments in accordance with an embodiment of the present invention;
fig. 2 is a schematic diagram of a target point cloud and a source point cloud according to an embodiment of the present invention, where (a) is target point cloud data and (b) is source point cloud data;
FIG. 3 is a flow chart of data preprocessing involved in an embodiment of the present invention;
FIG. 4 is a point distribution histogram of a point cloud along the Z direction according to an embodiment of the present invention;
FIG. 5 is a two-dimensional projection line segment extraction effect diagram according to an embodiment of the present invention;
FIG. 6 is an extracted point cloud segment and intersection point according to an embodiment of the present invention;
fig. 7 is a comparison of the registration result of the method according to the embodiment of the present invention with the original fine registration result, where (a) is the original fine registration result and (b) is the registration result of the method according to the present invention.
Detailed Description
Specific embodiments of a method for automatic registration of building point clouds based on two-dimensional projection line segments according to the present invention are described in detail below with reference to the accompanying drawings.
< example >
As shown in fig. 1, the method for automatically registering building point clouds based on two-dimensional projection line segments provided in the embodiment includes the following steps:
step 1. In order to verify the proposed method, the automatic registration algorithm proposed by the present invention is demonstrated with real point cloud data, as shown in fig. 2, apartment point cloud data in an open point cloud dataset, i.e. the Apartment point cloud data in a data set of open point cloud, i.e. the Apartment point cloud is scanned by using a FARO Focus 3D X330 HDR ground three-dimensional laser scanner, 7 stations are scanned in total and well registered, and scan_00 and scan_01 are selected for experiments, and scan_01 is taken as a target point cloud P t Scan_00 is used as source point cloud P s . Because the original data is well precisely registered, in order to verify the effectiveness of the automatic point cloud registration method based on the two-dimensional projection line segments, rotation and translation on an X-Y plane are applied to scan_00:
Figure BDA0003889903380000111
its inverse transform matrix is:
Figure BDA0003889903380000112
with a target point cloud P t For example, as shown in FIG. 3, for P t According to epsilon 1 Spatially uniform downsampling with a pitch of =0.01m to obtain a downsampled point cloud
Figure BDA0003889903380000113
A total of 626373 points, calculate +.>
Figure BDA0003889903380000114
The minimum value in the Z direction is-25.94 m, the maximum value is-23.20 m, the downsampling point cloud is divided into 27 parts in the Z direction according to the step length of 0.1m, the average point number of each layer is 23199, and the Z direction is countedAs shown in fig. 4, 4 peaks respectively correspond to the ground, bed, corridor top and ceiling in the target point cloud, which is significantly larger than the average point, the four point clouds are removed, and the rest point clouds may have some small horizontal planes, and can be filtered again to obtain filtered point clouds->
Figure BDA0003889903380000115
Projecting the filtered point cloud to an X-Y plane to obtain a projected two-dimensional point cloud +.>
Figure BDA0003889903380000116
Similarly, repeating the above process for the source point cloud to obtain a projected two-dimensional point cloud +.>
Figure BDA0003889903380000117
Step 2, extracting by using improved random sample consensus algorithm (RANSAC)
Figure BDA0003889903380000118
and />
Figure BDA0003889903380000119
Setting the iteration number as epsilon when executing RANSAC algorithm when extracting the two-dimensional point cloud line segments 2 =200, set a length threshold ε 3 =0.01m as the furthest distance from the point to the straight line, a length threshold ε is set 4 =0.05m as radius of neighbor search at clustering, set a numerical threshold ε ； =50 as the minimum number of points contained in a pass line segment, a length threshold epsilon is set 6 =0.3m as the shortest length of one qualified line segment, assuming that the extracted two-dimensional point cloud line segment sets are L respectively t and Ls For L t and Ls The point cloud segments in (a) are arranged in descending order according to the length, as shown in fig. 5.
Step 3, calculating L respectively t and Ls Qualified intersection points corresponding to the line segments and the affiliated virtual points thereof, only the intersection points of the line segments are shown in FIG. 6And the partial intersection points are displayed in an enlarged manner, wherein 55 red points (the gray scale of the point in the box indicated by the arrow on the upper left side of the figure corresponds to red) are used for representing the utilization of L t The obtained intersection point is 30 blue points (the gray scale of the point indicated by the arrow at the lower left side corresponds to blue), and L is used s The intersection point is obtained.
Step 4. Comparing theta t Any angle of
Figure BDA0003889903380000121
(1≤i≤n 1 ，n 1 Is I t Number of subsets) and θ s Any angle->
Figure BDA0003889903380000122
(1≤j≤n 2 ，n 2 Is P s Number of medium pass intersections)>
Figure BDA0003889903380000123
If |delta theta|<1 DEG, then consider P t and Ps The corresponding pair of intersection points and the virtual points attached to the intersection points are possibly in the same name point relationship, and the three pairs of intersection points are respectively +.>
Figure BDA0003889903380000124
And->
Figure BDA0003889903380000125
And (3) with
Figure BDA0003889903380000126
And->
Figure BDA0003889903380000127
Because of the correspondence between undefined line segments, there is also a correspondence case where homonymous points are +.>
Figure BDA0003889903380000128
And->
Figure BDA0003889903380000129
And->
Figure BDA00038899033800001210
And->
Figure BDA00038899033800001211
The conversion parameters calculated for the different combinations of the three pairs of points differ considerably, so that the two cases need to be considered separately for the combination +.>
Figure BDA00038899033800001212
And->
Figure BDA00038899033800001213
Figure BDA00038899033800001214
And->
Figure BDA00038899033800001215
And->
Figure BDA00038899033800001216
For example, a rotation matrix R, a translation vector T and a conversion matrix A between two coordinate systems are calculated; likewise, for combinations->
Figure BDA00038899033800001217
And->
Figure BDA00038899033800001218
And->
Figure BDA00038899033800001219
And->
Figure BDA00038899033800001220
A rotation matrix R, a translation vector T and a transformation matrix a between the two coordinate systems are also calculated.
Step 5, utilizing the combination
Figure BDA00038899033800001221
And->
Figure BDA00038899033800001222
And->
Figure BDA00038899033800001223
And->
Figure BDA00038899033800001224
Conversion matrix A will P s Unification to P t In the coordinate system, next, for I s Each intersection of->
Figure BDA00038899033800001225
(1≤i≤n 2 ，n 2 Is P s Number of qualified intersections), from I t Finding the nearest point to it +.>
Figure BDA00038899033800001226
And calculates the distance d between the two points if d < ε 7( wherein ε7 For the length threshold, the length threshold can be set according to practical conditions, and generally 0.01m is taken, then the counter C is used 1 Element->
Figure BDA00038899033800001227
The value of the K is increased by 1 (1 is less than or equal to k is less than or equal to n) 3 ，n 3 For all +.>
Figure BDA00038899033800001228
And->
Figure BDA00038899033800001229
And->
Figure BDA00038899033800001230
And->
Figure BDA00038899033800001231
Form combination and angle difference less than threshold 1 DEG total number of pairs, +>
Figure BDA00038899033800001232
Uniformly set to 0), and the registration residual of each combination is also calculated at the same time of counting. According to this method, the sum of +.>
Figure BDA00038899033800001233
And (3) with
Figure BDA00038899033800001234
And->
Figure BDA00038899033800001235
And->
Figure BDA00038899033800001236
Form-combined counter C 2 And each combined registration residual, for C 1 and C2 The elements of the group are sorted from large to small, the element with the largest accumulated value and the corresponding intersection point are found, more than one group of matching points (each group comprises three pairs of homonymous points) are matched according to the condition, and C 1 and C2 If there are several pairs of qualified paired combinations, singular value decomposition is used to solve parameters, and the obtained rotation and translation parameters are used to solve P s And P t Aligned on the X-Y plane, and finally, calculating the translation delta Z of the two point clouds along the Z-axis direction, wherein a simple mode is to calculate the average elevation h of the ground flat position of the two point clouds 1 and h2 The translation amount can be regarded as h 1 And h 2 The difference can also be measured by using the GPS of the scanner, and the final registration result can be obtained after the rotation of the X-Y plane and the translation of the Z direction.
In this embodiment, in order to compare the accuracy of the final registration result with the case of multiple pairs, in step 4, the intersection point and each additional virtual point (736 homonymous point combinations exist in each of the two combination forms, and 1650 combinations are added) are obtained through traversal, and only one or more pairs of intersection points and each additional virtual point need to be obtained in the actual calculation process, and all possible homonymous points do not need to be obtained through traversal comparison of all angles.
Specifically, in steps 4 and 5: co-determination of
Figure BDA0003889903380000132
And->
Figure BDA0003889903380000133
And->
Figure BDA0003889903380000134
And->
Figure BDA0003889903380000135
The total number of pairs combined in form and having an angle difference smaller than the threshold value of 1 DEG totals 736 pairs, likewise +.>
Figure BDA0003889903380000136
And->
Figure BDA0003889903380000137
And->
Figure BDA0003889903380000138
And->
Figure BDA0003889903380000139
The total number of pairs of formal combinations is also 736 pairs, if no angular constraint is imposed, the total number of pairs of a single combination is 1650 pairs; table 1 shows two counters C 1 and C2 Statistics of C 1 and C2 The numbers in a column represent the logarithm of homonymous intersections (pairs of homonymous intersections) with the same total number of inliers, if a number 1 represents 1 pair of homonymous intersections, and 4 represents 4 pairs of homonymous intersections. It can be seen that the corresponding homonymous intersections with the number of inliers 25, 24 and 23 are most likely to be used for registration, i.e. at least 18 pairs of homonymous intersections are found that can be used for registration.
TABLE 1
Figure BDA0003889903380000131
/>
Figure BDA0003889903380000141
In order to verify the reliability of the algorithm of the invention, the first 5 pairs of intersection points and the virtual points attached to the intersection points are simply taken for testing, so that the method can automatically find out a plurality of correct pairing combinations, and the coordinate values of the 5 pairs of intersection points are shown in a table 2.
TABLE 2
Figure BDA0003889903380000142
Firstly, 5 pairs of intersection points and virtual points attached to the intersection points are used for calculating a conversion matrix respectively:
Figure BDA0003889903380000143
Figure BDA0003889903380000144
Figure BDA0003889903380000145
Figure BDA0003889903380000151
Figure BDA0003889903380000152
/>
registration is then performed using the 5 different transformation matrices obtained above, respectively, and fig. 7 shows the use of transformation matrix a 1 Is compared to the original fine registration result, wherein the red part (light grey part in the figure) represents the target point cloud,blue parts (dark grey parts in the figure) represent source point clouds, and it can be seen that the automatic registration method provided by the invention can achieve a good rough registration effect. For further quantitative evaluation of the registration results, precision assessment was performed using formula (VIII) and formula (IX):
the rotation error is calculated as follows:
Figure BDA0003889903380000153
in the formula ,Rtrue Representing the actual rotation matrix, R representing the calculated rotation matrix.
The translation error is defined as follows:
Figure BDA0003889903380000154
table 3 shows the comparison of registration accuracy of 5 different combinations, and it can be seen that the rotation errors are all within 1.5 degrees, the translation error in the X direction is less than 0.02m, the translation error in the Y direction is less than 0.005m, and the total translation error is less than 0.02m, so that the automatic registration method based on the two-dimensional projection line segment provided by the invention can automatically calculate a plurality of groups of optimal pairing combinations, and has high accuracy and reliability.
TABLE 3 Table 3
Figure BDA0003889903380000155
Figure BDA0003889903380000161
The above embodiments are merely illustrative of the technical solutions of the present invention. The method of automatic registration of building point clouds based on two-dimensional projection line segments according to the present invention is not limited to what has been described in the above embodiments, but is subject to the scope defined by the claims. Any modifications, additions or equivalent substitutions made by those skilled in the art based on this embodiment are within the scope of the invention as claimed in the claims.

Claims (8)

1. The method for automatically registering the building point cloud based on the two-dimensional projection line segments is characterized by comprising the following steps of:
step 1, setting the target point cloud as P t The point cloud to be registered is P s For P t and Ps All are carried out in sequence: the pretreatment of spatial uniform downsampling, filtering the main horizontal plane in the Z direction and projecting to the X-Y plane to obtain a two-dimensional point cloud
Figure FDA0003889903370000011
and />
Figure FDA0003889903370000012
Step 2, setting the farthest distance from the point to the straight line, the radius of neighbor searching during clustering, the minimum number of points contained in the qualified line segment and the shortest length of the qualified line segment as extraction requirements, and respectively taking the points and the radius of the neighbor searching, the minimum number of points contained in the qualified line segment and the shortest length of the qualified line segment as extraction requirements
Figure FDA0003889903370000013
and />
Figure FDA0003889903370000014
Extracting two-dimensional point cloud line segments meeting the extraction requirements to form a two-dimensional point cloud line segment set L t and Ls For L t and Ls The point cloud line segments in (a) are arranged in descending order according to the length;
step 3, acquiring qualified intersection points, attached virtual points and included angles of the two-dimensional point cloud line segment sets: for L t Firstly, calculating the longest point cloud line segment
Figure FDA0003889903370000015
Straight line and L t Any point cloud segment left->
Figure FDA0003889903370000016
Angle formed by the straight line>
Figure FDA0003889903370000017
If it is
Figure FDA0003889903370000018
When the intersection point is qualified, the intersection point of the two straight lines and the attached two virtual points are calculated and compared: let the intersection point of two intersecting straight lines be p inters (c 1 ,c 2 ) Assuming that a virtual point on a straight line having a length r from the intersection point is p (x, y), the following relationship exists:
Figure FDA0003889903370000019
substituting the first polynomial in equation (7) into the second polynomial, the value of x can be found as:
Figure FDA00038899033700000110
wherein ,
Figure FDA00038899033700000111
each straight line is calculated to obtain two virtual points, and then the two virtual points are calculated to p for each straight line cent Distance d of (2) 1 and d2 Distance p on each straight line cent The nearest virtual point is the virtual point attached to the qualified intersection point, p cent The mass center of the two intersecting point cloud line segments;
once two line segments are compared, marking is carried out to avoid repeated comparison, and similarly, any two line segments are compared to obtain L t All qualified intersection points in the list form an intersection point set I t Each combination ofTwo virtual points are attached to the intersection point of the grid, and the included angle formed by each qualified intersection point and the attached two virtual points is calculated to form an included angle set theta t The method comprises the steps of carrying out a first treatment on the surface of the By using the L t The same method obtains L s Intersection point set I of all qualified intersection points in the list s Two virtual points attached to each qualified intersection point and a corresponding included angle set theta s ；
Step 4. Comparing theta t Any angle of
Figure FDA0003889903370000021
And theta s Any angle->
Figure FDA0003889903370000022
Difference (I) of->
Figure FDA0003889903370000023
n 1 Is I t The number of subsets of (1) j n 2 ，n 2 Is P s Number of pass-through intersections if |Δθ|<1 DEG, the corresponding pair of intersection points and the virtual points attached to the intersection points are considered to be the same name points, and the corresponding same name point relationship is the first +.>
Figure FDA0003889903370000024
And->
Figure FDA0003889903370000025
Figure FDA0003889903370000026
And->
Figure FDA0003889903370000027
Figure FDA0003889903370000028
And->
Figure FDA0003889903370000029
Or second +.>
Figure FDA00038899033700000210
And->
Figure FDA00038899033700000211
Figure FDA00038899033700000212
And->
Figure FDA00038899033700000213
Figure FDA00038899033700000214
And->
Figure FDA00038899033700000215
Two combination forms; respectively considering the two homonymous point combination relations, and solving conversion parameters corresponding to each homonymous point combination: a rotation matrix R, a translation vector T and a conversion matrix A between two coordinate systems;
step 5, setting a counter C 1 Corresponding to the first form of
Figure FDA00038899033700000216
And->
Figure FDA00038899033700000217
Figure FDA00038899033700000218
And->
Figure FDA00038899033700000219
Figure FDA00038899033700000220
And->
Figure FDA00038899033700000221
Is combined with the homonymy point, and p is converted by using a conversion matrix A of the homonymy point combination s Unification to p t In the coordinate system, next, for I s Each intersection point of (1) t Find out that it is nearest and belongs to I t Crossing of->
Figure FDA00038899033700000222
And calculates the distance d between these two points, once d < the length threshold ε 7 Two points are taken as a pair of homonymous intersection points, < ->
Figure FDA00038899033700000223
As a qualified inner point and to count the counter C 1 The number of inner points of (1) is increased by 1, and I is traversed according to the previous process s After each intersection point of (a) a counter C is obtained 1 The sum of all the accumulated values of the number of the inner points and the homonymous intersection point pair corresponding to each qualified inner point; in this way, a counter C is obtained for the second form of homonymous point combination 2 The sum of the cumulative values of the number of all the qualified interior points and the homonymous intersection point pair corresponding to each qualified interior point; from C 1 and C2 Finding out the homonymy intersection point pair corresponding to the counter meeting the maximum accumulated value of the number of the inner points or exceeding a certain threshold value requirement as the intersection point pair meeting the registration condition, if a plurality of groups of homonymy intersection point pairs meeting the registration condition exist, using singular value decomposition to solve the conversion parameter for integral registration, and using the obtained rotation matrix and translation vector to carry out P s And P t Aligned on the X-Y plane, and then calculates the translation delta Z, P of the two point clouds along the Z-axis direction s And P t And after the alignment of the X-Y plane and the translation in the Z direction, a final registration result is obtained.
2. The method for automatic registration of building point clouds based on two-dimensional projection line segments according to claim 1, wherein:
in step 4, the conversion parameter solving method is as follows: for the first form of homonymous point combination
Figure FDA0003889903370000031
And->
Figure FDA0003889903370000032
And (3) with
Figure FDA0003889903370000033
And->
Figure FDA0003889903370000034
First, a rotation matrix R is calculated:
Figure FDA0003889903370000035
in the formula ,
Figure FDA0003889903370000036
the translation vector T is:
Figure FDA0003889903370000037
the transformation matrix a between the two coordinate systems is:
Figure FDA0003889903370000038
and solving the conversion parameters of the second form homonymous point combination by adopting the same method.
3. The method for automatic registration of building point clouds based on two-dimensional projection line segments according to claim 1, wherein:
in step 4, at least three pairs of intersection points which may be homonymous points and virtual points attached to the intersection points are obtained; each pair of intersection points and virtual points attached to the intersection points are obtained and correspond to two combination forms, and corresponding two groups of conversion parameters are obtained through solving;
in step 5, two types of counter C with the same name point combination are obtained for the two groups of conversion parameters corresponding to each pair of intersection points obtained in step 4 1 and C2 The method comprises the steps of carrying out a first treatment on the surface of the The sub-counter corresponding to the kth group of homonymous point combinations in the first form homonymous point combination is marked as an element
Figure FDA0003889903370000041
n 3 The number of all homonymous point combinations contained in the homonymous point combinations in the first form; the second form employs the same method; then from all C 1 and C2 And finding out homonymous intersection point pairs corresponding to elements with the largest cumulative value of the number of the inner points or exceeding a certain threshold value as homonymous intersection point pairs meeting registration conditions for registration.
4. The method for automatic registration of building point clouds based on two-dimensional projection line segments according to claim 1, wherein:
in step 5, if there are multiple groups of homonymous intersection point pairs meeting registration conditions, solving the coordinate transformation relationship, and setting
Figure FDA0003889903370000042
and />
Figure FDA0003889903370000043
Respectively P t and Ps The u-th pair of homonymous intersection sets, the rotation matrix is R, the translation vector is T, and then:
Figure FDA0003889903370000044
/>
order the
Figure FDA0003889903370000045
and />
Figure FDA0003889903370000046
Centroid of +.>
Figure FDA0003889903370000047
and />
Figure FDA0003889903370000048
Will->
Figure FDA0003889903370000049
and />
Figure FDA00038899033700000410
Centralizing, including:
Figure FDA00038899033700000411
the total error equation at the time of coordinate conversion is:
Figure FDA00038899033700000412
wherein n is the logarithm of the homonymy point;
the total error equation after centering is:
Figure FDA00038899033700000413
when the value of delta is minimum, namely corresponding to the rotation matrix R and the translation vector T, the delta is appropriately transformed:
Figure FDA00038899033700000414
therefore, finding the minimum value of Δ is equivalent to finding the maximum value of Δ' which is:
Figure FDA0003889903370000051
wherein ,
Figure FDA0003889903370000052
singular value decomposition is performed on H, with:
[U ∑ V]＝SVD(H) (5-7)
the rotation matrix R can be expressed as:
R＝VU T (5-8)
after calculating the value of the rotation matrix R, the translation vector T can be found from equation (12) as:
Figure FDA0003889903370000053
the transformation matrix a between the two coordinate systems is:
Figure FDA0003889903370000054
using the rotation parameter obtained to determine r s And r t Aligned on the X-Y plane, and then, the translation amount delta Z of the two point clouds along the Z-axis direction is calculated.
5. The method for automatic registration of building point clouds based on two-dimensional projection line segments according to claim 1, wherein:
the method for obtaining the translation amount delta Z comprises the following steps: counting the average elevation h of the flat ground of the two point clouds 1 and h2 Will h 1 And h 2 The difference is taken as DeltaZ; or using the difference in elevation measured by the scanner's own GPS as deltaz.
6. The method for automatic registration of building point clouds based on two-dimensional projection line segments according to claim 1, wherein:
in step 5, when the plurality of groups of homonymous intersection point pairs meeting the registration condition are utilized for integral registration, homonymous intersection point pairs are uniformly selected in a detection area, and then the conversion parameters are obtained for the selected homonymous intersection point pairs by singular value decomposition.
7. The method for automatic registration of building point clouds based on two-dimensional projection line segments according to claim 1, wherein:
in step 5, registration is performed based on the same-name intersection point pair corresponding to the counter with the largest number of internal points, and if there are a plurality of counters with the largest number of internal points, one counter is selected randomly.
8. The method for automatic registration of building point clouds based on two-dimensional projection line segments according to claim 1, wherein:
wherein, in step 5, ε 7 Take 0.01m.
CN202211256855.1A 2022-10-14 2022-10-14 Building point cloud automatic registration method based on two-dimensional projection line segments Pending CN116091562A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211256855.1A CN116091562A (en) 2022-10-14 2022-10-14 Building point cloud automatic registration method based on two-dimensional projection line segments

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211256855.1A CN116091562A (en) 2022-10-14 2022-10-14 Building point cloud automatic registration method based on two-dimensional projection line segments

Publications (1)

Publication Number Publication Date
CN116091562A true CN116091562A (en) 2023-05-09

Family

ID=86206978

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211256855.1A Pending CN116091562A (en) 2022-10-14 2022-10-14 Building point cloud automatic registration method based on two-dimensional projection line segments

Country Status (1)

Country Link
CN (1) CN116091562A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN118089746A (en) * 2024-04-26 2024-05-28 嘉兴新生纪智能科技有限公司 Laser repositioning method, system and storage medium for indoor scene

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN118089746A (en) * 2024-04-26 2024-05-28 嘉兴新生纪智能科技有限公司 Laser repositioning method, system and storage medium for indoor scene

Similar Documents

Publication Publication Date Title
CN109410256B (en) Automatic high-precision point cloud and image registration method based on mutual information
CN111815757B (en) Large member three-dimensional reconstruction method based on image sequence
CN109410321B (en) Three-dimensional reconstruction method based on convolutional neural network
CN109887015B (en) Point cloud automatic registration method based on local curved surface feature histogram
Yi et al. Urban building reconstruction from raw LiDAR point data
CN104331699B (en) A kind of method that three-dimensional point cloud planarization fast search compares
Xu et al. Reconstruction of scaffolds from a photogrammetric point cloud of construction sites using a novel 3D local feature descriptor
CN103632366B (en) A kind of parameter identification method of ellipse target
CN111524168B (en) Point cloud data registration method, system and device and computer storage medium
CN106651942A (en) Three-dimensional rotation and motion detecting and rotation axis positioning method based on feature points
CN111080684A (en) Point cloud registration method for point neighborhood scale difference description
CN108921939A (en) A kind of method for reconstructing three-dimensional scene based on picture
CN105574527A (en) Quick object detection method based on local feature learning
CN114972459B (en) Point cloud registration method based on low-dimensional point cloud local feature descriptor
CN108388902B (en) Composite 3D descriptor construction method combining global framework point and local SHOT characteristics
CN111553292A (en) Rock mass structural plane identification and occurrence classification method based on point cloud data
CN103295239A (en) Laser-point cloud data automatic registration method based on plane base images
CN111028292A (en) Sub-pixel level image matching navigation positioning method
CN116379915A (en) Building mapping method, device, system and storage medium
Jiang et al. Learned local features for structure from motion of uav images: A comparative evaluation
Yuan et al. 3D point cloud recognition of substation equipment based on plane detection
CN108898629B (en) Projection coding method for enhancing aerial luggage surface texture in three-dimensional modeling
CN116452604B (en) Complex substation scene segmentation method, device and storage medium
CN112581511A (en) Three-dimensional reconstruction method and system based on approximate vertical scanning point cloud rapid registration
CN106228593B (en) A kind of image dense Stereo Matching method

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