From: Sandro Santilli <strk@keybit.net>
Date: Tue, 5 Jan 2016 16:39:55 +0000 (+0000)
Subject: Rewrite code to split a line by a (multi)point to improve robustness
X-Git-Tag: 2.3.0beta1~300
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=e354c14cfdac8940d99c1f620e4a39ee3feecddf;p=postgis

Rewrite code to split a line by a (multi)point to improve robustness

Fixes #3401 including unit and regress test for it.

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

diff --git a/liblwgeom/cunit/cu_split.c b/liblwgeom/cunit/cu_split.c
index 285491e9d..71ef560e9 100644
--- a/liblwgeom/cunit/cu_split.c
+++ b/liblwgeom/cunit/cu_split.c
@@ -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);
 }
 
 
diff --git a/liblwgeom/lwgeom_geos_split.c b/liblwgeom/lwgeom_geos_split.c
index 9e520aca7..1bf332475 100644
--- a/liblwgeom/lwgeom_geos_split.c
+++ b/liblwgeom/lwgeom_geos_split.c
@@ -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 ) {
diff --git a/topology/test/regress/st_modedgesplit.sql b/topology/test/regress/st_modedgesplit.sql
index 4684965e8..537090eea 100644
--- a/topology/test/regress/st_modedgesplit.sql
+++ b/topology/test/regress/st_modedgesplit.sql
@@ -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');
diff --git a/topology/test/regress/st_modedgesplit_expected b/topology/test/regress/st_modedgesplit_expected
index 263e9c311..62500b92a 100644
--- a/topology/test/regress/st_modedgesplit_expected
+++ b/topology/test/regress/st_modedgesplit_expected
@@ -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