Multibeam data processing: Adding and deleting vertices in a Delaunay triangulation

Gert Brouns 
Department of Geography, Ghent University, Belgium
Prof Alain De Wulf
Department of Geography, Ghent University, Belgium
Dr Denis Constales
Department of Mathematics, Ghent University, Belgium

ABSTRACT

Multibeam echosounder measurements serve to make a digital terrain model of the seafloor. The Delaunay triangulation is a widely appreciated and investigated mathematical model to represent the relief of such a terrain. However, the need for editing functionalities for the triangulated model rises.

A method for deleting a single vertex from a Delaunay triangulation is proposed. It is mathematically proved to be sufficient to calculate the conforming Delaunay triangulation of the polygonal hole that is created after eliminating a vertex to obtain the Delaunay triangulation of the reduced data set. The method has been extended for deleting a set of vertices.

A new method to replace part of an existing triangulation by new data is proposed; the method is adapted from the merge step in the divide-and-conquer algorithm.

Introduction

Obtaining an accurate model of the seafloor is a major concern in dredging works. Recent hydrographic surveying tools, especially the multibeam echosounder, yield a very dense sampling of the seafloor and this immense amount of data needs to be processed, according to some constraints imposed by the client, to form an accurate terrain model.

Two main categories of digital terrain models (DTM) exist: grid models and triangulated irregular networks (TIN).

Grid models represent the surface with a regular grid, each cell of which has a height value attributed to it. This happens by interpolating the measured values according to specific algorithms that range from local linear interpolation, over weighting inversely proportional to the square of the distance between grid cell centre and a number of neighbouring points, to polynomial interpolation. Triangulated irregular networks do not use interpolation; they utilise the discrete measurements as data points. This paper briefly depicts the main advantages and disadvantages of both methods.

Advantages and disadvantages of the different Digital Terrain Models

Interpolation

As outlined in the introduction, grid representation models always use some kind of interpolation to attribute values to the grid cells. By definition, interpolation does not exploit the original, unevenly spaced measurements to represent the terrain. Either values are attributed to grid cells in which no measurement has been taken, or if several measurements fall into one grid cell, they are clustered to result in only one grid cell value. Also inversely, if extracting the height of a measured spot from the grid model, the original measurement value is not reproduced. This means there is a double loss in precision and it is very hard to exert control on the link between the measurements and the derived grid model and therefore on the quality of the latter.

Smoothing

Interpolation results in a smoothing of the data. Sharp features will be rounded off and a certain level of detail is always lost. Polynomial interpolation, where a polynomial function is fitted through the data set can, in theory exactly describe the surface if enough polynomial terms are calculated. ‘Enough’ means a polynomial term, a point. Once the number of points n is getting over a thousand, it will become increasingly difficult to calculate these terms as systems involving matrices of dimensions n x n need to be solved. For multibeam data sets, which typically contain several hundreds of thousands to a few million points, this will clearly become quite inefficient. A solution is to calculate less terms and use a least squares fitting, but this inevitably will cause smoothing again. Smoothing has two important consequences in the specific field of marine engineering: contractors such as dredging companies always create additional relief on the seafloor. If, due to the measurement processing techniques, this artificial relief with often sharp discontinuities is smoothed, it will always be to the contractor’s disadvantage. On one hand, upper and lower limits of artificial slopes will be vaguely defined in both construction and as built charts, on the other hand volumes will systematically be underestimated in both dredging and land reclamation activities.

Relevance of data

Grid models usually employ rectangular representation areas. This generally does not reflect the real shape of the surveyed area, so extrapolation of the measurements or a dummy value (zero height) is sometimes used to fill up the gaps in the grid. In triangulated models, it is visually easier to perceive the quality and relevance of the information near the borders as firstly the whole area will be bordered by a convex hull, and secondly some triangles near the border will easily be interpreted as ‘filling triangles’ due to their huge area or elongated shape, from which it is readily (as well visually as mathematically) deduced that they contain no relevant information. In a triangulated model, the relevance of each point is clear: each point represents a measured value. In the grid model, the relationship between measurements and grid values is totally absent and interpretation becomes increasingly difficult.

Non-homogeneously spaced data

