moved here from GRASS Addons

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@34517 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/v.delaunay2/Makefile

@@ -0,0 +1,12 @@
+MODULE_TOPDIR = ../..
+
+PGM=v.delaunay
+
+LIBES     = $(VECTLIB) $(VECTLIB_REAL) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+ 
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd    

vector/v.delaunay2/data_types.h

@@ -0,0 +1,34 @@
+#ifndef DATA_TYPES_H
+#define DATA_TYPES_H
+
+#include <stdlib.h>
+
+#ifndef MY_NULL
+#define MY_NULL  0
+#endif
+#define TRUE  1
+#define FALSE 0
+
+typedef enum
+{ left, right } side;
+
+typedef unsigned char boolean;
+
+struct vertex
+{
+    double x, y, z;
+    struct edge *entry_pt;
+};
+
+struct edge
+{
+    struct vertex *org;
+    struct vertex *dest;
+    struct edge *onext;
+    struct edge *oprev;
+    struct edge *dnext;
+    struct edge *dprev;
+};
+
+extern struct vertex *sites;
+#endif

vector/v.delaunay2/defs.h

@@ -0,0 +1,11 @@
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+
+#ifdef MAIN
+struct Cell_head Window;
+BOUND_BOX Box;
+#else
+extern struct Cell_head Window;
+extern BOUND_BOX Box;
+#endif

vector/v.delaunay2/edge.c

@@ -0,0 +1,126 @@
+#include  "data_types.h"
+#include  "memory.h"
+#include  "edge.h"
+
+/* 
+ *  Construct an edge from vertices v1, v2 and add it to rings of edges e1, e2
+ */
+struct edge *join(struct edge *e1, struct vertex *v1,
+		  struct edge *e2, struct vertex *v2, side s)
+{
+
+    struct edge *new_edge;
+
+    /* v1, v2 - vertices to be joined.
+       e1, e2 - edges to which v1, v2 belong to */
+
+    new_edge = create_edge(v1, v2);
+
+    if (s == left) {
+	if (ORG(e1) == v1)
+	    splice(OPREV(e1), new_edge, v1);
+	else
+	    splice(DPREV(e1), new_edge, v1);
+	splice(e2, new_edge, v2);
+    }
+    else {
+	splice(e1, new_edge, v1);
+	if (ORG(e2) == v2)
+	    splice(OPREV(e2), new_edge, v2);
+	else
+	    splice(DPREV(e2), new_edge, v2);
+    }
+    return new_edge;
+}
+
+/* 
+ *  Remove an edge.
+ */
+void delete_edge(struct edge *e)
+{
+    struct vertex *u, *v;
+
+    /* Save destination and origin. */
+    u = ORG(e);
+    v = DEST(e);
+
+    /* Set entry points. */
+    if (u->entry_pt == e)
+	u->entry_pt = ONEXT(e);
+    if (v->entry_pt == e)
+	v->entry_pt = DNEXT(e);
+
+    /* Four edge references need adjustment */
+    if (ORG(ONEXT(e)) == u)
+	OPREV(ONEXT(e)) = OPREV(e);
+    else
+	DPREV(ONEXT(e)) = OPREV(e);
+
+    if (ORG(OPREV(e)) == u)
+	ONEXT(OPREV(e)) = ONEXT(e);
+    else
+	DNEXT(OPREV(e)) = ONEXT(e);
+
+    if (ORG(DNEXT(e)) == v)
+	OPREV(DNEXT(e)) = DPREV(e);
+    else
+	DPREV(DNEXT(e)) = DPREV(e);
+
+    if (ORG(DPREV(e)) == v)
+	ONEXT(DPREV(e)) = DNEXT(e);
+    else
+	DNEXT(DPREV(e)) = DNEXT(e);
+
+    free_edge(e);
+}
+
+
+ /*  Add an edge to a ring of edges. */
+void splice(struct edge *a, struct edge *b, struct vertex *v)
+{
+    struct edge *next;
+
+    /* b must be the unnattached edge and a must be the previous 
+       ccw edge to b. */
+
+    if (ORG(a) == v) {
+	next = ONEXT(a);
+	ONEXT(a) = b;
+    }
+    else {
+	next = DNEXT(a);
+	DNEXT(a) = b;
+    }
+
+    if (ORG(next) == v)
+	OPREV(next) = b;
+    else
+	DPREV(next) = b;
+
+    if (ORG(b) == v) {
+	ONEXT(b) = next;
+	OPREV(b) = a;
+    }
+    else {
+	DNEXT(b) = next;
+	DPREV(b) = a;
+    }
+}
+
+ /*  Create a new edge and initialize it */
+struct edge *create_edge(struct vertex *v1, struct vertex *v2)
+{
+    struct edge *new_edge;
+
+    new_edge = get_edge();
+
+    DNEXT(new_edge) = DPREV(new_edge) = ONEXT(new_edge) = OPREV(new_edge) =
+	new_edge;
+    ORG(new_edge) = v1;
+    DEST(new_edge) = v2;
+    if (v1->entry_pt == MY_NULL)
+	v1->entry_pt = new_edge;
+    if (v2->entry_pt == MY_NULL)
+	v2->entry_pt = new_edge;
+    return new_edge;
+}

