Elevation calculation method and system based on Delaunay triangulation network

Abstract

The invention discloses an elevation calculation method and system based on Delaunay triangulation networks, and the method comprises the steps: obtaining a plurality of elevation point sets to be calculated, and carrying out the triangulation of the elevation point sets, and obtaining a plurality of Delaunay triangulation networks; the Delaunay triangulation network is projected into a two-dimensional plane to obtain a triangle set, triangles in the triangle set are split to obtain a triangle undirected graph, and the undirected graph comprises a vertex set and an edge set; obtaining a to-be-queried point, and searching a first point location closest to the to-be-queried point in the undirected graph by using a KD tree; and based on the to-be-queried point and the first point location, determining a triangle containing the to-be-queried point, and based on a parameter equation of a plane where the triangle is located, obtaining an elevation value corresponding to the to-be-queried point. According to the method, the elevation value in a large-scale scene can be quickly and accurately calculated.

Description

Technical Field

[0001] The present application relates to the field of elevation estimation, and in particular to an elevation calculation method and system based on a Delaunay triangulation. Background Art

[0002] In CAD drawing operations, it is necessary to generate a Delaunay triangulation based on specified points, then determine the triangle where the specified point is located, and then calculate the elevation of the specified point in the triangle. In this process, determining the triangle where the specified point is located is very time-consuming.

[0003] The existing technology generally estimates the elevation of a specified point through a linear search method and a jump and walk method. The linear search method needs to traverse each triangle in the triangulation network and check whether the specified point is inside the triangle. This method is time-consuming and is not suitable for scenes with large triangulation networks. The jump and walk method needs to randomly select a starting triangle and determine whether the specified point is inside the triangle in turn. This method is related to the distribution of Delaunay triangles. If the point set is unevenly distributed or the starting point is not well selected, a large number of triangles need to be visited. Therefore, when processing scenes with large triangulation networks, it takes relatively long time.

[0004] Application Contents

[0005] The present application provides an elevation calculation method based on the Delaunay triangulation to quickly and accurately calculate elevation values ​​in large-scale scenes.

[0006] In a first aspect, the present application provides an elevation calculation method based on a Delaunay triangulation network, comprising:

[0007] Acquire a plurality of elevation point sets to be calculated, and triangulate the elevation point sets to obtain a plurality of Delaunay triangulations;

[0008] Projecting the Delaunay triangulation network onto a two-dimensional plane to obtain a triangle set, and splitting each triangle in the triangle set to obtain a triangle undirected graph, wherein the undirected graph includes a vertex set and an edge set;

[0009] Obtain a point to be queried, and use a KD tree to search for the first point in the undirected graph that is closest to the point to be queried;

[0010] Based on the point to be queried and the first point position, a triangle containing the point to be queried is determined, and based on a parametric equation of the plane where the triangle is located, an elevation value corresponding to the point to be queried is obtained.

[0011] In the embodiment of the present application, by obtaining several elevation point sets to be calculated and triangulating the elevation point sets to obtain several Delaunay triangulations, a Delaunay triangulation can be quickly generated, which is convenient for the subsequent determination of the triangle where the query point is located; by projecting the Delaunay triangulation onto a two-dimensional plane to obtain a triangle set, and splitting each triangle in the triangle set to obtain a triangle undirected graph, the triangulation network in the three-dimensional space can be projected onto a two-dimensional plane, which is convenient for splitting each vertex set and edge set in each triangle, thereby facilitating the subsequent determination of the triangle containing the query point. Triangle and elevation calculation; by obtaining the point to be queried, and using the KD tree to search for the first point position in the undirected graph that is closest to the point to be queried, the first point position that is closest to the point to be queried can be quickly obtained, which is convenient for the subsequent determination of the triangle containing the point to be queried and the elevation calculation; based on the point to be queried and the first point position, the triangle containing the point to be queried is determined, and based on the parametric equation of the plane where the triangle is located, the elevation value corresponding to the point to be queried is obtained, and calculation can be performed based on the determined parametric equation of the triangle containing the point to be queried, so that the elevation value in large-scale scenes can be quickly and accurately calculated.

[0012] Furthermore, the step of obtaining a plurality of elevation point sets to be calculated is as follows:

[0013] Collect several elevation points using RTK or total station;

[0014] Importing a plurality of the elevation points into CAD, a plurality of elevation point sets to be calculated are obtained.

[0015] In this way, several elevation points are collected in advance through RTK or total station, and imported into CAD to obtain several sets of elevation points to be calculated, which facilitates the subsequent construction of Delaunay triangulation.

[0016] Furthermore, the use of the KD tree to search for the first point in the undirected graph that is closest to the query point is specifically as follows:

[0017] Compare the point to be queried with the root node, if the value of the point to be queried in the current dimension is less than the value of the root node and the root node has no subtree, determine that the root node is a second point, and calculate a first distance between the second point and the point to be queried;

[0018] Search and calculate the second distance between the third point in the parent node corresponding to the root node and the point to be queried that is the shortest, and continuously update the second point based on the second distance until there is no smaller distance, and determine that the second point is the first point that is closest to the point to be queried.

[0019] In this way, by using the KD tree to search the first point in the undirected graph that is closest to the point to be queried, the first point that is closest to the point to be queried can be quickly acquired, which facilitates the subsequent determination of the triangle containing the point to be queried and the elevation calculation.

[0020] Further, the determining of a triangle including the point to be queried based on the point to be queried and the first point position is specifically:

[0021] Taking the first point as the starting point, draw a first ray through the point to be queried;

[0022] Extracting a plurality of candidate points adjacent to the first point from the undirected graph, using the first candidate point and the second candidate point as a second ray and a third ray passing through the first point, respectively, and calculating the angles between the first ray and the second ray and the first ray, respectively, to obtain a first angle and a second angle;

[0023] Compare the first angle with the second angle, and if the second angle is smaller than the first angle, determine the second candidate point corresponding to the second angle as the optimal candidate point;

