v.voronoi: add tesselation of areas

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

vector/v.voronoi/defs.h

@@ -3,5 +3,6 @@ extern struct Cell_head Window;
 extern struct bound_box Box;
 extern struct Map_info In, Out;
 extern int Type;
-extern int All;
 extern int Field;
+extern int in_area;
+extern double segf;

vector/v.voronoi/main.c

@@ -99,10 +99,10 @@ int main(int argc, char **argv)
     int i;
     int **cats, *ncats, nfields, *fields;
     struct {
-	struct Flag *line, *table;
+	struct Flag *line, *table, *area;
     } flag;
     struct {
-	struct Option *in, *out, *field;
+	struct Option *in, *out, *field, *segf;
     } opt;
     struct GModule *module;
     struct line_pnts *Points;
@@ -132,6 +132,18 @@ int main(int argc, char **argv)
     
     opt.out = G_define_standard_option(G_OPT_V_OUTPUT);
 
+    opt.segf = G_define_option();
+    opt.segf->type = TYPE_DOUBLE;
+    opt.segf->key = "segment";
+    opt.segf->answer = "0.25";
+    opt.segf->label = _("Factor to control boundary interpolation");
+    opt.segf->description = _("Applies to input areas only. Smaller values produce smoother output but can cause numerical instability.");
+
+    flag.area = G_define_flag();
+    flag.area->key = 'a';
+    flag.area->description =
+	_("Create Voronoi diagram for input areas");
+
     flag.line = G_define_flag();
     flag.line->key = 'l';
     flag.line->description =
@@ -147,7 +159,12 @@ int main(int argc, char **argv)
     else
 	Type = GV_BOUNDARY;
 
-    All = 0;
+    in_area = flag.area->answer;
+    segf = atof(opt.segf->answer);
+    if (segf < GRASS_EPSILON) {
+	segf = 0.25;
+	G_warning(_("Option '%s' is too small, set to %g"), opt.segf->key, segf);
+    }
 
     Points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
@@ -167,7 +184,10 @@ int main(int argc, char **argv)
     freeinit(&sfl, sizeof(struct Site));
 
     G_message(_("Reading features..."));
-    readsites();
+    if (in_area)
+	readbounds();
+    else
+	readsites();
 
     Vect_open_new(&Out, opt.out->answer, 0);
 
@@ -483,6 +503,16 @@ int main(int argc, char **argv)
 	}
     }
 
+    /* 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);
+    }
+    */
+
     /* build clean topology */
     Vect_build_partial(&Out, GV_BUILD_NONE);
     Vect_build(&Out);

vector/v.voronoi/sw_defs.h

@@ -103,7 +103,7 @@ int PQinitialize(void);
 int scomp(const void *, const void *);
 struct Site *nextone(void);
 int readsites(void);
-struct Site *readone(void);
+int readbounds(void);
 
 /* sw_memory.c */
 int freeinit(struct Freelist *, int);

vector/v.voronoi/sw_main.c

@@ -1,6 +1,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/vector.h>
 #include <grass/glocale.h>
@@ -31,8 +32,11 @@ struct Cell_head Window;
 struct bound_box Box;
 struct Map_info In, Out;
 int Type;
-int All;
 int Field;
+int in_area;
+double segf;
+
+int nsites_alloc;
 
 /* sort sites on y, then x, coord */
 int scomp(const void *v1, const void *v2)
@@ -103,12 +107,48 @@ void removeDuplicates()
 
 }
 
