Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2013 / Article
Special Issue

Intelligent Techniques for Simulation and Modelling

View this Special Issue

Research Article | Open Access

Volume 2013 |Article ID 523729 |

Huijie Zhang, Zhiqiang Ma, Yaxin Liu, Xinting He, Yun Ma, "A New Skeleton Feature Extraction Method for Terrain Model Using Profile Recognition and Morphological Simplification", Mathematical Problems in Engineering, vol. 2013, Article ID 523729, 16 pages, 2013.

A New Skeleton Feature Extraction Method for Terrain Model Using Profile Recognition and Morphological Simplification

Academic Editor: William Guo
Received10 Jul 2013
Accepted12 Aug 2013
Published24 Sep 2013


It is always difficul to reserve rings and main truck lines in the real engineering of feature extraction for terrain model. In this paper, a new skeleton feature extraction method is proposed to solve these problems, which put forward a simplification algorithm based on morphological theory to eliminate the noise points of the target points produced by classical profile recognition. As well all know, noise point is the key factor to influence the accuracy and efficiency of feature extraction. Our method connected the optimized feature points subset after morphological simplification; therefore, the efficiency of ring process and pruning has been improved markedly, and the accuracy has been enhanced without the negative effect of noisy points. An outbranching concept is defined, and the related algorithms are proposed to extract sufficient long trucks, which is capable of being consistent with real terrain skeleton. All of algorithms are conducted on many real experimental data, including GTOPO30 and benchmark data provided by PPA to verify the performance and accuracy of our method. The results showed that our method precedes PPA as a whole.

1. Introduction

Terrain model is an important component in the actual engineering about 3D geographic information system and is applied widely in the field of spatial information visualization, including battleground simulation, flight simulation training, and terrain deformation simulation. However, it is a big challenge to construct a terrain model, because terrain model usually contains massive data, and its feature is rather complex. A valid approach of building terrain model is to extract the key feature of terrain instead of using all of the terrain data. The aim is not only to guarantee the accuracy of terrain model but also to decrease the data scale.

Over the years, many skeleton feature extraction methods have been proposed, which are mainly based on terrain flow simulation [13] or geometrical gradient analysis [46]. The moving-window algorithm is an early classical algorithm, which takes elevation map as 2D image [4]. Although this algorithm can extract the basic outline feature of terrain, it does not consider the morphological information but just connects simply the extreme value of elevation in window into feature lines. Therefore, the accuracy is very low. Another classical algorithm is deterministic eight nodes (D8) proposed by O’Callaghan and Mark, which obtains the ridge lines by simulating the procedure of the surface water flow from 8 directions [7]. A series of improved algorithms are proposed based on D8 [811]. However, since drainage accumulation matrixes and flow direction matrixes are complex so as to take a lot of time overhead to calculate, the efficiency is not acceptable in many application fields.

Furthermore, for some special terrain, such as depression and flat, flow simulation algorithms are not capable of extracting the feature lines correctly [12]. The main reason is that flow simulation algorithm depends on the outlet to establish the flow direction matrixes. Nevertheless, it is possible that there is no outlet for depression and flat. Hence, depression filling algorithm is a relative valid approach to solve this problem [1315]. However, it may increase the area of flat and thus suffer from the new problems [16]. In addition, depression filling also leads to extra time consumption. Barnes proposed the simplified algorithm to improve the efficiency [15]. For terrain feature of flat, researchers also proposed the related strategies, including graph structure [17], elevation increments method [18], and radial basis function [19].

Because of the above limitations, most flow simulation algorithms cannot guarantee both the precision of the feature extraction and the global morphology of the feature lines. Chang et al. raise the Profile recognition and the Polygon-breaking Algorithm (PPA) [20]. Its advantage is not to be influenced by depression and flat and to be capable of extracting the global skeleton feature in most instances.

In recent years, PPA is applied in various terrain modeling, including terrain synthesis [21] and terrain deformation [2228], because of its fine characteristics mentioned. Zhou et al. utilized PPA to extract the skeleton feature of terrain and realized the terrain synthesis [21]. Gain improved the terrain synthesis and proposed the method based on user’s sketch map [29]. In addition, there are other applications of PPA in terrain visualization, including multiresolution terrain deformation [22] and enhanced terrain synthesis based on GPU [30]. In order to deal with global large scale data, the related work, such as artificial intelligence planning [31] and the differential evolution algorithm [32, 33], has also referential significance.

Nevertheless, PPA cannot reserve the ring feature, and the accuracy and stability of algorithm need to be improved further. In this paper, we proposed a new method, Morphological Profile recognition and Polygon breaking (MPPA). Its primary advantage is that MPPA is capable of recognizing the stable rings and reserving them. Moreover, MPPA defined the outbranching in order to remain the sufficient long truck lines, which is consistent with real terrain feature. Furthermore, the noise points have been eliminated by our simplification algorithm based on morphology; therefore, the whole performance of MPPA has been improved greatly.

2. Target Recognition and Feature Simplification

