AGLOBAL “ NATURAL ” GRID MODELING METHOD BASED ON TERRAIN TOPOLOGICAL STRUCTURE

In this paper, a Global “Natural” Grid (GNG) modeling method based on the terrain topological structure is proposed. GNG partitions the earth surface into “natural” cells by using the natural terrain feature elements, which may be more suitable for exploring and understanding some natural phenomena in the environment monitoring, hydrological analysis, species distribution, etc. This modeling method mainly consists of three steps: firstly, a variable-resolution geometric representation of global DEMs based on the restricted quadtree is constructed; secondly, the feature elements are extracted and an original GNG model is built; finally, a multi-resolution representation of GNG is established by recursively applying topological simplification. In the end, an experiment is done to validate the correctness and feasibility of this modeling method.


INTRODUCTION
Global "Natural" Grid (GNG) partitions the earth surface by using the natural terrain feature elements, such as feature points (pits, peaks and saddles) and feature lines (such as ridges and ravines).It decomposes the earth surface into "natural" cells rather than "regular" cells by the traditional Discrete Global Grid (DGG), such as Quadrangle Grid [Seong 2005], Quaternary Triangular Mesh [Chen et al. 2003], or Hexagon Grid [Sahr 2008] and Voronoi/TIN grid [Mostafavi et al. 2003].In fact, GNG can directly reveal the topological structure information related to the level sets of the height function defined on the earth surface and supports a knowledge-based approach to analyze, visualize and understand the earth surface behavior (in space and time) through a global perspective.So, as a topological representation, GNG model may be more suitable for exploring and understanding those natural phenomena in the environment monitoring, hydrological analysis, species distribution, meteorological service, climate change, and epidemic promulgation, etc.   The basis of terrain topological representation is Morse theory [Milnor 1963, Smale 1960].Morse theory can describe the shape of a surface by those fundamental elements (critical points, critical lines, and their relationships) that effectively act as landmarks on the surface and achieve a representative coverage of the entire surface, called Morse (or Morse-Smale) complex.In recent years, lots of algorithms for terrain morphology based on Morse Complex have been presented (see [Vatali et al. 2012;Biasotti et al. 2008] for more details).However, most of them and the related applications are limited to local terrain datasets and based on 2D planar mode.According to our knowledge, there are few reports on how to construct a multi-resolution topological representation for the * Corresponding author at: Geoscience and Surveying Engineering College, China University of Mining and Technology (Beijing), Beijing 100086, China.Tel: +86 13810439630.
whole earth surface on spherical mode.In this paper, a GNG modeling method is proposed based on the terrain topological structure.
The remainder of this paper is organized as follows.In Section 2, some basic notions on Morse Theory are introduced in briefly.In Section 3, a modeling method of the global "nature" grid is presented in details.Next, an experiment is designed and the experimental results are analyzed.In the end, some concluding remarks and future works are given.

MORSE THEORY
Morse theory is a powerful tool to capture the morphology of manifold on which a scalar function f is defined.It decomposes the shape into cells in which the gradient flow is uniform, and encodes the adjacencies of these cells in a complex which represents the topology as well as the geometric properties of the gradient of f.
Let f be a C 2 -differentiable real-valued function defined over a domain D in the 2D space.A point p of D is a critical point of f if and only if the gradient of f vanishes on p. Function f is said to be a Morse function when all its critical points are non-degenerate (i.e., if and only if the Hessian matrix Hesspf of the second derivatives of f at p is non-singular: the determinant of Hesspf is not zero ).This implies that the critical points of f are isolated.The number of negative eigenvalues of Hesspf is called the index of a critical point p.There are three types of non-degenerate critical point p: a minimum (pit), a saddle (pass), or a maximum (peak) if and only if p has index 0, 1 or 2, respectively.Points which are not critical are said to be regular.
An integral line of f, going through a point p, is a maximal path which is everywhere tangent to the gradient of f.Integral lines that converge to (originate from) a critical point p of index i form an i-cell ((2-i)-cell) called the descending (ascending) cell of p.The collection of all descending cells form a Euclidean cell complex, called the descending Morse complex (Figure 1a), and the collection of all ascending cells form also a Euclidean cell complex, called the ascending Morse complex International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-4/W2, 2013ISPRS WebMGS 2013& DMGIS 2013, 11 -12 November 2013, Xuzhou, Jiangsu, China Topics: Global Spatial Grid & Cloud-based Services (Figure 1b), which is dual with respect to the descending Morse one.A Morse function f satisfying the Morse-Smale condition, that is, descending and ascending Morse complexes intersect only transversally is called a Morse-Smale function.This means that the intersection of the edges of the descending and ascending cells is a saddle point.In this case, we can define a complex, called the Morse-Smale complex, whose cells are the connected components of the intersections of the descending and ascending cells (Figure 1c).In terrain analysis, the separatrix starting from saddle to peak is called the radge line (or radge), while the one starting from saddle to pit is called the ravine line (or ravine).Essentially, the surface network [Pfaltz 1976, Schneider andWood 2004] consisted of ridges and ravines is the critical net.According to Morse theory, it is impossible that valley and ridge line intersect at regular points.In this paper, both Morse (or Morse-Smale) complex and critical net (or surface network) for representing the terrain morphology are called the terrain morphological model.

