An automatic completion method of underground pipeline three-dimensional point cloud data

Abstract

The present application relates to the technical field of three-dimensional modeling, in particular to a kind of automatic completion method of underground pipeline three-dimensional point cloud data, comprising the following steps: according to underground pipeline three-dimensional point cloud data, the boundary ring vertex coordinates of the missing area inside point cloud are identified, the Euclidean distance sum between boundary ring vertex is calculated, the average curvature numerical value of boundary vertex is counted, the normal vector deviation variance value of the least square fitting plane where boundary is calculated, and the missing boundary geometric feature triplets are generated.The present application constructs local geometry frame sequence and projects adjacent points into the law plane to generate dynamic cross-section profile parameter set, breaks through the limitation that traditional method assumes pipeline as standard cylinder, can capture and restore local deformation or non-circular cross-section features that pipeline occurs along path, the generated mesh edge is stitched with original point cloud, and the integrity of underground pipeline model after reconstruction in geometric form and topological structure is ensured.

Description

Technical Field

[0001] This invention relates to the field of 3D modeling technology, and in particular to an automatic completion method for 3D point cloud data of underground pipelines. Background Technology

[0002] 3D modeling technology uses computer graphics, computational geometry, and data processing algorithms to transform the geometry, spatial location, and topology of real-world objects into digital 3D models.

[0003] Current technologies for processing 3D point cloud data of underground pipelines typically rely on general surface reconstruction algorithms or fitting methods based on simple geometric primitives. When faced with pipeline data exhibiting large-area occlusion or breaks, they often struggle to accurately distinguish between localized sparsity of surface data and actual missing pipeline structures. This can lead to long-distance pipeline breaks being misidentified as closed end faces or simple holes, resulting in incorrect sealing phenomena in the repaired model and making it difficult to reconstruct the continuous route of the pipeline. Therefore, improvements are needed. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and propose an automatic completion method for three-dimensional point cloud data of underground pipelines.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: an automatic completion method for three-dimensional point cloud data of underground pipelines, comprising the following steps: Based on the three-dimensional point cloud data of underground pipelines, the coordinates of the boundary ring vertices of the missing regions inside the point cloud are identified, the sum of the Euclidean distances between the vertices of the boundary rings is calculated, the average curvature value of the boundary vertices is calculated, the variance of the normal vector deviation of the least squares fitting plane where the boundary is located is calculated, the missing boundary geometric feature triplet is generated, and the boundary complexity classification index is calculated based on the missing boundary geometric feature triplet. Local contraction iterative operation is performed on the three-dimensional point cloud data of underground pipelines to calculate the geometric median coordinates of local point cloud clusters, generate a discrete central skeleton point set, perform interpolation operation on the discrete central skeleton point set, apply tangent vector continuity constraint conditions, predict and connect the skeleton points of the fracture area according to the tangent vector direction, and generate a continuous smooth skeleton curve. Calculate the tangent vector, principal normal vector, and secondary normal vector at each step length position on the continuous smooth skeleton curve, generate a local geometric frame sequence, project the coordinates of adjacent skeleton points in the three-dimensional point cloud data of underground pipelines onto the normal plane of the local geometric frame sequence, and generate a dynamic cross-sectional profile parameter set. The boundary complexity classification index is called to perform logical branch judgment, generate an initial mesh on the pipeline surface, and stitch the edge vertices of the initial mesh on the pipeline surface with the three-dimensional point cloud data of the underground pipeline to generate a complete three-dimensional model of the underground pipeline.

[0006] Preferably, the step of obtaining the boundary complexity classification index is as follows: Based on the 3D point cloud data of underground pipelines, spatial neighborhood limitation and sequential traversal are performed on the points around the missing area. A continuous boundary point sequence is extracted according to the connection relationship between points. The spatial distance and orientation consistency of the first and last points are verified. When the closure judgment condition is met, the sequence structure is fixed and the coordinates of the vertices of the missing boundary ring are generated. Based on the coordinates of the vertices of the missing boundary ring, the Euclidean distances between adjacent vertices are calculated sequentially and accumulated to obtain the total Euclidean distance. At the same time, the discrete curvature of each vertex is calculated based on the geometric changes of the local neighborhood of the vertex and the average curvature value is obtained. Plane fitting is performed on all vertices and the fluctuation of the normal vector deviating from the fitting plane is statistically analyzed to obtain the normal vector deviation variance value, forming a triplet of geometric features of the missing boundary. Calculate the boundary complexity classification index based on the missing boundary geometric feature triples.

[0007] Preferably, the steps for obtaining the discrete center skeleton point set are as follows: Based on the 3D point cloud data of underground pipelines, local point cloud clusters are divided according to a fixed neighborhood radius and a point-to-point adjacency list is established. The coordinates of each point in the local point cloud cluster are updated according to the number of iterations and synchronously shrunk to the central region of the local point cloud cluster. For local point cloud clusters whose coordinate displacement is less than the convergence threshold after each iteration, the iteration is stopped, and a stable local point cloud cluster coordinate sequence is generated. Based on the shrinking and stable local point cloud cluster coordinate sequence, all point coordinate components are extracted for each local point cloud cluster and the geometric median coordinates are calculated. The geometric median coordinates of each local point cloud cluster are sorted in spatial connectivity order. Fault segments are marked at positions where the distance between adjacent geometric median coordinates exceeds the fault determination threshold, and a discrete central skeleton point set is generated.

[0008] Preferably, the step of obtaining the continuous smooth skeleton curve is as follows: Based on the discrete central skeleton point set, a cubic polynomial interpolation segment is established for the interval between adjacent skeleton points and the interpolation segment coefficients are solved. The second derivative value at each interpolation node is calculated and converted into a tangent direction constraint. The skeleton point coordinates are extrapolated in the fracture section according to the tangent direction constraint and the endpoints of adjacent interpolation segments are connected to generate a continuous smooth skeleton curve.