The original data of terrain model is usually stored as the form of Digital Elevation Model (DEM). It is a large scale data set. It is necessary to extract the significant data which contains the skeleton structure information of the terrain from the massive data set. In our work, profile recognition strategy is adopted to recognize the initial feature points as candidate target points, which are similar to the target in PPA. Its advantage is that it can almost cover most of the feature points.

PPA connected all of the candidate target points into polygon areas and then disassembled those closed polygons in order to extract the ridge axes. Owing to there are a lot of noisy points, the efficiency of algorithm is really low. Moreover, its principle of removal of noisy points is the inexistence of the closed polygon, which did not accord with the original ridge morphology. Different from PPA, our method proposed a new strategy to refine the initial feature points by morphological simplification algorithm, which can eliminate a lot of noisy points in large quantity. And then we utilized the simplified feature points, that is, subset of initial target points, to exact the skeleton feature lines.

2.1. Profile Recognition

The profile recognition is a classical method that extracts the candidate target points from the original DEM data set [20]. For the adjacent points on the eight directions of the current central point, they are divided into four couples as (North, South), (East, West), (Northeast, Southwest), (Northwest, Southeast), as shown in Figure 1. Each couple represents a profile of terrain model, which has been identified by mean of different shapes in Figure 1.

Let the number of points on each profile be the profile length as , which is closely related to recognition result of feature points. In this paper, we set as 5, that is to say, there are two adjacent points at each side of the current point in each profile. Afterwards, the current point is determined whether it is a candidate target points according to the relationship between the elevation data of the current point with that of its adjacent points on different profiles. The definition of the candidate feature points set is as follows.

Definition 1. Let be an original set of elevation points, let be current point, let be the number of adjacent points for each profile, and let be the set of adjacent points on the th profile; then is a candidate feature point, if the following constraints are satisfied:(1)on all the profiles, there exist at least two points and ; their elevations are lower than that of ;(2)and should be at the different side on the same profile,where is profile length, are coordinates of current points, and is the sequence number of profile.

By using Definition 1, a set of candidate target points will be extracted, which almost covered the whole ridge feature of the terrain. It is important to note that the constraints are loose enough so as to prevent from missing the real target points. As a result, it is a data set with a lot of redundant noise points, but its whole feature information has been reserved. Figure 2 shows the set of candidate target points which we have obtained by using profile recognition.

2.2. Morphological Simplification

Morphology has presented its superior properties in the field of image processing. In our work, we adopted the morphological idea to solve the simplification problem of candidate target points. In this section, we defined a suit of binary coding sequence and then proposed a new simplification strategy for candidate target points according to the terrain morphology features. It is noteworthy that efficiency of our skeleton extraction method MPPA will be improved greatly, because the feature points are connected into polygon strips after the candidate target points have been simplified by our strategy.

Let be a binary coding sequence, and the mapping relationship between the bit code and the position of elevation point is presented by the following: where is the spatial position of the elevation point .

Obviously, it is easy to implement Formula (1). Nevertheless, we can determine the corresponding spatial position of each bit in according to the mapping relationship. Moreover, there exists a corresponding decimal value for each , and thus each can represent a kind of spatial morphology of terrain, named morphological coding. Since consists of nine bits, the value range of is . Hence, morphological coding can stand for 512 types of morphology of terrain at most. In our work, we defined 48 kinds of morphologies to simplify the candidate target points, shown in Figure 3. Figure 3(a) described the mapping relationship between and . In Figure 3(b), each morphological coding stands for a kind of morphology, which has been denoted as . However there are only 12 types of morphologies illustrating in Figure 3(b). And the other ones can be deduced by rotating these instances.

For each candidate target point, we detect the terrain morphology feature according to the relationship between it with its 8 adjacent points. If the morphology feature matched any morphological coding above-mentioned, then the current point should be removed from the set of candidate target points. The morphology feature of the current point is calculated by where the value range of is 0 or 1 and if is a candidate target point; otherwise .

The detail steps of simplification algorithm about candidate target points are described in Algorithm 1.

Input: The set of candidate target points.
Output: The optimized set of feature points.
1. Assign the 512 bits of contiguous space MTC to
  store the state of morphological coding.
2. Initialize the designed 48 kinds of morphologies in
 MTC into 1, other bits should be set into 0.
3. If there is no change in , then go to step 4;
 otherwise do the following loop:
3.1 For each , calculate its morphology
   feature using formula (2);
3.2 Evaluate whether matches one of 48
     kinds of morphologies;
3.3 Eliminate from set , if MTC
    equals to 1;
4. Algorithm end.

According to Algorithm 1, it is not hard to see that we take the morphology feature of the current point as the index in the contiguous space of MTC. Therefore, searching procedure is not necessary to recognize the morphology of terrain and then determine whether the current point should be eliminated from the set of candidate target points. As a result, the efficiency of algorithm is very high. In our experiment for real data, the outermost iteration is approximate 6 times to complete the simplification procedure.

2.3. Connection of Feature Points

