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

Authors: Huijie Zhang; Zhiqiang Ma; Yaxin Liu; Xinting He; Yun Ma
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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.


---

## Body

## 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 [1–3] or geometrical gradient analysis [4–6]. 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 2×2 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 [8–11]. 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 [13–15]. 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 [22–28], 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.Figure 1
Profile recognition.Let the number of points on each profile be the profile length aslp, which is closely related to recognition result of feature points. In this paper, we set lp 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.
LetV be an original set of elevation points, let vr,c∈V be current point, let m=⌊lp/2⌋ be the number of adjacent points for each profile, and let Pr,cn={vr±Δ,c±Δ∣0≤Δ≤m} be the set of adjacent points on the nth profile; then vr,c∈V is a candidate feature point, if the following constraints are satisfied:(1)
on all the profiles, there exist at least two pointsp1∈Pr,cn and p2∈Pr,cn; their elevations are lower than that of vr,c;(2)
p1andp2 should be at the different side on the same profile,wherelp is profile length, (r,c) are coordinates of current points, and n=1,2,3,4 is the sequence number of profile.By using Definition1, 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.Figure 2
Target points.


### 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.LetS=(bi=0or1∣i=0,1,2,…,8) be a binary coding sequence, and the mapping relationship between the bit code and the position of elevation point is presented by the following:
(1)r=⌊i3⌋+1,c=mod(i,3)+1,i=(r-1)*3+c-1,
where (r,c) is the spatial position of the elevation point vr,c.Obviously, it is easy to implement Formula (1). Nevertheless, we can determine the corresponding spatial position of each bit in S according to the mapping relationship. Moreover, there exists a corresponding decimal value dS=∑i=082bi for each S, and thus each dScan represent a kind of spatial morphology of terrain, named morphological coding. Since S consists of nine bits, the value range of dS is 20~29. Hence, morphological coding dS 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 (r,c) and dS. In Figure 3(b), each morphological coding stands for a kind of morphology, which has been denoted as dS. However there are only 12 types of morphologies illustrating in Figure 3(b). And the other ones can be deduced by rotating these instances.Morphological coding.
(a)
Binary sequence
(b)
A case of morphologiesFor 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 featureMF(r,c)of the current point is calculated by
(2)MF(r,c)=∑Δr=0,-1,1,Δc=0,-1,1vr+Δr,c+Δc*2i,
where the value range of vr+Δr,c+Δc is 0 or 1 and vr,c=1 if vr,c is a candidate target point; otherwise vr,c=0.The detail steps of simplification algorithm about candidate target points are described in Algorithm1.Algorithm 1:Simplification of target points.
Input: The setVf of candidate target points.
Output: The optimized setVf 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 inVf, then go to step 4;
otherwise do the following loop:
3.1 For each vr,c∈Vf, calculate its morphology
feature MF(r,c) using formula (2);
3.2 Evaluate whether MF(r,c) matches one of 48
kinds of morphologies;
3.3 Eliminate vr,c from set Vf, if MTC[MF(r,c)]
equals to 1;
4. Algorithm end.According to Algorithm1, it is not hard to see that we take the morphology feature MF(r,c) 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 Figure4, 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.Figure 4
Direction coefficient.If there is edge from current pointvr,c to the point in direction 2, at the same time there is edge from vr,c+1 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.


## 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.Figure 1
Profile recognition.Let the number of points on each profile be the profile length aslp, which is closely related to recognition result of feature points. In this paper, we set lp 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.
LetV be an original set of elevation points, let vr,c∈V be current point, let m=⌊lp/2⌋ be the number of adjacent points for each profile, and let Pr,cn={vr±Δ,c±Δ∣0≤Δ≤m} be the set of adjacent points on the nth profile; then vr,c∈V is a candidate feature point, if the following constraints are satisfied:(1)
on all the profiles, there exist at least two pointsp1∈Pr,cn and p2∈Pr,cn; their elevations are lower than that of vr,c;(2)
p1andp2 should be at the different side on the same profile,wherelp is profile length, (r,c) are coordinates of current points, and n=1,2,3,4 is the sequence number of profile.By using Definition1, 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.Figure 2
Target points.