Due to the nature of triangulated irregular networks, they will easily represent non-homogeneously sampled areas. If an area contains important features to be displayed, a denser sampling may be carried out. Around this core area, coarser sampling will suffice to show the general environment in which the survey takes place. A triangulated model, using the discrete measurements to represent the relief, will easily deal with this kind of terrain. Grid models would either have to fill the entire grid, taking lots of superfluous memory to store the data, or would require an extremely complex data structure to represent a variable density grid. The latter would eliminate the main advantage of grid models over triangulated models, which is their simplicity of implementation and subsequent manipulation.

Delaunay triangulation

In summary, it can be stated that thanks to the absence of interpolation, the consequent absence of smoothing, the clear relevance of data in the model and the natural representation of non-homogeneously sampled areas and in spite of its more difficult construction and ulterior manipulation, a triangulated irregular network is preferred to a regular grid model in situations where the highest level of accuracy of the digital terrain model is required.

For additional quality reasons, the kind of triangulated irregular network is chosen to be a Delaunay triangulation; meeting the Delaunay requirement (see below) means the triangulated model is entirely reproducible, which is a quality-control asset for the client.

Editing functionalities in grids and triangulated models

The initially obtained model very often needs to be updated after construction, when new information becomes available. In grid models, implementation of updating algorithms is fairly easy thanks to the evenly spaced grid cells, unlike in triangulated networks. Little basic work in that domain has been done for triangulated networks so far.

Updating may either mean adding points to the diagram or, on the contrary, deleting vertices from it. In either process, different routes can be followed; it is not only possible to insert or delete vertices one by one, it may also be necessary and quicker to delete or add a whole set of vertices at once or replace an old set of vertices by a new one.

After an initial recapitulation of some definitions, this paper addresses treats deletion of vertices from an existing Delaunay triangulation. Results thereof will be used in the section that handles the problem of inserting and replacing vertices.

Introduction to Delaunay triangulations

The definitions of a triangulation, of a Delaunay triangulation and those of some specific types of Delaunay triangulation are given below.

Definitions [1], [2]

A general triangulation of a set of vertices V is a graph that connects the vertices of V by open edges that are mutually empty, that is, they do not intersect each other and they meet in their open endpoints at the vertices of V.  The result is a subdivision of the plane in a number of disjoint, open triangles, called faces of the triangulation. Faces, edges and vertices together fill the convex polygon that is enclosed by the convex hull, the latter consisting of edges and points on the boundary of the polygon. The triangulation and the outer region of the polygon fill the entire plane.

Fig. 1: A general triangulation of a set of vertices

A triangulation is called a Delaunay triangulation when no vertex lies within the circumscribing circle of any triangle of the triangulation.

Fig. 2: Circumscribing circles

Fig. 3: The Delaunay triangulation of the given set of vertices

A constrained Delaunay triangulation is a triangulation of a set of points and straight-line segments that tries to respect the Delaunay property, but as the segments are preserved, the Delaunay property does not necessarily hold in all triangles near the segments. It is usually formed by creating a Delaunay triangulation of the set of vertices, subsequently inserting the straight-line segments and locally modifying the triangulation to eliminate edges that cross this segment.

Fig. 4: A constrained Delaunay triangulation

A conforming Delaunay triangulation is a constrained Delaunay triangulation where the Delaunay property does hold in all triangles. It is generally generated from a constrained Delaunay triangulation by adding points to the triangulation (called Steiner points) on the segments and hence by creating new triangles.

Fig. 5: A conforming Delaunay triangulation

Deleting vertices from a Delaunay triangulation

Deletion of a single vertex

What results after deletion of a vertex is not necessarily a convex polygon (see the counter example in Figure 6), but a star shaped polygon. A star shaped polygon is by definition any polygon of which all points in the interior are visible from at least one point in its interior[1], so it can be concave as well. It is trivial to show that the resulting polygon after deletion of a vertex is star shaped as all edges and vertices of the polygon can be seen from the position of the deleted vertex, so at least one point exists that fulfils the requirement in the definition.

Fig 6a: Vertex to be deleted                     Fig 6b: Influenced area                     Fig 6c: New triangulation

Fig. 6: Deletion of a vertex creates a star shaped polygonal hole that needs to be retriangulated.

Two lemmas and a corollary are stated and proved in order to show that it is sufficient to construct the conforming Delaunay triangulation of the hole that is created when removing a vertex, to obtain the triangulation of the reduced data set.