Morphology simplification is significant to improve the efficiency of MPPA. The main reason is that the number of optimized feature points is less than that of the initial candidate target points. And then we connect the optimized set of feature points into the polygon strips by the following rules. (1)For each feature point, connect with 4 adjacent feature points, shown in Figure 4, and 1, 2, 3, and 4 represent the direction coefficients.(2)Eliminate the edge, which the average of elevation of its two end points is lower, if there are two edges to be cross.

If there is edge from current point to the point in direction 2, at the same time there is edge from to the point in its direction 3, the cross phenomenon will appear, shown in Figure 4. The Rule (2) is a valid method to solve the above problems. Finally, the optimized feature segments composed by triangles and edges are obtained.

3. Ring Process

The most skeleton features of terrain are presented by tree topology. Hence, PPA breaks the polygon segments until there are not the closed polygons. Nevertheless, the real terrain contains a lot of atoll textures. As a result, tree topology is not enough to cover all of the skeleton features of real terrain. Therefore, we utilize graphic structure to store the skeleton features. However, not all of the rings are expected in the final result. Those single polygons and those sufficiently small rings should be disassembled, while big rings that can stand for the certain terrain feature are reserved by our method MPPA.

3.1. Polygon Disassembling

In MPPA, we take those rings that do not contain any nontarget points as single polygons. Figures 5(a)5(d) described the 4 kinds of triangles that need to be disassembled. Let be the current point; we just consider the triangles in the 4 directions as shown in Figure 4, since those in other directions can be processed by its left-top points.

For those triangles that accord with morphologies in Figure 5, we will break the edge with lower elevation. However, the undesirable parallelograms will appear, when two triangles shared a common edge, shown as dot lines in Figures 5(e)5(h). Hence, the edge shared by triangles will be remained (denoted by 2 in Figure 5), while the outer edges will be removed (denoted by 1 in Figure 5) in that case. The detailed procedure is described by Algorithm 2.

Input: Feature segments set FS with triangles and
  single edges.
Output: Feature segments set FS without any little triangles.
1. Traverse all of the feature segments:
 1.1 If there is not any case of 4 kinds of triangles of
  Figures 5(a)5(d) into set T, go to step 5; otherwise
  go to step 1.2;
 1.2 For , evaluate whether its three edges are shared;
 1.3 Denote the edge as 2, if it has been shared by
   two triangles; otherwise denote it as 1 and put it
   into set E.
2. Sort edges in set E in ascending order.
3. Remove the edge from FS and E, if
 satisfies the conditions as shown in Figures 5(a)5(d).
4. Go to step 1.
5. Algorithm end.

3.2. Judgment of Connected Domain

Besides breaking polygons, MPPA focuses on the ring process. In fact, there are few small rings in real terrain. It is reasonable that MPPA breaks the small rings but remains the big rings. Hence, it is important to evaluate the size of ring. The filling intensity is proposed in this paper to resolve this problem.

Definition 2. Let and be feature points set and nonfeature points set, respectively, and denote as a common label, if is enclosed in a closed area by the feature points of . The count of with the same label is called filling intensity.

By Definition 2, it is not difficult to find that filling intensity is the count of the interconnected points. Afterwards, we proposed the algorithm to judge the connected domain and calculate the filling intensity.

In Algorithm 3, each point is traversed once at most; therefore it does not bring big time overhead. Nevertheless, it is necessary to process the rings and reserve the long branches of skeleton features. As you can see, all of the nonfeature points have been assigned in the certain connected domain, parts of which have closed into rings, shown in Figure 6.

Input: Non-feature points set .
Output: The connected domain matrix DOM, and its
   filling intensity , .
1. Set , .
2. Set and DOM , ,
 where is the scale of original data.
3. Find one point that is neither target point, nor
 marked point from , let as the
  current non-feature point, and mark it, go to
  step 3.1; there is no such point go to step 4:
 3.1 Initialize array A;
 3.2 Set and ;
 3.3 Set and DOM ;
 3.4 Find one point that is neither target point, nor
    the marked point in eight adjacent points of
    current point, let be the point which is
    found at this time, if find the point go to
    step 3.6; if there is no such point, ,
     , go to step 3.5.
 3.5 If , go to step 4.
 3.6 Mark , , , = ,
     go to step 3.3.
4. Algorithm end.

3.3. Eliminate the Small Rings

For these small rings whose filling intensity is smaller than the criterion, MPPA consider it as the inauthentic terrain features. Hence, it is an important task to eliminate the sufficiently small rings. In Section 3.2, we have denoted the signs and calculated the filling intensity for each connected domain, which made it relative easy to remove the insignificant small rings.

By polygon breaking in Section 3.1, the basic skeleton features of terrain have come into being in substance. Afterwards, the main task is to denote those segments, as borders of small rings, and then delete them. The following rules are proposed to evaluate whether those segments should be eliminated. (1)Filling intensity of area where any adjacent point on both sides of the current segment belongs is less than the criterion.(2)Points on both sides of the current segment belong to different domains according to matrix.

