PRO triangulate_extend, xc, yc, ext, tr, ntr
; Description: This module generates an extended Delaunay triangulation of a set of unique points
; in two dimensions. As a first step, the module generates the Delaunay triangulation
; of the set of points with x and y coordinates "xc" and "yc", respectively, which
; is the unique set of disjointed triangles formed such that the circumcircle of any
; triangle contains no vertices from any other triangle.
; The Delaunay set of triangles is then extended as follows. For each point P, the
; module determines the set of points Sp that can be connected to P via a maximum
; number "ext" of edges of the Delaunay triangulation and forms all possible triangles
; containing P from the points in Sp. These triangles are appended to the Delaunay
; triangulation.
; Now let Q be a point from Sp, but different from P, and let Sq be the set of
; points that can be connected to Q via a maximum number "ext" of edges of the
; Delaunay triangulation. Also let R be a point from Sq that is not in Sp. Then, for
; each point P, the module forms all possible triangles from the points P, Q and R,
; and appends them to the Delaunay triangulation.
; This algorithm for generating an extended Delaunay triangulation guarantees that
; each triangle is unique. The development of this module was based on ideas in the
; paper by Pal & Bakos (Pal A. & Bakos G.A., 2006, PASP, 118, 1474). However, their
; description of how to generate an extended Delaunay triangulation should generate
; more triangles for a particular extension level than our algorithm, although in
; reality our algorithm generates more triangles for a particular extension level
; than the numbers of triangles reported in Pal & Bakos 2006. This implies that their
; implementation of the extended Delaunay triangulation is not as thorough as the one
; implemented in our algorithm.
; For the final set of triangles, the indices of each of the three vertices for each
; triangle are returned via the array "tr", and the number of triangles that have been
; formed are returned via the parameter "ntr".
;
; Input Parameters:
;
; xc - FLOAT/DOUBLE VECTOR - A one-dimensional vector of x coordinates for the set of points.
; yc - FLOAT/DOUBLE VECTOR - A one-dimensional vector of y coordinates for the set of points.
; ext - INTEGER/LONG - The extension level of the extended Delaunay triangulation (Default
; value = 0). The value "0" corresponds to a Delaunay triangulation. A
; positive integer indicates that for each point P, an extended triangulation
; is generated as described above. A negative value for this parameter is
; reset to the value "0".
;
; N.B: The coordinate list must not contain duplicate entries, otherwise the extended triangle
; generation algorithm will fail.
;
; Output Parameters:
;
; tr - LONG ARRAY - A 3 by "ntr" array of LONGS where each triplet of numbers represents the
; indices of the points corresponding to the three vertices of each triangle.
; ntr - LONG - The number of triangles that have been formed.
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
; 01/08/2008 - Module created (dmb)
;Check that "xc" and "yc" are one-dimensional vectors of the correct number type, and that they
;the same number of elements, which must be at least three in any case
;Check that "ext" is a non-negative number of the correct type, and reset the value to "0" if it
;is out of range
;Generate the Delaunay triangulation of the set of points
;If the value of "ext" is zero, then return with the Delaunay triangulation
;For each point
;Determine the group of points connected to the current point by a maximum of "ext" edges of the
;Delaunay triangulation
;Save the group of points for later use
;Determine the set of unique sizes of the groups of points connected to each point by a maximum
;of "ext" edges of the Delaunay triangulation
;Sort the unique group sizes into descending order
;If there is only one group size that is the same as the number of input points, then the extended
;Delaunay triangulation is equivalent to the full triangulation of the input points. Consequently,
;in this case, return the full triangulation.
;For each point P
;If there is more than one point in the set of points Sp connected to P by a maximum of "ext" edges of
;the Delaunay triangulation
;If necessary, set up a triangle array to hold new triangles
;For each point Q, different from P, in the set Sp
;Extract the point Q
;If there are more than two points in Sp
;Form all possible triangles containing P and Q from the points in Sp
;If there is more than one point in the set of points Sq connected to Q by a maximum of "ext" edges
;of the Delaunay triangulation
;For each point R, different from Q, in the set Sq
;Extract the point R
;If the point R is not in Sp
;Save the triangle formed by P, Q and R
;Store the set of triangles generated for point P
;Remove P from all groups in which it is present
;Extract the point Q for which P is in Sq
;If there is at least one point left in Sq
;Determine the subscripts of the points in Sq that are different to P
;If P is in Sq but there is also at least one point in Sq that is different to P
;Save only the points in Sq that are different to P
;If there are no points in Sq that are different to P