[0009] Preferably, the step of obtaining the local geometric frame sequence is as follows: Based on the continuous smooth skeleton curve, the coordinate sequence of skeleton sampling points is obtained by sampling sequentially with a fixed discrete step size. For each skeleton sampling point, the difference direction is calculated and normalized by selecting adjacent sampling points before and after each sampling point to obtain the tangent vector. The principal normal vector is calculated and normalized based on the change of the tangent vector. The binormal vector is obtained by performing a cross product between the tangent vector and the principal normal vector and normalizing it, thus generating a local geometric frame sequence.

[0010] Preferably, the steps for obtaining the dynamic cross-sectional profile parameter set are as follows: Based on the local geometric frame sequence, the coordinates of neighboring points of the underground pipeline three-dimensional point cloud data are selected at each skeleton sampling point location and a set of neighboring points is established. The coordinates of the neighboring points are multiplied by the tangent vector, principal normal vector and subnormal vector respectively to obtain local coordinate components. The tangent vector direction coordinate components are discarded and the principal normal vector direction and subnormal vector direction coordinate components are retained to generate the normal plane projection point coordinate sequence. Based on the coordinate sequence of the projection points on the normal plane, the horizontal and vertical components of the projection point coordinates are extracted for each skeleton sampling point. The residuals from the projection points to the elliptical trajectory are calculated point by point and accumulated to obtain the sum of squared distances. The major axis and minor axis parameters of the ellipse are used as variables to be determined and the parameter values ​​are corrected according to the iterative update rules until the sum of squared distances meets the convergence criteria, thereby generating a dynamic cross-sectional profile parameter set.

[0011] Preferably, the step of obtaining the initial mesh on the pipeline surface is as follows: Read the boundary complexity classification index value and the complexity threshold value, compare their sizes to generate branch identifiers. When the branch identifier is a scanning branch, lock the dynamic cross-sectional contour parameter set and the continuous smooth skeleton curve. When the branch identifier is a triangulation branch, lock the coordinates of the missing boundary ring vertex and generate a repair path branch identifier. When the repair path branch is identified as a scan branch, a cross-section placement sequence is generated according to the discrete step size of the continuous smooth skeleton curve. The cross-section vertices of the dynamic cross-section contour parameter set are sequentially placed to the position of each cross-section placement sequence and a vertex index table is established. Vertices with the same index are matched along adjacent cross-section positions and triangular patches are constructed. When the repair path branch is identified as a triangulation branch, the coordinates of the missing boundary ring vertices are projected onto the fitting plane and a plane coordinate sequence is generated to form the initial mesh of the pipeline surface.

[0012] Preferably, the steps for obtaining the complete three-dimensional model of the underground pipeline are as follows: Extract the initial mesh edge vertices of the pipeline surface and establish an edge loop index sequence. Search for the coordinates of the nearest neighbor points in the underground pipeline 3D point cloud data and establish a corresponding relationship table. Connect the edge vertices and nearest neighbor point coordinates one by one according to the corresponding relationship table and complete the stitched triangular facets to generate a complete underground pipeline 3D model.

[0013] Compared with the prior art, the advantages and positive effects of the present invention are as follows: This invention quantifies the geometric complexity of missing regions by identifying the coordinates of the boundary ring vertices within a point cloud and calculating the sum of Euclidean distances, average curvature, and variance of the normal vector deviation. This generates a boundary complexity classification index, which is used for logical branching to distinguish between simple surface holes and complex structural fractures. Skeleton-based scanning reconstruction is employed for high-complexity fracture regions, while planar triangulation is used for low-complexity regions. This effectively improves computational efficiency while ensuring model repair accuracy. The invention utilizes local contraction iteration and geometric median coordinate calculation to generate a discrete center skeleton point set, suppressing random noise interference on the point cloud surface and extracting points located within pipelines. Based on the robust skeleton points of the core, cubic polynomial interpolation is used and tangent vector continuity constraints are applied. This not only fills in the skeleton at the fracture point, but also ensures the smooth transition of the predicted curve in the tangent direction, avoiding unnatural angles in the pipeline route. A local geometric frame sequence is constructed and the neighboring points are projected into the normal plane to generate a dynamic cross-sectional profile parameter set. This breaks through the limitation of traditional methods that assume the pipeline is a standard cylinder. It can capture and restore the local deformation or non-circular cross-sectional features of the pipeline along the path. The generated mesh edges are stitched with the original point cloud, realizing a smooth transition from incomplete data to a watertight model, and ensuring the integrity of the reconstructed underground pipeline model in terms of geometry and topology. Attached Figure Description

[0014] Figure 1 This is a schematic diagram of the steps of the present invention. Detailed Implementation

[0015] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0016] Please see Figure 1 This invention provides a technical solution: an automatic completion method for three-dimensional point cloud data of underground pipelines, comprising the following steps: Based on the 3D point cloud data of underground pipelines, the coordinates of the boundary ring vertices of the missing regions inside the point cloud are identified, the sum of the Euclidean distances between the vertices of the boundary rings is calculated, the average curvature value of the boundary vertices is calculated, the variance of the normal vector deviation of the least squares fitting plane where the boundary is located is calculated, the missing boundary geometric feature triplet is generated, and the boundary complexity classification index is calculated based on the missing boundary geometric feature triplet. Local contraction iterative calculation is performed on the three-dimensional point cloud data of underground pipelines to calculate the geometric median coordinates of local point cloud clusters, generate a discrete central skeleton point set, perform interpolation calculation on the discrete central skeleton point set, apply tangent vector continuity constraints, predict and connect the skeleton points of the fracture area according to the direction of the tangent vector, and generate a continuous smooth skeleton curve. Calculate the tangent vector, principal normal vector, and secondary normal vector at each step length position on the continuous smooth skeleton curve, generate a local geometric frame sequence, project the coordinates of adjacent skeleton points in the three-dimensional point cloud data of underground pipelines onto the normal plane of the local geometric frame sequence, and generate a dynamic cross-sectional profile parameter set. The boundary complexity classification index is used to make logical branch judgments, generate the initial mesh of the pipeline surface, and stitch the edge vertices of the initial mesh of the pipeline surface with the 3D point cloud data of the underground pipeline to generate a complete 3D model of the underground pipeline.

