- Email: [email protected]
A procedure for automatically correcting invalid flat triangles occurring in triangulated contour data
PII:
Computers & Geosciences Vol. 24, No. 2, pp. 141±150, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0098-3004/98 $19.00 + 0.00 S0098-3004(97)00088-5
A PROCEDURE FOR AUTOMATICALLY CORRECTING INVALID FLAT TRIANGLES OCCURRING IN TRIANGULATED CONTOUR DATA J. MARK WARE University of Glamorgan, Pontypridd, Mid-Glamorgan, CF37 1DL, U.K. (e-mail: [email protected]) (Received 14 March 1997; accepted 1 June 1997) AbstractÐConstrained Delaunay triangulation is used widely as a basis for triangulating contour data. However, this method frequently produces triangulations that include ¯at triangles (i.e. where the three vertices of a triangle lie on the same contour). In many situations, such triangles fail to re¯ect the true shape of the surface being modelled. This paper presents details of a procedure, termed ¯at triangle corrector, which automatically corrects invalid ¯at triangles and produces a more accurate surface approximation. A detailed description of the procedure is provided in the form of a commentary and pseudocode listing. Experimental results on digitised contour data demonstrate that the approach is successful in correcting all ¯at triangles in a series of triangulations. # 1998 Elsevier Science Ltd. All rights reserved Key Words: Constrained Delaunay triangulation, Shelf triangles, Horizontal triangles, Surface modelling.
INTRODUCTION
Triangulating contour data using a conventional constrained Delaunay triangulation (CDT) approach will lead to solutions that include ¯at triangles (sometimes referred to as shelf or horizontal triangles). A ¯at triangle is de®ned as a triangle having its three vertices on the same contour. Flat triangles can exist either within single, closed contours or between pairs of adjacent contours (Fig. 1). In the situation of closed contours (Fig. 1A), ¯at triangles will sometimes be legitimate in that they represent areas where the real-world surface is truly ¯at (e.g. a lake). Alternatively, ¯at triangles occur in and around areas which represent peaks or troughs, resulting in truncated peaks and ¯attened troughs. Flat triangles that occur between pairs of adjacent contours (Fig. 1B) introduce regions of zero slope, thus introducing breaks, or arti®cial terraces, to the natural surface gradient. This paper presents a procedure, called ¯at triangle corrector (FTC), that seeks to eliminate all invalid ¯at triangles. The procedure is intended as a re®nement for use in conjunction with any existing CDT procedure. The development of FTC is justi®ed by the popularity of triangulation methods in GIS systems. Some previously suggested approaches to overcoming the ¯at triangle problem are reviewed. Next, details of a data structure suited to storing a CDT are provided. The FTC procedure, which takes the CDT (typically containing many ¯at triangles) as input, is then described and illus141
trated with pseudo code. Experimental results demonstrate the utility of the method on real-world data. PREVIOUS SOLUTIONS
Invalid ¯at triangles normally occur as a result of inappropriate data sampling. For example, the over-sampling of contours will often cause the average separating distance between adjacent points on individual contour lines to be much less than the average separating distance between adjacent contours. Because of the way in which constrained Delaunay triangulation works (i.e. points that are close together tend to be connected by triangle edges), this will, in general, lead to a higher probability of ¯at triangles occurring than if the average separating distance between points and the average separating distance between contours were more evenly matched. Inappropriate data sampling also occurs due to a failure to digitise and triangulate critical features, such as breaklines, ridges, valley ¯oors, spot heights, peaks and pits. Several approaches to overcoming the problem of ¯at triangles have been suggested previously. They can grouped into three classes: preventive measures, corrective procedures and speci®c solutions. Preventive measures One possible preventive approach to reducing the problem of ¯at triangles is to generalise the con-
142
J. M. Ware
Figure 1. Flat triangles (shaded) produced by constrained Delaunay triangulation. (A) Flat triangles within single, closed contour. (B) Flat triangles between adjacent pair of contours.
tours prior to triangulation (as suggested in E.S.R.I., 1991). This usually causes the average separating distance between points to become more equal to the average separating distance between contours, thus reducing the probability of ¯at triangles occurring. However, generalising contours has the detrimental eect of reducing the amount of contour detail. In extreme situations topological errors are introduced in the form of intersecting contours. Another preventive strategy is to ensure that all available critical features are included as constraints within the triangulation. The presence of such features will, in general, greatly reduce the chances of ¯at triangles occurring. The problem with this approach is that critical feature data are not always available. Also note that for both these preventive techniques, the removal of all ¯at triangles cannot be guaranteed. Corrective procedures Corrective procedures involve the swapping of edges and/or the insertion of new points subsequent to the initial triangulation. Brandli (1992) highlights two types of problem which can result from nonconstrained Delaunay triangulation of contour data, and provides corrective procedures. The ®rst situation is that of ¯at triangles. Such triangles are corrected by an interactive editor which, in addition to highlighting problem areas, allows for the swapping of edges and the insertion of new points. Initially, an attempt is made to eliminate ¯at regions, made up from one or more triangles, by swapping the edges which delimit them with a swap function. Given two triangles t1 (de®ned, say, by vertices {v1,v2,v3}) and t2 (de®ned {v1,v2,v4}), sharing a common edge {v1,v2}, this function updates triangulation topology such that t1 and t2 are subsequently de®ned {v2, v3, v4} and {v1, v3, v4} respectively (i.e. their common edge has been swapped and is now {v3,v4}). The swap function can only be applied to triangle pairs that, when combined, form a quadrilateral with no internal angle greater than or equal to 1808. In the event of an internal angle equalling or exceeding 1808, an
additional point must be inserted between the contours. The elevation of this point is set either by linear or cubic interpolation. The complete process is illustrated in Figure 2. The second situation described by Brandli is where individual triangle edges cross one or more contour lines. This problem, which does not arise when using constrained Delaunay triangulation, is overcome by the insertion of additional points. Vozenilek (1994) provides a corrective solution based solely on the insertion of what are referred to as ``help'' points. Initially, a CDT is created from contour lines, stream lines and spot heights. Flat triangles are then identi®ed, and their constituent edges extracted into a unique list of edges. Each edge is then inspected in turn. If an edge is not a contour edge (i.e. its endpoints are not neighbours on a contour) then the edge is classed as de®ning a problem region. The problem region is eliminated by adding a help point at the midpoint of the edge. Speci®c solutions The preventive measures and corrective procedures detailed previously are geared for use with conventional constrained and non-constrained Delaunay triangulation. There are, however, some non-Delaunay triangulation algorithms which have been developed speci®cally for triangulating contour data (for example, Fuchs, Keden and Uselton, 1977; Tipper, 1977; Christiansen and Sederberg, 1978; Ganthapy and Dennehy, 1982; Ekoule, Peyrin and Odet, 1991). A particular solution, of relevance here, is provided by Christensen (1987), who again highlights the contour-crossing edges and ¯at triangle problems associated with Delaunay triangulation. The ®rst of these problems is overcome by dividing the triangulation process up into a series of local triangulations performed between adjacent pairs of contours. The problem of ¯at triangles is solved by adding medial axis de®ning points between adjacent contour pairs. Medial axis points are given an elevation equal to the mean of the two contour elevations. The triangulation process obeys a simple
A procedure for automatically correcting invalid ¯at triangles
143
Figure 2. Correcting ¯at triangles. (A) Flat triangles (shaded). (B), (D), (E) and (F) Correction by means of edge swapping. (C) Correction by inserting point.
rule that ensures each triangle contains at least one contour point and at least one medial axis point, and thus ¯at triangles are avoided.
A TIN DATA STRUCTURE
This section gives a brief description of the TIN data structure adopted in the work presented in subsequent sections. Interested readers are referred to De Floriani (1987) and Woo (1985) for more comprehensive reviews of TIN data structures. A TIN can be considered as having three primary topological components, namely, vertices, edges and triangles. A data structure suited to encoding a TIN can be regarded as the combination of these basic components and a set of adjacency relations (De Floriani, 1987). This present paper assumes a triangle-based scheme in which each triangle is de®ned primarily by pointers to its three vertices and three adjacent triangles. This data structure is preferred to edge-based, point-based and alternative trianglebased structures because of its explicitness (i.e. triangles are explicitly de®ned, as are the neighbouring triangles of each triangle), making the FTC functions described in the next section more ecient and, incidentally, simpler to characterize. Additionally, each triangle de®nition records attribute information in the form of three Boolean ¯ags.
Each ¯ag relates to a particular edge of the triangle and indicates if that edge is a contour edge (¯ag is set to 1) or non-contour edge (¯ag is set to 0). The data structure is illustrated in Figure 3. Included in this ®gure are two C-like structure de®nitions which represent a triangulation of type triangulation as a linked-list of triangles of type triangle; this is the representation adopted in the functions described in the next section.
THE FLAT TRIANGLE CORRECTOR
This section gives details of the proposed FTC method for correcting invalid ¯at triangles. The procedure is applied to an existing CDT of the contour data. This initial triangulation is constructed using an algorithm based on the techniques of De Floriani and Puppo (1988), De Floriani (1989) and Heller (1990), and stored using the triangle-based data structure. Pseudo-code descriptions of the main FTC functions are provided the Appendix. Constructing a ¯at region Invalid ¯at triangles occur either in isolation or as contiguous sets. By visually studying triangulated contour data, it is observed that any invalid ¯at area, referred to as a Flat_Region, is bounded by two or more contour edges and, poss-
144
J. M. Ware
Figure 3. The triangle-based TIN data structure. (A) A sample CDT. (B) Triangle-based representation of sample triangulation. (C) C-like structure de®nition for storing CDT as linked-list of triangles.
ibly, by a single, non-¯at triangle tv made up of three non-contour edges. The absence of the non¯at triangle tv signi®es that Flat_Region lies within a single, closed contour; if tv exists then Flat_Region lies between an adjacent pair of contours (Fig. 4). The function FTC (FlatTriangle Corrector) involves a search for ¯at triangles through the triangulation T. A triangle is identi®ed as being ¯at if each of its vertices has the same elevation (z) value. Whenever a ¯at triangle t is found, the function ResolveFlat Triangle is called, taking t and T as its input$. ResolveFlatTriangle begins by ®nding Flat_Region. Making use of the function FindFlatRegion, the region is constructed by moving outward from t (using the adjacent triangle information associated with each triangle stored within T), adding ¯at triangles to Flat_Region as and when they are encountered. The movement outward along a particular path is terminated if either a contour edge is found or if a non-¯at triangle (tv) is encountered.
Having constructed Flat_Region, the function ResolveFlatTriangle proceeds as follows. If a non-¯at triangle tv has been found, then the function ResolveBetween Contours is called. Alternatively (i.e. tv=NULL), the function ResolveWithin Contour is invoked. Correcting between a pair of contours ResolveBetweenContours closely follows the approach suggested by Brandli (1992) in that an initial attempt is made to correct Flat_Region simply by edge swapping; if this does not succeed fully, additional points are added. To begin, the function FindAdjacent Triangle is called. This function, taking triangle tv, de®ned {v1,v2,v3}, and the list of triangles Flat_Region as its input, searches through Flat_Region looking for the triangle with which tv is adjacent (of which there will always be one, and one only). This triangle, ta, de®ned {v1,v2,v4}, is removed from Flat_Region and is subsequently returned to the calling procedure. It now follows that provided
A procedure for automatically correcting invalid ¯at triangles
145
Figure 4. (A) In situation of adjacent contour pair, Flat_Region is bounded on one side by contour edges (e1, e2, e3 and e4) and on other by non-¯at triangle tv. (B) For single, closed contour, Flat_Region is bounded only by contour edges (e1, e2, e3, e4, e5, e6 and e7).
the quadrilateral (Q) formed by tv and ta is convex (i.e. the maximum internal angle is less than 1808), then ta can be corrected by swapping the edge which it shares with tv. This is achieved by calling SwapCommonEdge and results in two non-¯at triangles tv and ta, de®ned {v2, v3, v4} and {v1, v3, v4} respectively. One of these triangles will have three non-contour edges; Swap CommonEdge ensures that it is always tv. In this way the process of ®nding ta and swapping the common edge can be repeated until eventually, provided Q is always convex, Flat_Region is emptied. If at any stage a non-convex quadrilateral Q is formed by tv and ta, it signi®es the need to insert a new point. This point, p, is positioned at the midpoint of the edge common to tv and ta, and is assigned an elevation using cubic interpolation. New points are added using AddPointToTriangulation. This function is based on the point insertion process described by De Floriani (1989). Having inserted p, some or all of the triangles in Flat_Region will have been corrected. In either situation, ResolveBetweenContours is exited; any ¯at triangles which remain are guaranteed to be encountered at a later stage, either during the main
FTC search or during the construction of a subsequent Flat_Region. Note that this edge-swapping approach is relatively ecient in limiting the number of new points inserted (see Table 1). Alternative methods, based solely on insertion of points, similar to those suggested by Vozenilek (1994) and Christensen (1987), could easily be implemented. Evaluation of which approach is in fact better will be the subject of future qualitative assessment experiments. Correcting within a single contour If a non-¯at triangle is not encountered, ResolveWithinContour is called. This function involves the addition of a single new point p at the centre of Flat_Region. The centre of Flat_Region is de®ned here as being the centre of its largest possible inscribed circle. As with ResolveBetweenContours, insertion of p will correct some or all of the triangles in Flat_Region; any ¯at triangles which remain are encountered at a later stage. It should be noted that FTC as described here does not dierentiate between closed contour situations which involve legitimate ¯at triangles (e.g. a lake) and those where ¯at triangles represent error (e.g. a truncated peak). In both situations, FTC
Table 1. Details of experimental data, CDTs prior to application of FTC and CDTs subsequent to application of FTC CDT prior to application CDT subsequent to application of FTC of FTC
Source data Data set 1 2 3 4
No. of contours 57 74 36 43
No. of points 4484 5562 2928 3384
No. of edges 4433 5494 2890 3337
Total no. of No. of ¯at Total no. of No. of points No. of ¯at triangles triangles triangles inserted triangles 8962 11118 5850 6766
1019 1059 1529 1048
9000 11162 5920 6970
19 32 35 14
0 0 0 0
146
J. M. Ware
Figure 7. Part of CDT generated from data set 2. (A) Before application of FTC. (B) After application of FTC.
Figure 5. Part of CDT generated from data set 1. (A) Before application of FTC. (B) After application of FTC.
Figure 6. Part of CDT generated from data set 1. (A) Before application of FTC. (B) After application of FTC.
Figure 8. Part of CDT generated from data set 3. (A) Before application of FTC. (B) After application of FTC.
A procedure for automatically correcting invalid ¯at triangles
147
Figure 9. Part of CDT generated from data set 3. (A) Before application of FTC. (B) After application of FTC.
Figure 11. Part of CDT generated from data set 4. (A) Before application of FTC. (B) After application of FTC.
Figure 10. Part of CDT generated from data set 3. (A) Before application of FTC. (B) After application of FTC.
Figure 12. Part of CDT generated from data set 4. (A) Before application of FTC. (B) After application of FTC.
148
J. M. Ware
attempts to correct the triangulation by calling ResolveWithinContour. In order to avoid an attempt at correcting a legitimate Flat_Region, it is necessary to have access to contour attribute data giving detail of what classes of geographical phenomena lie adjacent to individual contours.
EXPERIMENTAL RESULTS
The FTC procedure has been implemented in C on a Sparc workstation. Tests have been carried on CDTs generated from 1:10 000 scale UK Ordnance Survey digitised contour data. Table 1 gives details relating to source data, CDTs and experimental results. Note that in each experiment all ¯at triangles were successfully corrected. Parts of the triangulations, prior to (A) and subsequent to (B) application of FTC, are illustrated in Figures 5±12.
FINAL REMARKS
The paper outlines a procedure, FTC, for automatically correcting ¯at triangles occurring in triangulations generated from contour data. The procedure, having been tested on a number of data sets, has been shown to correct all ¯at triangles successfully in a given triangulation. The procedure, as it stands, provides users with a useful tool for improving models produced by existing CDT algorithms. Future work will concentrate on producing a qualitative assessment of the FTC approach. This paper has assumed that correcting ¯at triangles does indeed lead to an improved triangulated model; this assumption needs to be validated by, for example, a comprehensive visibility analysis study (see Kidner, 1996; Ware, Kidner and Herbert, 1996; for related work). Such a study will measure the eectiveness of Delaunay, constrainedDelaunay, non-Delaunay and FTC methods by comparing predicted results with ®eld observations, and analyse the alternative FTC edge-swapping and point-insertion strategies discussed here. A range of interpolation methods, including the cubic technique currently used, will be tested to estimate the elevation value for new points. AcknowledgmentsÐThe digital contour data used in experiments described in this paper has been provided by the UK Ordnance Survey.
REFERENCES Brandli, M. (1992) A triangulation-based method for geomorphological surface interpolation from contour lines. Proceedings 3rd European Conference and Exhibition on GIS (EGIS `92), Vol. 1, Munich, pp. 691±700. Christensen, A. H. J. (1987) Fitting a triangulation to contour lines. Proceedings 8th International Symposium of Computer Assisted Cartography, Baltimore, Maryland, pp. 57±67. Christiansen, H. N. and Sederberg, T. W. (1978) Conversion of complex contour line de®nitions into polygonal element mosaics. ACM Computer Graphics 12(3), 187±192. De Floriani, L. (1987) Surface representations based on triangular grids. Visual Computer 3(1), 27±50. De Floriani, L. (1989) A pyramidal data structure for triangle-based surface description. IEEE Computer Graphics and Applications, March 1989, pp. 67±78. De Floriani, L. and Puppo, E. (1988) Constrained Delaunay triangulation for multiresolution surface description. Proceedings 9th IEEE International Conference on Pattern Recognition, CS Press, Rome, 1988, pp. 566±569. Ekoule, A. B., Peyrin, F. C. and Odet, L. (1991) A triangulation algorithm from arbitrary shaped multiple planar contours. ACM Transactions on Graphics 10(2), 182±199. E.S.R.I. (1991) Surface modelling with TIN: surface analysis and display, ARC/INFO User's Guide. Environmental Systems Research Institute, Redlands, California. Fuchs, H., Kedem, Z. M. and Uselton, S. P. (1977) Optimal surface reconstruction from planar contours. Communications of the ACM 20(10), 693±702. Ganthapy, S. and Dennehy, T. G. (1982) A new general triangulation method for planar contours. ACM Computer Graphics 13(3), 311±319. Heller, M. (1990) Triangulation algorithms for adaptive terrain modelling. Proceedings 4th International Symposium on Spatial Data Handling, Vol. 1, pp. 163± 174. Kidner, D. B. (1996) Site selection and visibility analysis for a wind farm development: a problem for GIS? Proceedings 1st International Conference on GIS in Urban, Regional and Environmental Planning, Samos, Greece, pp. 220±237. Tipper, J. C. (1977) A method and Fortran program for the computerised reconstruction of three dimensional objects from serial sections. Computers & Geosciences 3(4), 579±599. Vozenilek, V. (1994) Generating surface models using elevations digitised from topographical maps. Proceedings 5th European Conference and Exhibition on GIS, Vol. 1, Utrecht, pp. 972±982. Ware, J. M., Kidner, D. B. and Herbert, M. J. (1996) Terrain in perspective: DEMs or TINs. Proceedings 2nd Joint European Conference on Geographic Information, Barcelona, pp. 388±397. Woo, T. C. (1985) A combinational analysis of boundary data structure schemata. IEEE Computer Graphics and Applications 5(3), 19±27.
A procedure for automatically correcting invalid ¯at triangles APPENDIX Pseudo-code for FTC struct triangulation {t:triangle_pointer; next:triangulation_pointer } struct triangle {vertex[1..3]:vertex_pointer; adjacent[1..3]:triangle_pointer; edge_¯ag[1..3]:Boolean } FTC(IN - T:triangulation, T_head:triangulation_pointer) current = T_head While current$NULL If (FlatTriangle(current.t)) ResolveFlatTriangle(T,T_head,current.t) Else current = current.next Endif Endwhile End FTC ResolveFlatTriangle(IN - T:triangulation, T_head:triangulation_pointer, t:triangle) FindFlatRegion(T, t, Flat_Region, head, tv) If (tv $ NULL) ResolveBetweenContours(T,T_head,tv) Else ResolveWithinContour(T, tv) Endif End ResolveFlatTriangle FindFlatRegion(IN OUT -
T:triangulation, t:triangle; Flat_Region:triangulation, FR_head:triangulation_pointer, tv:triangle)
tv = NULL FR_head = NULL TQ_head = NULL Add(t,Triangle_Queue,TQ_head) current = TQ_head While (current$NULL) ta = current.t Add(ta,Flat_Region,FR_head) For i =1 to 3 If (ta.edge_¯ag[i] =0) If (FlatTriangle(ta.adjacent[i])) If (ta.adjacent[i] ( Triangle_Queue) and (ta.adjacent[i] ( Flat_Region) Add(ta.adjacent[i],Triangle_Queue,TQ_head) Endif Else t v = ta Endif Endif Endfor Endwhile End FindFlatRegion
149
150
J. M. Ware
ResolveBetweenContours(IN - T, Flat_Region:triangulation, T_head,FR_head:triangulation_pointer, tv:triangle) Swapped = True While (Swapped) and (FR_head$NULL) ta = FindAdjacentTriangle(tv, Flat_Region) Swapped = SwapCommonEdge(tv,ta)) If (Swapped) Delete(ta,Flat_Region) t v = ta Else p = FindMidPoint(tv,ta) InterpolateElevation(p) AddPointToTriangulation(p,T,T_Head) Endif Endwhile End ResolveBetweenContours ResolveWithinContour(IN -
T, Flat_Region:triangulation, T_head,FR_head:triangulation_pointer) p = FindCentreOfRegion(Flat_Region, FR_head) InterpolateElevation(p) AddPointToTriangulation(p,T,T_Head) End ResolveWithinContour
Computers & Geosciences Vol. 24, No. 2, pp. 141±150, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0098-3004/98 $19.00 + 0.00 S0098-3004(97)00088-5
A PROCEDURE FOR AUTOMATICALLY CORRECTING INVALID FLAT TRIANGLES OCCURRING IN TRIANGULATED CONTOUR DATA J. MARK WARE University of Glamorgan, Pontypridd, Mid-Glamorgan, CF37 1DL, U.K. (e-mail: [email protected]) (Received 14 March 1997; accepted 1 June 1997) AbstractÐConstrained Delaunay triangulation is used widely as a basis for triangulating contour data. However, this method frequently produces triangulations that include ¯at triangles (i.e. where the three vertices of a triangle lie on the same contour). In many situations, such triangles fail to re¯ect the true shape of the surface being modelled. This paper presents details of a procedure, termed ¯at triangle corrector, which automatically corrects invalid ¯at triangles and produces a more accurate surface approximation. A detailed description of the procedure is provided in the form of a commentary and pseudocode listing. Experimental results on digitised contour data demonstrate that the approach is successful in correcting all ¯at triangles in a series of triangulations. # 1998 Elsevier Science Ltd. All rights reserved Key Words: Constrained Delaunay triangulation, Shelf triangles, Horizontal triangles, Surface modelling.
INTRODUCTION
Triangulating contour data using a conventional constrained Delaunay triangulation (CDT) approach will lead to solutions that include ¯at triangles (sometimes referred to as shelf or horizontal triangles). A ¯at triangle is de®ned as a triangle having its three vertices on the same contour. Flat triangles can exist either within single, closed contours or between pairs of adjacent contours (Fig. 1). In the situation of closed contours (Fig. 1A), ¯at triangles will sometimes be legitimate in that they represent areas where the real-world surface is truly ¯at (e.g. a lake). Alternatively, ¯at triangles occur in and around areas which represent peaks or troughs, resulting in truncated peaks and ¯attened troughs. Flat triangles that occur between pairs of adjacent contours (Fig. 1B) introduce regions of zero slope, thus introducing breaks, or arti®cial terraces, to the natural surface gradient. This paper presents a procedure, called ¯at triangle corrector (FTC), that seeks to eliminate all invalid ¯at triangles. The procedure is intended as a re®nement for use in conjunction with any existing CDT procedure. The development of FTC is justi®ed by the popularity of triangulation methods in GIS systems. Some previously suggested approaches to overcoming the ¯at triangle problem are reviewed. Next, details of a data structure suited to storing a CDT are provided. The FTC procedure, which takes the CDT (typically containing many ¯at triangles) as input, is then described and illus141
trated with pseudo code. Experimental results demonstrate the utility of the method on real-world data. PREVIOUS SOLUTIONS
Invalid ¯at triangles normally occur as a result of inappropriate data sampling. For example, the over-sampling of contours will often cause the average separating distance between adjacent points on individual contour lines to be much less than the average separating distance between adjacent contours. Because of the way in which constrained Delaunay triangulation works (i.e. points that are close together tend to be connected by triangle edges), this will, in general, lead to a higher probability of ¯at triangles occurring than if the average separating distance between points and the average separating distance between contours were more evenly matched. Inappropriate data sampling also occurs due to a failure to digitise and triangulate critical features, such as breaklines, ridges, valley ¯oors, spot heights, peaks and pits. Several approaches to overcoming the problem of ¯at triangles have been suggested previously. They can grouped into three classes: preventive measures, corrective procedures and speci®c solutions. Preventive measures One possible preventive approach to reducing the problem of ¯at triangles is to generalise the con-
142
J. M. Ware
Figure 1. Flat triangles (shaded) produced by constrained Delaunay triangulation. (A) Flat triangles within single, closed contour. (B) Flat triangles between adjacent pair of contours.
tours prior to triangulation (as suggested in E.S.R.I., 1991). This usually causes the average separating distance between points to become more equal to the average separating distance between contours, thus reducing the probability of ¯at triangles occurring. However, generalising contours has the detrimental eect of reducing the amount of contour detail. In extreme situations topological errors are introduced in the form of intersecting contours. Another preventive strategy is to ensure that all available critical features are included as constraints within the triangulation. The presence of such features will, in general, greatly reduce the chances of ¯at triangles occurring. The problem with this approach is that critical feature data are not always available. Also note that for both these preventive techniques, the removal of all ¯at triangles cannot be guaranteed. Corrective procedures Corrective procedures involve the swapping of edges and/or the insertion of new points subsequent to the initial triangulation. Brandli (1992) highlights two types of problem which can result from nonconstrained Delaunay triangulation of contour data, and provides corrective procedures. The ®rst situation is that of ¯at triangles. Such triangles are corrected by an interactive editor which, in addition to highlighting problem areas, allows for the swapping of edges and the insertion of new points. Initially, an attempt is made to eliminate ¯at regions, made up from one or more triangles, by swapping the edges which delimit them with a swap function. Given two triangles t1 (de®ned, say, by vertices {v1,v2,v3}) and t2 (de®ned {v1,v2,v4}), sharing a common edge {v1,v2}, this function updates triangulation topology such that t1 and t2 are subsequently de®ned {v2, v3, v4} and {v1, v3, v4} respectively (i.e. their common edge has been swapped and is now {v3,v4}). The swap function can only be applied to triangle pairs that, when combined, form a quadrilateral with no internal angle greater than or equal to 1808. In the event of an internal angle equalling or exceeding 1808, an
additional point must be inserted between the contours. The elevation of this point is set either by linear or cubic interpolation. The complete process is illustrated in Figure 2. The second situation described by Brandli is where individual triangle edges cross one or more contour lines. This problem, which does not arise when using constrained Delaunay triangulation, is overcome by the insertion of additional points. Vozenilek (1994) provides a corrective solution based solely on the insertion of what are referred to as ``help'' points. Initially, a CDT is created from contour lines, stream lines and spot heights. Flat triangles are then identi®ed, and their constituent edges extracted into a unique list of edges. Each edge is then inspected in turn. If an edge is not a contour edge (i.e. its endpoints are not neighbours on a contour) then the edge is classed as de®ning a problem region. The problem region is eliminated by adding a help point at the midpoint of the edge. Speci®c solutions The preventive measures and corrective procedures detailed previously are geared for use with conventional constrained and non-constrained Delaunay triangulation. There are, however, some non-Delaunay triangulation algorithms which have been developed speci®cally for triangulating contour data (for example, Fuchs, Keden and Uselton, 1977; Tipper, 1977; Christiansen and Sederberg, 1978; Ganthapy and Dennehy, 1982; Ekoule, Peyrin and Odet, 1991). A particular solution, of relevance here, is provided by Christensen (1987), who again highlights the contour-crossing edges and ¯at triangle problems associated with Delaunay triangulation. The ®rst of these problems is overcome by dividing the triangulation process up into a series of local triangulations performed between adjacent pairs of contours. The problem of ¯at triangles is solved by adding medial axis de®ning points between adjacent contour pairs. Medial axis points are given an elevation equal to the mean of the two contour elevations. The triangulation process obeys a simple
A procedure for automatically correcting invalid ¯at triangles
143
Figure 2. Correcting ¯at triangles. (A) Flat triangles (shaded). (B), (D), (E) and (F) Correction by means of edge swapping. (C) Correction by inserting point.
rule that ensures each triangle contains at least one contour point and at least one medial axis point, and thus ¯at triangles are avoided.
A TIN DATA STRUCTURE
This section gives a brief description of the TIN data structure adopted in the work presented in subsequent sections. Interested readers are referred to De Floriani (1987) and Woo (1985) for more comprehensive reviews of TIN data structures. A TIN can be considered as having three primary topological components, namely, vertices, edges and triangles. A data structure suited to encoding a TIN can be regarded as the combination of these basic components and a set of adjacency relations (De Floriani, 1987). This present paper assumes a triangle-based scheme in which each triangle is de®ned primarily by pointers to its three vertices and three adjacent triangles. This data structure is preferred to edge-based, point-based and alternative trianglebased structures because of its explicitness (i.e. triangles are explicitly de®ned, as are the neighbouring triangles of each triangle), making the FTC functions described in the next section more ecient and, incidentally, simpler to characterize. Additionally, each triangle de®nition records attribute information in the form of three Boolean ¯ags.
Each ¯ag relates to a particular edge of the triangle and indicates if that edge is a contour edge (¯ag is set to 1) or non-contour edge (¯ag is set to 0). The data structure is illustrated in Figure 3. Included in this ®gure are two C-like structure de®nitions which represent a triangulation of type triangulation as a linked-list of triangles of type triangle; this is the representation adopted in the functions described in the next section.
THE FLAT TRIANGLE CORRECTOR
This section gives details of the proposed FTC method for correcting invalid ¯at triangles. The procedure is applied to an existing CDT of the contour data. This initial triangulation is constructed using an algorithm based on the techniques of De Floriani and Puppo (1988), De Floriani (1989) and Heller (1990), and stored using the triangle-based data structure. Pseudo-code descriptions of the main FTC functions are provided the Appendix. Constructing a ¯at region Invalid ¯at triangles occur either in isolation or as contiguous sets. By visually studying triangulated contour data, it is observed that any invalid ¯at area, referred to as a Flat_Region, is bounded by two or more contour edges and, poss-
144
J. M. Ware
Figure 3. The triangle-based TIN data structure. (A) A sample CDT. (B) Triangle-based representation of sample triangulation. (C) C-like structure de®nition for storing CDT as linked-list of triangles.
ibly, by a single, non-¯at triangle tv made up of three non-contour edges. The absence of the non¯at triangle tv signi®es that Flat_Region lies within a single, closed contour; if tv exists then Flat_Region lies between an adjacent pair of contours (Fig. 4). The function FTC (FlatTriangle Corrector) involves a search for ¯at triangles through the triangulation T. A triangle is identi®ed as being ¯at if each of its vertices has the same elevation (z) value. Whenever a ¯at triangle t is found, the function ResolveFlat Triangle is called, taking t and T as its input$. ResolveFlatTriangle begins by ®nding Flat_Region. Making use of the function FindFlatRegion, the region is constructed by moving outward from t (using the adjacent triangle information associated with each triangle stored within T), adding ¯at triangles to Flat_Region as and when they are encountered. The movement outward along a particular path is terminated if either a contour edge is found or if a non-¯at triangle (tv) is encountered.
Having constructed Flat_Region, the function ResolveFlatTriangle proceeds as follows. If a non-¯at triangle tv has been found, then the function ResolveBetween Contours is called. Alternatively (i.e. tv=NULL), the function ResolveWithin Contour is invoked. Correcting between a pair of contours ResolveBetweenContours closely follows the approach suggested by Brandli (1992) in that an initial attempt is made to correct Flat_Region simply by edge swapping; if this does not succeed fully, additional points are added. To begin, the function FindAdjacent Triangle is called. This function, taking triangle tv, de®ned {v1,v2,v3}, and the list of triangles Flat_Region as its input, searches through Flat_Region looking for the triangle with which tv is adjacent (of which there will always be one, and one only). This triangle, ta, de®ned {v1,v2,v4}, is removed from Flat_Region and is subsequently returned to the calling procedure. It now follows that provided
A procedure for automatically correcting invalid ¯at triangles
145
Figure 4. (A) In situation of adjacent contour pair, Flat_Region is bounded on one side by contour edges (e1, e2, e3 and e4) and on other by non-¯at triangle tv. (B) For single, closed contour, Flat_Region is bounded only by contour edges (e1, e2, e3, e4, e5, e6 and e7).
the quadrilateral (Q) formed by tv and ta is convex (i.e. the maximum internal angle is less than 1808), then ta can be corrected by swapping the edge which it shares with tv. This is achieved by calling SwapCommonEdge and results in two non-¯at triangles tv and ta, de®ned {v2, v3, v4} and {v1, v3, v4} respectively. One of these triangles will have three non-contour edges; Swap CommonEdge ensures that it is always tv. In this way the process of ®nding ta and swapping the common edge can be repeated until eventually, provided Q is always convex, Flat_Region is emptied. If at any stage a non-convex quadrilateral Q is formed by tv and ta, it signi®es the need to insert a new point. This point, p, is positioned at the midpoint of the edge common to tv and ta, and is assigned an elevation using cubic interpolation. New points are added using AddPointToTriangulation. This function is based on the point insertion process described by De Floriani (1989). Having inserted p, some or all of the triangles in Flat_Region will have been corrected. In either situation, ResolveBetweenContours is exited; any ¯at triangles which remain are guaranteed to be encountered at a later stage, either during the main
FTC search or during the construction of a subsequent Flat_Region. Note that this edge-swapping approach is relatively ecient in limiting the number of new points inserted (see Table 1). Alternative methods, based solely on insertion of points, similar to those suggested by Vozenilek (1994) and Christensen (1987), could easily be implemented. Evaluation of which approach is in fact better will be the subject of future qualitative assessment experiments. Correcting within a single contour If a non-¯at triangle is not encountered, ResolveWithinContour is called. This function involves the addition of a single new point p at the centre of Flat_Region. The centre of Flat_Region is de®ned here as being the centre of its largest possible inscribed circle. As with ResolveBetweenContours, insertion of p will correct some or all of the triangles in Flat_Region; any ¯at triangles which remain are encountered at a later stage. It should be noted that FTC as described here does not dierentiate between closed contour situations which involve legitimate ¯at triangles (e.g. a lake) and those where ¯at triangles represent error (e.g. a truncated peak). In both situations, FTC
Table 1. Details of experimental data, CDTs prior to application of FTC and CDTs subsequent to application of FTC CDT prior to application CDT subsequent to application of FTC of FTC
Source data Data set 1 2 3 4
No. of contours 57 74 36 43
No. of points 4484 5562 2928 3384
No. of edges 4433 5494 2890 3337
Total no. of No. of ¯at Total no. of No. of points No. of ¯at triangles triangles triangles inserted triangles 8962 11118 5850 6766
1019 1059 1529 1048
9000 11162 5920 6970
19 32 35 14
0 0 0 0
146
J. M. Ware
Figure 7. Part of CDT generated from data set 2. (A) Before application of FTC. (B) After application of FTC.
Figure 5. Part of CDT generated from data set 1. (A) Before application of FTC. (B) After application of FTC.
Figure 6. Part of CDT generated from data set 1. (A) Before application of FTC. (B) After application of FTC.
Figure 8. Part of CDT generated from data set 3. (A) Before application of FTC. (B) After application of FTC.
A procedure for automatically correcting invalid ¯at triangles
147
Figure 9. Part of CDT generated from data set 3. (A) Before application of FTC. (B) After application of FTC.
Figure 11. Part of CDT generated from data set 4. (A) Before application of FTC. (B) After application of FTC.
Figure 10. Part of CDT generated from data set 3. (A) Before application of FTC. (B) After application of FTC.
Figure 12. Part of CDT generated from data set 4. (A) Before application of FTC. (B) After application of FTC.
148
J. M. Ware
attempts to correct the triangulation by calling ResolveWithinContour. In order to avoid an attempt at correcting a legitimate Flat_Region, it is necessary to have access to contour attribute data giving detail of what classes of geographical phenomena lie adjacent to individual contours.
EXPERIMENTAL RESULTS
The FTC procedure has been implemented in C on a Sparc workstation. Tests have been carried on CDTs generated from 1:10 000 scale UK Ordnance Survey digitised contour data. Table 1 gives details relating to source data, CDTs and experimental results. Note that in each experiment all ¯at triangles were successfully corrected. Parts of the triangulations, prior to (A) and subsequent to (B) application of FTC, are illustrated in Figures 5±12.
FINAL REMARKS
The paper outlines a procedure, FTC, for automatically correcting ¯at triangles occurring in triangulations generated from contour data. The procedure, having been tested on a number of data sets, has been shown to correct all ¯at triangles successfully in a given triangulation. The procedure, as it stands, provides users with a useful tool for improving models produced by existing CDT algorithms. Future work will concentrate on producing a qualitative assessment of the FTC approach. This paper has assumed that correcting ¯at triangles does indeed lead to an improved triangulated model; this assumption needs to be validated by, for example, a comprehensive visibility analysis study (see Kidner, 1996; Ware, Kidner and Herbert, 1996; for related work). Such a study will measure the eectiveness of Delaunay, constrainedDelaunay, non-Delaunay and FTC methods by comparing predicted results with ®eld observations, and analyse the alternative FTC edge-swapping and point-insertion strategies discussed here. A range of interpolation methods, including the cubic technique currently used, will be tested to estimate the elevation value for new points. AcknowledgmentsÐThe digital contour data used in experiments described in this paper has been provided by the UK Ordnance Survey.
REFERENCES Brandli, M. (1992) A triangulation-based method for geomorphological surface interpolation from contour lines. Proceedings 3rd European Conference and Exhibition on GIS (EGIS `92), Vol. 1, Munich, pp. 691±700. Christensen, A. H. J. (1987) Fitting a triangulation to contour lines. Proceedings 8th International Symposium of Computer Assisted Cartography, Baltimore, Maryland, pp. 57±67. Christiansen, H. N. and Sederberg, T. W. (1978) Conversion of complex contour line de®nitions into polygonal element mosaics. ACM Computer Graphics 12(3), 187±192. De Floriani, L. (1987) Surface representations based on triangular grids. Visual Computer 3(1), 27±50. De Floriani, L. (1989) A pyramidal data structure for triangle-based surface description. IEEE Computer Graphics and Applications, March 1989, pp. 67±78. De Floriani, L. and Puppo, E. (1988) Constrained Delaunay triangulation for multiresolution surface description. Proceedings 9th IEEE International Conference on Pattern Recognition, CS Press, Rome, 1988, pp. 566±569. Ekoule, A. B., Peyrin, F. C. and Odet, L. (1991) A triangulation algorithm from arbitrary shaped multiple planar contours. ACM Transactions on Graphics 10(2), 182±199. E.S.R.I. (1991) Surface modelling with TIN: surface analysis and display, ARC/INFO User's Guide. Environmental Systems Research Institute, Redlands, California. Fuchs, H., Kedem, Z. M. and Uselton, S. P. (1977) Optimal surface reconstruction from planar contours. Communications of the ACM 20(10), 693±702. Ganthapy, S. and Dennehy, T. G. (1982) A new general triangulation method for planar contours. ACM Computer Graphics 13(3), 311±319. Heller, M. (1990) Triangulation algorithms for adaptive terrain modelling. Proceedings 4th International Symposium on Spatial Data Handling, Vol. 1, pp. 163± 174. Kidner, D. B. (1996) Site selection and visibility analysis for a wind farm development: a problem for GIS? Proceedings 1st International Conference on GIS in Urban, Regional and Environmental Planning, Samos, Greece, pp. 220±237. Tipper, J. C. (1977) A method and Fortran program for the computerised reconstruction of three dimensional objects from serial sections. Computers & Geosciences 3(4), 579±599. Vozenilek, V. (1994) Generating surface models using elevations digitised from topographical maps. Proceedings 5th European Conference and Exhibition on GIS, Vol. 1, Utrecht, pp. 972±982. Ware, J. M., Kidner, D. B. and Herbert, M. J. (1996) Terrain in perspective: DEMs or TINs. Proceedings 2nd Joint European Conference on Geographic Information, Barcelona, pp. 388±397. Woo, T. C. (1985) A combinational analysis of boundary data structure schemata. IEEE Computer Graphics and Applications 5(3), 19±27.
A procedure for automatically correcting invalid ¯at triangles APPENDIX Pseudo-code for FTC struct triangulation {t:triangle_pointer; next:triangulation_pointer } struct triangle {vertex[1..3]:vertex_pointer; adjacent[1..3]:triangle_pointer; edge_¯ag[1..3]:Boolean } FTC(IN - T:triangulation, T_head:triangulation_pointer) current = T_head While current$NULL If (FlatTriangle(current.t)) ResolveFlatTriangle(T,T_head,current.t) Else current = current.next Endif Endwhile End FTC ResolveFlatTriangle(IN - T:triangulation, T_head:triangulation_pointer, t:triangle) FindFlatRegion(T, t, Flat_Region, head, tv) If (tv $ NULL) ResolveBetweenContours(T,T_head,tv) Else ResolveWithinContour(T, tv) Endif End ResolveFlatTriangle FindFlatRegion(IN OUT -
T:triangulation, t:triangle; Flat_Region:triangulation, FR_head:triangulation_pointer, tv:triangle)
tv = NULL FR_head = NULL TQ_head = NULL Add(t,Triangle_Queue,TQ_head) current = TQ_head While (current$NULL) ta = current.t Add(ta,Flat_Region,FR_head) For i =1 to 3 If (ta.edge_¯ag[i] =0) If (FlatTriangle(ta.adjacent[i])) If (ta.adjacent[i] ( Triangle_Queue) and (ta.adjacent[i] ( Flat_Region) Add(ta.adjacent[i],Triangle_Queue,TQ_head) Endif Else t v = ta Endif Endif Endfor Endwhile End FindFlatRegion
149
150
J. M. Ware
ResolveBetweenContours(IN - T, Flat_Region:triangulation, T_head,FR_head:triangulation_pointer, tv:triangle) Swapped = True While (Swapped) and (FR_head$NULL) ta = FindAdjacentTriangle(tv, Flat_Region) Swapped = SwapCommonEdge(tv,ta)) If (Swapped) Delete(ta,Flat_Region) t v = ta Else p = FindMidPoint(tv,ta) InterpolateElevation(p) AddPointToTriangulation(p,T,T_Head) Endif Endwhile End ResolveBetweenContours ResolveWithinContour(IN -
T, Flat_Region:triangulation, T_head,FR_head:triangulation_pointer) p = FindCentreOfRegion(Flat_Region, FR_head) InterpolateElevation(p) AddPointToTriangulation(p,T,T_Head) End ResolveWithinContour