+int addsite(double x, double y, double z, int id)
+{
+    if (nsites >= nsites_alloc) {
+	nsites_alloc += 100;
+	sites =
+	    (struct Site *)G_realloc(sites,
+				     (nsites_alloc) * sizeof(struct Site));
+    }
+    sites[nsites].coord.x = x;
+    sites[nsites].coord.y = y;
+    sites[nsites].coord.z = z;
+
+    sites[nsites].sitenbr = id;
+    sites[nsites].refcnt = 0;
+
+    if (nsites > 0) {
+	if (xmin > sites[nsites].coord.x)
+	    xmin = sites[nsites].coord.x;
+	if (xmax < sites[nsites].coord.x)
+	    xmax = sites[nsites].coord.x;
+	if (ymin > sites[nsites].coord.y)
+	    ymin = sites[nsites].coord.y;
+	if (ymax < sites[nsites].coord.y)
+	    ymax = sites[nsites].coord.y;
+    }
+    else {
+	xmin = xmax = sites[nsites].coord.x;
+	ymin = ymax = sites[nsites].coord.y;
+    }
+
+    nsites++;
+
+    return nsites;
+}
+
 /* read all sites, sort, and compute xmin, xmax, ymin, ymax */
 int readsites(void)
 {
     int nlines, ltype;
     struct line_pnts *Points;
     struct line_cats *Cats;
+    double z = 0.;
     
     Points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
@@ -129,40 +169,16 @@ int readsites(void)
 	if (!(ltype & GV_POINTS))
 	    continue;
 	/* G_percent(Vect_get_next_line_id(&In), nlines, 2); */
-		
-	if (!All) {
-	    if (!Vect_point_in_box(Points->x[0], Points->y[0], 0.0, &Box))
-		continue;
-	}
 
-	sites[nsites].coord.x = Points->x[0];
-	sites[nsites].coord.y = Points->y[0];
+	if (!Vect_point_in_box(Points->x[0], Points->y[0], 0.0, &Box))
+	    continue;
+
 	if (mode3d) {
 	    G_debug(3, "Points->z[0]: %f", Points->z[0]);
-	    sites[nsites].coord.z = Points->z[0];
-	}
-	else
-	    sites[nsites].coord.z = 0.0;
-
-	sites[nsites].sitenbr = nsites;
-	sites[nsites].refcnt = 0;
-
-	if (nsites > 1) {
-	    if (xmin > sites[nsites].coord.x)
-		xmin = sites[nsites].coord.x;
-	    if (xmax < sites[nsites].coord.x)
-		xmax = sites[nsites].coord.x;
-	    if (ymin > sites[nsites].coord.y)
-		ymin = sites[nsites].coord.y;
-	    if (ymax < sites[nsites].coord.y)
-		ymax = sites[nsites].coord.y;
-	}
-	else {
-	    xmin = xmax = sites[nsites].coord.x;
-	    ymin = ymax = sites[nsites].coord.y;
+	    z = Points->z[0];
 	}
 
-	nsites++;
+	addsite(Points->x[0], Points->y[0], z, nsites);
     }
 
     if (nsites < 2) {
@@ -186,18 +202,278 @@ int readsites(void)
     return 0;
 }
 
