Tessellations and triangulations

<< Click to Display Table of Contents >>

Navigation:  Building Blocks of Spatial Analysis > Geometric and Related Operations >

Tessellations and triangulations

A region may be divided into a set of non-overlapping areas in many ways. As we noted in Table 1‑1 a (regular or irregular) tessellation of a plane involves the subdivision of the plane into polygonal tiles that completely cover it. Within GIS regular tiles are almost always either square or rectangular, and form a (continuous) grid structure. Regular triangular and hexagonal grids are also possible in the plane but are rarely implemented in software packages. GIS does make extensive use of irregular triangular tessellations, both for division of plane regions and as an efficient means of representing surfaces.

Delaunay Triangulation

Given a set of points in the plane, P, lines may be drawn between these points to create a complete set of non-overlapping triangles (a triangulation) with the outer boundary being the convex hull of the point set. It has long been known that triangulation provides a secure method for locating points on the Earth’s surface by field survey. It has also been found that a desirable characteristic of such triangulations is that long thin triangles with very acute internal angles are to be avoided in order to provide the best framework for measurement and analysis.

To ensure the triangulation will have the best chance of meeting the characteristics desired, a construction rule was devised by the mathematician Delaunay. This states that three points form a Delaunay Triangulation if and only if (iff) a circle which passes through all three points contains no other points in the set. Figure 4‑32, below, illustrates this process for a set of spot height locations in a small area south of Edinburgh in Scotland (GB Ordnance Survey tile NT04). Each location is represented by an (x,y) pair. Many programs support the creation of the triangulation required, for example the delaunay(x,y) function in MATLab where x, and y are nx1 vectors containing the coordinates of all the points, or (as here) the Grid Data operation in Surfer, using the option to export the triangulation used as a base map. The circles illustrate how the circle circumscribing each triangle has no other points from the set within it. Worboys and Duckham (2004, pp 202-207) provide a useful discussion of triangulation algorithms.

A number of other observations should be made about this process:

this triangulation is unique in the sense that it is the only triangulation of the point set that satisfies the rule described above. However, with some point set arrangements (e.g. four points at the corners of a square or rectangle) multiple Delaunay triangulations are possible, hence to this extent it is not unique

long thin triangles with acute angles do occur using this procedure, although the number and severity of these is minimized compared to other triangulations

the center of the circumscribing circle (marked with a here) may lie outside the triangle in question, as is the case with one of the circles illustrated

this procedure is by no means the only, or necessarily the best method of triangulation for all purposes. For example, when modeling a surface there may be important breaks of surface continuity or surveyed transects which it is highly desirable to include within the dataset, which are not captured by simple Delaunay-type procedures (as noted in Section 4.2.13, Boundaries and zone membership)

If all points included in the Delaunay triangulation of a set P also lie on the convex hull of P the result is known as a furthest point Delaunay triangulation. The set of lines (edges) that belongs to all possible furthest site Delaunay triangulations is known as a Delaunay diagram, and has important applications in network analysis (see further, Chapter 7, Network and Location Analysis).

Figure 4‑32 Delaunay triangulation of spot height locations


TINs — Triangulated irregular networks

Regular lattices such as grids provide a very inefficient means of storing details of surface topography. Areas that have similar slope could be better represented and stored using a series of inter-connected triangles forming an irregular mesh or network (a TIN). The nodes or vertices of these triangles can be placed closer together where landscape detail is more complex and further apart where surfaces are simpler. TINs may be designed or derived. Designed TINs arise in engineering and surveying models, whereas derived TINs are typically programmatically generated from grid-based datasets (digital elevation models or DEMs).

Although such TINs could be generated by selecting vertices using key feature extraction and surface modeling, a simpler and faster approach is to systematically divide the entire region into more and more triangles based on simple subdivision rules. For example, if you start with a rectangular region and define vertices at each of the four corners, two triangulations are possible: (i) a line connecting the NE-SW corners creates one pair of triangles; or (ii) a line connecting the NW-SE corners creates a second pair. For each triangle a statistic can be computed, such as the maximum deviation in elevation of the triangle from the original DEM. The arrangement that produces the smallest maximum absolute deviation (a minimax criterion) can then be selected as the preferred triangulation. Each of the two initial triangles can then be subdivided further. For example, a point can be placed at the centroid of each triangle which creates a further three triangles, or 6 in total. The minimax criterion test can be applied for each of the subdivisions of the two original triangles and the process proceeds iteratively until the minimax criterion is less than a pre-defined threshold value (e.g. 5 meters) or the number of triangles generated has reached a maximum value, also pre-defined. This kind of procedure is used in many GIS packages — for example TNTMips describes this process as adaptive densification. TINs and DEMs are discussed in more detail in Chapter 6 (see for example, Figure 6‑5).

Voronoi/Thiessen polygons

Given a set of points in the plane, there exists an associated set of regions surrounding these points, such that all locations within any given region are closer to one of the points than to any other point. These regions may be regarded as the dual of the point set, and are known as proximity polygons, Voronoi polygons or Thiessen regions. These various terms are equivalent, as is the term Dirichlet cell, all being derived from the names of their proponents: Dirichlet (in 1850); Voronoi (in 1909); and Thiessen (in 1912). Detailed explanations and discussion of Voronoi polygons is provided in Boots (CATMOG 45) and Okabe et al. (2000).