Definitions

Let v be a vertex of the set of vertices V that does not belong to the convex hull of the Delaunay triangulationD(V) of V. Let D(V’) be the Delaunay triangulation of the set of vertices V’ = V/{v}.

Let R be the set of edges z on the boundary of the star shaped polygon that results after deleting v from D(V).

Let P be the set of the endpoints of the edges in R.

The following figures illustrate the above definitions:

Fig. 7: Delaunay triangulation of set V

Fig. 8: Delaunay triangulation of set V’

Fig. 9: Set of edges of the star shaped polygonal hole

Lemma 1

When a vertex v is deleted from the Delaunay triangulation of a set of vertices V, only those triangles containing v will become invalid in the triangulation of the set V/{v}.

Proving that at least those triangles having v as a vertex become invalid is trivial; when a vertex of a triangle is taken away, the triangle has necessarily to be removed (see Figure 10). This proves that at most those triangles become invalid, which is equivalent to proving that any triangle in the Delaunay triangulation that does not contain that vertex remains a valid triangle.

Fig. 10: Illustration Lemma 1

Proof

From the definition of a Delaunay triangulation, it follows that all circumscribing circles of the existing triangles are empty of vertices. When deleting a vertex no new vertices are added to the triangulation, so all existing circles remain empty of vertices. Hence, for all existing triangles, the Delaunay condition remains true.

Lemma 2

All triangles that are in the Delaunay triangulation D(V’) of V’=V/{v} that are not in D(V) are in the Delaunay triangulation D(P).

Fig. 11: Illustration Lemma 2

The shaded triangles in D(V’) (Figure 11) do not occur in D(V) (Figure 10). As can be seen in  Figure 11, all these triangles are present in D(P). It is now possible to proceed to a formal proof of Lemma 2 that is illustrated with data set V of Figure 7.

Proof

Fig. 12: Proof Lemma 2 – part 1

Each edge z of R is an edge that belongs either to a triangle of V that does not contain v (edges 1, 2 and 3 in Figure 12) or to the convex hull of V (edges 4 and 5 in Figure 12). In the first case, from lemma 1, z is also an edge of a triangle of V’=V/{v}. In the second case, z is on the convex hull of V/{v} as the convex hulls of V and V/{v} are identical, which is true because v is not on the convex hull of V. So all edges of R belong to D(V’).

Fig. 13: Proof Lemma 2 – part 2

If t is a triangle of D(V’) that does not exist in D(V) (triangles T1, T2 and T3 in Figure 13), then it lies in the convex polygon enclosed by the edges in R, and as it is a Delaunay triangle, no point of V’=V/{v} will be in its circumscribing circle (Figure 13). This means also that no point of PÌV falls inside that circle and hence, t will be a triangle of D(P).

Corollary

If v is not on the convex hull of the Delaunay triangulation D(V) then the Delaunay triangulation D(V’) of V’=V/{v} falls apart in the conforming Delaunay triangulation of the convex hull of V minus that of R, and the conforming Delaunay triangulation of R.

Fig. 14: Illustration corollary

Proof

The proof is immediate from lemmas 1 and 2, as one can omit the triangles in the Delaunay triangulation D(P) of P that are not contained inside the polygon defined by the edges in R, if any.

The Corollary implies that the star shaped polygonal hole created in D(V) can simply be filled with the triangles of its conforming Delaunay triangulation to obtain the Delaunay triangulation of the new data set V/{v}.

Construction of the triangles in the polygonal hole

Two different approaches exist to obtain the Delaunay triangulation of the polygonal hole that is created when deleting a vertex in a Delaunay triangulation: triangulating the polygon as a geometric primitive or constructing the conforming Delaunay triangulation from the vertex set containing the vertices of the polygon.

A number of triangulation algorithms have been devised to construct the (general) triangulation of a simple polygon (a simple polygon’s boundary’s constituting edges do, by definition[1(pp.3&11)], not intersect. As this condition is less restrictive than for star shaped polygons, the latter are a subset of the former and properties of simple polygons are also valid for star shaped polygons)[3]. The resulting triangulation can be locally transformed to a Delaunay triangulation by swapping edges. As the average number of triangles to which a vertex belongs is six or less[1(p.161)], the average number of vertices on the boundary of the polygonal hole is six also, so time complexity of the algorithm, which is a measure for its performance in terms of speed, is not a real concern here. An algorithm the time complexity of which is quadratically dependent on the size of the vertex set will be acceptable.

