TRIANGULATE

The TRIANGULATE procedure constructs a Delaunay triangulation of a planar set of points. Delaunay triangulations are very useful for the interpolation, analysis, and visual display of irregularly-gridded data. In most applications, after the irregularly gridded data points have been triangulated, the function TRIGRID is invoked to interpolate surface values to a regular grid. TRIANGULATE uses double precision for all computations.

Since Delaunay triangulations have the property that the circumcircle of any triangle in the triangulation contains no other vertices in its interior, interpolated values are only computed from nearby points.

TRIANGULATE can, optionally, return the adjacency list that describes, for each node, the adjacent nodes in the Delaunay triangulation. With this list, the Voronoi polygon (the polygon described by the set of points which are closer to that node than to any other node) can be computed for each node. This polygon contains the area influenced by its associated node. Tiling of the region in this manner is also called Dirichlet, Wigner-Seithz, or Thiessen tessellation.

The grid returned by the TRIGRID function can be input to various routines such as SURFACE, TV, and CONTOUR. See the description of TRIGRID for an example.

TRIANGULATE and TRIGRID can also be used to perform gridding and interpolation over the surface of a sphere. The interpolation is C₁ continuous, meaning that the result is continuous over both the function value and its first derivative. This feature is ideal for interpolating an irregularly-sampled dataset over part or all of the surface of the earth (or other (spherical) celestial bodies). Extrapolation outside the convex hull of sample points is also supported. To perform spherical gridding, you must include the FVALUE and SPHERE keywords described below.

Examples

Here we perform the triangulation on a dataset and display the results:

x = [96,171,107,153,150,51,194,92]

y = [183,185,253,306,232,267,272,395]

triangulate, x, y, triang

p = plot(x, y, 'or', sym_size=1.5, /sym_filled)

for i=0,triang.dim[1]-1 do $

  o = polyline(x[triang[*,i]], y[triang[*,i]], /data, conn=[4,0,1,2,0])

t = text(x, y, '  '+([0:x.length]).toString(), /data, clip=0, font_size=14)

For more examples using the TRIANGULATE routine, see the TRIGRID function.

Syntax

TRIANGULATE, X, Y, Triangles [, B] [, CONNECTIVITY=variable] [, SPHERE=variable [, /DEGREES]] [, FVALUE=variable] [, REPEATS=variable] [, TOLERANCE=value]

Arguments

X

An array that contains the X coordinates of the points to be triangulated.

Y

An array that contains the Y coordinates of the points to be triangulated. Parameters X and Y must have the same number of elements.

Triangles

A named variable that, on exit, contains the list of triangles in the Delaunay triangulation of the points specified by the X and Y arguments. Triangles is a longword array dimensioned (3, number of triangles), where Triangles[0, i], Triangles[1, i], and Triangles[2, i] contain the indices of the vertices of the i-th triangle (i.e., X[Tr[*, i]] and Y[Triangles[*, i]] are the X and Y coordinates of the vertices of the i-th triangle).

B

An optional, named variable that, upon return, contains a list of the indices of the boundary points in counterclockwise order.

Keywords

CONNECTIVITY

Set this keyword to a named variable in which the adjacency list for each of the N nodes (xy point) is returned. The list has the following form:

Each element i, 0 ≤i < N, contains the starting index of the connectivity list for node i within the list array. To obtain the adjacency list for node i, extract the list elements from LIST[i] to LIST[i+1]-1.

The adjacency list is ordered in the counter-clockwise direction. The first item on the list of boundary nodes is the subscript of the node itself. For interior nodes, the list contains the subscripts of the adjacent nodes in counter-clockwise order.

For example, the call:

TRIANGULATE, X, Y, CONNECTIVITY = LIST

returns the adjacency list in the variable LIST.

The subscripts of the nodes adjacent to X[i] and Y[i] are contained in the array:

LIST[LIST[i] : LIST[i+1]-1]

DEGREES

Set this keyword to indicate that the X and Y arguments contain longitude and latitude coordinates specified in degrees. This keyword is only effective if the SPHERE keyword is specified. If DEGREES is not set, X and Y are assumed to be specified in radians when a spherical triangulation is performed.

FVALUE

Set this keyword to a named variable that contains sample values for each longitude/latitude point in a spherical triangulation. On output, the elements of FVALUE are rearranged to correspond to the new ordering of X and Y (as described in the SPHERE keyword, below). This reordered array can be passed to TRIGRID to complete the interpolation.

REPEATS

Set this keyword to a named variable to return a (2, n) list of the indices of duplicated points. That is, for each i,

X[REPEATS[0,i]] = X[REPEATS[1,i]]

and

Y[REPEATS[0,i]] = Y[REPEATS[1,i]]

Note: Use the GRID_INPUT procedure to handle repeated points (duplicate locations).

SPHERE

Set this keyword to a named variable to indicate that spherical gridding should be performed. In this case, X and Y are interpreted as longitude and latitude in radians (or degrees if /DEGREES is set).

The results from the spherical triangulation are returned in the SPHERE variable as a structure that contains the following fields:

XYZ - A double array of length N containing the longitude/latitude values converted to xyz coordinates on the sphere.
IEND - An array of length N containing the index within IADJ of the end of each modified adjacency list.
IADJ - An array containing the indices of the neighbors in the adjacency sets. For example, for node k=0, the indices of the neighbors are stored in IADJ[0]...IADJ[IEND[0]]. For node k > 0, the neighbor indices are stored in IADJ[IEND[k–1]+1]...IADJ[IEND[k]]. Note that node k is a boundary node if and only if IADJ[IEND[k]] = 0.

This structure can then be passed to TRIGRID to perform spherical gridding.

Note: The X and Y parameters are converted to double precision and are rearranged to match the spherical triangulation.

TOLERANCE

Set this keyword to the tolerance to be used when determining whether points are colinear. The default is zero, which assumes that none of the points are colinear.

Tip: TRIANGULATE is typically used for irregular grids, where none of the points are colinear. For regular or semi-regular grids, you can use the TOLERANCE value to avoid creating artificial triangles along the boundary, and to triangulate the appropriate points within the grid. In this case, an appropriate tolerance value might be 10^-6 * Max(X,Y) (or 10^-12 for double precision), where Max(X,Y) is the maximum absolute value of the X and Y arrays.

Version History

Pre-4.0

Introduced

Resources and References

The spherical gridding technique used in IDL is based on the paper “Interpolation of Data on the Surface of a Sphere”, R. Renka, Oak Ridge National Laboratory Report ORNL/CSD-108, 1982.
The TRIANGULATE procedure uses the divide-and-conquer method described in “Two Algorithms for Constructing a Delaunay Triangulation,” D.T. Lee and B.J. Schachter, Int. J. of Computer and Information Sci., Vol. 9, 1980, pp. 219-242.
The storage mechanism for the SPHERE keyword is described in Cline, A.K.; Renka, R.L. "A storage-efficient method for construction of a Thiessen triangulation." Rocky Mountain J. Math. 14 (1984), no. 1, 119-140.