A GNG MODELLING METHOD
Above all, Morse theory investigates the relation between the properties of a function defined on a manifold and the properties of the manifold itself.Fortunately, the earth surface is manifold, where Morse theory can be used directly.In general, a GNG modeling method mainly consists of three steps: firstly, a variable-resolution geometric representation is constructed based on the restricted quadtree; after that, the feature elements are extracted and an original GNG model is constructed; finally, the topological simplification and refine operations are applied to the original GNG model recursively and a multi-resolution representation of GNG is established.

A Variable-resolution Geometric Representation
Since the data volume of the global DEMs is huge, a variable-resolution representation of the earth surface is more appropriate to achieve fine expression in land area while rough expression in marine region (about 71% of the whole earth surface).However, some cracks or T-junctions may occur, which is adverse to the extraction of the feature elements.In this section, based on the restricted quadtree [David 2002], a triangulated global diamond grid is adopted to construct a seamless expression of the global variable-resolution DEMs without long and narrow triangles.

Global Diamond Subdivision:
The main principle of the global diamond subdivision is as follows: an octahedron embedded in a spheroid is selected as the framework of the spherical grid partitioning, of which the projection edges coincide with the 0, ± 90, ( ± )180 meridian respectively.In this way, four original facets, called basic diamonds, can be constructed first of all (as shown in Figure 2).Similar to the planetary quadtree tiling, each facet is redivided into four child-diamonds recursively until some desired resolution is satisfied.As a result, the whole earth surface is simulated and the diamond subdivision is encoded by a quadtree.
Figure 2. The global diamond grid [Zhao et al. 2007] As this diamond grid is very similar to the regular square grid, it gets easy in the spatial operations especially in the neighborhood search and conductive to simulate a variable-resolution earth model encoded by the restricted quadtree.
To construct the seamless expression of the global variable-resolution DEMs, the earth surface is uniformly partitioned to a certain level at first.After that, the maximum elevation difference of each diamond vertexes is calculated and compared with the given threshold, and then according to the algorithm of restricted quadtree tiling [David 2002], the diamond node is determined to be recursively refined, until the maximum elevation difference does not exceed the threshold value any more.As a result, the representation level of the flat ocean region is progressively shallower than that of the fluctuant land area, where the level between two adjacent diamonds exceeds one at most.Finally, leaf node is triangulated according to the adjacencies of the quadtree nodes.

Vertex Calculation:
In this section, we mainly address the problem of the vertex coordinates of the child diamond calculated directly from the father node.In global diamond gird, the latitude of any child node vertex is easy to be obtained by averaging the latitude values of the father node vertexes.However, the longitude calculation of the child node vertexes is complicated.For brevity, as shown in Figure 3, the calculation with the central latitude M0>0 is only introduced in details to get the longitude of the child nodes from the father diamond, while other conditions are similar and skipped over.As shown in Figure 3b, let L1, L2 is the west and east longitude value of a father diamond (remarked in red in Figure 3a), and L3, L4, L5, L6 is the longitude of child nodes to be calculated.On the M0 latitude line, the child diamond number N on its central latitude line is 90/(L2-L1) at current level while twice at the next level at the same latitude line.Therefore, the child diamond number is 2N-1, 2N+1 on the latitude line M1 and M2 respectively since M0>0.So, L3, L4 can be formulated as follows: Here, L0 is the west longitude of the basic diamond and k is an integer of which the value is between 0 and 2N-1 and can be determined according to the following inequality: Analogously, L5 and L6 can be obtained.Finally, all the vertexes of each child diamond can be located in this top-down calculation, while the ones of the root node, that is, the basic diamond are known after all.According to the geographic coordinate, the evaluation of any vertex can be obtained from the database of the global DEMs.

