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

diff --git a/liblwgeom/lwsegmentize.c b/liblwgeom/lwsegmentize.c
index 61cf59de587f2f124ce3492e94a1261868b4a0dd..660a3867a4a2c8baac22d4f62a1e18d2f5b2216e 100644 (file)
--- 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 6a26d45cd8ab98a01a8b9a56966b25c37f771937..adc98a4999d85900dc73b203c7b42265ad3c4c19 100644 (file)
--- 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 6d00b48d921dff6a652c9d30c20b93a35349f1f5..cc3ff0162f2b8b6bf8ae555cec4141efb1ef5b0b 100644 (file)
--- 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))