Complex terrain multi-resolution modeling method based on lod and triangulated irregular network folding

Abstract

The application relates to the technical field of geographic information processing, and discloses a complex terrain multi-resolution modeling method based on LOD and triangulation net folding, which comprises the following steps: firstly, constructing an irregular triangulation net and extracting a landform skeleton network, then carrying out spatial block division according to sampling distribution density and elevation change rate; for each block, calculating an entanglement sensitivity index representing the response of ridge-valley spatial tightness to the change of analysis granularity, and constructing a comprehensive evaluation system containing a geometric deviation component and a skeleton structure maintaining component according to the index; finally, performing controlled grid folding optimization based on the evaluation system, and performing weighted parameterized suture on the boundaries of adjacent blocks by using the above sensitivity index. The application can accurately perceive and retain the key landform skeleton in a complex mountainous scene, avoid structural topological errors in the simplification process, and ensure the spatial structure consistency and splicing smoothness of the multi-resolution model at different levels.

Description

Technical Field

[0001] This invention relates to the field of geographic information processing technology, and more specifically, to a multi-resolution modeling method for complex terrain based on LOD and triangulation folding. Background Technology

[0002] In applications such as digital earth, disaster simulation, and military simulation, real-time rendering of large-scale, high-precision 3D terrain environments is a core requirement. Due to the massive amount of terrain data, multi-resolution modeling techniques are widely used in industry to balance rendering efficiency and visual accuracy. The typical construction process usually includes: first, building a basic irregular triangular mesh; second, using spatial indexing structures such as quadtrees to manage terrain in blocks; and finally, generating models with different levels of detail based on viewpoint distance using algorithms such as edge folding or triangle shrinking. In existing terrain simplification and modeling methods, the criteria for determining whether a mesh needs to be folded or simplified typically rely on two main indicators: the distribution density of sampling points and the degree of geometric undulation of the terrain surface (such as root mean square error of elevation, curvature, or roughness). The basic logic is to retain more meshes in areas with significant undulations and drastically reduce the meshes in flat areas to achieve data lightweighting.

[0003] However, the conventional methods based on geometric undulations have significant limitations when dealing with complex mountainous scenes containing evolutionary traces such as river capture, wind gaps, and watershed migration. In these special geomorphic scenes, a ridge-valley role reversal phenomenon often occurs, where the original ancient river valley has transformed into a watershed passage due to geological evolution. Such areas are often spatially characterized as gently sloping wide valleys or low-lying saddle zones, with ridge lines and valley lines being extremely close and parallel in spatial location.

[0004] Because current technologies primarily allocate mesh resources based on the degree of geometric undulation, they often misclassify low-undulation ridge-valley zones as simple flat areas, leading to intensive mesh simplification. This approach easily causes the spatial structure of ridges and valleys to become stuck or incorrectly merged, resulting in positional drift or broken connections of key geomorphic framework structures (such as watershed boundaries and runoff paths) in the simplified model. This not only disrupts the topological consistency of the terrain but also leads to structural misalignments and visual abrupt changes that cannot be eliminated by conventional crack repair methods when switching between different levels of detail, severely impacting the realism of the representation and the accuracy of the analysis of complex terrain environments. Summary of the Invention

[0005] This invention provides a multi-resolution modeling method for complex terrain based on LOD and triangular mesh folding, which solves the technical problems mentioned in the background art.

[0006] A multi-resolution modeling method for complex terrain based on LOD and triangulation folding includes:

[0007] A basic irregular triangular network is constructed and a geomorphic skeleton network is extracted. The geomorphic skeleton network includes a set of valley lines based on the classification of water catchment parameters and a set of ridge lines based on the classification of watershed diffusion parameters.

[0008] Based on the sampling point distribution density and the rate of change of terrain elevation, the basic irregular triangular network is recursively divided into quadtrees to obtain multiple terrain blocks;

[0009] For each terrain block, calculate the ridge-valley spatial entanglement sensitivity, which characterizes the degree of response of the spatial proximity aggregation value from the valley line set to the ridge line set within the block to the change of spatial analysis granularity.

[0010] A simplified comprehensive error evaluation index for triangulation is constructed. The comprehensive error evaluation index includes a geometric deviation component and a skeleton structure preservation component weighted by the spatial entanglement sensitivity of the ridge valley. The skeleton structure preservation component is used to quantify the degree of deviation of the simplification operation from the distance field from the vertex to the terrain skeleton network.

[0011] The retention complexity of terrain blocks is calculated based on the ridge and valley spatial entanglement sensitivity to determine the number of target grids, and triangular mesh folding optimization is performed in combination with the comprehensive error evaluation index to generate a multi-resolution terrain model.

[0012] The beneficial effects of this invention include: by introducing a sensitivity index that characterizes the spatial entanglement of ridge and valley sets throughout the multi-resolution modeling process, it overcomes the defect of conventional methods that rely solely on sampling distribution and elevation undulations, resulting in the loss of the geomorphic skeleton structure; by deeply integrating this index into the grid simplification evaluation system and the block boundary stitching process, it effectively suppresses topological drift and misalignment in complex geomorphic scenarios such as river capture or wind gaps, achieving a significant reduction in data volume while accurately maintaining the spatial continuity and visual realism of the key mountain geomorphic skeleton, and significantly improving the structural stability and seamless stitching of the 3D terrain model at different display levels. Attached Figure Description

[0013] Figure 1 This is a flowchart of the multi-resolution modeling method for complex terrain based on LOD and triangular mesh folding of the present invention. Detailed Implementation

[0014] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0015] like Figure 1 As shown, a multi-resolution modeling method for complex terrain based on LOD and triangular mesh folding includes:

[0016] A basic irregular triangular network is constructed and a geomorphic skeleton network is extracted. The geomorphic skeleton network includes a set of valley lines based on the classification of water catchment parameters and a set of ridge lines based on the classification of watershed diffusion parameters.

[0017] Based on the sampling point distribution density and the rate of change of terrain elevation, the basic irregular triangular network is recursively divided into quadtrees to obtain multiple terrain blocks;

[0018] For each terrain block, calculate the ridge-valley spatial entanglement sensitivity, which characterizes the degree of response of the spatial proximity aggregation value from the valley line set to the ridge line set within the block to the change of spatial analysis granularity.

[0019] A simplified comprehensive error evaluation index for triangulation is constructed. The comprehensive error evaluation index includes a geometric deviation component and a skeleton structure preservation component weighted by the spatial entanglement sensitivity of the ridge valley. The skeleton structure preservation component is used to quantify the degree of deviation of the simplification operation from the distance field from the vertex to the terrain skeleton network.

[0020] The retention complexity of terrain blocks is calculated based on the ridge and valley spatial entanglement sensitivity to determine the number of target grids, and triangular mesh folding optimization is performed in combination with the comprehensive error evaluation index to generate a multi-resolution terrain model.

[0021] In a preferred embodiment, constructing a basic irregular triangular mesh and extracting a terrain skeleton network includes:

[0022] Terrain sampling point set The basic irregular triangulation is obtained by performing Delaunay triangulation on the plane. ;

[0023] Computing each vertex in a basic irregular triangular mesh Convergence accumulation :

[0024]

[0025] in, As vertices The representative area, The flow direction points to the vertex. The set of upstream adjacent vertices, upstream adjacent vertex Flowing to the top Multi-flow allocation weights;

[0026] Computing each vertex in a basic irregular triangular mesh cumulative diffusion :

[0027]

[0028] in, As vertices The set of adjacent vertices flowing downstream. As vertices Flowing downstream adjacent vertices Dual allocation weights;

[0029] Extract valley line set With ridge line assembly The valley line set is composed of the confluence accumulation. The ridge set consists of connected paths exceeding a set threshold, and is composed of diffusion accumulation. Connectivity paths exceeding a set threshold constitute the network.

