v.voronoi: abuse it for skeletons

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

vector/v.voronoi/Makefile

@@ -8,7 +8,7 @@ SWEEP_OBJS = sw_edgelist.o sw_geometry.o sw_heap.o sw_main.o sw_memory.o sw_outp
 
 v_voronoi_OBJS = vo_main.o vo_extend.o vo_write.o $(SWEEP_OBJS)
 
-LIBES = $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB)
+LIBES = $(VECTORLIB) $(DIG2LIB) $(DBMILIB) $(GISLIB) $(MATHLIB)
 DEPENDENCIES = $(VECTORDEP) $(DBMIDEP) $(GISDEP)
 EXTRA_INC = $(VECT_INC)
 EXTRA_CFLAGS = $(VECT_CFLAGS)

vector/v.voronoi/clean_topo.c

@@ -0,0 +1,149 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include "sw_defs.h"
+#include "defs.h"
+
+int clean_topo(void)
+{
+    int type, line, nlines;
+    int area, nareas;
+    int verbose;
+    int err_boundaries, err_centr_out, err_centr_dupl, err_nocentr;
+    struct bound_box box;
+    double snap_thresh;
+
+    /* cleaning part 1: count errors */
+    G_message(_("Searching for topology errors..."));
+    verbose = G_verbose();
+    G_set_verbose(0);
+    Vect_build_partial(&Out, GV_BUILD_CENTROIDS);
+    G_set_verbose(verbose);
+    err_boundaries = err_centr_out = err_centr_dupl = err_nocentr = 0;
+    nlines = Vect_get_num_lines(&Out);
+    for (line = 1; line <= nlines; line++) {
+
+	if (!Vect_line_alive(&Out, line))
+	    continue;
+
+	type = Vect_get_line_type(&Out, line);
+	if (type == GV_BOUNDARY) {
+	    int left, right;
+
+	    Vect_get_line_areas(&Out, line, &left, &right);
+
+	    if (left == 0 || right == 0) {
+		G_debug(3, "line = %d left = %d right = %d", line, 
+			left, right);
+		err_boundaries++;
+	    }
+	}
+	if (type == GV_CENTROID) {
+	    area = Vect_get_centroid_area(&Out, line);
+	    if (area == 0)
+		err_centr_out++;
+	    else if (area < 0)
+		err_centr_dupl++;
+	}
+    }
+
+    err_nocentr = 0;
+    nareas = Vect_get_num_areas(&Out);
+    for (area = 1; area <= nareas; area++) {
+	if (!Vect_area_alive(&Out, area))
+	    continue;
+	line = Vect_get_area_centroid(&Out, area);
+	if (line == 0)
+	    err_nocentr++;
+    }
+
+    /* cleaning part 2: snap */
+    /* TODO: adjust snapping treshold to ULP */
+    Vect_get_map_box(&Out, &box);
+    snap_thresh = fabs(box.W);
+    if (snap_thresh < fabs(box.E))
+	snap_thresh = fabs(box.E);
+    if (snap_thresh < fabs(box.N))
+	snap_thresh = fabs(box.N);
+    if (snap_thresh < fabs(box.S))
+	snap_thresh = fabs(box.S);
+    snap_thresh = d_ulp(snap_thresh);
+    
+    if (err_nocentr || err_centr_dupl || err_centr_out) {
+	int nmod;
+
+	G_important_message(_("Cleaning output topology"));
+	Vect_snap_lines(&Out, GV_BOUNDARY, snap_thresh, NULL);
+	do {
+	    Vect_break_lines(&Out, GV_BOUNDARY, NULL);
+	    Vect_remove_duplicates(&Out, GV_BOUNDARY, NULL);
+	    nmod =
+		Vect_clean_small_angles_at_nodes(&Out, GV_BOUNDARY, NULL);
+	} while (nmod > 0);
+
+	G_message(_("Removing dangles..."));
+	Vect_remove_dangles(&Out, GV_BOUNDARY, -1.0, NULL);
+	G_message(_("Removing bridges..."));
+	Vect_remove_bridges(&Out, NULL, NULL, NULL);
+
+	err_boundaries = 0;
+	nlines = Vect_get_num_lines(&Out);
+	for (line = 1; line <= nlines; line++) {
+
+	    if (!Vect_line_alive(&Out, line))
+		continue;
+
+	    type = Vect_get_line_type(&Out, line);
+	    if (type == GV_BOUNDARY) {
+		int left, right;
+
+		Vect_get_line_areas(&Out, line, &left, &right);
+
+		if (left == 0 || right == 0) {
+		    G_debug(3, "line = %d left = %d right = %d", line, 
+			    left, right);
+		    err_boundaries++;
+		}
+	    }
+	}
+    }
+    /* cleaning part 3: remove remaining incorrect boundaries */
+    if (err_boundaries) {
+	G_important_message(_("Removing incorrect boundaries from output"));
+	nlines = Vect_get_num_lines(&Out);
+	for (line = 1; line <= nlines; line++) {
+
+	    if (!Vect_line_alive(&Out, line))
+		continue;
+
+	    type = Vect_get_line_type(&Out, line);
+	    if (type == GV_BOUNDARY) {
+		int left, right;
+
+		Vect_get_line_areas(&Out, line, &left, &right);
+
+		/* &&, not ||, no typo */
+		if (left == 0 && right == 0) {
+		    G_debug(3, "line = %d left = %d right = %d", line, 
+			    left, right);
+		    Vect_delete_line(&Out, line);
+		}
+	    }
+	}
+    }
+
+    /* this is slow:
+    if (in_area) {
+	if (Type == GV_LINE)
+	    G_message(_("Merging lines ..."));
+	else
+	    G_message(_("Merging boundaries ..."));
+	Vect_merge_lines(&Out, Type, NULL, NULL);
+    }
+    */
+
+    return 1;
+}