[0024] The optimal candidate point is used as a starting point to continue drawing rays through the point to be queried until a duplicate point appears, and a triangle containing the point to be queried is obtained based on the rays between the duplicate points.

[0025] In this way, by determining a triangle containing the point to be queried, and performing calculations based on the parameter equation of the determined triangle containing the point to be queried, the elevation value in a large-scale scene can be calculated quickly and accurately.

[0026] Furthermore, the parametric equation of the plane where the triangle is located is specifically:

[0027] A(x-x1)+B(y-y1)+C(z-z1)=0;

[0028] In the formula, A is ((y2-y1)(z3-z1)-(z2-z1)(y3-y1)), B is ((x2-x1)(z3-z1)-(z2-z1)(x3-x1)), C is ((x2-x1)(y3-y1)-(y2-y1)(x3-x1)), x1, y1 and z1 are the positions of the first point.

[0029] In this way, by performing calculation based on the determined parametric equation of the triangle containing the point to be queried, the elevation value in a large-scale scene can be calculated quickly and accurately.

[0030] In a second aspect, the present application provides an elevation calculation system based on a Delaunay triangulation network, comprising: a first acquisition module, a first processing module, a second acquisition module, and a second processing module;

[0031] The first acquisition module is used to acquire a plurality of elevation point sets to be calculated, and triangulate the elevation point sets to obtain a plurality of Delaunay triangulations;

[0032] The first processing module is used to project the Delaunay triangulation network into a two-dimensional plane to obtain a triangle set, and split each triangle in the triangle set to obtain a triangle undirected graph, wherein the undirected graph includes a vertex set and an edge set;

[0033] The second acquisition module is used to acquire the point to be queried, and use the KD tree to search the first point in the undirected graph that is closest to the point to be queried;

[0034] The second processing module is used to determine a triangle containing the point to be queried based on the point to be queried and the first point position, and obtain an elevation value corresponding to the point to be queried based on a parametric equation of the plane where the triangle is located.

[0035] In the embodiment of the present application, by obtaining several elevation point sets to be calculated and triangulating the elevation point sets to obtain several Delaunay triangulations, a Delaunay triangulation can be quickly generated, which is convenient for the subsequent determination of the triangle where the query point is located; by projecting the Delaunay triangulation onto a two-dimensional plane to obtain a triangle set, and splitting each triangle in the triangle set to obtain a triangle undirected graph, the triangulation network in the three-dimensional space can be projected onto a two-dimensional plane, which is convenient for splitting each vertex set and edge set in each triangle, thereby facilitating the subsequent determination of the triangle containing the query point. Triangle and elevation calculation; by obtaining the point to be queried, and using the KD tree to search for the first point position in the undirected graph that is closest to the point to be queried, the first point position that is closest to the point to be queried can be quickly obtained, which is convenient for the subsequent determination of the triangle containing the point to be queried and the elevation calculation; based on the point to be queried and the first point position, the triangle containing the point to be queried is determined, and based on the parametric equation of the plane where the triangle is located, the elevation value corresponding to the point to be queried is obtained, and calculation can be performed based on the determined parametric equation of the triangle containing the point to be queried, so that the elevation value in large-scale scenes can be quickly and accurately calculated.

[0036] Furthermore, the first acquisition module includes: a collection unit and an import unit;

[0037] The acquisition unit is used to acquire a number of elevation points using RTK or a total station;

[0038] The importing unit is used to import the plurality of elevation points into CAD to obtain a plurality of elevation point sets to be calculated.

[0039] Furthermore, the second acquisition module includes: a comparison unit and an update unit;

[0040] The comparison unit is configured to compare the to-be-queried point with a root node, and if the value of the to-be-queried point in the current dimension is less than the value of the root node and the root node has no subtree, determine that the root node is a second point position, and calculate a first distance between the second point position and the to-be-queried point;

[0041] The updating unit is used to search for and calculate the second distance between the third point position in the parent node corresponding to the root node and the point to be queried, and continuously update the second point position based on the second distance until there is no smaller distance, and determine that the second point position is the first point position closest to the point to be queried.

[0042] Further, the second processing module includes: a third processing unit, a calculation unit, a second comparison unit and a construction unit;

[0043] The third processing unit is used to make a first ray passing through the to-be-queried point from the first point as a starting point;

[0044] The calculation unit is used to extract a plurality of candidate points adjacent to the first point from the undirected graph, use the first candidate point and the second candidate point as a second ray and a third ray to pass through the first point, and calculate the angles between the first ray and the second ray and the first ray, respectively, to obtain a first angle and a second angle;

[0045] The second comparison unit is used to compare the first angle with the second angle, and if the second angle is smaller than the first angle, determine that the second candidate point corresponding to the second angle is the optimal candidate point;

[0046] The construction unit is used to take the optimal candidate point as a starting point and continue to make rays passing through the point to be queried until a repeated point appears, and obtain a triangle containing the point to be queried based on the rays between the repeated points.

[0047] Furthermore, the parametric equation of the plane where the triangle is located is specifically:

[0048] A(x-x1)+B(y-y1)+C(z-z1)=0;

[0049] In the formula, A is ((y2-y1)(z3-z1)-(z2-z1)(y3-y1)), B is ((x2-x1)(z3-z1)-(z2-z1)(x3-x1)), C is ((x2-x1)(y3-y1)-(y2-y1)(x3-x1)), x1, y1 and z1 are the positions of the first point. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 It is a flowchart of an embodiment of an elevation calculation method based on Delaunay triangulation provided by the present application;

[0051] Figure 2 It is a schematic diagram for determining the optimal candidate point provided by the present application;

[0052] Figure 3 It is a structural schematic diagram of an embodiment of an elevation calculation system based on Delaunay triangulation provided by the present application. DETAILED DESCRIPTION