Figure 7 described the distributing instances of points on both sides of the current segment, which bold line shows the current segment detected, green point indicates the start point of segment, and the yellow point indicates the end point of segment. The other points denote the detected points on both sides of the current segment. In order to distinguish the connected domain, we use to mark the different domains. If the above rules are satisfied, the current segment should be denoted and put them into a set. In terms of this strategy, those segments are traversed and processed in the same way. Finally, the segment with the lowest average elevation of its two end points should be eliminated. As a result, the undesirable small rings are disassembled.

3.4. Evaluate and Reserve Stable Big Rings

It is important to recognize accurate skeleton features with ring and prevent extracting wrong features; therefore not all the big rings should be reserved. It is a key problem to recognize the stable rings and remain them in the skeleton features. At the same time, unstable big rings still need to be disassembled.

By using the same method as in Section 3.2, we find out the big rings and determine the segments which belong to those big rings. By contrast with Rule (1) in Section 3.3, the opposite rule is adopted as the evaluation principle of big rings, that is to say, the filling intensity is bigger than the criterion. However, Rule (2) in Section 3.3 is also suitable for the judgment of big rings.

Nevertheless, these rules are not sufficient to judge the stability of big rings. An additional and significant rule is that the segments that comprised the big rings should be composed by the thick points. It is important to note that the above judgment is completed by original target points before the morphology simplification is performed. If the big ring contained some single segments composed by sparse target points, there is reason to believe that the big ring is unstable and this phenomenon is aroused by noise points.

In addition, the reasons for using the original target points to judge the density of segments on big rings are in order to avoid the effect of omitted noise points during the procedure of morphology simplification. Algorithm 4 described the detailed steps to evaluate and disassemble the unstable big rings. It is noteworthy that reduction of single branches around rings should have been done before Algorithm 4.

Input: Segments set segres without polygons and
  small rings.
Output: Segments set segbr and segments set
   segres that has been disassembled.
1. Traverse segments set segres. If all the segments
 have been traversed, then go to step 5.
2. Let as the current segment:
 2.1 Find the adjacent points of as Figure 7, and
   evaluate the filling intensities of these
   points, if , then go to step 1;
 2.2 Evaluate whether satisfies the Rule (2) in
   Section 3.3.
3. If the points (Figure 7) on the both sides of are
 non-target points according to the original target
 points set, and then:
 3.1 Disassemble the rings and remove from
 3.2 Otherwise, denote as the segment of the
    stable ring and put it into segbr.
4. If there is no any change for the segments set
 segres, then go to step 5.
5. Algorithm end.

4. Customizable Pruning of Skeleton Features

Although the skeleton features have been extracted by the above stages, there are some redundant little branches, shown in Figure 8. Moreover, the skeleton features that only contain the trunk or the branches with specified step length are usually required in many application fields. However, there are few extraction methods that can provide the personalized strategy of pruning.

4.1. Recognize Outbranch Points

In our MPPA, we proposed the customizable pruning strategy to solve the above problems. Since the branches whose step lengths are equal to one are bound to be undesirable, they should be pruned. In Figure 8, the blue lines described the little branches, whereas, those long branches, especially for main lines, should be remained, even by pruning for many times. Therefore, we proposed the method to recognize the outbranches and reserve them.

The premise of recognizing the outbranches is to find the branch points. In order to clearly express the concept of branch points, the formalized describing is as follows.

Definition 3. Let be a feature point, and let be the number of segments that are connected directly with ; then is also a branch point, if .

Based on Definition 3, the outbranch point is defined as follows.

Definition 4. Let be a branch point, and let be the number of branch points that are connected directly with ; then is also an outbranch point, if and only if .

As shown in Figure 9, the branches that are denoted by red squares are outbranches. It is very important to reserve the main feature lines. According to Definition 4, the outbranches must be located at the end of the skeleton features. If this kind of branches can be preserved, then the main feature lines will not be influenced during the period of the repeated pruning. Algorithm 5 described the procedure to recognize the outbranches.

Input: Array with branch points signs.
Output: Array with out-branching points signs.
1. Let as the count of the adjacent branch points of
 the current point.
2. For each branch point in , do:
 2.1 Set ;
 2.2 For each branch connected with , do:
  2.2.1 Find the end point of each branch
    connected with , according to the
    direction coefficients of ;
  2.2.2 If is an end point that has not any
    direction coefficient, then go to 2.2;
    otherwise, if is an out-branching
    point then set , and go to 2.2.
 2.3 If is equal to 1, denote as out-branch
   point in .
3. Algorithm end.

4.2. Prereserve Outbranches

In fact, the outbranching points are the border of feature lines. Moreover, the skeleton feature of real terrain is expected to be as long as possible and be consecutive. The main aim of recognizing outbranching points is to reserve the long and continuous feature lines. If all of branches are treated with by the uniform principle, then the feature lines will become shorter. Hence, we proposed the concept of the outbranch based on the outbranching point.

Definition 5. Let be an outbranch points set, and let be a branch in direction of end point of ; then is an outbranch, if is the longest branch in the end segments that connected with .

