From: Sandro Santilli Date: Tue, 28 Feb 2012 17:17:39 +0000 (+0000) Subject: TopoGeo_addPoint: use a more functional tolerance when snapping (#1613) X-Git-Tag: 2.0.0beta1~26 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=9f7ad76ab4b78c0e659688c9554a3ab459b469b6;p=postgis 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 6146c9d3b..3d78ec3fa 100644 --- 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 623436147..564e9b4e9 100644 --- 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 61fd16e39..0ed92a18b 100644 --- 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