[0053] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0054] It should be understood that the step numbers used in this article are only for the convenience of description and are not intended to limit the order in which the steps are executed.

[0055] It should be understood that the terms used in the present specification are only for the purpose of describing specific embodiments and are not intended to limit the present invention. As used in the present specification and the appended claims, the singular forms "a", "an" and "the" are intended to include plural forms unless the context clearly indicates otherwise.

[0056] The terms “include” and “comprising” indicate the presence of described features, integers, steps, operations, elements and / or components, but do not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or combinations thereof.

[0057] The term "and / or" means and includes any and all possible combinations of one or more of the associated listed items.

[0058] In CAD drawing, after constructing the Delaunay triangulation, it is necessary to locate the triangle to which the specified point belongs to calculate its elevation, which is a time-consuming process. Existing elevation estimation technologies include linear search and jump and walk methods. The linear search method requires traversing all triangles, which is inefficient and not suitable for large-scale triangulation. The jump and walk method has a random starting point and relies on the Delaunay distribution. When the point set is uneven or the starting point is not well selected, it takes a relatively long time.

[0059] The Delaunay triangulation (also known as the Delaunay triangulation) is a special triangulation network that follows a series of strict geometric criteria to construct a set of adjacent and non-overlapping triangles.

[0060] Embodiment 1

[0061] Please refer to Figure 1 , Figure 1 It is a flowchart of an embodiment of an elevation calculation method based on Delaunay triangulation provided by the present application, including steps S1 to S4;

[0062] Step S1, obtaining a plurality of elevation point sets to be calculated, and triangulating the elevation point sets to obtain a plurality of Delaunay triangulations;

[0063] It can be understood that several elevation points can be collected using RTK or total station; and by importing several of the elevation points into CAD, several sets of elevation points to be calculated can be obtained, where the elevation points are point information in a three-dimensional space, including x-axis coordinates, y-axis coordinates and z-axis coordinates.

[0064] In this way, several elevation points are collected in advance through RTK or total station, and imported into CAD to obtain several sets of elevation points to be calculated, which facilitates the subsequent construction of Delaunay triangulation.

[0065] After obtaining several elevation point sets, it is necessary to preprocess the elevation point sets and select the Delaunay triangulation algorithm to triangulate the preprocessed elevation point sets to obtain several Delaunay triangulation networks; wherein the triangulation algorithm includes but is not limited to the Bowyer-Watson algorithm, the Lawson local optimization algorithm (LOP), etc., for example but not limitation, taking the Bowyer-Watson algorithm for triangulation as an example, first, it is necessary to construct a sufficiently large super triangle as the initial triangular network, the vertices of which include all the points in the elevation point set; secondly, the points in the elevation point set are inserted into the triangular mesh one by one, for each newly inserted point, the influencing triangle containing the point is found, and the common edges of the influencing triangles are deleted, and the new point is connected to the vertices of the influencing triangle to form three new triangles. Again, perform an empty circumscribed circle check on each newly formed triangle to determine whether it meets the Delaunay condition. If the circumscribed circle of a triangle contains other elevation points, it means that the triangle does not meet the Delaunay condition. At this time, it is necessary to optimize the triangle by exchanging diagonals until all triangles meet the empty circumscribed circle condition; then, repeat the above steps of inserting points, constructing triangles, and detecting and optimizing empty circumscribed circles until all elevation points are inserted into the triangular mesh; finally, delete the triangles related to the super triangle to obtain the final Delaunay triangulation.

[0066] It should be noted that data preprocessing methods include but are not limited to removing duplicate points, processing noise points and filling missing data, and these data preprocessing methods can be achieved through existing technologies, but data preprocessing methods are not the focus of this application, so they will not be expanded here.

[0067] It should be noted that the Bowyer-Watson algorithm is an incremental algorithm for generating Delaunay triangulation. Delaunay triangulation is a special triangulation method that generates a triangulated network that meets the Delaunay condition by gradually inserting points and rebuilding the triangular mesh until all points are added. The Lawson Local Optimization Procedure (LOP) is an algorithm used to optimize the triangulated network structure during the triangulation process to ensure that it meets the Delaunay triangulation characteristics. It optimizes triangles that do not meet the Delaunay condition by exchanging diagonals.

[0068] Step S2, projecting the Delaunay triangulation network onto a two-dimensional plane to obtain a triangle set, and splitting each triangle in the triangle set to obtain a triangle undirected graph, wherein the undirected graph includes a vertex set and an edge set;

[0069] It is understandable that when a Delaunay triangulation is obtained, it is necessary to simplify the Delaunay triangulation into a two-dimensional plane for further analysis. Specifically, first, it is necessary to predetermine a suitable projection method to better retain the structural characteristics of the three-dimensional Delaunay triangulation; second, project each point in the Delaunay triangulation into a two-dimensional plane according to the selected projection method, and during the projection process, it is necessary to calculate the coordinates of each point on the two-dimensional plane and retain the relative position relationship between the points; finally, re-triangulate the two-dimensional plane according to the projected point set to obtain a set of triangles in the two-dimensional plane, wherein the triangulation method has been detailed in step S1, so it will not be repeated here.

[0070] When a triangle set on a two-dimensional plane is obtained, it is necessary to extract all unique vertices from the triangle set and collect all vertices to form a vertex set V, V = {v1, v2, ..., v n}, where v i Indicates the i-th vertex. Next, traverse each triangle in the triangle set and add the three edges connected to each triangle to the edge set E. Each edge consists of two vertices, indicating that the two vertices are adjacent in the triangle. For each triangle Δ ijk , whose vertices are v i , v j , v k , its three sides (v i ,v j ),(v i ,v k ),(v j ,v k ) is added to the set of non-phase edges E, that is, (v i ,v j )∈E, if and only if (v j ,v i )∈E. Finally, the vertex set V and edge set E are used to construct an undirected graph G=(V,E).

