[GRASS-SVN] r34516 - in grass/branches/develbranch_6/vector: .
v.delaunay2
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Nov 26 16:34:54 EST 2008
Author: neteler
Date: 2008-11-26 16:34:54 -0500 (Wed, 26 Nov 2008)
New Revision: 34516
Added:
grass/branches/develbranch_6/vector/v.delaunay2/
grass/branches/develbranch_6/vector/v.delaunay2/Makefile
grass/branches/develbranch_6/vector/v.delaunay2/data_types.h
grass/branches/develbranch_6/vector/v.delaunay2/defs.h
grass/branches/develbranch_6/vector/v.delaunay2/description.html
grass/branches/develbranch_6/vector/v.delaunay2/edge.c
grass/branches/develbranch_6/vector/v.delaunay2/edge.h
grass/branches/develbranch_6/vector/v.delaunay2/geom_primitives.h
grass/branches/develbranch_6/vector/v.delaunay2/geometry.c
grass/branches/develbranch_6/vector/v.delaunay2/geometry.h
grass/branches/develbranch_6/vector/v.delaunay2/in_out.c
grass/branches/develbranch_6/vector/v.delaunay2/in_out.h
grass/branches/develbranch_6/vector/v.delaunay2/main.c
grass/branches/develbranch_6/vector/v.delaunay2/memory.c
grass/branches/develbranch_6/vector/v.delaunay2/memory.h
grass/branches/develbranch_6/vector/v.delaunay2/v_delaunay_spearfish60_archsites.png
Log:
moved here from GRASS Addons
Added: grass/branches/develbranch_6/vector/v.delaunay2/Makefile
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/Makefile (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/Makefile 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/data_types.h
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/data_types.h (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/data_types.h 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/data_types.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/defs.h
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/defs.h (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/defs.h 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/defs.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/description.html
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/description.html (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/description.html 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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>
Added: grass/branches/develbranch_6/vector/v.delaunay2/edge.c
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/edge.c (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/edge.c 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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;
+}
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/edge.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/edge.h
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/edge.h (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/edge.h 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/edge.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/geom_primitives.h
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/geom_primitives.h (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/geom_primitives.h 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/geom_primitives.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/geometry.c
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/geometry.c (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/geometry.c 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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;
+ }
+ }
+}
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/geometry.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/geometry.h
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/geometry.h (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/geometry.h 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/geometry.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/in_out.c
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/in_out.c (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/in_out.c 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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;
+}
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/in_out.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/in_out.h
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/in_out.h (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/in_out.h 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/in_out.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/main.c
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/main.c (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/main.c 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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);
+}
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/memory.c
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/memory.c (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/memory.c 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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;
+}
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/memory.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/memory.h
===================================================================
--- grass/branches/develbranch_6/vector/v.delaunay2/memory.h (rev 0)
+++ grass/branches/develbranch_6/vector/v.delaunay2/memory.h 2008-11-26 21:34:54 UTC (rev 34516)
@@ -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
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/memory.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/branches/develbranch_6/vector/v.delaunay2/v_delaunay_spearfish60_archsites.png
===================================================================
(Binary files differ)
Property changes on: grass/branches/develbranch_6/vector/v.delaunay2/v_delaunay_spearfish60_archsites.png
___________________________________________________________________
Name: svn:mime-type
+ image/png
More information about the grass-commit
mailing list