A high-fidelity multiresolution DEM model for Earth systems

The impact of topography on Earth systems variability is well recognised. As numerical simulations evolved to incorporate broader scales and finer processes, accurately assimilating or transforming the topography to produce more exact land-atmosphere-ocean interactions, has proven to be quite challenging. Numerical schemes of Earth 10 systems often use empirical parameterisation at sub-grid scale with downscaling to express topographic endogenous processes, or rely on insecure point interpolation to induce topographic forcing, which creates bias and input uncertainties. DEM (digital elevation model) generalisation provides more sophisticated systematic topographic transformation, but existing methods are often difficult to be incorporated because of unwarranted grid quality. Meanwhile, approaches over discrete sets often employ heuristic approximating which are generally not best 15 performed. Based on DEM generalisation, this article proposes a high-fidelity multiresolution DEM model with guaranteed grid quality for Earth systems. The generalised DEM surface is initially approximated as a triangulated irregular network (TIN) via selected feature points and possible input features. The TIN surface is then optimised through an energy-minimised centroidal Voronoi tessellation (CVT). By devising a robust discrete curvature as density function and exact geometry clipping as energy reference, the developed curvature CVT (cCVT) converges, the 20 generalised surface evolves to a further approximation to the original DEM surface, and the points with the dual triangles become spatially equalised with the curvature distribution, exhibiting a quasi-uniform high quality and an adaptive variable-resolution. The cCVT model was then evaluated on real LiDAR-derived DEM datasets and compared to the classical heuristic model. The experimental results show that the cCVT multiresolution model outperforms classical heuristic DEM generalisations in terms of both surface approximation precision and surface morphology 25 retention.

should be considered.(Now we realized that the curvature of the Earth itself might not be of so concerned).We modified this sentence to follow its initial intention as: "Furthermore, considering the topography over a wider range may require re-implementing cCVT on the geographical coordinate base instead of on the currently used projection coordinate base."

[List of Changes]
1. P2L18."Although there are many situations where dynamic conditions are stressed as stronger impacts on predictions."According to Comment #1 2. P12L30."Evaluation of the cCVT multiresolution DEM model on Earth and environmental systems over wide-ranged domains and scales is a topic for future studies."According to Comment#2 3. P12L31."Furthermore, considering the topography over a wider range may require re-implementing cCVT on the geographical coordinate base instead of on the currently used projection coordinate base.""According to Comment#3 1 Introduction

Topography in Earth systems
Topography is one of the main factors controlling processes operating at or near the surface layer of the planet (Florinsky and Pankratov, 2015;Wilson and Gallant, 2000).With the success of Earth and environment systems with these scale diversified processes, there exist persistent demands for extending their utility to new and expanding scopes (Ringler et al., 2008;Tarolli, 2014;Wilson, 2012), as exemplified by lapse-rate controlled functional plant distributions (Ke et al., 2012), orographic forcing imposed on oceanic and atmospheric dynamics (Nunalee et al., 2015;Brioude et al., 2012;Hughes et al., 2015), topographic dominated flood inundations (Bilskie et al., 2015;Hunter et al., 2007), and many other geomorphological (Wilson, 2012), soil (Florinsky and Pankratov, 2015), and ecological (Leempoel et al., 2015) examples from Earth systems.However, as numerical simulation systems evolved to incorporate broader scales and finer processes to produce more exact predictions (Ringler et al., 2011;Weller et al., 2016;Wilson, 2012;Zarzycki et al., 2014), how to accurately assimilate or transform the fine resolution topography has proven to be a quite difficult task (Bilskie et al., 2015;Chen et al., 2015;Tarolli, 2014).Earth and environmental simulations usually adopt sub-grid schemes to express topography heterogeneous processes (Fiddes and Gruber, 2014;Kumar et al., 2012;Wilby and Wigley, 1997).The sub-grid schemes are designed for the empirical parameterisation rather than accurate topography representation, which often leads to mixed-up uncertainties and bias of endogenous variability (Jiménez and Dudhia, 2013;Nunalee et al., 2015).However, underresolved representation could be improved by variable-resolution enhancement, and bias of simulations can be justified by more fidelity topography transformation (Nunalee et al., 2015;Ringler et al., 2011).Topography is also commonly treated as a static boundary layer in dynamics simulations, where different interpolation strategies and mesh refinement techniques are used to convey terrain variation (Guba et al., 2014;Kesserwani and Liang, 2012;Nikolos and Delis, 2009;Weller et al., 2016).But a mesh constructed from interpolated vertices does not necessarily comply with the terrain relief, and elevation errors are frequently reported as an input uncertainty (Bilskie and Hagen, 2013;Hunter et al., 2007;Nunalee et al., 2015;Wilson and Gallant, 2000).Although there are many situations where dynamic conditions are stressed as stronger impacts on predictions (Cea and Bladé;Budd et al., 2015), the underlying topography is still very important due to its increasingly improved fidelity to the Earth's surface (Bates, 2012;Tarolli, 2014), and a sophisticated topography transformation would be beneficial to reduce discrepancies arisen from physical inconsistencies (Chen et al., 2015;Glover, 1999;Ringler et al., 2011).