[0017] The steps to obtain the boundary complexity classification index are as follows: Based on the 3D point cloud data of underground pipelines, spatial neighborhood limitation and sequential traversal are performed on the points around the missing area. A continuous boundary point sequence is extracted according to the connection relationship between points. The spatial distance and orientation consistency of the first and last points are verified. When the closure judgment condition is met, the sequence structure is fixed and the coordinates of the vertices of the missing boundary ring are generated. Based on the coordinates of the vertices of the missing boundary ring, the Euclidean distances between adjacent vertices are calculated sequentially and summed to obtain the total Euclidean distance. At the same time, the discrete curvature of each vertex is calculated based on the geometric changes of the local neighborhood of the vertex and the average curvature value is obtained. Plane fitting is performed on all vertices and the fluctuation of the normal vector deviating from the fitted plane is statistically analyzed to obtain the normal vector deviation variance value, forming a triplet of geometric features of the missing boundary. Based on the missing boundary geometric feature triples, the boundary complexity classification index is calculated using the following formula: ; in, As a classification index for boundary complexity, This is the sum of the Euclidean distances between the vertices of the missing boundary ring. This is the average result of the discrete curvature at the vertices of the missing boundary ring. This represents the variance of the normal vector deviation after fitting a plane to the vertices of the missing boundary loop. The dimensionless scaling parameter acting on the distance-curvature product term. Let be the dimensionless proportionality parameter acting on the squared term of the deviation from the normal vector. Let be the dimensionless scaling parameter acting on the interaction term between distance curvature and normal vector deviation. This represents the output of the hyperbolic tangent function with the variance of the normal vector deviation as the independent variable.