[0030] It should be noted that the terrain sampling point set is a collection of multiple three-dimensional coordinate points on the terrain surface, which can be obtained through LiDAR scanning aerial photogrammetry or ground measurement. The basic irregular triangular mesh (BITM) is a triangular mesh model formed by delaunay triangulation of the terrain sampling point set. Vertices are the basic nodes that form triangles in the BITM. Elevation values ​​are the vertical height data of the terrain corresponding to the vertices, which can be obtained through LiDAR scanning aerial photogrammetry or ground measurement. The confluence accumulation is the weighted sum of the vertex's own representative area and the confluence accumulation of all upstream adjacent vertices flowing towards that vertex, used to characterize the degree of water flow convergence at the vertex. The representative area is the area covered by each vertex in its planar projection, usually taken as the area of ​​the corresponding Voronoi polygon. The upstream adjacent vertex set is the set of adjacent vertices on the terrain surface whose water flow direction points towards the current vertex. The multi-flow direction allocation weight is the proportion of water flow allocated from upstream adjacent vertices to the current vertex, preferably a value between 0 and 1, based on the terrain slope and aspect relationship. When using the D-infinity method for allocation, the weight is positively correlated with the slope between adjacent vertices. The diffusion accumulation is a weighted sum of the area represented by the vertex itself and the diffusion accumulation of all its downstream adjacent vertices, used to characterize the degree of water diffusion at the vertex. The downstream adjacent vertex set is the set of vertices on the terrain surface whose water flow direction originates from the current vertex. The dual allocation weight is the proportion of water flow allocated from the current vertex to its downstream adjacent vertices, preferably a value between 0 and 1, based on the same criteria as the multi-flow allocation weight, ensuring duality with the confluence accumulation calculation logic. The valley set is the set of connected paths whose confluence accumulation exceeds a set threshold, used to characterize gully lines in the terrain. The ridge set is the set of connected paths whose diffusion accumulation exceeds a set threshold, used to characterize ridge lines in the terrain. The confluence accumulation threshold is the critical value for filtering the valley set, preferably the 80th percentile of the total confluence accumulation, determined by statistically analyzing the distribution of total confluence accumulation across the entire region to retain connected paths in high-convergence areas. The diffusion accumulation threshold is the critical value for filtering the ridge line set. The preferred value is the 80th percentile of the global diffusion accumulation. The basis for this value is to retain the connected paths in high diffusion areas by statistically analyzing the global diffusion accumulation distribution.

[0031] It should be noted that the confluence accumulation and diffusion accumulation employ a dual calculation logic. The confluence accumulation is calculated based on the forward flow direction, while the diffusion accumulation is calculated based on the reverse flow direction, with their calculation forms being completely symmetrical. This design ensures consistency in the extraction logic of valley lines and ridge lines, avoiding imbalances in the skeleton network caused by a single extraction method. For example, in a mountainous terrain, water flows from the ridge to the valley. The confluence accumulation gradually increases at the valley apex, while the diffusion accumulation gradually increases at the ridge apex. Dual calculations can accurately capture the core routes of both valleys and ridges simultaneously. The extraction of both valley line sets and ridge line sets uses a dual screening rule combining accumulation thresholds and connected paths, rather than relying solely on a single accumulation value. This rule excludes isolated vertices with high accumulation, ensuring that the extracted skeleton network has continuous topographical significance. For example, if an isolated vertex with high confluence accumulation exists in a certain area, it will not be included in the valley line set because it does not meet the connected path condition, preventing breaks or abnormal branches in the skeleton network.

[0032] It should be noted that the specific allocation rule for multi-flow direction weighting adopts the D-infinity method. This method allocates weights based on the slope angle between the current vertex and its upstream adjacent vertices. The larger the slope angle, the larger the weight, and the sum of the weights of all upstream adjacent vertices is 1. For example, if the current vertex has three upstream adjacent vertices with corresponding slope angles of 30 degrees, 45 degrees, and 60 degrees, the weight allocation ratio is 0.2, 0.3, and 0.5, with a total of 1. The specific allocation rule for dual weighting is consistent with the multi-flow direction weighting, using the D-infinity method. It allocates weights based on the slope angle between the current vertex and its downstream adjacent vertices, ensuring duality with the flow allocation logic for calculating the confluence accumulation. For example, if the current vertex has two downstream adjacent vertices with corresponding slope angles of 30 degrees and 60 degrees, the weight allocation ratio is 0.3 and 0.7, with a total of 1. The specific method for setting the confluence accumulation threshold is as follows: calculate the confluence accumulation of all vertices in the entire domain, sort them from smallest to largest, and take the value corresponding to the 80th percentile of the sorted value as the threshold. For example, if the entire domain has 1000 vertices, and the value of the 800th vertex after sorting the cumulative diffusion amount is 500 square meters, then the cumulative diffusion amount threshold is set to 500 square meters. The specific method for setting the cumulative diffusion amount threshold is the same as the cumulative diffusion amount threshold: calculate the cumulative diffusion amount of all vertices in the entire domain, sort them from smallest to largest, and take the value corresponding to the 80th percentile of the sorted value as the threshold. For example, if the entire domain has 1000 vertices, and the value of the 800th vertex after sorting the cumulative diffusion amount is 400 square meters, then the cumulative diffusion amount threshold is set to 400 square meters.

[0033] It should be noted that, based on the analysis of water flow direction, dual calculations of confluence and diffusion accumulation are used to extract valley lines and ridge lines from the terrain, forming a complete geomorphic skeleton network. This design fully utilizes the correlation between the hydrological characteristics and geometric morphology of the terrain, ensuring that the extracted skeleton network conforms to the actual terrain evolution patterns. Its beneficial effects include good geometric adaptability of the constructed basic irregular triangular network, preserving key terrain features. The extracted valley and ridge line sets have clear physical meaning and continuity, providing reliable skeleton support for subsequent multi-resolution modeling and preventing the loss of core terrain structures during subsequent simplification. Simultaneously, the dual calculation logic ensures the integrity and consistency of the skeleton network, improving the stability and reliability of the entire modeling method.

[0034] In a preferred embodiment, the basic irregular triangular network is recursively partitioned using a quadtree based on the sampling point distribution density and the rate of change of terrain elevation, resulting in multiple terrain blocks, including:

[0035] Calculate the current terrain tiles Partition evaluation value :

[0036]

[0037] in, The sampling point distribution density for the current terrain patch. This is a median reference value for the distribution density of sampling points across the entire region. The rate of change of terrain elevation for the current terrain segment. This is the median reference value for the rate of change of topographic elevation across the entire region. Density weighting coefficient, This is the weighting coefficient for the rate of change of elevation;

[0038] Check the partition evaluation value Does the recursion termination condition meet?

[0039]

[0040] in, This is the preset partition stop threshold;

[0041] If the termination condition is not met, the current terrain block is subdivided into quadtrees; if the termination condition is met, the current terrain block is output.

[0042] It should be noted that the current terrain tile is the terrain region unit being processed during the quadtree partitioning process. The partition evaluation value is a quantitative indicator measuring the complexity of the current terrain tile, obtained by a weighted sum of the normalized sampling point distribution density and the normalized terrain elevation change rate. The sampling point distribution density is the density of sampling points within the current terrain tile, calculated by dividing the number of sampling points within the tile by the tile area. The median reference value for the global sampling point distribution density is the value in the middle after sorting the sampling point distribution densities of all terrain tiles within the entire terrain region from smallest to largest. The terrain elevation change rate is the standard deviation of the elevation values ​​of all vertices within the current terrain tile, used to characterize the degree of terrain undulation. The median reference value for the global terrain elevation change rate is the value in the middle after sorting the terrain elevation change rates of all terrain tiles within the entire terrain region from smallest to largest. The density weighting coefficient adjusts the influence of sampling point distribution density on the zoning evaluation value, with a preferred value of 0.4. This value is chosen to balance the influence of sampling point density and the rate of change of terrain elevation, ensuring that the zoning results simultaneously adapt to geometric density and terrain undulation characteristics. The elevation change rate weighting coefficient adjusts the influence of the rate of change of terrain elevation on the zoning evaluation value, with a preferred value of 0.6. This value follows the same criteria as the density weighting coefficient, and their sum is 1, ensuring the rationality of the evaluation value calculation. The zoning stop threshold is the critical value for determining whether the current terrain block should stop quadtree partitioning, with a preferred value of 1.5. This value is determined by statistically analyzing the zoning evaluation value distribution of the initial blocks across the entire region, avoiding excessive computation due to overly fine partitioning or loss of terrain details due to overly coarse partitioning.