A Modified Algorithm of Constructing GNG
The STD algorithm [Magillo et al. 2009] has the advantage of being simple, using only comparisons (and no floating-point computations).However, it may induce intersections between ridges and ravines at the regular points, which is not consistent with the Morse theory.Moreover, this algorithm is only suitable for local terrain datasets.Some improvements need to be done to prevent intersections and to make it suitable for the topological representation of the Global DEMs.In this section, a modified algorithm of constructing GNG is presented by taking advantage of the duality of Morse Complexes.The main principle is as followings: From the dual structure of the ascending and descending Morse complex (in Figure 1), it can be observed that the boundary of an ascending cell is consisted of ravines and saddles, and that the ridges starting from the boundary saddles to the central peak route within the cell.Dually, the boundary of a descending cell is consisted of ridges and saddles, the ravines starting from the boundary saddles to the central pit route within the cell.
The modified algorithm mainly consists of three steps.Firstly, one of Morse complex, descending or ascending Morse complex, is constructed by the STD algorithm [Magillo et al. 2009] after the topology relationship of each two adjacent vertexes and two adjacent triangles as well as vertex and triangle is established in the triangulated diamond grid.Here, all the marine areas are identified first of all and each is regarded as a pit region according to the common sense.Secondly, ridges (ravines), that is, the separatrix of the dual Morse complex is extracted in every constructed cell to prevent those intersections.The saddles belonging to a cell are classified into two groups: the boundary ones or the internal ones.So, a cell saddle whether located on the boundary or not need be determined.What's more, there may be some macro-saddles, that is, one feature line passes two saddles and the path between the two saddles (called macro-saddle line) plays the dual role of both ravine line and ridge line.Therefore, all the macro-saddles are identified and any of the macro-saddle lines is extracted according to the corresponding path established in the first stage.Furthermore, it should be noticed that the ravines starting from the saddles located on the boundary of the marine regions are not extracted and the end point of each ravine is a pit or a macro-saddle or a point of some marine area.Finally, regarding the feature lines extracted in the second stage as the boundary, the dual complex of the former is constructed and the topological representation of the earth surface, i.e. the original GNG is constructed.
Taking the construction of descending Morse complex in the first stage for example, the flow chart of constructing the original GNG is shown in Figure 4.

A Multi-resolution Representation of GNG
It is well known that the result and conclusion may vary as the scale changes when we observe and analyze a geography phenomenon.Therefore, a multi-resolution representation of the terrain morphology is critical for interactive analysis and exploration of the earth surface in order to maintain and analyze characteristic terrain features at different levels of abstraction.To this aim, topological simplification is indispensable.
The basic topological simplification operation is collapsing a maximum-saddle pair into a maximum, or a minimum-saddle pair into a minimum, so as to maintain the consistency of the underling complex.Taking the former (in Figure 5) for instance, the peak u and the saddle v are adjacent, that is, there is a ridge connecting u and v in the critical net.Therefore, u and v can be paired.Moreover, there is also a peak point w connecting with v. Since the persistence [Elebrunner 2001] of w is larger than that of u, the pair u-v is deleted and the ridges previously terminating to u is assigned to terminate to w after one topological simplification operation is proceed.In addition, as v is not a saddle any more, the two ravines connecting v are eliminated from the critical net.Symmetrically, the minimum-saddle pair also can be cancelled.However, if taking care of the presence of macro-saddles and also of multiple saddles in the real terrain surface, the simplification becomes more complicated (see [Danavora et al. 2007, Floriani et al. 2010] for more details).In summary, the co-separatrix of any two adjacent ascending or descending Morse cells will be cancelled, as well as the critical point pair with smaller persistence in every topological simplification operator.Certainly, the reverse operation, that is, refinement is also can be implemented according to the requirement of multi-resolution representation.The construction of a multi-resolution GNG mainly consists of three steps.Firstly, adjacent cells are paired.Secondly, and the persistence of each pair is calculated.Thirdly, simplification or refinement operations are applied and a multi-resolution GNG is achieved.

