v.voronoi for areas (#26)

v.voronoi

Fix for area skeletons and voronoi diagrams for areas.

Numerical stability of v.voronoi has been improved but is not perfect.

vector/v.voronoi/defs.h

@@ -1,7 +1,7 @@
 
 extern struct Cell_head Window;
 extern struct bound_box Box;
-extern struct Map_info In, Out;
+extern struct Map_info In, Out, Out2;
 extern int Type;
 extern int Field;
 extern int in_area;

vector/v.voronoi/main.c

@@ -202,6 +202,8 @@ int main(int argc, char **argv)
     Vect_region_box(&Window, &Box);
     Box.T = 0.5;
     Box.B = -0.5;
+    xcenter = (Box.W + Box.E) / 2.0;
+    ycenter = (Box.S + Box.N) / 2.0;
 
     freeinit(&sfl, sizeof(struct Site));
 
@@ -211,8 +213,17 @@ int main(int argc, char **argv)
     else
 	readsites();
 
-    if (Vect_open_new(&Out, opt.out->answer, 0) < 0)
-	G_fatal_error(_("Unable to create vector map <%s>"), opt.out->answer);
+    if (in_area) {
+	if (Vect_open_tmp_new(&Out, NULL, 0) < 0)
+	    G_fatal_error(_("Unable to create temporary vector map"));
+
+	if (Vect_open_new(&Out2, opt.out->answer, 0) < 0)
+	    G_fatal_error(_("Unable to create vector map <%s>"), opt.out->answer);
+    }
+    else {
+	if (Vect_open_new(&Out, opt.out->answer, 0) < 0)
+	    G_fatal_error(_("Unable to create vector map <%s>"), opt.out->answer);
+    }
 
     Vect_hist_copy(&In, &Out);
     Vect_hist_command(&Out);
@@ -299,6 +310,103 @@ int main(int argc, char **argv)
 	G_free(coor);
     }
 
+    if (in_area) {
+	int ltype;
+
+	/* Voronoi diagrams for areas */
+	/* write input boundaries with an area on both sides to output */
+	if (copybounds() > 0) {
+	    int cat;
+	    double x, y;
+	    double snap_thresh;
+	    int nmod;
+
+	    /* snap Voronoi edges to boundaries */
+	    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 = 2 * d_ulp(snap_thresh);
+
+	    /* TODO: snap only Voronoi edges */
+	    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);
+
+	    /* break lines */
+	    Vect_break_lines(&Out, GV_BOUNDARY, NULL);
+
+	    /* remove edge parts that are within an input area */
+	    nlines = Vect_get_num_lines(&Out);
+	    for (line = 1; line <= nlines; line++) {
+		if (!Vect_line_alive(&Out, line))
+		    continue;
+
+		ltype = Vect_get_line_type(&Out, line);
+		if (!(ltype & GV_BOUNDARY))
+		    continue;
+
+		Vect_read_line(&Out, Points, Cats, line);
+		cat = 0;
+		Vect_cat_get(Cats, 1, &cat);
+		if (cat != 1)
+		    continue;
+
+		Vect_line_prune(Points);
+		if (Points->n_points == 2) {
+		    x = (Points->x[0] + Points->x[1]) / 2.0;
+		    y = (Points->y[0] + Points->y[1]) / 2.0;
+		}
+		else {
+		    i = 1;
+		    while (Points->x[0] == Points->x[i] &&
+		           Points->y[0] == Points->y[i] &&
+			   i < Points->n_points - 1) {
+			i++;
+		    }
+		    i = Points->n_points / 2.;
+		    x = Points->x[i];
+		    y = Points->y[i];
+		}
+		if (Vect_find_area(&In, x, y) > 0)
+		    Vect_delete_line(&Out, line);
+	    }
+	}
+
+	/* copy result to final output */
+	Vect_reset_cats(Cats);
+
+	nlines = Vect_get_num_lines(&Out);
+
+	for (line = 1; line <= nlines; line++) {
+
+	    if (!Vect_line_alive(&Out, line))
+		continue;
+
+	    ltype = Vect_get_line_type(&Out, line);
+	    if (!(ltype & GV_BOUNDARY))
+		continue;
+
+	    Vect_read_line(&Out, Points, NULL, line);
+	    Vect_line_prune(Points);
+
+	    Vect_write_line(&Out2, GV_BOUNDARY, Points, Cats);
+	}
+
+	/* close temporary vector */
+	Vect_close(&Out);
+	/* is this a dirty hack ? */
+	Out = Out2;
+	Vect_build_partial(&Out, GV_BUILD_BASE);
+    }
+
     nfields = Vect_cidx_get_num_fields(&In);
     cats = (int **)G_malloc(nfields * sizeof(int *));
     ncats = (int *)G_malloc(nfields * sizeof(int));

vector/v.voronoi/sw_defs.h