[0043] It should be noted that the zoning evaluation value adopts a calculation logic of normalization combined with weighted summation. First, the sampling point distribution density and the rate of change of terrain elevation are normalized by dividing each by their respective global median reference value. Then, they are summed using weighted coefficients. This design eliminates the influence of differences in the dimensions of the two indicators, giving the evaluation value a unified comparison standard. For example, if the sampling point distribution density of a certain terrain block is 20 points per square meter, and the global median reference value is 10 points per square meter, the normalized value is 2. The rate of change of terrain elevation is 8 meters, and the global median reference value is 4 meters, the normalized value is 2. Using weighted coefficients of 0.4 and 0.6, the zoning evaluation value is 2 × 0.4 + 2 × 0.6 = 2. If it exceeds the threshold of 1.5, further subdivision is required. The termination condition of the quadtree recursive partitioning depends only on the single zoning evaluation value, without introducing other branch determination conditions. This design ensures the determinism of the partitioning; the same terrain data input will always yield the same partitioning result. For example, when dividing terrain into different areas, regardless of their geographical location, the division stops as long as the partition evaluation value is less than or equal to 1.5, thus ensuring the consistency of the partitioning rules.

[0044] It should be noted that the specific rules for determining the values ​​of the density weighting coefficient and the elevation change rate weighting coefficient are that their sum is fixed at 1, and the ratio can be adjusted according to the characteristics of the terrain data. When the density distribution of the sampling points in the terrain data varies greatly, the density weighting coefficient can be increased to 0.5, and the elevation change rate weighting coefficient can be decreased to 0.5; when the terrain undulations are more significant, the elevation change rate weighting coefficient can be kept at 0.6, and the density weighting coefficient at 0.4. For example, for plain areas, where the sampling point density difference is small and the terrain undulations are gentle, both can be adjusted to 0.5; for mountainous areas, where the terrain undulations are large, the ratio of 0.4 and 0.6 should be maintained. The specific method for setting the partition stopping threshold is as follows: calculate the partition evaluation value of all initial root node blocks in the entire region, and take 1.5 times their average value as the final threshold. For example, if the initial global domain has 10 root node blocks with partition evaluation values ​​of 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, and 2.8 respectively, and the average value is 1.96, then the threshold is set to 1.96 × 1.5 ≈ 2.94 to ensure that most initial blocks can be reasonably subdivided according to complexity.

[0045] It should be noted that, based on the sampling point distribution density and elevation change rate of the terrain, a partition evaluation value is constructed through normalization and weighted summation. Complex terrain is then divided into multiple suitable terrain blocks using a quadtree recursive partitioning method. By utilizing the global median reference value to eliminate the influence of extreme values, and by balancing the effects of sampling point distribution density and elevation change rate through weighting coefficients, the partitioning results are ensured to match the actual complexity of the terrain. This results in terrain block sizes that are appropriate for the terrain complexity, with more refined blocks for complex areas and coarser blocks for flat areas. This avoids the loss of detail or computational redundancy caused by uniform partitioning, and provides uniform and terrain-feature-fitting processing units for subsequent steps such as sensitivity calculation and triangulation folding, improving the efficiency and stability of the entire modeling process.

[0046] In a preferred embodiment, for each terrain block, the ridge-valley spatial entanglement sensitivity, which characterizes the degree of response of the spatial proximity aggregation value from the valley line set to the ridge line set within the block to changes in spatial analysis granularity, is calculated, including:

[0047] Get current terrain tiles Valley line set subset subset of the ridge set ;

[0048] exist Upper sampling point Calculate the distance from each sampling point to shortest distance ;

[0049] Calculate the analysis granularity of the current terrain patch in any space. Spatial proximity aggregation value :

[0050]

[0051] in, This represents the sampling step size along the valley line;

[0052] Calculate the ridge-valley spatial entanglement sensitivity of the current terrain patch. :

[0053]

[0054] in, For a fixed granularity scaling ratio, The basic analysis granularity for the current terrain partitioning is defined as follows:

[0055]

[0056] in, The area of ​​the current terrain is divided into blocks. This represents the number of basic vertices within the block. It is a proportionality constant.

[0057] It should be noted that terrain segments are independent terrain processing units obtained after recursive partitioning using a quadtree. The valley set subset is the set of valley line segments extracted from the current terrain segment. The ridge set subset is the set of ridge line segments extracted from the current terrain segment. Discrete sampling points are discrete points extracted from the valley set subset at a fixed arc length step. Euclidean distance is the straight-line distance from each discrete sampling point to the ridge set subset. Spatial analysis granularity is the observation scale that measures spatial proximity, which can be understood as a detail resolution window. The spatial proximity aggregation value is a numerical value that quantifies the proximity between valley and ridge lines based on an exponential decay model, and is an approximate integral sum of the ratio of the negative distance of all discrete sampling points to the spatial analysis granularity. The sampling step size is a fixed arc length interval for extracting discrete sampling points along the valley line, preferably half of the basic analysis granularity. The value is chosen to control the computational load while ensuring sampling accuracy, so that the sampling points can fully reflect the valley morphology. Ridge-valley spatial entanglement sensitivity is a quantitative indicator characterizing the responsiveness of spatial proximity aggregation values ​​from valley line sets to ridge line sets within a block to changes in spatial analysis granularity. The granularity scaling ratio is a fixed constant used to scale the spatial analysis granularity, preferably set to 2. This value is chosen because a doubling of the scale stably reflects the response characteristics of spatial proximity aggregation values, avoiding sensitivity calculation distortion caused by excessively small or large scale spans. The basic analysis granularity is a benchmark scale proportional to the average vertex spacing within the terrain block, used to unify the reference standard for spatial analysis. The terrain block area is the projected area of ​​the current terrain block on a two-dimensional plane. The number of basic vertices within the block is the total number of vertices in the basic irregular triangular mesh within the current terrain block. The scaling constant is a coefficient that maps the average vertex spacing to the basic analysis granularity, preferably set to 1. This value is chosen because it directly correlates the average vertex spacing with the basic analysis granularity, simplifying calculations while ensuring scale consistency.

[0058] It should be noted that the spatial proximity aggregation value is calculated using an exponential decay model, rather than a simple distance mean or summation. This model highlights the contribution of nearby sampling points, making the aggregation value more accurately reflect the actual proximity between valleys and ridges. For example, within a certain terrain patch, some discrete sampling points are 1 meter away from the ridge, while others are 10 meters away. When the spatial analysis granularity is 5 meters, the exponential term corresponding to the sampling point at a distance of 1 meter is 0.818, and that at a distance of 10 meters is 0.135. The aggregation value effectively highlights the proximity characteristic. The ridge-valley spatial entanglement sensitivity is calculated using logarithmic differencing, focusing on the response of the aggregation value to changes in spatial analysis granularity, rather than the proximity value at a single scale. This design can capture the scale dependence of valley-ridge entanglement. For example, in a conventional terrain patch, when the granularity increases from 5 meters to 10 meters, the logarithmic difference of the aggregation value is 0.3, and the sensitivity is 0.3; while in a patch with strong ridge-valley entanglement, the logarithmic difference is 0.8 and the sensitivity is 0.8 under the same granularity change, clearly distinguishing the two types of terrain. The granularity of the basic analysis is tied to the average vertex spacing within the terrain tiles, establishing a correlation between sensitivity calculation and the resolution of the terrain tiles themselves, thus avoiding subjectivity in setting the analysis scale. For example, for a tile with an average vertex spacing of 4 meters, the granularity of the basic analysis is 4 meters, ensuring that the spatial analysis scale matches the accuracy of terrain discretization.