Multiresolution DEM model
Systematic scale transformation of topographic data has long been studied under terrain generalisation, where precise surface approximation and terrain structural feature retention have both been pursued (Ai and Li, 2010;Chen et al., 2015;Guilbert et al., 2014;Jenny et al., 2011;Weibel, 1992;Zhou and Chen, 2011).Triangulated irregular networks (TIN) are generally chosen as a substitute for the regularly spaced grids (RSG), and terrain feature points (critical points or salient points from some significance metric) are selected for constructing the network.Triangular networks are used for their adaptiveness to locally enhanced multiresolution schemes.Critical points or salient points are selected because they can effectively improve the approximation precision (Heckbert and Garland, 1997;Zakšek and Podobnikar, 2005;Zhou and Chen, 2011).As surface approximation precision and terrain feature retention are competitive for the redistribution of feature points, DEM (digital elevation model) generalisation is differentiated from terrain generalisation by its emphasis on surface approximation as a whole, with the aim of providing precise surface interpolation (Guilbert et al., 2014).Terrain generalisation emphasises geomorphology or landform depiction, where map generalisation measures (such as abstracting, smoothing) are drawn to produce progressive data reduction, with the effect that the main relief features are strongly stressed while non-structural details are massively suppressed (Ai and Li, 2010;Guilbert et al., 2014;Jenny et al., 2011).Since the static topographic layers are commonly composed directly from DEM datasets for diverse simulation interests, maintaining precise surface approximation for rigorous boundary conditions is often more important than 'sparse' geomorphology representation.While DEM datasets are usually used interchangeably with topography or terrain in Earth systems, we will use DEMs and topography indiscriminately hereafter.Existing DEM generalisations can be catalogued into two broad classes, namely, heuristic refinements and smoothfittings, according to differences in surface approximation strategy.The first class of approaches is due to the computational feasibility consideration, for selecting a TIN surface that best approximates the original DEM surface from exhaustively enumerating (of triangular combinations) requires exponential time (Chen and Li, 2012;Heckbert and Garland, 1997).It thus forces existing researches to employ some heuristic strategy, in which insertion (or deletion) refinements on feature points are adopted, to find a sub-optimal approximation that is computationally practical.In each insertion (or deletion), rearranging the entire existing grid to obtain a better approximation is also computationally prohibitive and thus not adopted (Chen et al., 2015;Heckbert and Garland, 1997;Lee, 1991); this may result in the clustering of feature points (Chen and Li, 2012).Among those existing heuristic approaches, trenching the pre-extracted terrain features (drainage streamlines, for example) into the TIN surface seems quite appealing (compound method) (Chen and Zhou, 2012;Zhou and Chen, 2011), but the quality of the generalised TIN surface cannot be guaranteed, and the existence of sliver triangles makes it difficult to be incorporated with simulation numerical stability (Kim et al., 2014;Weller et al., 2009).The second class of approaches recognised the TIN surface constructed from locally computed feature points as not-well approximation to the original DEM surface.Much research thus considered global approximation instead of relying on an elaborate feature point selection scheme, such as bi-linear, bi-quadric, multiquadric, Kriging, or general radial base function-based fitting (Aguilar et al., 2005;Chen et al., 2015;Schneider, 2005;Shi et al., 2005).The proposed multi-quadric method (Chen et al., 2012;Chen et al., 2015), for example, well approximates the original DEM surface with a high-order smooth surface, and the smooth surface provides a kind of rejection mechanism to cure the feature point clustering problem.However, the high-order radial base function is computationally expensive when a broad scenario is involved (Chen et al., 2015;Mitášová and Hofierka, 1993).In brief, existing DEM transformations are neither well performing with respect to loyalty to the original terrain surface nor easily incorporated by the numerical schemes.

