Add a vertex-snap-tolerance parameter to ptarray_substring

This is aimed at improving robustness for ST_Split, which now
uses an hard-coded tolerance of 1e-14 (see #1311)

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

diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index aef35d35a2a92d4df24bbeff0e02a8ceb941d5b9..4d7624597939780c317585f41e188db3526f18c7 100644 (file)
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -818,7 +818,14 @@ extern int ptarray_isclosedz(const POINTARRAY *pa);
 extern void ptarray_longitude_shift(POINTARRAY *pa);
 extern void ptarray_reverse(POINTARRAY *pa);
 extern POINTARRAY* ptarray_flip_coordinates(POINTARRAY *pa);
-extern POINTARRAY *ptarray_substring(POINTARRAY *pa, double d1, double d2);
+
+/**
+ * @d1 start location (distance from start / total distance)
+ * @d2   end location (distance from start / total distance)
+ * @param tolerance snap to vertices at locations < tolerance away from given ones
+ */
+extern POINTARRAY *ptarray_substring(POINTARRAY *pa, double d1, double d2,
+                                                          double tolerance);
 
 
 /**

diff --git a/liblwgeom/lwgeom_geos_split.c b/liblwgeom/lwgeom_geos_split.c
index 33dcabea0d28f39593afec07882247be7c6eb7af..b550a48f757b61bae8cf64e23a2c69c07f6edf80 100644 (file)
--- a/liblwgeom/lwgeom_geos_split.c
+++ b/liblwgeom/lwgeom_geos_split.c
@@ -163,6 +163,7 @@ lwline_split_by_point_to(const LWLINE* lwline_in, const LWPOINT* blade_in,
        POINT2D pt;
        POINTARRAY* pa1;
        POINTARRAY* pa2;
+       double vstol; /* vertex snap tolerance */
 
        /* Possible outcomes:
         *
@@ -193,9 +194,9 @@ lwline_split_by_point_to(const LWLINE* lwline_in, const LWPOINT* blade_in,
        }
 
        /* There is a real intersection, let's get two substrings */
-
-       pa1 = ptarray_substring(lwline_in->points, 0, loc);
-       pa2 = ptarray_substring(lwline_in->points, loc, 1);
+       vstol = 1e-14; /* TODO: take this as parameter ? */
+       pa1 = ptarray_substring(lwline_in->points, 0, loc, vstol);
+       pa2 = ptarray_substring(lwline_in->points, loc, 1, vstol);
 
        /* NOTE: I've seen empty pointarrays with loc != 0 and loc != 1 */
        if ( pa1->npoints == 0 || pa2->npoints == 0 ) {

diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 059d95092208bdc541f20800fd797421c3071e7d..bd2992d9a04846c5590259fdd7b9cb39b259a960 100644 (file)
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -668,7 +668,7 @@ ptarray_force_dims(const POINTARRAY *pa, int hasz, int hasm)
 }
 
 POINTARRAY *
-ptarray_substring(POINTARRAY *ipa, double from, double to)
+ptarray_substring(POINTARRAY *ipa, double from, double to, double tolerance)
 {
        POINTARRAY *dpa;
        POINT4D pt;
@@ -725,13 +725,7 @@ ptarray_substring(POINTARRAY *ipa, double from, double to)
 
                        LWDEBUG(3, " Before start");
 
-                       /*
-                        * Didn't reach the 'from' point,
-                        * nothing to do
-                        */
-                       if ( from > tlength + slength ) goto END;
-
-                       else if ( from == tlength + slength )
+                       if ( fabs ( from - ( tlength + slength ) ) <= tolerance )
                        {
 
                                LWDEBUG(3, "  Second point is our start");
@@ -739,12 +733,12 @@ ptarray_substring(POINTARRAY *ipa, double from, double to)
                                /*
                                 * Second point is our start
                                 */
-                               ptarray_append_point(dpa, &p2, LW_TRUE);
+                               ptarray_append_point(dpa, &p2, LW_FALSE);
                                state=1; /* we're inside now */
                                goto END;
                        }
 
-                       else if ( from == tlength )
+                       else if ( fabs(from - tlength) <= tolerance )
                        {
 
                                LWDEBUG(3, "  First point is our start");
@@ -752,7 +746,7 @@ ptarray_substring(POINTARRAY *ipa, double from, double to)
                                /*
                                 * First point is our start
                                 */
-                               ptarray_append_point(dpa, &p1, LW_TRUE);
+                               ptarray_append_point(dpa, &p1, LW_FALSE);
 
                                /*
                                 * We're inside now, but will check
@@ -761,6 +755,12 @@ ptarray_substring(POINTARRAY *ipa, double from, double to)
                                state=1;
                        }
 
+                       /*
+                        * Didn't reach the 'from' point,
+                        * nothing to do
+                        */
+                       else if ( from > tlength + slength ) goto END;
+
                        else  /* tlength < from < tlength+slength */
                        {
 
@@ -771,9 +771,10 @@ ptarray_substring(POINTARRAY *ipa, double from, double to)
                                 * second point
                                 */
                                dseg = (from - tlength) / slength;
+
                                interpolate_point4d(&p1, &p2, &pt, dseg);
 
-                               ptarray_append_point(dpa, &pt, LW_TRUE);
+                               ptarray_append_point(dpa, &pt, LW_FALSE);
 
                                /*
                                 * We're inside now, but will check
@@ -788,20 +789,10 @@ ptarray_substring(POINTARRAY *ipa, double from, double to)
 
                        LWDEBUG(3, " Inside");
 
-                       /*
-                        * Didn't reach the 'end' point,
-                        * just copy second point
-                        */
-                       if ( to > tlength + slength )
-                       {
-                               ptarray_append_point(dpa, &p2, LW_FALSE);
-                               goto END;
-                       }
-
                        /*
                         * 'to' point is our second point.
                         */
-                       else if ( to == tlength + slength )
+                       if ( fabs(to - ( tlength + slength ) ) <= tolerance )
                        {
 
                                LWDEBUG(3, " Second point is our end");
@@ -814,7 +805,7 @@ ptarray_substring(POINTARRAY *ipa, double from, double to)
                         * 'to' point is our first point.
                         * (should only happen if 'to' is 0)
                         */
-                       else if ( to == tlength )
+                       else if ( fabs(to - tlength) <= tolerance )
                        {
 
                                LWDEBUG(3, " First point is our end");
@@ -824,6 +815,16 @@ ptarray_substring(POINTARRAY *ipa, double from, double to)
                                break; /* substring complete */
                        }
 
+                       /*
+                        * Didn't reach the 'end' point,
+                        * just copy second point
+                        */
+                       else if ( to > tlength + slength )
+                       {
+                               ptarray_append_point(dpa, &p2, LW_FALSE);
+                               goto END;
+                       }
+
                        /*
                         * 'to' point falls on this segment
                         * Interpolate and break.

diff --git a/postgis/lwgeom_functions_analytic.c b/postgis/lwgeom_functions_analytic.c
index 6c6bb6f32ffc5f09b2ad446d64235fa3487eff9e..4354a71956b802eef6cda4e4b535d734f91ac780 100644 (file)
--- a/postgis/lwgeom_functions_analytic.c
+++ b/postgis/lwgeom_functions_analytic.c
@@ -761,7 +761,7 @@ Datum LWGEOM_line_substring(PG_FUNCTION_ARGS)
 
                ipa = iline->points;
 
-               opa = ptarray_substring(ipa, from, to);
+               opa = ptarray_substring(ipa, from, to, 0);
 
                if ( opa->npoints == 1 ) /* Point returned */
                        olwgeom = (LWGEOM *)lwpoint_construct(iline->srid, NULL, opa);
@@ -827,7 +827,7 @@ Datum LWGEOM_line_substring(PG_FUNCTION_ARGS)
                                subto = (to - minprop) / (maxprop - minprop);
 
 
-                       opa = ptarray_substring(subline->points, subfrom, subto);
+                       opa = ptarray_substring(subline->points, subfrom, subto, 0);
                        if ( opa && opa->npoints > 0 )
                        {
                                if ( opa->npoints == 1 ) /* Point returned */