[0059] It should be noted that the discrete sampling points are extracted as follows: along the centerline of the valley subset, starting from the starting point, points are sequentially extracted according to the set sampling step size until the end point. If the last segment is less than the sampling step size, the end point is included. For example, if the valley subset is 20 meters long and the sampling step size is 4 meters, a total of 5 discrete sampling points are extracted, with a spacing of 4 meters between each point. The Euclidean distance is calculated as follows: for each discrete sampling point, the shortest straight-line distance to all line segments of the ridge subset is calculated, and the minimum value is taken as the Euclidean distance of that sampling point. For example, if the distances from a sampling point to the three line segments of the ridge subset are 3 meters, 2.5 meters, and 4 meters respectively, the final Euclidean distance is taken as 2.5 meters. The specific calculation method for the granularity of the basic analysis is as follows: first, calculate the average vertex spacing within the terrain block, that is, the square root of the terrain block area divided by the square root of the number of basic vertices within the block, and then multiply by the proportionality constant to obtain the granularity of the basic analysis. For example, if the terrain block area is 10,000 square meters, the number of basic vertices within the block is 100, the average vertex spacing is 10 meters, and the scale constant is 1, then the basic analysis granularity is 10 meters.

[0060] It should be noted that the spatial proximity relationship between valley lines and ridge lines in complex terrain is scale-dependent. By extracting a subset of the skeleton within a terrain block, sampling at a fixed step size and calculating the Euclidean distance, an exponential decay model is used to obtain the spatial proximity aggregation value. Then, logarithmic difference is used to quantify the response of the aggregation value to changes in the analysis granularity, thus obtaining the ridge-valley spatial entanglement sensitivity. This design, combined with quadtree partitioning, can provide a dedicated quantification index of skeleton entanglement features for each block. Its beneficial effect is that the sensitivity calculation closely matches the actual structure of the terrain block, accurately capturing ridge-valley entanglement features that are difficult to detect by conventional methods. This provides a reliable basis for subsequent comprehensive error evaluation and triangulation folding optimization, enabling subsequent steps to specifically protect sensitive skeleton structures and improve the rationality and stability of multi-resolution modeling.

[0061] In a preferred embodiment, the geometric deviation component includes:

[0062] For the local vertex set Fold to new apex Candidate folding operations Calculate the geometric deviation components :

[0063]

[0064] in, For the new peak The three-dimensional coordinates The original vertices in the local vertex set The three-dimensional coordinates These are the vertical weighting coefficients, and .

[0065] It should be noted that candidate folding operations are vertex merging operations that may be performed during the triangulation process, including edge folding or triangle folding. The local vertex set is a group of adjacent vertices to be merged in the candidate folding operation. The new vertex is a single vertex generated after merging the local vertex set in the candidate folding operation. The geometric deviation component is a quantitative indicator measuring the degree of change in local terrain geometry caused by the candidate folding operation. The x-coordinate of the new vertex is its planar coordinate value in the horizontal x-direction. The y-coordinate of the new vertex is its planar coordinate value in the horizontal y-direction. The z-coordinate of the new vertex is its corresponding vertical elevation value. The x-coordinate of the original vertex is the planar coordinate value in the horizontal x-direction of each original vertex in the local vertex set. The y-coordinate of the original vertex is the planar coordinate value in the horizontal y-direction of each original vertex in the local vertex set. The z-coordinate of the original vertex is its corresponding vertical elevation value. The vertical weighting coefficient is a coefficient used to amplify the weight of errors in the vertical elevation direction. The preferred value is 2. The reason for this value is that elevation features have a more significant impact on terrain morphology in terrain modeling. By weighting, the elevation accuracy can be protected first, and the distortion of terrain undulation caused by folding can be avoided.

[0066] It should be noted that the geometric deviation component uses a calculation logic combining the squared error of the horizontal coordinates with the weighted squared error of the vertical axis, rather than calculating the error using equal weights on the three-dimensional coordinates. This design fully considers the characteristics of terrain data. Horizontal coordinates mainly reflect planar position, while vertical elevation directly determines the terrain undulation. Weighting highlights the importance of elevation error. For example, for a local vertex set of a mountainous terrain, the total horizontal coordinate error after folding is 3, and the total vertical elevation error is 2. Calculated using a vertical weighting coefficient of 2, the geometric deviation component is 3 + 2 × 2 = 7. If calculated using equal weights, it would be 5. The weighted design can more strictly constrain the distortion in the elevation direction. The vertical weighting coefficient is set to a fixed value greater than 1, rather than being dynamically adjusted, ensuring the consistency and determinism of the error calculation rules. This design avoids fluctuations in error evaluation standards caused by dynamic coefficients. For example, regardless of whether the terrain block is gentle or steep, a uniform vertical weighting coefficient of 2 is used, making the folding operation errors of different areas comparable.

[0067] It should be noted that the method for generating the 3D coordinates of the new vertex uses a least-squares fitting algorithm. Taking the 3D coordinates of the local vertex set as input, it solves for the optimal coordinates of the new vertex by minimizing the geometric deviation components from all original vertices to the new vertex. For example, if the local vertex set contains three vertices with coordinates (1,2,3), (2,3,4), and (3,4,5), the least-squares fitting calculation yields the new vertex coordinates as (2,3,4), where the geometric deviation components are minimized. The vertical weighting coefficient ranges from 1.5 to 3, with a preferred value of 2. This value is determined through statistical analysis of numerous terrain modeling cases, finding that this range strikes a balance between preserving elevation features and controlling computational complexity. When high terrain elevation accuracy is required, a value of 3 can be used; when computational efficiency needs to be considered, a value of 1.5 can be used. For example, in high-precision terrain modeling scenarios, a vertical weighting coefficient of 3 is used to strengthen elevation error constraints; in conventional scenarios, a value of 2 is used to balance accuracy and efficiency.

[0068] It should be noted that, considering the crucial role of elevation features in terrain modeling, this design separates the error calculations of horizontal coordinates and vertical elevation, introduces a vertical weighting coefficient to amplify the weight of elevation errors, and constructs a geometric deviation component to quantify the impact of folding operations on terrain geometry. This design aligns with the essential characteristics of terrain data, and its beneficial effects include prioritizing the preservation of terrain elevation information, reducing distortion of terrain undulations during triangulation folding, and ensuring that the simplified mesh retains the key geometric features of the original terrain. Simultaneously, the clearly defined error calculation rules ensure the consistency and repeatability of folding operation evaluation, providing a reliable quantitative basis for geometric errors in the subsequent construction of comprehensive error evaluation indicators, and improving the rationality and stability of triangulation folding optimization.

[0069] In a preferred embodiment, the skeleton structure retains components, including:

[0070] For the local vertex set Fold to new apex Candidate folding operations Calculate the skeleton structure preservation components :

[0071]

[0072] in, and These are the sets of ridge lines from the vertices before and after the collapse to the current terrain tile. The shortest distance;

[0073] As vertices The skeleton influence weight is defined as:

[0074]

[0075] in, As vertices Set of valley lines for current terrain segmentation The shortest distance, The skeleton affects the attenuation scale, and The value is equal to the basic analysis granularity of the current terrain patch. .

[0076] It should be noted that the candidate folding operation is the vertex merging operation to be performed during the triangulation process, including edge folding or triangle folding. The local vertex set is a group of adjacent vertices to be merged in the candidate folding operation. The new vertex is a single vertex generated after merging the local vertex set in the candidate folding operation. The skeleton structure preservation component is a quantitative indicator measuring the impact of the candidate folding operation on the stability of the terrain skeleton network. The shortest distance from the original vertex to the ridge set is the minimum straight-line distance from each original vertex in the local vertex set to the subset of ridges within the terrain block. The shortest distance from the original vertex to the valley set is the minimum straight-line distance from each original vertex in the local vertex set to the subset of valleys within the terrain block. The skeleton influence weight of a vertex is a numerical value that quantifies the degree of influence of the vertex on the terrain skeleton network, constructed based on an exponential decay function, and consists of the sum of distance decay terms to ridges and valleys. The skeleton influence decay scale is a characteristic scale in the exponential decay function that controls the weight decay rate. It is preferably set to equal the basic analysis granularity of the terrain block, and its value is determined by matching the block resolution to ensure that the weight calculation conforms to the terrain discretization accuracy. The shortest distance from the new vertex to the ridge set after folding is the minimum straight-line distance from the newly generated vertex to the subset of ridge sets within the terrain tile. The subset of ridge sets within a terrain tile is a set of ridge segments extracted from the current terrain tile. The subset of valley sets within a terrain tile is a set of valley segments extracted from the current terrain tile.