[0071] It should be noted that the projection method includes but is not limited to orthogonal projection, perspective projection and isometric projection, etc. Projecting the Delaunay triangulation into a two-dimensional plane is not the focus of this application, so it will not be expanded here.

[0072] It should be noted that an undirected graph is a concept in graph theory, in which the edge has no direction, that is, the edge connecting two vertices can be regarded as bidirectional. In an undirected graph, vertices and edges together constitute the basic elements of the graph, and the edge set generally defines the connection relationship between vertices. A triangular undirected graph can be understood as a special undirected graph, which contains multiple triangular structures formed by three vertices connected end to end by three edges.

[0073] Step S3, obtaining a point to be queried, and using a KD tree to search for the first point in the undirected graph that is closest to the point to be queried;

[0074] It can be understood that after obtaining the triangular undirected graph, a KD tree can be established based on the vertex set V to speed up the search efficiency of the nearest point, where the process of constructing the KD tree is as follows: (1) arbitrarily select a vertex from the vertex set V as the root node; (2) In each layer of the KD tree, it is necessary to cyclically select different dimensions as the partitioning dimensions; (3) For the selected partitioning dimension, calculate the median of all vertices on the dimension, and use the median as the root node to divide the remaining points into two left and right subtrees; (4) Repeat steps (2) and (3) for the vertices of the left and right subtrees respectively, and recursively construct the left subtree and the right subtree until the subtree is empty.

[0075] After the KD tree is constructed, the query point p can be obtained. i , and use the KD tree to search the undirected graph G with p i Euclidean distance to the nearest point v s The specific query process is as follows: First, select the point p to be queried i The first dimension is the current dimension d, and the query point is compared with the same dimension d t The root node p c Compare, if the query point p i The value at the current dimension is less than the root node p c When the value of , enter the root node p c Continue searching the left subtree of the node; otherwise, go to the root node p c Continue searching the right subtree of the node, and repeat the search until the root node p c If there is no subtree, then the root node p is determined c is a second point that is closer, and a first distance d between the second point and the point to be queried is calculated. c ; Secondly, search and calculate the root node p c The corresponding parent node p t The third point in the query point has the smallest second distance d a , and based on the second distance d aThe second point position is continuously updated until there is no smaller distance, and the second point position is determined to be the first point position v closest to the query point. s By way of example and not limitation, if the second distance d a <First distance d c , then calculate the parent node p t and the query point p i The distance between t , if d t <d c , then its parent node p t For the node with the closest distance, update the second point, i.e. p c =p t , d c =d t , recursively search for p t The left and right subtrees of are repeatedly recursively searched to determine the point closest to the query point.

[0076] In this way, by using the KD tree to search the first point in the undirected graph that is closest to the point to be queried, the first point that is closest to the point to be queried can be quickly acquired, which facilitates the subsequent determination of the triangle containing the point to be queried and the elevation calculation.

[0077] Step S4: based on the point to be queried and the first point position, determine a triangle containing the point to be queried, and based on a parametric equation of the plane where the triangle is located, obtain an elevation value corresponding to the point to be queried.

[0078] It can be understood that after determining the first point position closest to the query point, the next path node v can be searched. p (that is, the optimal candidate point below) to determine the triangle containing the query point. Specifically, first, take the first point v s Make the first ray l as the starting point s Passing the query point p i ; Secondly, randomly extract the first point v from the undirected graph G s The adjacent vertex set V includes several candidate points. The first candidate point v is randomly selected from the vertex set V. p , and take the first candidate point v p Draw the second ray l as the starting point p After the first point v s , where the initial optimal angle d is set p =2π, traverse the vertex set V and select the second candidate point v k , and take the second candidate point v k Make the third ray l as the starting point k After the first point vs , and calculate the second ray l p and the third ray l k With the first ray l s The first angle d is obtained p and the second angle d k ; Again, the first angle d p and the second angle d k For comparison, if the second angle d k is smaller than the first angle d p , then determine the second angle d k The corresponding second candidate point v k is the optimal candidate point, even if v p =v k , and make d p =d k , otherwise keep v p is the optimal candidate point; finally, the optimal candidate point v p Push it into the path queue, and continue to make rays passing through the query point as the starting point, repeat the above steps until a duplicate point appears in the path queue, and take the last 3 points of the path queue, that is, form a triangle loop based on the rays between the duplicate points, and finally get the triangle containing the query point p i The triangle ΔS, where the schematic diagram for determining the optimal candidate point is as follows Figure 2 shown.

[0079] In this way, by determining a triangle containing the point to be queried, and performing calculations based on the parameter equation of the determined triangle containing the point to be queried, the elevation value in a large-scale scene can be calculated quickly and accurately.

[0080] When the query point p is determined i The triangle ΔS can get the query point p i Specifically, find the parametric equation f(x, y, z) of the plane where ΔS is located. The process of obtaining the parametric equation is as follows: (1) Let the three points of triangle ΔS be a(x1, y1, z1), b(x2, y2, z2), and c(x3, y3, z3), and find the two vectors: The cross product of the two vectors can be used to obtain the normal vector of the plane where ΔS is located, and the plane normal vector is brought into the query point p i The coordinate x i ,y i Get the elevation value z i .in, right Expand the determinant to get: according to The plane parametric equation f(x,y,z) can be written as:

[0081] A(x-x1)+B(y-y1)+C(z-z1)=0;

[0082] In the formula, A is ((y2-y1)(z3-z1)-(z2-z1)(y3-y1)), B is ((x2-x1)(z3-z1)-(z2-z1)(x3-x1)), C is ((x2-x1)(y3-y1)-(y2-y1)(x3-x1)), x1, y1 and z1 are the positions of the first point.

[0083] Therefore, the query point p i (x i ,y i ) into f(x,y,z) to get the corresponding elevation value z i , by repeating this method, the elevation values ​​corresponding to several candidate points in the vertex set V can be calculated.