@@ -57,6 +57,7 @@ extern struct Site *bottomsite;
 extern int nedges;
 extern struct Freelist efl;
 extern double xmin, xmax, ymin, ymax, deltax, deltay;
+extern double xcenter, ycenter;
 extern struct Freelist hfl;
 extern struct Halfedge *ELleftend, *ELrightend;
 extern int ELhashsize;
@@ -111,6 +112,7 @@ int scomp(const void *, const void *);
 struct Site *nextone(void);
 int readsites(void);
 int readbounds(void);
+int copybounds(void);
 
 /* sw_memory.c */
 int freeinit(struct Freelist *, int);

vector/v.voronoi/sw_main.c

@@ -19,6 +19,7 @@ struct Site *bottomsite;
 int nedges;
 struct Freelist efl;
 double xmin, xmax, ymin, ymax, deltax, deltay;
+double xcenter, ycenter;
 struct Freelist hfl;
 struct Halfedge *ELleftend, *ELrightend;
 int ELhashsize;
@@ -30,7 +31,7 @@ int PQmin;
 
 struct Cell_head Window;
 struct bound_box Box;
-struct Map_info In, Out;
+struct Map_info In, Out, Out2;
 int Type;
 int Field;
 int in_area;
@@ -138,6 +139,9 @@ int addsite(double x, double y, double z, int id)
 	ymin = ymax = sites[nsites].coord.y;
     }
 
+    sites[nsites].coord.x -= xcenter;
+    sites[nsites].coord.y -= ycenter;
+
     nsites++;
 
     return nsites;
@@ -259,7 +263,7 @@ int readbounds(void)
 
     Vect_set_constraint_type(&In, GV_BOUNDARY);
     Vect_set_constraint_field(&In, Field);
-    
+
     l = 0;
     maxdist = 0;
     for (line = 1; line <= nlines; line++) {
@@ -272,7 +276,7 @@ int readbounds(void)
 	    continue;
 
 	area_id = 0;
-	if (n_areas(line, &area_id) != 1)
+	if (n_areas(line, &area_id) < 1)
 	    continue;
 
 	Vect_read_line(&In, Points, Cats, line);
@@ -299,8 +303,14 @@ int readbounds(void)
 	    continue;
 
 	area_id = 0;
-	if (n_areas(line, &area_id) != 1)
-	    continue;
+	if (skeleton) {
+	    if (n_areas(line, &area_id) < 1)
+		continue;
+	}
+	else {
+	    if (n_areas(line, &area_id) != 1)
+		continue;
+	}
 
 	Vect_read_line(&In, Points, Cats, line);
 	Vect_line_prune(Points);
@@ -484,3 +494,44 @@ int readbounds(void)
 
     return 0;
 }
+
+int copybounds(void)
+{
+    int line, nlines, ltype;
+    int area_id;
+    int ncopied;
+    struct line_pnts *Points;
+    struct line_cats *Cats;
+
+    Points = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+    
+    nlines = Vect_get_num_lines(&In);
+    ncopied = 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;
+
+	if (n_areas(line, &area_id) != 2)
+	    continue;
+
+	Vect_read_line(&In, Points, Cats, line);
+	Vect_line_prune(Points);
+
+	Vect_reset_cats(Cats);
+	Vect_cat_set(Cats, 1, 2);
+	Vect_write_line(&Out, GV_BOUNDARY, Points, Cats);
+	ncopied++;
+    }
+
+    Vect_destroy_line_struct(Points);
+    Vect_destroy_cats_struct(Cats);
+
+    return ncopied;
+}

vector/v.voronoi/v.voronoi.html

@@ -16,6 +16,12 @@ will extract only the center line.
 
 <h2>NOTES</h2>
 
+<em>v.voronoi</em> suffers from numerical instability, results can 
+sometimes contain many artefacts. When creating Voronoi diagrams for areas 
+or skeletons for areas, it is highly recommended to simplify the areas first 
+with <em><a href="v.generalize.html">v.generalize</a></em>.
+
+<p>
 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>),

vector/v.voronoi/vo_extend.c

@@ -48,7 +48,7 @@ int extend_line(double s, double n, double w, double e,
 	}
 
 	/* south */