vector/v.delaunay2/edge.h

@@ -0,0 +1,23 @@
+#ifndef EDGE_H
+#define EDGE_H
+
+#define ORG(e)   ((e)->org)
+#define DEST(e)  ((e)->dest)
+#define ONEXT(e) ((e)->onext)
+#define OPREV(e) ((e)->oprev)
+#define DNEXT(e) ((e)->dnext)
+#define DPREV(e) ((e)->dprev)
+
+#define OTHER_VERTEX(e,p) (ORG(e) == p ? DEST(e) : ORG(e))
+#define NEXT(e,p)         (ORG(e) == p ? ONEXT(e) : DNEXT(e))
+#define PREV(e,p)         (ORG(e) == p ? OPREV(e) : DPREV(e))
+
+#define SAME_EDGE(e1,e2) (e1 == e2)
+
+struct edge *join(struct edge *e1, struct vertex *v1,
+		  struct edge *e2, struct vertex *v2, side s);
+void delete_edge(struct edge *e);
+void splice(struct edge *a, struct edge *b, struct vertex *v);
+struct edge *create_edge(struct vertex *v1, struct vertex *v2);
+
+#endif

vector/v.delaunay2/geom_primitives.h

@@ -0,0 +1,23 @@
+#ifndef GEOM_PRIMITIVES_H
+#define GEOM_PRIMITIVES_H
+
+#define CREATE_VECTOR(p1, p2, dx, dy) (dx = p2->x - p1->x, dy = p2->y - p1->y)
+
+#define DOT_PRODUCT_2V(a1, a2, b1, b2) (a1 * b1 + a2 * b2)
+
+#define CROSS_PRODUCT_2V(a1, a2, b1, b2) (a1 * b2 - a2 * b1)
+/*
+   cross-product around p2
+   ((p2->x - p1->x) * (p3->y - p2->y) - (p2->y - p1->y) * (p3->x - p2->x))
+
+   around p1
+ */
+#define CROSS_PRODUCT_3P(p1, p2, p3) ((p2->x - p1->x) * (p3->y - p1->y) - (p2->y - p1->y) * (p3->x - p1->x))
+
+/* predicate testing if p3 is to the left of the line through p1 and p2 */
+#define LEFT_OF(p1, p2, p3) (CROSS_PRODUCT_3P(p1, p2, p3) > 0)
+
+/* predicate testing if p3 is to the right of the line through p1 and p2 */
+#define RIGHT_OF(p1, p2, p3) (CROSS_PRODUCT_3P(p1, p2, p3) < 0)
+
+#endif

vector/v.delaunay2/geometry.c

