Fix angle reminder computation

Closes #3772 again (see comment:6)
Add specific testcase

git-svn-id: http://svn.osgeo.org/postgis/trunk@15465 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/cunit/cu_lwstroke.c b/liblwgeom/cunit/cu_lwstroke.c
index 48d9850f60df5f5996fc764d081c9e7a332e8f88..5d73d9cd9e255daf4d9d383c816cb460b1635eac 100644 (file)
--- a/liblwgeom/cunit/cu_lwstroke.c
+++ b/liblwgeom/cunit/cu_lwstroke.c
@@ -20,6 +20,8 @@
 #include "liblwgeom_internal.h"
 #include "cu_tester.h"
 
+/* #define SKIP_TEST_RETAIN_ANGLE 1 */
+
 
 static LWGEOM* lwgeom_from_text(const char *str)
 {
@@ -108,21 +110,25 @@ static void test_lwcurve_linearize(void)
        ASSERT_STRING_EQUAL(str, "LINESTRING(30 70,74 96,126 96,170 70,196 26,200 0)");
        lwfree(str);
        lwgeom_free(out);
+
        /* 3 segment per quadrant - symmetric */
        out = lwcurve_linearize(in, 3, toltype, LW_LINEARIZE_FLAG_SYMMETRIC);
        str = lwgeom_to_text(out, 2);
        ASSERT_STRING_EQUAL(str, "LINESTRING(30 70,70 96,116 98,158 80,190 46,200 0)");
        lwfree(str);
        lwgeom_free(out);
+
+#ifndef SKIP_TEST_RETAIN_ANGLE
        /* 3 segment per quadrant - symmetric/retain_angle */
        out = lwcurve_linearize(in, 3, toltype,
                LW_LINEARIZE_FLAG_SYMMETRIC |
                LW_LINEARIZE_FLAG_RETAIN_ANGLE
        );
        str = lwgeom_to_text(out, 2);
-       ASSERT_STRING_EQUAL(str, "LINESTRING(30 70,62 92,114 100,160 80,192 38,200 0)");
+       ASSERT_STRING_EQUAL(str, "LINESTRING(30 70,40 80,86 100,138 92,180 60,200 14,200 0)");
        lwfree(str);
        lwgeom_free(out);
+#endif /* SKIP_TEST_RETAIN_ANGLE */
 
        lwgeom_free(in);
 
@@ -170,13 +176,34 @@ static void test_lwcurve_linearize(void)
        ASSERT_STRING_EQUAL(str, "LINESTRING(0 0,50 86,150 86,200 0)");
        lwfree(str);
        lwgeom_free(out);
+
+#ifndef SKIP_TEST_RETAIN_ANGLE
        /* Maximum of 20 units of difference, symmetric/retain angle */
        out = lwcurve_linearize(in, 20, toltype,
                LW_LINEARIZE_FLAG_SYMMETRIC |
                LW_LINEARIZE_FLAG_RETAIN_ANGLE
        );
        str = lwgeom_to_text(out, 2);
-       ASSERT_STRING_EQUAL(str, "LINESTRING(0 0,40 80,160 80,200 0)");
+       ASSERT_STRING_EQUAL(str, "LINESTRING(0 0,4 28,100 100,196 28,200 0)");
+       lwfree(str);
+       lwgeom_free(out);
+#endif /* SKIP_TEST_RETAIN_ANGLE */
+
+       lwgeom_free(in);
+
+       /*
+        * Clockwise ~90 degrees south-west to north-west quadrants
+        * starting at ~22 degrees west of vertical line
+        *
+        * See https://trac.osgeo.org/postgis/ticket/3772
+        */
+       toltype = LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION;
+       in = lwgeom_from_text("CIRCULARSTRING(71.96 -65.64,22.2 -18.52,20 50)");
+
+       /* 4 units of max deviation - symmetric */
+       out = lwcurve_linearize(in, 4, toltype, LW_LINEARIZE_FLAG_SYMMETRIC);
+       str = lwgeom_to_text(out, 2);
+       ASSERT_STRING_EQUAL(str, "LINESTRING(72 -66,34 -38,16 4,20 50)");
        lwfree(str);
        lwgeom_free(out);
 
@@ -227,15 +254,18 @@ static void test_lwcurve_linearize(void)
        ASSERT_STRING_EQUAL(str, "LINESTRING(0 0,50 86,150 86,200 0)");
        lwfree(str);
        lwgeom_free(out);
+
+#ifndef SKIP_TEST_RETAIN_ANGLE
        /* Maximum of 70 degrees, symmetric/retain angle */
        out = lwcurve_linearize(in, 70 * M_PI / 180, toltype,
                LW_LINEARIZE_FLAG_SYMMETRIC |
                LW_LINEARIZE_FLAG_RETAIN_ANGLE
        );
        str = lwgeom_to_text(out, 2);
-       ASSERT_STRING_EQUAL(str, "LINESTRING(0 0,42 82,158 82,200 0)");
+       ASSERT_STRING_EQUAL(str, "LINESTRING(0 0,6 34,100 100,194 34,200 0)");
        lwfree(str);
        lwgeom_free(out);
+#endif /* SKIP_TEST_RETAIN_ANGLE */
 
        lwgeom_free(in);
 

diff --git a/liblwgeom/lwstroke.c b/liblwgeom/lwstroke.c
index b0101300dddcf3230a7b2a6d079884d137542290..60739184e668c5333c86414f76334c4a84216e18 100644 (file)
--- a/liblwgeom/lwstroke.c
+++ b/liblwgeom/lwstroke.c
@@ -31,7 +31,7 @@
 
 #include "liblwgeom_internal.h"
 
-/* #define POSTGIS_DEBUG_LEVEL 4 */
+/*#define POSTGIS_DEBUG_LEVEL 3*/
 
 #include "lwgeom_log.h"
 
@@ -147,6 +147,7 @@ lwarc_linearize(POINTARRAY *to,
        LWDEBUG(2, "lwarc_linearize called.");
 
        radius = lw_arc_center(t1, t2, t3, &center);
+       LWDEBUGF(2, " center is POINT(%.15g %.15g) - radius:%g", center.x, center.y, radius);
        p2_side = lw_segment_side(t1, t3, t2);
 
        /* Matched start/end points imply circle */
@@ -179,7 +180,7 @@ lwarc_linearize(POINTARRAY *to,
                        return -1;
                }
                increment = fabs(M_PI_2 / perQuad);
-               LWDEBUGF(2, "lwarc_linearize: perQuad:%d, increment:%g", perQuad, increment);
+               LWDEBUGF(2, "lwarc_linearize: perQuad:%d, increment:%g (%g degrees)", perQuad, increment, increment*180/M_PI);
 
        }}
        else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_DEVIATION )
