Rewrite code to split a line by a (multi)point to improve robustness
authorSandro Santilli <strk@keybit.net>
Tue, 5 Jan 2016 16:39:55 +0000 (16:39 +0000)
committerSandro Santilli <strk@keybit.net>
Tue, 5 Jan 2016 16:39:55 +0000 (16:39 +0000)
Fixes #3401 including unit and regress test for it.

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

liblwgeom/cunit/cu_split.c
liblwgeom/lwgeom_geos_split.c
topology/test/regress/st_modedgesplit.sql
topology/test/regress/st_modedgesplit_expected

index 285491e9d48fa8d6370b1087b3ea370c7f8331a4..71ef560e90f058688a648b9629cf89d4dd617ab3 100644 (file)
@@ -246,6 +246,32 @@ static void test_lwgeom_split(void)
        lwgeom_free(ret);
        lwgeom_free(geom);
        lwgeom_free(blade);
+
+  /* See #3401 -- robustness issue */
+  geom = lwgeom_from_wkt("LINESTRING(-180 0,0 0)", LW_PARSER_CHECK_NONE);
+  CU_ASSERT(geom != NULL);
+  blade = lwgeom_from_wkt("POINT(-20 0)", LW_PARSER_CHECK_NONE);
+  ret = lwgeom_split(geom, blade);
+  CU_ASSERT(ret != NULL);
+       {
+               LWCOLLECTION *split = lwgeom_as_lwcollection(ret);
+               LWLINE *l1, *l2;
+               POINT2D pt;
+               CU_ASSERT(split != NULL);
+               l1 = lwgeom_as_lwline(split->geoms[0]);
+               CU_ASSERT(l1 != NULL);
+               getPoint2d_p(l1->points, 1, &pt);
+               ASSERT_DOUBLE_EQUAL(pt.x, -20);
+               ASSERT_DOUBLE_EQUAL(pt.y, 0);
+               l2 = lwgeom_as_lwline(split->geoms[1]);
+               CU_ASSERT(l2 != NULL);
+               getPoint2d_p(l2->points, 0, &pt);
+               ASSERT_DOUBLE_EQUAL(pt.x, -20);
+               ASSERT_DOUBLE_EQUAL(pt.y, 0);
+       }
+       lwgeom_free(ret);
+       lwgeom_free(geom);
+       lwgeom_free(blade);
 }
 
 
index 9e520aca7a75355d951beb2cb70e5908a919a1fb..1bf332475e9b5fd98dcbe0d5c48995452fd305b0 100644 (file)
@@ -22,7 +22,8 @@
  *
  **********************************************************************/
 
-
+#include "../postgis_config.h"
+/*#define POSTGIS_DEBUG_LEVEL 4*/
 #include "lwgeom_geos.h"
 #include "liblwgeom_internal.h"
 