@@ -0,0 +1,287 @@
+#include <stddef.h>
+#include "data_types.h"
+#include "memory.h"
+#include "geometry.h"
+#include "geom_primitives.h"
+#include "edge.h"
+
+static void find_lowest_cross_edge(struct edge *r_cw_l, struct vertex *s,
+				   struct edge *l_ccw_r, struct vertex *u,
+				   struct edge **l_lower,
+				   struct vertex **org_l_lower,
+				   struct edge **r_lower,
+				   struct vertex **org_r_lower);
+
+static void merge(struct edge *r_cw_l, struct vertex *s,
+		  struct edge *l_ccw_r, struct vertex *u,
+		  struct edge **l_tangent);
+
+void divide(struct vertex *sites_sorted[], unsigned int l, unsigned int r,
+	    struct edge **l_ccw, struct edge **r_cw)
+{
+
+    unsigned int n;
+    unsigned int split;
+    struct edge *l_ccw_l, *r_cw_l, *l_ccw_r, *r_cw_r, *l_tangent;
+    struct edge *a, *b, *c;
+    double c_p;
+
+    n = r - l + 1;
+    if (n == 2) {
+	/* Base case #1 - 2 sites in region. Construct an edge from 
+	   two sites in the region       */
+	*l_ccw = *r_cw = create_edge(sites_sorted[l], sites_sorted[r]);
+    }
+    else if (n == 3) {
+	/* Base case #2 - 3 sites. Construct a triangle or two edges */
+	a = create_edge(sites_sorted[l], sites_sorted[l + 1]);
+	b = create_edge(sites_sorted[l + 1], sites_sorted[r]);
+	splice(a, b, sites_sorted[l + 1]);
+	c_p =
+	    CROSS_PRODUCT_3P(sites_sorted[l], sites_sorted[l + 1],
+			     sites_sorted[r]);
+
+	if (c_p > 0.0) {
+	    /* Create a triangle */
+	    c = join(a, sites_sorted[l], b, sites_sorted[r], right);
+	    *l_ccw = a;
+	    *r_cw = b;
+	}
+	else if (c_p < 0.0) {
+	    /* Create a triangle */
+	    c = join(a, sites_sorted[l], b, sites_sorted[r], left);
+	    *l_ccw = c;
+	    *r_cw = c;
+	}
+	else {
+	    /* Cross product is zero. Sites are located on a line. 
+	       Triangle cannot be created */
+	    *l_ccw = a;
+	    *r_cw = b;
+	}
+    }
+    else if (n > 3) {
+	/* Recursive case. Continue splitting */
+
+	/* Splitting point */
+	split = (l + r) / 2;
+
+	/* Divide into two halves */
+	divide(sites_sorted, l, split, &l_ccw_l, &r_cw_l);
+	divide(sites_sorted, split + 1, r, &l_ccw_r, &r_cw_r);
+
+	/* Merge the two triangulations */
+	merge(r_cw_l, sites_sorted[split], l_ccw_r, sites_sorted[split + 1],
+	      &l_tangent);
+
+	/* The lower tangent added by merge may have invalidated 
+	   l_ccw_l or r_cw_r. Update them if necessary. */
+	if (ORG(l_tangent) == sites_sorted[l])
+	    l_ccw_l = l_tangent;
+	if (DEST(l_tangent) == sites_sorted[r])
+	    r_cw_r = l_tangent;
+
+	/* Update leftmost ccw edge and rightmost cw edge */
+	*l_ccw = l_ccw_l;
+	*r_cw = r_cw_r;
+    }
+}
+
+/*
+ *  Find the lowest cross edge of the two triangulations
+ */
+static void find_lowest_cross_edge(struct edge *r_cw_l, struct vertex *s,
+				   struct edge *l_ccw_r, struct vertex *u,
+				   struct edge **l_lower,
+				   struct vertex **org_l_lower,
+				   struct edge **r_lower,
+				   struct vertex **org_r_lower)
+{
+    struct edge *l, *r;
+    struct vertex *o_l, *o_r, *d_l, *d_r;
+    boolean ready;
+
+    l = r_cw_l;
+    r = l_ccw_r;
+    o_l = s;
+    d_l = OTHER_VERTEX(l, s);
+    o_r = u;
+    d_r = OTHER_VERTEX(r, u);
+    ready = FALSE;
+
+    while (!ready)
+	/* left_of */
+	if (LEFT_OF(o_l, d_l, o_r)) {
+	    l = PREV(l, d_l);
+	    o_l = d_l;
+	    d_l = OTHER_VERTEX(l, o_l);
+	    /* right_of */
+	}
+	else if (RIGHT_OF(o_r, d_r, o_l)) {
+	    r = NEXT(r, d_r);
+	    o_r = d_r;
+	    d_r = OTHER_VERTEX(r, o_r);
+	}
+	else
+	    ready = TRUE;
+
+    *l_lower = l;
+    *r_lower = r;
+    *org_l_lower = o_l;
+    *org_r_lower = o_r;
+}
+
+/* 
+ *  The most time-expensive function, most of the work gets done here.
+ */
+static void merge(struct edge *r_cw_l, struct vertex *s,
+		  struct edge *l_ccw_r, struct vertex *u,
+		  struct edge **l_tangent)
+{
+    struct edge *base, *l_cand, *r_cand;
+    struct vertex *org_base, *dest_base;
+    long double u_l_c_o_b, v_l_c_o_b, u_l_c_d_b, v_l_c_d_b;
+    long double u_r_c_o_b, v_r_c_o_b, u_r_c_d_b, v_r_c_d_b;
+
+    /* cross product */
+    long double c_p_l_cand, c_p_r_cand;
+
+    /* dot product */
+    long double d_p_l_cand, d_p_r_cand;
+    boolean above_l_cand, above_r_cand, above_next, above_prev;
+    struct vertex *dest_l_cand, *dest_r_cand;
+    long double cot_l_cand, cot_r_cand;
+    struct edge *l_lower, *r_lower;
+    struct vertex *org_r_lower, *org_l_lower;
+
+    /* Create first cross edge by joining lower common tangent */
+    find_lowest_cross_edge(r_cw_l, s, l_ccw_r, u, &l_lower, &org_l_lower,
+			   &r_lower, &org_r_lower);
+    base = join(l_lower, org_l_lower, r_lower, org_r_lower, right);
+    org_base = org_l_lower;
+    dest_base = org_r_lower;
+
+    /* Need to return lower tangent. */
+    *l_tangent = base;
+
+    /* The merge loop */
+    while (TRUE) {
+	/* Initialise edges l_cand and r_cand */
+	l_cand = NEXT(base, org_base);
+	r_cand = PREV(base, dest_base);
+	dest_l_cand = OTHER_VERTEX(l_cand, org_base);
+	dest_r_cand = OTHER_VERTEX(r_cand, dest_base);
+
+	/* Vectors used for above and modified IN_CIRCLE tests
+	   u/v left/right candidate origin/destination */
+	CREATE_VECTOR(dest_l_cand, org_base, u_l_c_o_b, v_l_c_o_b);
+	CREATE_VECTOR(dest_l_cand, dest_base, u_l_c_d_b, v_l_c_d_b);
+	CREATE_VECTOR(dest_r_cand, org_base, u_r_c_o_b, v_r_c_o_b);
+	CREATE_VECTOR(dest_r_cand, dest_base, u_r_c_d_b, v_r_c_d_b);
+
+	/* Above tests. */
+	c_p_l_cand =
+	    CROSS_PRODUCT_2V(u_l_c_o_b, v_l_c_o_b, u_l_c_d_b, v_l_c_d_b);
+	c_p_r_cand =
+	    CROSS_PRODUCT_2V(u_r_c_o_b, v_r_c_o_b, u_r_c_d_b, v_r_c_d_b);
+	above_l_cand = c_p_l_cand > 0.0;
+	above_r_cand = c_p_r_cand > 0.0;
+
+	/* Terminate merge loop. No valid sites left in L or R.
+	   The top-most cross-edge have already been added. */
+	if (!above_l_cand && !above_r_cand)
+	    break;
+
+	/* Move to next l_cand ccw, delete the old l_cand edge, until the 
+	   in_circle test gets invalid. */
+	if (above_l_cand) {
+	    double u_n_o_b, v_n_o_b, u_n_d_b, v_n_d_b;
+	    double c_p_next, d_p_next, cot_next;
+	    struct edge *next;
+	    struct vertex *dest_next;
+
+	    d_p_l_cand =
+		DOT_PRODUCT_2V(u_l_c_o_b, v_l_c_o_b, u_l_c_d_b, v_l_c_d_b);
+	    cot_l_cand = d_p_l_cand / c_p_l_cand;
+
+	    while (TRUE) {
+		next = NEXT(l_cand, org_base);
+		dest_next = OTHER_VERTEX(next, org_base);
+		CREATE_VECTOR(dest_next, org_base, u_n_o_b, v_n_o_b);
+		CREATE_VECTOR(dest_next, dest_base, u_n_d_b, v_n_d_b);
+		c_p_next =
+		    CROSS_PRODUCT_2V(u_n_o_b, v_n_o_b, u_n_d_b, v_n_d_b);
+		above_next = c_p_next > 0.0;
+
+		if (!above_next)
+		    break;	/* Terminate loop. */
+
+		d_p_next = DOT_PRODUCT_2V(u_n_o_b, v_n_o_b, u_n_d_b, v_n_d_b);
+		cot_next = d_p_next / c_p_next;
+
+		if (cot_next > cot_l_cand)
+		    break;	/* Terminate loop. */
+
+		delete_edge(l_cand);
+		l_cand = next;
+		cot_l_cand = cot_next;
+	    }
+	}
+
+	/* Essentially the same done for r_cand symmetrically. */
+	/* Move to prev r_cand cw, delete the old r_cand edge, until the 
+	   in_circle test gets invalid. */
+	if (above_r_cand) {
+	    double u_p_o_b, v_p_o_b, u_p_d_b, v_p_d_b;
+	    double c_p_prev, d_p_prev, cot_prev;
+	    struct edge *prev;
+	    struct vertex *dest_prev;
+
+	    d_p_r_cand =
+		DOT_PRODUCT_2V(u_r_c_o_b, v_r_c_o_b, u_r_c_d_b, v_r_c_d_b);
+	    cot_r_cand = d_p_r_cand / c_p_r_cand;
+
+	    while (TRUE) {
+		prev = PREV(r_cand, dest_base);
+		dest_prev = OTHER_VERTEX(prev, dest_base);
+		CREATE_VECTOR(dest_prev, org_base, u_p_o_b, v_p_o_b);
+		CREATE_VECTOR(dest_prev, dest_base, u_p_d_b, v_p_d_b);
+		c_p_prev =
+		    CROSS_PRODUCT_2V(u_p_o_b, v_p_o_b, u_p_d_b, v_p_d_b);
+		above_prev = c_p_prev > 0.0;
+
+		if (!above_prev)
+		    break;	/* Terminate. */
+
+		d_p_prev = DOT_PRODUCT_2V(u_p_o_b, v_p_o_b, u_p_d_b, v_p_d_b);
+		cot_prev = d_p_prev / c_p_prev;
+
+		if (cot_prev > cot_r_cand)
+		    break;	/* Terminate. */
+
+		delete_edge(r_cand);
+		r_cand = prev;
+		cot_r_cand = cot_prev;
+	    }
+	}
+
+	/*
+	   Add a successive L-R cross edge from base 
+	   If both candidates are valid, choose the one with smallest circumcircle. 
+	   Set the base to the new L-R cross edge.
+	 */
+	dest_l_cand = OTHER_VERTEX(l_cand, org_base);
+	dest_r_cand = OTHER_VERTEX(r_cand, dest_base);
+	if (!above_l_cand ||
+	    (above_l_cand && above_r_cand && cot_r_cand < cot_l_cand)) {
+	    /* Connect to the right */
+	    base = join(base, org_base, r_cand, dest_r_cand, right);
+	    dest_base = dest_r_cand;
+	}
+	else {
+	    /* Connect to the left */
+	    base = join(l_cand, dest_l_cand, base, dest_base, right);
+	    org_base = dest_l_cand;
+	}
+    }
+}