[0077] It should be noted that the skeleton influence weights employ a double exponential decay design, incorporating both the distances from vertices to ridges and valleys, rather than relying solely on a single skeleton element. This design comprehensively considers the relationship between vertices and the terrain skeleton network, making the weights more closely aligned with the actual terrain structure. For example, if a vertex is 3 meters from a ridge and 4 meters from a valley, with a base analysis granularity of 5 meters, the corresponding decay terms are 0.549 and 0.449 respectively, resulting in a skeleton influence weight of 0.998, reflecting the strong influence of the skeleton network on that vertex. The skeleton influence decay scale is forcibly bound to the base analysis granularity of the terrain sub-section, rather than being set independently. This design establishes a correlation between weight calculation and terrain sub-section resolution, avoiding subjectivity in weighting. For example, for a sub-section with a base analysis granularity of 6 meters, the decay scale is fixed at 6 meters, ensuring consistent weight calculation rules across different sub-sections and improving the comparability of error evaluations. The skeleton structure maintains the absolute change in the distance from vertices to ridges before and after component focusing folding, directly quantifying the stability impact of the skeleton network, rather than indirectly deriving it through geometric errors. For example, if the distance from a vertex to the ridge line is 5 meters before folding and 7 meters after folding, the absolute change is 2 meters. By combining the influence weight of the skeleton, the disturbance of the skeleton caused by folding can be accurately measured.

[0078] It should be noted that the method for calculating the shortest distance from a vertex to the ridge / valley set is as follows: decompose the ridge / valley subset into several line segments, calculate the perpendicular distance from the vertex to each line segment, and take the minimum of all perpendicular distances as the shortest distance. For example, if the perpendicular distances from a vertex to the three line segments of the ridge subset are 4 meters, 3.2 meters, and 5 meters, the shortest distance is 3.2 meters. The method for calculating the shortest distance from the new vertex to the ridge set after folding is as follows: first, obtain the 3D coordinates of the new vertex through least-squares fitting, and then use the same line segment decomposition-perpendicular distance calculation method as the original vertex to obtain the shortest distance from the new vertex to the ridge subset. For example, if the new vertex coordinates are (5,6,7), the ridge subset is decomposed into two line segments, and the calculated perpendicular distances are 2.8 meters and 3.5 meters, respectively. The shortest distance is 2.8 meters. The specific steps for calculating the skeleton influence weight are as follows: First, calculate the shortest distance from the vertex to the ridge and valley lines respectively. Then, divide each distance by the skeleton influence attenuation scale, take the negative value, and calculate the exponential function. Finally, add the results of the two exponential functions to obtain the skeleton influence weight. For example, if the distance from the ridge line is 2 meters, the distance from the valley line is 3 meters, and the attenuation scale is 5 meters, the weight is exp(-2 / 5) + exp(-3 / 5) ≈ 0.670 + 0.549 = 1.219.

[0079] It should be noted that, based on the supporting role of the geomorphic skeleton network in terrain modeling, by calculating the shortest distance from vertices to ridges and valleys, a skeleton influence weight bound to the block resolution is constructed. This quantifies the perturbation of the skeleton network's distance field by folding operations, forming skeleton structure preservation components. This design aligns with the stability requirements of the skeleton network, and its beneficial effect is that it can accurately capture the impact of folding operations on the core geomorphic structure, avoiding skeleton network distortion during simplification. Simultaneously, the explicit calculation rules ensure consistent evaluation, and a comprehensive error index is constructed in conjunction with the geometric deviation component. This allows the triangulation folding optimization to both consider geometric morphology and protect the skeleton structure, improving the structural integrity and reliability of the multi-resolution model.

[0080] In a preferred embodiment, calculating the preservation complexity of terrain patches based on the ridge-valley spatial entanglement sensitivity to determine the target grid number includes:

[0081] Calculate the current terrain tiles Preservation complexity :

[0082]

[0083] in, These represent the sampling point distribution density and the rate of change of terrain elevation for the current terrain segment, respectively. This is the corresponding global reference value. Sensitivity to entanglement in ridge valley space. For the corresponding exponential weighting constant, This is the sensitivity amplification factor;

[0084] Calculate the target grid number of the current terrain tile. :

[0085]

[0086] in, and These represent the minimum and maximum number of triangles allowed within a terrain segment, respectively. Rounded to the nearest integer.

[0087] It should be noted that terrain segments are independent terrain processing units obtained after recursive partitioning using a quadtree. Preservation complexity is a quantitative indicator that integrates the geometric features and entanglement features of the terrain segments, and is composed of the product of the density factor, elevation change rate factor, and structural risk factor. Sampling point distribution density is the density of sampling points within a terrain segment, calculated by dividing the number of sampling points within the segment by the segment area. The global reference value for sampling point distribution density is the median of the sampling point distribution density of all terrain segments within the entire terrain region. The elevation change rate is the standard deviation of the elevation values ​​of all vertices within a terrain segment, used to characterize the degree of terrain undulation. The global reference value for the elevation change rate is the median of the elevation change rate of all terrain segments within the entire terrain region. Ridge-valley spatial entanglement sensitivity is an indicator characterizing the responsiveness of the spatial proximity aggregation values ​​of valley lines and ridge lines within a segment to changes in analysis granularity. The density exponential weighting constant is a constant that adjusts the influence of the density factor on preservation complexity, preferably set to 1. This value is chosen to balance the effects of density with other factors, avoiding a single factor dominating the complexity calculation. The exponential weighting constant for the terrain elevation change rate is a constant that adjusts the influence of the elevation change rate factor on the retention complexity. A value of 1 is preferred, and its selection is consistent with the exponential weighting constant for density, ensuring a balanced contribution of geometric features. The exponential weighting constant for the structural risk factor is a constant that adjusts the influence of the structural risk factor on the retention complexity. A value of 1.5 is preferred, and its selection is based on strengthening the influence of skeleton entanglement features, enabling high-sensitivity blocks to achieve higher retention complexity. The sensitivity amplification factor is a coefficient that amplifies the contribution of ridge and valley spatial entanglement sensitivity to the structural risk factor. A value of 2 is preferred, and its selection is based on highlighting the differences in skeleton entanglement between different blocks, avoiding weak influence due to excessively low sensitivity values. The target mesh number is the number of triangles that the terrain block needs to retain at a specific level of detail. The minimum number of triangles is the lower bound of the allowed number of triangles within the terrain block, and its selection is based on ensuring the basic structural integrity of the block mesh, avoiding oversimplification. The maximum number of triangles is the upper limit of the number of triangles allowed within a terrain tile. The preferred value is 1000. The basis for this value is to control the amount of computation and avoid low modeling efficiency caused by too many meshes.

[0088] It should be noted that the retained complexity employs a three-factor product design, integrating the sampling point distribution density (geometric density), the rate of change of terrain elevation (geometric undulation), and the ridge-valley spatial entanglement sensitivity (structural risk), rather than determining complexity from a single dimension. This design comprehensively reflects the actual complexity of the blocks. For example, for a mountain block with a density ratio of 2, an elevation change rate ratio of 1.8, a ridge-valley spatial entanglement sensitivity of 1.2, a sensitivity amplification factor of 2, and a structural risk factor of 3.4, the retained complexity is 2 × 1.8 × 3.4^1.5 ≈ 2 × 1.8 × 6.24 ≈ 22.46, comprehensively reflecting both geometric and structural characteristics. The target grid number uses a nonlinear normalized mapping, and this nonlinear mapping avoids the uncontrolled grid number caused by extreme values ​​in the retained complexity. For example, when the retention complexity is 10, the number of target grids after mapping is 100 + (1000 - 100) × 10 / (10 + 1) ≈ 918; when the retention complexity is 0.5, the number of grids after mapping is 100 + 900 × 0.5 / 1.5 = 400, ensuring that the number of grids in blocks with different complexities is within a reasonable range.