The red lines in Figure 10 showed the longest branch of each outbranching point. Those longest outbranches are preserved, so that these significant branches will not be influenced during the pruning, that is to say, the main trunk of skeleton features are still long enough, although most of insignificant branches have been eliminated.

4.3. Process of Pruning

By the process of rings and the preservation of outbranch mentioned before, our MPPA method can extract more precise skeleton features of real terrain. At the same time, it is a customizable method, which can determine the pruning length according to user’s different requirements. The detailed steps are presented in Algorithm 6.

Input: The current segments set segres, ring
  segments set segbr, and out branch sign
  array .
Output: The result segments set segres with the
   final skeleton features.
1. Let be branch points array produced in Section
 4.1 and let loop = 1 as loop times.
2. Call Algorithm 5 to update the out-branching set
of the current segments set segres.
3. Let lengthcriterion be criterion of branch length.
4. Search the unchecked directions of adjacent points
 of , until reach the end points of this
5. For each branch at the beginning of :
 5.1 Let ;
 5.2 Calculate the length of ;
 5.3 If , then eliminate and
   the branch ;
 5.4 Check direction coefficient of , if there is the
   direction that has not be traversed, then go to
   step 4;
 5.5 .
6. Retrieval the stable ring according to ring segments
 set segbr.
7. Retrieval the out branch according to
 out-branching sign array .
8. Algorithm end.

By Algorithm 6, the final personalized skeleton features have been extracted, as shown in Figure 11. As we can see, the final skeleton features have reserved the longer main ridge lines, and those insignificant little branches have been eliminated, as shown in Figure 11(a). Moreover, it is noteworthy that a stable ring has been remained in the final skeleton features, which are in accordance with the real terrain features in the natural world. It is not hard to find that the screenshot in Figure 11(b) is the known Crater Lake. Obviously, the ring skeleton feature has been extracted clearly.

5. Experimental Results and Discussion

To test validity of the results extracted by our MPPA, many experiments have been conducted. Since PPA is applied widely in synthesis and deformation of terrain [22, 27, 30], we compare our method MPPA with the classical method PPA. The experimental platform is a PC with a 3.0 GHz CPU, DDR3 1600 MHz 4 GB main memory, and Nvidia Geforce GTX 560 (1 GB) graphic card, and the experiment program is coded in VS2008.

Our MPPA is conducted based on GTOPO30, which contains the global elevation data. Its resolution is 30 arc second, that is, approximately 1 km. In order to compare MPPA with PPA, we download the benchmark data set ex100 and ex200 that is provided by PPA. There are 10000 points in ex100 and 40000 points in ex200.

Without loss of generality, we select the data L1T7, L1T15, and L1T16 randomly by our another roaming system of terrain visualization for GTOPO30. Their longitude and latitude are (97.458, 0.258), (41.458, 33.858), and (99.458, 31.458), respectively. Furthermore, the scale of data is about 16900 points. In addition, we also select G1, G2, and G3 to verify the effect of skeleton features extraction of MPPA. These three blocks of data will be described in Section 5.3.

In addition, Crater Lake has been used as our experimental data, which is of 10 m resolution with elevation points. In our experiment, Crater Lake is simplified as data set with points. The aim is to verify the effect of ring process of MPPA.

5.1. Performance Contrast Analysis

Whether for PPA or MPPA, profile recognition method is utilized to detect the target points. The main advantage is that it can cover all of terrain features because of its super loose judgment condition. Nevertheless, a lot of noise points will come into being.

The advantage of MPPA is that it first simplifies the points to obtain the subset of the points with more precise features, before those target points are connected into segments. However, PPA connected all of the target points to construct the polygon segments. Hence, the number of original segments of MPPA is less than that of PPA greatly. Table 1 described the number contrast of original segments between PPA and MPPA. This means that subsequent algorithms of MPPA are conducted on the little subset of segments, compared with the segments set of PPA. There is reason to believe that MPPA is of relatively high efficiency.

DEM dataNumber of segmentsSimplification ratio


Certainly, the fine efficiency of MPPA is at the price of the time overhead of simplification algorithm. In fact, the simplification algorithm of MPPA replaced with the polygon-breaking algorithm of PPA. Afterwards, we analyze and contrast the two methods.

MPPA proposed the simplification method of noise points based on the morphology. By defining the 48 kinds of morphological coding, the noise points are eliminated largely. From the perspective of simplification algorithm performance, each target point in MPPA need to be traversed only once to judge whether it should be reserved. If the feature morphology composed by the target points and its adjacent points satisfied anyone of the 48 kinds of morphologies, the target point will be eliminated.

Assume that the scale of elevation data is , and then on the surface, the simplification algorithm of MPPA needs to calculate the times. Nevertheless, our method does not need to match the 48 kinds of morphologies but calculate a morphological coding value and directly index the designed morphological signs. The signs have been built in advance by creating a space with 512 bits, equal to 64 bytes. It is really little space overhead. So, the calculating times of the simplification algorithm of MPPA is less than the original data scale . Obviously, the time complexity of the simplification method is .