-	nx = (c - b * s) / a;
+	nx = (c - b * (s - ycenter)) / a + xcenter;
 	if (Vect_point_in_box(nx, s, 0.0, &Box) &&
 	    ((nx > x && knownPointAtLeft) || (nx <= x && !knownPointAtLeft)))
 	{
@@ -58,7 +58,7 @@ int extend_line(double s, double n, double w, double e,
 	}
 
 	/* north */
-	nx = (c - b * n) / a;
+	nx = (c - b * (n - ycenter)) / a + xcenter;
 	if (Vect_point_in_box(nx, n, 0.0, &Box) &&
 	    ((nx > x && knownPointAtLeft) || (nx <= x && !knownPointAtLeft)))
 	{
@@ -69,7 +69,7 @@ int extend_line(double s, double n, double w, double e,
 
 	if (knownPointAtLeft) {
 	    /* east */
-	    ny = (c - a * e) / b;
+	    ny = (c - a * (e - xcenter)) / b + ycenter;
 	    if (Vect_point_in_box(e, ny, 0.0, &Box)) {
 		*c_x = e;
 		*c_y = ny;
@@ -78,7 +78,7 @@ int extend_line(double s, double n, double w, double e,
 	}
 	else {
 	    /* west */
-	    ny = (c - a * w) / b;
+	    ny = (c - a * (w - xcenter)) / b + ycenter;
 	    if (Vect_point_in_box(w, ny, 0.0, &Box)) {
 		*c_x = w;
 		*c_y = ny;

vector/v.voronoi/vo_write.c

@@ -12,6 +12,7 @@ static int write_skeleton(struct line_pnts *Points)
     int i, area1, area2;
     static struct line_pnts *APoints = NULL;
     static struct line_cats *Cats = NULL;
+    int ret;
 
     if (!APoints) {
 	APoints = Vect_new_line_struct();
@@ -22,21 +23,51 @@ static int write_skeleton(struct line_pnts *Points)
 	return 0;
     if ((area2 = Vect_find_area(&In, Points->x[1], Points->y[1])) == 0)
 	return 0;
-    if (area1 != area2)
-	return 0;
 
+    ret = 1;
     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;
+	ret = 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 (ret) {
+	Vect_get_area_points(&In, area1, APoints);
 	if (Vect_line_check_intersection(Points, APoints, 0))
-	    return 0;
+	    ret = 0;
+    }
+
+    if (ret) {
+	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)) {
+		ret = 0;
+		break;
+	    }
+	}
     }
 
+
+    if (!ret && area1 != area2) {
+	ret = 1;
+	if (!Vect_get_area_centroid(&In, area2))
+	    ret = 0;
+	if (ret) {
+	    Vect_get_area_points(&In, area2, APoints);
+	    if (Vect_line_check_intersection(Points, APoints, 0))
+		ret = 0;
+	}
+
+	if (ret) {
+	    for (i = 0; i < Vect_get_area_num_isles(&In, area2); i++) {
+		Vect_get_isle_points(&In, Vect_get_area_isle(&In, area2, i), APoints);
+		if (Vect_line_check_intersection(Points, APoints, 0)) {
+		    ret = 0;
+		    break;
+		}
+	    }
+	}
+    }
+    if (!ret)
+	return 0;
+
     Vect_get_area_cats(&In, area1, Cats);
     Vect_write_line(&Out, GV_LINE, Points, Cats);
 
@@ -66,16 +97,18 @@ int write_ep(struct Edge *e)
     if (!Points) {
 	Points = Vect_new_line_struct();
 	Cats = Vect_new_cats_struct();
+	if (in_area)
+	    Vect_cat_set(Cats, 1, 1);
     }
 
     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;
-	y1 = e->ep[le]->coord.y;
-	x2 = e->ep[re]->coord.x;
-	y2 = e->ep[re]->coord.y;
+	x1 = e->ep[le]->coord.x + xcenter;
+	y1 = e->ep[le]->coord.y + ycenter;
+	x2 = e->ep[re]->coord.x + xcenter;
+	y2 = e->ep[re]->coord.y + ycenter;
 
 	if (!Vect_point_in_box(x1, y1, 0.0, &Box) ||
 	    !Vect_point_in_box(x2, y2, 0.0, &Box)) {
@@ -98,19 +131,19 @@ int write_ep(struct Edge *e)
 	int knownPointAtLeft = -1;
 
 	if (e->ep[le] != NULL) {
-	    x1 = e->ep[le]->coord.x;
-	    y1 = e->ep[le]->coord.y;
+	    x1 = e->ep[le]->coord.x + xcenter;
+	    y1 = e->ep[le]->coord.y + ycenter;
 	    knownPointAtLeft = 1;
 	}
 	else if (e->ep[re] != NULL) {
-	    x1 = e->ep[re]->coord.x;
-	    y1 = e->ep[re]->coord.y;
+	    x1 = e->ep[re]->coord.x + xcenter;
+	    y1 = e->ep[re]->coord.y + ycenter;
 	    knownPointAtLeft = 0;
 	}
 
 	if (knownPointAtLeft == -1) {
-	    x2 = (e->reg[le]->coord.x + e->reg[re]->coord.x) / 2.0;
-	    y2 = (e->reg[le]->coord.y + e->reg[re]->coord.y) / 2.0;
+	    x2 = (e->reg[le]->coord.x + e->reg[re]->coord.x) / 2.0 + xcenter;
+	    y2 = (e->reg[le]->coord.y + e->reg[re]->coord.y) / 2.0 + ycenter;
 	    knownPointAtLeft = 0;
 	    if (!extend_line(Box.S, Box.N, Box.W, Box.E,
 			     e->a, e->b, e->c, x2, y2, &x1, &y1,