Aims and contributions
The purpose of this article is to devise a multiresolution DEM model that optimises surface approximation and guarantees grid quality that can be easily incorporated into the simulation systems.Multiresolution is an effective paradigm to model scale diversity (Du et al., 2010;Guba et al., 2014;Ringler et al., 2011;Weller et al., 2016).Amongst a number of promising approaches, we are especially fascinated by the centroidal Voronoi tessellations (CVTs) method as an intuitive way to redistribute samples with a designated function (Du et al., 1999;Du et al., 2010;Ringler et al., 2011) to develop an optimised surface transformation method to realise multiresolution terrain models.CVT is essentially a two-step optimisation loop, i.e., spatial domain equalisation from Voronoi tessellation and property domain equalisation from barycentre computation (Du et al., 1999).To make this general optimisation method work for DEM transformation, we made the following contributions: (1) The generalised DEM surface is initially approximated by a triangular grid constructed from selected feature points.
The selection of feature points have important morphological structures embedded (in the form of serialised points sequences), for computed (such as D8 flow algorithm) or auxiliary input morphological lines has been proved to have significant influence on the quality of the transformed DEM surface (Zhou and Chen, 2011).The proposed method keeps the structural lines in the optimisation loop and makes it different to existing CVT implementations where stationary points are commonly not considered.
(2) For the discrete TIN surface, we compute robust mean curvature on each facet.The attached curvature acts as a frequency distribution.In this discrete spatial domain and frequency domain, the CVT loops and makes sample facets equalised from both domains.Spatial equalisation warrants a quasi-uniform grid quality, while the curvature domain equalisation warrants adaptive distribution conforming to the terrain relief.It is thus a totally different DEM generalisation approach, and we called it curvature CVT (cCVT).
(3) Existing CVT implementations often undertake a clustering approach.However, clustering over discrete sets suffers from numerical issues such as zigzag boundaries, invalid cluster cells (Valette et al., 2008) and limited grid quality.
By devising an exact geometry clipping technique, this article develops a dedicated CVT algorithm for DEM transformation that helps to improve or avoid the numeric problems listed above.
The cCVT works on discrete sets but has a global optimisation mechanism.It promises an optimised surface approximation and quality grid, which can be used to build a high-fidelity multiresolution terrain model.From this terrain model, reliable surface variables can be estimated under a coupled system, or improved computational mesh can be constructed and refined to possible dynamic conditions.

Organisation of the article
The rest of the article is organised as follows.In Section 2, the theory behind CVT for optimised DEM surfaces is introduced, techniques for incorporating DEM generalisation principles and fast convergence are presented, and the differences between the cCVT implementation and classical clustering approach are discussed.In Section 3, the cCVT model is tested with real LiDAR-derived terrain datasets.Section 4 discusses some considerations, comparable results, possible causes, and interpretations of the cCVT model.Finally, Section 5 briefly presents a short conclusion and outlook.

Definition
Centroidal Voronoi tessellation is a space tessellation for each Voronoi cell's geometrical centre (in the spatial domain) that coincides with its barycentre from the abstract property domain (Du et al., 1999).Here, the property domain is analogous to the frequency domain.For the vertex set t {  }  1 in Ω ⊂  3 , the Voronoi tessellation graph is defined as That is, a Voronoi cell   is the set of points whose distance to   is less than that to any other vertices.| • | is the Euclidean norm.Every vertex and its corresponding dual cell commonly have some intensity scalar  attached from some abstract property domain, which is called a density function.The total potential energy of the Voronoi graph V of a terrain surface can be computed by summing up the potential energies of all cells   : . (2) The energy minimiser is which minimises the surface's total potential, since   ̅ = (̅ ,  ̅, ̅ ) satisfies In other words, when   coincides with the barycentre  ̅  , each cell's potential effect on the property domain (gravity) becomes equalised to a stable energy state.