Nevertheless, PPA starts with an end point of each segment to traverse all the branches that connected with it directly or indirectly, until return its end point of , or there is no path to go. It is noteworthy that the segments set of PPA is so big because it is constructed by fully connecting all of the original target points with a lot of noise points. And a mass of target points have been traversed for many times repeatedly. At worst, the time complexity of the polygon breaking method is . Even though the less points need to be traversed with process of the polygon breaking, its average time complexity is approximate .

Obviously, our simplification method can eliminate the noise points rapidly. The following chat presents the change of the noise point number with the increasing of iteration times, as shown in Figure 23.

5.2. Accuracy Analysis and Contrast

Except for the difference mentioned in Section 5.1, the effects of feature extraction by MPPA differ from those by PPA, which can be presented in the following three aspects.

First, since the principle of polygon breaking of PPA is that there is no closed polygon, it is impossible to extract the skeleton feature with rings, such as crater. Nevertheless, our method MPPA is accomplished in ring process. It can recognize the stable rings, reserve them, and at the same time eliminate unstable rings and small rings.

Furthermore, during the procedure of polygon breaking of PPA, it traverses all of the branches blindly and then reaches the other end point or backdates the jump-off point. Obviously, it does not consider any morphological information about skeleton features, which might lead to some error phenomenon. Figure 12 described the renderings contrast between PPA and MPPA. In order to present the visualized effect of feature extraction, we render the terrain model as the form of gray scale image, in which the light area shows the higher elevation, and the dark area shows the lower elevation.

As we can see, some area with lower elevation has been taken as skeleton feature by PPA, shown as the part B denoted by blue rectangle in Figure 12(a). It is not hard to see that the both sides of the part B are bright; therefore, this area is a valley, which should not be contained in the skeleton features. By MPPA, this area has been removed from the skeleton feature segments, shown in Figure 12(b). At the same time, the area A, denoted by white rectangle in Figure 12(a), ought to be reserved; however, PPA breaks the key feature segments.

The last but the most important, for some terrain elevation data, the pruning effect of MPPA is different from that of PPA. Since PPA connected all of points including a lot of noise points into polygon segments, there are so many insignificant little branches, after polygon breaking. Figure 13(a) showed the case of little branches, which might influence the effect of pruning. Whereas Figure 13(b) is the effect of feature segments connected by MPPA after target points simplification. It is not difficult to find that the trunk line extracted by the two methods is consistent, which also verified the validity of MPPA.

It is noteworthy that the feature segments of MPPA have few little branches. It is beneficial to prune accurately. Figure 14(b) showed the effect of MPPA after feature segments have been pruned, which is consistent with the effect of Figure 13(b). Nevertheless, as shown in Figure 14(a), PPA has omitted the trunk lines, which is inconsistent with its results before pruning, shown in Figure 13(a).

In our experiment, the step length of segment is 15, and that of segment is 5. Since the short branch should be pruned preferentially, ought to be eliminated. However, PPA has removed , which lost the accuracy of skeleton features. MPPA reserved and then connected it with segment ; therefore, a long truck line was formed. Obviously, these skeleton features that extracted by our method are more close to the real terrain features. Compared with PPA, MPPA reserved the long trunk lines and is of more precise.

5.3. Stability of Algorithm and Rendering Results

By executing the program of PPA, we find the various results for the same data. However, the results of MPPA are consistent though the experiment is executed for many times repeatedly. Three real elevation data are adopted to conduct the experiments, including ex100, ex200B0, and ex200B1. The results are shown in Figures 15, 16, and 17, where the two left screenshots are produced by PPA, and the right results are extracted by MPPA.

Apparently, these results from the three data sources demonstrate that many details of feature lines extracted by PPA are different even for the same data, which have been denoted by blue rectangle in Figures 1517. Hence, MPPA is of good stability, compared with PPA.

In order to further verify the effects of MPPA, we conduct MPPA on many real elevation data, including ex200B2, ex200B3, and other three blocks from GTOPO30.

Figure 18 described the three screenshots, including the original terrain, skeleton features, and longest truck line. About 10000 points are contained, and the longest truck feature has been denoted by blue line whose total length of steps is 202. Similarly, the results from data source ex200B3 have been presented in Figure 19. The length of the longest truck line is 210.

Other data sources are from GTOPO30, whose resolution is 1 km and the scale is to contain points. Their longitude and latitude are about , , and , and their screenshots are shown in Figures 20, 21, and 22. Obviously, the skeleton features with long truck lines are extracted well. The length of the longest branch is 201, 231, and 264, respectively, as shown in Figures 2022. By the above experiments and analysis for the experimental results, MPPA is of fine performance and accuracy of feature extraction.

6. Conclusions

At present, most features extraction methods of terrain skeleton are not competent to extract the ring features. Furthermore, the trunk line is not long enough to present the real terrain. Moreover, the accuracy of skeleton features for some terrain is not very high because of the negative effect of noisy points. Aimed at the above problems, a new method of skeleton feature extraction is proposed to eliminate the noise points, to reserve the truck line, and to recognize and remain the stable ring, which is significant in terrain modeling engineering.