[0084] In this way, by performing calculation based on the determined parametric equation of the triangle containing the point to be queried, the elevation value in a large-scale scene can be calculated quickly and accurately.

[0085] The following is a time complexity comparison of the present application with the linear search method and the jump and walk method, to prove that the present application has achieved better technical effects than the prior art, that is, it can quickly and accurately calculate the elevation values ​​in large-scale scenes. For the present application, ignoring the time-consuming process of Delaunay triangulation, assuming that the scale of the Delaunay triangulation is m, and the scale of the points to be queried is n, in the preprocessing stage, the algorithm takes O(m) to traverse the triangle, and the process of building an undirected graph uses a red-black tree data structure to store edges (that is, to record edge sets). Each triangle has 3 edges stored twice, that is, it is inserted into the red-black tree 6 times, and the overall complexity is O(mlog m), so the complexity in the KD tree building stage is O(mlog m); and in the query stage, the KD tree single query complexity is O(log m), taking the nearest point as the starting point, the path length of the path loop does not exceed 4, so the overall time complexity of querying n points is O(nlog m). Therefore, the overall complexity of the present application is: Under the same problem scale, the time complexity of the linear search method is: O(nm); and under the same problem scale, the time complexity of the jump and walk method is:

[0086] In the embodiment of the present application, by obtaining several elevation point sets to be calculated and triangulating the elevation point sets to obtain several Delaunay triangulations, a Delaunay triangulation can be quickly generated, which is convenient for the subsequent determination of the triangle where the query point is located; by projecting the Delaunay triangulation onto a two-dimensional plane to obtain a triangle set, and splitting each triangle in the triangle set to obtain a triangle undirected graph, the triangulation network in the three-dimensional space can be projected onto a two-dimensional plane, which is convenient for splitting each vertex set and edge set in each triangle, thereby facilitating the subsequent determination of the triangle containing the query point. Triangle and elevation calculation; by obtaining the point to be queried, and using the KD tree to search for the first point position in the undirected graph that is closest to the point to be queried, the first point position that is closest to the point to be queried can be quickly obtained, which is convenient for the subsequent determination of the triangle containing the point to be queried and the elevation calculation; based on the point to be queried and the first point position, the triangle containing the point to be queried is determined, and based on the parametric equation of the plane where the triangle is located, the elevation value corresponding to the point to be queried is obtained, and calculation can be performed based on the determined parametric equation of the triangle containing the point to be queried, so that the elevation value in large-scale scenes can be quickly and accurately calculated.

[0087] Embodiment 2

[0088] Please refer to Figure 3 , Figure 3 It is a structural schematic diagram of an embodiment of an elevation calculation system based on a Delaunay triangulation provided by the present application, comprising a first acquisition module 01, a first processing module 02, a second acquisition module 03 and a second processing module 04;

[0089] The first acquisition module 01 is used to acquire a number of elevation point sets to be calculated, and triangulate the elevation point sets to obtain a number of Delaunay triangulation networks; the first acquisition module 01 includes: an acquisition unit and an import unit; the acquisition unit is used to acquire a number of elevation points using RTK or a total station; the import unit is used to import the several elevation points into CAD to obtain a number of elevation point sets to be calculated.

[0090] It can be understood that several elevation points can be collected using RTK or total station; and by importing several of the elevation points into CAD, several sets of elevation points to be calculated can be obtained, where the elevation points are point information in a three-dimensional space, including x-axis coordinates, y-axis coordinates and z-axis coordinates.

[0091] In this way, several elevation points are collected in advance through RTK or total station, and imported into CAD to obtain several sets of elevation points to be calculated, which facilitates the subsequent construction of Delaunay triangulation.

[0092] After obtaining several elevation point sets, it is necessary to preprocess the elevation point sets and select the Delaunay triangulation algorithm to triangulate the preprocessed elevation point sets to obtain several Delaunay triangulation networks; wherein the triangulation algorithm includes but is not limited to the Bowyer-Watson algorithm, the Lawson local optimization algorithm (LOP), etc., for example but not limitation, taking the Bowyer-Watson algorithm for triangulation as an example, first, it is necessary to construct a sufficiently large super triangle as the initial triangular network, the vertices of which include all the points in the elevation point set; secondly, the points in the elevation point set are inserted into the triangular mesh one by one, for each newly inserted point, the influencing triangle containing the point is found, and the common edges of the influencing triangles are deleted, and the new point is connected to the vertices of the influencing triangle to form three new triangles. Again, perform an empty circumscribed circle check on each newly formed triangle to determine whether it meets the Delaunay condition. If the circumscribed circle of a triangle contains other elevation points, it means that the triangle does not meet the Delaunay condition. At this time, it is necessary to optimize the triangle by exchanging diagonals until all triangles meet the empty circumscribed circle condition; then, repeat the above steps of inserting points, constructing triangles, and detecting and optimizing empty circumscribed circles until all elevation points are inserted into the triangular mesh; finally, delete the triangles related to the super triangle to obtain the final Delaunay triangulation.

[0093] It should be noted that data preprocessing methods include but are not limited to removing duplicate points, processing noise points and filling missing data, and these data preprocessing methods can be achieved through existing technologies, but data preprocessing methods are not the focus of this application, so they will not be expanded here.

[0094] It should be noted that the Bowyer-Watson algorithm is an incremental algorithm for generating Delaunay triangulation. Delaunay triangulation is a special triangulation method that generates a triangulated network that meets the Delaunay condition by gradually inserting points and rebuilding the triangular mesh until all points are added. The Lawson Local Optimization Procedure (LOP) is an algorithm used to optimize the triangulated network structure during the triangulation process to ensure that it meets the Delaunay triangulation characteristics. It optimizes triangles that do not meet the Delaunay condition by exchanging diagonals.

[0095] The first processing module 02 is used to project the Delaunay triangulation network into a two-dimensional plane to obtain a triangle set, and split each triangle in the triangle set to obtain a triangle undirected graph, wherein the undirected graph includes a vertex set and an edge set;

