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

[Sandro Santilli]

Author: strk
Date: 2018-04-24 06:26:49 -0700 (Tue, 24 Apr 2018)
New Revision: 16553

Modified:
   trunk/liblwgeom/cunit/cu_lwstroke.c
   trunk/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.
References #4031 for trunk (2.5.0dev) - to be backported

Modified: trunk/liblwgeom/cunit/cu_lwstroke.c
===================================================================
--- trunk/liblwgeom/cunit/cu_lwstroke.c	2018-04-24 11:52:20 UTC (rev 16552)
+++ trunk/liblwgeom/cunit/cu_lwstroke.c	2018-04-24 13:26:49 UTC (rev 16553)
@@ -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: trunk/liblwgeom/lwstroke.c
===================================================================
--- trunk/liblwgeom/lwstroke.c	2018-04-24 11:52:20 UTC (rev 16552)
+++ trunk/liblwgeom/lwstroke.c	2018-04-24 13:26:49 UTC (rev 16553)
@@ -207,13 +207,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);
 	}}