Enhanced speed _ST_DWithin(g,g,d) that returns as soon as g and g are within d of...

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

diff --git a/lwgeom/box2d.c b/lwgeom/box2d.c
index 1da3bb091381d415d2ccf590c21e736eae94e23b..9e0fe8b4b64e35c06b9d71b57ce27a394ed362c9 100644 (file)
--- 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 0ae0c0406b7a24ff87a0db0547429c10c938f726..c5c0f79aa73e3913e6676763de56428f5b2a6788 100644 (file)
--- 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 94bda2dbdec11b076fb584a2e0df99e39ab65e5b..21b78a19a49ba88eab86bee988c5b2b6171ca915 100644 (file)
--- 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; i<polygon->nrings; 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; i<polygon->nrings; 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 228e4e10510d61943a5847626773bcba1c71cbcf..00fb95b4c726388c4c7f31591e56b6c8e25fc8b0 100644 (file)
--- 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 ad99987e25c3340fe79c5eca92a672540dac3a01..b67d46797a7fc7ee5879b0c2bffd2ab2edb20265 100644 (file)
--- 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 6ecaea6b7abaeca66a6e02dc1ffea1400dfb6589..ddf69d0bd465f5a29db89bc5d5c5a5d26a6411a8 100644 (file)
--- 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 546aad2106792d1220bd05c96a590a052c9a2e60..a6d493e871c9b3621deb0bbf5b5135ee9b0c47e5 100644 (file)
--- 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)
 {