vector/v.voronoi/defs.h

@@ -5,4 +5,5 @@ extern struct Map_info In, Out;
 extern int Type;
 extern int Field;
 extern int in_area;
+extern int skeleton;
 extern double segf;

vector/v.voronoi/main.c

@@ -7,9 +7,10 @@
  *               Andrea Aime <aaime libero.it>
  *               Radim Blazek <radim.blazek gmail.com> (GRASS 5.1 v.voronoi) 
  *               Glynn Clements <glynn gclements.plus.com>,  
- *               Markus Neteler <neteler itc.it>
+ *               Markus Neteler <neteler itc.it>,
+ *               Markus Metz
  * PURPOSE:      produce a Voronoi diagram using vector points
- * COPYRIGHT:    (C) 1993-2006, 2001 by the GRASS Development Team
+ * COPYRIGHT:    (C) 1993-2006, 2001, 2013 by the GRASS Development Team
  *
  *               This program is free software under the GNU General Public
  *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -51,6 +52,8 @@ typedef struct
     double x, y;
 } COOR;
 
+
+
 int cmp(void *a, void *b)
 {
     COOR *ca = (COOR *) a;
@@ -99,22 +102,20 @@ int main(int argc, char **argv)
     int i;
     int **cats, *ncats, nfields, *fields;
     struct {
-	struct Flag *line, *table, *area;
-    } flag;
-    struct {
-	struct Option *in, *out, *field, *smooth;
+	struct Option *in, *out, *field, *smooth, *thin;
     } opt;
+    struct {
+	struct Flag *line, *table, *area, *skeleton;
+    } flag;
     struct GModule *module;
     struct line_pnts *Points;
     struct line_cats *Cats;
-    struct bound_box box;
     int node, nnodes;
     COOR *coor;
     int verbose;
     int ncoor, acoor;
-    int line, nlines, type, ctype, area, nareas;
-    int err_boundaries, err_centr_out, err_centr_dupl, err_nocentr;
-    double snap_thresh;
+    int line, nlines, type, ctype;
+    double thresh;
 
     G_gisinit(argv[0]);
 
@@ -137,13 +138,27 @@ int main(int argc, char **argv)
     opt.smooth->key = "smoothness";
     opt.smooth->answer = "0.25";
     opt.smooth->label = _("Factor for output smoothness");
-    opt.smooth->description = _("Applies to input areas only. Smaller values produce smoother output but can cause numerical instability.");
+    opt.smooth->description = _("Applies to input areas only. "
+                                "Smaller values produce smoother output but can cause numerical instability.");
+
+    opt.thin = G_define_option();
+    opt.thin->type = TYPE_DOUBLE;
+    opt.thin->key = "thin";
+    opt.thin->answer = "-1";
+    opt.thin->label = _("Maximum dangle length of skeletons");
+    opt.thin->description = _("Applies only to skeleton extraction. "
+                                "Default = -1 will extract the center line.");
 
     flag.area = G_define_flag();
     flag.area->key = 'a';
     flag.area->description =
 	_("Create Voronoi diagram for input areas");
 
+    flag.skeleton = G_define_flag();
+    flag.skeleton->key = 's';
+    flag.skeleton->description =
+	_("Extract skeletons for input areas");
+
     flag.line = G_define_flag();
     flag.line->key = 'l';
     flag.line->description =
@@ -165,6 +180,11 @@ int main(int argc, char **argv)
 	segf = 0.25;
 	G_warning(_("Option '%s' is too small, set to %g"), opt.smooth->key, segf);
     }
+    thresh = atof(opt.thin->answer);
+    
+    skeleton = flag.skeleton->answer;
+    if (skeleton)
+	Type = GV_LINE;
 
     Points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
@@ -184,7 +204,7 @@ int main(int argc, char **argv)
     freeinit(&sfl, sizeof(struct Site));
 
     G_message(_("Reading features..."));
-    if (in_area)
+    if (in_area || skeleton)
 	readbounds();
     else
 	readsites();
@@ -200,72 +220,79 @@ int main(int argc, char **argv)
     plot = 0;
     debug = 0;
 
-    G_message(_("Processing Voronoi triangulation..."));
+    G_message(_("Voronoi triangulation for %d points..."), nsites);
     voronoi(nextone);
+    G_message(_("Writing edges..."));
+    vo_write();
 
-    /* Close free ends by current region */
     verbose = G_verbose();
     G_set_verbose(0);
     Vect_build_partial(&Out, GV_BUILD_BASE);
     G_set_verbose(verbose);
 
