[GRASS-SVN] r57844 - grass/branches/releasebranch_6_4/lib/gis

[svn_grass at osgeo.org]

Author: mmetz
Date: 2013-09-26 03:30:53 -0700 (Thu, 26 Sep 2013)
New Revision: 57844

Modified:
   grass/branches/releasebranch_6_4/lib/gis/plot.c
Log:
libgis: backport fix FPE error tolerance for plotting

Modified: grass/branches/releasebranch_6_4/lib/gis/plot.c
===================================================================
--- grass/branches/releasebranch_6_4/lib/gis/plot.c	2013-09-26 10:27:51 UTC (rev 57843)
+++ grass/branches/releasebranch_6_4/lib/gis/plot.c	2013-09-26 10:30:53 UTC (rev 57844)
@@ -485,16 +485,20 @@
     }
 
     /* check if perimeter has odd number of points */
-    if (np % 2)
+    if (np % 2) {
+	G_debug(1, "Weird internal error: perimeter has odd number of points");
 	return OUT_OF_SYNC;
+    }
 
     /* sort the edge points by col(x) and then by row(y) */
     qsort(P, np, sizeof(POINT), &edge_order);
 
     /* plot */
     for (i = 1; i < np; i += 2) {
-	if (P[i].y != P[i - 1].y)
+	if (P[i].y != P[i - 1].y) {
+	    G_debug(1, "Weird internal error: edge leaves row");
 	    return OUT_OF_SYNC;
+	}
 	row_fill(P[i].y, P[i - 1].x + shift1, P[i].x + shift1);
     }
     if (window.proj == PROJECTION_LL) {	/* now do wrap-around, part 2 */
@@ -622,8 +626,10 @@
     }				/* for() */
 
     /* check if perimeter has odd number of points */
-    if (np % 2)
+    if (np % 2) {
+	G_debug(1, "Weird internal error: perimeter has odd number of points");
 	return OUT_OF_SYNC;
+    }
 
     /* sort the edge points by col(x) and then by row(y) */
     qsort(P, np, sizeof(POINT), &edge_order);
@@ -631,8 +637,10 @@
     /* plot */
     for (j = 0; j < rings; j++) {
 	for (i = 1; i < np; i += 2) {
-	    if (P[i].y != P[i - 1].y)
+	    if (P[i].y != P[i - 1].y) {
+		G_debug(1, "Weird internal error: edge leaves row");
 		return OUT_OF_SYNC;
+	    }
 	    row_fill(P[i].y, P[i - 1].x + shift1[j], P[i].x + shift1[j]);
 	}
 	if (window.proj == PROJECTION_LL) {	/* now do wrap-around, part 2 */
@@ -661,17 +669,26 @@
 static int edge(double x0, double y0, double x1, double y1)
 {
     register double m;
-    double dy, x;
+    double x, d;
     int ystart, ystop;
+    int exp;
 
+    /* tolerance to avoid FPE */
+    d = GRASS_EPSILON;
+    if (y0 != y1) {
+	if (fabs(y0) > fabs(y1))
+	    d = fabs(y0);
+	else
+	    d = fabs(y1);
 
-    /* tolerance to avoid FPE */
-    dy = y0 - y1;
-    if (fabs(dy) < 1e-10)
+	d = frexp(d, &exp);
+	exp -= 53;
+	d = ldexp(d, exp);
+    }
+
+    if (fabs(y0 - y1) < d)
 	return 1;
 
-    m = (x0 - x1) / dy;
-
     if (y0 < y1) {
 	ystart = iceil(y0);
 	ystop = ifloor(y1);
@@ -687,6 +704,7 @@
     if (ystart > ystop)
 	return 1;		/* does not cross center line of row */
 
+    m = (x0 - x1) / (y0 - y1);
     x = m * (ystart - y0) + x0;
     while (ystart <= ystop) {
 	if (!edge_point(x, ystart++))