[postgis-tickets] r14549 - Rewrite code to split a line by a (multi)point to improve robustness
Sandro Santilli
strk at keybit.net
Tue Jan 5 08:27:58 PST 2016
Author: strk
Date: 2016-01-05 08:27:58 -0800 (Tue, 05 Jan 2016)
New Revision: 14549
Modified:
branches/2.2/NEWS
branches/2.2/liblwgeom/cunit/cu_split.c
branches/2.2/liblwgeom/lwgeom_geos_split.c
branches/2.2/topology/test/regress/st_modedgesplit.sql
branches/2.2/topology/test/regress/st_modedgesplit_expected
Log:
Rewrite code to split a line by a (multi)point to improve robustness
References #3401 for 2.2 branch.
Includes unit and regress test.
Modified: branches/2.2/NEWS
===================================================================
--- branches/2.2/NEWS 2016-01-04 19:52:21 UTC (rev 14548)
+++ branches/2.2/NEWS 2016-01-05 16:27:58 UTC (rev 14549)
@@ -22,6 +22,7 @@
- #3390, Compilation under Alpine Linux 3.2
gives an error when compiling the postgis and postgis_topology extension
- #3393, ST_Area NaN for some polygons
+ - #3401, Improve ST_Split robustness on 32bit systems
- #3404, ST_ClusterWithin crashes backend
- #3407, Fix crash on splitting a face or an edge
defining multiple TopoGeometry objects
Modified: branches/2.2/liblwgeom/cunit/cu_split.c
===================================================================
--- branches/2.2/liblwgeom/cunit/cu_split.c 2016-01-04 19:52:21 UTC (rev 14548)
+++ branches/2.2/liblwgeom/cunit/cu_split.c 2016-01-05 16:27:58 UTC (rev 14549)
@@ -249,6 +249,32 @@
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);
}
Modified: branches/2.2/liblwgeom/lwgeom_geos_split.c
===================================================================
--- branches/2.2/liblwgeom/lwgeom_geos_split.c 2016-01-04 19:52:21 UTC (rev 14548)
+++ branches/2.2/liblwgeom/lwgeom_geos_split.c 2016-01-05 16:27:58 UTC (rev 14549)
@@ -35,6 +35,8 @@
*
**********************************************************************/
+#include "../postgis_config.h"
+/*#define POSTGIS_DEBUG_LEVEL 4*/
#include "lwgeom_geos.h"
#include "liblwgeom_internal.h"
@@ -220,11 +222,13 @@
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:
*
@@ -240,30 +244,69 @@
*/
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 ) {
Modified: branches/2.2/topology/test/regress/st_modedgesplit.sql
===================================================================
--- branches/2.2/topology/test/regress/st_modedgesplit.sql 2016-01-04 19:52:21 UTC (rev 14548)
+++ branches/2.2/topology/test/regress/st_modedgesplit.sql 2016-01-05 16:27:58 UTC (rev 14549)
@@ -148,3 +148,13 @@
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');
Modified: branches/2.2/topology/test/regress/st_modedgesplit_expected
===================================================================
--- branches/2.2/topology/test/regress/st_modedgesplit_expected 2016-01-04 19:52:21 UTC (rev 14548)
+++ branches/2.2/topology/test/regress/st_modedgesplit_expected 2016-01-05 16:27:58 UTC (rev 14549)
@@ -42,3 +42,9 @@
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
More information about the postgis-tickets
mailing list