-    ncoor = 0;
-    acoor = 100;
-    coor = (COOR *) G_malloc(sizeof(COOR) * acoor);
-
-    nnodes = Vect_get_num_nodes(&Out);
-    for (node = 1; node <= nnodes; node++) {
-	double x, y;
-
-	if (Vect_get_node_n_lines(&Out, node) < 2) {	/* add coordinates */
-	    Vect_get_node_coor(&Out, node, &x, &y, NULL);
+    if (skeleton) {
+	G_message(_("Thin skeletons ..."));
+	thin_skeleton(thresh);
+    }
+    else {
+	/* Close free ends by current region */
+	ncoor = 0;
+	acoor = 100;
+	coor = (COOR *) G_malloc(sizeof(COOR) * acoor);
+
+	nnodes = Vect_get_num_nodes(&Out);
+	for (node = 1; node <= nnodes; node++) {
+	    double x, y;
+
+	    if (Vect_get_node_n_lines(&Out, node) < 2) {	/* add coordinates */
+		Vect_get_node_coor(&Out, node, &x, &y, NULL);
+
+		if (ncoor == acoor - 5) {	/* always space for 5 region corners */
+		    acoor += 100;
+		    coor = (COOR *) G_realloc(coor, sizeof(COOR) * acoor);
+		}
 
-	    if (ncoor == acoor - 5) {	/* always space for 5 region corners */
-		acoor += 100;
-		coor = (COOR *) G_realloc(coor, sizeof(COOR) * acoor);
+		coor[ncoor].x = x;
+		coor[ncoor].y = y;
+		ncoor++;
 	    }
+	}
 
-	    coor[ncoor].x = x;
-	    coor[ncoor].y = y;
-	    ncoor++;
+	/* Add region corners */
+	coor[ncoor].x = Box.W;
+	coor[ncoor].y = Box.S;
+	ncoor++;
+	coor[ncoor].x = Box.E;
+	coor[ncoor].y = Box.S;
+	ncoor++;
+	coor[ncoor].x = Box.E;
+	coor[ncoor].y = Box.N;
+	ncoor++;
+	coor[ncoor].x = Box.W;
+	coor[ncoor].y = Box.N;
+	ncoor++;
+
+	/* Sort */
+	qsort(coor, ncoor, sizeof(COOR), (void *)cmp);
+
+	/* add last (first corner) */
+	coor[ncoor].x = Box.W;
+	coor[ncoor].y = Box.S;
+	ncoor++;
+
+	for (i = 1; i < ncoor; i++) {
+	    if (coor[i].x == coor[i - 1].x && coor[i].y == coor[i - 1].y)
+		continue;		/* duplicate */
+
+	    Vect_reset_line(Points);
+	    Vect_append_point(Points, coor[i].x, coor[i].y, 0.0);
+	    Vect_append_point(Points, coor[i - 1].x, coor[i - 1].y, 0.0);
+	    Vect_write_line(&Out, Type, Points, Cats);
 	}
-    }
 
-    /* Add region corners */
-    coor[ncoor].x = Box.W;
-    coor[ncoor].y = Box.S;
-    ncoor++;
-    coor[ncoor].x = Box.E;
-    coor[ncoor].y = Box.S;
-    ncoor++;
-    coor[ncoor].x = Box.E;
-    coor[ncoor].y = Box.N;
-    ncoor++;
-    coor[ncoor].x = Box.W;
-    coor[ncoor].y = Box.N;
-    ncoor++;
-
-    /* Sort */
-    qsort(coor, ncoor, sizeof(COOR), (void *)cmp);
-
-    /* add last (first corner) */
-    coor[ncoor].x = Box.W;
-    coor[ncoor].y = Box.S;
-    ncoor++;
-
-    for (i = 1; i < ncoor; i++) {
-	if (coor[i].x == coor[i - 1].x && coor[i].y == coor[i - 1].y)
-	    continue;		/* duplicate */
-
-	Vect_reset_line(Points);
-	Vect_append_point(Points, coor[i].x, coor[i].y, 0.0);
-	Vect_append_point(Points, coor[i - 1].x, coor[i - 1].y, 0.0);
-	Vect_write_line(&Out, Type, Points, Cats);
+	G_free(coor);
     }
 
-    G_free(coor);
-
-    /* Copy input points as centroids */
     nfields = Vect_cidx_get_num_fields(&In);
     cats = (int **)G_malloc(nfields * sizeof(int *));
     ncats = (int *)G_malloc(nfields * sizeof(int));
@@ -278,7 +305,8 @@ int main(int argc, char **argv)
 	fields[i] = Vect_cidx_get_field_number(&In, i);
     }
 
-    if (flag.line->answer)
+    /* Copy input points */
+    if (Type == GV_LINE)
 	ctype = GV_POINT;
     else
 	ctype = GV_CENTROID;
@@ -298,8 +326,8 @@ int main(int argc, char **argv)
 	if (!Vect_point_in_box(Points->x[0], Points->y[0], 0.0, &Box))
 	    continue;
 
-	Vect_write_line(&Out, ctype, Points, Cats);
-
+	if (!skeleton)
+	    Vect_write_line(&Out, ctype, Points, Cats);
 
 	for (i = 0; i < Cats->n_cats; i++) {
 	    int f, j;
@@ -382,136 +410,10 @@ int main(int argc, char **argv)
 	}
     }
 
-
     Vect_close(&In);
-
-    /* cleaning part 1: count errors */
-    G_message(_("Searching for topology errors..."));
-    G_set_verbose(0);
-    Vect_build_partial(&Out, GV_BUILD_CENTROIDS);
-    G_set_verbose(verbose);
-    err_boundaries = err_centr_out = err_centr_dupl = err_nocentr = 0;
-    nlines = Vect_get_num_lines(&Out);
-    for (line = 1; line <= nlines; line++) {
-
-	if (!Vect_line_alive(&Out, line))
-	    continue;
-
-	type = Vect_get_line_type(&Out, line);
-	if (type == GV_BOUNDARY) {
-	    int left, right;
-
-	    Vect_get_line_areas(&Out, line, &left, &right);
-
-	    if (left == 0 || right == 0) {
-		G_debug(3, "line = %d left = %d right = %d", line, 
-			left, right);
-		err_boundaries++;
-	    }
-	}
-	if (type == GV_CENTROID) {
-	    area = Vect_get_centroid_area(&Out, line);
-	    if (area == 0)
-		err_centr_out++;
-	    else if (area < 0)
-		err_centr_dupl++;
-	}
-    }
-
-    err_nocentr = 0;
-    nareas = Vect_get_num_areas(&Out);
-    for (area = 1; area <= nareas; area++) {
-	if (!Vect_area_alive(&Out, area))
-	    continue;
-	line = Vect_get_area_centroid(&Out, area);
-	if (line == 0)
-	    err_nocentr++;
-    }
-
-    /* cleaning part 2: snap */
-    /* TODO: adjust snapping treshold to ULP */
-    Vect_get_map_box(&Out, &box);
-    snap_thresh = fabs(box.W);
-    if (snap_thresh < fabs(box.E))
-	snap_thresh = fabs(box.E);
-    if (snap_thresh < fabs(box.N))
-	snap_thresh = fabs(box.N);
-    if (snap_thresh < fabs(box.S))
-	snap_thresh = fabs(box.S);
-    snap_thresh = d_ulp(snap_thresh);
     
-    if (err_nocentr || err_centr_dupl || err_centr_out) {
-	int nmod;
-
-	G_important_message(_("Cleaning output topology"));
-	Vect_snap_lines(&Out, GV_BOUNDARY, snap_thresh, NULL);
-	do {
-	    Vect_break_lines(&Out, GV_BOUNDARY, NULL);
-	    Vect_remove_duplicates(&Out, GV_BOUNDARY, NULL);
-	    nmod =
-		Vect_clean_small_angles_at_nodes(&Out, GV_BOUNDARY, NULL);
-	} while (nmod > 0);
-
-	G_message(_("Removing dangles..."));
-	Vect_remove_dangles(&Out, GV_BOUNDARY, -1.0, NULL);
-	G_message(_("Removing bridges..."));
-	Vect_remove_bridges(&Out, NULL, NULL, NULL);
-
-	err_boundaries = 0;
-	nlines = Vect_get_num_lines(&Out);
-	for (line = 1; line <= nlines; line++) {
-
-	    if (!Vect_line_alive(&Out, line))
-		continue;
-
-	    type = Vect_get_line_type(&Out, line);
-	    if (type == GV_BOUNDARY) {
-		int left, right;
-
-		Vect_get_line_areas(&Out, line, &left, &right);
-
-		if (left == 0 || right == 0) {
-		    G_debug(3, "line = %d left = %d right = %d", line, 
-			    left, right);
-		    err_boundaries++;
-		}
-	    }
-	}
-    }
-    /* cleaning part 3: remove remaining incorrect boundaries */
-    if (err_boundaries) {
-	G_important_message(_("Removing incorrect boundaries from output"));
-	nlines = Vect_get_num_lines(&Out);
-	for (line = 1; line <= nlines; line++) {
-
-	    if (!Vect_line_alive(&Out, line))
-		continue;
-
-	    type = Vect_get_line_type(&Out, line);
-	    if (type == GV_BOUNDARY) {
-		int left, right;
-
-		Vect_get_line_areas(&Out, line, &left, &right);
-
-		/* &&, not ||, no typo */
-		if (left == 0 && right == 0) {
-		    G_debug(3, "line = %d left = %d right = %d", line, 
-			    left, right);
-		    Vect_delete_line(&Out, line);
-		}
-	    }
-	}
-    }
-
-    /* this is slow:
-    if (in_area) {
-	if (Type == GV_LINE)
-	    G_message(_("Merging lines ..."));
-	else
-	    G_message(_("Merging boundaries ..."));
-	Vect_merge_lines(&Out, Type, NULL, NULL);
-    }
-    */
+    if (Type == GV_BOUNDARY)
+	clean_topo();
 
     /* build clean topology */
     Vect_build_partial(&Out, GV_BUILD_NONE);
@@ -519,5 +421,6 @@ int main(int argc, char **argv)
     Vect_close(&Out);
 
     G_done_msg(" ");
+
     exit(EXIT_SUCCESS);
 }

