Point cloud integrity improvement method and system based on Delaunay triangle and minimum curvature

Abstract

The invention discloses a Delaunay triangle and curvature minimization-based point cloud integrity improvement method and system, and the method comprises the following steps: converting a point cloud under a three-dimensional rectangular coordinate system into a spherical coordinate system, and then projecting the point cloud to a two-dimensional plane coordinate system to form a two-dimensional plane point cloud; carrying out Delaunay triangularization, and then projecting back to a three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular mesh; the side length, the area and the unit normal vector of each triangular mesh are calculated, and abnormal triangular meshes are removed according to the side length; searching boundary points of the point cloud according to the topological relation of the triangular mesh edges; searching the boundary points of the missing part according to the corresponding relation of the boundary points; filling a plurality of points in the missing part based on spherical projection; and calculating the curvature of the filling point, and performing iteration to enable the curvature to be minimum so as to correct the position of the filling point. According to the method, the boundary points of the point cloud can be effectively detected, the filling points are generated, and the missing parts are repaired, so that the integrity and the accuracy of the point cloud are improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of point cloud processing, and in particular relates to a method and system for improving point cloud integrity based on Delaunay triangle and minimum curvature. Background Art

[0002] As an important part of transportation infrastructure, the health of bridges is directly related to traffic safety and smoothness. With the development of science and technology, 3D laser scanning technology has been widely used in bridge inspection due to its advantages such as high precision and non-contact measurement. The point cloud data of bridge components obtained by 3D laser scanners can accurately reflect the geometric shape and surface characteristics of the bridge.

[0003] However, in practical applications, there are some problems with bridge component point clouds. On the one hand, point clouds have the characteristic of large density variations, which makes it impossible to rely solely on spatial distance when establishing the topological relationship between points. No matter what boundary detection method is used, it is challenging. On the other hand, the original point cloud model may have the problem of missing edges, which needs to be repaired to generate a complete point cloud model.

[0004] Most of the existing point cloud repair methods are aimed at studying point clouds of small size, arbitrarily complex surfaces, and uniform point density. For point clouds of bridge components that are often in the millions or even tens of millions, the efficiency of many algorithms is affected and they are not suitable for processing such large-scale point clouds.

[0005] Therefore, it is urgent to solve the above problems. Summary of the invention

[0006] Purpose of the invention: The purpose of the present invention is to provide a point cloud integrity improvement method based on Delaunay triangulation and minimum curvature, which can effectively detect the boundary points of the point cloud, generate filling points and repair the missing parts to improve the integrity and accuracy of the point cloud.

[0007] The second object of the present invention is to provide a point cloud integrity improvement system based on Delaunay triangulation and minimum curvature.

[0008] Technical solution: To achieve the above objectives, the present invention discloses a point cloud integrity improvement method based on Delaunay triangle and minimum curvature, comprising the following steps:

[0009] S100, transforming the point cloud in the three-dimensional rectangular coordinate system into a spherical coordinate system, and then projecting it onto a two-dimensional plane coordinate system to form a two-dimensional plane point cloud;

[0010] S200, performing Delaunay triangulation on the two-dimensional plane point cloud, and then projecting it back to a three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular mesh;

[0011] S300, calculating the side length, area and unit normal vector of each triangular mesh, and eliminating abnormal triangular meshes according to the side length;

[0012] S400, searching for boundary points of the point cloud according to the topological relationship of the triangular mesh edges;

[0013] S500, searching for boundary points of the missing part according to the corresponding relationship of the boundary points;

[0014] S600, filling a number of points in the missing part based on spherical projection;

[0015] S700, calculating the curvature of the filling point, and iterating to minimize the curvature to correct the position of the filling point.

[0016] Optionally, the step S100 specifically includes the following steps: according to the principle of a three-dimensional laser scanner, the point cloud can be converted into a spherical coordinate system with the instrument as the origin; referring to the idea of ​​parametric surfaces in classical differential geometry, the horizontal angle α and the vertical angle β of the spherical coordinate system are used as the horizontal axis and the vertical axis of the two-dimensional rectangular coordinate system, thereby projecting the three-dimensional point cloud onto a two-dimensional plane.

[0017] Optionally, the step S200 specifically includes the following steps: performing Delaunay triangulation on the two-dimensional plane point cloud based on the Bowyer-Watson method, and then projecting it back to a three-dimensional rectangular coordinate system to achieve two-dimensional grid division of the three-dimensional point cloud and form a spatial two-dimensional triangular mesh, each of which contains three vertices and a connection relationship between them; specifically, first generating an initial grid from a given point set, and then adding the remaining points to the initial grid in turn, finding a triangle in the original triangular mesh whose circumscribed circle contains the newly inserted point, called the influencing triangle of the newly inserted point, deleting the common edges of the triangle, connecting the newly inserted point with all the vertices of the influencing triangle to generate a new triangle, and optimizing according to the Delaunay standard until all points are added; after performing Delaunay triangulation on the two-dimensional point cloud, projecting it back to a three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular mesh, achieving two-dimensional grid division of the three-dimensional point cloud.

[0018] Optionally, the step S300 specifically includes the following steps: calculating the side length, area and unit normal vector of each triangular mesh by a formula, setting the side length threshold to be less than 0.5m, and the ratio of the longest side to the shortest side of the triangular mesh cannot exceed 25, and deleting abnormal meshes;

[0019] e1=(x2,y2,z2)-(x3,y3,z3)

[0020] e2=(x1,y1,z1)-(x3,y3,z3)

[0021] e3=(x2,y2,z2)-(x1,y1,z1)

[0022]

[0023] Where e1, e2, e3 are the three side lengths of the triangular mesh, S is the area of ​​the triangular mesh, and N is the unit normal vector of the triangular mesh.

[0024] Optionally, the step S500 specifically includes the following steps: for missing points caused by occlusion relationships, the boundary points are divided into several groups according to the boundary lines they belong to, based on the connection relationship of the edge grid, and corresponding points are searched for based on the angle difference in these groups; if the average distance between the corresponding points is less than a certain threshold, it can be considered that the points of the two groups constitute a pair of obstacle boundaries and missing boundaries, where the group farthest from the instrument is the missing boundary, and the boundary points of the missing part are found.

[0025] Optionally, the step S600 specifically includes the following steps: according to the scanning principle of the three-dimensional laser scanner, the acquisition of the point cloud is based on the spherical coordinate system, and the angle at which the instrument emits the laser is not affected by the target object, but is only related to the state of the instrument itself; the spherical coordinate system has three parameters, the horizontal angle α, the vertical angle β and the distance S. When S is ignored, the point cloud becomes a sphere with uniform density in the spherical coordinate system. The missing parts on the spherical surface are interpolated according to the scanning resolution to complete the preliminary filling of the missing parts.

