First of all, the origin and first base plane should be determined. Assuming a set of discrete points n and a vector set, V _i is composed of vector lines that connect all the discrete points to the origin (refer to Figure 1). P(i) ( 1 ≤ i ≤ n ) can be considered as an origin. Then, the product of each vector pair is calculated. The pair with maximal absolute values are designated as vectors OA and OB, where O is the origin, A is the initial anchor point, B is the initial floating-point, and OAB is designated as the first base plane.

The point dataset should be sequenced next. The random 3D discrete points are sequenced according to the spatial proximity principle. The initial anchor point A is the first point in the point column, and the following new point nearest to the current point must be searched. The initial floating-point B is taken as the final point of the point column, and O is not included in the point column.

Finally, the feature points should be selected. The distance between points and the first base plane is calculated. When the maximum distance is less than the specified threshold, all points must be deleted. When the maximum distance exceeds the specified threshold, the sequenced point dataset can be divided into two-point datasets for the threshold, and then the above processes are repeated for each data set. (Fei et al., 2006).

There is a natural coupling relationship between rivers and landforms (Shu, 2012). In simplifying river vector data and DEM based on the 3D_DP algorithm, the curved shape of the rivers in the horizontal direction cannot be retained. Particularly, the points of the river with different bending sizes are eliminated under the control of a certain threshold in the local scale with flat areas.

The BAI is proposed to improve the 3D_DP algorithm; it can identify the bending size of the river shape. BAI is used to adjust the distance between point and plane in the 3D_DP algorithm, and the pseudo distance between point and plane is obtained. The pseudo distance between point and plane is compared with the preset threshold value to determine the selection of river points afterward. Finally, the river and DEM’s overall shape and main features are preserved.

The merged point dataset was sequenced based on the improved 3D_DP algorithm according to the importance of each point. The aim of sequencing is to obtain a point dataset to support the pre-processing of LOD modelling. In the improved 3D_DP algorithm, the threshold value and the BAI determine the selected point’s importance level. When the threshold value in the 3D_DP algorithm is set to zero, we regard all the points as a cut-off point involving simplification. The reorganization of the dataset can be accomplished based on the sequencing order of the point. The importance of the point was then assigned to its attribute.

Furthermore, the river points are sequenced according to the pseudo distance between the point-and plane with the bending adjustment index. Some specific points have to be retained, such as the points derived from boundary extraction. Since these points do not participate in the subsequent simplification process, they are not required to undergo the sequencing process. The attribute (importance) of the points is defined as the most critical points. The intersections of river endpoints and tributaries are retained, but they can be deleted at a selected level if necessary.

The simplification of the LOD model is to obtain a simplified model with minimal deviations. In this algorithm, the point was deleted on the basis of the descending order of the importance. An error evaluation method is then employed to quantitatively assess the error arising from the simplified model (Li et al., 2007).

Hausdorff distance is widely applied to measure the degree of similarity between two datasets. It can be used to represent and evaluate the magnitude of the simplification error of the LOD model (Cao, 2013). Here and in what follows, we describe the process of error evaluation.

Given two point datasets: M = m 1 , . . . m p is the original model, N = n 1 , . . . n p is the simplified model. The N is dynamically composed of point sets in which points were deleted on the basis of the descending order of the importance.

The Hausdorff distance d _E represents the distance between m, an element of M, to N, given by Equation 3: d E m , N = min n ∈ N d m , n (3) where, d m , n is the Euclidean distance from point m to point n (See Equation 4): d m , n = m x − n x 2 + m y − n y 2 + m z − n z 2 (4)

The absolute value of the height difference between point m and point n can be calculated as Equation 5 h m , n = H m − H n (5)

The vertical spatial error h _s based on the one-way Hausdorff distance is defined as Equation 6: h S M , N = a v g m ∈ M h m , N (6) where h _s is the vertical spatial error of the model derived from the simplification process. When h _s is less than the threshold value derived from the screen projection error, the new dataset (N) can replace the original dataset (M) to conclude the simplification process.

The vertical spatial error is computed during the data preprocessing stage to reduce the computational load in the real-time display stage.

To avoid the occurrence of spatial conflicts such as “river climbing” (Shu, 2012), we treat the river line segment as a mandatory constraint to generate a Delaunay network, which is constructed to achieve unified LOD modelling of the river elements and DEM. The points are deleted from the original river-constrained Delaunay triangulation network based on the order of importance. Therefore, the relevant triangles and the river line segments as constraint edges are deleted for each deleted vertex. Geometric gaps derived from point deletion have to be triangulated locally. Based on the ear-cutting algorithm proposed by Kong, the algorithm adds the judgment of constrained edges to complete the local triangulation operation of polygons (Wibowo et al., 2019; Zhai et al., 2010).

The polygon holes that need to be triangulated are convex polygons or concave polygons. Convex polygons are more accessible to implement than concave polygons. Taking the concave polygon as an example, we present the polygon triangulation algorithm designed and used in this algorithm. As shown in Figure 4, a concave heptagon (seven-sided polygon) consists of vertices A, B, C, D, E, F, and G, and CX and XA are river line segments with constraint edges, then the point X has to be deleted.

In this algorithm, the diagonal of the concave polygon refers to the line of any two non-adjacent vertices, which does not intersect with polygon’s edges and is located within the polygon. For example, AF is not considered as a diagonal because it is not located within the polygon. The workflow of the algorithm is as follows:

(1) To begin with, we determine whether constraint edges for river segments exist. When a constrained river line segment exists, the line segment AC has to be connected after deleting the point X. The line AC acts as a constrained edge, dividing the polygon into two polygons, ABC and CDEFGA. When there is no river segment in the polygon, the procedure of the algorithm is executed directly.

The vertical spatial error from the simplification can be displayed on the two-dimensional screen. And the error derived from projecting the image into a 2D perspective is called screen projection error. The calculation of the screen projection error is directly related to the line of sight direction, the distance of view, and the angle of intersection. The accurate estimate demands enormous computational resources and remarkably reduces the algorithm’s efficiency. Therefore, the algorithm proposed and adopted a simplified approach (Tang and Dou, 2023).

As shown in Figure 5, taking the vertical axis projected onto the screen (elevation direction) as an example, the model spatial error of the model simplified by the calculated in the elevation direction can be directly used as the vertical axis direction error (OE in Figure 5) to participate in the subsequent calculation. When the distance between the view plane is certain and the center line of sight is horizontal, the size of the model space error OE (corresponding to 3.2 model space error h S ) on the projection plane after perspective projection will generate the largest projection error on the screen.

This maximum projection error can be used as the judgment condition, the model space error h S (the model space error of elevation direction calculated in 3.2) is directly used to participate in the calculation, which can simplify the calculation and improve the efficiency. Based on the properties of similar triangles, the calculation of the corresponding relationship is as Equation 7: B C = f S × O E (7)

The screen projection error BC is equal to δ and the model spatial error O E is equal to h S . The calculation of δ is as Equation 8: δ = f S × h S (8) where δ is screen projection error; f is the distance from the observation point to the nearsighted surface, namely the focal length); S is the distance from the observation point to the object; and, O E is the model space error, which corresponds to h S calculated in Section 3.2.

From Equation 8, the three parameters related to screening projection error are the focal length f, the model spatial error O E , and the distance S. Obviously, the value of f is given and the value of S can be determined by the position of the viewpoint in the real-time display stage. Once the tolerance of the screen projection error δ is known, we may deduce the tolerance of the model spatial error, which can then be used as a control parameter to simplify the LOD modelling. Different line-of-sight conditions correspond to different spatial error limits of models, which also fit LOD models with different accuracy requirements.

Assume that the number of pixels is taken as the threshold value of screen projection error,but the δ is in unit of distance in Equation 8, so the pixels have to be converted into unit of distance. Suppose P is the number of pixels in the projector plane occupied by BC (Figure 5), and H is the number of pixels occupied by the screen (AM in Figure 5). The corresponding relationship is illustrated as Equation 9: δ A C = P H / 2 (9)

From Equations 8, 9, the following relationship can be inferred in Equation 10: δ = f S × h S = 2 P f ⁡ tan α 2 H (10)

The model space error h S can be derived from Formula 10. The number of pixel P can be designated as the limitation of the screen projection error according to the actual visual distance, and then the limit of the model space error h _lim can be calculated by hs (see Equation 11). h lim = 2 P S ⁡ tan α 2 H (11) where S is the line-of-sight, a is the field-of-view angle, P is the pixel number, and H is the pixel resolution of the screen projection plane. In this paper, a single pixel or several pixels can be specified as the threshold for the screen projection error.

In the real-time dynamic display stage of the model, the visual distance between the viewpoint and the center of each block can be obtained in real time, and the h _lim can be deduced in real time according to a specified value of the screen projection error (namely,the number of pixel P). Then, according to the deletion sequence of the arranged points of each block, the model spatial error h S is compared with h _lim . When h S ≤ h _lim , all the points in front of this point are deleted. When h S > h _lim , this point and those that follow are preserved.

The displayed number of points, river segments, and triangles can be determined based on the above point dataset which is generated on the basis of h S and h _lim . As the model is simplified in the data preprocessing stage, the changes in the topological structure of any point, river segment, and triangle caused by point deletion have been recorded by river variables and triangle variables, which were defined in the data structure.

The algorithm is realized in C# programming and OpenGL rendering (Zeng and Chen, 2018; Richard and Nicholas, 2010). The computer hardware used in the experiment is configured with Intel Core i5 CPU, 2.5 GHz frequency, 8G memory, and 2G graphics card, and the software environment was Win7 operating system. In this experiment, we used a display device with a screen resolution of 1920*1,080 for rendering.

The original DEM is derived from ASTER GDEM with a 1-arc second spatial resolution, WGS84 geographic coordinate, and GeoTIFF format. Vector dataset are selected from the vector dataset of river network (GADM V1, http://www.gadm.org/) of Baoding and its surrounding areas, a region which is located in Hebei Province, China, with about 220 × 220 km². The topography of the study site extends from northwest to southeast with various landforms including flat, rugged, and hilly terrains. The northwest part is dominated by mountains, and the southeast region is mainly composed of plains. The elevation ranges from 8 to 2,893 m. In this study, we eliminated scattered and tiny rivers while retaining five complete tree-like river systems, consisting of 131 river segments. A total of 13,997 DEM points and 3,911 river points are extracted from the original DEM and river network dataset, respectively. Figure 6 is the interface of the experimental prototype system written by the authors.

During the real-time rendering of the data, we set a limitation value (The number of pixel) for the screen projection error to meet visual requirements. That is to say, when deleting a specified set of points under a certain sight distance, the screen projection error of model must be less than or equal to the specified number of pixel, then the simplified model can be used to replace the original model.