Instead of treating the star shaped polygon as a geometric primitive and triangulating it as such, probably the neatest solution is to simply construct the constrained Delaunay triangulation of that polygon, as suggested in corollary. In contrast to the first method, this is done by treating the points of the polygon as a new vertex set that needs to be triangulated. The edges of the conformed Delaunay triangulation will subsequently be constrained to the edges of the (star shaped) polygonal hole.

Special case: deleting a vertex on the convex hull

The conforming Delaunay triangulation will also deal with the deletion of a vertex on the convex hull of the initial triangulation. This is illustrated in Figure 15. In such a case, an open polygon is created, a case that was explicitly not treated in the preceding lemmas and corollary.

15a:                                              15b:                                                 15c:
           Original triangulation                     Half open polygondefined by R                    Updated Delaunay triangulation

                       Fig. 15: Deleting a vertex on the convex hull

Generally, the open side of the polygon defined by R will be closed by joining the beginning and the end of the polyline R. This line will be a new edge on the convex hull of D(V’), connecting the two vertices where the convex hull of D(V) was cut off due to the elimination of more than one triangle bordering the convex hull. Sometimes, however, the extremities of the open polygon will not be ‘visible’ to each other (as is the case in Figure 15); a line connecting both extremities creates a non-simple polygon if the line intersects one or more of the edges of R. In such a case, the polygon is closed by connecting two vertices on R which are as far away from each other as possible, resulting in a simple polygon and one or more dangling edges that will not be a part of it. Once these dangling edges are identified and the polygon created, filling the latter with triangles can be done as proved above. The dangling edges will, together with the closing line, complete the new convex hull of D(V’) (Figure 15c).

Deletion of a set of vertices

If a set of vertices needs to be deleted, it is possible to proceed slightly differently than a one-by-one approach: all vertices and corresponding triangles can be deleted at once and only then the resulting hole or holes have to be triangulated. There is a major concern, however: the resulting polygon is not necessarily star-shaped. It could well be that deleted areas are ‘hidden’ for each other by subsisting triangles of the original Delaunay triangulation.

Again, a conforming Delaunay triangulation will bring an elegant solution here. This means constructing the Delaunay triangulation of the vertices on the boundary of the simple polygon with its boundary edges as mandatory edges, eliminating all triangles that fall outside the simple polygon and inserting this triangulated polygon into the polygonal hole of the original triangulation.

Inserting vertices into an existing Delaunay triangulation

Inserting vertices one by one

Inserting individual vertices into an existing Delaunay triangulation is the basic operation of the incremental algorithm[4], [5].

 Fig. 16: Inserting a vertex in a Delaunay triangulation

The update step can be described as follows: the triangle containing the new vertex is identified, all neighbouring triangles are checked whether or not their circumscribing circle contains the new vertex, if so they are marked and subsequently eliminated (as they do not meet the Delaunay requirement any more)  and finally new Delaunay triangles are constructed.

The time complexity for adding one vertex to an existing triangulation is linearly dependent on n, the number of vertices in the triangulation. In the worst case, all triangles must be checked and the number of operations that need to be performed is linearly related to the size of the triangulation. Usually, of course, performance is much better as only a limited number of triangles in the neighbourhood of the inserted vertex will be affected.

Inserting a block of vertices – practical applications

A novel method to replace a block of vertices in an existing triangulation with new data is presented.

Replacing old data by new data is of particular interest in hydrography. Assume a detailed bathymetric chart of a port is available and needs to be updated on a regular basis. Therefore, a survey vessel, equipped with a multibeam echosounder, is going out regularly. A large area of the port cannot be densely sampled in a one-day campaign but not the entire chart area can be updated yet, as this requires much more than one day’s work. Nevertheless, this recent information is invaluable for port authorities to know, for instance, fairway depths, so the existing chart should be updated as soon as possible.

Another closely related area of application is port maintenance and dredging. Dredged areas are resurveyed and the old data for these areas is replaced in the existing depth chart. In that way, the chart is kept as up to date as possible, allowing for an active intervention in the planning of the dredging activities.

