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 PDFInfo
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 68
- 230000001186 cumulative effect Effects 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 32
- 238000013519 translation Methods 0.000 claims description 30
- 238000006243 chemical reaction Methods 0.000 claims description 26
- 230000009466 transformation Effects 0.000 claims description 10
- 238000000605 extraction Methods 0.000 claims description 8
- 238000000354 decomposition reaction Methods 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 7
- 238000001914 filtration Methods 0.000 claims description 2
- 238000001514 detection method Methods 0.000 claims 1
- 238000007781 pre-processing Methods 0.000 abstract description 3
- 230000000694 effects Effects 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 241000238097 Callinectes sapidus Species 0.000 description 1
- 241000350052 Daniellia ogea Species 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 239000013256 coordination polymer Substances 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011158 quantitative evaluation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
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 cloudAndextracting 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
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 and />
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 segmentStraight line and L t Any point cloud segment left->Angle formed by the straight line>If it isWhen 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: />
Substituting the first polynomial in equation (7) into the second polynomial, the value of x can be found as:
wherein ,
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 ofAnd theta s Any angle->Difference (I) of->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 +.>And->And->And (3) withOr second +.>And->And->And->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 ofAnd->And->And->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->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, < ->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 combinationAnd (3) withAnd->And->First, a rotation matrix R is calculated:
the translation vector T is:
the transformation matrix a between the two coordinate systems is:
solving the second form homonymous point combination by adopting the same methodAnd->And->And->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 element1≤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 and />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:
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:
where n is the logarithm of the homonymy point.
The total error equation after centering is:
when the value of delta is minimum, namely corresponding to the rotation matrix R and the translation vector T, the delta is appropriately transformed:
therefore, finding the minimum value of Δ is equivalent to finding the maximum value of Δ' which is:
wherein ,
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:
the transformation matrix a between the two coordinate systems is:
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 and />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 +.> and />Next, the +.A. are counted separately according to a certain step size> and />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 +.> and />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> and />For-> and />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, and />Is projected to the X-Y plane, so as to obtain the projected two-dimensional point cloud +.> and />
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) and />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->For example, random from->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:
assume a plane psi c Through the spatial point P (x c ,y c ,z c ) And byAs 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:
based on the fact CP ζab, t can be expressed as:
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:
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 processedOne point p of (a) i A k-nearest neighbor search is performed to find the corresponding neighbor point (denoted +.>) D is then aver The following can be calculated:
wherein K isI 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 fromRemoving; will->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 +.>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 +.>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:
its inverse transform matrix is:
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 cloudA total of 626373 points, calculate +.>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->Projecting the filtered point cloud to an X-Y plane to obtain a projected two-dimensional point cloud +.>Similarly, repeating the above process for the source point cloud to obtain a projected two-dimensional point cloud +.>
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(1≤i≤n 1 ,n 1 Is I t Number of subsets) and θ s Any angle->(1≤j≤n 2 ,n 2 Is P s Number of medium pass intersections)>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 +.>And->And (3) withAnd->Because of the correspondence between undefined line segments, there is also a correspondence case where homonymous points are +.>And->And->And->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 +.>And-> And->And->For example, a rotation matrix R, a translation vector T and a conversion matrix A between two coordinate systems are calculated; likewise, for combinations->And->And->And->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 combinationAnd->And->And->Conversion matrix A will P s Unification to P t In the coordinate system, next, for I s Each intersection of->(1≤i≤n 2 ,n 2 Is P s Number of qualified intersections), from I t Finding the nearest point to it +.>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->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 +.>And->And->And->Form combination and angle difference less than threshold 1 DEG total number of pairs, +>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 +.>And (3) withAnd->And->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 ofAnd->And->And->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 +.>And->And->And->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
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
Firstly, 5 pairs of intersection points and virtual points attached to the intersection points are used for calculating a conversion matrix respectively:
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:
in the formula ,Rtrue Representing the actual rotation matrix, R representing the calculated rotation matrix.
The translation error is defined as follows:
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
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 and />
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 and />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 segmentStraight line and L t Any point cloud segment left->Angle formed by the straight line>If it isWhen 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:
substituting the first polynomial in equation (7) into the second polynomial, the value of x can be found as:
wherein ,
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 ofAnd theta s Any angle->Difference (I) of->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 +.>And-> And-> And->Or second +.>And-> And-> And->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 ofAnd-> And-> And->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->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, < ->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 combinationAnd->And (3) withAnd->First, a rotation matrix R is calculated:
the translation vector T is:
the transformation matrix a between the two coordinate systems is:
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 elementn 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 and />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:
the total error equation at the time of coordinate conversion is:
wherein n is the logarithm of the homonymy point;
the total error equation after centering is:
when the value of delta is minimum, namely corresponding to the rotation matrix R and the translation vector T, the delta is appropriately transformed:
therefore, finding the minimum value of Δ is equivalent to finding the maximum value of Δ' which is:
wherein ,
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:
the transformation matrix a between the two coordinate systems is:
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.
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)
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 |
-
2022
- 2022-10-14 CN CN202211256855.1A patent/CN116091562A/en active Pending
Cited By (1)
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 |