From 8ccff1d3dd47a64d860474ba6c34f5ebd80c8d73 Mon Sep 17 00:00:00 2001
From: Sandro Santilli <strk@keybit.net>
Date: Mon, 12 Aug 2013 18:11:36 +0000
Subject: [PATCH] Require all arc edges to form the same angle (#2423)

Note: gives expected result for #183

git-svn-id: http://svn.osgeo.org/postgis/trunk@11778 b70326c6-7e19-0410-871a-916f4a2858ee
---
 liblwgeom/lwsegmentize.c | 43 ++++++++++++++++++++++++++++++++++++++--
 regress/tickets.sql      |  4 ++++
 regress/tickets_expected |  3 ++-
 3 files changed, 47 insertions(+), 3 deletions(-)

diff --git a/liblwgeom/lwsegmentize.c b/liblwgeom/lwsegmentize.c
index 61cf59de5..660a3867a 100644
--- a/liblwgeom/lwsegmentize.c
+++ b/liblwgeom/lwsegmentize.c
@@ -467,6 +467,28 @@ lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad)
 	return ogeom;
 }
 
+/**
+ * Return ABC angle in radians
+ * TODO: move to lwalgorithm
+ */
+static double
+lw_arc_angle(const POINT2D *a, const POINT2D *b, const POINT2D *c)
+{
+  POINT2D ab, cb;
+
+  ab.x = b->x - a->x;
+  ab.y = b->y - a->y;
+
+  cb.x = b->x - c->x;
+  cb.y = b->y - c->y;
+
+  double dot = (ab.x * cb.x + ab.y * cb.y); /* dot product */
+  double cross = (ab.x * cb.y - ab.y * cb.x); /* cross product */
+
+  double alpha = atan2(cross, dot);
+
+  return alpha;
+}
 
 /**
 * Returns LW_TRUE if b is on the arc formed by a1/a2/a3, but not within
@@ -495,7 +517,17 @@ static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D
 	{
 		int a2_side = lw_segment_side(t1, t3, t2);
 		int b_side  = lw_segment_side(t1, t3, tb);
-		
+		double angle1 = lw_arc_angle(t1, t2, t3);
+		double angle2 = lw_arc_angle(t2, t3, tb);
+
+		/* Is the angle similar to the previous one ? */
+		diff = fabs(angle1 - angle2);
+		LWDEBUGF(4, " angle1: %g, angle2: %g, diff:%g", angle1, angle2, diff);
+		if ( diff > EPSILON_SQLMM ) 
+		{
+			return LW_FALSE;
+		}
+
 		/* Is the point b on the same side of a1/a3 as the mid-point a2 is? */
 		/* If not, it's in the unbounded part of the circle, so it continues the arc, return true. */
 		if ( b_side != a2_side )
@@ -550,6 +582,7 @@ pta_desegmentize(POINTARRAY *points, int type, int srid)
 {
 	int i = 0, j, k;
 	POINT4D a1, a2, a3, b;
+	POINT4D first;
 	char *edges_in_arcs;
 	int found_arc = LW_FALSE;
 	int current_arc = 1;
@@ -592,6 +625,8 @@ pta_desegmentize(POINTARRAY *points, int type, int srid)
 		getPoint4d_p(points, i  , &a1);
 		getPoint4d_p(points, i+1, &a2);
 		getPoint4d_p(points, i+2, &a3);
+		memcpy(&first, &a1, sizeof(POINT4D));
+
 		for( j = i+3; j < num_edges+1; j++ )
 		{
 			LWDEBUGF(4, "i=%d, j=%d", i, j);
@@ -612,6 +647,10 @@ pta_desegmentize(POINTARRAY *points, int type, int srid)
 				current_arc++;
 				break;
 			}
+
+			memcpy(&a1, &a2, sizeof(POINT4D));
+			memcpy(&a2, &a3, sizeof(POINT4D));
+			memcpy(&a3,  &b, sizeof(POINT4D));
 		}
 		/* Jump past all the edges that were added to the arc */
 		if ( found_arc )
@@ -623,7 +662,7 @@ pta_desegmentize(POINTARRAY *points, int type, int srid)
 			arc_edges = j - 1 - i;
 			num_quadrants = 1; /* silly guess, TODO: compute */
 			LWDEBUGF(4, "arc defined by %d edges found", arc_edges);
-			if ( a1.x == b.x && a1.y == b.y ) {
+			if ( first.x == b.x && first.y == b.y ) {
 				LWDEBUG(4, "arc is a circle");
 				num_quadrants = 4;
 			}
diff --git a/regress/tickets.sql b/regress/tickets.sql
index 6a26d45cd..adc98a499 100644
--- a/regress/tickets.sql
+++ b/regress/tickets.sql
@@ -833,5 +833,9 @@ SELECT '#2412', ST_AsText(ST_LineToCurve('LINESTRING(0 0,10 0,20 0)'));
 
 SELECT '#2420', ST_AsText(ST_LineToCurve('LINESTRING(0 0,10 0,10 10,0 10,0 0)'));
 
+SELECT '#2423', ST_AsText(ST_SnapToGrid(ST_CurveToLine(ST_LineToCurve(
+  ST_Intersection(ST_Buffer(ST_Point(0,0),10),ST_MakeEnvelope(-10,0,10,10))
+), 4), 1e-5));
+
 -- Clean up
 DELETE FROM spatial_ref_sys;
diff --git a/regress/tickets_expected b/regress/tickets_expected
index 6d00b48d9..cc3ff0162 100644
--- a/regress/tickets_expected
+++ b/regress/tickets_expected
@@ -40,7 +40,7 @@ NOTICE:  No points or linestrings in input array
 #179a|
 NOTICE:  No points or linestrings in input array
 #179b|
-#183|CIRCULARSTRING(0 0,0.5 1.2071067812,0 1)
+#183|COMPOUNDCURVE(CIRCULARSTRING(0 0,0.4653039147 1.2062550401,1 0),(1 0,0 1))
 #210a|
 NOTICE:  No points or linestrings in input array
 #210b|
@@ -248,3 +248,4 @@ ERROR:  invalid GML representation
 #2415.2|MULTISURFACE(CURVEPOLYGON(CIRCULARSTRING(10 0,15 1,20 0,18 5,20 10,10 10,10 0)))
 #2412|LINESTRING(0 0,10 0,20 0)
 #2420|LINESTRING(0 0,10 0,10 10,0 10,0 0)
+#2423|POLYGON((-10 0,-9.2388 3.82683,-7.07107 7.07107,-3.82683 9.2388,0 10,3.82683 9.2388,7.07107 7.07107,9.2388 3.82683,10 0,-10 0))
-- 
2.50.1