vector/v.delaunay2/geometry.h

@@ -0,0 +1,6 @@
+#ifndef GEOMETRY_H
+#define GEOMETRY_H
+
+void divide(struct vertex *sites_sorted[], unsigned int l, unsigned int r,
+	    struct edge **l_ccw, struct edge **r_cw);
+#endif

vector/v.delaunay2/in_out.c

@@ -0,0 +1,207 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+#include "data_types.h"
+#include "memory.h"
+#include "edge.h"
+#include "geom_primitives.h"
+
+void output_edges(struct vertex *sites_sorted[], unsigned int n, int mode3d,
+		  int Type, struct Map_info map_out)
+{
+    struct edge *e_start, *e;
+    struct vertex *u, *v;
+    unsigned int i;
+
+    static struct line_pnts *Points = NULL;
+    static struct line_cats *Cats = NULL;
+
+    if (!Points) {
+	Points = Vect_new_line_struct();
+	Cats = Vect_new_cats_struct();
+    }
+    double x1, y1, z1, x2, y2, z2;
+
+    for (i = 0; i < n; i++) {
+	u = sites_sorted[i];
+	e_start = e = u->entry_pt;
+	do {
+	    v = OTHER_VERTEX(e, u);
+	    if (u < v) {
+		x1 = sites[u - sites].x;
+		y1 = sites[u - sites].y;
+		x2 = sites[v - sites].x;
+		y2 = sites[v - sites].y;
+
+		Vect_reset_line(Points);
+
+		if (mode3d) {
+		    z1 = sites[u - sites].z;
+		    z2 = sites[v - sites].z;
+		    Vect_append_point(Points, x1, y1, z1);
+		    Vect_append_point(Points, x2, y2, z2);
+		}
+		else {
+		    Vect_append_point(Points, x1, y1, 0.0);
+		    Vect_append_point(Points, x2, y2, 0.0);
+		}
+		Vect_write_line(&map_out, Type, Points, Cats);
+	    }
+	    e = NEXT(e, u);
+	} while (!SAME_EDGE(e, e_start));
+    }
+}
+
+/* Print the ring of triangles about each vertex. */
+
+void output_triangles(struct vertex *sites_sorted[], unsigned int n,
+		      int mode3d, int Type, struct Map_info map_out)
+{
+    struct edge *e_start, *e, *next;
+    struct vertex *u, *v, *w;
+    unsigned int i;
+    struct vertex *temp;
+
+    struct line_pnts *Points = Vect_new_line_struct();
+    struct line_cats *Cats = Vect_new_cats_struct();
+
+    double x1, y1, z1, x2, y2, z2, x3, y3, z3;
+
+    for (i = 0; i < n; i++) {
+	u = sites_sorted[i];
+	e_start = e = u->entry_pt;
+	do {
+	    v = OTHER_VERTEX(e, u);
+	    if (u < v) {
+		next = NEXT(e, u);
+		w = OTHER_VERTEX(next, u);
+		if (u < w)
+		    if (SAME_EDGE(NEXT(next, w), PREV(e, v))) {
+			/* Triangle. */
+			if (v > w) {
+			    temp = v;
+			    v = w;
+			    w = temp;
+			}
+			x1 = sites[u - sites].x;
+			y1 = sites[u - sites].y;
+			x2 = sites[v - sites].x;
+			y2 = sites[v - sites].y;
+			x3 = sites[w - sites].x;
+			y3 = sites[w - sites].y;
+
+			if (mode3d) {
+			    z1 = sites[u - sites].z;
+			    z2 = sites[v - sites].z;
+			    z3 = sites[w - sites].z;
+			    Vect_reset_line(Points);
+			    Vect_append_point(Points, x1, y1, z1);
+			    Vect_append_point(Points, x2, y2, z2);
+			    Vect_write_line(&map_out, Type, Points, Cats);
+
+			    Vect_reset_line(Points);
+			    Vect_append_point(Points, x2, y2, z2);
+			    Vect_append_point(Points, x3, y3, z3);
+			    Vect_write_line(&map_out, Type, Points, Cats);
+
+			    Vect_reset_line(Points);
+			    Vect_append_point(Points, x3, y3, z3);
+			    Vect_append_point(Points, x1, y1, z1);
+			    Vect_write_line(&map_out, Type, Points, Cats);
+			}
+			else {
+			    Vect_reset_line(Points);
+			    Vect_append_point(Points, x1, y1, 0.0);
+			    Vect_append_point(Points, x2, y2, 0.0);
+			    Vect_write_line(&map_out, Type, Points, Cats);
+
+			    Vect_reset_line(Points);
+			    Vect_append_point(Points, x2, y2, 0.0);
+			    Vect_append_point(Points, x3, y3, 0.0);
+			    Vect_write_line(&map_out, Type, Points, Cats);
+
+			    Vect_reset_line(Points);
+			    Vect_append_point(Points, x3, y3, 0.0);
+			    Vect_append_point(Points, x1, y1, 0.0);
+			    Vect_write_line(&map_out, Type, Points, Cats);
+			}
+		    }
+	    }
+	    /* Next edge around u. */
+	    e = NEXT(e, u);
+	} while (!SAME_EDGE(e, e_start));
+    }
+}
+
+void remove_duplicates(struct vertex *list[], unsigned int *size)
+{
+    int n = *size - 1;
+    int left = 0;
+    int right = 1;
+    int shift = 0;
+    int empty = 0;
+
+    if (n > 0)
+	for (; right < n; left++, right++)
+	    if (list[left]->x == list[right]->x &&
+		list[left]->y == list[right]->y) {
+		if (shift == 0) {
+		    shift = 1;
+		    empty = right;
+		}
+		(*size)--;
+	    }
+	    else if (shift == 1)
+		list[empty++] = list[right];
+}
+
+/* returns number of sites read */
+int read_sites(int mode3d, int complete_map, struct Map_info map_in,
+	       BOUND_BOX Box)
+{
+    int nlines, line, allocated, nsites;
+    struct line_pnts *Points;
+
+    Points = Vect_new_line_struct();
+    nlines = Vect_get_num_lines(&map_in);
+    alloc_sites(nlines);
+    allocated = nlines;
+
+    nsites = 0;
+    for (line = 1; line <= nlines; line++) {
+	int type;
+
+	type = Vect_read_line(&map_in, Points, NULL, line);
+	if (!(type & GV_POINTS))
+	    continue;
+	if (!complete_map) {
+	    if (!Vect_point_in_box(Points->x[0], Points->y[0], 0.0, &Box))
+		continue;
+	}
+	sites[nsites].x = Points->x[0];
+	sites[nsites].y = Points->y[0];
+	if (mode3d) {
+	    G_debug(3, "Points->z[0]: %f", Points->z[0]);
+	    sites[nsites].z = Points->z[0];
+	}
+	/* Initialise entry edge vertices. */
+	sites[nsites].entry_pt = MY_NULL;
+
+	nsites += 1;
+	/* number 100 was arbitrarily chosen /*
+	   /*        if (nsites == allocated && line != nlines){
+	   allocated += 100;
+	   realloc_sites(allocated);
+	   }
+	 */
+    }
+    if (nsites != nlines)
+	realloc_sites(nsites);
+    alloc_edges(nsites);
+
+    return nsites;
+}