Lloyd Relaxation
The most classical energy minimisation process of centroidal Voronoi tessellation is expressed by Lloyd's Relaxation (Lloyd, 1982).The main idea of this algorithm is to first tessellate the surface, and then perform density integration over the area to find a 'gravity' barycentre for each tessellated cell, which is used as the new site for the iteration.The pseudo code of this procedure is shown below.
We follow Lloyd's elegant idea.The barycentre of a 2-dimensional Voronoi cell may fall outside this surface patch, so an additional calculation may be needed.Du et al. suggested projecting the barycentre onto a nearest facet and using instead the constrained projection point for the new site (Du et al., 2003).Others suggest quadric interpolations over all the facets of the cluster for more accurate site calculations (Valette et al., 2008).

Clustering CVTs
Lloyd's Relaxation requires Voronoi tessellation on a discrete 2-dimensional surface, but direct Voronoi tessellation on a piecewise smooth surface requires costly geodesic computation and may be challenged by complicated numerical issues (Cabello et al., 2009;Kimmel and Sethian, 1998).Du et al. suggested that CVT could be realised through some clustering approaches (Cohen-Steiner et al., 2004;Du et al., 1999;Du et al., 2003), that is, using the associated property as density function to cluster facets and then find the clustered cells' barycentres to create new clustering sites.Through this heuristic iteration, the new sites along with the new tessellations compose better and better approximation to the original surface, with their spatial distribution conforming to the pre-defined density function.The clustering approach avoids geodesic tessellation by direct facet combination, which is computationally light.The greatest expenditure then comes from global distance computation for identifying every cell to its cluster centre.However, this k-means like clustering over discrete facets suffers from some numeric issues, such as zigzag cluster cell boundaries -since no geodesic Voronoi tessellation was used, and invalid clusters due to disconnected set of facets (Valette and Chassery, 2004;Valette et al., 2008).Furthermore, the key to the quick clustering algorithm is that it avoids generating new sites (to avoid surface reconstruction) and relies on existing sites (or facets).Thus, the generated grid cells may not be well qualified.

Curvature as density function
Terrain surface critical points such as peaks, pits, and saddles are treated as gravity equilibria and key elements depicting the surface geometry in the large (Banchoff, 1967;Milnor, 1963); a further extension of the critical points on a secondorder surface derivatives will describe a more fundamental terrain geometry shape (Jenny et al., 2011;Kennelly, 2008).When constructing a generalised DEM surface, these feature points are commonly used as a base set, and additional input points, or pre-extracted terrain structures are embedded for further approximation (Guilbert et al., 2014;Zakšek and Podobnikar, 2005;Zhou and Chen, 2011).The additional input points or pre-extracted terrain structures of interest are also commonly required in numerical simulation setups for cross-checking or validation purposes.Based on these observations, and considering requirements of the CVT variational framework, this article proposes a feature point based scheme (including boundary points, feature points, and pre-extracted structural points of interest) as initial Voronoi sites.For optimised spatial distribution of these sample points, we calculate a robust discrete mean curvature as density function, which is based on the recognition of curvature's flexibility in capturing shape characteristics and capability in conducting shape evolution (Banchoff, 1967;Kennelly, 2008;Pan et al., 2012).Curvature's ability to flexibly describe terrain morphology has been appreciated by many researchers.For example, P. J. Kennelly noted that, compared to the results of the flow accumulation model, curvature based delineation of drainage networks is not limited to one pixel thickness and requires no depression filling (Kennelly, 2008).The robust discrete curvature calculation is referred to Meyer et al. (2003).