vector/v.voronoi/skeleton.c

@@ -0,0 +1,472 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include "defs.h"
+
+static int next_dist(int line, int side, double mf)
+{
+    double d, dist, nextdist, totaldist;
+    int nextline, nlines;
+    int node;
+    static struct line_pnts *Points = NULL;
+
+    G_debug(3, "next_dist()");
+
+    if (!Points)
+	Points = Vect_new_line_struct();
+
+    Vect_read_line(&Out, Points, NULL, abs(line));
+    dist = Vect_line_length(Points);
+
+    if (line < 0)
+	Vect_get_line_nodes(&Out, -line, &node, NULL);
+    else
+	Vect_get_line_nodes(&Out, line, NULL, &node);
+
+    nlines = Vect_get_node_n_lines(&Out, node);
+
+    if (nlines == 1)
+	return 1;
+
+    nextdist = totaldist = 0.;
+    while (nlines > 1) {
+	nextline = dig_angle_next_line(&(Out.plus), -line, side, GV_LINE, NULL);
+	Vect_read_line(&Out, Points, NULL, abs(nextline));
+	d = Vect_line_length(Points);
+	nextdist += d;
+	totaldist += d;
+
+	if (nextline < 0)
+	    Vect_get_line_nodes(&Out, -nextline, &node, NULL);
+	else
+	    Vect_get_line_nodes(&Out, nextline, NULL, &node);
+
+	nlines = Vect_get_node_n_lines(&Out, node);
+	if (nlines > 2)
+	    nextdist = 0.;
+	line = nextline;
+    }
+
+    if (totaldist > nextdist && dist > nextdist)
+	return (totaldist < mf * dist);
+
+    return (dist > nextdist);
+}
+
+static int loop_test(int line, int node, struct line_pnts *Points,
+                     double l, double mf)
+{
+    int line1, line2, nextline;
+    int n1, n2, nspikes, nout;
+    double l1, minl, maxl;
+    
+    if (Vect_get_node_n_lines(&Out, node) != 3)
+	return 0;
+	
+    line1 = line2 = 0;
+    line1 = Vect_get_node_line(&Out, node, 0);
+    if (abs(line1) == abs(line))
+	line1 = 0;
+    line2 = Vect_get_node_line(&Out, node, 1);
+    if (abs(line2) == abs(line))
+	line2 = 0;
+    if (line1 == 0) {
+	line1 = line2;
+	line2 = 0;
+    }
+    if (line2 == 0)
+	line2 = Vect_get_node_line(&Out, node, 2);
+    
+    if (abs(line1) == abs(line2))
+	return 1;
+
+    nextline = dig_angle_next_line(&(Out.plus), -line, GV_LEFT, GV_LINE, NULL);
+    line1 = nextline;
+    
+    nspikes = 1;
+    nout = 0;
+    minl = 0;
+    maxl = 0;
+    do {
+	if (nextline < 0)
+	    Vect_get_line_nodes(&Out, -nextline, &n1, NULL);
+	else
+	    Vect_get_line_nodes(&Out, nextline, NULL, &n1);
+
+	if (Vect_get_node_n_lines(&Out, n1) == 1)
+	    return 0;
+
+	if (n1 != node && Vect_get_node_n_lines(&Out, n1) > 2) {
+	    nspikes++;
+	    line2 = dig_angle_next_line(&(Out.plus), -nextline, GV_LEFT, GV_LINE, NULL);
+
+	    if (line2 < 0)
+		Vect_get_line_nodes(&Out, -line2, &n2, NULL);
+	    else
+		Vect_get_line_nodes(&Out, line2, NULL, &n2);
+
+	    if (Vect_get_node_n_lines(&Out, n2) == 1) {
+		Vect_read_line(&Out, Points, NULL, abs(line2));
+		l1 = Vect_line_length(Points);
+		if (minl == 0 || minl > l1) {
+		    minl = l1;
+		}
+		if (maxl < l1) {
+		    maxl = l1;
+		}
+	    }
+	    else
+		nout++;
+	}
+	
+	nextline = dig_angle_next_line(&(Out.plus), -nextline, GV_RIGHT, GV_LINE, NULL);
+	
+    } while (abs(nextline) != abs(line1));
+
+    if (minl == 0)
+	minl = l;
+    if (maxl == 0)
+	maxl = mf * l;
+
+    nspikes -= nout;
+
+    if (mf > 1) {
+	return (nspikes < 3 || (nspikes == 3 && (l > minl || mf * l > maxl)));
+
+	if (l > minl)
+	    return 1;
+	if (nspikes == 3 && l < minl)
+	    return (mf * l > maxl);
+	else
+	    return (nspikes < 3);
+    }
+
+    return (nspikes < 3 || (nspikes > 2 && l > minl));
+}
+
+static int break_loop(int line, int node, struct line_pnts *Points)
+{
+    int line1, line2, firstline, nextline;
+    int n1;
+    double l1, l2;
+    
+    if (Vect_get_node_n_lines(&Out, node) != 3)
+	return 0;
+	
+    line1 = line2 = 0;
+    line1 = Vect_get_node_line(&Out, node, 0);
+    if (abs(line1) == abs(line))
+	line1 = 0;
+    line2 = Vect_get_node_line(&Out, node, 1);
+    if (abs(line2) == abs(line))
+	line2 = 0;
+    if (line1 == 0) {
+	line1 = line2;
+	line2 = 0;
+    }
+    if (line2 == 0)
+	line2 = Vect_get_node_line(&Out, node, 2);
+    
+    if (abs(line1) == abs(line2))
+	return 1;
+
+    nextline = dig_angle_next_line(&(Out.plus), -line, GV_LEFT, GV_LINE, NULL);
+    firstline = nextline;
+
+    do {
+	if (nextline < 0)
+	    Vect_get_line_nodes(&Out, -nextline, &n1, NULL);
+	else
+	    Vect_get_line_nodes(&Out, nextline, NULL, &n1);
+
+	if (Vect_get_node_n_lines(&Out, n1) == 1)
+	    return 0;
+
+	nextline = dig_angle_next_line(&(Out.plus), -nextline, GV_RIGHT, GV_LINE, NULL);
+	
+    } while (abs(nextline) != abs(firstline));
+
+    if (abs(nextline) != abs(firstline)) {
+	G_warning("no loop ???");
+	return 0;
+    }
+
+    Vect_read_line(&Out, Points, NULL, abs(line1));
+    l1 = Vect_line_length(Points);
+
+    Vect_read_line(&Out, Points, NULL, abs(line2));
+    l2 = Vect_line_length(Points);
+
+    if (l1 > l2)
+	Vect_delete_line(&Out, abs(line1));
+    else
+	Vect_delete_line(&Out, abs(line2));
+
+    Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
+
+    return 1;
+}
+
+static int length_test(int line, int node, struct line_pnts *Points,
+                       double l, double mf)
+{
+    int line1, line2, nlines;
+    int n1;
+    double l1, l2;
+    
+    if (Vect_get_node_n_lines(&Out, node) != 3)
+	return 0;
+
+    line1 = Vect_get_node_line(&Out, node, 0);
+    line2 = 0;
+    if (abs(line1) == abs(line))
+	line1 = 0;
+    line2 = Vect_get_node_line(&Out, node, 1);
+    if (abs(line2) == abs(line))
+	line2 = 0;
+    if (line1 == 0) {
+	line1 = line2;
+	line2 = 0;
+    }
+    if (line2 == 0)
+	line2 = Vect_get_node_line(&Out, node, 2);
+
+    if (line1 < 0)
+	Vect_get_line_nodes(&Out, -line1, &n1, NULL);
+    else
+	Vect_get_line_nodes(&Out, line1, NULL, &n1);
+
+    nlines = Vect_get_node_n_lines(&Out, n1);
+    if (nlines > 1)
+	return 0;
+
+    if (line2 < 0)
+	Vect_get_line_nodes(&Out, -line2, &n1, NULL);
+    else
+	Vect_get_line_nodes(&Out, line2, NULL, &n1);
+
+    nlines = Vect_get_node_n_lines(&Out, n1);
+    if (nlines > 1)
+	return 0;
+
+    Vect_read_line(&Out, Points, NULL, abs(line1));
+    l1 = Vect_line_length(Points);
+
+    Vect_read_line(&Out, Points, NULL, abs(line2));
+    l2 = Vect_line_length(Points);
+
+    if (l1 > mf * l2) {
+	if (mf * l < l1 && l < l2)
+	    return 0;
+    }
+    if (l2 > mf * l1) {
+	if (mf * l < l2 && l < l1)
+	    return 0;
+    }
+
+    return (mf * l > l1 || mf * l > l2);
+}
+
+/* thin the skeletons */
+int thin_skeleton(double thresh)
+{
+    int i;
+    int node, n1, n2;
+    struct line_pnts *Points;
+    struct ilist *list;
+    double l, minl, maxl;
+    int nlines, line, minline, line2;
+    int counter = 1;
+    double morphof = 1.6180339887;
+
+    Points = Vect_new_line_struct();
+    list = Vect_new_list();
+
+    if (thresh < 0)
+	morphof = 1;
+
+    Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
+
+    while (TRUE) {
+	G_verbose_message(_("Pass %d"), counter++);
+	Vect_reset_list(list);
+
+	for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
+	    if (!Vect_node_alive(&Out, node))
+		continue;
+	    nlines = Vect_get_node_n_lines(&Out, node);
+	    
+	    if (nlines > 1)
+		continue;
+
+	    line = Vect_get_node_line(&Out, node, 0);
+	    if (line < 0)
+		Vect_get_line_nodes(&Out, -line, &n1, NULL);
+	    else
+		Vect_get_line_nodes(&Out, line, NULL, &n1);
+
+	    nlines = Vect_get_node_n_lines(&Out, n1);
+	    
+	    if (nlines < 3)
+		continue;
+
+	    Vect_read_line(&Out, Points, NULL, abs(line));
+	    minline = line;
+	    minl = Vect_line_length(Points);
+
+	    if (nlines == 3) {
+		if (loop_test(line, n1, Points, minl, morphof))
+		    continue;
+		if (length_test(line, n1, Points, minl, morphof))
+		    continue;
+	    }
+
+	    maxl = 0;
+	    for (i = 0; i < nlines; i++) {
+		line2 = Vect_get_node_line(&Out, n1, i);
+		if (abs(line2) == abs(minline) || abs(line2) == abs(line))
+		    continue;
+		if (line2 < 0)
+		    Vect_get_line_nodes(&Out, -line2, &n2, NULL);
+		else
+		    Vect_get_line_nodes(&Out, line2, NULL, &n2);
+		if (Vect_get_node_n_lines(&Out, n2) > 1)
+		    continue;
+		Vect_read_line(&Out, Points, NULL, abs(line2));
+		l = Vect_line_length(Points);
+		if (minl > l) {
+		    minl = l;
+		    minline = line2;
+		}
+		if (maxl == 0 || maxl < l) {
+		    maxl = l;
+		}
+	    }
+	    if (thresh < 0) {
+		G_ilist_add(list, minline);
+	    }
+	    else  {
+		if (minl < thresh)
+		    G_ilist_add(list, minline);
+	    }
+	}
+	if (list->n_values == 0)
+	    break;
+	nlines = 0;
+	for (i = 0; i < list->n_values; i++) {
+	    line = list->value[i];
+	    if (Vect_line_alive(&Out, abs(line))) {
+		if (next_dist(line, GV_RIGHT, morphof))
+		    continue;
+		if (next_dist(line, GV_LEFT, morphof))
+		    continue;
+		Vect_delete_line(&Out, abs(line));
+		nlines++;
+	    }
+	}
+	if (nlines == 0)
+	    break;
+	else
+	    Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
+    }
+
+    if (thresh >= 0)
+	return 0;
+
+    for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
+	if (!Vect_node_alive(&Out, node))
+	    continue;
+	nlines = Vect_get_node_n_lines(&Out, node);
+	
+	if (nlines > 1)
+	    continue;
+
+	line = Vect_get_node_line(&Out, node, 0);
+	if (line < 0)
+	    Vect_get_line_nodes(&Out, -line, &n1, NULL);
+	else
+	    Vect_get_line_nodes(&Out, line, NULL, &n1);
+
+	nlines = Vect_get_node_n_lines(&Out, n1);
+	
+	if (nlines == 3) {
+	    break_loop(line, n1, Points);
+	}
+    }
+
+    while (TRUE) {
+	G_verbose_message(_("Pass %d"), counter++);
+	Vect_reset_list(list);
+
+	for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
+	    if (!Vect_node_alive(&Out, node))
+		continue;
+	    nlines = Vect_get_node_n_lines(&Out, node);
+	    
+	    if (nlines > 1)
+		continue;
+
+	    line = Vect_get_node_line(&Out, node, 0);
+	    if (line < 0)
+		Vect_get_line_nodes(&Out, -line, &n1, NULL);
+	    else
+		Vect_get_line_nodes(&Out, line, NULL, &n1);
+
+	    nlines = Vect_get_node_n_lines(&Out, n1);
+	    
+	    if (nlines < 3)
+		continue;
+
+	    Vect_read_line(&Out, Points, NULL, abs(line));
+	    minline = line;
+	    minl = Vect_line_length(Points);
+
+	    if (nlines == 3) {
+		if (break_loop(line, n1, Points))
+		    continue;
+	    }
+
+	    for (i = 0; i < nlines; i++) {
+		line2 = Vect_get_node_line(&Out, n1, i);
+		if (abs(line2) == abs(minline) || abs(line2) == abs(line))
+		    continue;
+		if (line2 < 0)
+		    Vect_get_line_nodes(&Out, -line2, &n2, NULL);
+		else
+		    Vect_get_line_nodes(&Out, line2, NULL, &n2);
+		if (Vect_get_node_n_lines(&Out, n2) > 1)
+		    continue;
+		Vect_read_line(&Out, Points, NULL, abs(line2));
+		l = Vect_line_length(Points);
+		if (minl > l) {
+		    minl = l;
+		    minline = line2;
+		}
+	    }
+	    if (thresh < 0 || minl < thresh)
+		G_ilist_add(list, minline);
+	}
+	if (list->n_values == 0)
+	    break;
+	nlines = 0;
+	for (i = 0; i < list->n_values; i++) {
+	    line = list->value[i];
+	    if (Vect_line_alive(&Out, abs(line))) {
+		if (next_dist(line, GV_RIGHT, morphof))
+		    continue;
+		if (next_dist(line, GV_LEFT, morphof))
+		    continue;
+		Vect_delete_line(&Out, abs(line));
+		nlines++;
+	    }
+	}
+	if (nlines == 0)
+	    break;
+	else
+	    Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
+    }
+
+    return 0;
+}