In both cases, the hydrographer wants to fit his newly collected data into the existing chart. Obviously the depth information gathered in that day’s campaign will prevail over the old depth information. Areas which have been left unsounded will keep their old depth values. Removing the old data and seamlessly inserting the new one in an optimal and rigorous manner will be the challenge here.

Inserting a block of vertices – theoretical elaboration

Inserting a block of vertices as a whole is only possible if this block is to replace a part of the existing triangulation. If old data needs to be kept in place and the new vertices are only additional information and not replacing information, the incremental algorithm as described above has to be employed.

The method now proposed is based on the merge step of the divide-and-conquer algorithm as presented in Shewchuck[6], see also [7] and [8]. This divide-and-conquer algorithm, intended to construct a Delaunay triangulation of a set of vertices V in a very fast way, repeatedly cuts V in half until a base case is attained and then recursively stitches together the intermediate triangulations (Figure 17). Since this algorithm is an adaptation to fulfil specific needs, a few problems that do not occur in the original algorithm have to be overcome.

1.   Points of the two sets occur in the same spatial area. As the new data set prevails over the old one, the old data needs to be partially eliminated. A hole will have to be cut into the old data set in which the new data set fits.

2.   Whereas in the original divide-and-conquer algorithm, the data sets are separated by a vertical or horizontal line, the separating barrier is now curvilinear.

3.   The start edge for creating new triangles in the old algorithm is the lower common tangent connecting both convex hulls. Here another starting criterion needs to be found.

                                           Fig 17a: Original data set (left) and new data set (right) Fig 17b: Affected  triangles

                                   Fig 17c: Curvilinear separating barrier Fig 17d: Merge step     Fig 17e: New triangulation

Fig. 17: Inserting a new Delaunay triangulation in an old one.

The major difficulty is situated in the first stage, when creating a hole in the old triangulation.

If the spatial resolution of the new triangulation is coarser than that of the old one, instead of creating a hole in the old triangulation, a set of small holes will be created around the places where a new point is located, the neighbouring new point being too far to create a link between the two open spaces. This is illustrated in Figure 18.

Fig. 18: Detecting affected triangles when replacing old data (left) by new (right) Creation of ‘holes’ in the old triangulation

Taking the convex hull as a cutting line does not seem a good option, as the new area might be non-convex. A better idea is to eliminate all sliver triangles and border triangles. Both kinds of triangles are represented in Figure 19.

Fig. 19: Sliver and border triangles

Sliver triangles are the long, elongated triangles that appear near the convex hull of most Delaunay triangulations. They can be detected by their small minimum angle. The numerical value of that angle can be arbitrarily fixed depending on the needs, for instance to 25°.

Border triangles also deviate from a mean value; they have a comparably large edge length to the mean and do occur when the sampled area is, for instance, banana-shaped as in Figure 19.

Experience has shown that using a maximum area value to eliminate border triangles does not work very well. It is a much worse criterion than the proposed maximum edge length condition.

What results after eliminating sliver and border triangles is the concave hull. This hull is mathematically much less well defined than the convex hull, as it depends entirely upon the criteria used to eliminate triangles from the triangulation.

In Figures 20 to 22 an example is given of the difference between using the convex hull and the concave hull as cutting line.

Fig. 20: Replacing old by new data: detecting affected triangles

Fig. 21: Cut and merged triangulations using the concave cutting criterion

Fig. 22: Cut and merged triangulations using the convex cutting criterion

This concave cutting line is a good alternative to merely detecting affected triangles when the spatial resolution of the new data set is coarser, as is the case in Figure 22. However, no real-world tests have been carried out on that method yet.

Practical applications and benefits of THE TIN

Triangulated irregular networks (TIN) are a favourite scheme to construct a digital terrain model (DTM) from a seafloor measured at discrete spots, especially in marine engineering, where the highest accuracy of the digital terrain model is required. The alternative method, overlaying the sampled points with a regular grid, has four important drawbacks, as treated extensively earlier:

•    Grids use interpolation of the measured data, causing loss of information and loss of precision of the resulting terrain model.

•    Grids result always in a smoother representation of the relief than the reality. Especially for volume calculations, this is disadvantageous.