Improved CVT implementation from approximation
The Lloyd Relaxation demonstrates an effective way for heuristic approximation.To follow this elegant approximation, an edge bi-sectioning based dual operation (Du et al., 2010) approach is utilised.That is: from the sample points, an initial TIN surface is constructed.We compute its dual mesh and take this space tessellation as approximated Voronoi cells.The approximated Voronoi tessellation is then optimised within the cCVT iteration.But different to clustering approach, we use each approximated Voronoi cell to (vertically) clip the original dense DEM surface, called referring patch.From this exactly clipped referring patch, we compute accurate energy estimations for the new sites.The global clipping computation is localised using a kd-tree structure.The localisation and accurate referring energy computation makes the cCVT iteration converge fast.The efficiency of the cCVT approximation as a whole is comparable to that of the elegant clustering approach (also has kd-tree fast location embedded).We go no further with the complexity analysis, but provide an implementation of the classical clustering algorithm with the same settings as the cCVT method in the attachment.The pseudo-code of this improved cCVT iteration is described as follows.Here, we illustrate this algorithm using a numerical mountain model.The analytic equation is (5) It has two peaks, two saddles and a pit.We rasterise it with a 4949=2401 regular grid (Figure 1, left).For effectiveness, we set the transformation scale to 0.1 (Ratio = 0.1), that is, there are approximately 240 points left.The sample set includes 56 boundary points, 5 critical points, and an additional 169 random points for visual saturation purposes (Figure 1, left; the red points are randomly generated points, the blue points are boundary points, and the green points are critical points).Relief feature points are always abundant in a real terrain dataset, so additional random points are rarely needed.A robust mean curvature estimation is computed on the original high-resolution surface oriPd (Figure 1, right), from which we can clearly distinguish critical points as peaks, saddles, and pits.The initially approximated TIN surface from the sample set is shown in Figure 2 (left), and its generated dual mesh is shown in Figure 2 (right), which corresponds to step 3 of Algorithm2. Figure 3 shows the dual cell of the sample points, which is the key idea of the cCVT approximation.Figure 4 and Figure 5 show the algorithm steps 4.3.1 and 4.3.2,respectively, where the exact clipping is completed on the original DEM surface.Figure 6 and Figure 7 show the final computation on the reference patch of the first sample point, which corresponds to the algorithm steps 4.4.1 and 4.4.2. Figure 8 exhibits the result of the first iteration, compared to that of the final iteration, with the initial sample points included (top).A comparison of the constructed approximate TIN grids of the initial state and final state is illustrated in the middle, while the curvature distribution that represents the terrain feature comparison is illustrated at the bottom.
The results show how the embedded stationary points (control points and boundary points), feature points, and random points are spatially equalised (Figure 8).Additionally, the cCVT generated a variable-resolution terrain grid (middle right); the convergent TIN grid exhibited nearly uniform high quality, and the convergence process generally resembled Lloyd's Relaxation (Figure 9).Notably, the direct reference on the original DEM surface is realised by the exact geometry clipping, which linearly interpolated the high-resolution surface.This clipping technique has several important benefits: it guarantees accurate energy estimation, it avoids the generation of invalid clustering cells or zigzagging cells, and it promises exact site position calculation, which will result in improved grid quality.