[0096] It is understandable that when a Delaunay triangulation is obtained, it is necessary to simplify the Delaunay triangulation into a two-dimensional plane for further analysis. Specifically, first, it is necessary to predetermine a suitable projection method to better retain the structural characteristics of the three-dimensional Delaunay triangulation; second, project each point in the Delaunay triangulation into a two-dimensional plane according to the selected projection method, and during the projection process, it is necessary to calculate the coordinates of each point on the two-dimensional plane and retain the relative position relationship between the points; finally, re-triangulate the two-dimensional plane according to the projected point set to obtain a set of triangles in the two-dimensional plane, wherein the triangulation method has been detailed in the first acquisition module, so it will not be repeated here.

[0097] When a triangle set on a two-dimensional plane is obtained, it is necessary to extract all unique vertices from the triangle set and collect all vertices to form a vertex set V, V = {v1, v2, ..., v n}, where v i Indicates the i-th vertex. Next, traverse each triangle in the triangle set and add the three edges connected to each triangle to the edge set E. Each edge consists of two vertices, indicating that the two vertices are adjacent in the triangle. For each triangle Δ ijk , whose vertices are v i , v j , v k , its three sides (v i ,v j ),(v i ,v k ),(v j ,v k ) is added to the set of non-phase edges E, that is, (v i ,v j )∈E, if and only if (v j ,v i )∈E. Finally, the vertex set V and edge set E are used to construct an undirected graph G=(V,E).

[0098] It should be noted that the projection method includes but is not limited to orthogonal projection, perspective projection and isometric projection, etc. Projecting the Delaunay triangulation into a two-dimensional plane is not the focus of this application, so it will not be expanded here.

[0099] It should be noted that an undirected graph is a concept in graph theory, in which the edge has no direction, that is, the edge connecting two vertices can be regarded as bidirectional. In an undirected graph, vertices and edges together constitute the basic elements of the graph, and the edge set generally defines the connection relationship between vertices. A triangular undirected graph can be understood as a special undirected graph, which contains multiple triangular structures formed by three vertices connected end to end by three edges.

[0100] The second acquisition module 03 is used to acquire the point to be queried, and use the KD tree to search the first point in the undirected graph that is closest to the point to be queried;

[0101] It can be understood that after obtaining the triangular undirected graph, a KD tree can be established based on the vertex set V to speed up the search efficiency of the nearest point, where the process of constructing the KD tree is as follows: (1) arbitrarily select a vertex from the vertex set V as the root node; (2) In each layer of the KD tree, it is necessary to cyclically select different dimensions as the partitioning dimensions; (3) For the selected partitioning dimension, calculate the median of all vertices on the dimension, and use the median as the root node to divide the remaining points into two left and right subtrees; (4) Repeat steps (2) and (3) for the vertices of the left and right subtrees respectively, and recursively construct the left subtree and the right subtree until the subtree is empty.

[0102] Specifically, the second acquisition module 03 includes: a comparison unit and an update unit; the comparison unit is used to compare the point to be queried with the root node. If the value of the point to be queried in the current dimension is less than the value of the root node and the root node has no subtree, the root node is determined to be the second point position, and the first distance between the second point position and the point to be queried is calculated; the update unit is used to search and calculate the second distance between the third point position in the parent node corresponding to the root node and the point to be queried with the shortest distance, and continuously update the second point position based on the second distance until there is no smaller distance, and determine that the second point position is the first point position closest to the point to be queried.

[0103] After the KD tree is constructed, the query point p can be obtained. i , and use the KD tree to search the undirected graph G with p i Euclidean distance to the nearest point v s The specific query process is as follows: First, select the point p to be queried i The first dimension is the current dimension d, and the query point is compared with the same dimension d t The root node p c Compare, if the query point p i The value at the current dimension is less than the root node p c When the value of , enter the root node pc Continue searching the left subtree of the node; otherwise, go to the root node p c Continue searching the right subtree of the node, and repeat the search until the root node p c If there is no subtree, then the root node p is determined c is a second point that is closer, and a first distance d between the second point and the point to be queried is calculated. c ; Secondly, search and calculate the root node p c The corresponding parent node p t The third point in the query point has the smallest second distance d a , and based on the second distance d a The second point position is continuously updated until there is no smaller distance, and the second point position is determined to be the first point position v closest to the query point. s By way of example and not limitation, if the second distance d a <First distance d c , then calculate the parent node p t and the query point p i The distance between t , if d t <d c , then its parent node p t For the node with the closest distance, update the second point, i.e. p c =p t , d c =d t , recursively search for p t The left and right subtrees of are repeatedly recursively searched to determine the point closest to the query point.

[0104] In this way, by using the KD tree to search the first point in the undirected graph that is closest to the point to be queried, the first point that is closest to the point to be queried can be quickly acquired, which facilitates the subsequent determination of the triangle containing the point to be queried and the elevation calculation.

[0105] The second processing module 04 is used to determine a triangle containing the point to be queried based on the point to be queried and the first point position, and obtain an elevation value corresponding to the point to be queried based on a parametric equation of the plane where the triangle is located.

[0106] The second processing module 04 includes: a third processing unit, a calculation unit, a second comparison unit and a construction unit; the third processing unit is used to make a first ray passing through the point to be queried with the first point as the starting point; the calculation unit is used to extract a number of candidate points adjacent to the first point from the undirected graph, and make a second ray and a third ray passing through the first point with the first candidate point and the second candidate point respectively, and calculate the angles between the first ray and the second ray and the first ray respectively to obtain the first angle and the second angle; the second comparison unit is used to compare the first angle with the second angle, and if the second angle is smaller than the first angle, determine that the second candidate point corresponding to the second angle is the optimal candidate point; the construction unit is used to continue making rays passing through the point to be queried with the optimal candidate point as the starting point until a repeated point appears, and obtain a triangle containing the point to be queried based on the rays between the repeated points.