vector/v.voronoi/sw_defs.h

@@ -66,6 +66,12 @@ extern struct Halfedge *PQhash;
 extern int PQcount;
 extern int PQmin;
 
+/* clean_topo.c */
+int clean_topo(void);
+
+/* skeleton.c */
+int thin_skeleton(double);
+
 /* sw_edgelist.c */
 int ELinitialize(void);
 struct Halfedge *HEcreate(struct Edge *, int);
@@ -119,5 +125,6 @@ int extend_line(double, double, double, double, double, double, double,
 		double, double, double *, double *, int);
 
 /* vo_write.c */
+int vo_write(void);
 int write_ep(struct Edge *);
 

vector/v.voronoi/sw_main.c

@@ -34,6 +34,7 @@ struct Map_info In, Out;
 int Type;
 int Field;
 int in_area;
+int skeleton;
 double segf;
 
 int nsites_alloc;

vector/v.voronoi/sw_voronoi.c

@@ -1,4 +1,5 @@
 #include <stdlib.h>
+#include <grass/gis.h>
 #include "sw_defs.h"
 
 /* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
@@ -14,6 +15,7 @@ int voronoi(struct Site *(*nextsite) (void))
     int pm;
     struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
     struct Edge *e;
+    int counter = 0;
 
     PQinitialize();
     bottomsite = (*nextsite) ();
@@ -29,6 +31,8 @@ int voronoi(struct Site *(*nextsite) (void))
 	    (newsite->coord.y == newintstar.y && 
 	    newsite->coord.x < newintstar.x))) {	/* new site is smallest */
 
