TopoGeo_addPoint: use a more functional tolerance when snapping (#1613)

All of this looks like magic but it isn't.
I actually wonder if ST_ModEdgeSplit and ST_ModEdgesSplit and ST_Split
itself should do this internally, and if in doing so we wouldn't need
to do it from higher levels. It doesn't indeed feel comfortable to do
all this noise on such an high level.

Anyway this commit adds a now-passing regression test for the topology
building issue and that's A Good Thing.

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

diff --git a/topology/sql/populate.sql.in.c b/topology/sql/populate.sql.in.c
index 6146c9d3bd91037dc52549bf9502ac0b9795a1c9..3d78ec3fa3cae0fc88c3d9d9d4391bcf4b793c99 100644 (file)
--- a/topology/sql/populate.sql.in.c
+++ b/topology/sql/populate.sql.in.c
@@ -643,6 +643,7 @@ DECLARE
   sql text;
   prj GEOMETRY;
   snapedge GEOMETRY;
+  snaptol FLOAT8;
 BEGIN
 
   -- 0. Check arguments
@@ -694,11 +695,29 @@ BEGIN
 #ifdef POSTGIS_TOPOLOGY_DEBUG
       RAISE DEBUG ' Snapping edge to contain closest point';
 #endif
-      -- Snap the edge geometry to the projected point
-      -- The tolerance is an arbitrary number.
-      -- How much would be enough to ensure any projected point is within
-      -- that distance ? ST_Distance() may very well return 0.
-      snapedge := ST_Snap(rec.geom, prj, 1e-12); -- arbitrarely low number
+      -- The tolerance must be big enough for snapping to happen
+      -- and small enough to snap only to the projected point.
+      -- Unfortunately ST_Distance returns 0 because it also uses
+      -- a projected point internally, so we need another way.
+      --
+      -- A pragmatic test conducted using algoritm shown here:
+      -- http://stackoverflow.com/questions/7408407/generate-next-largest-or-smallest-representable-floating-point-number-without-bi
+      -- showed the "tolerance" growing by an order of magnitude proportionally
+      -- with the order of magnitude of the input, starting with something like
+      -- 3.5527136788005009294e-15 for the starting value of 9.0
+      --
+      -- TODO: make this simpler
+      --
+      snaptol := 3.6 * power(10,  - ( 15 - log(10.0, coalesce(nullif(greatest(abs(ST_X(prj)), abs(ST_Y(prj)))::numeric, 0), 1)) ));
+#ifdef POSTGIS_TOPOLOGY_DEBUG
+      RAISE DEBUG 'Tolerance for snapping to point % = %', ST_AsText(prj), snaptol;
+#endif
+      snapedge := ST_Snap(rec.geom, prj, snaptol);
+#ifdef POSTGIS_TOPOLOGY_DEBUG
+      IF NOT ST_Contains(snapedge, prj) THEN -- or if equal ?
+        RAISE WARNING 'Edge within % distance from node still does not contain the node after snapping to it with tolerance %', tolerance, snaptol;
+      END IF;
+#endif
       PERFORM ST_ChangeEdgeGeom(atopology, rec.edge_id, snapedge);
     END IF;
     id := topology.ST_ModEdgeSplit(atopology, rec.edge_id, prj);

diff --git a/topology/test/regress/topogeo_addlinestring.sql b/topology/test/regress/topogeo_addlinestring.sql
index 623436147afa49f453fcf8507c08831f2e09f721..564e9b4e9274bed200dd6b277614ed0766627a67 100644 (file)
--- a/topology/test/regress/topogeo_addlinestring.sql
+++ b/topology/test/regress/topogeo_addlinestring.sql
@@ -119,6 +119,13 @@ SELECT check_changes();
 SELECT 'iso_ex_2segs', TopoGeo_addLineString('city_data', 'LINESTRING(37 20, 43 19, 41 16)');
 SELECT check_changes();
 
+-- See http://trac.osgeo.org/postgis/attachment/ticket/1613
 
+SELECT '#1613.1', TopoGeo_addLineString('city_data', 'LINESTRING(556267.562954 144887.066638, 556267 144887.4)') ORDER BY 2;
+SELECT check_changes();
+SELECT '#1613.2', TopoGeo_addLineString('city_data', 'LINESTRING(556250 144887, 556267 144887.07, 556310.04 144887)') ORDER BY 2;
+SELECT check_changes();
+
+-- Cleanups
 DROP FUNCTION check_changes();
 SELECT DropTopology('city_data');

diff --git a/topology/test/regress/topogeo_addlinestring_expected b/topology/test/regress/topogeo_addlinestring_expected
index 61fd16e3989a4d71cf0cce9799ffb19d806f0a19..0ed92a18b38309013f94b88390f579a926dc517f 100644 (file)
--- a/topology/test/regress/topogeo_addlinestring_expected
+++ b/topology/test/regress/topogeo_addlinestring_expected
@@ -95,4 +95,17 @@ N|45||POINT(22 37)
 E|49|sn44|en4
 E|50|sn4|en45
 iso_ex_2segs|28
+#1613.1|51
+N|46||POINT(556267.6 144887)
+N|47||POINT(556267 144887.4)
+E|51|sn46|en47
+#1613.2|53
+#1613.2|54
+N|48||POINT(556250 144887)
+N|49||POINT(556267.6 144887)
+N|50||POINT(556310 144887)
+E|51|sn46|en49
+E|52|sn49|en47
+E|53|sn48|en49
+E|54|sn49|en50
 Topology 'city_data' dropped