From 55a1ecae41063e73e2d27ebe5afbb0db0bf54925 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Wed, 28 May 2008 23:03:11 +0000 Subject: [PATCH] Enhanced speed _ST_DWithin(g,g,d) that returns as soon as g and g are within d of each other, rather than using distance naively. Change ST_DWithin to use enhanced op. (Issue 20) git-svn-id: http://svn.osgeo.org/postgis/trunk@2796 b70326c6-7e19-0410-871a-916f4a2858ee --- lwgeom/box2d.c | 4 +-- lwgeom/liblwgeom.h | 1 + lwgeom/lwgeom_functions_analytic.c | 26 ++++++++--------- lwgeom/lwgeom_functions_basic.c | 45 ++++++++++++++++++++++++++++++ lwgeom/lwgeom_svg.c | 4 +-- lwgeom/lwpostgis.sql.in.c | 8 +++++- lwgeom/measures.c | 30 ++++++++++++-------- 7 files changed, 89 insertions(+), 29 deletions(-) diff --git a/lwgeom/box2d.c b/lwgeom/box2d.c index 1da3bb091..9e0fe8b4b 100644 --- a/lwgeom/box2d.c +++ b/lwgeom/box2d.c @@ -37,12 +37,12 @@ box2d_same(BOX2DFLOAT4 *box1, BOX2DFLOAT4 *box2) (box1->xmin==box2->xmin) && (box1->ymax==box2->ymax) && (box1->ymin==box2->ymin)); -#if 0 // changed to exact equality +#if 0 return(FPeq(box1->xmax, box2->xmax) && FPeq(box1->xmin, box2->xmin) && FPeq(box1->ymax, box2->ymax) && FPeq(box1->ymin, box2->ymin)); -#endif // 0 +#endif } BOX2DFLOAT4 * diff --git a/lwgeom/liblwgeom.h b/lwgeom/liblwgeom.h index 0ae0c0406..c5c0f79aa 100644 --- a/lwgeom/liblwgeom.h +++ b/lwgeom/liblwgeom.h @@ -946,6 +946,7 @@ extern double distance2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2); extern double distance2d_line_poly(LWLINE *line, LWPOLY *poly); extern int azimuth_pt_pt(POINT2D *p1, POINT2D *p2, double *ret); extern double lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2); +extern double lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance); extern int lwgeom_pt_inside_circle(POINT2D *p, double cx, double cy, double rad); extern int32 lwgeom_npoints(uchar *serialized); extern char ptarray_isccw(const POINTARRAY *pa); diff --git a/lwgeom/lwgeom_functions_analytic.c b/lwgeom/lwgeom_functions_analytic.c index 94bda2dbd..21b78a19a 100644 --- a/lwgeom/lwgeom_functions_analytic.c +++ b/lwgeom/lwgeom_functions_analytic.c @@ -794,7 +794,7 @@ Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS) PG_LWGEOM *out_geom = NULL; LWGEOM *out_lwgeom; gridspec grid; - //BOX3D box3d; + /* BOX3D box3d; */ if ( PG_ARGISNULL(0) ) PG_RETURN_NULL(); datum = PG_GETARG_DATUM(0); @@ -859,7 +859,7 @@ Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS) out_lwgeom->bbox = box3d_to_box2df(&box3d); } -#endif // 0 +#endif /* 0 */ #if VERBOSE elog(NOTICE, "SnapToGrid made a %s", lwgeom_typename(TYPE_GETTYPE(out_lwgeom->type))); @@ -880,7 +880,7 @@ Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS) PG_LWGEOM *out_geom = NULL; LWGEOM *out_lwgeom; gridspec grid; - //BOX3D box3d; + /* BOX3D box3d; */ POINT4D offsetpoint; if ( PG_ARGISNULL(0) ) PG_RETURN_NULL(); @@ -965,7 +965,7 @@ Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS) out_lwgeom->bbox = box3d_to_box2df(&box3d); } -#endif // 0 +#endif /* 0 */ #if VERBOSE elog(NOTICE, "SnapToGrid made a %s", lwgeom_typename(TYPE_GETTYPE(out_lwgeom->type))); @@ -1018,7 +1018,7 @@ Datum LWGEOM_line_substring(PG_FUNCTION_ARGS) opa = ptarray_substring(ipa, from, to); - if ( opa->npoints == 1 ) // Point returned + if ( opa->npoints == 1 ) /* Point returned */ olwgeom = (LWGEOM *)lwpoint_construct(iline->SRID, NULL, opa); else olwgeom = (LWGEOM *)lwline_construct(iline->SRID, NULL, opa); @@ -1366,8 +1366,8 @@ int point_in_polygon_deprecated(LWPOLY *polygon, LWPOINT *point) /* assume bbox short-circuit has already been attempted */ ring = polygon->rings[0]; - //root = createTree(ring); - //if(point_in_ring(root, &pt) != 1) + /* root = createTree(ring); */ + /* if(point_in_ring(root, &pt) != 1) */ if(point_in_ring_deprecated(polygon->rings[0], &pt) != 1) { #ifdef PGIS_DEBUG @@ -1379,8 +1379,8 @@ int point_in_polygon_deprecated(LWPOLY *polygon, LWPOINT *point) for(i=1; inrings; i++) { ring = polygon->rings[i]; - //root = createTree(ring); - //if(point_in_ring(root, &pt) != -1) + /* root = createTree(ring); */ + /* if(point_in_ring(root, &pt) != -1) */ if(point_in_ring_deprecated(polygon->rings[i], &pt) != -1) { #ifdef PGIS_DEBUG @@ -1447,8 +1447,8 @@ int point_outside_polygon_deprecated(LWPOLY *polygon, LWPOINT *point) /* assume bbox short-circuit has already been attempted */ ring = polygon->rings[0]; - //root = createTree(ring); - //if(point_in_ring(root, &pt) == -1) + /* root = createTree(ring); */ + /* if(point_in_ring(root, &pt) == -1) */ if(point_in_ring_deprecated(ring, &pt) == -1) { #ifdef PGIS_DEBUG @@ -1460,8 +1460,8 @@ int point_outside_polygon_deprecated(LWPOLY *polygon, LWPOINT *point) for(i=1; inrings; i++) { ring = polygon->rings[i]; - //root = createTree(ring); - //if(point_in_ring(root, &pt) == 1) + /* root = createTree(ring); */ + /* if(point_in_ring(root, &pt) == 1) */ if(point_in_ring_deprecated(ring, &pt) == 1) { #ifdef PGIS_DEBUG diff --git a/lwgeom/lwgeom_functions_basic.c b/lwgeom/lwgeom_functions_basic.c index 228e4e105..00fb95b4c 100644 --- a/lwgeom/lwgeom_functions_basic.c +++ b/lwgeom/lwgeom_functions_basic.c @@ -34,6 +34,7 @@ Datum LWGEOM_summary(PG_FUNCTION_ARGS); Datum LWGEOM_npoints(PG_FUNCTION_ARGS); Datum LWGEOM_nrings(PG_FUNCTION_ARGS); Datum LWGEOM_area_polygon(PG_FUNCTION_ARGS); +Datum LWGEOM_dwithin(PG_FUNCTION_ARGS); Datum postgis_uses_stats(PG_FUNCTION_ARGS); Datum postgis_autocache_bbox(PG_FUNCTION_ARGS); Datum postgis_scripts_released(PG_FUNCTION_ARGS); @@ -1629,6 +1630,50 @@ Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS) PG_RETURN_FLOAT8(mindist); } +/* Minimum 2d distance between objects in geom1 and geom2. */ +PG_FUNCTION_INFO_V1(LWGEOM_dwithin); +Datum LWGEOM_dwithin(PG_FUNCTION_ARGS) +{ + PG_LWGEOM *geom1; + PG_LWGEOM *geom2; + double mindist, tolerance; + +#ifdef PROFILE + profstart(PROF_QRUN); +#endif + + geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + tolerance = PG_GETARG_FLOAT8(2); + + if( tolerance < 0 ) { + elog(ERROR,"Tolerance cannot be less than zero\n"); + PG_RETURN_NULL(); + } + + if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) + { + elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); + PG_RETURN_NULL(); + } + + mindist = lwgeom_mindistance2d_recursive_tolerance( + SERIALIZED_FORM(geom1), + SERIALIZED_FORM(geom2), + tolerance + ); + +#ifdef PROFILE + profstop(PROF_QRUN); + profreport("dist",geom1, geom2, NULL); +#endif + + PG_FREE_IF_COPY(geom1, 0); + PG_FREE_IF_COPY(geom2, 1); + + PG_RETURN_BOOL(tolerance >= mindist); +} + /* * Maximum 2d distance between linestrings. * Returns NULL if given geoms are not linestrings. diff --git a/lwgeom/lwgeom_svg.c b/lwgeom/lwgeom_svg.c index ad99987e2..b67d46797 100644 --- a/lwgeom/lwgeom_svg.c +++ b/lwgeom/lwgeom_svg.c @@ -280,7 +280,7 @@ print_svg_path_abs(char *result, POINTARRAY *pa, int precision, int polygonRing) { getPoint2d_p(pa, u, &pt); - //close PATH with 'Z' for polygon rings if last point equals first point + /* close PATH with 'Z' for polygon rings if last point equals first point */ if(u > 0 && u == (pa->npoints - 1) && polygonRing) { POINT2D firstPoint; @@ -333,7 +333,7 @@ print_svg_path_rel(char *result, POINTARRAY *pa, int precision, int polygonRing) if(u == (pa->npoints - 1) && polygonRing) { - //close PATH with 'z' if last point equals first point + /* close PATH with 'z' if last point equals first point */ POINT2D firstPoint; getPoint2d_p(pa, 0, &firstPoint); if(pt.x == firstPoint.x && pt.y == firstPoint.y) diff --git a/lwgeom/lwpostgis.sql.in.c b/lwgeom/lwpostgis.sql.in.c index 6ecaea6b7..ddf69d0bd 100644 --- a/lwgeom/lwpostgis.sql.in.c +++ b/lwgeom/lwpostgis.sql.in.c @@ -4430,10 +4430,16 @@ CREATEFUNCTION ST_Touches(geometry,geometry) AS 'SELECT $1 && $2 AND _ST_Touches($1,$2)' LANGUAGE 'SQL' _IMMUTABLE; -- WITH (iscachable); +-- Availability: 1.3.4 +CREATEFUNCTION _ST_DWithin(geometry,geometry,float8) + RETURNS boolean + AS '@MODULE_FILENAME@', 'LWGEOM_dwithin' + LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable); + -- Availability: 1.2.2 CREATEFUNCTION ST_DWithin(geometry, geometry, float8) RETURNS boolean - AS 'SELECT $1 && ST_Expand($2,$3) AND $2 && ST_Expand($1,$3) AND ST_Distance($1, $2) < $3' + AS 'SELECT $1 && ST_Expand($2,$3) AND $2 && ST_Expand($1,$3) AND _ST_DWithin($1, $2, $3)' LANGUAGE 'SQL' _IMMUTABLE; -- WITH (iscachable); -- Deprecation in 1.2.3 diff --git a/lwgeom/measures.c b/lwgeom/measures.c index 546aad210..a6d493e87 100644 --- a/lwgeom/measures.c +++ b/lwgeom/measures.c @@ -47,7 +47,7 @@ int pt_in_ring_2d(POINT2D *p, POINTARRAY *ring) #ifdef PGIS_DEBUG lwnotice("pt_in_ring_2d called with point: %g %g", p->x, p->y); - //printPA(ring); + /* printPA(ring); */ #endif /* loop through all edges of the polygon */ @@ -91,10 +91,10 @@ double distance2d_pt_pt(POINT2D *p1, POINT2D *p2) return sqrt ( hside*hside + vside*vside ); - // the above is more readable - //return sqrt( - // (p2->x-p1->x) * (p2->x-p1->x) + (p2->y-p1->y) * (p2->y-p1->y) - // ); + /* the above is more readable + return sqrt( + (p2->x-p1->x) * (p2->x-p1->x) + (p2->y-p1->y) * (p2->y-p1->y) + ); */ } /*distance2d from p to line A->B */ @@ -510,7 +510,7 @@ double distance2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2) poly2->rings[j]); if ( d <= 0 ) return 0.0; - // mindist is -1 when not yet set + /* mindist is -1 when not yet set */ if (mindist > -1) mindist = LW_MIN(mindist, d); else mindist = d; #ifdef PGIS_DEBUG @@ -676,6 +676,12 @@ lwnotice("in lwgeom_polygon_perimeter (%d rings)", poly->nrings); double lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2) +{ + return lwgeom_mindistance2d_recursive_tolerance( lw1, lw2, 0.0 ); +} + +double +lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance) { LWGEOM_INSPECTED *in1, *in2; int i, j; @@ -688,13 +694,13 @@ lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2) { uchar *g1 = lwgeom_getsubgeometry_inspected(in1, i); int t1 = lwgeom_getType(g1[0]); - double dist=0; + double dist=tolerance; /* it's a multitype... recurse */ if ( t1 >= 4 ) { - dist = lwgeom_mindistance2d_recursive(g1, lw2); - if ( dist == 0 ) return 0.0; /* can't be closer */ + dist = lwgeom_mindistance2d_recursive_tolerance(g1, lw2, tolerance); + if ( dist <= tolerance ) return tolerance; /* can't be closer */ if ( mindist == -1 ) mindist = dist; else mindist = LW_MIN(dist, mindist); continue; @@ -779,7 +785,7 @@ lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2) } else /* it's a multitype... recurse */ { - dist = lwgeom_mindistance2d_recursive(g1, g2); + dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance); } if (mindist == -1 ) mindist = dist; @@ -791,7 +797,7 @@ lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2) #endif - if (mindist <= 0.0) return 0.0; /* can't be closer */ + if (mindist <= tolerance) return tolerance; /* can't be closer */ } @@ -802,6 +808,8 @@ lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2) return mindist; } + + int lwgeom_pt_inside_circle(POINT2D *p, double cx, double cy, double rad) { -- 2.50.1