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

It is always difficul to reserve rings andmain truck lines in the real engineering of feature extraction for terrainmodel. In this paper, a new skeleton feature extraction method is proposed to solve these problems, which put forward a simplification algorithm based onmorphological theory to eliminate the noise points of the target points produced by classical profile recognition. Aswell 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.


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][2][3] or geometrical gradient analysis [4][5][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][9][10][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][14][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], Mathematical Problems in Engineering 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 Polygonbreaking 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][23][24][25][26][27][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.

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.

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 V , ∈  be current point, let  = ⌊  /2⌋ be the number of adjacent points for each profile, and let   , = {V ±Δ,±Δ | 0 ≤ Δ ≤ } be the set of adjacent points on the th profile; then V , ∈ V is a candidate feature point, if the following constraints are satisfied: (1) on all the profiles, there exist at least two points  1 ∈   , and  2 ∈   , ; their elevations are lower than that of V , ; (2)  1 and  2 should be at the different side on the same profile, where   is profile length, (, ) are coordinates of current points, and  = 1, 2, 3, 4 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.

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  = (  = 0 or 1 |  = 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: where (, ) is the spatial position of the elevation point V , .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   = ∑ 8 =0 2   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 2 0 ∼ 2 9 .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 MF (,) of the current point is calculated by where the value range of V +Δ,+Δ is 0 or 1 and The detail steps of simplification algorithm about candidate target points are described in Algorithm 1.
According to Algorithm 1, it is not hard to see that we take the morphology feature MF (,) 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.

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 V , to the point in direction 2, at the same time there is edge from V ,+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.

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 Mathematical Problems in Engineering    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.

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 V , 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.

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 V  ∈   as a common label, if V  is enclosed in a closed area by the feature points of   .The count of V  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.

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 DOM 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 DOM(, ) =   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.

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.

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.

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 V , be a feature point, and let  be the number of segments that are connected directly with V , ; then V , is also a branch point, if  > 2.
Based on Definition 3, the outbranch point is defined as follows.
Definition 4. Let V , be a branch point, and let  be the number of branch points that are connected directly with V , ; then V , is also an outbranch point, if and only if  = 1.
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.

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.

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.
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 Input: Segments set seg res without polygons and small rings.Output: Segments set seg br and segments set seg res that has been disassembled.1. Traverse segments set seg res .If all the segments have been traversed, then go to step 5. 2. Let  ∈ seg res as the current segment: 2.1 Find the adjacent points of  as Figure 7, and evaluate the filling intensities ins of these points, if ins < criterion, then go to step 1; 2.2 Evaluate whether  satisfies the Rule (2) in Section 3.  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.

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 4123×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.

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.
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 polygonbreaking 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 48 ×  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 ( 2 ).Even though the less points need to be traversed with process of the polygon breaking, its average time complexity is approximate ( 3/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 Figure 23.

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

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 15-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.
Figure 18 described the three screenshots, including the original terrain, skeleton features, and longest truck line.

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

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 seg res ; 3.2 Otherwise, denote  as the segment of the stable ring and put it into seg br .4. If there is no any change for the segments set seg res , then go to step 5. 5. Algorithm end.Algorithm 4: Evaluate and process the big rings.

Figure 12 :
Figure 12: Contrast of accuracy (the top right corner of ex200B3).
(b)  is the effect of feature segments connected by MPPA after target

Figure 13 :
Figure 13: Contrast of segments before pruning (the left of ex100).

Figure 14 :
Figure 14: Contrast of segments after pruning (the left of ex100).

Figure 23 :
Figure 23: Eliminating rate of noise points of MPPA.
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 V , ∈   , calculate its morphology feature MF (,) using formula (2); 3.2 Evaluate whether MF (,) matches one of 48 kinds of morphologies; 3.3 Eliminate V , from set   , if MTC[MF (,) ] equals to 1; 4. Algorithm end.Algorithm 1: Simplification of target points.
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).

Table 1 :
Contrast of segments between PPA and MPPA.