Figure 4‑33A shows an example of Voronoi region creation, in this case using ArcGIS. The point set from which the map was generated is also shown and as can be seen, points are frequently close to polygon edges. Indeed, in some applications it may be desirable to exclude points that are closer together than some predefined tolerance level. Note that in this example Voronoi regions at the edges of the sample extend to infinity, and are shown as bounded (e.g. by the MBR of the point set plus 10%). Figure 4‑33B shows the same point set plus the associated Delaunay triangulation once more, created within the MATLab mathematical modeling package, which (correctly) limits the set of regions generated to those that can be unambiguously computed. Figure 4‑33B enables us to see exactly how Voronoi regions are created. There are two steps: the first is the triangulation process, connecting the point set as described earlier; and the second is bisection of the triangulation lines of adjacent points to define Voronoi polygon boundaries.

There are many reasons an analyst may wish to generate Voronoi polygons from a point set, but a common requirement is to assign values to regions or cells based on the attributes of the point set. The most common assignment is that of the individual point, or nearest-neighbor model. Other options are offered by some GIS packages in much the same way as is applied to grid and image files (in which each cell is surrounded by 4 or 8 immediate neighbors, depending on the model used). Note that the MATLab implementation does not extend to the MBR or other bounding region (e.g. rectangle) defined by the user. A selection of (a) MATLab and (b) GRASS functions that support operations discussed in this section is shown in Table 4‑7.

Table 4‑7 Selected MATLab/GRASS planar geometric analysis functions



(a) convhull; (b) v.hull

Compute the convex hull of a point set, optionally providing its area

(a) delaunay; (b) v.delaunay

Delaunay triangulation

(a) dsearch

Nearest point search of Delaunay triangulation. Identifies the vertex of a Delaunay triangle that is closest to a specified point

(a) inpolygon

True (1) for points inside polygonal region, False (0) for points outside, and assigned 0.5 if points lie on the boundary

(a) polyarea; (b) area_poly2.c

Area of a polygon

(a) rectint

Area of intersection for two or more rectangles

(a) tsearch

Closest triangle search. Identifies the Delaunay triangle that includes a specified point

(a) voronoi; (b) v.voronoi

Voronoi diagram creation

Figure 4‑33 Voronoi regions generated in ArcGIS and MATLab





Assignment facilities provided within ArcGIS include: local mean value — the cell is assigned the average value of the cell and its direct neighbors; and local modal value — determined by computing the frequency distribution of all adjacent cells in 5 classes and then assigning the cell the class value from this set that is the most common amongst the cell and its immediate neighbors. These assignments represent local smoothing procedures. Local variation rather than values may also be computed. For example:

local entropy value (equivalent to Shannon’s Landscape Diversity Index, SHDI — see Section 5.3.4, Landscape Metrics) — in this model the adjacent cells are again grouped into 5 sets (smart quantiles or natural breaks classes) and the number in each class is represented as a proportion of the total, pi, i=1‑5, from which Shannon’s entropy statistic is computed (Table 1‑3) and this value is assigned to the cell — values will range from 0 to log2(5); and

local standard deviation and local inter-quartile range — these statistics are computed for the cell and its immediate neighbors and the value assigned to the selected cell

Voronoi polygons may also be generated in grid or raster datasets (perhaps better described as Voronoi patches). Since the definition of a Voronoi polygon in the plane is that it consists of all locations closer to selected points than to any other points, such regions can be generated by a form of spreading algorithm. In this case we start with a number of points (the source points, sometimes called the target points) and spread outwards in all directions at a constant rate, like ripples spreading across the surface of a pond when a stone has been thrown into the middle. The process stops when the ripple (spreading) from another point is encountered. The rectangular grid structure and the number of search directions (typically 8) affect the form of the resulting cells (Figure 4‑34).

With the raster representation Voronoi patch shape and boundary definition are not as clear cut as in vector operations. Within GIS packages processing activities of this type are usually carried out using techniques described as cost distance, or accumulated cost distance tools (e.g. in ArcGIS Spatial Analyst). In image processing the procedures adopt a somewhat different set of algorithms, based on scanning rather than spreading, which are known as distance transforms. The latter may be exact (in Euclidean distance terms) or approximate (see further, Section 4.4.3, Distance Transforms). Furthermore, with both approaches the space across which the scan or spread procedure is conducted may be homogenous (i.e. treated as uniform in all directions and at all points) or non-homogenous, containing areas of varying difficulty to traverse, varying cost, and/or including obstacles.

Figure 4‑34 Voronoi cells for a homogeneous grid using a 3x3 distance transform


Okabe and colleagues have developed a software package, SANET, which can be installed as an ArcGIS add-in for “Spatial Analysis on a Network”. The SANET software includes a network-based Voronoi tool, very similar to some of the network analysis tools provided in ArcGIS itself and in network-oriented packages such as TransCAD and Cube. An illustration of its facilities is shown in Figure 4‑35. In this example, showing the Shibuya district of Tokyo, the point set (off-network sites which are Bank premises) have been linked to the nearest network edge, and then non-directed network Voronoi regions generated from these locations (new network nodes) using shortest path distances. If one-way streets are taken into account separate inward- and outward-directed Voronoi regions are computed (see Okabe and Satoh, 2008, for more details). In each case these regions derived using unweighted calculations. Weighted calculations, for example taking into account the cost or time of travel across network links results in additional region definitions.

Figure 4‑35 Network-based Voronoi regions — Shibuya district, Tokyo


Data courtesy: Profs A Okabe, K Okunuki and S Shiode (University of Tokyo, Japan)