[GRASS-SVN] r32667 - in grass-addons/vector/voronoi: . demo
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Aug 10 10:00:11 EDT 2008
Author: pkelly
Date: 2008-08-10 10:00:11 -0400 (Sun, 10 Aug 2008)
New Revision: 32667
Added:
grass-addons/vector/voronoi/demo/
grass-addons/vector/voronoi/demo/Makefile
grass-addons/vector/voronoi/demo/README
grass-addons/vector/voronoi/demo/data_types.h
grass-addons/vector/voronoi/demo/edge.c
grass-addons/vector/voronoi/demo/edge.h
grass-addons/vector/voronoi/demo/error.c
grass-addons/vector/voronoi/demo/error.h
grass-addons/vector/voronoi/demo/geom_primitives.h
grass-addons/vector/voronoi/demo/geometry.c
grass-addons/vector/voronoi/demo/geometry.h
grass-addons/vector/voronoi/demo/in_out.c
grass-addons/vector/voronoi/demo/in_out.h
grass-addons/vector/voronoi/demo/main.c
grass-addons/vector/voronoi/demo/memory.c
grass-addons/vector/voronoi/demo/memory.h
Log:
Interim demo code (not yet a GRASS module and still with
bugs to be fixed and optimisations to be made) from Martin
Pavlovsky
Added: grass-addons/vector/voronoi/demo/Makefile
===================================================================
--- grass-addons/vector/voronoi/demo/Makefile (rev 0)
+++ grass-addons/vector/voronoi/demo/Makefile 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,37 @@
+CC = gcc
+CFLAGS = -O0 -g3 $$(pkg-config --cflags --libs cairo)
+LDFLAGS = -lm
+INCS =\
+ data_types.h\
+ edge.h\
+ error.h\
+ geometry.h\
+ geom_primitives.h\
+ in_out.h\
+ memory.h\
+
+SRC =\
+ edge.c\
+ error.c\
+ in_out.c\
+ main.c\
+ memory.c\
+ geometry.c
+
+
+OBJS =\
+ edge.o\
+ error.o\
+ in_out.o\
+ main.o\
+ memory.o\
+ geometry.o
+
+TEXT = Makefile README $(INCS) $(SRC)
+
+guibas-stolfi : $(OBJS)
+ $(CC) -o $@ $(CFLAGS) $(OBJS) $(LDFLAGS)
+
+clean:
+ rm -f guibas-stolfi $(OBJS)
+
Added: grass-addons/vector/voronoi/demo/README
===================================================================
--- grass-addons/vector/voronoi/demo/README (rev 0)
+++ grass-addons/vector/voronoi/demo/README 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,39 @@
+ Delaunay Triangulation Demonstration Code
+ (Guibas-Stolfi algorithm with winged edge data structure)
+ =========================================================
+
+Source code from Martin Pavlovsky, July 2008 as part of interim evaluation
+in the Google Summer of Code project "Reimplementation of v.voronoi and
+v.delaunay modules in the Vector library of GRASS GIS using more efficient
+algorithms".
+
+Description from the author:
+
+Date: Sun, 13 Jul 2008 18:57:46 +0200
+From: "Martin Pavlovsky" <mpavlovsky at gmail.com>
+To: "Paul Kelly" <paul-grass at stjohnspoint.co.uk>
+Subject: UPDATE: SVGAlib replaced with cairo, freezes are no longer present
+
+[...]
+Source code in the attachement is just the implementation of the
+guibas-stolfi algorithm with winged-edge data structure, not an actual GRASS
+module. I will start working on the module itself as soon as I feel that the
+algorithm is robust and stable enough.
+
+Since I understand that it is almost impossible for a human to check if a
+triangulation of lets say 1000 or more sites is correct by just reading text
+output of the edges of the triangulation (e.g. pairs of the form:
+{edge_origin edge_destination}) I decided to use graphic library Cairo
+to visualize the triangulation.
+
+You can compile the sources using makefile included in the directory. Then
+you run the executable guibas-stolfi.
+
+The program generates the Delaunay-triangulation of N randomly generated
+sites(this is just for testing and presentation purposes). After starting
+the program, you are prompted to enter number of sites. After you run the
+program, you can find the visual output in file "dt.png".
+
+You may notice that triangulation for a higher number of sites (>500) are
+sometimes not correct(merging of two halves is not complete occasionaly). I
+am working on correcting this bug.
Added: grass-addons/vector/voronoi/demo/data_types.h
===================================================================
--- grass-addons/vector/voronoi/demo/data_types.h (rev 0)
+++ grass-addons/vector/voronoi/demo/data_types.h 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,27 @@
+#include <stdlib.h>
+
+#ifndef NULL
+#define NULL 0
+#endif
+#define TRUE 1
+#define FALSE 0
+
+typedef enum {left, right} side;
+
+typedef unsigned char boolean;
+
+struct vertex {
+ double x,y;
+ 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;
Added: grass-addons/vector/voronoi/demo/edge.c
===================================================================
--- grass-addons/vector/voronoi/demo/edge.c (rev 0)
+++ grass-addons/vector/voronoi/demo/edge.c 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,117 @@
+#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 == NULL)
+ v1->entry_pt = new_edge;
+ if (v2->entry_pt == NULL)
+ v2->entry_pt = new_edge;
+ return new_edge;
+}
Added: grass-addons/vector/voronoi/demo/edge.h
===================================================================
--- grass-addons/vector/voronoi/demo/edge.h (rev 0)
+++ grass-addons/vector/voronoi/demo/edge.h 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,19 @@
+#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);
+
Added: grass-addons/vector/voronoi/demo/error.c
===================================================================
--- grass-addons/vector/voronoi/demo/error.c (rev 0)
+++ grass-addons/vector/voronoi/demo/error.c 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,7 @@
+#include <stdlib.h>
+#include <stdio.h>
+
+void problem(char *message){
+ fprintf(stderr, message);
+ exit(EXIT_FAILURE);
+}
Added: grass-addons/vector/voronoi/demo/error.h
===================================================================
--- grass-addons/vector/voronoi/demo/error.h (rev 0)
+++ grass-addons/vector/voronoi/demo/error.h 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,2 @@
+
+void problem(char *message);
\ No newline at end of file
Added: grass-addons/vector/voronoi/demo/geom_primitives.h
===================================================================
--- grass-addons/vector/voronoi/demo/geom_primitives.h (rev 0)
+++ grass-addons/vector/voronoi/demo/geom_primitives.h 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,15 @@
+
+#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))
+
+
Added: grass-addons/vector/voronoi/demo/geometry.c
===================================================================
--- grass-addons/vector/voronoi/demo/geometry.c (rev 0)
+++ grass-addons/vector/voronoi/demo/geometry.c 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,255 @@
+#include "data_types.h"
+#include "memory.h"
+#include "geometry.h"
+#include "geom_primitives.h"
+#include "edge.h"
+#include <stddef.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 edge references to be returned */
+ *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)
+ if (CROSS_PRODUCT_3P(o_l, d_l, o_r) > 0.0) {
+ l = PREV(l, d_l);
+ o_l = d_l;
+ d_l = OTHER_VERTEX(l, o_l);
+ } else if (CROSS_PRODUCT_3P(o_r, d_r, o_l) < 0.0) {
+ 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;
+ double u_l_c_o_b, v_l_c_o_b, u_l_c_d_b, v_l_c_d_b;
+ double u_r_c_o_b, v_r_c_o_b, u_r_c_d_b, v_r_c_d_b;
+ /* cross product */
+ double c_p_l_cand, c_p_r_cand;
+ /* dot product */
+ 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;
+ 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 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 for above and "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;
+ if (!above_l_cand && !above_r_cand)
+ break; /* Terminate loop */
+
+ /* Advance 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;
+ }
+ }
+
+ /* Do the same for r_cand symmetrically. */
+ 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;
+ }
+ }
+
+ /*
+ * Now add a cross edge from base to either l_cand or r_cand.
+ * If both are valid choose on the basis of the in_circle test .
+ * Advance base and whichever candidate was chosen.
+ */
+ 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;
+ }
+ }
+}
Added: grass-addons/vector/voronoi/demo/geometry.h
===================================================================
--- grass-addons/vector/voronoi/demo/geometry.h (rev 0)
+++ grass-addons/vector/voronoi/demo/geometry.h 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,3 @@
+
+void divide(struct vertex *sites_sorted[], unsigned int l, unsigned int r,
+ struct edge **l_ccw, struct edge **r_cw);
Added: grass-addons/vector/voronoi/demo/in_out.c
===================================================================
--- grass-addons/vector/voronoi/demo/in_out.c (rev 0)
+++ grass-addons/vector/voronoi/demo/in_out.c 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,137 @@
+#include "data_types.h"
+#include "memory.h"
+#include "edge.h"
+#include "geom_primitives.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <cairo.h>
+#include <math.h>
+
+cairo_surface_t *surface;
+cairo_t *cr;
+
+static void print_edges(unsigned int n);
+static void print_triangles(unsigned int n);
+
+void read_points(unsigned int n){
+ int i;
+ for (i = 0; i < n; i++)
+ if (scanf("%f %f", &sites[i].x, &sites[i].y) != 2)
+ problem("Error reading sites.");
+
+}
+
+void print_results(unsigned int n, char option){
+ /* Print either triangles or edges */
+ if (option == 't')
+ print_triangles(n);
+ else
+ print_edges(n);
+}
+
+/*
+ * Print the ring of edges about each vertex.
+ */
+static void print_edges(unsigned int n){
+ struct edge *e_start, *e;
+ struct vertex *u, *v;
+ unsigned int i;
+
+ for (i = 0; i < n; i++) {
+ u = &sites[i];
+ e_start = e = u->entry_pt;
+ do{
+ v = OTHER_VERTEX(e, u);
+ if (u < v)
+ if (printf("%d %d\n", u - sites, v - sites) == EOF)
+ problem("Error printing results\n");
+ e = NEXT(e, u);
+ } while (!SAME_EDGE(e, e_start));
+ }
+}
+
+/* draw triangulation */
+void draw_sites(int n){
+ int i;
+
+ surface = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, 1280, 800);
+
+ cr = cairo_create(surface);
+/* for (i = 0; i < n; i++) {
+ cairo_arc(cr, sites[i].x, sites[i].y, 1, 0, 2 * M_PI);
+ cairo_stroke (cr);
+ }*/
+}
+
+void draw(unsigned int n){
+ struct edge *e_start, *e;
+ struct vertex *u, *v;
+ unsigned int i;
+
+ draw_sites(n);
+ cairo_set_line_width (cr, 1);
+
+ int x1, y1, x2, y2;
+ for (i = 0; i < n; i++) {
+ u = &sites[i];
+ e_start = e = u->entry_pt;
+ do{
+ v = OTHER_VERTEX(e, u);
+ if (u < v) {
+ x1 = (int) sites[u - sites].x;
+ y1 = (int) sites[u - sites].y;
+ x2 = (int) sites[v - sites].x;
+ y2 = (int) sites[v - sites].y;
+
+ cairo_move_to (cr, x1, y1);
+ cairo_line_to (cr, x2, y2);
+ cairo_stroke (cr);
+ /* printf("e1=%d e2=%d\n", u - sites, v - sites);
+ printf("%3d %3d %3d %3d \n", x1, y1, x2, y2);
+ */
+ }
+ e = NEXT(e, u);
+ } while (!SAME_EDGE(e, e_start));
+ }
+ cairo_destroy (cr);
+ cairo_surface_write_to_png (surface, "dt.png");
+ cairo_surface_destroy (surface);
+}
+
+/*
+ * Print the ring of triangles about each vertex.
+ */
+static void print_triangles(unsigned int n)
+{
+ struct edge *e_start, *e, *next;
+ struct vertex *u, *v, *w;
+ unsigned int i;
+ struct vertex *temp;
+
+ for (i = 0; i < n; i++) {
+ u = &sites[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;
+ }
+ if (printf("%d %d %d\n", u - sites, v - sites, w - sites) == EOF)
+ problem("Error printing results\n");
+ }
+ }
+ /* Next edge around u. */
+ e = NEXT(e, u);
+ } while (!SAME_EDGE(e, e_start));
+ }
+}
Added: grass-addons/vector/voronoi/demo/in_out.h
===================================================================
--- grass-addons/vector/voronoi/demo/in_out.h (rev 0)
+++ grass-addons/vector/voronoi/demo/in_out.h 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,3 @@
+void read_points(unsigned int n);
+void print_results(unsigned int n, char option);
+void draw(unsigned int n);
\ No newline at end of file
Added: grass-addons/vector/voronoi/demo/main.c
===================================================================
--- grass-addons/vector/voronoi/demo/main.c (rev 0)
+++ grass-addons/vector/voronoi/demo/main.c 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,64 @@
+#include "data_types.h"
+#include "memory.h"
+#include "edge.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+int compare(const struct vertex **p1, const struct vertex **p2);
+
+int main(int argc, char *argv[]){
+
+ unsigned int i;
+ unsigned int n;
+ struct edge *l_cw, *r_ccw;
+ struct vertex **sites_sorted;
+
+ printf("Enter number of sites: ");
+ if (scanf("%d", &n) != 1)
+ problem("Number of sites could not be read.");
+ if (n <= 0)
+ problem("Number of sites has to be positive integer.");
+
+ alloc_memory(n);
+ //read_points(n);
+
+ srand((unsigned int)time(0));
+ for (i = 0; i < n; i++){
+ sites[i].x = rand() % 1260 + 10;
+ sites[i].y = rand() % 780 + 10;
+ }
+
+ /* Initialise entry edge vertices. */
+ for (i = 0; i < n; i++)
+ sites[i].entry_pt = NULL;
+ /* Sort. */
+ sites_sorted = (struct vertex **) malloc((unsigned)n*sizeof(struct vertex *));
+ if (sites_sorted == NULL)
+ problem("Not enough memory\n");
+ for (i = 0; i < n; i++)
+ sites_sorted[i] = sites + i;
+
+ qsort(sites_sorted, n, sizeof(struct vertex *), (void *)compare);
+
+// for (i = 0; i < n; i++)
+// printf("(%3d) x=%3.0f y=%3.f\n", i, sites[i].x, sites[i].y);
+
+ /* Triangulate. */
+ divide(sites_sorted, 0, n-1, &l_cw, &r_ccw);
+
+ free((char *)sites_sorted);
+
+ draw(n);
+
+ free_memory();
+
+ exit(EXIT_SUCCESS);
+ return 0; /* To keep lint quiet. */
+}
+
+/* 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);
+}
Added: grass-addons/vector/voronoi/demo/memory.c
===================================================================
--- grass-addons/vector/voronoi/demo/memory.c (rev 0)
+++ grass-addons/vector/voronoi/demo/memory.c 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,45 @@
+#include "data_types.h"
+#include <stdlib.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;
+
+ /* Point storage. */
+ sites = (struct vertex *) calloc(n, sizeof(struct vertex));
+ if (sites == NULL)
+ problem("Not enough memory\n");
+
+ /* Edges. */
+ /* Euler's formula - at most 3n edges on a set of n sites */
+ n_free_e = 3 * n;
+ edges = e = (struct edge *) calloc(n_free_e, sizeof(struct edge));
+ if (edges == NULL)
+ problem("Not enough memory\n");
+ free_list_e = (struct edge **) calloc(n_free_e, sizeof(struct edge *));
+ if (free_list_e == NULL)
+ problem("Not enough memory\n");
+ for (i = 0; i < n_free_e; i++, e++)
+ free_list_e[i] = e;
+}
+
+void free_memory(){
+ free(sites);
+ free(edges);
+ free(free_list_e);
+}
+
+struct edge *get_edge(){
+ if (n_free_e < 1) problem("All allocated edges have been used\n");
+ return (free_list_e[--n_free_e]);
+}
+
+void free_edge(struct edge *e){
+ free_list_e[n_free_e++] = e;
+}
Added: grass-addons/vector/voronoi/demo/memory.h
===================================================================
--- grass-addons/vector/voronoi/demo/memory.h (rev 0)
+++ grass-addons/vector/voronoi/demo/memory.h 2008-08-10 14:00:11 UTC (rev 32667)
@@ -0,0 +1,5 @@
+
+void alloc_memory(unsigned int n);
+void free_memory();
+struct edge *get_edge();
+void free_edge(struct edge *e);
More information about the grass-commit
mailing list