+	    G_percent(counter++, nsites, 2);
+
 	    lbnd = ELleftbnd(&(newsite->coord));
 	    rbnd = ELright(lbnd);
 	    bot = rightreg(lbnd);
@@ -95,12 +99,7 @@ int voronoi(struct Site *(*nextsite) (void))
 	    break;
     }
 
-    for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd)) {
-	e = lbnd->ELedge;
-	write_ep(e);
-    }
-    
-    /* TODO: free memory */
+    G_percent(1, 1, 1);
 
     return 0;
 }

vector/v.voronoi/v.voronoi.html

@@ -1,61 +1,69 @@
 <h2>DESCRIPTION</h2>
 
-<em>v.voronoi</em> uses an existing vector points map (<b>input</b>) to create
-a Voronoi diagram (Thiessen polygons) in a new vector map (<b>output</b>).
-<p>The bounds of the output map are limited by the current region.
-(see <em>g.region</em>)
+<em>v.voronoi</em> creates a Voronoi diagram (Thiessen polygons) from 
+points or centroids. 
+<p>The bounds of the output map are limited by the current region 
+(see <em><a href="g.region.html">g.region</a></em>).
+
 <p>
-<br>
-Voronoi diagram and Delaunay triangulation example:
-<center>
-<img src="v_voronoi_delaunay.png" border="1"><br>
-<table border="0" width="590">
-<tr><td><center>
-<i>Delaunay Triangulation (left pane), Voronoi diagram (center pane),
-and both (right pane)</i>
-</center></td></tr>
-</table>
-</center>
+The <em>-a</em> flag can be used to create a Voronoi diagram for areas.
 