vector/v.delaunay2/in_out.h

@@ -0,0 +1,11 @@
+#ifndef IN_OUT_H
+#define IN_OUT_H
+
+int read_sites(int mode3d, int complete_map, struct Map_info map_in,
+	       BOUND_BOX Box);
+void output_edges(struct vertex *sites_sorted[], unsigned int n, int mode3d,
+		  int Type, struct Map_info map_out);
+void output_triangles(struct vertex *sites_sorted[], unsigned int n,
+		      int mode3d, int Type, struct Map_info map_out);
+void remove_duplicates(struct vertex *list[], unsigned int *size);
+#endif

vector/v.delaunay2/main.c

@@ -0,0 +1,164 @@
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+#include "data_types.h"
+#include "memory.h"
+#include "geometry.h"
+#include "edge.h"
+#include "in_out.h"
+
+int compare(const struct vertex **p1, const struct vertex **p2);
+
+int main(int argc, char *argv[])
+{
+
+    /* GRASS related variables */
+    char *mapset;
+    struct Map_info map_in, map_out;
+    struct Cell_head Window;
+    BOUND_BOX Box;
+    struct GModule *module;
+    struct Flag *reg_flag, *line_flag;
+    struct Option *in_opt, *out_opt;
+
+    struct line_pnts *Points;
+    struct line_cats *Cats;
+    int nareas, area;
+
+    int Type;
+    int complete_map;
+    int mode3d;
+
+    /* ---------------------- */
+
+    unsigned int i;
+    unsigned int n;
+    struct edge *l_cw, *r_ccw;
+    struct vertex **sites_sorted;
+
+    /* GRASS related manipulations */
+    G_gisinit(argv[0]);
+    module = G_define_module();
+    module->keywords = _("vector");
+    module->description = _("Creates a Delaunay triangulation from an input "
+			    "vector map containing points or centroids.");
+
+    in_opt = G_define_standard_option(G_OPT_V_INPUT);
+    out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+
+    reg_flag = G_define_flag();
+    reg_flag->key = 'r';
+    reg_flag->description = _("Use only points in current region");
+
+    line_flag = G_define_flag();
+    line_flag->key = 'l';
+    line_flag->description =
+	_("Output triangulation as a graph (lines), not areas");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    if (line_flag->answer)
+	Type = GV_LINE;
+    else
+	Type = GV_BOUNDARY;
+
+    complete_map = reg_flag->answer ? 0 : 1;
+
+    Points = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+
+    /* open files */
+    if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
+	G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
+
+    Vect_set_open_level(2);
+    Vect_open_old(&map_in, in_opt->answer, mapset);
+
+    /* check if we have a 3D input points map */
+    mode3d = Vect_is_3d(&map_in);
+
+    if (mode3d) {
+	if (0 > Vect_open_new(&map_out, out_opt->answer, 1))
+	    G_fatal_error(_("Unable to create vector map <%s>"),
+			  out_opt->answer);
+    }
+    else if (0 > Vect_open_new(&map_out, out_opt->answer, 0))
+	G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
+
+    Vect_hist_copy(&map_in, &map_out);
+    Vect_hist_command(&map_out);
+
+    /* initialize working region */
+    G_get_window(&Window);
+    G_percent(0, 100, 1);
+    Vect_region_box(&Window, &Box);
+
+    n = read_sites(mode3d, complete_map, map_in, Box);
+
+    /* Sort. */
+    sites_sorted =
+	(struct vertex **)G_malloc((unsigned)n * sizeof(struct vertex *));
+    if (sites_sorted == MY_NULL)
+	G_fatal_error(_("Not enough memory."));
+    for (i = 0; i < n; i++)
+	sites_sorted[i] = sites + i;
+    qsort(sites_sorted, n, sizeof(struct vertex *), (void *)compare);
+
+    remove_duplicates(sites_sorted, &n);
+
+    /* Triangulate. */
+    divide(sites_sorted, 0, n - 1, &l_cw, &r_ccw);
+
+    output_edges(sites_sorted, n, mode3d, Type, map_out);
+
+    free((char *)sites_sorted);
+    free_memory();
+
+    Vect_close(&map_in);
+
+    if (Type == GV_BOUNDARY) {
+	Vect_build_partial(&map_out, GV_BUILD_AREAS);
+	nareas = Vect_get_num_areas(&map_out);
+	G_debug(3, "nareas = %d", nareas);
+	/*  Assign centroid to each area */
+	for (area = 1; area <= nareas; area++) {
+	    double x, y, z, angle, slope;
+	    int ret;
+
+	    G_percent(area, nareas, 2);
+	    Vect_reset_line(Points);
+	    Vect_reset_cats(Cats);
+	    ret = Vect_get_point_in_area(&map_out, area, &x, &y);
+	    if (ret < 0) {
+		G_warning(_("Unable to calculate area centroid"));
+		continue;
+	    }
+	    ret = Vect_tin_get_z(&map_out, x, y, &z, &angle, &slope);
+	    G_debug(3, "area centroid z: %f", z);
+	    if (ret < 0) {
+		G_warning(_("Unable to calculate area centroid z coordinate"));
+		continue;
+	    }
+	    Vect_append_point(Points, x, y, z);
+	    Vect_cat_set(Cats, 1, area);
+	    Vect_write_line(&map_out, GV_CENTROID, Points, Cats);
+	}
+    }
+    Vect_build(&map_out);
+    Vect_close(&map_out);
+
+    exit(EXIT_SUCCESS);
+}
+
+/* compare first according to x-coordinate, if equal then y-coordinate */
+int compare(const struct vertex **p1, const struct vertex **p2)
+{
+    if ((*p1)->x == (*p2)->x)
+	return ((*p1)->y < (*p2)->y);
+    return ((*p1)->x < (*p2)->x);
+}

