[postgis-tickets] r14535 - Use recursive snapping to improve predictability

[Sandro Santilli]

Author: strk
Date: 2015-12-31 08:16:53 -0800 (Thu, 31 Dec 2015)
New Revision: 14535

Modified:
   trunk/liblwgeom/lwgeom_topo.c
   trunk/topology/test/regress/topogeo_addlinestring.sql
   trunk/topology/test/regress/topogeo_addlinestring_expected_newsnap
   trunk/topology/test/regress/topogeo_addlinestring_expected_oldsnap
Log:
Use recursive snapping to improve predictability

Fixes geometry-intersects-edge exception when snapping twice
to the same pointset. See #3412.

Includes automated testcase for both old and new geos snap
(3.3.8- and 3.3.9+)

Modified: trunk/liblwgeom/lwgeom_topo.c
===================================================================
--- trunk/liblwgeom/lwgeom_topo.c	2015-12-31 16:15:43 UTC (rev 14534)
+++ trunk/liblwgeom/lwgeom_topo.c	2015-12-31 16:16:53 UTC (rev 14535)
@@ -428,6 +428,43 @@
  *
  ************************************************************************/
 
+static LWGEOM *
+_lwt_toposnap(LWGEOM *src, LWGEOM *tgt, double tol)
+{
+  LWGEOM *tmp = src;
+  LWGEOM *tmp2;
+  int changed;
+  int iterations = 0;
+
+  int maxiterations = lwgeom_count_vertices(tgt);
+
+  /* GEOS snapping can be unstable */
+  /* See https://trac.osgeo.org/geos/ticket/760 */
+  do {
+    LWGEOM *tmp3;
+    tmp2 = lwgeom_snap(tmp, tgt, tol);
+    ++iterations;
+    changed = ( lwgeom_count_vertices(tmp2) != lwgeom_count_vertices(tmp) );
+#if GEOS_NUMERIC_VERSION < 30309
+    /* Up to GEOS-3.3.8, snapping could duplicate points */
+    if ( changed ) {
+      tmp3 = lwgeom_remove_repeated_points( tmp2, 0 );
+      lwgeom_free(tmp2);
+      tmp2 = tmp3;
+      changed = ( lwgeom_count_vertices(tmp2) != lwgeom_count_vertices(tmp) );
+    }
+#endif /* GEOS_NUMERIC_VERSION < 30309 */
+    LWDEBUGF(2, "After iteration %d, geometry changed ? %d (%d vs %d vertices)", iterations, changed, lwgeom_count_vertices(tmp2), lwgeom_count_vertices(tmp));
+    if ( tmp != src ) lwgeom_free(tmp);
+    tmp = tmp2;
+  } while ( changed && iterations <= maxiterations );
+
+  LWDEBUGF(1, "It took %d/%d iterations to properly snap",
+              iterations, maxiterations);
+
+  return tmp;
+}
+
 static void
 _lwt_release_faces(LWT_ISO_FACE *faces, int num_faces)
 {
@@ -5144,7 +5181,7 @@
       -- a projected point internally, so we need another way.
       */
       snaptol = _lwt_minTolerance(prj);
-      snapedge = lwgeom_snap(g, prj, snaptol);
+      snapedge = _lwt_toposnap(g, prj, snaptol);
       snapline = lwgeom_as_lwline(snapedge);
 
       LWDEBUGF(1, "Edge snapped with tolerance %g", snaptol);
@@ -5561,7 +5598,7 @@
       LWDEBUGF(1, "Snapping noded, with srid=%d "
                   "to interesecting edges, with srid=%d",
                   noded->srid, iedges->srid);
-      snapped = lwgeom_snap(noded, iedges, tol);
+      snapped = _lwt_toposnap(noded, iedges, tol);
       lwgeom_free(noded);
       LWDEBUGG(1, snapped, "Snapped");
       LWDEBUGF(1, "Diffing snapped, with srid=%d "
@@ -5634,7 +5671,7 @@
 
       /* TODO: consider snapping once against all elements
        *      (rather than once with edges and once with nodes) */
-      tmp = lwgeom_snap(noded, inodes, tol);
+      tmp = _lwt_toposnap(noded, inodes, tol);
       lwgeom_free(noded);
       noded = tmp;
       LWDEBUGG(1, noded, "Node-snapped");

Modified: trunk/topology/test/regress/topogeo_addlinestring.sql
===================================================================
--- trunk/topology/test/regress/topogeo_addlinestring.sql	2015-12-31 16:15:43 UTC (rev 14534)
+++ trunk/topology/test/regress/topogeo_addlinestring.sql	2015-12-31 16:16:53 UTC (rev 14535)
@@ -276,3 +276,17 @@
  ORDER BY 1;
 SELECT 't3280.end', topology.DropTopology('bug3280');
 
+-- See http://trac.osgeo.org/postgis/ticket/3412
+SELECT 't3412.start', CreateTopology('bug3412', 0, 0.001) > 0;
+SELECT 't3412.L1', TopoGeo_AddLinestring('bug3412',
+'LINESTRING(
+599671.37 4889664.32,
+599665.52 4889680.18,
+599671.37 4889683.4,
+599671.37 4889781.87
+)'
+::geometry, 0);
+SELECT 't3412.L2', TopoGeo_AddLinestring('bug3412',
+'0102000000020000003AB42BBFEE4C22410010C5A997A6524167BB5DBDEE4C224117FE3DA85FA75241'
+::geometry, 0);
+SELECT 't3412.end', DropTopology('bug3412');

Modified: trunk/topology/test/regress/topogeo_addlinestring_expected_newsnap
===================================================================
--- trunk/topology/test/regress/topogeo_addlinestring_expected_newsnap	2015-12-31 16:15:43 UTC (rev 14534)
+++ trunk/topology/test/regress/topogeo_addlinestring_expected_newsnap	2015-12-31 16:16:53 UTC (rev 14535)
@@ -183,3 +183,10 @@
 t3280|L1b4
 t3280|L1b2
 t3280.end|Topology 'bug3280' dropped
+t3412.start|t
+t3412.L1|1
+t3412.L2|2
+t3412.L2|3
+t3412.L2|5
+t3412.L2|4
+t3412.end|Topology 'bug3412' dropped

Modified: trunk/topology/test/regress/topogeo_addlinestring_expected_oldsnap
===================================================================
--- trunk/topology/test/regress/topogeo_addlinestring_expected_oldsnap	2015-12-31 16:15:43 UTC (rev 14534)
+++ trunk/topology/test/regress/topogeo_addlinestring_expected_oldsnap	2015-12-31 16:16:53 UTC (rev 14535)
@@ -182,3 +182,10 @@
 t3280|L1b4
 t3280|L1b2
 t3280.end|Topology 'bug3280' dropped
+t3412.start|t
+t3412.L1|1
+t3412.L2|2
+t3412.L2|3
+t3412.L2|5
+t3412.L2|4
+t3412.end|Topology 'bug3412' dropped