[0089] It should be noted that the range of values ​​for the exponential weighting constants is as follows: the exponential weighting constants for density and elevation change rate range from 0.8 to 1.2, while the exponential weighting constants for structural risk factors range from 1.2 to 1.8. When the geometric features of the terrain data dominate, the exponential weighting constants for density and elevation change rate can be increased to 1.2; when the skeleton structure features are more critical, the exponential weighting constant for structural risk factors can be increased to 1.8. For example, in high-precision modeling scenarios, the exponential weighting constant for structural risk factors is set to 1.8 to strengthen skeleton protection. The specific calculation steps for nonlinear normalized mapping are as follows: first, calculate the ratio of the retention complexity to (retention complexity + 1) to obtain the normalization coefficient; then multiply this coefficient by (maximum number of triangles - minimum number of triangles); finally, add the minimum number of triangles, and round the result to the nearest integer. For example, if the retention complexity is 3, the minimum number of triangles is 100, the maximum is 1000, the normalization coefficient is 3 / 4 = 0.75, and the target grid number is 100 + 900 × 0.75 = 775. The minimum and maximum number of triangles are determined based on the area of ​​each terrain segment. The minimum number of triangles per square kilometer is 100, and the maximum is 1000. For every square kilometer increase or decrease in the segment area, the minimum and maximum number of triangles increase or decrease by 100 respectively. For example, for a segment of 0.5 square kilometers, the minimum number of triangles is 50, and the maximum is 500; for a segment of 2 square kilometers, the minimum is 200, and the maximum is 2000.

[0090] It should be noted that, based on the geometric features of terrain blocks (sampling point density, elevation change rate) and skeleton structure features (susceptibility to spatial entanglement of ridges and valleys), complexity is preserved through the product of three factors, and then nonlinear normalization is applied to map to a reasonable range of grid numbers to determine the target grid number. This design balances geometric density, terrain undulation, and skeleton entanglement risk. Its beneficial effect is that the target grid number accurately matches the actual complexity of the blocks, avoiding oversimplification or undersimplification caused by a single geometric feature. At the same time, the explicit calculation rules ensure the consistency and repeatability of the target grid number for different blocks, providing a scientific basis for subsequent triangulation folding optimization. This allows the simplification process to control computational load while protecting key skeleton structures, improving the overall quality and stability of the multi-resolution model.

[0091] In a preferred embodiment, triangulation folding optimization is performed in conjunction with the comprehensive error evaluation index to generate a multi-resolution terrain model, including:

[0092] Define candidate folding operations Comprehensive error evaluation index :

[0093]

[0094] in, For geometric deviation components, To maintain the components of the skeletal structure, Sensitivity to entanglement in ridge valley space. Here is the weighting constant. This is the sensitivity amplification factor;

[0095] The boundary constraint energy component is defined as follows:

[0096]

[0097] in, The distance from the vertex to the block boundary. The boundary influence scale is set to a value equal to the granularity of the basic analysis. ;

[0098] Iterative selection Perform a fold until a block is formed. Number of triangular facets Equal to the number of target grids ;

[0099] generate Layer of detail model , of which Number of target grids at each level Defined as:

[0100]

[0101] in, This represents the number of triangles in the original block.

[0102] It should be noted that the candidate folding operation is the vertex merging operation to be performed during the triangulation process, including edge folding or triangle folding. The comprehensive error evaluation index is a quantitative indicator measuring the impact of candidate folding operations on the terrain, composed of a weighted sum of geometric deviation components, skeleton structure preservation components, and boundary constraint energy components. The geometric deviation component is a quantitative indicator measuring the change in local terrain geometry caused by the candidate folding operation, obtained by summing the squared errors of horizontal coordinates and the weighted squared errors of vertical elevation. The skeleton structure preservation component is a quantitative indicator measuring the impact of candidate folding operations on the stability of the geomorphic skeleton network, obtained by multiplying the vertex skeleton influence weights by the changes in ridge distances before and after folding. The boundary constraint energy component is a quantitative indicator measuring the degree of damage to terrain block boundaries caused by candidate folding operations, obtained by multiplying the exponential decay term of the vertex-to-boundary distance by the square of the folding displacement. The ridge-valley spatial entanglement sensitivity is an indicator characterizing the responsiveness of the spatial proximity aggregation values ​​of valley lines and ridge lines within a terrain block to changes in analysis granularity. The weighting constant for the skeleton structure preservation component is a constant that adjusts the proportion of the skeleton structure preservation component in the overall error. A value of 0.5 is preferred, chosen to balance the protection priorities of geometry and skeleton structure, avoiding a single component dominating the evaluation. The weighting constant for the boundary constraint energy component is a constant that adjusts the proportion of the boundary constraint energy component in the overall error. A value of 0.3 is preferred, chosen to moderately strengthen boundary protection and prevent distortion of block boundaries due to folding. The sensitivity amplification factor is a coefficient that amplifies the influence of ridge-valley spatial entanglement sensitivity on the skeleton structure preservation component. A value of 2 is preferred, ensuring the sensitivity's discriminative power for different blocks. The distance from a vertex to a block boundary is the shortest straight-line distance from each vertex in the local vertex set to the edge of the terrain block. The boundary influence scale is a characteristic scale that controls the decay rate of the boundary constraint energy component. A value equal to the basic analysis granularity of the terrain block is preferred, chosen to match the block resolution, ensuring that boundary protection conforms to the terrain discretization accuracy. A new vertex is a single vertex generated after merging the local vertex set by candidate folding operations. An original vertex is the initial vertex to be merged in the local vertex set. The optimal folding operation is the candidate folding operation with the smallest comprehensive error evaluation index value. The number of triangle faces within a terrain tile is the total number of faces in the triangulation of the current terrain tile. The target mesh number is the number of triangle faces that the terrain tile needs to retain at a specific level of detail. The maximum level of detail is the total number of levels included in the multi-resolution model, preferably 5, chosen to balance multi-resolution requirements and computational efficiency, avoiding excessive resource consumption due to too many levels. The k-th level of detail model is the triangulation model of the terrain tile at the k-th level of detail. The target mesh number at the k-th level is the number of triangle faces that the terrain tile needs to retain at the k-th level of detail. The number of triangles in the original tile is the total number of initial triangle faces in the basic irregular triangulation of the terrain tile.

[0103] It should be noted that the comprehensive error evaluation index adopts a three-component design combining geometric deviation, DPSS weighted skeleton preservation, and boundary constraints. The skeleton preservation component is weighted by multiplying the sensitivity amplification factor with the ridge-valley spatial entanglement sensitivity, achieving dynamic protection for highly entangled blocks. For example, if a highly sensitive block has a ridge-valley spatial entanglement sensitivity of 1.5, a sensitivity amplification factor of 2, and a skeleton preservation component weight of 0.5, its weighted contribution is 0.5 × (1 + 2 × 1.5) = 2, higher than the 0.5 × (1 + 2 × 0.3) = 0.8 for a low-sensitivity block, thus providing targeted protection for the skeleton structure. The boundary constraint energy component adopts an exponential decay design; the closer the vertex is to the boundary, the greater the constraint weight, and the boundary influence scale is bound to the granularity of the basic analysis, ensuring that boundary protection is adapted to the block resolution. For example, if the basic analysis granularity of a certain block is 5 meters, the distance from the vertex near the boundary to the boundary is 1 meter, the attenuation term is exp(-1 / 5) = 0.818, and the constraint weight is relatively large; the distance from the vertex far from the boundary is 10 meters, the attenuation term is 0.135, and the constraint weight is relatively small, achieving differentiated boundary protection. The number of target meshes in the multi-resolution detail pyramid decreases exponentially to ensure smooth transitions between levels, and all levels use the same comprehensive error evaluation index to ensure consistency of models at different resolutions. For example, if the original block has 1600 triangles, the number of targets in the first layer is 400, and the number in the second layer is 100, the simplification ratio between levels is consistent to avoid abrupt morphological changes.