vector/v.delaunay2/memory.c

@@ -0,0 +1,88 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+#include "data_types.h"
+
+struct vertex *sites;
+static struct edge *edges;
+static struct edge **free_list_e;
+
+static unsigned int n_free_e;
+
+void alloc_memory(unsigned int n)
+{
+    struct edge *e;
+    int i;
+
+    /* Sites storage. */
+    sites = (struct vertex *)G_calloc(n, sizeof(struct vertex));
+    if (sites == MY_NULL)
+	G_fatal_error(_("Not enough memory."));
+
+    /* Edges. Euler's formula - at most 3n edges on a set of n sites */
+    n_free_e = 3 * n;
+    edges = e = (struct edge *)G_calloc(n_free_e, sizeof(struct edge));
+    if (edges == MY_NULL)
+	G_fatal_error(_("Not enough memory."));
+
+    free_list_e = (struct edge **)G_calloc(n_free_e, sizeof(struct edge *));
+    if (free_list_e == MY_NULL)
+	G_fatal_error(_("Not enough memory."));
+    for (i = 0; i < n_free_e; i++, e++)
+	free_list_e[i] = e;
+}
+
+void alloc_sites(unsigned int n)
+{
+    /* Sites storage. */
+    sites = (struct vertex *)G_calloc(n, sizeof(struct vertex));
+    if (sites == MY_NULL)
+	G_fatal_error(_("Not enough memory."));
+}
+
+void realloc_sites(unsigned int n)
+{
+    /* Sites storage. */
+    sites = (struct vertex *)G_realloc(sites, n * sizeof(struct vertex));
+    if (sites == MY_NULL)
+	G_fatal_error(_("Not enough memory."));
+}
+
+void alloc_edges(unsigned int n)
+{
+    struct edge *e;
+    int i;
+
+    /* Edges. Euler's formula - at most 3n edges on a set of n sites */
+    n_free_e = 3 * n;
+    edges = e = (struct edge *)G_calloc(n_free_e, sizeof(struct edge));
+    if (edges == MY_NULL)
+	G_fatal_error(_("Not enough memory."));
+
+    free_list_e = (struct edge **)G_calloc(n_free_e, sizeof(struct edge *));
+    if (free_list_e == MY_NULL)
+	G_fatal_error(_("Not enough memory."));
+    for (i = 0; i < n_free_e; i++, e++)
+	free_list_e[i] = e;
+}
+
+
+void free_memory()
+{
+    G_free(sites);
+    G_free(edges);
+    G_free(free_list_e);
+}
+
+struct edge *get_edge()
+{
+    if (n_free_e < 1)
+	G_fatal_error(_("All allocated edges have been used."));
+    return (free_list_e[--n_free_e]);
+}
+
+void free_edge(struct edge *e)
+{
+    free_list_e[n_free_e++] = e;
+}

