[postgis-tickets] r16579 - Fix bug in lwgeom_median and avoid division by zero

[Raul]

Author: algunenano
Date: 2018-05-21 03:59:22 -0700 (Mon, 21 May 2018)
New Revision: 16579

Modified:
   branches/2.4/NEWS
   branches/2.4/liblwgeom/lwgeom_median.c
Log:
Fix bug in lwgeom_median and avoid division by zero

References #3997
Closes https://github.com/postgis/postgis/pull/244


Modified: branches/2.4/NEWS
===================================================================
--- branches/2.4/NEWS	2018-05-21 10:41:32 UTC (rev 16578)
+++ branches/2.4/NEWS	2018-05-21 10:59:22 UTC (rev 16579)
@@ -11,6 +11,7 @@
            (Paul Ramsey)
   - #3980, delay freeing input until processing complete (lucasvr)
   - #4090, PG 11 support (Paul Ramsey, Raúl Marín)
+  - #3997, fix bug in lwgeom_median and avoid division by zero (Raúl Marín)
 
 
 PostGIS 2.4.4

Modified: branches/2.4/liblwgeom/lwgeom_median.c
===================================================================
--- branches/2.4/liblwgeom/lwgeom_median.c	2018-05-21 10:41:32 UTC (rev 16578)
+++ branches/2.4/liblwgeom/lwgeom_median.c	2018-05-21 10:59:22 UTC (rev 16579)
@@ -85,8 +85,8 @@
 		double dx = 0;
 		double dy = 0;
 		double dz = 0;
-		double r_inv;
-		POINT3D alt;
+		double d_sqr;
+
 		for (i = 0; i < npoints; i++)
 		{
 			if (distances[i] > 0)
@@ -93,17 +93,20 @@
 			{
 				dx += (points[i].x - curr->x) / distances[i];
 				dy += (points[i].y - curr->y) / distances[i];
-				dz += (points[i].y - curr->z) / distances[i];
+				dz += (points[i].z - curr->z) / distances[i];
 			}
 		}
 
-		r_inv = 1.0 / sqrt ( dx*dx + dy*dy + dz*dz );
+		d_sqr = sqrt (dx*dx + dy*dy + dz*dz);
 
-		alt.x = FP_MAX(0, 1.0 - r_inv)*next.x + FP_MIN(1.0, r_inv)*curr->x;
-		alt.y = FP_MAX(0, 1.0 - r_inv)*next.y + FP_MIN(1.0, r_inv)*curr->y;
-		alt.z = FP_MAX(0, 1.0 - r_inv)*next.z + FP_MIN(1.0, r_inv)*curr->z;
+		if (d_sqr > DBL_EPSILON)
+		{
+			double r_inv = 1.0 / d_sqr;
 
-		next = alt;
+			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;
+		}
 	}
 
 	delta = distance3d_pt_pt(curr, &next);