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

[Raul]

Author: algunenano
Date: 2018-05-21 04:03:01 -0700 (Mon, 21 May 2018)
New Revision: 16580

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

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


Modified: branches/2.3/NEWS
===================================================================
--- branches/2.3/NEWS	2018-05-21 10:59:22 UTC (rev 16579)
+++ branches/2.3/NEWS	2018-05-21 11:03:01 UTC (rev 16580)
@@ -5,6 +5,7 @@
 
   - #4071, ST_ClusterKMeans crash on NULL/EMPTY fixed (Darafei Praliaskouski)
   - #4074, Disable interrupt tests by default (Regina Obe)
+  - #3997, fix bug in lwgeom_median and avoid division by zero (Raúl Marín)
 
 PostGIS 2.3.7
 2018/04/06

Modified: branches/2.3/liblwgeom/lwgeom_median.c
===================================================================
--- branches/2.3/liblwgeom/lwgeom_median.c	2018-05-21 10:59:22 UTC (rev 16579)
+++ branches/2.3/liblwgeom/lwgeom_median.c	2018-05-21 11:03:01 UTC (rev 16580)
@@ -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);