[0107] It can be understood that after determining the first point position closest to the query point, the next path node v can be searched. p (that is, the optimal candidate point below) to determine the triangle containing the query point. Specifically, first, take the first point v s Make the first ray l as the starting point s Passing the query point p i ; Secondly, randomly extract the first point v from the undirected graph G s The adjacent vertex set V includes several candidate points. The first candidate point v is randomly selected from the vertex set V. p , and take the first candidate point v p Draw the second ray l as the starting point p After the first point v s , where the initial optimal angle d is set p =2π, traverse the vertex set V and select the second candidate point v k , and take the second candidate point v k Make the third ray l as the starting point k After the first point v s , and calculate the second ray l p and the third ray l k With the first ray l s The first angle d is obtained p and the second angle d k ; Again, the first angle d p and the second angle d k For comparison, if the second angle d k is smaller than the first angle d p , then determine the second angle d kThe corresponding second candidate point v k is the optimal candidate point, even if v p =v k , and make d p =d k , otherwise keep v p is the optimal candidate point; finally, the optimal candidate point v p Push it into the path queue, and continue to make rays passing through the query point as the starting point, repeat the above steps until a duplicate point appears in the path queue, and take the last 3 points of the path queue, that is, form a triangle loop based on the rays between the duplicate points, and finally get the triangle containing the query point p i The triangle ΔS, where the schematic diagram for determining the optimal candidate point is as follows Figure 2 shown.

[0108] In this way, by determining a triangle containing the point to be queried, and performing calculations based on the parameter equation of the determined triangle containing the point to be queried, the elevation value in a large-scale scene can be calculated quickly and accurately.

[0109] When the query point p is determined i The triangle ΔS can get the query point p i Specifically, find the parametric equation f(x, y, z) of the plane where ΔS is located. The process of obtaining the parametric equation is as follows: (1) Let the three points of triangle ΔS be a(x1, y1, z1), b(x2, y2, z2), and c(x3, y3, z3), and find the two vectors: The cross product of the two vectors can be used to obtain the normal vector of the plane where ΔS is located, and the plane normal vector is brought into the query point p i The coordinate x i ,y i Get the elevation value z i .in, right Expand the determinant to get: according to The plane parametric equation f(x,y,z) can be written as:

[0110] A(x-x1)+B(y-y1)+C(z-z1)=0;

[0111] In the formula, A is ((y2-y1)(z3-z1)-(z2-z1)(y3-y1)), B is ((x2-x1)(z3-z1)-(z2-z1)(x3-x1)), C is ((x2-x1)(y3-y1)-(y2-y1)(x3-x1)), x1, y1 and z1 are the positions of the first point.

[0112] Therefore, the query point p i (x i ,y i ) into f(x,y,z) to get the corresponding elevation value z i , by repeating this method, the elevation values ​​corresponding to several candidate points in the vertex set V can be calculated.

[0113] In this way, by performing calculation based on the determined parametric equation of the triangle containing the point to be queried, the elevation value in a large-scale scene can be calculated quickly and accurately.

[0114] The following is a time complexity comparison of the present application with the linear search method and the jump and walk method, to prove that the present application has achieved better technical effects than the prior art, that is, it can quickly and accurately calculate the elevation values ​​in large-scale scenes. For the present application, ignoring the time-consuming process of Delaunay triangulation, assuming that the scale of the Delaunay triangulation is m, and the scale of the points to be queried is n, in the preprocessing stage, the algorithm takes O(m) to traverse the triangle, and the process of building an undirected graph uses a red-black tree data structure to store edges (that is, to record edge sets). Each triangle has 3 edges stored twice, that is, it is inserted into the red-black tree 6 times, and the overall complexity is O(mlog m), so the complexity in the KD tree building stage is O(mlog m); and in the query stage, the KD tree single query complexity is O(log m), taking the nearest point as the starting point, the path length of the path loop does not exceed 4, so the overall time complexity of querying n points is O(nlog m). Therefore, the overall complexity of the present application is: Under the same problem scale, the time complexity of the linear search method is: O(nm); and under the same problem scale, the time complexity of the jump and walk method is:

[0115] In the embodiment of the present application, by obtaining several elevation point sets to be calculated and triangulating the elevation point sets to obtain several Delaunay triangulations, a Delaunay triangulation can be quickly generated, which is convenient for the subsequent determination of the triangle where the query point is located; by projecting the Delaunay triangulation onto a two-dimensional plane to obtain a triangle set, and splitting each triangle in the triangle set to obtain a triangle undirected graph, the triangulation network in the three-dimensional space can be projected onto a two-dimensional plane, which is convenient for splitting each vertex set and edge set in each triangle, thereby facilitating the subsequent determination of the triangle containing the query point. Triangle and elevation calculation; by obtaining the point to be queried, and using the KD tree to search for the first point position in the undirected graph that is closest to the point to be queried, the first point position that is closest to the point to be queried can be quickly obtained, which is convenient for the subsequent determination of the triangle containing the point to be queried and the elevation calculation; based on the point to be queried and the first point position, the triangle containing the point to be queried is determined, and based on the parametric equation of the plane where the triangle is located, the elevation value corresponding to the point to be queried is obtained, and calculation can be performed based on the determined parametric equation of the triangle containing the point to be queried, so that the elevation value in large-scale scenes can be quickly and accurately calculated.

[0116] The specific embodiments described above further describe the purpose, technical solutions and beneficial effects of the present application in detail. It should be understood that the above description is only a specific embodiment of the present application and is not intended to limit the scope of protection of the present application. It is particularly pointed out that for those skilled in the art, any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the scope of protection of the present application.

Claims

