[postgis-tickets] r16366 - lwgeom_median: Avoid division by zero

[Paul Ramsey]

Author: pramsey
Date: 2018-01-26 07:18:40 -0800 (Fri, 26 Jan 2018)
New Revision: 16366

Modified:
   trunk/liblwgeom/cunit/cu_algorithm.c
   trunk/liblwgeom/lwgeom_median.c
Log:
lwgeom_median: Avoid division by zero
(Raúl Marín Rodríguez)
Closes #3997
Closes https://github.com/postgis/postgis/pull/198


Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c	2018-01-26 15:10:58 UTC (rev 16365)
+++ trunk/liblwgeom/cunit/cu_algorithm.c	2018-01-26 15:18:40 UTC (rev 16366)
@@ -1122,8 +1122,9 @@
 	LWPOINT* expected_result = NULL;
 	POINT4D actual_pt;
 	POINT4D expected_pt;
+	const double tolerance = FP_TOLERANCE / 10.0;
 
-	LWPOINT* result = lwgeom_median(g, FP_TOLERANCE / 10.0, iter_count, fail_if_not_converged);
+	LWPOINT* result = lwgeom_median(g, tolerance, iter_count, fail_if_not_converged);
 	int passed = LW_FALSE;
 
 	if (expected != NULL)
@@ -1143,10 +1144,10 @@
 		passed = passed && (lwgeom_has_z((LWGEOM*) expected_result) == lwgeom_has_z((LWGEOM*) result));
 		if (!lwgeom_is_empty((LWGEOM*) result))
 		{
-			passed = passed && FP_EQUALS(actual_pt.x, expected_pt.x);
-			passed = passed && FP_EQUALS(actual_pt.y, expected_pt.y);
-			passed = passed && (!lwgeom_has_z((LWGEOM*) expected_result) || FP_EQUALS(actual_pt.z, expected_pt.z));
-			passed = passed && (!lwgeom_has_m((LWGEOM*) expected_result) || FP_EQUALS(actual_pt.m, expected_pt.m));
+			passed = passed && abs(actual_pt.x - expected_pt.x) < tolerance;
+			passed = passed && abs(actual_pt.y - expected_pt.y) < tolerance;
+			passed = passed && (!lwgeom_has_z((LWGEOM*) expected_result) || abs(actual_pt.z - expected_pt.z) < tolerance);
+			passed = passed && (!lwgeom_has_m((LWGEOM*) expected_result) || abs(actual_pt.m - expected_pt.m) < tolerance);
 		}
 		if (!passed)
 			printf("median_test input %s (parsed %s) expected %s got %s\n", input, lwgeom_to_ewkt(g), lwgeom_to_ewkt((LWGEOM*) expected_result), lwgeom_to_ewkt((LWGEOM*) result));
@@ -1229,6 +1230,19 @@
 	/* Unsupported geometry type */
 	do_median_test("POLYGON((1 0,0 1,1 2,2 1,1 0))", NULL, LW_TRUE, 1000);
 	do_median_test("POLYGON((1 0,0 1,1 2,2 1,1 0))", NULL, LW_FALSE, 1000);
+
+	/* Intermediate point included in the set */
+	do_median_test("MULTIPOINT ZM ("
+			"(0 0 20000 0.5),"
+			"(0 0 59000 0.5),"
+			"(0 48000 -20000 1.3),"
+			"(0 -48000 -20000 1.3),"
+			"(0 -3000 -3472.22222222222262644208967685699462890625 1),"
+			"(0 3000 3472.22222222222262644208967685699462890625 1),"
+			"(0 0 -1644.736842105263121993630193173885345458984375 1),"
+			"(0 0 1644.736842105263121993630193173885345458984375 1)"
+			")",
+		"POINT (0 0 0)", LW_TRUE, 296);
 }
 
 static void test_point_density(void)

Modified: trunk/liblwgeom/lwgeom_median.c
===================================================================
--- trunk/liblwgeom/lwgeom_median.c	2018-01-26 15:10:58 UTC (rev 16365)
+++ trunk/liblwgeom/lwgeom_median.c	2018-01-26 15:18:40 UTC (rev 16366)
@@ -94,7 +94,7 @@
 			double dx = 0;
 			double dy = 0;
 			double dz = 0;
-			double r_inv;
+			double d_sqr;
 			for (i = 0; i < npoints; i++)
 			{
 				if (distances[i] > DBL_EPSILON)
@@ -105,11 +105,15 @@
 				}
 			}
 
-			r_inv = 1.0 / sqrt ( dx*dx + dy*dy + dz*dz );
-
-			next.x = FP_MAX(0, 1.0 - r_inv)*next.x + FP_MIN(1.0, r_inv)*curr->x;
-			next.y = FP_MAX(0, 1.0 - r_inv)*next.y + FP_MIN(1.0, r_inv)*curr->y;
-			next.z = FP_MAX(0, 1.0 - r_inv)*next.z + FP_MIN(1.0, r_inv)*curr->z;
+			d_sqr = sqrt(dx*dx + dy*dy + dz*dz);
+			/* Avoid division by zero if the intermediate point is the median */
+			if (d_sqr > DBL_EPSILON)
+			{
+				double r_inv = FP_MAX(0, 1.0 / d_sqr);
+				next.x = (1.0 - r_inv)*next.x + r_inv*curr->x;
+				next.y = (1.0 - r_inv)*next.y + r_inv*curr->y;
+				next.z = (1.0 - r_inv)*next.z + r_inv*curr->z;
+			}
 		}
 
 		delta = distance3d_pt_pt(curr, &next);