Add box-wise ORDER BY for KNN (#701)

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

diff --git a/postgis/gserialized_gist_2d.c b/postgis/gserialized_gist_2d.c
index f53e47b4b97faeff0cd61c96b34f97616427d08a..436708566025dd592baa9715b717bad63ae2eeda 100644 (file)
--- a/postgis/gserialized_gist_2d.c
+++ b/postgis/gserialized_gist_2d.c
@@ -80,12 +80,13 @@ Datum gserialized_overleft_2d(PG_FUNCTION_ARGS);
 Datum gserialized_overright_2d(PG_FUNCTION_ARGS);
 Datum gserialized_overabove_2d(PG_FUNCTION_ARGS);
 Datum gserialized_overbelow_2d(PG_FUNCTION_ARGS);
-Datum gserialized_boxdistance_2d(PG_FUNCTION_ARGS);
+Datum gserialized_distance_box_2d(PG_FUNCTION_ARGS);
+Datum gserialized_distance_centroid_2d(PG_FUNCTION_ARGS);
 
 /*
 ** true/false test function type
 */
-typedef bool (*box2df_predicate)(BOX2DF *a, BOX2DF *b);
+typedef bool (*box2df_predicate)(const BOX2DF *a, const BOX2DF *b);
 
 
 #if POSTGIS_DEBUG_LEVEL > 0
@@ -234,7 +235,7 @@ static inline void box2df_validate(BOX2DF *b)
        return;
 }
 
-static bool box2df_overlaps(BOX2DF *a, BOX2DF *b)
+static bool box2df_overlaps(const BOX2DF *a, const BOX2DF *b)
 {
        if ( (a->xmin > b->xmax) || (b->xmin > a->xmax) ||
             (a->ymin > b->ymax) || (b->ymin > a->ymax) )
@@ -245,7 +246,7 @@ static bool box2df_overlaps(BOX2DF *a, BOX2DF *b)
        return TRUE;
 }
 
-static bool box2df_contains(BOX2DF *a, BOX2DF *b)
+static bool box2df_contains(const BOX2DF *a, const BOX2DF *b)
 {
        if ( (a->xmin > b->xmin) || (a->xmax < b->xmax) ||
             (a->ymin > b->ymin) || (a->ymax < b->ymax) )
@@ -256,13 +257,13 @@ static bool box2df_contains(BOX2DF *a, BOX2DF *b)
        return TRUE;
 }
 
-static bool box2df_within(BOX2DF *a, BOX2DF *b)
+static bool box2df_within(const BOX2DF *a, const BOX2DF *b)
 {
        POSTGIS_DEBUG(5, "entered function");
        return box2df_contains(b,a);
 }
 
-static bool box2df_equals(BOX2DF *a, BOX2DF *b)
+static bool box2df_equals(const BOX2DF *a, const BOX2DF *b)
 {
        if ( (a->xmin != b->xmin) || (a->xmax != b->xmax) ||
             (a->ymin != b->ymin) || (a->ymax != b->ymax) )
@@ -273,49 +274,49 @@ static bool box2df_equals(BOX2DF *a, BOX2DF *b)
        return TRUE;
 }
 
-static bool box2df_overleft(BOX2DF *a, BOX2DF *b)
+static bool box2df_overleft(const BOX2DF *a, const BOX2DF *b)
 {
        /* a.xmax <= b.xmax */
        return a->xmax <= b->xmax;
 }
 
-static bool box2df_left(BOX2DF *a, BOX2DF *b)
+static bool box2df_left(const BOX2DF *a, const BOX2DF *b)
 {
        /* a.xmax < b.xmin */
        return a->xmax < b->xmin;
 }
 
-static bool box2df_right(BOX2DF *a, BOX2DF *b)
+static bool box2df_right(const BOX2DF *a, const BOX2DF *b)
 {
        /* a.xmin > b.xmax */
        return a->xmin > b->xmax;
 }
 
-static bool box2df_overright(BOX2DF *a, BOX2DF *b)
+static bool box2df_overright(const BOX2DF *a, const BOX2DF *b)
 {
        /* a.xmin >= b.xmin */
        return a->xmin >= b->xmin;
 }
 
