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