[GRASS-SVN] r57268 - grass/trunk/vector/v.voronoi

[svn_grass at osgeo.org]

Author: mmetz
Date: 2013-07-25 04:14:07 -0700 (Thu, 25 Jul 2013)
New Revision: 57268

Modified:
   grass/trunk/vector/v.voronoi/main.c
   grass/trunk/vector/v.voronoi/sw_geometry.c
Log:
v.voronoi: fix #957, #1682

Modified: grass/trunk/vector/v.voronoi/main.c
===================================================================
--- grass/trunk/vector/v.voronoi/main.c	2013-07-25 10:30:27 UTC (rev 57267)
+++ grass/trunk/vector/v.voronoi/main.c	2013-07-25 11:14:07 UTC (rev 57268)
@@ -187,7 +187,7 @@
     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);
@@ -410,7 +410,7 @@
     if (err_nocentr || err_centr_dupl || err_centr_out) {
 	int nmod;
 
-	G_important_message(_("Output needs topological cleaning"));
+	G_important_message(_("Cleaning output topology"));
 	Vect_snap_lines(&Out, GV_BOUNDARY, 1e-7, NULL);
 	do {
 	    Vect_break_lines(&Out, GV_BOUNDARY, NULL);

Modified: grass/trunk/vector/v.voronoi/sw_geometry.c
===================================================================
--- grass/trunk/vector/v.voronoi/sw_geometry.c	2013-07-25 10:30:27 UTC (rev 57267)
+++ grass/trunk/vector/v.voronoi/sw_geometry.c	2013-07-25 11:14:07 UTC (rev 57268)
@@ -1,5 +1,6 @@
 #include <stdlib.h>
 #include <math.h>
+#include <grass/gis.h>
 #include "sw_defs.h"
 
 int geominit(void)
@@ -32,12 +33,22 @@
     newedge->ep[0] = (struct Site *)NULL;
     newedge->ep[1] = (struct Site *)NULL;
 
-    dx = s2->coord.x - s1->coord.x;
-    dy = s2->coord.y - s1->coord.y;
+    if (s1->coord.x < s2->coord.x || 
+        (s1->coord.x == s2->coord.x && s1->coord.y < s2->coord.y)) {
+	dx = s2->coord.x - s1->coord.x;
+	dy = s2->coord.y - s1->coord.y;
+	newedge->c =
+	    s1->coord.x * dx + s1->coord.y * dy + (dx * dx + dy * dy) * 0.5;
+    }
+    else {
+	dx = s1->coord.x - s2->coord.x;
+	dy = s1->coord.y - s2->coord.y;
+	newedge->c =
+	    s2->coord.x * dx + s2->coord.y * dy + (dx * dx + dy * dy) * 0.5;
+    }
+
     adx = dx > 0 ? dx : -dx;
     ady = dy > 0 ? dy : -dy;
-    newedge->c =
-	s1->coord.x * dx + s1->coord.y * dy + (dx * dx + dy * dy) * 0.5;
     if (adx > ady) {
 	newedge->a = 1.0;
 	newedge->b = dy / dx;
@@ -55,12 +66,27 @@
     return (newedge);
 }
 
+/* single precision ULP */
+static double d_ulp(double d)
+{
+    int exp;
 
+
+    if (d == 0)
+	return GRASS_EPSILON;
+
+    d = frexp(d, &exp);
+    exp -= 22;
+    d = ldexp(d, exp);
+
+    return d;
+}
+
 struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2)
 {
     struct Edge *e1, *e2, *e;
     struct Halfedge *el;
-    double d, xint, yint;
+    double d, dt, xint, yint;
     int right_of_site;
     struct Site *v;
 
@@ -72,9 +98,21 @@
 	return ((struct Site *)NULL);
 
     d = e1->a * e2->b - e1->b * e2->a;
-    if (-1.0e-10 < d && d < 1.0e-10)
+    if (fabs(e1->a * e2->b) > fabs(e1->b * e2->a)) {
+	dt = fabs(e1->a * e2->b);
+    }
+    else
+	dt = fabs(e1->b * e2->a);
+
+    if (dt != dt)
 	return ((struct Site *)NULL);
 
+    dt = d_ulp(dt);
+    G_debug(4, "dt = %g", dt);
+
+    if (-dt < d && d < dt)
+	return ((struct Site *)NULL);
+
     xint = (e1->c * e2->b - e2->c * e1->b) / d;
     yint = (e2->c * e1->a - e1->c * e2->a) / d;