-static bool box2df_overbelow(BOX2DF *a, BOX2DF *b)
+static bool box2df_overbelow(const BOX2DF *a, const BOX2DF *b)
 {
        /* a.ymax <= b.ymax */
        return a->ymax <= b->ymax;
 }
 
-static bool box2df_below(BOX2DF *a, BOX2DF *b)
+static bool box2df_below(const BOX2DF *a, const BOX2DF *b)
 {
        /* a.ymax < b.ymin */
        return a->ymax < b->ymin;
 }
 
-static bool box2df_above(BOX2DF *a, BOX2DF *b)
+static bool box2df_above(const BOX2DF *a, const BOX2DF *b)
 {
        /* a.ymin > b.ymax */
        return a->ymin > b->ymax;
 }
 
-static bool box2df_overabove(BOX2DF *a, BOX2DF *b)
+static bool box2df_overabove(const BOX2DF *a, const BOX2DF *b)
 {
        /* a.ymin >= b.ymin */
        return a->ymin >= b->ymin;
@@ -325,7 +326,7 @@ static bool box2df_overabove(BOX2DF *a, BOX2DF *b)
 * Calculate the centroid->centroid distance between the boxes.
 * We return the square distance to avoid a call to sqrt.
 */
-static double box2df_distance_leaf(BOX2DF *a, BOX2DF *b)
+static double box2df_distance_leaf_centroid(const BOX2DF *a, const BOX2DF *b)
 {
     /* The centroid->centroid distance between the boxes */
     double a_x = (a->xmax + a->xmin) / 2.0;
@@ -343,7 +344,7 @@ static double box2df_distance_leaf(BOX2DF *a, BOX2DF *b)
 * between the boxes.
 * We return the square distance to avoid a call to sqrt.
 */
-static double box2df_distance_node(BOX2DF *node, BOX2DF *query)
+static double box2df_distance_node_centroid(const BOX2DF *node, const BOX2DF *query)
 {
     BOX2DF q;
     double qx, qy;
@@ -413,6 +414,61 @@ static double box2df_distance_node(BOX2DF *node, BOX2DF *query)
     return d;
 }
 
+/* Quick distance function */
+static inline double pt_distance(double ax, double ay, double bx, double by)
+{
+       return sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
+}
+
+/**
+* Calculate the box->box distance.
+*/
+static double box2df_distance(const BOX2DF *a, const BOX2DF *b)
+{
+    /* Check for overlap */
+    if ( box2df_overlaps(a, b) )
+        return 0.0;
+
+    if ( box2df_left(a, b) )
+    {
+        if ( box2df_above(a, b) )
+                       return pt_distance(a->xmax, a->ymin, b->xmin, b->ymax);
+               if ( box2df_below(a, b) )
+                       return pt_distance(a->xmax, a->ymax, b->xmin, b->ymin);
+               else
+                       return b->xmin - a->xmax;
+       }
+       if ( box2df_right(a, b) )
+       {
+        if ( box2df_above(a, b) )
+                       return pt_distance(a->xmin, a->ymin, b->xmax, b->ymax);
+               if ( box2df_below(a, b) )
+                       return pt_distance(a->xmin, a->ymax, b->xmax, b->ymin);
+               else
+                       return a->xmin - b->xmax;
+       }
+       if ( box2df_above(a, b) )
+       {
+               if ( box2df_left(a, b) )
+                       return pt_distance(a->xmax, a->ymin, b->xmin, b->ymax);
+               if ( box2df_right(a, b) )
+                       return pt_distance(a->xmin, a->ymin, b->xmax, b->ymax);
+               else
+                       return a->ymin - b->ymax;
+       }
+       if ( box2df_below(a, b) )
+       {
+               if ( box2df_left(a, b) )
+                       return pt_distance(a->xmax, a->ymax, b->xmin, b->ymin);
+               if ( box2df_right(a, b) )
+                       return pt_distance(a->xmin, a->ymax, b->xmax, b->ymin);
+               else
+                       return b->ymin - a->ymax;
+       }
+       
+       return INFINITY;
+}
+
 
 /**
 * Peak into a #GSERIALIZED datum to find the bounding box. If the
@@ -500,12 +556,32 @@ gserialized_datum_predicate_2d(Datum gs1, Datum gs2, box2df_predicate predicate)
 * GiST 2-D Index Operator Functions
 */
 
-PG_FUNCTION_INFO_V1(gserialized_boxdistance_2d);
-Datum gserialized_boxdistance_2d(PG_FUNCTION_ARGS)
+PG_FUNCTION_INFO_V1(gserialized_distance_centroid_2d);
+Datum gserialized_distance_centroid_2d(PG_FUNCTION_ARGS)
+{
+       BOX2DF b1, b2;
+       Datum gs1 = PG_GETARG_DATUM(0);
+       Datum gs2 = PG_GETARG_DATUM(1);    
+       
+       POSTGIS_DEBUG(3, "entered function");
+
+       /* Must be able to build box for each argument (ie, not empty geometry). */
+       if ( (gserialized_datum_get_box2df_p(gs1, &b1) == LW_SUCCESS) &&
+            (gserialized_datum_get_box2df_p(gs2, &b2) == LW_SUCCESS) )
+       {           
+               double distance = box2df_distance_leaf_centroid(&b1, &b2);
+               POSTGIS_DEBUGF(3, "got boxes %s and %s", box2df_to_string(&b1), box2df_to_string(&b2));
+               PG_RETURN_FLOAT8(distance);
+       }
+       PG_RETURN_FLOAT8(INFINITY);
+}
+
+PG_FUNCTION_INFO_V1(gserialized_distance_box_2d);
+Datum gserialized_distance_box_2d(PG_FUNCTION_ARGS)
 {
        BOX2DF b1, b2;
-    Datum gs1 = PG_GETARG_DATUM(0);
-    Datum gs2 = PG_GETARG_DATUM(1);    
+       Datum gs1 = PG_GETARG_DATUM(0);
+       Datum gs2 = PG_GETARG_DATUM(1);    
        
        POSTGIS_DEBUG(3, "entered function");
 
@@ -513,11 +589,11 @@ Datum gserialized_boxdistance_2d(PG_FUNCTION_ARGS)
        if ( (gserialized_datum_get_box2df_p(gs1, &b1) == LW_SUCCESS) &&
             (gserialized_datum_get_box2df_p(gs2, &b2) == LW_SUCCESS) )
        {           
-        double distance = box2df_distance_leaf(&b1, &b2);
+               double distance = box2df_distance(&b1, &b2);
                POSTGIS_DEBUGF(3, "got boxes %s and %s", box2df_to_string(&b1), box2df_to_string(&b2));
-        PG_RETURN_FLOAT8(distance);
+               PG_RETURN_FLOAT8(distance);
        }
-    PG_RETURN_FLOAT8(INFINITY);
+       PG_RETURN_FLOAT8(INFINITY);
 }
 
 PG_FUNCTION_INFO_V1(gserialized_same_2d);
@@ -865,19 +941,23 @@ Datum gserialized_gist_consistent_2d(PG_FUNCTION_ARGS)
 ** represents the distance to the index entry; for an internal tree node, the
 ** result must be the smallest distance that any child entry could have.
 ** 
+** Strategy 13 = centroid-based distance tests
+** Strategy 14 = box-based distance tests
 */
 PG_FUNCTION_INFO_V1(gserialized_gist_distance_2d);
 Datum gserialized_gist_distance_2d(PG_FUNCTION_ARGS)
 {
        GISTENTRY *entry = (GISTENTRY*) PG_GETARG_POINTER(0);
        BOX2DF query_box;
+       BOX2DF *entry_box;
        StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
        double distance;
 
        POSTGIS_DEBUG(4, "[GIST] 'distance' function called");
 
-    /* We are using '13' as the gist distance strategy number */
-    if ( strategy != 13 ) {
+    /* We are using '13' as the gist distance-betweeen-centroids strategy number 
+    *  and '14' as the gist distance-between-boxes strategy number */
+    if ( strategy != 13 && strategy != 14 ) {
         elog(ERROR, "unrecognized strategy number: %d", strategy);
                PG_RETURN_FLOAT8(INFINITY);
        }
@@ -889,18 +969,26 @@ Datum gserialized_gist_distance_2d(PG_FUNCTION_ARGS)
                PG_RETURN_FLOAT8(INFINITY);
        }
 
+       /* Get the entry box */
+    entry_box = (BOX2DF*)DatumGetPointer(entry->key);
+       
+       /* Box-style distance test */
+       if ( strategy == 14 )
+       {
+               distance = (double)box2df_distance(entry_box, &query_box);
+               PG_RETURN_FLOAT8(distance);
+       }
+
        /* Treat leaf node tests different from internal nodes */
        if (GIST_LEAF(entry))
        {
-           /* Calculate distance to centroid for leaves */
-        BOX2DF *leaf_box = (BOX2DF*)DatumGetPointer(entry->key);
-        distance = (double)box2df_distance_leaf(leaf_box, &query_box);
+           /* Calculate distance to leaves */
+               distance = (double)box2df_distance_leaf_centroid(entry_box, &query_box);
        }
        else
        {
-           /* Calculate distance to nearest corner for internal nodes */
-        BOX2DF *node_box = (BOX2DF*)DatumGetPointer(entry->key);
-        distance = (double)box2df_distance_node(node_box, &query_box);
+           /* Calculate distance for internal nodes */
+               distance = (double)box2df_distance_node_centroid(entry_box, &query_box);
        }
 
        PG_RETURN_FLOAT8(distance);

diff --git a/postgis/postgis.sql.in.c b/postgis/postgis.sql.in.c
index c385a44d51a3e1377771c2421db252684475b794..7b1e22f31bf51498dbe2852d3c488c4846e1676f 100644 (file)
--- a/postgis/postgis.sql.in.c
+++ b/postgis/postgis.sql.in.c
@@ -515,7 +515,7 @@ CREATE TYPE box2df (
 );\r
 \r
 -- Availability: 2.0.0\r
-CREATE OR REPLACE FUNCTION geometry_gist_boxdistance_2d(internal,geometry,int4) \r
+CREATE OR REPLACE FUNCTION geometry_gist_distance_2d(internal,geometry,int4) \r
        RETURNS float8 \r
        AS 'MODULE_PATHNAME' ,'gserialized_gist_distance_2d'\r
        LANGUAGE 'C';\r
@@ -603,16 +603,27 @@ CREATE OPERATOR ~= (
 );\r
 \r
 -- Availability: 2.0.0\r
-CREATE OR REPLACE FUNCTION geometry_boxdistance(geometry, geometry) \r
+CREATE OR REPLACE FUNCTION geometry_distance_centroid(geometry, geometry) \r
        RETURNS float8 \r
-       AS 'MODULE_PATHNAME' ,'gserialized_boxdistance_2d'\r
+       AS 'MODULE_PATHNAME' ,'gserialized_distance_centroid_2d'\r
+       LANGUAGE 'C' IMMUTABLE STRICT;\r
+\r
+-- Availability: 2.0.0\r
+CREATE OR REPLACE FUNCTION geometry_distance_box(geometry, geometry) \r
+       RETURNS float8 \r
+       AS 'MODULE_PATHNAME' ,'gserialized_distance_box_2d'\r
        LANGUAGE 'C' IMMUTABLE STRICT;\r
 \r
 #if POSTGIS_PGSQL_VERSION >= 91\r
 CREATE OPERATOR <-> (\r
-    LEFTARG = geometry, RIGHTARG = geometry, PROCEDURE = geometry_boxdistance,\r
+    LEFTARG = geometry, RIGHTARG = geometry, PROCEDURE = geometry_distance_centroid,\r
     COMMUTATOR = '<->'\r
 );\r
+\r
+CREATE OPERATOR <#> (\r
+    LEFTARG = geometry, RIGHTARG = geometry, PROCEDURE = geometry_distance_box,\r
+    COMMUTATOR = '<#>'\r
+);\r
 #endif\r
 \r
 -- Availability: 2.0.0\r
@@ -753,7 +764,8 @@ CREATE OPERATOR CLASS gist_geometry_ops_2d
        OPERATOR        12       |&> ,\r
 #if POSTGIS_PGSQL_VERSION >= 91\r
        OPERATOR        13       <-> FOR ORDER BY pg_catalog.float_ops,\r
-       FUNCTION        8        geometry_gist_boxdistance_2d (internal, geometry, int4),\r
+       OPERATOR        14       <#> FOR ORDER BY pg_catalog.float_ops,\r
+       FUNCTION        8        geometry_gist_distance_2d (internal, geometry, int4),\r
 #endif\r
        FUNCTION        1        geometry_gist_consistent_2d (internal, geometry, int4),\r
        FUNCTION        2        geometry_gist_union_2d (bytea, internal),\r