## 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.LetS=(bi=0or1∣i=0,1,2,…,8) be a binary coding sequence, and the mapping relationship between the bit code and the position of elevation point is presented by the following:
(1)r=⌊i3⌋+1,c=mod(i,3)+1,i=(r-1)*3+c-1,
where (r,c) is the spatial position of the elevation point vr,c.Obviously, it is easy to implement Formula (1). Nevertheless, we can determine the corresponding spatial position of each bit in S according to the mapping relationship. Moreover, there exists a corresponding decimal value dS=∑i=082bi for each S, and thus each dScan represent a kind of spatial morphology of terrain, named morphological coding. Since S consists of nine bits, the value range of dS is 20~29. Hence, morphological coding dS 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 (r,c) and dS. In Figure 3(b), each morphological coding stands for a kind of morphology, which has been denoted as dS. However there are only 12 types of morphologies illustrating in Figure 3(b). And the other ones can be deduced by rotating these instances.Morphological coding.
(a)
Binary sequence
(b)
A case of morphologiesFor 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 featureMF(r,c)of the current point is calculated by
(2)MF(r,c)=∑Δr=0,-1,1,Δc=0,-1,1vr+Δr,c+Δc*2i,
where the value range of vr+Δr,c+Δc is 0 or 1 and vr,c=1 if vr,c is a candidate target point; otherwise vr,c=0.The detail steps of simplification algorithm about candidate target points are described in Algorithm1.Algorithm 1:Simplification of target points.
Input: The setVf of candidate target points.
Output: The optimized setVf 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 inVf, then go to step 4;
otherwise do the following loop:
3.1 For each vr,c∈Vf, calculate its morphology
feature MF(r,c) using formula (2);
3.2 Evaluate whether MF(r,c) matches one of 48
kinds of morphologies;
3.3 Eliminate vr,c from set Vf, if MTC[MF(r,c)]
equals to 1;
4. Algorithm end.According to Algorithm1, it is not hard to see that we take the morphology feature MF(r,c) 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 Figure4, 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.Figure 4
Direction coefficient.If there is edge from current pointvr,c to the point in direction 2, at the same time there is edge from vr,c+1 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. Figures5(a)–5(d) described the 4 kinds of triangles that need to be disassembled. Let vr,c 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.Triangles breaking.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)For those triangles that accord with morphologies in Figure5, 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.Algorithm 2:Polygon disassembling.
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
Figures5(a)–5(d) into set T, go to step 5; otherwise
go to step 1.2;
1.2 Fort∈T, 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 setE.
2. Sort edges in setE in ascending order.
3. Remove the edgee∈E from FS and E, if e
satisfies the conditions as shown in Figures5(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.
LetVf and Vn be feature points set and nonfeature points set, respectively, and denote vn∈Vn as a common label, if vn is enclosed in a closed area by the feature points of Vf. The count of vn with the same label is called filling intensity.By Definition2, 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 Algorithm3, 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.Algorithm 3:Connected domain judgment.
Input: Non-feature points setVn.
Output: The connected domain matrix DOM, and its
filling intensityInsi, i=1,2,3,….
1. SetInsi=0, i=1,2,3,….
2. Seti=0 and DOM(r,c)=0, 0≤r,c≤m,
wherem is the scale of original data.
3. Find one point that is neither target point, nor
marked point fromVn, let vr,c∈Vn as the
current non-feature point,i=i+1 and mark it, go to
step 3.1; there is no such point go to step 4:
3.1 Initialize arrayA;
3.2 Setk=1 and A[k]=vr,c;
3.3 SetInsi=Insi+1 and DOM(r,c)=i;
3.4 Find one point that is neither target point, nor
the marked point in eight adjacent points of
current point, letvrr,cc be the point which is
found at this time, if find the point go to
step 3.6; if there is no such point,k=k-1,
vr,c=A[k], go to step 3.5.
3.5 Ifk=1, go to step 4.
3.6 Markvrr,cc, k=k+1, A[k]=vrr,cc,vr,c = vrr,cc,
go to step 3.3.
4. Algorithm end.Figure 6
Connected domain.


### 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 Section3.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 Section3.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 toDOM matrix.Figure7 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 DOM(r,c)=ni 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.Evaluation of segment eliminating.
(a)
(b)
(c)
(d)


### 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 Section3.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. Algorithm4 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.Algorithm 4:Evaluate and process the big rings.
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. Lets∈segres as the current segment:
2.1 Find the adjacent points ofs as Figure 7, and
evaluate the filling intensitiesins of these
points, ifins<criterion, then go to step 1;
2.2 Evaluate whethers satisfies the Rule (2) in
Section3.3.
3. If the points (Figure7) on the both sides of s are
non-target points according to the original target
points set, and then:
3.1 Disassemble the rings and removes from
segres;
3.2 Otherwise, denotes 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.


## 3.1. Polygon Disassembling

In MPPA, we take those rings that do not contain any nontarget points as single polygons. Figures5(a)–5(d) described the 4 kinds of triangles that need to be disassembled. Let vr,c 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.Triangles breaking.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)For those triangles that accord with morphologies in Figure5, 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.Algorithm 2:Polygon disassembling.
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
Figures5(a)–5(d) into set T, go to step 5; otherwise
go to step 1.2;
1.2 Fort∈T, 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 setE.
2. Sort edges in setE in ascending order.
3. Remove the edgee∈E from FS and E, if e
satisfies the conditions as shown in Figures5(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.
LetVf and Vn be feature points set and nonfeature points set, respectively, and denote vn∈Vn as a common label, if vn is enclosed in a closed area by the feature points of Vf. The count of vn with the same label is called filling intensity.By Definition2, 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 Algorithm3, 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.Algorithm 3:Connected domain judgment.
Input: Non-feature points setVn.
Output: The connected domain matrix DOM, and its
filling intensityInsi, i=1,2,3,….
1. SetInsi=0, i=1,2,3,….
2. Seti=0 and DOM(r,c)=0, 0≤r,c≤m,
wherem is the scale of original data.
3. Find one point that is neither target point, nor
marked point fromVn, let vr,c∈Vn as the
current non-feature point,i=i+1 and mark it, go to
step 3.1; there is no such point go to step 4:
3.1 Initialize arrayA;
3.2 Setk=1 and A[k]=vr,c;
3.3 SetInsi=Insi+1 and DOM(r,c)=i;
3.4 Find one point that is neither target point, nor
the marked point in eight adjacent points of
current point, letvrr,cc be the point which is
found at this time, if find the point go to
step 3.6; if there is no such point,k=k-1,
vr,c=A[k], go to step 3.5.
3.5 Ifk=1, go to step 4.
3.6 Markvrr,cc, k=k+1, A[k]=vrr,cc,vr,c = vrr,cc,
go to step 3.3.
4. Algorithm end.Figure 6
Connected domain.


## 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 Section3.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 Section3.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 toDOM matrix.Figure7 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 DOM(r,c)=ni 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.Evaluation of segment eliminating.
(a)
(b)
(c)
(d)


## 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 Section3.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. Algorithm4 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.Algorithm 4:Evaluate and process the big rings.
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. Lets∈segres as the current segment:
2.1 Find the adjacent points ofs as Figure 7, and
evaluate the filling intensitiesins of these
points, ifins<criterion, then go to step 1;
2.2 Evaluate whethers satisfies the Rule (2) in
Section3.3.
3. If the points (Figure7) on the both sides of s are
non-target points according to the original target
points set, and then:
3.1 Disassemble the rings and removes from
segres;
3.2 Otherwise, denotes 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 Figure8. 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.Figure 8
Little branches.


### 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 Figure8, 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.
Letvr,c be a feature point, and let n be the number of segments that are connected directly with vr,c; then vr,c is also a branch point, if n>2.Based on Definition3, the outbranch point is defined as follows.Definition 4.
Letvr,c be a branch point, and let n be the number of branch points that are connected directly with vr,c; then vr,c is also an outbranch point, if and only if n=1.As shown in Figure9, 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.Algorithm 5:Recognize the outbranches.
Input: ArrayAb with branch points signs.
Output: ArrayAo with out-branching points signs.
1. Letn as the count of the adjacent branch points of
the current point.
2. For each branch pointpb in Ab, do:
2.1 Setn=0;
2.2 For each branch connected withpb, do:
2.2.1 Find the end point pa of each branch
connected with pb, according to the
direction coefficients of pb;
2.2.2 If pa is an end point that has not any
direction coefficient, then go to 2.2;
otherwise, if pa is an out-branching
point then set n=n+1, and go to 2.2.
2.3 Ifn is equal to 1, denote pb as out-branch
point inAo.
3. Algorithm end.Figure 9
Branch points.


### 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.
LetSb be an outbranch points set, and let eb be a branch in direction of end point of b∈Sb; then eb is an outbranch, if eb is the longest branch in the end segments that connected with b.The red lines in Figure10 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.Figure 10
Outbranching.


### 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 Algorithm6.Algorithm 6:Repeated pruning of MPPA.
Input: The current segments set segres, ring
segments set segbr, and out branch sign
array Ao.
Output: The result segments set segres with the
final skeleton features.
1. LetAb be branch points array produced in Section
4.1 and let loop = 1 as loop times.
2. Call Algorithm5 to update the out-branching set
Ao of the current segments set segres.
3. Let lengthcriterion be criterion of branch length.
4. Search the unchecked directions of adjacent points
ofp∈Ab, until reach the end points of this
branch.
5. For each brancheb at the beginning of p:
5.1 Letlengthcriterion=lengthcriterion+loop;
5.2 Calculate the lengthn of eb;
5.3 Ifn<lengthcriterion, then eliminate p and
the brancheb;
5.4 Check direction coefficient ofp, if there is the
direction that has not be traversed, then go to
step 4;
5.5loop=loop+1.
6. Retrieval the stable ring according to ring segments
set segbr.
7. Retrieval the out branch according to
out-branching sign arrayAo.
8. Algorithm end.By Algorithm6, 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.Results of skeleton features.
(a)
Ex100
(b)
Crater lake


## 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 Figure8, 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.
Letvr,c be a feature point, and let n be the number of segments that are connected directly with vr,c; then vr,c is also a branch point, if n>2.Based on Definition3, the outbranch point is defined as follows.Definition 4.
Letvr,c be a branch point, and let n be the number of branch points that are connected directly with vr,c; then vr,c is also an outbranch point, if and only if n=1.As shown in Figure9, 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.Algorithm 5:Recognize the outbranches.
Input: ArrayAb with branch points signs.
Output: ArrayAo with out-branching points signs.
1. Letn as the count of the adjacent branch points of
the current point.
2. For each branch pointpb in Ab, do:
2.1 Setn=0;
2.2 For each branch connected withpb, do:
2.2.1 Find the end point pa of each branch
connected with pb, according to the
direction coefficients of pb;
2.2.2 If pa is an end point that has not any
direction coefficient, then go to 2.2;
otherwise, if pa is an out-branching
point then set n=n+1, and go to 2.2.
2.3 Ifn is equal to 1, denote pb as out-branch
point inAo.
3. Algorithm end.Figure 9
Branch points.


## 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.
LetSb be an outbranch points set, and let eb be a branch in direction of end point of b∈Sb; then eb is an outbranch, if eb is the longest branch in the end segments that connected with b.The red lines in Figure10 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.Figure 10
Outbranching.


## 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 Algorithm6.Algorithm 6:Repeated pruning of MPPA.
Input: The current segments set segres, ring
segments set segbr, and out branch sign
array Ao.
Output: The result segments set segres with the
final skeleton features.
1. LetAb be branch points array produced in Section
4.1 and let loop = 1 as loop times.
2. Call Algorithm5 to update the out-branching set
Ao of the current segments set segres.
3. Let lengthcriterion be criterion of branch length.
4. Search the unchecked directions of adjacent points
ofp∈Ab, until reach the end points of this
branch.
5. For each brancheb at the beginning of p:
5.1 Letlengthcriterion=lengthcriterion+loop;
5.2 Calculate the lengthn of eb;
5.3 Ifn<lengthcriterion, then eliminate p and
the brancheb;
5.4 Check direction coefficient ofp, if there is the
direction that has not be traversed, then go to
step 4;
5.5loop=loop+1.
6. Retrieval the stable ring according to ring segments
set segbr.
7. Retrieval the out branch according to
out-branching sign arrayAo.
8. Algorithm end.By Algorithm6, 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.Results of skeleton features.
(a)
Ex100
(b)
Crater lake


## 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 Section5.3.In addition, Crater Lake has been used as our experimental data, which is of 10 m resolution with4123×4207 elevation points. In our experiment, Crater Lake is simplified as data set with 128×128 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. Table1 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.Table 1
Contrast of segments between PPA and MPPA.
DEM data
Number of segments
Simplification ratio
PPA
MPPA
ex100
8032
1968
75.50%
ex200
41374
11876
71.30%
L1T7
16208
6138
62.13%
L1T15
22928
12405
45.90%
L1T16
23608
8687
63.20%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 isn, and then on the surface, the simplification algorithm of MPPA needs to calculate the 48×n 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 n. Obviously, the time complexity of the simplification method is O(n).Nevertheless, PPA starts with an end point of each segments to traverse all the branches that connected with it directly or indirectly, until return its end point of s, 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 O(n2). Even though the less points need to be traversed with process of the polygon breaking, its average time complexity is approximate O(n3/2).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 Figure23.


### 5.2. Accuracy Analysis and Contrast

Except for the difference mentioned in Section5.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. Figure12 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.Contrast of accuracy (the top right corner of ex200B3).
(a)
PPA
(b)
MPPAAs 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 Figure12(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. Figure13(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.Contrast of segments before pruning (the left of ex100).
(a)
PPA
(b)
MPPAIt is noteworthy that the feature segments of MPPA have few little branches. It is beneficial to prune accurately. Figure14(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).Contrast of segments after pruning (the left of ex100).
(a)
PPA
(b)
MPPAIn our experiment, the step length of segmentoa is 15, and that of segment ob is 5. Since the short branch should be pruned preferentially, ob ought to be eliminated. However, PPA has removed oa, which lost the accuracy of skeleton features. MPPA reserved oa and then connected it with segment oc; 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 Figures15, 16, and 17, where the two left screenshots are produced by PPA, and the right results are extracted by MPPA.Contrast of stability (ex100).
(a)
One of results of PPA
(b)
Another result of PPA
(c)
MPPAContrast of stability (ex200B0).
(a)
One of results of PPA
(b)
Another result of PPA
(c)
MPPAContrast of stability (ex200B1).
(a)
One of results of PPA
(b)
Another result of PPA
(c)
MPPAApparently, 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 Figures15–17. 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.Figure18 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.Screenshots extracted by MPPA (right below of ex200B2).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchScreenshots extracted by MPPA (left below of ex200B3).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchOther data sources are from GTOPO30, whose resolution is 1 km and the scale is to contain93*75=6975 points. Their longitude and latitude are about G1(137.502,69.995), G2(138.278,70.620), and G3(140.628,67.495), 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 20–22. By the above experiments and analysis for the experimental results, MPPA is of fine performance and accuracy of feature extraction.Screenshots extracted by MPPA (G1).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchScreenshots extracted by MPPA (G2).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchScreenshots extracted by MPPA (G3).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchFigure 23
Eliminating rate of noise points 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. Table1 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.Table 1
Contrast of segments between PPA and MPPA.
DEM data
Number of segments
Simplification ratio
PPA
MPPA
ex100
8032
1968
75.50%
ex200
41374
11876
71.30%
L1T7
16208
6138
62.13%
L1T15
22928
12405
45.90%
L1T16
23608
8687
63.20%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 isn, and then on the surface, the simplification algorithm of MPPA needs to calculate the 48×n 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 n. Obviously, the time complexity of the simplification method is O(n).Nevertheless, PPA starts with an end point of each segments to traverse all the branches that connected with it directly or indirectly, until return its end point of s, 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 O(n2). Even though the less points need to be traversed with process of the polygon breaking, its average time complexity is approximate O(n3/2).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 Figure23.


## 5.2. Accuracy Analysis and Contrast

Except for the difference mentioned in Section5.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. Figure12 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.Contrast of accuracy (the top right corner of ex200B3).
(a)
PPA
(b)
MPPAAs 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 Figure12(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. Figure13(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.Contrast of segments before pruning (the left of ex100).
(a)
PPA
(b)
MPPAIt is noteworthy that the feature segments of MPPA have few little branches. It is beneficial to prune accurately. Figure14(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).Contrast of segments after pruning (the left of ex100).
(a)
PPA
(b)
MPPAIn our experiment, the step length of segmentoa is 15, and that of segment ob is 5. Since the short branch should be pruned preferentially, ob ought to be eliminated. However, PPA has removed oa, which lost the accuracy of skeleton features. MPPA reserved oa and then connected it with segment oc; 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 Figures15, 16, and 17, where the two left screenshots are produced by PPA, and the right results are extracted by MPPA.Contrast of stability (ex100).
(a)
One of results of PPA
(b)
Another result of PPA
(c)
MPPAContrast of stability (ex200B0).
(a)
One of results of PPA
(b)
Another result of PPA
(c)
MPPAContrast of stability (ex200B1).
(a)
One of results of PPA
(b)
Another result of PPA
(c)
MPPAApparently, 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 Figures15–17. 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.Figure18 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.Screenshots extracted by MPPA (right below of ex200B2).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchScreenshots extracted by MPPA (left below of ex200B3).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchOther data sources are from GTOPO30, whose resolution is 1 km and the scale is to contain93*75=6975 points. Their longitude and latitude are about G1(137.502,69.995), G2(138.278,70.620), and G3(140.628,67.495), 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 20–22. By the above experiments and analysis for the experimental results, MPPA is of fine performance and accuracy of feature extraction.Screenshots extracted by MPPA (G1).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchScreenshots extracted by MPPA (G2).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchScreenshots extracted by MPPA (G3).
(a)
Original Terrain
(b)
Skeleton feature results
(c)
Longest branchFigure 23
Eliminating rate of noise points of MPPA.


## 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.



---

*Source: 523729-2013-09-24.xml*