[0104] It should be noted that the weight constants for the skeleton structure preservation component and the boundary constraint energy component range as follows: the weight constant for the skeleton structure preservation component ranges from 0.3 to 0.7, and the weight constant for the boundary constraint energy component ranges from 0.2 to 0.4. The sum of the weights of the two and the geometric deviation component is 1. When the terrain skeleton structure is more critical, the skeleton weight is 0.7 and the boundary weight is 0.2; when the boundary integrity requirement is high, the skeleton weight is 0.3 and the boundary weight is 0.4. For example, in complex mountain modeling, the skeleton weight is 0.7 to strengthen the protection of the core structure. The specific calculation steps for the boundary constraint energy component are as follows: First, calculate the shortest distance from the vertex to the block boundary. Then, divide this distance by the boundary influence scale, take the negative value, and calculate the exponential function. Next, calculate the squared displacement of the new vertex compared to the original vertex. Finally, multiply the result of the exponential function by the squared displacement, and sum the products of all vertices to obtain the boundary constraint energy component. For example, if a vertex is 3 meters from the boundary, the boundary influence scale is 5 meters, and the squared displacement is 4, the contribution of this vertex is exp(-3 / 5)×4≈2.156. The optimal folding operation selection and update process is as follows: Initially, all possible candidate folding operations are generated. The comprehensive error evaluation index value of each operation is calculated. The operation with the smallest index value is selected for execution. After execution, candidate operations related to this operation are deleted, new candidate operations are generated, and the index value is calculated. This process is repeated until the number of triangles reaches the target mesh number. For example, if there are 100 candidate operations initially, after executing the optimal operation, 20 related operations are deleted, 15 new operations are generated, and the remaining 95 candidate operations enter the next round of selection. The calculation method for the target mesh number in the k-th layer is as follows: Multiply the number of triangles in the original block by 2 to the power of -2k, and round the result to the nearest integer. For example, if the original number of triangles is 2000, when k=1 in the first layer, the target number is 2000×2^(-2)=500; when k=2 in the second layer, the target number is 2000×2^(-4)=125.

[0105] It should be noted that the simplification of triangular meshes by folding must consider geometric shape, skeleton structure, and boundary integrity. A comprehensive error evaluation index incorporating these three components is constructed. The skeleton protection strength is dynamically adjusted based on the entanglement sensitivity of ridge and valley spaces. The optimal folding operation is selected to progressively simplify to the target mesh number, and then a multi-resolution detail pyramid is generated according to an exponentially decreasing rule. This design aligns with the core requirements of multi-resolution modeling. Its beneficial effect is that it balances geometric accuracy, skeleton integrity, and boundary stability during the simplification process, avoiding model distortion caused by single-dimensional optimization. A unified evaluation standard is used for different levels to ensure the consistency and smooth transition of the multi-resolution model, providing terrain models adapted to different accuracy requirements for subsequent applications. At the same time, the deterministic operation process improves the repeatability and reliability of the modeling method.

[0106] In a preferred embodiment, the boundary stitching optimization of adjacent terrain blocks is also included, as follows:

[0107] Divide adjacent terrain into blocks and Shared Boundary Parameterization ;

[0108] Set boundary resampling step size ,in These represent the basic analysis granularity for the two sub-blocks;

[0109] At parameter points Generate stitch boundary points :

[0110]

[0111] in, Two blocks are respectively in the parameters Interpolated coordinates of the original boundary points at the location;

[0112] The stitch weight is defined as:

[0113]

[0114]

[0115] in, The entanglement sensitivity of the ridge-valley space in the two segments are respectively. This is the sensitivity amplification factor.

[0116] It should be noted that adjacent terrain blocks are two geographically adjacent independent terrain processing units after quadtree recursive partitioning. The shared boundary curve is the common edge line between two adjacent terrain blocks. The parameterized variable is the variable used for unified parameterization of the shared boundary curve, with a value ranging from 0 to 1. The boundary resampling step size is a fixed interval for resampling the shared boundary curve, taking the smaller value among the basic analysis granularities of the two adjacent terrain blocks. The basic analysis granularity is a benchmark scale proportional to the average vertex spacing within the terrain block, used as a reference standard for unified spatial analysis. Parameter points are equidistant scattered points generated on the shared boundary curve according to the boundary resampling step size. The stitching boundary point is a common boundary point generated after fusing the coordinates of the original boundary points of two adjacent terrain blocks. The interpolated coordinates of the original boundary points are the coordinate values ​​obtained by interpolation at the shared boundary parameter point positions of the adjacent terrain blocks. The stitching weight is a value that adjusts the contribution of the two adjacent terrain blocks to the stitching boundary point, obtained through a linear transformation of their respective ridge-valley spatial entanglement sensitivity. Ridge-valley spatial entanglement sensitivity is an index characterizing the degree to which the spatial proximity aggregation value of valley lines and ridge lines within a terrain segment responds to changes in analysis granularity. The sensitivity amplification factor is a coefficient that amplifies the influence of ridge-valley spatial entanglement sensitivity on stitch weights; a value of 2 is preferred to ensure the sensitivity's discriminative power across different segments.

[0117] It should be noted that the stitching weights are designed based on the ridge-valley spatial entanglement sensitivity. A linear transformation is used to give higher-sensitivity blocks greater weight, rather than using equal weights or simple distance weighting. This design allows entanglement-sensitive blocks to dominate the morphology during boundary stitching, avoiding structural mismatches. For example, block one has a ridge-valley spatial entanglement sensitivity of 1.6, block two has 0.8, and the sensitivity amplification factor is 2. The stitching weights are 4.2 and 2.6 respectively. Block one has a higher proportion in the boundary morphology, protecting its sensitive skeletal structure. The shared boundary uses a unified parameterization based on arc length, rather than distance or angle parameterization, ensuring a one-to-one correspondence between the boundary points of the two blocks. For example, if the shared boundary curve length is 20 meters, the parameterized variables are uniformly distributed from 0 to 1. A parameter value of 0.25 corresponds to a position 5 meters from the starting point. Both block one and block two have corresponding interpolated coordinates at this position, ensuring stitching alignment. The boundary resampling step size is taken as the smaller value between the two block-based basic analysis granularities. This ensures that the resampling accuracy is adapted to the finer block resolution and avoids the loss of boundary details due to an excessively large sampling step size. For example, if the basic analysis granularity of block one is 6 meters and that of block two is 4 meters, the resampling step size is set to 4 meters to match the finer resolution of block two.

[0118] It should be noted that the specific algorithm for the unified parameterization of the shared boundary curve adopts the natural parameterization method. Using the arc length of the shared boundary curve as a benchmark, the arc length from the start point to the end point is normalized into a parameterized variable of 0 to 1, with the parameter value proportional to the arc length. For example, if the total arc length of the shared boundary curve is 25 meters, a parameter value of 0.4 corresponds to a position 10 meters from the start point, and a parameter value of 0.8 corresponds to a position 20 meters from the start point. The calculation method for the interpolated coordinates of the original boundary points uses cubic spline interpolation. Using the original vertex coordinates on the adjacent terrain block boundaries as control points, the interpolation function is solved at the set parameter point positions to obtain the corresponding coordinate values. For example, if the original vertex coordinates on the block boundary are (10,20,5), (20,30,6), and (30,40,7), the interpolated coordinates at the parameter value of 0.3 are (16,26,5.5) obtained through cubic spline interpolation. The number of parameter points generated is determined by the ratio of the total length of the shared boundary curve to the boundary resampling step size. If the ratio is not an integer, it is rounded up to ensure that the entire shared boundary is covered. For example, if the length of the shared boundary curve is 18 meters and the boundary resampling step size is 4 meters, the ratio is 4.5, which is rounded up to 5 parameter points with a spacing of 4 meters. The last parameter point covers the remaining 2 meters.

[0119] It should be noted that, based on the differences in the entanglement features of the skeletons of adjacent terrain blocks, the shared boundary is parameterized with a unified arc length. A finer granularity of basic analysis is used as the resampling step size, and stitching weights are constructed by combining the spatial entanglement sensitivity of ridges and valleys. This generates stitching boundary points that fuse the features of both sides, which are then embedded into the mesh through constrained triangulation. This design aligns with the boundary integrity requirements of multi-resolution modeling. Its beneficial effect is that it allows ridge-valley entanglement-sensitive blocks to dominate the morphology during boundary stitching, reducing structural mismatches. Unified parameterization and cubic spline interpolation ensure smooth and continuous boundaries, while constrained triangulation achieves seamless embedding. This forms a closed loop with the aforementioned sensitivity calculation and triangulation folding optimization, improving the boundary consistency and structural integrity of the entire multi-resolution model.

[0120] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.

Claims