[0018] Specifically, based on the 3D point cloud data of underground pipelines, the search range is defined in the established octree index structure. The neighborhood search radius is set to three times the accuracy of the scanning equipment, for example, 0.03 meters. Using this radius, the scattered point cloud at the edge of the missing hole is spatially confined and sequentially traversed. During the traversal, the eigenvalues ​​of the covariance matrix of each neighborhood point set are calculated, and the eigenvalue ratio is used ( To identify boundary feature points, points with a ratio exceeding a preset difference of 0.3 are marked as candidate boundary points. Based on the Euclidean distance principle, a connection relationship is established between candidate points. A doubly linked list structure is used to store and connect broken edge fragments, thereby extracting a continuous sequence of boundary points with their beginning and end connected. Then, the closure of the first and last points of the sequence is checked. The Euclidean distance between the first and last points in the three-dimensional Cartesian coordinate system is calculated and compared with the average density spacing of the point cloud. If the distance is less than 1.5 times the average spacing (e.g., less than 0.015 meters), the dot product of the tangent vectors at the first and last points is calculated. If the dot product value is greater than 0.98, it indicates that the tangent direction has high collinearity, which satisfies the closure judgment condition. At this time, the linked list structure is locked to fix the sequence, and the sorted ordered point set is output to generate the coordinates of the vertices of the missing boundary loop. Based on the coordinates of the missing boundary ring vertices, the quantization extraction process of geometric feature parameters is initiated. Data for each vertex on the boundary ring is read sequentially. The Euclidean distance between adjacent vertices is calculated segment by segment using the distance formula between two points. The lengths of all calculated line segments are stored in an accumulator register for summation until the entire closed loop is traversed, thus obtaining the total Euclidean distance. Simultaneously, with each boundary vertex as the center and a radius of 0.05 meters, local neighborhood points are searched. The least squares method is used to fit a local osculating circle, and the reciprocal of the osculating circle radius is calculated as the discrete curvature value of that vertex, denoted by . The curvature values ​​of all vertices on the ring are summed and divided by the total number of vertices to form the average curvature value. Then, the coordinates of all vertices on the boundary ring are selected to construct a data matrix. The eigenvectors of the data matrix are calculated by principal component analysis (PCA). The vector corresponding to the smallest eigenvalue is selected as the normal vector of the best-fit plane. The dot product of the local normal vector of each boundary vertex and the normal vector of the fitting plane is calculated one by one. The dispersion of all dot product results relative to the mean is statistically analyzed. The fluctuation of the normal vector deviating from the fitting plane is calculated by the variance formula, that is, the dimensionless normal vector deviation variance value is obtained. The sum of distances, average curvature, and variance values ​​obtained above are packaged and stored to form a missing boundary geometric feature triplet. The boundary complexity classification index calculation formula achieves precise classification of the difficulty of repairing pipeline missing areas by nonlinearly coupling macroscopic scale features with microscopic geometric features. The first term... The dimensionless "total spin" characteristic is constructed using the product of length and curvature, reflecting the overall spatial curvature of the missing region; the second term The squared term of variance is used to amplify the sensitivity of surface unevenness, enabling the identification of minute bumps and undulations; the third term... By introducing the hyperbolic tangent function as a saturation activation factor, an interactive response mechanism between scale and roughness is established to prevent exponential explosion when processing large-scale and extremely rough damaged data, thus ensuring the robustness of the evaluation index under different pipe diameters.

[0019] The sum of the Euclidean distances between the vertices of the missing boundary ring represents the physical length of the perimeter of the missing region, in meters (m). This parameter is used to characterize the physical size of the hole to be repaired. The value is derived from the cumulative distance of the boundary ring vertex sequence in the previous steps. The specific steps are: read the ordered list of boundary vertices stored in memory. Iterate through the list to calculate the number of adjacent points. and The L2 norm distance between them. For example, when inspecting a damaged section of a 500mm diameter pipe, 85 boundary points were extracted. After the program accumulated and calculated, the boundary perimeter was found to be 1.85 meters, i.e. .

[0020] This is the average result of the discrete curvature of the vertices of the missing boundary ring, representing the average degree of curvature of the missing edge, in cubic meters (m). This parameter describes the complexity of the geometry; a higher curvature indicates more severe surface bending. The values ​​are derived from statistics on the local geometric properties of the boundary points. The specific steps for obtaining these values ​​are as follows: For each boundary point... Get its neighborhood At each point, construct the covariance matrix and solve for the eigenvalues. Then, estimate the local curvature using the surface variability formula. Finally, the average value is calculated. In the above example, the mean curvature at each point is calculated, and the result is 2.4, that is... Note that... (m) and ( After multiplying, the units cancel each other out, and the result is a dimensionless value.

[0021] The variance of the normal vector deviation after fitting the plane to the vertices of the missing boundary ring represents the undulation dispersion of the missing region relative to the ideal plane, and is a dimensionless value. The numerical values ​​are obtained by analyzing the consistency of the normal vectors. The specific steps are as follows: First, calculate the normal vector of the best-fit plane of the boundary loop. Then calculate each vertex Local normal vector Calculate the dot product of the two. (Results are between 0 and 1), calculate the variance of all dot product results. In the example, due to slight metal tearing and flanging at the damaged edge, the normal vector directions are inconsistent, resulting in a calculated variance of 0.08. .

[0022] This is a dimensionless scaling parameter acting on the distance-curvature product term, used to adjust the contribution of the "total spin" feature to the overall complexity index. This parameter is set based on regression analysis of a historical repair database. The specific steps for obtaining this parameter are: collecting 1000 historical pipeline repair cases and extracting the values ​​from each case... The value was used as the independent variable, and the corresponding normalized value of grid generation time was used as the dependent variable. Linear regression analysis was performed, and the regression coefficients were extracted and normalized. Based on fitting with a large amount of experimental data, this parameter was set to 0.45, i.e. .

[0023] This is a dimensionless scaling parameter acting on the squared term of the normal vector deviation, used to adjust the weight of surface roughness in the complexity assessment. This parameter is obtained through an expert system scoring method. The specific steps are as follows: A sample set containing different roughness levels is established; three senior reverse engineers are invited to score the repair difficulty of the samples (0-10 points); and the scoring results are compared with the sample's... Correlation analysis was performed on the values, and the Pearson correlation coefficient was calculated as the basis for weight setting. After analysis, this parameter was set to 0.35, i.e. .

[0024] This is a dimensionless scaling parameter acting on the interaction term between distance curvature and normal vector deviation, used to balance the bias caused by the extreme nature of a single feature. This parameter is obtained through a grid search method. The specific steps are: setting a search interval with a step size of 0.05, constructing a loss function, and iterating through parameter combinations on the validation set to find the parameter value corresponding to the minimum loss. After optimization, this parameter is set to 0.20, i.e. .

[0025] This represents the output of the hyperbolic tangent function with the variance of the normal vector deviation as the independent variable. This function is used to convert dimensionless variance values... Mapped to The interval serves to smooth and normalize the numerical values. The specific execution process is as follows: The values ​​obtained from the aforementioned calculations... The value is directly substituted into the hyperbolic tangent function formula for calculation. The formula is as follows: For example, when hour, .

[0026] Calculations based on parameters: Substitute the above parameters into the formula: , , , , , .

[0027] The calculation process is as follows: The first step is to calculate the product of distance and curvature (dimension cancellation process: ): ; The second step is to calculate the squared term of the normal vector deviation: ; The third step is to calculate the interaction items: First, calculate the hyperbolic tangent: ; Then calculate the complete interaction item: ; Step 4: Calculate the overall index. ; Calculated boundary complexity classification index The value is 2.071. This numerical result indicates that the area with missing underground pipelines is at a medium-to-high complexity level (the preset threshold is 2.0). A value exceeding the threshold means that the area not only has a large spatial span and curvature (due to...). The region is dominated by a certain degree of surface undulation. Based on this result, the system will automatically determine that this region is not suitable for a simple planar triangulation algorithm. Instead, it should trigger the "scan branch" logic, calling a scan line algorithm based on dynamic cross-sectional contours for mesh repair. This ensures the geometric accuracy and surface smoothness of the final generated model, avoiding creases or holes at bends.

[0028] The steps for obtaining the discrete central skeleton point set are as follows: Based on the 3D point cloud data of underground pipelines, local point cloud clusters are divided according to a fixed neighborhood radius and a point-to-point adjacency list is established. The coordinates of each point in the local point cloud cluster are updated according to the number of iterations and synchronously shrunk to the central region of the local point cloud cluster. For local point cloud clusters whose coordinate displacement is less than the convergence threshold after each iteration, the iteration is stopped, and a stable local point cloud cluster coordinate sequence is generated. Based on the coordinate sequence of the shrinking and stable local point cloud clusters, all point coordinate components are extracted for each local point cloud cluster and the geometric median coordinates are calculated. The geometric median coordinates of each local point cloud cluster are sorted in spatial connectivity order. Fault segments are marked at positions where the distance between adjacent geometric median coordinates exceeds the fault determination threshold, and a discrete central skeleton point set is generated.

[0029] Specifically, based on the 3D point cloud data of underground pipelines, the local neighborhood radius needs to be determined first through statistical analysis. One thousand sampling points are randomly selected from the point cloud data, and the Euclidean distance between each sampling point and its nearest neighbor is calculated. All calculated distance values ​​are summed and divided by the number of sampling points to obtain the average resolution of the point cloud. Set a fixed neighborhood radius Five to eight times the average resolution, for example, calculated Meters, then set Meters, using the KD-Tree data structure to construct a spatial index for the point cloud, with each point as the center, within a radius Within a spherical region, nearest neighbors are searched, and their indices are stored in a point-to-point adjacency list. Then, an iterative shrinking process is initiated, with a maximum of 30 iterations. In each iteration, each point in the point cloud is processed. Calculate its new coordinate position The calculation formula is: ,in, For the updated point coordinates, For point The set of neighboring points in the adjacency list For the first in the neighborhood The coordinates of the neighboring points, This represents the Euclidean distance between two points. This is a Gaussian kernel function used to assign weights based on distance; the closer the points, the greater the weight. After updating the coordinates of all points in one round, it calculates the displacement of each point before and after the update. Set a convergence threshold The threshold is 0.001 meters. This threshold is set based on the accuracy requirements of pipeline reconstruction, meaning the maximum allowable positional error is 1 millimeter, and the statistical displacement is less than the convergence threshold. The number of point cloud clusters, if the displacement of all points within a certain local point cloud cluster is less than If the cluster stops iterating, then the next round of updates continues until the convergence condition is met or the maximum number of iterations is reached. Finally, all converged point coordinates are reorganized according to the original index to generate a shrinking and stable local point cloud cluster coordinate sequence. Based on the coordinate sequence of the shrunk and stable local point cloud clusters, traverse each shrunk local point cloud cluster and extract the three-dimensional coordinate components of all points within the cluster. The mean value of each component is calculated using the averaging method and used as the geometric median coordinate of the cluster. The calculation formula is: ,in, Indicates the first The geometric median coordinates of a local point cloud cluster. This represents the number of points within the cluster. For the first in the cluster After obtaining the geometric median of all clusters by analyzing the coordinate components of each point, the point with the smallest Z-axis coordinate is selected as the starting seed point. A nearest neighbor search strategy is then used to find the unvisited median point closest to the current seed point as the next connection node. This process is iterated until all median points are traversed, forming an initial spatial connectivity order. Subsequently, statistical analysis is performed on the distances between adjacent geometric median points in the connected sequence to establish a distance set. Calculate the average value of the set. and standard deviation Set a fracture detection threshold The setup process is as follows: [Take...] As a baseline length, considering that occlusion during pipeline scanning can lead to data loss, a certain range of jumps is allowed, therefore, it is set to... For example, calculated rice, Rice, then Meters, checking the distance between adjacent geometric median coordinates one by one. ,like If the data is broken at a certain location, the location is marked as the start and end boundary of the broken segment. The connected sequence is split at the break point, noise points and isolated short segments are removed, the point set on the main path is retained, and a discrete center skeleton point set is generated.

[0030] The steps for obtaining a continuous smooth skeleton curve are as follows: Based on the discrete central skeleton point set, a cubic polynomial interpolation segment is established for the interval between adjacent skeleton points and the interpolation segment coefficients are solved. The second derivative value at each interpolation node is calculated and converted into a tangent direction constraint. The skeleton point coordinates are extrapolated in the fracture section according to the tangent direction constraint and the endpoints of adjacent interpolation segments are connected to generate a continuous smooth skeleton curve.

[0031] Specifically, based on the discrete center skeleton point set, the skeleton point set is divided into several continuous intervals, and for each interval, adjacent skeleton points... and Construct a cubic polynomial interpolation function between them The function expression is: ,in, The normalization parameter has a range of values. , Let be the coefficient vector of the interpolation segment to be determined. Indicates the corresponding parameters on the interpolation curve To determine the coordinates of the points, endpoint coordinate constraints are used to solve for the coefficients. and and the tangent vector constraint at the endpoints and ,in and Given the tangent vector at the skeleton point, the second derivative value at each interpolation node is calculated using the central difference method. The continuity condition of the second derivative is used to smoothly adjust the tangent direction. Specifically, the tangent direction is calculated at the previous point. Current point and the last point The rate of change of the constructed vector is used to obtain the initial tangent vector. Then, the tangent direction is optimized by minimizing the objective function using the second derivative. The optimized tangent direction is then substituted into the polynomial equation system as a constraint to solve for the coefficients. to In the identified fracture section, the tangent direction at the fracture point is used... Perform extrapolation prediction and set the extrapolation step size. Calculate the coordinates of the predicted points based on the average spacing of the skeleton points. ,in This is the last skeletal point before fracture. Let be the unit tangent vector at that point. Using a step size scalar, the extrapolation continues until the distance to the endpoint of the next skeleton segment is less than the extrapolation step size. The extrapolated points are then smoothly connected to the original skeleton points. Finally, all interpolated segments and the repaired broken segments are merged to generate a continuous and smooth skeleton curve.

[0032] The steps for obtaining the local geometric frame sequence are as follows: Based on the continuous smooth skeleton curve, the coordinate sequence of skeleton sampling points is obtained by sampling sequentially according to a fixed discrete step size. For each skeleton sampling point, the difference direction is calculated and normalized by selecting the adjacent sampling points before and after it to obtain the tangent vector. The principal normal vector is calculated and normalized based on the change of the tangent vector. The binormal vector is obtained by performing a cross product between the tangent vector and the principal normal vector and normalizing it. A local geometric frame sequence is generated.

[0033] Specifically, based on the continuous smooth skeleton curve, the discrete sampling step size parameter is set. This parameter is set based on the average point spacing of the three-dimensional point cloud data of underground pipelines. The setting principle is Should be greater than To eliminate minor noise interference, while ensuring complete capture of geometric features by using a radius less than one-tenth of the pipeline's minimum bending radius, such as measured... Rice, setting Meters, sampling is performed at equal intervals along the arc length of the curve at this step size to obtain samples containing... The coordinate sequence of skeleton sampling points of discrete points ,in The total number of sampling points, for each skeleton sampling point in the sequence. Select its preceding sampling points With subsequent sampling points The tangent direction vector at this point is calculated using the central difference algorithm. The calculation formula is as follows: ,in, For the first Unit tangent vector at each skeleton sampling point This is the three-dimensional coordinate vector of the next adjacent sampling point. The three-dimensional coordinate vector of the previous adjacent sampling point. The Euclidean norm operation is used to normalize vectors. Next, the rate of change of the tangent vector is calculated to determine the direction of the principal normal. The formula is: ,in, For the first The unit principal normal vector at each point Let be the tangent vector at the next point. Let be the tangent vector at the previous point, indicating the inward direction of the curve's curvature. Using the right-hand rule, the cross product is performed to calculate the binormal vector. The formula is: ,in, For the first The unit binormal vector at each point Let be the tangent vector at the current point. Let be the principal normal vector of the current point. The cross product operator represents vectors. The calculation result needs to be normalized and verified to ensure that the three vectors are pairwise orthogonal and have a magnitude of 1. The three calculated vectors are combined to construct the Frenet frame. The calculation is completed by traversing all sampling points to generate a local geometric frame sequence.

[0034] The steps for obtaining the dynamic cross-sectional profile parameter set are as follows: Based on the local geometric frame sequence, the coordinates of neighboring points of the underground pipeline 3D point cloud data are selected at each skeleton sampling point location and a set of neighboring points is established. The coordinates of the neighboring points are multiplied by the tangent vector, principal normal vector and subnormal vector respectively to obtain the local coordinate components. The tangent vector direction coordinate components are discarded and the principal normal vector direction and subnormal vector direction coordinate components are retained to generate the normal plane projection point coordinate sequence. Based on the coordinate sequence of the projection points on the normal plane, the horizontal and vertical components of the projection point coordinates are extracted for each skeleton sampling point. The residuals from the projection points to the elliptical trajectory are calculated point by point and accumulated to obtain the sum of squared distances. The major axis and minor axis parameters of the ellipse are used as variables to be determined and the parameter values ​​are corrected according to the iterative update rules until the sum of squared distances meets the convergence criteria, thus generating a dynamic cross-sectional profile parameter set.

[0035] Specifically, the local neighborhood search radius is set based on the local geometric frame sequence. This radius value is set with reference to the maximum theoretical radius of the pipeline design. And add redundancy margin, setting the formula as follows: ,in For the search radius, This is the maximum radius value for the pipeline standard. This is the redundancy factor, with a value ranging from 1.2 to 1.5. For example, if the maximum radius of the pipeline is known to be 0.5 meters, then... Then the calculation yields meters, for each skeleton sampling point With the center of the ball, Using a radius, a spatial index is constructed in the 3D point cloud data of underground pipelines, and a range query is performed to extract all point cloud data points falling within the range of the sphere, thus establishing a set of neighboring points. ,in Given the total number of points in the set, iterate through each neighboring point in the set. Construct relative position vectors from skeleton points to neighboring points. ,in For the first The nearest neighbor point relative to the first The displacement vector of each skeleton point The coordinates of the nearest points Let be the coordinates of the skeleton point. Project this relative position vector onto the three axes of the local frame, and decompose the coordinate components using the dot product operation. The calculation formula is as follows: ,in, Let be the coordinate components of the projection point along the principal normal direction. Let be the coordinate components of the projection point in the direction of the binormal. Let be the coordinate components of the projection point along the tangent direction. The first For each skeleton point, the unit principal normal, binormal, and tangent vectors are defined, and a tangential offset threshold is set. It is half the sampling step size, that is Rice, only retain The points are chosen to ensure that the projection points are located near the actual pipeline cross-section, and the tangential direction coordinate components are discarded. Only the principal normal component is retained. component of the binormal direction As two-dimensional plane coordinates, all the retained two-dimensional coordinates are organized according to the skeleton point index to generate a sequence of normal plane projection point coordinates; Based on the coordinate sequence of the projection points on the normal plane, the set of two-dimensional projection points corresponding to each skeleton sampling point is... Perform ellipse fitting and initialize the major axis parameters of the ellipse. With minor axis parameters The initial value is set as the average distance from the centroid of the point set to all points, and the objective function of the distance residual is constructed. The calculation formula is: ,in, Let be the sum of squared distances from all projected points to the fitted ellipse. This represents the total number of projection points involved in the fitting process. The index variable of the point. For the first The x-coordinates of the projection points For the first The ordinates of the projection points The major axis parameter of the ellipse to be optimized is... The minor axis parameter of the ellipse to be optimized. For the first The angles of each point in the polar coordinate system, and the denominator term. Indicates the angle The polar radius length corresponding to the elliptical trajectory is used to correct the parameters according to the gradient descent method and the iteration rules. The updated gradient of the major axis parameter is calculated using the following formula: ,in, For the updated major axis parameter values, This represents the major axis parameter value for the current iteration step. The learning rate is set to 0.01. For the objective function with respect to the major axis parameter Similarly, update the minor axis parameters using the partial derivatives. After each iteration, calculate the change in the objective function value. ,in, This represents the absolute value of the difference between the function values ​​in the two iterations. For the current number The numerical value of the objective function obtained from the next iteration. For the previous round The objective function value obtained from the next iteration is used as the convergence criterion, which is that the change is less than the convergence threshold. threshold The method for obtaining this information is based on the root mean square error requirement of the pipeline cross-section fitting. (For example, 0.001 meters), the formula for calculating the threshold is: ,in For points, Using the preset root mean square error standard, if 100 points participate in the fitting in a certain iteration, then ,when Stop iteration when the time is right and output the final optimized result. and Values ​​are used to generate a dynamic cross-sectional profile parameter set.

[0036] The steps for obtaining the initial mesh on the pipeline surface are as follows: Read the boundary complexity classification index value and the complexity threshold value, compare their values ​​to generate branch identifiers. When scanning a branch, lock the dynamic cross-sectional contour parameter set and the continuous smooth skeleton curve. When triangulating a branch, lock the coordinates of the missing boundary ring vertex and generate a repair path branch identifier. When the repair path branch is identified as a scan branch, a cross-section placement sequence is generated according to the discrete step size of the continuous smooth skeleton curve. The cross-section vertices of the dynamic cross-section profile parameter set are placed sequentially to the position of each cross-section placement sequence and a vertex index table is established. Vertices with the same index are matched along adjacent cross-section positions and triangular patches are constructed. When the repair path branch is identified as a triangulation branch, the coordinates of the missing boundary ring vertices are projected onto the fitting plane and a plane coordinate sequence is generated to form the initial mesh of the pipeline surface.

[0037] Specifically, the boundary complexity classification index value and the complexity threshold value are read. First, a set of boundary complexity classification indices calculated in the previous steps are extracted. Simultaneously, it reads the complexity determination threshold from the preset configuration parameters. The threshold The settings were not arbitrarily chosen, but rather derived from statistical analysis of a large number of historical pipeline repair cases. The specific steps involved collecting data from two thousand samples of different types of underground pipeline damage. Based on the difficulty of manual repair, the samples were categorized into "simple planar type" and "complex curved surface type," and the parameters for each sample were calculated. The values ​​are then plotted and a frequency distribution histogram is generated. A Bayesian classifier is used to calculate the intersection of the two class distribution curves, or to calculate the numerical point that minimizes the classification error rate. The calculation formula is as follows: ,in The segmentation threshold variable is used for traversal. The probability of misclassification at this threshold, for example, calculated for "simple" samples. The values ​​are concentrated in the range of 0.5 to 1.8, indicating "complex" samples. The values ​​are concentrated in the range of 2.2 to 4.5, and the optimal split point between the two is calculated to be 2.0, i.e., set... Next, establish the numerical comparator logic to compare the current input... Value and Perform a floating-point comparison, if ,For example If the missing region is determined to have high curvature or non-planar features, a more adaptive scanline algorithm is required. In this case, the system state is latched, and the previously generated dynamic cross-sectional profile parameter set is stored. The sequence and continuous smooth skeleton curve data are loaded into the cache, and the logical branch identifier variable is set. A value of 1 represents a "scan branch", and vice versa. ,For example If the missing region is approximated as a plane, a more computationally efficient planar partitioning algorithm is suitable. In this case, the original vertex coordinate list of the missing boundary loop is locked, unnecessary skeleton and section parameters are released to save memory, and the logical branch identifier variable is set. The value is set to 0, which represents "triangulation branch". The final output of this state variable generates the repair path branch identifier. When a repair path branch is identified as a scan branch, a cross-section placement sequence is generated based on the discrete step size of the continuous smooth skeleton curve. When a branch identified as a scan branch is detected, the resampling point set of the continuous smooth skeleton curve is read. For each skeleton point Read the corresponding local geometric frame and dynamic cross-sectional profile parameters Set the discrete angle step size of the cross section The degree is 10 degrees, meaning each cross-section generates 36 vertices. The position of each vertex in the local coordinate system is calculated, and then mapped to the world coordinate system using a coordinate transformation matrix. The calculation formula is as follows: ,in, For the first On the first cross section The three-dimensional world coordinates of each vertex. The coordinates of the current skeleton point. Principal normal vector, It is the binormal vector. The vertex index is... and To determine the major and minor axis parameters of the ellipse at that location, after generating all vertices, a global vertex index table is created. Following the topological rules of "top left - bottom left - bottom right" and "top left - bottom right - top right," adjacent cross-sections are then indexed. and The corresponding vertices are connected to construct a triangular patch index. If a branch is detected as a triangulation branch, the coordinate sequence of the missing boundary loop vertices is read, the least-squares fitting plane of the sequence is calculated, and the plane normal vector is obtained. and center point A projection matrix is ​​constructed to project all three-dimensional boundary points onto the two-dimensional plane, resulting in a two-dimensional coordinate set. The constrained Delaunay triangulation algorithm is invoked, using the boundary edges as forced constraint segments to generate a non-overlapping triangular mesh topology in the 2D plane. Then, the generated 2D mesh nodes are mapped back to 3D space using inverse projection transformation. The calculation formula is as follows: ,in and The basis vectors of the local planar coordinate system are used to integrate the generated geometric patch data to form the initial mesh of the pipeline surface.

[0038] The steps to obtain a complete 3D model of underground pipelines are as follows: Extract the initial mesh edge vertices on the pipeline surface and establish an edge loop index sequence. Search for the coordinates of the nearest neighbor points in the underground pipeline 3D point cloud data and establish a correspondence table. Connect the edge vertices and nearest neighbor point coordinates one by one according to the correspondence table and complete the stitched triangular facets to generate a complete underground pipeline 3D model.

[0039] Specifically, the initial mesh edge vertices on the pipeline surface are extracted and an edge loop index sequence is established. The edge structure of all triangular facets in the generated initial mesh is traversed, and the reference count of each edge is counted. Edges with a reference count of only 1 are identified as open edges. These edges are sorted according to their geometric connectivity, and an ordered edge loop index sequence is constructed. Set the distance threshold for nearest neighbor search. The threshold is set to three times the average spacing of the point cloud, for example, 0.06 meters. The KD-Tree algorithm is used to search for and match points in the original underground pipeline 3D point cloud data. Each vertex The closest boundary point in space Establish a point-to-point mapping table For each pair of matching points in the mapping table, the consistency of their normal vectors is checked. If the angle between the two normal vectors is less than 45 degrees, it is considered a valid stitching pair. For valid point pairs, a "zipper" stitching strategy is adopted, that is, stitching at the vertices of the mesh edge. , Boundary points of the point cloud , Constructing transition triangles between them involves calculating the quadrilaterals. Given the length of the diagonal, choose the shorter diagonal to divide the quadrilateral into two triangles. For example, if... Then construct a triangle and This process fills the tiny gaps between the newly generated mesh and the original point cloud, and finally updates the normal vector information of all vertices in a unified manner, merges the vertex buffer and the index buffer, and generates a complete 3D model of the underground pipeline.

[0040] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for automatically completing three-dimensional point cloud data of underground pipelines, characterized in that, Includes the following steps: Based on the three-dimensional point cloud data of underground pipelines, the coordinates of the boundary ring vertices of the missing regions inside the point cloud are identified, the sum of the Euclidean distances between the vertices of the boundary rings is calculated, the average curvature value of the boundary vertices is calculated, the variance of the normal vector deviation of the least squares fitting plane where the boundary is located is calculated, the missing boundary geometric feature triplet is generated, and the boundary complexity classification index is calculated based on the missing boundary geometric feature triplet. Local contraction iterative operation is performed on the three-dimensional point cloud data of underground pipelines to calculate the geometric median coordinates of local point cloud clusters, generate a discrete central skeleton point set, perform interpolation operation on the discrete central skeleton point set, apply tangent vector continuity constraint conditions, predict and connect the skeleton points of the fracture area according to the tangent vector direction, and generate a continuous smooth skeleton curve. Calculate the tangent vector, principal normal vector, and secondary normal vector at each step length position on the continuous smooth skeleton curve, generate a local geometric frame sequence, project the coordinates of adjacent skeleton points in the three-dimensional point cloud data of underground pipelines onto the normal plane of the local geometric frame sequence, and generate a dynamic cross-sectional profile parameter set. The boundary complexity classification index is called to perform logical branch judgment, generate an initial mesh on the pipeline surface, and stitch the edge vertices of the initial mesh on the pipeline surface with the three-dimensional point cloud data of the underground pipeline to generate a complete three-dimensional model of the underground pipeline.

2. The method for automatically completing three-dimensional point cloud data of underground pipelines according to claim 1, characterized in that, The steps for obtaining the boundary complexity classification index are as follows: Based on the 3D point cloud data of underground pipelines, spatial neighborhood limitation and sequential traversal are performed on the points around the missing area. A continuous boundary point sequence is extracted according to the connection relationship between points. The spatial distance and orientation consistency of the first and last points are verified. When the closure judgment condition is met, the sequence structure is fixed and the coordinates of the vertices of the missing boundary ring are generated. Based on the coordinates of the vertices of the missing boundary ring, the Euclidean distances between adjacent vertices are calculated sequentially and accumulated to obtain the total Euclidean distance. At the same time, the discrete curvature of each vertex is calculated based on the geometric changes of the local neighborhood of the vertex and the average curvature value is obtained. Plane fitting is performed on all vertices and the fluctuation of the normal vector deviating from the fitting plane is statistically analyzed to obtain the normal vector deviation variance value, forming a triplet of geometric features of the missing boundary. Calculate the boundary complexity classification index based on the missing boundary geometric feature triples.

3. The method for automatically completing three-dimensional point cloud data of underground pipelines according to claim 1, characterized in that, The steps for obtaining the discrete central skeleton point set are as follows: Based on the 3D point cloud data of underground pipelines, local point cloud clusters are divided according to a fixed neighborhood radius and a point-to-point adjacency list is established. The coordinates of each point in the local point cloud cluster are updated according to the number of iterations and synchronously shrunk to the central region of the local point cloud cluster. For local point cloud clusters whose coordinate displacement is less than the convergence threshold after each iteration, the iteration is stopped, and a stable local point cloud cluster coordinate sequence is generated. Based on the shrinking and stable local point cloud cluster coordinate sequence, all point coordinate components are extracted for each local point cloud cluster and the geometric median coordinates are calculated. The geometric median coordinates of each local point cloud cluster are sorted in spatial connectivity order. Fault segments are marked at positions where the distance between adjacent geometric median coordinates exceeds the fault determination threshold, and a discrete central skeleton point set is generated.

4. The method for automatically completing three-dimensional point cloud data of underground pipelines according to claim 1, characterized in that, The steps for obtaining the continuous smooth skeleton curve are as follows: Based on the discrete central skeleton point set, a cubic polynomial interpolation segment is established for the interval between adjacent skeleton points and the interpolation segment coefficients are solved. The second derivative value at each interpolation node is calculated and converted into a tangent direction constraint. The skeleton point coordinates are extrapolated in the fracture section according to the tangent direction constraint and the endpoints of adjacent interpolation segments are connected to generate a continuous smooth skeleton curve.

5. The method for automatically completing three-dimensional point cloud data of underground pipelines according to claim 1, characterized in that, The steps for obtaining the local geometric frame sequence are as follows: Based on the continuous smooth skeleton curve, the coordinate sequence of skeleton sampling points is obtained by sampling sequentially with a fixed discrete step size. For each skeleton sampling point, the difference direction is calculated and normalized by selecting adjacent sampling points before and after each sampling point to obtain the tangent vector. The principal normal vector is calculated and normalized based on the change of the tangent vector. The binormal vector is obtained by performing a cross product between the tangent vector and the principal normal vector and normalizing it, thus generating a local geometric frame sequence.

6. The method for automatically completing three-dimensional point cloud data of underground pipelines according to claim 1, characterized in that, The steps for obtaining the dynamic cross-sectional profile parameter set are as follows: Based on the local geometric frame sequence, the coordinates of neighboring points of the underground pipeline three-dimensional point cloud data are selected at each skeleton sampling point location and a set of neighboring points is established. The coordinates of the neighboring points are multiplied by the tangent vector, principal normal vector and subnormal vector respectively to obtain local coordinate components. The tangent vector direction coordinate components are discarded and the principal normal vector direction and subnormal vector direction coordinate components are retained to generate the normal plane projection point coordinate sequence. Based on the coordinate sequence of the projection points on the normal plane, the horizontal and vertical components of the projection point coordinates are extracted for each skeleton sampling point. The residuals from the projection points to the elliptical trajectory are calculated point by point and accumulated to obtain the sum of squared distances. The major axis and minor axis parameters of the ellipse are used as variables to be determined and the parameter values ​​are corrected according to the iterative update rules until the sum of squared distances meets the convergence criteria, thereby generating a dynamic cross-sectional profile parameter set.

7. The method for automatically completing three-dimensional point cloud data of underground pipelines according to claim 1, characterized in that, The steps for obtaining the initial mesh on the pipeline surface are as follows: Read the boundary complexity classification index value and the complexity threshold value, compare their sizes to generate branch identifiers. When the branch identifier is a scanning branch, lock the dynamic cross-sectional contour parameter set and the continuous smooth skeleton curve. When the branch identifier is a triangulation branch, lock the coordinates of the missing boundary ring vertex and generate a repair path branch identifier. When the repair path branch is identified as a scan branch, a cross-section placement sequence is generated according to the discrete step size of the continuous smooth skeleton curve. The cross-section vertices of the dynamic cross-section contour parameter set are sequentially placed to the position of each cross-section placement sequence and a vertex index table is established. Vertices with the same index are matched along adjacent cross-section positions and triangular patches are constructed. When the repair path branch is identified as a triangulation branch, the coordinates of the missing boundary ring vertices are projected onto the fitting plane and a plane coordinate sequence is generated to form the initial mesh of the pipeline surface.

8. The method for automatically completing three-dimensional point cloud data of underground pipelines according to claim 1, characterized in that, The steps for obtaining the complete three-dimensional model of the underground pipeline are as follows: Extract the initial mesh edge vertices of the pipeline surface and establish an edge loop index sequence. Search for the coordinates of the nearest neighbor points in the underground pipeline 3D point cloud data and establish a corresponding relationship table. Connect the edge vertices and nearest neighbor point coordinates one by one according to the corresponding relationship table and complete the stitched triangular facets to generate a complete underground pipeline 3D model.

Images