QHULL

The QHULL procedure constructs convex hulls, Delaunay triangulations, and Voronoi diagrams of a set of points of 2-dimensions or higher. It uses and is based on the program QHULL, which is described in Barber, Dobkin and Huhdanpaa, "The Quickhull Algorithm for Convex Hulls," ACM Transactions on Mathematical Software, Vol. 22, No 4, December 1996, Pages 469-483.

Syntax

Keywords: [, BOUNDS=variable ] [, CONNECTIVITY=variable ] [, /DELAUNAY ] [, SPHERE=variable ] [, VDIAGRAM=array ] [, VNORMALS=array ] [, VVERTICES=array ]

Arguments

An input argument providing an nd-by-np array containing the locations of np points, in nd dimensions. The number of dimensions, nd, must be greater than or equal to 2.

V0, V1, V2, ..., V(N-1)

Input vectors of dimension np-by-1 elements each containing the i-th coordinate of np points in nd dimensions. A maximum of seven input vectors may be specified.

An nd1-by-nt array containing the indices of either the convex hull (nd1 is equal to nd), or the Delaunay triangulation (nd1 is equal to nd+1) of the input points.

Keywords

BOUNDS

Set this keyword equal to a named variable that will contain the indices of the vertices of the convex hull.

CONNECTIVITY

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

Each element i, 0 £ i < np, contains the starting index of the connectivity list (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 not ordered. To obtain the connectivity list, either the DELAUNAY or SPHERE keywords must also be specified.

For example, to perform a spherical triangulation, use the following procedure call:

QHULL, lon, lat, CONNECTIVITY = list, SHPERE = sphere

which returns the adjacency list in the variable list. The subscripts of the nodes adjacent to lon[i] and lat[i] are contained in the array: list[list[i] :list[i+1] - 1].

DELAUNAY

Performs a Delaunay triangulation and returns the vertex indices of the resulting polyhedra; otherwise, the convex hull of the data are returned.

SPHERE

Computes the Delaunay triangulation of the points which lie on the surface of a sphere. The V0 argument contains the longitude, in degrees, and V1 contains the latitude, in degrees, of each point.

VDIAGRAM

For two-dimensional arrays, VDIAGRAM is a 4-by-nv integer array. For each Voronoi ridge, i, VDIAGRAM[0:1, i] contains the index of the two input points the ridge bisects. VDIAGRAM [2:3, i] contains the indices within VVERTICES of the Voronoi vertices. In the case of an unbounded half-space, VDIAGRAM[2, i] is set to a negative index, j, indicating that the corresponding Voronoi ridge is unbounded, and that the equation for the ridge is contained in VNORMAL[*, -j-1], and starts at Voronoi vertex [3, i].

For three-dimensional or higher dimensional arrays, VDIAGRAM is returned as a connectivity vector. This vector has the form [n, v0, v1, i0, i1, ..., in-3], where n is the number of points needed to describe that particular Voronoi ridge, v0 and v1 contain the indices for the two input points that the ridge bisects, and i0...in -3 contain the indices within VVERTICES of the Voronoi vertices. In the case of an unbounded half-space, VDIAGRAM[i] is set to a negative index, j, indicating that the corresponding Voronoi ridge is unbounded, and that the equation for the ridge is contained in VNORMAL[*, -j-1], and starts at VDIAGRAM[i+1].

VNORMALS

When specified, this keyword returns the normals of each Voronoi ridge that is unbounded. The normals consist of a (nd+1)-by-nu array, where nd is the number of dimensions, and nu is the number of unbounded vertices. The first nd elements in each row contain the equation for the unbounded ridge, normalized by the last element in each row. See the preceding description of VDIAGRAM for details.

VVERTICES

Examples

pro ex_qhull 
 
    ; Create a collection of random points. 
    n = 20 
    seed = 15 
    x = RANDOMU(seed, n) 
    y = RANDOMU(seed, n) 
 
    ; Construct the Delaunay triangulation 
    ; and the Voronoi diagram. 
    QHULL, x, y, triangle, /DELAUNAY, $ 
        VDIAGRAM=vdiag, VVERTICES=vvert, VNORM=vnorm 
 
    ; Plot our input points. 
    PLOT, [-0.1, 1.1], [-0.1, 1.1], /NODATA, $ 
        XSTYLE=4, YSTYLE=4 
    PLOTS, x, y, PSYM=4 
 
    ; Plot the Voronoi diagram. 
    for i=0,N_ELEMENTS(vdiag[2,*])-1 do begin 
        j = vdiag[2,i] + 1 
        ; Bounded or unbounded? 
        if (j gt 0) then begin   ; Bounded. 
            PLOTS, vvert[*, vdiag[2:3,i]], PSYM=-5 
        endif else begin 
            ; Unbounded, retrieve normal equation. 
            ; Vnorm[2,*] contains the normalization parameter. 
            offset = [vnorm[1,-j], -vnorm[0,-j]]*vnorm[2,-j] 
            xy1 = vvert[*, vdiag[3,i]] 
            ; Reverse the normal direction if necessary. 
            if (xy1[1] lt xy1[0]) then offset = -offset 
            PLOTS, [[xy1], [xy1 + 5*offset]] 
        endelse 
    endfor 
 
end