-/* read one site */
-struct Site *readone(void)
+/* valid boundary: one area with centroid */
+int n_areas(int line, int *aid)
 {
-    struct Site *s;
+    int larea, rarea, ncentroids;
+    
+    ncentroids = 0;
+    Vect_get_line_areas(&In, line, &larea, &rarea);
+    
+    if (larea < 0)
+	larea = Vect_get_isle_area(&In, -larea);
+    if (larea > 0) {
+	if (Vect_get_area_centroid(&In, larea) > 0) {
+	    ncentroids++;
+	    *aid = larea;
+	}
+    }
+    if (rarea < 0)
+	rarea = Vect_get_isle_area(&In, -rarea);
+    if (rarea > 0) {
+	if (Vect_get_area_centroid(&In, rarea) > 0) {
+	    ncentroids++;
+	    *aid = rarea;
+	}
+    }
+    
+    return ncentroids;
+}
+
+/* read all boundaries, sort, and compute xmin, xmax, ymin, ymax */
+int readbounds(void)
+{
+    int line, nlines, ltype, node, nnodes;
+    struct line_pnts *Points;
+    struct line_cats *Cats;
+    int i;
+    int area_id;
+    double x, y, z, x1, y1, z1;
+    double maxdist, sdist, l, dx, dy, dz;
+    int nconnected;
+    struct ilist *linelist, *arealist;
+    
+    Points = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+    
+    nlines = Vect_get_num_lines(&In);
 
-    s = (struct Site *)getfree(&sfl);
-    s->refcnt = 0;
-    s->sitenbr = siteidx;
-    siteidx++;
+    nsites = 0;
+    nsites_alloc = nlines * 2;
+    sites = (struct Site *)G_malloc(nsites_alloc * sizeof(struct Site));
 
-    if (scanf("%lf %lf", &(s->coord.x), &(s->coord.y)) == EOF)
-	return ((struct Site *)NULL);
+    Vect_set_constraint_type(&In, GV_BOUNDARY);
+    Vect_set_constraint_field(&In, Field);
+    
+    l = 0;
+    maxdist = 0;
+    for (line = 1; line <= nlines; line++) {
+
+	if (!Vect_line_alive(&In, line))
+	    continue;
+	ltype = Vect_get_line_type(&In, line);
+
+	if (!(ltype & GV_BOUNDARY))
+	    continue;
+
+	area_id = 0;
+	if (n_areas(line, &area_id) != 1)
+	    continue;
+
+	Vect_read_line(&In, Points, Cats, line);
+	Vect_line_prune(Points);
+	
+	l += Vect_line_length(Points);
+	nsites += Points->n_points;
+    }
+    if (nsites)
+	maxdist = segf * l / nsites;
+    G_verbose_message("Maximum segment length: %g", maxdist);
+    
+    nsites = 0;
+    z = 0.;
+    z1 = 0;
+    dz = 0;
+    for (line = 1; line <= nlines; line++) {
+
+	if (!Vect_line_alive(&In, line))
+	    continue;
+	ltype = Vect_get_line_type(&In, line);
+
+	if (!(ltype & GV_BOUNDARY))
+	    continue;
+
+	area_id = 0;
+	if (n_areas(line, &area_id) != 1)
+	    continue;
+
+	Vect_read_line(&In, Points, Cats, line);
+	Vect_line_prune(Points);
+
+	if (nsites + Points->n_points > nsites_alloc) {
+	    nsites_alloc = nsites + Points->n_points;
+	    sites =
+		(struct Site *)G_realloc(sites,
+					 (nsites_alloc) * sizeof(struct Site));
+	}
+
+	for (i = 0; i < Points->n_points; i++) {
+	    if (!Vect_point_in_box(Points->x[i], Points->y[i], 0.0, &Box))
+		continue;
+
+	    x = Points->x[i];
+	    y = Points->y[i];
+	    if (mode3d) {
+		G_debug(3, "Points->z[i]: %f", Points->z[i]);
+		z = Points->z[i];
+	    }
+
+	    if (i > 0 && i < Points->n_points - 1)
+		addsite(x, y, z, area_id);
+	    
+	    /* densify */
+	    if (maxdist > 0 && i < Points->n_points - 1) {
+		dx = Points->x[i + 1] - Points->x[i];
+		dy = Points->y[i + 1] - Points->y[i];
+		if (mode3d)
+		    dz = Points->z[i + 1] - Points->z[i];
+		l = sqrt(dx * dx + dy * dy);
+		
+		if (l > maxdist) {
+		    int n = ceil(l / maxdist) + 0.5;
+		    double step = l / n;
+		    
+		    while (--n) {
+			sdist = (step * n) / l;
+			x1 = x + sdist * dx;
+			y1 = y + sdist * dy;
+			if (mode3d)
+			    z1 = z + sdist * dz;
+
+			addsite(x1, y1, z1, area_id);
+		    }
+		}
+	    }
+	}
+    }
+    
+    /* process nodes */
+    nnodes = Vect_get_num_nodes(&In);
+    linelist = Vect_new_list();
+    arealist = Vect_new_list();
+    
+    for (node = 1; node <= nnodes; node++) {
+	Vect_get_node_coor(&In, node, &x, &y, &z);
+	if (!mode3d)
+	    z = 0.;
+
+	if (!Vect_point_in_box(x, y, 0.0, &Box))
+	    continue;
+
+	/* count number of connected boundaries
+	 * must be >= 2 */
+	/* count number of valid boundaries */
+	nlines = Vect_get_node_n_lines(&In, node);
+	nconnected = 0;
+	Vect_reset_list(linelist);
+	Vect_reset_list(arealist);
+	
+	for (i = 0; i < nlines; i++) {
+	    line = Vect_get_node_line(&In, node, i);
+	    
+	    ltype = Vect_get_line_type(&In, abs(line));
+	    
+	    if (!(ltype & GV_BOUNDARY))
+		continue;
+	    
+	    if (n_areas(abs(line), &area_id) == 1) {
+		Vect_list_append(linelist, line);
+		Vect_list_append(arealist, area_id);
+		nconnected++;
+	    }
+	}
+	
+	if (arealist->n_values == 1) {
+	    
+	    area_id = arealist->value[0];
+	    addsite(x, y, z, area_id);
+	}
+	else if (arealist->n_values > 1) {
+	    /* displacement */
+	    double displace;
+	    
+	    if (fabs(x) > fabs(y))
+		displace = d_ulp(fabs(x));
+	    else
+		displace = d_ulp(fabs(y));
+	    
+	    displace *= 2;
+	    
+	    for (i = 0; i < linelist->n_values; i++) {
 
-    return (s);
+		line = linelist->value[i];
+		
+		if (n_areas(abs(line), &area_id) != 1)
+		    G_fatal_error(_("All boundaries in the list should be valid"));
+
+		Vect_read_line(&In, Points, Cats, abs(line));
+		Vect_line_prune(Points);
+		
+		if (Points->n_points < 2)
+		    G_fatal_error(_("Boundary is degenerate"));
+
+		if (line < 0) {
+		    dx = Points->x[Points->n_points - 2] - Points->x[Points->n_points - 1]; 
+		    dy = Points->y[Points->n_points - 2] - Points->y[Points->n_points - 1]; 
+		}
+		else {
+		    dx = Points->x[1] - Points->x[0]; 
+		    dy = Points->y[1] - Points->y[0]; 
+		}
+		l = sqrt(dx * dx + dy * dy);
+		if (displace > l)
+		    displace = l;
+	    }
+
+	    for (i = 0; i < linelist->n_values; i++) {
+
+		line = linelist->value[i];
+		
+		if (n_areas(abs(line), &area_id) != 1)
+		    G_fatal_error(_("All boundaries in the list should be valid"));
+
+		Vect_read_line(&In, Points, Cats, abs(line));
+		Vect_line_prune(Points);
+		
+		if (Points->n_points < 2)
+		    G_fatal_error(_("Boundary is degenerate"));
+
+		if (line < 0) {
+		    dx = Points->x[Points->n_points - 2] - Points->x[Points->n_points - 1]; 
+		    dy = Points->y[Points->n_points - 2] - Points->y[Points->n_points - 1]; 
+		}
+		else {
+		    dx = Points->x[1] - Points->x[0]; 
+		    dy = Points->y[1] - Points->y[0]; 
+		}
+		l = sqrt(dx * dx + dy * dy);
+		if (l > displace * 2) {
+		    sdist = displace / l;
+		    x1 = x + sdist * dx;
+		    y1 = y + sdist * dy;
+		    z1 = 0;
+
+		    addsite(x1, y1, z1, area_id);
+		}
+	    }
+	}
+    }
+
+    if (nsites < 2) {
+	const char *name = Vect_get_full_name(&In);
+	Vect_close(&In);
+	G_fatal_error(_("Found %d vertices in <%s>, but at least 2 are needed"),
+	              nsites, name);
+    }
+
+
+    qsort(sites, nsites, sizeof(struct Site), scomp);
+    removeDuplicates();
+
+    Vect_destroy_line_struct(Points);
+    Vect_destroy_cats_struct(Cats);
+    Vect_destroy_list(linelist);
+    Vect_destroy_list(arealist);
+
+    return 0;
 }