vector/v.delaunay2/memory.h

@@ -0,0 +1,10 @@
+#ifndef MEMORY_H
+#define MEMORY_H
+void alloc_memory(unsigned int n);
+void free_memory();
+struct edge *get_edge();
+void free_edge(struct edge *e);
+void alloc_sites(unsigned int n);
+void realloc_sites(unsigned int n);
+void alloc_edges(unsigned int n);
+#endif

vector/v.delaunay2/v.delaunay2.html

@@ -0,0 +1,50 @@
+<H2>DESCRIPTION</H2>
+
+<EM>v.delaunay2</EM> uses an existing vector points map (<B>input</B>)
+to create a Delaunay triangulation vector map (<B>output</B>).
+<P>
+
+<BR>
+Delaunay triangulation example:
+<center>
+<img src=v_delaunay_spearfish60_archsites.png border=1><BR>
+<table border=0 width=590>
+<tr><td><center>
+<i>Delaunay Triangulation</i>
+</center></td></tr>
+</table>
+</center>
+
+
+<H2>EXAMPLE</H2>
+
+Commands used with the Spearfish dataset to create the above figure.
+<div class="code"><pre>
+  g.region n=4927250 s=4920000 w=588650 e=605850
+  v.delaunay2 -lr in=archsites out=arch_delaunay
+  d.vect map=arch_delaunay color=0:0:255
+</pre></div>
+
+
+<H2>REFERENCES</H2>
+<EM>Leonid Guibas and Jorge Stolfi, (1985). 
+Primitives for the 
+Manipulation of General Subdivisions and the Computation of
+Voronoi Diagrams, ACM Transactions on Graphics, Vol 4, No. 2, 
+April 1985, Pages 74-123
+</EM>
+
+<H2>SEE ALSO</H2>
+<EM>
+<A HREF="v.voronoi.html">v.voronoi</A>, 
+<A HREF="v.hull.html">v.hull</A>
+</EM>
+
+
+<H2>AUTHORS</H2>
+Martin Pavlovsky, Google Summer of Code 2008, Student<br>
+Paul Kelly, Mentor<br>
+
+
+<p>
+<i>Last changed: $Date: (Mon, 18 Aug 2008) $</i></p>