+<p>
+The <em>-s</em> flag can be used to extract the center line of areas or 
+skeletons of areas with <em>thin</em> >= 0. Smaller values for the 
+<em>thin</em> option will preserve more detail, while negative values 
+will extract only the center line.
 
 <h2>NOTES</h2>
 
 Voronoi diagrams may be used for nearest-neighbor flood filling.
 Give the centroids attributes (start with
 <em><a href="v.db.addcolumn.html">v.db.addcolumn</a></em>),
-then optionally convert to a raster map with
+then optionally convert the result to a raster map with
 <em><a href="v.to.rast.html">v.to.rast</a></em>.
 
 
 <h2>EXAMPLE</h2>
 
-Commands used with the Spearfish dataset to create the above figure.
+<h4>Voronoi diagram for points</h4>
+This example uses the hospitals in the North Carolina dataset.
+<div class="code"><pre>
+  g.region -p rast=elev_state_500m
+  v.voronoi in=hospitals out=hospitals_voronoi
+</pre></div>
+
+Result:
+<center>
+<img src="v_voronoi_points.png" border="0"><br>
+<i>Voronoi diagram for hospitals in North Carolina</i>
+</center>
+
+<h4>Voronoi diagram for areas</h4>
+This example uses urban areas in the North Carolina dataset.
 <div class="code"><pre>