Experimental datasets
Two sites with significant geomorphological characteristics were selected.Experimental site 1 is Mount St. Helens, located in Skamania County, W.A., USA.This mountain is an active volcano whose last eruption occurred in May 1980, and deep magma chambers have been observed recently (Hand, 2015).This site was selected for its typical mountain morphology along cone ridges and evident fluvial features downhill, where heavy pyroclastic materials and deposits are present.These two distinctively different terrain structures mingle together, posing challenges for DEM generalisation.The St. Helens dataset was selected from the Puget Sound LiDAR dataset (http://wagda.lib.washington.edu/data/type/elevation/lidar/st_helens/).This LiDAR dataset was collected in late 2002.The selected dataset is a 29243894 regular grid with a 3 m cell size and covers an area of approximately 102 km 2 .The elevation ranges from 855.32 to 2539.34 m.The image and hillshade views of these data are illustrated in Figure 10.Experimental site 2 is the Columbia Plateau, USA.This area has been labelled UTM zone 11, so we hereafter call it UTM11 (http://gis.ess.washington.edu/data/raster/tenmeter/).This LiDAR dataset was collected in 2009.The selected site is located on the border between Columbia County and Walla Walla County, WA.The southeastern corner is located in the Wenaha-Tucannon Wilderness, Umatilla National Forest.This area contains rugged basaltic ridges with steep canyon slopes at high elevations (average of 1700 m).The northwestern area is located near Dayton City, which is a vast agricultural and ranching area with relatively smoother morphology at low elevations (average of 500 m).This site is selected due to the coexistence of these two different prominent surface morphologies.If the generalisation scheme emphasises the high-elevation areas with sharp variations, the surface interpolation as a whole might be unbalanced, which may result in smoothing of the low-elevation areas.The selected UTM11 dataset is a 3875  3758 grid with a 10 m cell size and covers an area of over 1456 km 2 .The elevation ranges from 3533 to 19340 cm.The image and hillshade views of these data are shown in Figure 11.

Comparison method
As previously mentioned, DEM generalisation has long been studied in geoscience, and numerous methods have been proposed over time.One of the most classical approaches is the hierarchical insertion (or decimation) of feature points to construct a TIN grid under a destination scale.This type of heuristic feature point refinement (HFPR) performs very well in terms of surface approximation and terrain structure retention.For this reason, although HFPR methods generally cannot guarantee high-quality grids, these methods are suitable for comparison purposes.
A typical HFPR starts with four corner points from a dense DEM image and constructs a Delaunay triangular grid that contains two triangles.The rest of the points are weighted according to their distances from the triangular surface or other error criteria and then queued.The point with the highest priority in the queue is selected, and the grid is modified using incremental Delaunay triangulation.This process repeats until some error threshold is satisfied (Heckbert and Garland, 1997).Michael Garland provided a classical HFPR implementation (http://mgarland.org/software.html), and many other variants are available in GIS, meshing, and visualisation tool suites.

Quantitative comparisons
We performed the processes from Algorithm2 for the two experimental datasets, including triangulation and curvature estimation, boundary point extraction and marking, feature point extraction based on curvature significance and marking, and optimisation loop through cCVT.For effectiveness, the transformation Ratio was set to range from 5% to 0.1% point density (comparable to the 3.1% to 0.6% setting in (Zhou and Chen, 2011)).The accuracy of the surface approximation determines the final surface interpolation precision and is thus a basic quality comparison index.Here, we applied a statistical interpolation method to measure the surface approximation precision.From each triangle on the cCVT-generated quasi-uniform TIN grid, a random point is selected and a vertical line is introduced to intersect the original dense DEM surface and the HFPR generalised TIN at the same time.Error estimates for the surface approximation could be obtained from these intersection points.We computed the mean error, maximum error, and root mean squared error (RMSE) for this elevation interpolation (TIN interpolation); the results are listed in Table 1.Furthermore, we computed the aspect ratios of the triangles for both generalised TIN surfaces to measure the grid quality, which are also listed in Table 1.RMSEs with various transformation Ratios are listed in independent rows and columns in Table 1.From the results in Table 1, we can conclude that under the same resolution (point density), the transformed DEM surface obtained using the cCVT method is generally more precise than that obtained using the HFPR method.In all cases, surface approximation precision (compared to the original) decreased as the resolution coarsened.

Qualitative comparison
A qualitative index is usually measured from the aspect of terrain structure retention.According to the resulting TINs from the two experimental datasets, both the cCVT and HFPR methods performed well based on visual examination.However, upon closer inspection, the surface generated by cCVT has a smoother transition effect than that generated by HFPR (Figure 12).HFPR accumulated more samples around sharp features (c.f. Figure 13), and its surface exhibited clearer impressions because flat details were smoothed out.From visual examination, it may be concluded that under the same transformation conditions, HFPRs may exert a stronger generalisation effect than cCVTs.However, a stronger generalisation effect actually decreases the precision of the general approximation, which may result in structural distortion or misconfiguration.Figure 14 illustrates a closer examination of St. Helens.Some structural details on the original surface were recovered by the cCVTs but not by the HFPRs.This terrain structure loss occurred on both smooth areas and steep areas, as illustrated in Figure 14. Figure 15 illustrates similar structural detail losses by the HFPRs in the UTM11 dataset.Terrain structural features could also be measured from DEM derivatives such as the slope, aspect, and hydrological structural lines.Here, we used contours to compare the generalisation accuracy using experimental site 2. For the same configurations (80 m elevation increments), we generated contours for the original dense TIN (rendered in red), the cCVT-generated TIN (rendered in blue), and the HFPR-generated TIN (rendered in black), and overlaid the three sets of contours for comparison (Figure 16).The illustrations demonstrate that in most cases (Figure 16 b, c), the contours from the cCVT-generalised surface accurately conform with those from the original dense surface, while the contours from the HFPR-generalised surface generally did not, except for some cases in steeper areas with sharp curvature variations (Figure 16 d).This result can be explained by the HFPR's stronger accumulation of sample points near sharp features, which guaranteed an edging out, if we noticed that the inspection area d is much smaller than b or c.

Discussion
Topography transformation of DEM surfaces has been a deeply studied topic in geoscience, simplification techniques and generalisation principles are widely realised and adopted.Extracting terrain feature points and using these points to construct a generalised surface has proven to be one valuable approach; its success may be due to the feature points' strong capability to capture terrain structures.However, discrete surface as TIN grids that are constructed purely from feature points may not be best approximated to the original high-resolution surfaces.Take the mountain equation in Figure 1 as an example.It has at least two peaks, two saddles, and one pit close to the zero level.Assume that the scale transformation requires that only two critical points are left; selecting both peak points is more reasonable than selecting the pit point, even if the pit point has a stronger quantitative index (curvature in this case) than the peak points.This observation implies that if global surface interpolation precision is of more importantly demanded, a robust approach that has overall considerations on surface approximation and terrain feature retention should be adopted.Among those classical DEM generalisation approaches, heuristic feature point refinement is an outstanding example.As illustrated by Table 1, Figure 12, and Figure 16, HFPR methods perform excellently in terms of surface approximation and surface morphology retention.For the treatment of feature points, these methods use a heuristic strategy by incremental Delaunay triangulation, which considers the point with the largest error with respect to the constructed TIN.However, the impact of the insertion of a new feature point on the inserted feature points is not considered due to computational burden.That is, modifications are only taken place on the triangle where the point with largest vertical error locates on.As a result, feature points may cluster around reliefs with sharp variations, as shown in Figure 12.Too many feature points accumulating near sharp features means that relatively few feature points are present in flat areas, which will eventually lead to terrain structure misconfiguration, as shown in Figure 14, Figure 15, and Figure 16.Sometimes, this type of structural loss is unacceptable.For example, the terrain relief at high elevation under the studied scale (10 m cell spacing) in experimental site 2 exhibited a fiercer landform than at lower elevations.The accumulation of too many sample points in high-elevation areas may result in the distortion of the smooth anthropogenic terrain morphology in low-elevation areas.cCVT starts by constructing a terrain-adaptive variable-resolution grid.The cCVT iteration uses a robust mean curvature as density function which is based on the curvature's capability to characterise shapes and conduct shape evolution.Under the curvature guidance, the two-step optimisation (c.f.Algorithm1 in Section 2.2) loops both to spatial equalisation and frequency equalisation.The process of spatial equalising of feature points has been seldom considered by classical approaches, which may explain why cCVT generally prevails over HFPRs (c.f.Table 1).Notably, the triangles from the spatial equalisation exhibited a maximum aspect ratio that was less than 5.0, which implies that the constructed terrain grid satisfied the numerical stability requirement from classical finite element or finite volume computations.On the other hand, CVT is an approach within variational framework.The result of the iteration depends on the boundary conditions and initial conditions.Hence, this article employed a feature point scheme (with additional input points considered) as a relatively stationary initial condition to maintain algorithm stability.The requirement of embedding feature points of interest, along with consideration of avoiding the problematic k-means like clustering, prompted us to develop a non-clustering approach with an exact energy referring method.Experiments on ten million DEM points demonstrated that the exact clipping approach performed comparably to the elegant clustering approach.

Conclusions
In this article, a high-fidelity multiresolution DEM model was proposed.The variable-resolution with high-fidelity was 5 achieved by the developed curvature-based CVT.cCVT optimisation increases the precision of surface approximation compared to existing heuristic DEM generalisations, while the equalisation of feature points from the spatial domain guarantees a high-quality grid.
Multiresolution models are essential tools to incorporate more scales, while a high-fidelity generalised DEM model can be used to construct a concrete topographic layer from which fine endogenous or exogenous processes can be assessed 10 under proper scale conditions.Evaluation of the cCVT multiresolution DEM model on Earth and environmental systems over wide-ranged domains and scales is a topic for future studies.Furthermore, considering the topography over a wider range may require re-implementing cCVT on the geographical coordinate base instead of on the currently used projection coordinate base..
Use kt-tree for fast intersection computation of point ( ̅,  ̅,  ̅) and oriPd, with the result used as the projected nearest point; push it into the new candidate point set  ' ; } 4.5) Using {B, C} as constraints, Delaunay triangulate point set {B, C, F ' } and obtain reconstructed TIN'; 4.6) Compute E on TIN ' ; }