1. A multi-resolution modeling method for complex terrain based on LOD and triangular mesh folding, characterized in that, include: A basic irregular triangular network is constructed and a geomorphic skeleton network is extracted. The geomorphic skeleton network includes a set of valley lines based on the classification of water catchment parameters and a set of ridge lines based on the classification of watershed diffusion parameters. Constructing a basic irregular triangular network and extracting the topographic skeleton network includes: Perform Delaunay triangulation on the terrain sampling point set on a two-dimensional plane projection to generate a basic irregular triangular network, and assign an elevation value to each vertex; The flow accumulation of each vertex in the basic irregular triangular network is calculated based on flow direction analysis. The flow accumulation is the weighted sum of the area represented by the vertex itself and the flow accumulation of all upstream adjacent vertices whose flow direction points to the vertex. The diffusion accumulation of each vertex in the basic irregular triangular network is calculated based on the reverse flow analysis. The diffusion accumulation is the weighted sum of the area represented by the vertex itself and the diffusion accumulation of all adjacent vertices flowing downstream from the vertex. In the basic irregular triangular network, connected paths with high confluence accumulation are extracted into valley line sets, and connected paths with high diffusion accumulation are extracted into ridge line sets. Based on the sampling point distribution density and the rate of change of terrain elevation, the basic irregular triangular network is recursively divided into quadtrees to obtain multiple terrain blocks; For each terrain block, calculate the ridge-valley spatial entanglement sensitivity, which characterizes the degree of response of the spatial proximity aggregation value from the valley line set to the ridge line set within the block to the change of spatial analysis granularity. A simplified comprehensive error evaluation index for triangulation is constructed. The comprehensive error evaluation index includes a geometric deviation component and a skeleton structure preservation component weighted by the spatial entanglement sensitivity of the ridge valley. The skeleton structure preservation component is used to quantify the degree of deviation of the simplification operation from the distance field from the vertex to the terrain skeleton network. The retention complexity of terrain blocks is calculated based on the ridge and valley spatial entanglement sensitivity to determine the number of target grids, and triangular mesh folding optimization is performed in combination with the comprehensive error evaluation index to generate a multi-resolution terrain model.

2. The method for multi-resolution modeling of complex terrain based on LOD and triangular mesh folding according to claim 1, characterized in that, Based on the sampling point distribution density and the rate of change of terrain elevation, the basic irregular triangular network is recursively partitioned using a quadtree to obtain multiple terrain blocks, including: Starting from the global root node, calculate the partition evaluation value of the current terrain patch; The zoning evaluation value is composed of a weighted sum of the normalized sampling point distribution density and the normalized topographic elevation change rate; Determine whether the partition evaluation value is less than a preset partition stopping threshold; If the partition evaluation value is greater than the partition stopping threshold, the current terrain block is divided into four sub-blocks and the calculation is performed recursively. If the partition evaluation value is less than or equal to the partition stopping threshold, the partitioning is stopped and the current block is taken as the final terrain block.

3. The method for multi-resolution modeling of complex terrain based on LOD and triangular mesh folding according to claim 1, characterized in that, For each terrain patch, the spatial entanglement sensitivity of the ridge-valley line is calculated, which characterizes the degree of response of the spatial proximity aggregation value from the valley line set to the ridge line set within the patch to changes in spatial analysis granularity. This includes: Extract the valley line set subset and ridge line set subset within the current terrain segment; Discrete sampling points are extracted on the valley line subset with a fixed arc length step, and the Euclidean distance from each discrete sampling point to the ridge line subset is calculated. The spatial proximity aggregation value of the current terrain block is calculated based on the exponential decay model. The spatial proximity aggregation value is the integral approximation of the ratio of the negative distance of all discrete sampling points to the spatial analysis granularity. Determine the basic analysis granularity of the current terrain segment, wherein the basic analysis granularity is proportional to the average vertex spacing within the terrain segment; The ridge valley spatial entanglement sensitivity is calculated as the difference between the logarithm of the spatial neighbor aggregation value at the scaled spatial analysis granularity and the logarithm of the spatial neighbor aggregation value at the basic analysis granularity, divided by the logarithm of the granularity scaling ratio.

4. The method for multi-resolution modeling of complex terrain based on LOD and triangular mesh folding according to claim 1, characterized in that, Geometric deviation components include: For each candidate folding operation, identify the set of local vertices to be merged by the operation and the new vertices generated after the folding; Calculate the sum of the squared error in the horizontal coordinate direction and the weighted squared error in the vertical elevation direction for each vertex in the local vertex set; The sum is used as a geometric deviation component to measure the change in local geometry caused by the folding operation, wherein the error weighting coefficient in the vertical elevation direction is greater than the error weighting coefficient in the horizontal coordinate direction, so as to maintain the elevation characteristics of the terrain.

5. The method for multi-resolution modeling of complex terrain based on LOD and triangular mesh folding according to claim 1, characterized in that, The skeleton structure retains components, including: For each candidate folding operation, calculate the shortest distance from each vertex in the local vertex set involved to the ridge set and the valley set; The skeleton influence weight of each vertex is constructed based on the exponential decay function. The skeleton influence weight is composed of the exponential decay term of the shortest distance from the vertex to the ridge line set and the exponential decay term of the shortest distance to the valley line set. The characteristic scale of decay is consistent with the granularity of the basic analysis of terrain segmentation. Calculate the absolute change in the shortest distance from the vertices to the set of ridges before and after folding; The sum of the product of the skeleton influence weight and the absolute change over the local vertex set is determined as the skeleton structure preservation component that measures the impact of the folding operation on the stability of the terrain skeleton network.

6. The method for multi-resolution modeling of complex terrain based on LOD and triangular mesh folding according to claim 1, characterized in that, The retention complexity of terrain patches is calculated based on the ridge-valley spatial entanglement sensitivity to determine the target grid number, including: The retention complexity of the terrain patch is calculated, which is composed of the product of three factors: the first factor is a power function of the ratio of the distribution density of the terrain patch sampling points to the global reference value; the second factor is a power function of the ratio of the terrain elevation change rate of the terrain patch to the global reference value; and the third factor is a power function of the structural risk factor enhanced by the spatial entanglement sensitivity of the ridge and valley. Construct a mapping function for the target number of grid cells, and transform the retained complexity to a preset range of grid cell numbers through a nonlinear normalization mapping; The number of target grids for the current terrain patch at a specific level of detail is calculated based on the mapping function, and the result is rounded.

7. The method for multi-resolution modeling of complex terrain based on LOD and triangular mesh folding according to claim 1, characterized in that, Combining the aforementioned comprehensive error evaluation index, triangulation folding optimization is performed to generate a multi-resolution terrain model, including: For each terrain block, a set of candidate folding operations is established, and the comprehensive error evaluation index value of each candidate folding operation is calculated. The comprehensive error evaluation index value includes not only the geometric deviation component and the skeleton structure preservation component, but also the boundary constraint energy component used to prevent folding from destroying the block boundary. Perform a folding operation with the minimum comprehensive error evaluation index value, and update the associated candidate folding operations and their evaluation index values. Repeat this process until the number of triangles within the terrain tile reaches the target number of grids. A multi-resolution detail pyramid is constructed by setting the number of target grids at each level that decreases geometrically with the level index, and by performing a folding sequence using the same comprehensive error evaluation index, a series of terrain patch grid models from the original resolution to the lowest resolution are generated.

8. The method for multi-resolution modeling of complex terrain based on LOD and triangular mesh folding according to claim 1, characterized in that, This also includes boundary stitching optimization for adjacent terrain blocks, as follows: Extract the shared boundary curve of two adjacent terrain blocks and perform unified parameterization on the shared boundary curve based on arc length; Determine the boundary resampling step size, which is the smaller value among the basic analysis granularities of two adjacent terrain blocks; The stitching weights of two adjacent terrain blocks on a shared boundary are calculated. The stitching weights are obtained by linear transformation of the ridge-valley spatial entanglement sensitivity of each terrain block, so that the block with higher ridge-valley spatial entanglement sensitivity has a larger weight. Generate stitch boundary point coordinates, wherein the stitch boundary point coordinates are the weighted average of the boundary point coordinates of two adjacent terrain blocks at the same parameter position, and the weight is the stitch weight; By using constrained triangulation, the generated stitching boundary points are re-embedded into the mesh of two adjacent terrain blocks as common boundaries.

Images