@@ -192,7 +193,7 @@ lwarc_linearize(POINTARRAY *to,
                }
                halfAngle = acos( -tol / radius + 1 );
                increment = 2 * halfAngle;
-               LWDEBUGF(2, "lwarc_linearize: maxDiff:%g, radius:%g, halfAngle:%g, increment:%g", tol, radius, halfAngle, increment);
+               LWDEBUGF(2, "lwarc_linearize: maxDiff:%g, radius:%g, halfAngle:%g, increment:%g (%g degrees)", tol, radius, halfAngle, increment, increment*180/M_PI);
        }}
        else if ( tolerance_type == LW_LINEARIZE_TOLERANCE_TYPE_MAX_ANGLE )
        {
@@ -214,42 +215,47 @@ lwarc_linearize(POINTARRAY *to,
        a2 = atan2(p2->y - center.y, p2->x - center.x);
        a3 = atan2(p3->y - center.y, p3->x - center.x);
 
+       LWDEBUGF(2, "lwarc_linearize A1:%g A2:%g A3:%g",
+               a1*180/M_PI, a2*180/M_PI, a3*180/M_PI);
+
        if ( flags & LW_LINEARIZE_FLAG_SYMMETRIC )
-       {
+       {{
+               /* Calculate total arc angle, in radians */
+               double angle = a1 - a3;
+               if ( angle < 0 ) angle += M_PI * 2;
+               LWDEBUGF(2, "lwarc_linearize SYMMETRIC requested - total angle %g deg",
+                                angle * 180 / M_PI);
                if ( flags & LW_LINEARIZE_FLAG_RETAIN_ANGLE )
                {{
-                       /* Total arc angle, in radians */
-                       double angle = a3 - a1;
                        /* Number of steps */
-                       int steps = floor(angle / increment);
+                       int steps = trunc(angle / increment);
                        /* Angle reminder */
                        double angle_reminder = angle - ( increment * steps );
                        angle_shift = angle_reminder / 2.0;
 
-                       LWDEBUGF(2, "lwarc_linearize SYMMETRIC/RETAIN_ANGLE operation requested - "
-                                "A:%g - LINESTRING(%g %g,%g %g,%g %g) - R:%g",
-                                angle, p1->x, p1->y, center.x, center.y, p3->x, p3->y,
-                                angle_reminder);
+                       LWDEBUGF(2, "lwarc_linearize RETAIN_ANGLE operation requested - "
+                                "total angle %g, steps %d, increment %g, reminder %g",
+                                angle * 180 / M_PI, steps, increment * 180 / M_PI,
+                                angle_reminder * 180 / M_PI);
                }}
                else
                {{
-                       /* Total arc angle, in radians */
-                       double angle = fabs(a3 - a1);
                        /* Number of segments in output */
                        int segs = ceil(angle / increment);
                        /* Tweak increment to be regular for all the arc */
                        increment = angle/segs;
 
                        LWDEBUGF(2, "lwarc_linearize SYMMETRIC operation requested - "
-                                                       "A:%g - LINESTRING(%g %g,%g %g,%g %g) - S:%d - I:%g",
-                                                       angle, p1->x, p1->y, center.x, center.y, p3->x, p3->y,
-                                                       segs, increment);
+                                                       "total angle %g degrees - LINESTRING(%g %g,%g %g,%g %g) - S:%d -   I:%g",
+                                                       angle*180/M_PI, p1->x, p1->y, center.x, center.y, p3->x, p3->y,
+                                                       segs, increment*180/M_PI);
                }}
-       }
+       }}
 
        /* p2 on left side => clockwise sweep */
        if ( clockwise )
        {
+               LWDEBUG(2, " Clockwise sweep");
                increment *= -1;
                angle_shift *= -1;
                /* Adjust a3 down so we can decrement from a1 to a3 cleanly */
@@ -261,6 +267,7 @@ lwarc_linearize(POINTARRAY *to,
        /* p2 on right side => counter-clockwise sweep */
        else
        {
+               LWDEBUG(2, " Counterclockwise sweep");
                /* Adjust a3 up so we can increment from a1 to a3 cleanly */
                if ( a3 < a1 )
                        a3 += 2.0 * M_PI;
@@ -277,11 +284,18 @@ lwarc_linearize(POINTARRAY *to,
                clockwise = LW_FALSE;
        }
 
+       LWDEBUGF(2, "lwarc_linearize angle_shift:%g, increment:%g",
+               angle_shift * 180/M_PI, increment * 180/M_PI);
+
        /* Sweep from a1 to a3 */
        ptarray_append_point(pa, p1, LW_FALSE);
        ++points_added;
-       for ( angle = a1 + increment - angle_shift; clockwise ? angle > a3 : angle < a3; angle += increment )
+       if ( angle_shift ) angle_shift -= increment;
+       LWDEBUGF(2, "a1:%g (%g deg), a3:%g (%g deg), inc:%g, shi:%g, cw:%d",
+               a1, a1 * 180 / M_PI, a3, a3 * 180 / M_PI, increment, angle_shift, clockwise);
+       for ( angle = a1 + increment + angle_shift; clockwise ? angle > a3 : angle < a3; angle += increment )
        {
+               LWDEBUGF(2, " SA: %g ( %g deg )", angle, angle*180/M_PI);
                pt.x = center.x + radius * cos(angle);
                pt.y = center.y + radius * sin(angle);
                pt.z = interpolate_arc(angle, a1, a2, a3, p1->z, p2->z, p3->z);

diff --git a/regress/curvetoline_expected b/regress/curvetoline_expected
index 0b81e38e2306ba6f31e4198816cf7327fc484f82..115897e9b29f235baceaeacff5ed00d7ba4066a3 100644 (file)
--- a/regress/curvetoline_expected
+++ b/regress/curvetoline_expected
@@ -3,6 +3,6 @@ semicircle2|LINESTRING(0 0,72 -96,184 -54,200 0)
 semicircle2.sym|LINESTRING(0 0,50 -86,150 -86,200 0)
 semicircle3|LINESTRING(0 0,24 -64,82 -98,150 -86,194 -34,200 0)
 semicircle3.sym|LINESTRING(0 0,20 -58,70 -96,130 -96,180 -58,200 0)
-semicircle3.sym.ret|LINESTRING(0 0,14 -50,66 -94,134 -94,186 -50,200 0)
+semicircle3.sym.ret|LINESTRING(0 0,2 -18,36 -76,100 -100,164 -76,198 -18,200 0)
 multiarc1|LINESTRING(0 0,30 -70,100 -100,170 -70,200 0,258 142,400 200,542 142,600 0)
 multiarc1.maxerr20.sym|LINESTRING(0 0,50 -86,150 -86,200 0,258 142,400 200,542 142,600 0)