Figure 1
Figure 1 High-resolution grid of a mountain equation (5).Left: original grid oriPd in 49×49 resolution, rendered in mean curvature.The sample points were also rendered on oirPd; the blue points are boundary points, the green points are critical points, and the red points are random points.Right: robust mean curvature estimation.The saddle terrain features, peaks, and pits can be distinguished, compared to the depiction on the left.

Figure 2
Figure 2 Initial TIN surface (left) and its dual grid TD (right).The initial sample points on the dual grid are also rendered.

Figure A triangle
Figure A triangle Tj (rendered in semi-transparent blue) with its minimal bounding box's intersection part narrPd with the original grid on the approximate grid (left), the localised intersection part narrPd (middle), and the intersection part on the original grid oriPd (right).

Figure
Figure Exact clipping steps of narrPd with Ti.The sequence from left to right illustrates the edge clipped results.

Figure 6
Figure 6 Reference patch refj on the original DEM surface of an initial Voronoi cell, with the centre at ri (left).Right: refj on the original DEM surface oriPd.

Figure 7
Figure 7 Barycentre computation based on the reference patch refj; the original site is the white block in the circle, and the projection on oriPd of the newly computed site is depicted as the blue block in the circle.

Figure 8
Figure 8 Comparison of converged results.Top left: reconstructed TIN surface from one iteration, with the initial points presented.Top right: the converged TIN surface with the initial sites, after approximately 140 loops.Middle: the initial approximated TIN surface (left) and the final TIN surface (right).Bottom: curvature distribution on the original surface (left) and the generalised grid (right).

