[postgis-tickets] r16577 - Survive to big max deviation values passed to ST_CurveToLine

[Sandro Santilli]

Author: strk
Date: 2018-05-20 07:46:48 -0700 (Sun, 20 May 2018)
New Revision: 16577

Modified:
   branches/2.4/NEWS
   branches/2.4/liblwgeom/cunit/cu_lwstroke.c
   branches/2.4/liblwgeom/lwstroke.c
Log:
Survive to big max deviation values passed to ST_CurveToLine

When using "max-deviation" tolerance type, passing a tolerance
bigger than twice the radius of any arc resulted in entering
an infinite loop, only limited by availability of RAM.

This commit fixes the bug by being careful in what's fed to
acos()...

Includes a unit test.
Closes #4031 in 2.4 branch (2.4.5dev)

Modified: branches/2.4/NEWS
===================================================================
--- branches/2.4/NEWS	2018-05-18 09:10:03 UTC (rev 16576)
+++ branches/2.4/NEWS	2018-05-20 14:46:48 UTC (rev 16577)
@@ -3,9 +3,12 @@
 
   * Bug Fixes and Enchantments
 
+  - #4031, Survive to big MaxError tolerances passed to ST_CurveToLine
+           (Sandro Santilli)
   - #4071, ST_ClusterKMeans crash on NULL/EMPTY fixed (Darafei Praliaskouski)
   - #4079, ensure St_AsMVTGeom outputs CW oriented polygons (Paul Ramsey)
-  - #4070, use standard interruption error code on GEOS interruptions (Paul Ramsey)
+  - #4070, use standard interruption error code on GEOS interruptions
+           (Paul Ramsey)
   - #3980, delay freeing input until processing complete (lucasvr)
 
 

Modified: branches/2.4/liblwgeom/cunit/cu_lwstroke.c
===================================================================
--- branches/2.4/liblwgeom/cunit/cu_lwstroke.c	2018-05-18 09:10:03 UTC (rev 16576)
+++ branches/2.4/liblwgeom/cunit/cu_lwstroke.c	2018-05-20 14:46:48 UTC (rev 16577)
@@ -224,6 +224,17 @@
 	lwfree(str);
 	lwgeom_free(out);
 
+	/* max deviation bigger than twice the radius
+	 * we really only want to make sure NOT to enter
+	 * an infinite loop here.
+	 * See https://trac.osgeo.org/postgis/ticket/4031
+	 */
+	out = lwcurve_linearize(in, 500, toltype, LW_LINEARIZE_FLAG_SYMMETRIC);
+	str = lwgeom_to_text(out, 2);
+	ASSERT_STRING_EQUAL(str, "LINESTRING(20 50,72 -66)");
+	lwfree(str);
+	lwgeom_free(out);
+
 	lwgeom_free(in);
 
 	/***********************************************************

Modified: branches/2.4/liblwgeom/lwstroke.c
===================================================================
--- branches/2.4/liblwgeom/lwstroke.c	2018-05-18 09:10:03 UTC (rev 16576)
+++ branches/2.4/liblwgeom/lwstroke.c	2018-05-20 14:46:48 UTC (rev 16577)
@@ -206,13 +206,45 @@
 	}}
 	else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION )
 	{{
-		double halfAngle;
+		double halfAngle, maxErr;
 		if ( tol <= 0 )
 		{
 			lwerror("lwarc_linearize: max deviation must be bigger than 0, got %.15g", tol);
 			return -1;
 		}
-		halfAngle = acos( -tol / radius + 1 );
+
+		/*
+		 * Ref: https://en.wikipedia.org/wiki/Sagitta_(geometry)
+		 *
+		 * An arc "sagitta" (distance between middle point of arc and
+		 * middle point of corresponding chord) is defined as:
+		 *
+		 *   sagitta = radius * ( 1 - cos( angle ) );
+		 *
+		 * We want our sagitta to be at most "tolerance" long,
+		 * and we want to find out angle, so we use the inverse
+		 * formula:
+		 *
+		 *   tol = radius * ( 1 - cos( angle ) );
+		 *   1 - cos( angle ) =  tol/radius
+		 *   - cos( angle ) =  tol/radius - 1
+		 *   cos( angle ) =  - tol/radius + 1
+		 *   angle = acos( 1 - tol/radius )
+		 *
+		 * Constraints: 1.0 - tol/radius must be between -1 and 1
+		 * which means tol/radius must be between 0 and 2 times
+		 * the radius, which makes sense as you cannot have a
+		 * sagitta bigger than twice the radius!
+		 *
+		 */
+		maxErr = tol;
+		if ( maxErr > radius * 2 )
+		{
+			maxErr = radius * 2;
+			LWDEBUGF(2, "lwarc_linearize: tolerance %g is too big, "
+			            "using arc-max 2 * radius == %g", tol, maxErr);
+		}
+		halfAngle = acos( 1.0 - maxErr / radius );
 		increment = 2 * halfAngle;
 		LWDEBUGF(2, "lwarc_linearize: maxDiff:%g, radius:%g, halfAngle:%g, increment:%g (%g degrees)", tol, radius, halfAngle, increment, increment*180/M_PI);
 	}}