In this paper, we defined the candidate feature points and proposed a morphological simplification algorithm to eliminate the noise points rapidly, which is benefit to the subsequent ring process and pruning. The feature point subset that has been simplified contains less data, compared with the original target points; therefore, the efficiency of ring process and pruning algorithms has been improved. Moreover, the accuracy of extraction has also been enhanced without the negative effect of noise points.

In addition, a concept of filling intensity is defined to determine a big stable ring. For small and unstable rings, we proposed the disassembling algorithm. In order to keep the long truck line, a concept of outbranch is defined, and the pruning algorithm is proposed to reserve the outbranch feature lines and eliminate the little branches.

Finally, all algorithms have been tested on many real elevation data and benchmark data. The experimental results show that our MPPA outweighs PPA in the aspects of performance and accuracy. Moreover, we proposed the strategy of ring process and customized pruning rules.

In the future, we will extract the skeleton features on the huge set of terrain elevation and expand these theories and algorithms of MPPA in the field of 3D visualization, such as 3D terrain synthesis and terrain deformation.


The authors are grateful for the helpful comments and suggestions of the reviewers, which have improved the presentation. This work is supported by Research Fund for the Doctoral Program of Higher Education of China (no. 20100043120012) and by National Natural Science Foundation of China for Young Scholars (no. 41101434).


  1. D. Bavera, C. de Michele, and M. Pepe, “Melted snow volume control in the snowmelt runoff model using a snow water equivalent statistically based model,” Hydrological Processes, vol. 26, no. 22, pp. 3405–3415, 2012. View at: Publisher Site | Google Scholar
  2. I. C. Fuller and M. Marden, “Slope-channel coupling in steepland terrain: a field-based conceptual model from the Tarndale gully and fan, Waipaoa catchment, New Zealand,” Geomorphology, vol. 128, no. 3-4, pp. 105–115, 2011. View at: Publisher Site | Google Scholar
  3. M. Cavalli, S. Trevisani, F. Comiti, and L. Marchi, “Geomorphometric assessment of spatial sediment connectivity in small Alpine catchments,” Geomorphology, vol. 188, pp. 31–41, 2013. View at: Google Scholar
  4. T. K. Peucker and D. H. Douglas, “Detection of surface-specific points by local parallel processing of discrete terrain elevation data,” Computer Graphics and Image Processing, vol. 4, no. 4, pp. 375–387, 1975. View at: Google Scholar
  5. J. Wood, The Geomorphological Characterisation of Digital Elevation Models, University of Leicester, Leicester, UK, 1996.
  6. L. E. Band, “Topographic partition of watersheds with digital elevation models,” Water Resources Research, vol. 22, no. 1, pp. 15–24, 1986. View at: Google Scholar
  7. J. F. O'Callaghan and D. M. Mark, “The extraction of drainage networks from digital elevation data,” Computer Vision, Graphics, & Image Processing, vol. 28, no. 3, pp. 323–344, 1984. View at: Google Scholar
  8. J. Fairfield and P. Leymarie, “Drainage networks from grid digital elevation models,” Water Resources Research, vol. 27, no. 5, pp. 709–717, 1991. View at: Publisher Site | Google Scholar
  9. M. C. Costa-Cabral and S. J. Burges, “Digital elevation model networks (DEMON): a model of flow over hillslopes for computation of contributing and dispersal areas,” Water Resources Research, vol. 30, no. 6, pp. 1681–1692, 1994. View at: Google Scholar
  10. D. G. Tarboton, “A new method for the determination of flow directions and upslope areas in grid digital elevation models,” Water Resources Research, vol. 33, no. 2, pp. 309–319, 1997. View at: Google Scholar
  11. A. Meisels, S. Raizman, and A. Karnieli, “Skeletonizing a DEM into a drainage network,” Computers and Geosciences, vol. 21, no. 1, pp. 187–196, 1995. View at: Google Scholar
  12. R. Moussa, M. Voltz, and P. Andrieux, “Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events,” Hydrological Processes, vol. 16, no. 2, pp. 393–412, 2002. View at: Publisher Site | Google Scholar
  13. C. Gascuel-Odoux, P. Aurousseau, T. Doray et al., “Incorporating landscape features to obtain an object-oriented landscape drainage network representing the connectivity of surface flow pathways over rural catchments,” Hydrological Processes, vol. 25, no. 23, pp. 3625–3636, 2011. View at: Publisher Site | Google Scholar
  14. F. Cazorzi, G. D. Fontana, A. D. Luca, G. Sofia, and P. Tarolli, “Drainage network detection and assessment of network storage capacity in agrarian landscape,” Hydrological Processes, vol. 27, no. 4, pp. 541–553, 2013. View at: Publisher Site | Google Scholar
  15. R. Barnes, C. Lehman, and D. Mulla, “An efficient assignment of drainage direction over flat surfaces in raster digital elevation models,” Computers & Geosciences, 2013. View at: Publisher Site | Google Scholar
  16. F. Nardi, S. Grimaldi, and M. Santini, “Hydrogeomorphic properties of simulated drainage patterns using digital elevation models: the flat area issue/Propriétés hydro-géomorphologiques de réseaux de drainage simulés à partir de modèles numériques de terrain: la question des zones planes,” Hydrological Sciences Journal, vol. 53, no. 6, pp. 1176–1193, 2008. View at: Google Scholar
  17. R. Jones, “Algorithms for using a DEM for mapping catchment areas of stream sediment samples,” Computers and Geosciences, vol. 28, no. 9, pp. 1051–1060, 2002. View at: Publisher Site | Google Scholar
  18. R. Jana, T. V. Reshmidevi, P. S. Arun, and T. I. Eldho, “An enhanced technique in construction of the discrete drainage network from low-resolution spatial database,” Computers and Geosciences, vol. 33, no. 6, pp. 717–727, 2007. View at: Publisher Site | Google Scholar
  19. H. Zhang and G. Huang, “Building channel networks for flat regions in digital elevation models,” Hydrological Processes, vol. 23, no. 20, pp. 2879–2887, 2009. View at: Publisher Site | Google Scholar
  20. Y. Chang, G. Song, and S. Hsu, “Automatic extraction of ridge and valley axes using the profile recognition and polygon-breaking algorithm,” Computers and Geosciences, vol. 24, no. 1, pp. 83–93, 1998. View at: Google Scholar
  21. H. Zhou, J. Sun, G. Turk, and J. M. Rehg, “Terrain synthesis from digital elevation models,” IEEE Transactions on Visualization and Computer Graphics, vol. 13, no. 4, pp. 834–848, 2007. View at: Publisher Site | Google Scholar
  22. H. Hnaidi, E. Guérin, S. Akkouche, A. Peytavie, and E. Galin, “Feature based terrain generation using diffusion equation,” Computer Graphics Forum, vol. 29, no. 7, pp. 2179–2186, 2010. View at: Publisher Site | Google Scholar
  23. O. Mora, J. J. Mallorqui, and A. Broquetas, “Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images,” IEEE Transactions on Geoscience and Remote Sensing, vol. 41, no. 10, pp. 2243–2253, 2003. View at: Publisher Site | Google Scholar
  24. J. Crause, A. Flower, and P. Marais, “A system for real-time deformable terrain,” in Proceedings of the 2011 ACM the South African Institute of Computer Scientists and Information Technologists Conference on Knowledge, Innovation and Leadership in a Diverse, Multidisciplinary Environment, pp. 77–86, October 2011. View at: Google Scholar
  25. R. Westerteiger, T. Compton, and T. Bemadin, “Interactive retro-deformation of rerrain for reconstructing 3D fault displacements,” IEEE Transactions on Visualization and Computer Graphics, vol. 18, no. 12, pp. 2208–2215, 2012. View at: Google Scholar
  26. R. Xu, F. Li, and L. Liao, “A 3D non-regular military symbol deformation method along the battlefield terrain,” in Proceedings of the 3rd IEEE International Conference on Computer Science and Information Technology (ICCSIT '10), pp. 207–211, July 2010. View at: Publisher Site | Google Scholar
  27. R. Kooima, J. Leigh, A. Johnson, D. Roberts, M. Subbarao, and T. A. Defanti, “Planetary-scale terrain composition,” IEEE Transactions on Visualization and Computer Graphics, vol. 15, no. 5, pp. 719–733, 2009. View at: Publisher Site | Google Scholar
  28. O. Št'ava, B. Beneš, M. Brisbin, and J. Křivánek, “Interactive terrain modeling using hydraulic erosion,” in Proceedings of the 2008 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pp. 201–210, 2008. View at: Google Scholar
  29. J. Gain, P. Marais, and W. Straßer, “Terrain sketching,” in Proceedings of the 2009 Symposium on Interactive 3D Graphics and Games, pp. 31–38, March 2009. View at: Publisher Site | Google Scholar
  30. F. P. Tasse, J. Gain, and P. Marais, “Enhanced texture-based terrain synthesis on graphics Hardware,” Computer Graphics Forum, vol. 31, no. 6, pp. 1959–1972, 2012. View at: Google Scholar
  31. X. Li, J. Wang, J. Zhou, and M. Yin, “A perturb biogeography based optimization with mutation for global numerical optimization,” Applied Mathematics and Computation, vol. 218, no. 2, pp. 598–609, 2011. View at: Publisher Site | Google Scholar
  32. X. Li and M. Yin, “Application of differential evolution algorithm on self-potential data,” PloS ONE, vol. 7, no. 12, Article ID e51199, 2012. View at: Google Scholar
  33. X. Li and M. Yin, “An opposition-based differential evolution Algorithm for permutation flow shop scheduling based on diversity measure,” Advances in Engineering Software, vol. 55, pp. 10–31, 2013. View at: Google Scholar

Copyright © 2013 Huijie Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

Article of the Year Award: Outstanding research contributions of 2020, as selected by our Chief Editors. Read the winning articles.