Figure 9
Figure 9 Trajectories of point convergences.The red points indicate the initial sample set, and the trajectories show the convergence trends, with closer gaps between candidate points.The right side shows a close view of the convergence of two points.These trends imply that the cCVT's convergence complies with Lloyd Relaxation linear convergence.

Figure 10
Figure 10 Experimental Site 1: Mount St. Helens.Left: image view.Middle: hillshade view.It is a 2924x3894 grid with a 3 m cell size.

Figure 11
Figure 11 Experimental Site 2: UTM11 Zone.Left: image view.Middle: hillshade view.It is a 3875 x 3758 grid with a 10 m cell size.

Figure 12
Figure 12 Visual examination of St. Helens.Top: cCVT grid.Bottom: HFPR grid.The latter grid appears more rigid than the former, which implies a stronger generalisation effect.

Figure 13
Figure 13 Grid quality from an intuitive comparison.Top: cCVT-generalised grid with nearly uniform triangles.Bottom: HFPR generalised grid with irregular triangles.

Figure 14
Figure 14 Detail loss of the HFPR generalisation grid.The inspected area in the experimental site is bounded by the white rectangle (a); the magnified inspection area on the original dense TIN (b); the area on the cCVT-generalised TIN (c); and the area on the HFPR-generalised TIN (d).From the close-up view, it can be seen that the HFPR method generates a rougher grid than the CVTs.Thus, structural distortion or misconfiguration might be introduced.

Figure 15
Figure 15 Detail loss from the HFPR method in the UTM11 dataset.(a) The inspection area in the entire experimental site.The magnified view of the inspection area is shown on the original dense DEM (b), on the cCVTgeneralised TIN (c), and on the HFPR-generated TIN surface (d).The fold morphology in the white ellipse was recovered by the cCVT method but not by the HFPR method.

Figure 16
Figure 16 Comparison of the derived contour lines.The contours from the dense UTM11 dataset are shown in (a) and rendered in red.Three areas in the boxes were magnified to show the differences in the contour configurations.The blue contour lines are from the cCVT-generalised TIN surface, and the black lines are from the HFPR-generated TIN surface.According to the illustrations in (a) and (b), the contours that were derived from the cCVT grid (rendered in blue) are more in line with those from the original dense surface (rendered in red).The contours from the HFPR grid (rendered in black) may sometimes edge out on areas with steep slopes, as shown in (d), because the HFPR method accumulated a relatively abundant of sample points around these areas.
Compute its minimal bounding box BBoxj, quickly compute its intersection with oriPd using a kd-tree, obtain a narrowed reference geometry narrPd; 4.3.2) Compute exact intersection of Tj and narrPd, push the result into reference sets REF={refj}; } 4.4) For each refj in REF