[0026] Optionally, the step S700 specifically includes the following steps:

[0027] S701. Establish a local coordinate system:

[0028] In a three-dimensional rectangular coordinate system, some surfaces are difficult to represent. Therefore, it is more convenient to calculate by establishing a local coordinate system and representing the surface as a parametric surface. Assume:

[0029]

[0030] Where x, y, z are all differentiable functions of u and v, the domain (u, v) is a region D in two-dimensional space, (u, v) = (x(u, v), y(u, v), z(u, v) gives a mapping from two-dimensional space D to three-dimensional space. Under this mapping, the range of D constitutes a surface D in three-dimensional space.

[0031] A point P in the point cloud i (x i ,y i ,z i ) to establish a local coordinate system (u, v, w): point P i (x i ,yi ,z i ) is the origin, the w axis and point P i The directions of the normal vectors are consistent, the u-axis and the v-axis are orthogonal to each other in the tangent plane, and form a rectangular coordinate system with the w-axis;

[0032] Move the origin O of the original global coordinate system to P i At point , the transformation matrix T is:

[0033]

[0034] Rotate around the x-axis by an angle of α, the rotation matrix R x for:

[0035]

[0036] in

[0037] Rotate angle β around the Y axis, rotation matrix R y for:

[0038]

[0039] in

[0040] S702, quadratic parametric surface fitting: the local shape of any surface can be approximated by a quadratic surface,

[0041] S(u,v)=(u,v,w(u,v))

[0042] w(u,v)=au 2 +buv+cv 2 +eu+fu

[0043] When the coefficients e and f of the linear term are 0, the quadratic surface is simplified to a quadratic parabola:

[0044] w(u,v)=au 2 +buv+cv 2

[0045] The quadratic parabola is used to fit the point cloud locally; for any point P in the point cloud i , search for its neighboring point P j ∈N(P i ), P i The coordinate value in the local coordinate system is (u i ,v i ,w i ), its transformation matrix B is:

[0046]

[0047] Point p i The neighboring point p j The coordinates in the local coordinate system (u j ,v j ,w j ) into the formula:

[0048] w(u,v)=au 2 +buv+cv 2

[0049] Obtain an overdetermined set of equations, solve them using the least squares method, and obtain the parameters a, b, c of the best fitting quadratic parabola;

[0050] S703, curvature calculation: calculating the normal vector, the first basic form and the second basic form of the surface;

[0051] The surface is represented by a parametric equation: r(u,v)=(x(u,v),y(u,v),z(u,v)), where u and v are parameters. The normal vector of the surface at the point x(u,v),y(u,v),z(u,v) is:

[0052]

[0053] in and are two tangent vectors of the surface, and the cross product of the two gives the normal vector n;

[0054] The first basic form is denoted as ds 2 =Edu 2 +2Fdudv+Gdv 2 , where E = S u ·S u , F = S u ·S v , G=S v ·S v , S is the parametric equation of the quadratic parabola;

[0055] The second basic form is expressed as II = Ldu 2 +2Mdudv+Ndv 2 , where L = S uu n, M = S uu n, N = S uu n,

[0056] For point P i , whose coordinates in the local coordinate system are (0,0,0), we can find the quadratic parabola at point P i The first-order partial derivative at is:

[0057] Su | (0,0) =(1,0,w u )=(1,0,2cu+bu)=(1,0,0)

[0058] S u | (0,0) =(1,0,w v )=(1,0,2cv+bu)=(0,1,0)

[0059] The quadratic parabola at point P i The second-order partial derivatives and mixed partial derivatives at are:

[0060] S uu | (0,0) =(0,0,w uu )=(0,0,b)

[0061] S vv | (0,0) =(0,0,w vv )=(1,0,2c)

[0062] According to the formula K=4ac-b 2 , H=a+c calculates the Gaussian curvature K and the average curvature H;

[0063] S704, repairing the missing parts of the point cloud: After completing the point filling of the missing parts of the point cloud, a number of filling points are obtained: P t (Rcosθcosα, Rcosθsinα, Rsinθ), where the horizontal angle α and the vertical angle θ are known, and there is only one unknown quantity for the filling point, that is, the distance R from the point to the instrument; select the filling point closest to the reference point as the seed point, and t Select the R value of the nearest reference point as the initial value R0 of the iteration, iterate several steps in the positive and negative directions with a fixed step size, calculate the Gaussian curvature of each step, and select the R value that minimizes the absolute value of the Gaussian curvature. min As the R value of the seed point; after the position of the seed point is determined, the position of the seed point is used as a reference point for the remaining filling points to calculate the curvature, and then the next filling point closest to the reference point is selected as the seed point, until the positions of all filling points are determined to complete the repair of the missing part.

[0064] Based on the same inventive concept, the present invention discloses a point cloud integrity improvement system based on Delaunay triangulation and minimum curvature, including a two-dimensional conversion module for converting a point cloud in a three-dimensional rectangular coordinate system into a spherical coordinate system, and then projecting it onto a two-dimensional plane coordinate system to form a two-dimensional plane point cloud;

[0065] The triangulation module is used to perform Delaunay triangulation on the two-dimensional plane point cloud and then project it back to the three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular mesh;

[0066] The abnormal elimination module is used to calculate the side length, area and unit normal vector of each triangular mesh, and eliminate abnormal triangular meshes according to the side length;

[0067] The boundary query module is used to find the boundary points of the point cloud according to the topological relationship of the triangular mesh edge, and find the boundary points of the missing part according to the corresponding relationship of the boundary points;

[0068] A filling module is used to fill in a number of points in the missing part based on spherical projection;

[0069] The correction module calculates the curvature of the filling point and iteratively minimizes the curvature to correct the position of the filling point.

[0070] Optionally, the abnormal removal module calculates the side length, area and unit normal vector of each triangular mesh by a formula, sets the side length threshold to be less than 0.5m, and the ratio of the longest side to the shortest side of the triangular mesh cannot exceed 25, and deletes the abnormal mesh;

[0071] e1=(x2,y2,z2)-(x3,y3,z3)

[0072] e2=(x1,y1,z1)-(x3,y3,z3)

[0073] e3=(x2,y2,z2)-(x1,y1,z1)

[0074]

[0075] Where e1, e2, e3 are the three side lengths of the triangular mesh, S is the area of ​​the triangular mesh, and N is the unit normal vector of the triangular mesh.

[0076] Optionally, a local coordinate system is established in the correction module:

[0077] In a three-dimensional rectangular coordinate system, some surfaces are difficult to represent. Therefore, it is more convenient to calculate by establishing a local coordinate system and representing the surface as a parametric surface. Assume:

[0078]

[0079] Where x, y, z are all differentiable functions of u and v, the domain (u, v) is a region D in two-dimensional space, (u, v) = (x(u, v), y(u, v), z(u, v) gives a mapping from two-dimensional space D to three-dimensional space. Under this mapping, the range of D constitutes a surface D in three-dimensional space.

[0080] A point P in the point cloud i (x i ,y i ,z i ) to establish a local coordinate system (u, v, w): point P i (x i ,y i ,z i ) is the origin, the w axis and point P i The directions of the normal vectors are consistent, the u-axis and the v-axis are orthogonal to each other in the tangent plane, and form a rectangular coordinate system with the w-axis;

[0081] Move the origin O of the original global coordinate system to P i At point , the transformation matrix T is:

[0082]

[0083] Rotate around the x-axis by an angle of α, the rotation matrix R x for:

[0084]

[0085] in

[0086] Rotate angle β around the Y axis, rotation matrix R y for:

[0087]

[0088] in

[0089] Quadratic parametric surface fitting: The local shape of any surface can be approximated by a quadratic surface,

[0090] S(u,v)=(u,v,w(u,v))

[0091] w(u,v)=au 2 +buv+cv 2 +eu+fu

[0092] When the coefficients e and f of the linear term are 0, the quadratic surface is simplified to a quadratic parabola:

[0093] w(u,v)=au 2 +buv+cv 2

[0094] The quadratic parabola is used to fit the point cloud locally; for any point P in the point cloud i , search for its neighboring point P j ∈N(P i ), P iThe coordinate value in the local coordinate system is (u i ,v i ,w i ), its transformation matrix B is:

[0095]

[0096] Point p i The neighboring point p j The coordinates in the local coordinate system (u j ,v j ,w j ) into the formula:

[0097] w(u,v)=au 2 +buv+cv 2

[0098] Obtain an overdetermined set of equations, solve them using the least squares method, and obtain the parameters a, b, c of the best fitting quadratic parabola;

[0099] Curvature calculation: calculate the normal vector, first fundamental form and second fundamental form of the surface;

[0100] The surface is represented by a parametric equation: r(u,v)=(x(u,v),y(u,v),z(u,v)), where u and v are parameters. The normal vector of the surface at the point x(u,v),y(u,v),z(u,v) is:

[0101]

[0102] in and are two tangent vectors of the surface, and the cross product of the two gives the normal vector n;

[0103] The first basic form is denoted as ds 2 =Edu 2 +2Fdudv+Gdv 2 , where E = S u ·S u , F = S u ·S v , G=S v ·S v , S is the parametric equation of the quadratic parabola;

[0104] The second basic form is expressed as II = Ldu 2 +2Mdudv+Ndv 2 , where L = S uu n, M = S uu n, N = S uu n,

[0105] For point P i , whose coordinates in the local coordinate system are (0,0,0), we can find the quadratic parabola at point P i The first-order partial derivative at is:

[0106] S u | (0,0) =(1,0,w u )=(1,0,2cu+bu)=(1,0,0)

[0107] S u | (0,0) =(1,0,w v )=(1,0,2cv+bu)=(0,1,0)

[0108] The quadratic parabola at point P i The second-order partial derivatives and mixed partial derivatives at are:

[0109] S uu | (0,0) =(0,0,w uu )=(0,0,b)

[0110] S vv | (0,0) =(0,0,w vv )=(1,0,2c)

[0111] According to the formula K=4ac-b 2 , H=a+c calculates the Gaussian curvature K and the average curvature H;

[0112] Repair the missing parts of the point cloud: After completing the point filling of the missing parts of the point cloud, several filling points are obtained: P t (Rcosθcosα, Rcosθsinα, Rsinθ), where the horizontal angle α and the vertical angle θ are known, and there is only one unknown quantity for the filling point, that is, the distance R from the point to the instrument; select the filling point closest to the reference point as the seed point, and t Select the R value of the nearest reference point as the initial value R0 of the iteration, iterate several steps in the positive and negative directions with a fixed step size, calculate the Gaussian curvature of each step, and select the R value that minimizes the absolute value of the Gaussian curvature. min As the R value of the seed point; after the position of the seed point is determined, the position of the seed point is used as a reference point for the remaining filling points to calculate the curvature, and then the next filling point closest to the reference point is selected as the seed point, until the positions of all filling points are determined to complete the repair of the missing part.

[0113] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages: the present invention generates a Delaunay triangulated mesh based on spherical projection, thereby improving the accuracy of the mesh, and obtains the boundary of the missing part of the point cloud based on the triangulated mesh; generates filling points at the missing parts, and corrects the positions of the filling points according to the standard of minimum Gaussian curvature of the filling points, and finally completes the repair of the missing parts, thereby improving the integrity of the point cloud; the edge detection of the present invention is more adaptable to changes in point cloud density, and has a higher accuracy rate in edge recognition; for the missing parts of the point cloud on the surface of the bridge component, the repair effect is better, and the algorithm has high accuracy; the present invention is suitable for point cloud processing of various bridge components, such as concrete pier components and concrete variable-section box girder components, and has a better point cloud repair effect on these components. BRIEF DESCRIPTION OF THE DRAWINGS

[0114] Figure 1 A schematic diagram of a system framework in one embodiment of the present invention;

[0115] Figure 2 A schematic diagram of projecting a three-dimensional point cloud into a two-dimensional plane coordinate system in one embodiment of the present invention;

[0116] Figure 3 A schematic diagram of inserting a new node P using the Bowyer-Watson method in one embodiment of the present invention;

[0117] Figure 4 A schematic diagram of finding an influence triangle using the Bowyer-Watson method in one embodiment of the present invention;

[0118] Figure 5 A schematic diagram of deleting common edges using the Bowyer-Watson method in one embodiment of the present invention;

[0119] Figure 6 A schematic diagram of constructing a new triangle using the Bowyer-Watson method in one embodiment of the present invention;

[0120] Figure 7 A schematic diagram of an original triangular mesh generated based on a point cloud of a T-beam component in one embodiment of the present invention;

[0121] Figure 8 A schematic diagram of a deleted triangular mesh generated based on a point cloud of a T-beam component in one embodiment of the present invention;

[0122] Fig. 9 A schematic diagram of a missing boundary of a search point cloud in an embodiment of the present invention;

[0123] Fig.10 FIG. 4 is a schematic diagram of a local coordinate system in an embodiment of the present invention. DETAILED DESCRIPTION

[0124] The technical solution of the present invention is further described below in conjunction with the accompanying drawings.

[0125] Example 1

[0126] like Figure 1 As shown, a point cloud integrity improvement method based on Delaunay triangulation and minimum curvature includes the following steps:

[0127] S100, transforming the point cloud in the three-dimensional rectangular coordinate system into a spherical coordinate system, and then projecting it onto a two-dimensional plane coordinate system to form a two-dimensional plane point cloud;

[0128] Specifically, the following steps are included: According to the principle of 3D laser scanner, the point cloud can be converted to a spherical coordinate system with the instrument as the origin; referring to the idea of ​​parametric surface in classical differential geometry, the horizontal angle α and vertical angle β of the spherical coordinate system are used as the horizontal axis and vertical axis of the 2D rectangular coordinate system, so as to project the 3D point cloud onto a 2D plane, such as Figure 2 As shown in the figure, from the absolute distance of adjacent points, different degrees of deformation will occur after projection, but the relative positions of adjacent points will not change.

[0129] S200, performing Delaunay triangulation on the two-dimensional plane point cloud, and then projecting it back to a three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular mesh;

[0130] The specific steps include: Figure 3 to Figure 6 As shown, the two-dimensional plane point cloud is Delaunay triangulated based on the Bowyer-Watson method, and then projected back to the three-dimensional rectangular coordinate system to realize the two-dimensional grid division of the three-dimensional point cloud and form a spatial two-dimensional triangular grid, each of which contains three vertices and the connection relationship between them; specifically, an initial grid is first generated from a given point set, and then the remaining points are added to the initial grid in turn, and the triangle whose circumscribed circle contains the new inserted point is found in the original triangular grid, which is called the influence triangle of the new inserted point, the common edges of the triangle are deleted, and the new inserted point is connected with all the vertices of the influence triangle to generate a new triangle, and optimized according to the Delaunay standard until all points are added; after the two-dimensional point cloud is Delaunay triangulated, it is projected back to the three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular grid, and the two-dimensional grid division of the three-dimensional point cloud is realized.

[0131] S300, calculating the side length, area and unit normal vector of each triangular mesh, and eliminating abnormal triangular meshes according to the side length;

[0132] The specific steps include: calculating the side length, area and unit normal vector of each triangular mesh by formula, setting the side length threshold to be less than 0.5m, and the ratio of the longest side to the shortest side of the triangular mesh cannot exceed 25, deleting abnormal meshes, and the effects before and after deleting abnormal meshes are shown in the figure. Figure 7 and Figure 8 As shown;

[0133] e1=(x2,y2,z2)-(x3,y3,z3)

[0134] e2=(x1,y1,z1)-(x3,y3,z3)

[0135] e3=(x2,y2,z2)-(x1,y1,z1)

[0136]

[0137] Where e1, e2, e3 are the three side lengths of the triangular mesh, S is the area of ​​the triangular mesh, and N is the unit normal vector of the triangular mesh.

[0138] S400, finding the boundary points of the point cloud according to the topological relationship of the triangular mesh edges; through the characteristics of the triangular mesh, that is, an edge of the edge mesh is isolated and not connected to other meshes, and the two points connected by the edge are boundary points, so as to quickly find all the boundary points in the point cloud.

[0139] S500, searching for boundary points of the missing part according to the corresponding relationship of the boundary points;

[0140] like Fig. 9 As shown in the figure, for the missing problem caused by occlusion, according to the connection relationship of the edge grid, the boundary points are divided into several groups according to the boundary lines they belong to, and the corresponding points are searched on these groups by angle difference; if the average distance of the corresponding points is less than a certain threshold, it can be considered that the points of the two groups constitute a pair of obstacle boundaries and missing boundaries, and the group far away from the instrument is the missing boundary; for the boundary points generated by the resolution relationship, they have obvious linear characteristics and are difficult to form a closed loop, and most of them can be solved by splicing the clouds of multiple measuring stations, so this type of missing can be ignored during detection.

[0141] S600. Fill in a number of points at the missing parts based on spherical projection. According to the scanning principle of the 3D laser scanner, the acquisition of point cloud is based on the spherical coordinate system, and the angle at which the instrument emits laser is not affected by the target object, but is only related to the state of the instrument itself. The spherical coordinate system has three parameters, the horizontal angle α, the vertical angle β and the distance S. If S is ignored, the point cloud becomes a sphere with uniform density in the spherical coordinate system. The missing parts on the spherical surface are interpolated according to the scanning resolution to complete the preliminary filling of the missing parts.

[0142] S700, calculating the curvature of the filling point, and iterating to minimize the curvature to correct the position of the filling point;

[0143] The specific steps are:

[0144] S701. Establish a local coordinate system:

[0145] In a three-dimensional rectangular coordinate system, some surfaces are difficult to represent. Therefore, it is more convenient to calculate by establishing a local coordinate system and representing the surface as a parametric surface. Assume:

[0146]

[0147] Where x, y, z are all differentiable functions of u and v, the domain (u, v) is a region D in two-dimensional space, (u, v) = (x(u, v), y(u, v), z(u, v) gives a mapping from two-dimensional space D to three-dimensional space. Under this mapping, the range of D constitutes a surface D in three-dimensional space.

[0148] A point P in the point cloud i (x i ,y i ,z i ) to establish a local coordinate system (u, v, w): point P i (x i ,y i ,z i ) is the origin, the w axis and point P i The normal vector direction is consistent, the u-axis and the v-axis are orthogonal to each other in the tangent plane, and form a rectangular coordinate system with the w-axis. The local coordinate system is as shown in the attached figure. Fig.10 As shown;

[0149] Move the origin O of the original global coordinate system to P i At point , the transformation matrix T is:

[0150]

[0151] Rotate around the x-axis by an angle of α, the rotation matrix R x for:

[0152]

[0153] in

[0154] Rotate angle β around the Y axis, rotation matrix R y for:

[0155]

[0156] in

[0157] S702, quadratic parametric surface fitting: the local shape of any surface can be approximated by a quadratic surface,

[0158] S(u,v)=(u,v,w(u,v))

[0159] w(u,v)=au 2 +buv+cv 2 +eu+fu

[0160] When the coefficients e and f of the linear term are 0, the quadratic surface is simplified to a quadratic parabola:

[0161] w(u,v)=au 2 +buv+cv 2

[0162] The quadratic parabola is used to fit the point cloud locally; for any point P in the point cloud i , search for its neighboring point P j ∈N(P i ), P i The coordinate value in the local coordinate system is (u i ,v i ,w i ), its transformation matrix B is:

[0163]

[0164] Point p i The neighboring point p j The coordinates in the local coordinate system (u j ,v j ,w j ) into the formula:

[0165] w(u,v)=au 2 +buv+cv 2

[0166] The overdetermined equations are obtained and solved using the least squares method to obtain the parameters a, b, and c of the best fitting quadratic parabola.

[0167] S703, curvature calculation: calculating the normal vector, the first basic form and the second basic form of the surface;

[0168] The surface is represented by a parametric equation: r(u,v)=(x(u,v),y(u,v),z(u,v)), where u and v are parameters. The normal vector of the surface at the point x(u,v),y(u,v),z(u,v) is:

[0169]

[0170] in and are two tangent vectors of the surface, and the cross product of the two gives the normal vector n;

[0171] The first basic form of the line element is expressed as ds 2 =Edu 2 +2Fdudv+Gdv 2 Indicates that E = S u ·S u , F = S u ·S v , G=S v ·S v , S is the parametric equation of the quadratic parabola;

[0172] The second basic form is expressed as II = Ldu 2 +2Mdudv+Ndv 2 , where L = S uu n, M = S uu n, N = S uu n,

[0173] For point P i , whose coordinates in the local coordinate system are (0,0,0), we can find the quadratic parabola at point P i The first-order partial derivative at is:

[0174] S u | (0,0) =(1,0,w u )=(1,0,2cu+bu)=(1,0,0)

[0175] S u | (0,0) =(1,0,w v )=(1,0,2cv+bu)=(0,1,0)

[0176] The quadratic parabola at point P i The second-order partial derivatives and mixed partial derivatives at are:

[0177] S uu | (0,0) =(0,0,w uu )=(0,0,b)

[0178] S vv | (0,0) =(0,0,w vv )=(1,0,2c)

[0179] According to the formula K=4ac-b 2, H=a+c, the Gaussian curvature K and the average curvature H are calculated.

[0180] S704, repairing the missing parts of the point cloud: After completing the point filling of the missing parts of the point cloud, a number of filling points are obtained: P t (Rcosθcosα, Rcosθsinα, Rsinθ), where the horizontal angle α and the vertical angle θ are known, and there is only one unknown quantity for the filling point, that is, the distance R from the point to the instrument; select the filling point closest to the reference point as the seed point, and t Select the R value of the nearest reference point as the initial iteration value R0, iterate 1000 steps in the positive and negative directions with a step size of 0.1 mm, calculate the Gaussian curvature of each step, and select the R value that minimizes the absolute value of the Gaussian curvature. min As the R value of the seed point; after the position of the seed point is determined, the position of the seed point is used as a reference point for the remaining filling points to calculate the curvature, and then the next filling point closest to the reference point is selected as the seed point, until the positions of all filling points are determined to complete the repair of the missing part.

[0181] Example 2: Based on the same inventive concept, the present invention discloses a point cloud integrity improvement system based on Delaunay triangulation and minimum curvature, including a two-dimensional conversion module for converting a point cloud in a three-dimensional rectangular coordinate system into a spherical coordinate system, and then projecting it onto a two-dimensional plane coordinate system to form a two-dimensional plane point cloud.

[0182] According to the principle of 3D laser scanner, the point cloud can be converted to a spherical coordinate system with the instrument as the origin; referring to the idea of ​​parametric surface in classical differential geometry, the horizontal angle α and vertical angle β of the spherical coordinate system are used as the horizontal axis and vertical axis of the 2D rectangular coordinate system, so as to project the 3D point cloud onto the 2D plane, such as Figure 2 As shown in the figure, from the absolute distance of adjacent points, different degrees of deformation will occur after projection, but the relative positions of adjacent points will not change.

[0183] The triangulation module is used to perform Delaunay triangulation on the two-dimensional plane point cloud, and then project it back to the three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular mesh.

[0184] like Figure 3 to Figure 6As shown, the two-dimensional plane point cloud is Delaunay triangulated based on the Bowyer-Watson method, and then projected back to the three-dimensional rectangular coordinate system to realize the two-dimensional grid division of the three-dimensional point cloud and form a spatial two-dimensional triangular grid, each of which contains three vertices and the connection relationship between them; specifically, an initial grid is first generated from a given point set, and then the remaining points are added to the initial grid in turn, and the triangle whose circumscribed circle contains the new inserted point is found in the original triangular grid, which is called the influence triangle of the new inserted point, the common edges of the triangle are deleted, and the new inserted point is connected with all the vertices of the influence triangle to generate a new triangle, and optimized according to the Delaunay standard until all points are added; after the two-dimensional point cloud is Delaunay triangulated, it is projected back to the three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular grid, and the two-dimensional grid division of the three-dimensional point cloud is realized.

[0185] The abnormal elimination module is used to calculate the side length, area and unit normal vector of each triangular mesh and eliminate abnormal triangular meshes according to the side length.

[0186] Calculate the side length, area and unit normal vector of each triangular mesh, and remove abnormal triangular meshes according to the side length; calculate the side length, area and unit normal vector of each triangular mesh by formula, set the side length threshold to be less than 0.5m, and the ratio of the longest side to the shortest side of the triangular mesh cannot exceed 25, delete abnormal meshes, and the effects before and after deleting abnormal meshes are as follows: Figure 7 and Figure 8 As shown;

[0187] e1=(x2,y2,z2)-(x3,y3,z3)

[0188] e2=(x1,y1,z1)-(x3,y3,z3)

[0189] e3=(x2,y2,z2)-(x1,y1,z1)

[0190]

[0191] Where e1, e2, e3 are the three side lengths of the triangular mesh, S is the area of ​​the triangular mesh, and N is the unit normal vector of the triangular mesh.

[0192] The boundary query module is used to find the boundary points of the point cloud according to the topological relationship of the triangular mesh edges, and to find the boundary points of the missing parts according to the corresponding relationship of the boundary points; to find the boundary points of the point cloud according to the topological relationship of the triangular mesh edges; through the characteristics of the triangular mesh, that is, an edge of the edge mesh is isolated and not connected to other meshes, and the two points connected by the edge are boundary points, so as to quickly find all the boundary points in the point cloud.

[0193] Find the boundary points of the missing parts according to the corresponding relationship of the boundary points; Fig. 9 As shown in the figure, for the missing problem caused by occlusion, according to the connection relationship of the edge grid, the boundary points are divided into several groups according to the boundary lines they belong to, and the corresponding points are searched on these groups by angle difference; if the average distance of the corresponding points is less than a certain threshold, it can be considered that the points of the two groups constitute a pair of obstacle boundaries and missing boundaries, and the group far away from the instrument is the missing boundary; for the boundary points generated by the resolution relationship, they have obvious linear characteristics and are difficult to form a closed loop, and most of them can be solved by splicing the clouds of multiple measuring stations, so this type of missing can be ignored during detection.

[0194] The filling module is used to fill in a number of points in the missing parts based on spherical projection; fill in a number of points in the missing parts based on spherical projection; according to the scanning principle of the 3D laser scanner, the acquisition of point cloud is based on the spherical coordinate system, and the angle at which the instrument emits the laser is not affected by the target object, but is only related to the state of the instrument itself; the spherical coordinate system has three parameters, the horizontal angle α, the vertical angle β and the distance S. If S is ignored, the point cloud becomes a sphere with uniform density in the spherical coordinate system. The missing parts on the spherical surface are interpolated according to the scanning resolution to complete the preliminary filling of the missing parts.

[0195] The correction module calculates the curvature of the filling point and iteratively minimizes the curvature to correct the position of the filling point.

[0196] Correction module to establish a local coordinate system:

[0197] In a three-dimensional rectangular coordinate system, some surfaces are difficult to represent. Therefore, it is more convenient to calculate by establishing a local coordinate system and representing the surface as a parametric surface. Assume:

[0198]

[0199] Where x, y, z are all differentiable functions of u and v, the domain (u, v) is a region D in two-dimensional space, (u, v) = (x(u, v), y(u, v), z(u, v) gives a mapping from two-dimensional space D to three-dimensional space. Under this mapping, the range of D constitutes a surface D in three-dimensional space.

[0200] A point P in the point cloud i (x i ,y i ,z i ) to establish a local coordinate system (u, v, w): point P i (x i ,y i ,z i ) is the origin, the w axis and point P iThe direction of the normal vector is consistent, the u-axis and the v-axis are orthogonal to each other in the tangent plane, and form a rectangular coordinate system with the w-axis. The local coordinate system is as shown in the attached figure. Fig.10 shown.

[0201] Move the origin O of the original global coordinate system to P i At point , the transformation matrix T is:

[0202]

[0203] Rotate around the x-axis by an angle of α, the rotation matrix R x for:

[0204]

[0205] in

[0206] Rotate angle β around the Y axis, rotation matrix R y for:

[0207]

[0208] in

[0209] Quadratic parametric surface fitting: The local shape of any surface can be approximated by a quadratic surface,

[0210] S(u,v)=(u,v,w(u,v))

[0211] w(u,v)=au 2 +buv+cv 2 +eu+fu

[0212] When the coefficients e and f of the linear term are 0, the quadratic surface is simplified to a quadratic parabola:

[0213] w(u,v)=au 2 +buv+cv 2

[0214] The quadratic parabola is used to fit the point cloud locally; for any point P in the point cloud i , search for its neighboring point P j ∈N(P i ), P i The coordinate value in the local coordinate system is (u i ,v i ,w i ), its transformation matrix B is:

[0215]

[0216] Point p i The neighboring point pj The coordinates in the local coordinate system (u j ,v j ,w j ) into the formula:

[0217] w(u,v)=au 2 +buv+cv 2

[0218] The overdetermined equations are obtained and solved using the least squares method to obtain the parameters a, b, and c of the best fitting quadratic parabola.

[0219] Curvature calculation: calculate the normal vector, first fundamental form and second fundamental form of the surface;

[0220] The surface is represented by a parametric equation: r(u,v)=(x(u,v),y(u,v),z(u,v)), where u and v are parameters. The normal vector of the surface at the point x(u,v),y(u,v),z(u,v) is:

[0221]

[0222] in and are two tangent vectors of the surface, and the cross product of the two gives the normal vector n;

[0223] The first basic form of the line element is expressed as ds 2 =Edu 2 +2Fdudv+Gdv 2 Indicates that E = S u ·S u , F = S u ·S v , G=S v ·S v , S is the parametric equation of the quadratic parabola;

[0224] The second basic form is expressed as II = Ldu 2 +2Mdudv+Ndv 2 , where L = S uu n, M = S uu n, N = S uu n,

[0225] For point P i , whose coordinates in the local coordinate system are (0,0,0), we can find the quadratic parabola at point P i The first-order partial derivative at is:

[0226] S u | (0,0) =(1,0,wu )=(1,0,2cu+bu)=(1,0,0)

[0227] S u | (0,0) =(1,0,w v )=(1,0,2cv+bu)=(0,1,0)

[0228] The quadratic parabola at point P i The second-order partial derivatives and mixed partial derivatives at are:

[0229] S uu | (0,0) =(0,0,w uu )=(0,0,b)

[0230] S vv | (0,0) =(0,0,w vv )=(1,0,2c)

[0231] According to the formula K=4ac-b 2 , H=a+c, the Gaussian curvature K and the average curvature H are calculated.

[0232] Repair the missing parts of the point cloud: After completing the point filling of the missing parts of the point cloud, several filling points are obtained: P t (Rcosθcosα, Rcosθsinα, Rsinθ), where the horizontal angle α and the vertical angle θ are known, and there is only one unknown quantity for the filling point, that is, the distance R from the point to the instrument; select the filling point closest to the reference point as the seed point, and t Select the R value of the nearest reference point as the initial iteration value R0, iterate 1000 steps in the positive and negative directions with a step size of 0.1 mm, calculate the Gaussian curvature of each step, and select the R value that minimizes the absolute value of the Gaussian curvature. min As the R value of the seed point; after the position of the seed point is determined, the position of the seed point is used as a reference point for the remaining filling points to calculate the curvature, and then the next filling point closest to the reference point is selected as the seed point, until the positions of all filling points are determined to complete the repair of the missing part.

Claims

1. A point cloud integrity improvement method based on Delaunay triangulation and curvature minimization, characterized in that: The steps include: S100, transforming the point cloud in the three-dimensional rectangular coordinate system into a spherical coordinate system, and then projecting it onto a two-dimensional plane coordinate system to form a two-dimensional plane point cloud; S200, performing Delaunay triangulation on the two-dimensional plane point cloud, and then projecting it back to a three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular mesh; S300, calculating the side length, area and unit normal vector of each triangular mesh, and eliminating abnormal triangular meshes according to the side length; S400, searching for boundary points of the point cloud according to the topological relationship of the triangular mesh edges; S500, searching for boundary points of the missing part according to the corresponding relationship of the boundary points; S600, filling a number of points in the missing part based on spherical projection; S700, calculating the curvature of the filling point, and iterating to minimize the curvature to correct the position of the filling point.

2. The point cloud integrity improvement method based on Delaunay triangulation and curvature minimization according to claim 1, characterized in that: The step S100 specifically includes the following steps: according to the principle of the three-dimensional laser scanner, the point cloud can be converted into a spherical coordinate system with the instrument as the origin; referring to the idea of ​​parametric surfaces in classical differential geometry, the horizontal angle α and the vertical angle β of the spherical coordinate system are used as the horizontal axis and the vertical axis of the two-dimensional rectangular coordinate system, so as to project the three-dimensional point cloud onto a two-dimensional plane.

3. The point cloud integrity improvement method based on Delaunay triangulation and curvature minimization according to claim 1, characterized in that: The step S200 specifically includes the following steps: performing Delaunay triangulation on the two-dimensional plane point cloud based on the Bowyer-Watson method, and then projecting it back to a three-dimensional rectangular coordinate system to achieve two-dimensional grid division of the three-dimensional point cloud and form a spatial two-dimensional triangular grid, each of which contains three vertices and the connection relationship between them; specifically, first generating an initial grid from a given point set, and then adding the remaining points to the initial grid in turn, finding a triangle in the original triangular grid whose circumscribed circle contains the newly inserted point, which is called the influencing triangle of the newly inserted point, deleting the common edges of the triangle, connecting the new inserted point with all the vertices of the influencing triangle to generate a new triangle, and optimizing it according to the Delaunay standard until all points are added; after performing Delaunay triangulation on the two-dimensional point cloud, projecting it back to a three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular grid, achieving two-dimensional grid division of the three-dimensional point cloud.

4. The point cloud integrity improvement method based on Delaunay triangulation and curvature minimization according to claim 1, characterized in that: The step S300 specifically includes the following steps: calculating the side length, area and unit normal vector of each triangular mesh by formula, setting the side length threshold to be less than 0.5m, and the ratio of the longest side to the shortest side of the triangular mesh cannot exceed 25, and deleting abnormal meshes; e1=(x2,y2,z2)-(x3,y3,z3) e2=(x1,y1,z1)-(x3,y3,z3) e3=(x2,y2,z2)-(x1,y1,z1) Where e1, e2, e3 are the three side lengths of the triangular mesh, S is the area of ​​the triangular mesh, and N is the unit normal vector of the triangular mesh.

5. The point cloud integrity improvement method based on Delaunay triangulation and curvature minimization according to claim 1, characterized in that: The step S500 specifically includes the following steps: for missing points caused by occlusion, the boundary points are divided into several groups according to the boundary lines they belong to, based on the connection relationship of the edge grid, and corresponding points are searched for in these groups based on the angle difference; if the average distance of the corresponding points is less than a certain threshold, it can be considered that the points of the two groups constitute a pair of obstacle boundaries and missing boundaries, where the group farthest from the instrument is the missing boundary, and the boundary points of the missing part are found.

6. The point cloud integrity improvement method based on Delaunay triangulation and curvature minimization according to claim 1, characterized in that: The step S600 specifically includes the following steps: according to the scanning principle of the three-dimensional laser scanner, the acquisition of the point cloud is based on the spherical coordinate system, and the angle at which the instrument emits the laser is not affected by the target object, but is only related to the state of the instrument itself; the spherical coordinate system has three parameters, the horizontal angle α, the vertical angle β and the distance S. When S is ignored, the point cloud becomes a spherical surface with uniform density in the spherical coordinate system. The missing parts on the spherical surface are interpolated according to the scanning resolution to complete the preliminary filling of the missing parts.

7. The point cloud integrity improvement method based on Delaunay triangulation and curvature minimization according to claim 1, characterized in that: The step S700 specifically includes the following steps: S701. Establish a local coordinate system: Some surfaces are difficult to represent in a three-dimensional rectangular coordinate system, so it is more convenient to calculate by establishing a local coordinate system and representing the surface as a parametric surface. set up: Where x, y, z are all differentiable functions of u and v, the domain (u, v) is a region D in two-dimensional space, (u, v) = (x(u, v), y(u, v), z(u, v) gives a mapping from two-dimensional space D to three-dimensional space. Under this mapping, the range of D constitutes a surface D in three-dimensional space. A point P in the point cloud i (x i ,y i ,z i ) to establish a local coordinate system (u, v, w): point P i (x i ,y i ,z i ) is the origin, the w axis and point P i The directions of the normal vectors are consistent, the u-axis and the v-axis are orthogonal to each other in the tangent plane, and form a rectangular coordinate system with the w-axis; Move the origin O of the original global coordinate system to P i At point , the transformation matrix T is: Rotate around the x-axis by an angle of α, the rotation matrix R x for: in Rotate angle β around the Y axis, rotation matrix R y for: in S702, quadratic parametric surface fitting: the local shape of any surface can be approximated by a quadratic surface, S(u,v)=(u,v,w(u,v)) w(u,v)=au 2 +buv+cv 2 +eu+fu When the coefficients e and f of the linear term are 0, the quadratic surface is simplified to a quadratic parabola: w(u,v)=au 2 +buv+cv 2 The quadratic parabola is used to fit the point cloud locally; for any point P in the point cloud i , search for its neighboring point P j ∈N(P i ), P i The coordinate value in the local coordinate system is (u i ,v i ,w i ), its transformation matrix B is: Point p i The neighboring point p j The coordinates in the local coordinate system (u j ,v j ,w j ) into the formula: w(u,v)=au 2 +buv+cv 2 Obtain an overdetermined set of equations, solve them using the least squares method, and obtain the parameters a, b, c of the best fitting quadratic parabola; S703, curvature calculation: calculating the normal vector, the first basic form and the second basic form of the surface; The surface is represented by a parametric equation: r(u,v)=(x(u,v),y(u,v),z(u,v)), where u and v are parameters. The normal vector of the surface at the point x(u,v),y(u,v),z(u,v) is: in and are two tangent vectors of the surface, and the cross product of the two gives the normal vector n; The first basic form is denoted as ds 2 =Edu 2 +2Fdudv+Gdv 2 , where E = S u ·S u , F = S u ·S v , G=S v ·S v , S is the parametric equation of the quadratic parabola; The second basic form is expressed as II = Ldu 2 +2Mdudv+Ndv 2 , where L = S uu n, M = S uu n, N=S uu ·n, For point P i , whose coordinates in the local coordinate system are (0,0,0), we can find the quadratic parabola at point P i The first-order partial derivative at is: S u | (0,0) =(1,0,w u )=(1,0,2cu+bu)=(1,0,0) S u | (0,0) =(1,0,w v )=(1,0,2cv+this)=(0,1,0) The quadratic parabola at point P i The second-order partial derivatives and mixed partial derivatives at are: S uu | (0,0) =(0,0,w uu )=(0,0,b) S vv | (0,0) =(0,0,in vv )=(1,0,2c) According to the formula K=4ac-b 2 , H=a+c calculates the Gaussian curvature K and the average curvature H; S704, repairing the missing parts of the point cloud: After completing the point filling of the missing parts of the point cloud, a number of filling points are obtained: P t (Rcosθcosα, Rcosθsinα, Rsinθ), where the horizontal angle α and the vertical angle θ are known, and there is only one unknown quantity for the filling point, that is, the distance R from the point to the instrument; select the filling point closest to the reference point as the seed point, and t Select the R value of the nearest reference point as the initial value R0 of the iteration, iterate several steps in the positive and negative directions with a fixed step size, calculate the Gaussian curvature of each step, and select the R value that minimizes the absolute value of the Gaussian curvature. min As the R value of the seed point; after the position of the seed point is determined, the position of the seed point is used as a reference point for the remaining filling points to calculate the curvature, and then the next filling point closest to the reference point is selected as the seed point, until the positions of all filling points are determined to complete the repair of the missing part.

8. A point cloud integrity improvement system based on Delaunay triangulation and curvature minimization, characterized in that: It includes a two-dimensional conversion module, which is used to convert the point cloud in the three-dimensional rectangular coordinate system into the spherical coordinate system, and then project it onto the two-dimensional plane coordinate system to form a two-dimensional plane point cloud; The triangulation module is used to perform Delaunay triangulation on the two-dimensional plane point cloud and then project it back to the three-dimensional rectangular coordinate system to form a spatial two-dimensional triangular mesh; The abnormal elimination module is used to calculate the side length, area and unit normal vector of each triangular mesh, and eliminate abnormal triangular meshes according to the side length; The boundary query module is used to find the boundary points of the point cloud according to the topological relationship of the triangular mesh edge, and find the boundary points of the missing part according to the corresponding relationship of the boundary points; A filling module is used to fill in a number of points in the missing part based on spherical projection; The correction module calculates the curvature of the filling point and iteratively minimizes the curvature to correct the position of the filling point.

9. The point cloud integrity improvement system based on Delaunay triangulation and curvature minimization according to claim 8, characterized in that: The abnormal removal module calculates the side length, area and unit normal vector of each triangular mesh by formula, sets the side length threshold to be less than 0.5m, and the ratio of the longest side to the shortest side of the triangular mesh cannot exceed 25, and deletes the abnormal mesh; e1=(x2,y2,z2)-(x3,y3,z3) e2=(x1,y1,z1)-(x3,y3,z3) e3=(x2,y2,z2)-(x1,y1,z1) Where e1, e2, e3 are the three side lengths of the triangular mesh, S is the area of ​​the triangular mesh, and N is the unit normal vector of the triangular mesh.

10. The point cloud integrity improvement system based on Delaunay triangulation and curvature minimization according to claim 8, characterized in that: A local coordinate system is established in the correction module: Some surfaces are difficult to represent in a three-dimensional rectangular coordinate system, so it is more convenient to calculate by establishing a local coordinate system and representing the surface as a parametric surface. set up: Where x, y, z are all differentiable functions of u and v, the domain (u, v) is a region D in two-dimensional space, (u, v) = (x(u, v), y(u, v), z(u, v) gives a mapping from two-dimensional space D to three-dimensional space. Under this mapping, the range of D constitutes a surface D in three-dimensional space. A point P in the point cloud i (x i ,y i ,z i ) to establish a local coordinate system (u, v, w): point P i (x i ,y i ,z i ) is the origin, the w axis and point P i The directions of the normal vectors are consistent, the u-axis and the v-axis are orthogonal to each other in the tangent plane, and form a rectangular coordinate system with the w-axis; Move the origin O of the original global coordinate system to P i At point , the transformation matrix T is: Rotate around the x-axis by an angle of α, the rotation matrix R x for: in Rotate angle β around the Y axis, rotation matrix R y for: in Quadratic parametric surface fitting: The local shape of any surface can be approximated by a quadratic surface, S(u,v)=(u,v,w(u,v)) w(u,v)=au 2 +buv+cv 2 +eu+fu When the coefficients e and f of the linear term are 0, the quadratic surface is simplified to a quadratic parabola: w(u,v)=au 2 +buv+cv 2 The quadratic parabola is used to fit the point cloud locally; for any point P in the point cloud i , search for its neighboring point P j ∈N(P i ), P i The coordinate value in the local coordinate system is (u i ,v i ,w i ), its transformation matrix B is: Point p i The neighboring point p j The coordinates in the local coordinate system (u j ,v j ,w j ) into the formula: w(u,v)=au 2 +buv+cv 2 Obtain an overdetermined set of equations, solve them using the least squares method, and obtain the parameters a, b, c of the best fitting quadratic parabola; Curvature calculation: calculate the normal vector, first fundamental form and second fundamental form of the surface; The surface is represented by a parametric equation: r(u,v)=(x(u,v),y(u,v),z(u,v)), where u and v are parameters. The normal vector of the surface at the point x(u,v),y(u,v),z(u,v) is: in and are two tangent vectors of the surface, and the cross product of the two gives the normal vector n; The first basic form is denoted as ds 2 =Edu 2 +2Fdudv+Gdv 2 , where E = S u ·S u , F = S u ·S v , G=S v ·S v , S is the parametric equation of the quadratic parabola; The second basic form is expressed as II = Ldu 2 +2Mdudv+Ndv 2 , where L = S uu n, M = S uu n, N=S uu ·n, For point P i , whose coordinates in the local coordinate system are (0,0,0), we can find the quadratic parabola at point P i The first-order partial derivative at is: S u | (0,0) =(1,0,w u )=(1,0,2cu+bu)=(1,0,0) S u | (0,0) =(1,0,w v )=(1,0,2cv+this)=(0,1,0) The quadratic parabola at point P i The second-order partial derivatives and mixed partial derivatives at are: S uu | (0,0) =(0,0,w uu )=(0,0,b) S vv | (0,0) =(0,0,in vv )=(1,0,2c) According to the formula K=4ac-b 2 , H=a+c calculates the Gaussian curvature K and the average curvature H; Repair the missing parts of the point cloud: After completing the point filling of the missing parts of the point cloud, several filling points are obtained: P t (Rcosθcosα, Rcosθsinα, Rsinθ), where the horizontal angle α and the vertical angle θ are known, and there is only one unknown quantity for the filling point, that is, the distance R from the point to the instrument; select the filling point closest to the reference point as the seed point, and t Select the R value of the nearest reference point as the initial value R0 of the iteration, iterate several steps in the positive and negative directions with a fixed step size, calculate the Gaussian curvature of each step, and select the R value that minimizes the absolute value of the Gaussian curvature. min The R value as the seed point; After the position of the seed point is determined, the position of the seed point is used as a reference point for the remaining filling points to calculate the curvature, and then the next filling point closest to the reference point is selected as the seed point until the positions of all filling points are determined to complete the repair of the missing part.

Images