1. A method for calculating elevation based on Delaunay triangulation, characterized in that: include: Acquire a plurality of elevation point sets to be calculated, and triangulate the elevation point sets to obtain a plurality of Delaunay triangulations; Projecting the Delaunay triangulation network onto a two-dimensional plane to obtain a triangle set, and splitting each triangle in the triangle set to obtain a triangle undirected graph, wherein the undirected graph includes a vertex set and an edge set; Obtain a point to be queried, and use a KD tree to search for the first point in the undirected graph that is closest to the point to be queried; Based on the point to be queried and the first point position, a triangle containing the point to be queried is determined, and based on a parametric equation of a plane where the triangle is located, an elevation value corresponding to the point to be queried is obtained.

2. The elevation calculation method based on Delaunay triangulation according to claim 1 is characterized in that: The step of obtaining a plurality of elevation point sets to be calculated is as follows: Collect several elevation points using RTK or total station; Importing a plurality of the elevation points into CAD, a plurality of elevation point sets to be calculated are obtained.

3. The elevation calculation method based on Delaunay triangulation according to claim 1 is characterized in that: The method of using the KD tree to search the first point in the undirected graph that is closest to the query point is specifically as follows: Compare the point to be queried with the root node, if the value of the point to be queried in the current dimension is less than the value of the root node and the root node has no subtree, determine that the root node is a second point, and calculate a first distance between the second point and the point to be queried; Search and calculate the second distance between the third point in the parent node corresponding to the root node and the point to be queried that is the shortest, and continuously update the second point based on the second distance until there is no smaller distance, and determine that the second point is the first point that is closest to the point to be queried.

4. The elevation calculation method based on Delaunay triangulation according to claim 1 is characterized in that: The determining of a triangle including the point to be queried based on the point to be queried and the first point position is specifically: Taking the first point as the starting point, draw a first ray passing through the point to be queried; Extracting a plurality of candidate points adjacent to the first point from the undirected graph, using the first candidate point and the second candidate point as a second ray and a third ray passing through the first point, respectively, and calculating the angles between the first ray and the second ray and the first ray, respectively, to obtain a first angle and a second angle; Compare the first angle with the second angle, and if the second angle is smaller than the first angle, determine the second candidate point corresponding to the second angle as the optimal candidate point; The optimal candidate point is used as a starting point to continue drawing rays through the point to be queried until a duplicate point appears, and a triangle containing the point to be queried is obtained based on the rays between the duplicate points.

5. The elevation calculation method based on Delaunay triangulation according to claim 1 is characterized in that: The parametric equation of the plane where the triangle is located is specifically: A(x-x1)+B(y-y1)+C(z-z1)=0; In the formula, A is ((y2-y1)(z3-z1)-(z2-z1)(y3-y1)), B is ((x2-x1)(z3-z1)-(z2-z1)(x3-x1)), C is ((x2-x1)(y3-y1)-(y2-y1)(x3-x1)), x1, y1 and z1 are the positions of the first point.

6. An elevation calculation system based on Delaunay triangulation, characterized in that: include: A first acquisition module, a first processing module, a second acquisition module and a second processing module; The first acquisition module is used to acquire a plurality of elevation point sets to be calculated, and triangulate the elevation point sets to obtain a plurality of Delaunay triangulations; The first processing module is used to project the Delaunay triangulation network into a two-dimensional plane to obtain a triangle set, and split each triangle in the triangle set to obtain a triangle undirected graph, wherein the undirected graph includes a vertex set and an edge set; The second acquisition module is used to acquire the point to be queried, and use the KD tree to search the first point in the undirected graph that is closest to the point to be queried; The second processing module is used to determine a triangle containing the point to be queried based on the point to be queried and the first point position, and obtain an elevation value corresponding to the point to be queried based on a parametric equation of the plane where the triangle is located.

7. The elevation calculation system based on Delaunay triangulation according to claim 6 is characterized in that: The first acquisition module includes: a collection unit and an import unit; The acquisition unit is used to acquire a number of elevation points using RTK or a total station; The importing unit is used to import the plurality of elevation points into CAD to obtain a plurality of elevation point sets to be calculated.

8. The elevation calculation system based on Delaunay triangulation according to claim 6, characterized in that: The second acquisition module includes: a comparison unit and an update unit; The comparison unit is configured to compare the to-be-queried point with a root node, and if the value of the to-be-queried point in the current dimension is less than the value of the root node and the root node has no subtree, determine that the root node is a second point position, and calculate a first distance between the second point position and the to-be-queried point; The updating unit is used to search for and calculate the second distance between the third point position in the parent node corresponding to the root node and the point to be queried, and continuously update the second point position based on the second distance until there is no smaller distance, and determine that the second point position is the first point position closest to the point to be queried.

9. The elevation calculation system based on Delaunay triangulation according to claim 6, characterized in that: The second processing module includes: a third processing unit, a calculation unit, a second comparison unit and a construction unit; The third processing unit is used to make a first ray passing through the to-be-queried point from the first point as a starting point; The calculation unit is used to extract a plurality of candidate points adjacent to the first point from the undirected graph, use the first candidate point and the second candidate point as a second ray and a third ray to pass through the first point, and calculate the angles between the first ray and the second ray and the first ray, respectively, to obtain a first angle and a second angle; The second comparison unit is used to compare the first angle with the second angle, and if the second angle is smaller than the first angle, determine that the second candidate point corresponding to the second angle is the optimal candidate point; The construction unit is used to take the optimal candidate point as a starting point and continue to make rays passing through the point to be queried until a repeated point appears, and obtain a triangle containing the point to be queried based on the rays between the repeated points.

10. The elevation calculation system based on Delaunay triangulation according to claim 6, characterized in that: The parametric equation of the plane where the triangle is located is specifically: A(x-x1)+B(y-y1)+C(z-z1)=0; In the formula, A is ((y2-y1)(z3-z1)-(z2-z1)(y3-y1)), B is ((x2-x1)(z3-z1)-(z2-z1)(x3-x1)), C is ((x2-x1)(y3-y1)-(y2-y1)(x3-x1)), x1, y1 and z1 are the positions of the first point.

Images