-  g.region n=4927250 s=4919400 w=588650 e=594850
-  d.frame -c fr=one at=0,100,0,33.3333
-  d.frame -c fr=two at=0,100,33.3333,66.6667
-  d.frame -c fr=three at=0,100,66.6667,100
-
-  v.delaunay -lr in=archsites out=arch_delaunay
-  d.frame -s one
-  d.vect arch_delaunay
-  d.vect archsites color=red fcolor=red size=5 icon=basic/circle
-
-  v.voronoi -l in=archsites out=arch_voronoi
-  d.frame -s two
-  d.vect arch_voronoi type=line
-  d.vect archsites color=red fcolor=red size=5 icon=basic/circle
-
-  d.frame -s three
-  d.vect arch_voronoi type=line
-  d.vect arch_delaunay color=blue
-  d.vect archsites color=red fcolor=red size=5 icon=basic/circle
+  g.region -p n=162500 s=80000 w=727000 e=846000 res=500
+  v.voronoi in=urbanarea out=urbanarea_voronoi -a
 </pre></div>
 
+Result:
+<center>
+<img src="v_voronoi_areas.png" border="0"><br>
+<i>Voronoi diagram for urban areas in North Carolina</i>
+</center>
 
-<h2>BUGS</h2>
-Only attribute table of field 1 is copied. 
+<h4>Skeletons and center lines of areas</h4>
+This example uses urban areas in the North Carolina dataset.
+<div class="code"><pre>
+  g.region -p n=161000 s=135500 w=768500 e=805500 res=500
+  v.voronoi in=urbanarea out=urbanarea_centerline -s
+  v.voronoi in=urbanarea out=urbanarea_skeleton -s thin=2000
+</pre></div>
 
+Result:
+<center>
+<img src="v_voronoi_skeleton.png" border="0"><br>
+<i>Skeleton (blue) and center line (red) for urban areas in North Carolina</i>
+</center>
 
 <h2>REFERENCES</h2>
 <em>Steve J. Fortune,  (1987).  A Sweepline Algorithm for
@@ -73,6 +81,7 @@ Only attribute table of field 1 is copied.
 <h2>AUTHORS</h2>
 James Darrell McCauley, Purdue University<br>
 GRASS 5 update, improvements: <a href="mailto:aaime@libero.it">Andrea Aime</a>, Modena, Italy<br>
-GRASS 5.7 update: Radim Blazek
+GRASS 5.7 update: Radim Blazek<br>
+Markus Metz
 
 <p><i>Last changed: $Date$</i>

vector/v.voronoi/vo_write.c

@@ -7,6 +7,56 @@
 #include "sw_defs.h"
 #include "defs.h"
 
+static int write_skeleton(struct line_pnts *Points)
+{
+    int i, area1, area2;
+    static struct line_pnts *APoints = NULL;
+    static struct line_cats *Cats = NULL;
+
+    if (!APoints) {
+	APoints = Vect_new_line_struct();
+	Cats = Vect_new_cats_struct();
+    }
+
+    if ((area1 = Vect_find_area(&In, Points->x[0], Points->y[0])) == 0)
+	return 0;
+    if ((area2 = Vect_find_area(&In, Points->x[1], Points->y[1])) == 0)
+	return 0;
+    if (area1 != area2)
+	return 0;
+
+    if (!Vect_get_area_centroid(&In, area1))
+	return 0;
+    Vect_get_area_points(&In, area1, APoints);
+    if (Vect_line_check_intersection(Points, APoints, 0))
+	return 0;
+
+    for (i = 0; i < Vect_get_area_num_isles(&In, area1); i++) {
+	Vect_get_isle_points(&In, Vect_get_area_isle(&In, area1, i), APoints);
+	if (Vect_line_check_intersection(Points, APoints, 0))
+	    return 0;
+    }
+
+    Vect_get_area_cats(&In, area1, Cats);
+    Vect_write_line(&Out, GV_LINE, Points, Cats);
+
+    return 1;
+}
+
+int vo_write(void)
+{
+    struct Halfedge *lbnd;
+
+    for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd)) {
+	write_ep(lbnd->ELedge);
+    }
+
+    /* TODO: free memory */
+
+    return 1;
+}
+
+
 int write_ep(struct Edge *e)
 {
     static struct line_pnts *Points = NULL;
@@ -39,7 +89,10 @@ int write_ep(struct Edge *e)
 	Vect_reset_line(Points);
 	Vect_append_point(Points, x1, y1, 0.0);
 	Vect_append_point(Points, x2, y2, 0.0);
-	Vect_write_line(&Out, Type, Points, Cats);
+	if (skeleton)
+	    write_skeleton(Points);
+	else
+	    Vect_write_line(&Out, Type, Points, Cats);
     }
     else {
 	int knownPointAtLeft = -1;
@@ -80,7 +133,10 @@ int write_ep(struct Edge *e)
 	    Vect_reset_line(Points);
 	    Vect_append_point(Points, x1, y1, 0.0);
 	    Vect_append_point(Points, x2, y2, 0.0);
-	    Vect_write_line(&Out, Type, Points, Cats);
+	    if (skeleton)
+		write_skeleton(Points);
+	    else
+		Vect_write_line(&Out, Type, Points, Cats);
 	}
     }