vector/v.voronoi/sw_voronoi.c

@@ -25,8 +25,11 @@ int voronoi(int triangulate, struct Site *(*nextsite) (void))
 	if (!PQempty())
 	    newintstar = PQ_min();
 
-	if (newsite != (struct Site *)NULL && (PQempty()
-					       || newsite->coord.y < newintstar.y || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x))) {	/* new site is smallest */
+	if (newsite != (struct Site *)NULL && 
+	    (PQempty() || newsite->coord.y < newintstar.y || 
+	    (newsite->coord.y == newintstar.y && 
+	    newsite->coord.x < newintstar.x))) {	/* new site is smallest */
+
 	    out_site(newsite);
 	    lbnd = ELleftbnd(&(newsite->coord));
 	    rbnd = ELright(lbnd);
@@ -54,9 +57,8 @@ int voronoi(int triangulate, struct Site *(*nextsite) (void))
 		   temp->coord.y == newsite->coord.y);
 	    newsite = temp;
 	}
-	else if (!PQempty())
+	else if (!PQempty()) {
 	    /* intersection is smallest */
-	{
 	    lbnd = PQextractmin();
 	    llbnd = ELleft(lbnd);
 	    rbnd = ELright(lbnd);

vector/v.voronoi/vo_write.c

@@ -20,6 +20,9 @@ int write_ep(struct Edge *e)
 
     if (!triangulate) {
 	double x1, y1, x2, y2;
+	
+	if (in_area && e->reg[le]->sitenbr == e->reg[re]->sitenbr)
+	    return 0;
 
 	if (e->ep[le] != NULL && e->ep[re] != NULL) {	/* both end defined */
 	    x1 = e->ep[le]->coord.x;