•    Because of interpolation, the relevance of the grid cell values are unclear unlike in triangulated models where each value represents a measurement.

•    The grid model is not adaptive. Whereas TIN will naturally represent areas with detailed relief information with a denser triangle pattern than more coarsely sampled areas, grids will be much more inflexible to deal with varied levels of detail and would lose their main advantage, which is their ease of calculation and manipulation, if one tried to do so.

Utilising a Delaunay triangulation assures maximum quality for the resulting DTM and its derived products such as volume calculations. This is for instance important for both dredging contractor and client. The contractor has a powerful tool at his disposal as the theoretically backed proposed editing tools allow him to update and adapt the work charts with the most recent data flexibly, without loss of accuracy which would occur when using a grid or rough cut-and-paste techniques. The client has the assurance that rigorous and precise techniques are used to obtain the volume calculations.

Implications of the proposed techniques on multibeam data sets

With the proposed techniques, it will be possible to edit triangulated irregular networks resulting from multibeam echosounding data sets nearly as flexibly as it is to edit grids. Both manual editing (point-by-point insertion and elimination of vertices) and automated techniques (replacing areas with old data by new) are presented. Hence the multibeam data, in spite of their large volume, can finally be used for what they’re worth, with minimal loss of accuracy which is inherent when using grids.

Conclusions

A method for deleting a single vertex from a Delaunay triangulation has been proposed and it has been mathematically proved to be sufficient to calculate the conforming Delaunay triangulation of the polygonal hole created after eliminating a vertex to obtain the Delaunay triangulation of the reduced data set. The method has been extended for deleting a set of vertices.

The incremental Delaunay triangulation algorithm has been indicated to be the best method by which to insert new vertices one by one in an existing Delaunay triangulation. However, this method is not adapted to replace part of an existing triangulation with a new triangulated data set. To do that a new method, that is a generalisation of the merge step in the divide-and-conquer algorithm, has been proposed.

The proposed editing tools allow full exploitation of the gathered multibeam data. Whereas multibeam data is usually interpolated by assigning values to a regular grid, with inherent loss of accuracy, it is now possible to deal directly with the discrete measurements.

References:

1    O’Rourke, J., 1998. Computational Geometry in C, Cambridge University Press, 2nd edition.

2    Fleischmann, Peter, 2000. 5 Delaunay triangulation. www.iue.tuwien.ac.at/diss/fleischmann/diss_new/node41.html.

3    Kumar, S., Surface Triangulation: A Survey, University of North Carolina, 1996, 
      www.citeseer.nj.nec.com/kumar96surface.html

4    Su, P and Drysdale, R.L.S., 1996. A Comparison of Sequential Delaunay Triangulation Algorithms,
      www.cs.berkeley.edu/~jrs/mesh/present.html

5    Brouns, G., de Wulf, A., Constales, D. Delaunay triangulation algorithms useful for multibeam echosounding. 
      Submitted for publication.

6    Shewchuck, J.R., 1996. Triangulation Algorithms and Data Structures,
      www.cs.cmu.edu/~quake/tripaper/triangle2.html

7    L. Guibas and J. Stolfi, 1985. Primitives for the manipulation of general subdivisions and the computation 
      of Voronoi diagrams, ACM Transactions on Graphics, April 1985, 4(2):74-123

8    Rex A. Dwyer, 1987. A Faster Divide-and-Conquer Algorithm for Constructing Delaunay Triangulations, 
      Algorithmica 2(2):137-151, 1987

Gert Brouns

is a researcher at the Department of Geography at Ghent University. He previously held a two year fellowship at CERN in Geneva, Switzerland. He obtained an M.Sc. in Surveying at University College London in 1998 and graduated in 1997 with a degree in Geography – Land Surveying option, from Ghent University.

Alain De Wulf

has been Professor in Surveying and Geomatic Engineering in the Department of Geography at Ghent University since 1995. He obtained a Ph.D. in Engineering in 1993, and previously obtained an M.Sc. in Civil Engineering and one in Industrial Management. He also graduated with a degree in Computer Science from the same university.

Denis Constales

has been a Postdoctoral Researcher in the Department of Mathematical Analysis at Ghent University since 1995. He obtained a Ph.D. in Mathematics in 1989 and previously graduated with a degree in Mathematics from the same university.