EXPERIMENTS AND DISCUSSIONS
An experiment system with C++ and OSG (Open Source developed by using GTOPO30 data is designed to validate the feasibility and correctness of our modeling method.
Based on the restricted quadtree, a triangulated variable-resolution global diamond grid is constructed without cracks and T-junctions as well as thin and long triangles, as shown in Figure 6a.Here, we uniformly partition the earth surface into the diamond grid at level 3 first of all, and then each diamond node is refined recursively to the max level 7, if the maximum elevation difference of the diamond vertexes exceeds the given threshold 10.0 m.As a result, as shown in Figure 6a, the resolution of the land areas is higher than that of the marine regions and the level between two adjacent diamonds does not exceed one any more.
The original GNG (in Figure 6c) of the global DEMs is constructed by the modified STD algorithm and it can be regarded as the overlap of the global descending Morse complex (in Figure 6b) and the global ascending Morse complex (in Figure 6c).In hydrology, each descending cell is a catchment area and all the rivers converge into the sea.What's more, attention should be paid that some islands occurs in the marine areas and represent as isolated peaks (in Figure 6b), that is, no saddles and ravines connect them.Here, we regard these isolated peaks as the degeneration of the boundary of the marine descending cell, that is to say, all the ridges on the boundary of the marine descending cell degenerate into a peak point.Similarly, a ravine may degenerate into a saddle if the saddle is a boundary point of some marine region.As a result, these situations indicate that the Morse-Smale condition is not satisfied any more in the GNG.
By applying the topological simplification operators (or refine operators), a multi-resolution GNG is achieved, as shown in Figure 7. Gradually, 2000, 4000 and 6000 steps of simplification operators are averagely assigned to the global descending Morse complex and the global ascending Morse complex and the responding terrain morphology is shown in Figure 7a, 7b and 7c respectively.As a result, the main terrain features, that is, higher mountains and deeper canyons in the earth surface is maintained.Meanwhile, the primary global terrain structure is obtained through the multi-resolution representation of GNG.

CONCLUDING REMARKS
In this paper, a multi-resolution GNG model of the global DEMs is proposed.Decomposing the earth surface into "natural" cells by using the natural terrain feature elements, GNG directly reveals the topological structure information of the earth surface and supports a knowledge-based approach to analyze, visualize and understand the earth surface behavior (in space and time) through a global perspective.Some contributions of this paper are as followings: 1) To extract the natural feature elements, a variable-resolution geometric representation of the global DEMs without cracks, T-junctions and thin and long triangles is presented based on the restricted quadtree.
2) Taking advantage of the duality of the two Morse Complexes, a modified algorithm of constructing the original GNG is developed for the morphology of the earth surface.
3) By applying the simplification and refinement operators, How to improve the efficiency of modeling GNG is a issue that we consider in the future works.
Integral lines connecting saddle points to other critical points are called separatrices.All separatrices forms a network, called the critical net.It is the 1-skeleton of the Morse-Smale complex.Morse and Morse-Smale complex: (a) Descending Morse complex, (b) Ascending Morse complex, (c) Morse-Smale complex[Danavora et al. 2007]

International
Figure 3.The structure of the diamond grid

Figure 4 .
Figure 4.The flow chart of constructing the original GNG The basic operation of topological simplification: (a) before simplification, (b) after simplification Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-4/W2, 2013 ISPRS WebMGS 2013 & DMGIS 2013, 11 -12 November 2013, Xuzhou, Jiangsu, China Topics: Global Spatial Grid & Cloud-based Services a multi-resolution GNG model of the global DEMs is finally given, which is suitable and critical for exploring and understanding the natural phenomena.

Figure 6 .
Figure 6.The original geometric and topological representation of global DEMs: (a) the variable-resolution representation, (b) global descending Morse complex, (c) global ascending Morse complex, (d) the original GNG