@@ -208,11 +209,13 @@ int
 lwline_split_by_point_to(const LWLINE* lwline_in, const LWPOINT* blade_in,
                          LWMLINE* v)
 {
-       double loc, dist;
+       double mindist = -1;
        POINT4D pt, pt_projected;
+       POINT4D p1, p2;
+       POINTARRAY *ipa = lwline_in->points;
        POINTARRAY* pa1;
        POINTARRAY* pa2;
-       double vstol; /* vertex snap tolerance */
+       int i, nsegs, seg = -1;
 
        /* Possible outcomes:
         *
@@ -228,30 +231,69 @@ lwline_split_by_point_to(const LWLINE* lwline_in, const LWPOINT* blade_in,
         */
 
        getPoint4d_p(blade_in->point, 0, &pt);
-       loc = ptarray_locate_point(lwline_in->points, &pt, &dist, &pt_projected);
-
-       /* lwnotice("Location: %g -- Distance: %g", loc, dist); */
 
-       if ( dist > 0 )   /* TODO: accept a tolerance ? */
+       /* Find closest segment */
+       getPoint4d_p(ipa, 0, &p1);
+       nsegs = ipa->npoints - 1;
+       for ( i = 0; i < nsegs; i++ )
        {
-               /* No intersection */
-               return 0;
+               getPoint4d_p(ipa, i+1, &p2);
+               double dist;
+               dist = distance2d_pt_seg((POINT2D*)&pt, (POINT2D*)&p1, (POINT2D*)&p2);
+               LWDEBUGF(4, " Distance of point %g %g to segment %g %g, %g %g: %g", pt.x, pt.y, p1.x, p1.y, p2.x, p2.y, dist);
+               if (i==0 || dist < mindist )
+               {
+                       mindist = dist;
+                       seg=i;
+               }
+               p1 = p2;
        }
 
-       if ( loc == 0 || loc == 1 )
+       LWDEBUGF(3, "Closest segment: %d", seg);
+       LWDEBUGF(3, "mindist: %g", mindist);
+
+       /* No intersection */
+       if ( mindist > 0 ) return 0;
+
+       /* empty or single-point line, intersection on boundary */
+       if ( seg < 0 ) return 1;
+
+       /*
+        * We need to project the
+        * point on the closest segment.
+        */
+       getPoint4d_p(ipa, seg, &p1);
+       getPoint4d_p(ipa, seg+1, &p2);
+       closest_point_on_segment(&pt, &p1, &p2, &pt_projected);
+
+       LWDEBUGF(3, "Projected point:(%g %g), seg:%d, p1:(%g %g), p2:(%g %g)", pt_projected.x, pt_projected.y, seg, p1.x, p1.y, p2.x, p2.y);
+
+       /* When closest point == an endpoint, this is a boundary intersection */
+       if ( ( (seg == nsegs-1) && p4d_same(&pt_projected, &p2) ) ||
+            ( (seg == 0)       && p4d_same(&pt_projected, &p1) ) )
        {
-               /* Intersection is on the boundary */
                return 1;
        }
 
-       /* There is a real intersection, let's get two substrings */
+       /* This is an internal intersection, let's build the two new pointarrays */
 
-       /* Compute vertex snap tolerance based on line length
-        * TODO: take as parameter ? */
-       vstol = ptarray_length_2d(lwline_in->points) / 1e14;
+       pa1 = ptarray_construct_empty(FLAGS_GET_Z(ipa->flags), FLAGS_GET_M(ipa->flags), seg+2);
+       /* TODO: replace with a memcpy ? */
+       for (i=0; i<=seg; ++i)
+       {
+               getPoint4d_p(ipa, i, &p1);
+               ptarray_append_point(pa1, &p1, LW_FALSE);
+       }
+       ptarray_append_point(pa1, &pt_projected, LW_FALSE);
 
-       pa1 = ptarray_substring(lwline_in->points, 0, loc, vstol);
-       pa2 = ptarray_substring(lwline_in->points, loc, 1, vstol);
+       pa2 = ptarray_construct_empty(FLAGS_GET_Z(ipa->flags), FLAGS_GET_M(ipa->flags), ipa->npoints-seg);
+       ptarray_append_point(pa2, &pt_projected, LW_FALSE);
+       /* TODO: replace with a memcpy (if so need to check for duplicated point) ? */
+       for (i=seg+1; i<ipa->npoints; ++i)
+       {
+               getPoint4d_p(ipa, i, &p1);
+               ptarray_append_point(pa2, &p1, LW_FALSE);
+       }
 
        /* NOTE: I've seen empty pointarrays with loc != 0 and loc != 1 */
        if ( pa1->npoints == 0 || pa2->npoints == 0 ) {
index 4684965e8fd9edc1e086291b3fc2fb7a95567e08..537090eea08b292c3e19d2f114afb73318364b4d 100644 (file)
@@ -148,3 +148,13 @@ DROP TABLE t;
 DROP FUNCTION check_changes();
 SELECT DropTopology('city_data');
 DELETE FROM spatial_ref_sys where srid = 4326;
+
+-- See https://trac.osgeo.org/postgis/ticket/3401
+SELECT 't3401.start', CreateTopology('bug3401') > 1;
+SELECT 't3401.N1', ST_AddIsoNode('bug3401', 0, ST_MakePoint(-180, 0));
+SELECT 't3401.N2', ST_AddIsoNode('bug3401', 0, ST_MakePoint(0, 0));
+SELECT 't3401.L1', ST_AddIsoEdge('bug3401', 1, 2, 'LINESTRING(-180 0, 0 0)'::geometry);
+SELECT 't3401.valid_before', * FROM ValidateTopology('bug3401');
+SELECT 't3401.split', ST_ModEdgeSplit('bug3401', 1, ST_MakePoint(-20, 0));
+SELECT 't3401.valid_after', * FROM ValidateTopology('bug3401');
+SELECT 't3401.end', DropTopology('bug3401');
index 263e9c311cb36a232e61c1b4d98f80ed01b50abe..62500b92aa5d1e72b3acf2fe077d688838a9e2ab 100644 (file)
@@ -42,3 +42,9 @@ seq_reset|1|1|1
 robust.1|E1|N3
 robust.2|t|t
 Topology 'city_data' dropped
+t3401.start|t
+t3401.N1|1
+t3401.N2|2
+t3401.L1|1
+t3401.split|3
+t3401.end|Topology 'bug3401' dropped