v.voronoi bugfix for special cases

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

vector/v.voronoi/main.c

@@ -107,12 +107,14 @@ int main(int argc, char **argv)
     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;
 
     G_gisinit(argv[0]);
 
@@ -407,11 +409,22 @@ int main(int argc, char **argv)
     }
 
     /* 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, 1e-7, NULL);
+	Vect_snap_lines(&Out, GV_BOUNDARY, snap_thresh, NULL);
 	do {
 	    Vect_break_lines(&Out, GV_BOUNDARY, NULL);
 	    Vect_remove_duplicates(&Out, GV_BOUNDARY, NULL);
@@ -419,6 +432,11 @@ int main(int argc, char **argv)
 		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++) {

vector/v.voronoi/sw_defs.h

@@ -88,6 +88,7 @@ double dist(struct Site *, struct Site *);
 int makevertex(struct Site *);
 int deref(struct Site *);
 int ref(struct Site *);
+double d_ulp(double);
 
 /* sw_heap.c */
 int PQinsert(struct Halfedge *, struct Site *, double);

vector/v.voronoi/sw_geometry.c

@@ -67,14 +67,16 @@ struct Edge *bisect(struct Site *s1, struct Site *s2)
 }
 
 /* single precision ULP */
-static double d_ulp(double d)
+double d_ulp(double d)
 {
     int exp;
 
-
     if (d == 0)
 	return GRASS_EPSILON;
 
+    if (d < 0)
+	d = fabs(d);
+
     d = frexp(d, &exp);
     exp -= 22;
     d = ldexp(d, exp);

vector/v.voronoi/sw_main.c

@@ -126,7 +126,9 @@ int readsites(void)
 	if(ltype == -2)
 	    break;
 
-	G_percent(Vect_get_next_line_id(&In), nlines, 2);
+	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))
@@ -142,6 +144,9 @@ int readsites(void)
 	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;
@@ -160,10 +165,10 @@ int readsites(void)
 	nsites++;
     }
 
-    if (nsites < 3) {
+    if (nsites < 2) {
 	const char *name = Vect_get_full_name(&In);
 	Vect_close(&In);
-	G_fatal_error(_("Found %d points/centroids in <%s>, but at least 3 are needed"),
+	G_fatal_error(_("Found %d points/centroids in <%s>, but at least 2 are needed"),
 	              nsites, name);
     }
 
@@ -172,16 +177,6 @@ int readsites(void)
 	    (struct Site *)G_realloc(sites,
 				     (nsites) * sizeof(struct Site));
 
-    if (xmin == Box.W)
-	Box.W -= GRASS_EPSILON;
-    if (xmax == Box.E)
-	Box.E += GRASS_EPSILON;
-    if (ymin == Box.S)
-	Box.S -= GRASS_EPSILON;
-    if (ymax == Box.N)
-	Box.N += GRASS_EPSILON;
-
-
     qsort(sites, nsites, sizeof(struct Site), scomp);
     removeDuplicates();
 

vector/v.voronoi/sw_output.c

@@ -47,17 +47,9 @@ int out_bisector(struct Edge *e)
 
 int out_ep(struct Edge *e)
 {
-    static struct line_pnts *Points = NULL;
-    static struct line_cats *Cats = NULL;
-
-    if (!Points) {
-	Points = Vect_new_line_struct();
-	Cats = Vect_new_cats_struct();
-    }
-
-    if (!triangulate & plot)
+    if ((!triangulate) & plot)
 	clip_line(e);
-    if (!triangulate & !plot) {
+    if ((!triangulate) & !plot) {
 	/*
 	   fprintf (stdout,"e %d", e->edgenbr);
 	   fprintf (stdout," %d ", e->ep[le] != (struct Site *)NULL ? e->ep[le]->sitenbr : -1);
@@ -85,7 +77,7 @@ int out_vertex(struct Site *v)
 
 int out_site(struct Site *s)
 {
-    if (!triangulate & plot & !debug)
+    if ((!triangulate) & plot & (!debug))
 	circle(s->coord.x, s->coord.y, cradius);
 
     /*

vector/v.voronoi/vo_extend.c

@@ -33,7 +33,7 @@ int extend_line(double s, double n, double w, double e,
 {
     double nx, ny;		/* intersection coordinates */
 
-    if (x > w && x < e && y > s && y < n) {
+    if (x >= w && x <= e && y >= s && y <= n) {
 	/* vertical line? */
 	if (a == 0) {
 	    *c_x = knownPointAtLeft ? e : w;

vector/v.voronoi/vo_write.c

@@ -42,23 +42,43 @@ int write_ep(struct Edge *e)
 	    Vect_write_line(&Out, Type, Points, Cats);
 	}
 	else {
-	    int knownPointAtLeft;
+	    int knownPointAtLeft = -1;
 
 	    if (e->ep[le] != NULL) {
 		x1 = e->ep[le]->coord.x;
 		y1 = e->ep[le]->coord.y;
 		knownPointAtLeft = 1;
 	    }
-	    else {
+	    else if (e->ep[re] != NULL) {
 		x1 = e->ep[re]->coord.x;
 		y1 = e->ep[re]->coord.y;
 		knownPointAtLeft = 0;
 	    }
 
+	    if (knownPointAtLeft == -1) {
+
+		/*
+		G_warning("Dead edge!");
+		return 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;
+		knownPointAtLeft = 0;
+		if (!extend_line(Box.S, Box.N, Box.W, Box.E,
+				 e->a, e->b, e->c, x2, y2, &x1, &y1,
+				 knownPointAtLeft)) {
+				     
+		    G_warning("Undefined edge, unable to extend line");
+		    
+		    return 0;
+		}
+		G_debug(0, "x1 = %g, y1 = %g, x2 = %g, y2 = %g", x1, y1, x2, y2);
+		knownPointAtLeft = 1;
+	    }
 	    if (extend_line(Box.S, Box.N, Box.W, Box.E,
 			    e->a, e->b, e->c, x1, y1, &x2, &y2,
-			    knownPointAtLeft)
-		) {
+			    knownPointAtLeft)) {
+
 		/* Don't write zero length */
 		if (x1 == x2 && y1 == y2)
 		    return 0;