Added ST_DWithin(raster, raster) and regression tests. Ticket is #1922

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

diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c
index c0d3def85134637d2746bba942e98628f2f11c07..e418f3950c58202e1415d15440c136da13733cb8 100644 (file)
--- a/raster/rt_core/rt_api.c
+++ b/raster/rt_core/rt_api.c
@@ -11862,6 +11862,8 @@ int rt_raster_within_distance(
        lwgeom_free(surface1);
        lwgeom_free(surface2);
 
+       RASTER_DEBUGF(3, "(mindist, distance) = (%f, %f)", mindist, distance);
+
        /* if distance >= mindist, true */
        if (FLT_EQ(mindist, distance) || distance > mindist)
                *dwithin = 1;

diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c
index eaa40f436181aab4f3dcdc197dff1f980e92b969..87933b528e43f361effe68c330afdfdd4a06949d 100644 (file)
--- a/raster/rt_pg/rt_pg.c
+++ b/raster/rt_pg/rt_pg.c
@@ -288,6 +288,9 @@ Datum RASTER_covers(PG_FUNCTION_ARGS);
 /* determine if the first raster is covered by the second raster */
 Datum RASTER_coveredby(PG_FUNCTION_ARGS);
 
+/* determine if the two rasters are within the specified distance of each other */
+Datum RASTER_dwithin(PG_FUNCTION_ARGS);
+
 /* determine if two rasters are aligned */
 Datum RASTER_sameAlignment(PG_FUNCTION_ARGS);
 
@@ -10707,6 +10710,147 @@ Datum RASTER_coveredby(PG_FUNCTION_ARGS)
        PG_RETURN_BOOL(result);
 }
 
+/**
+ * See if the two rasters are within the specified distance of each other
+ */
+PG_FUNCTION_INFO_V1(RASTER_dwithin);
+Datum RASTER_dwithin(PG_FUNCTION_ARGS)
+{
+       const int set_count = 2;
+       rt_pgraster *pgrast[2];
+       int pgrastpos[2] = {-1, -1};
+       rt_raster rast[2] = {NULL};
+       uint32_t bandindex[2] = {0};
+       uint32_t hasbandindex[2] = {0};
+       double distance = 0;
+
+       uint32_t i;
+       uint32_t j;
+       uint32_t k;
+       uint32_t numBands;
+       int rtn;
+       int result;
+
+       for (i = 0, j = 0; i < set_count; i++) {
+               /* pgrast is null, return null */
+               if (PG_ARGISNULL(j)) {
+                       for (k = 0; k < i; k++) {
+                               rt_raster_destroy(rast[k]);
+                               PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
+                       }
+                       PG_RETURN_NULL();
+               }
+               pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
+               pgrastpos[i] = j;
+               j++;
+
+               /* raster */
+               rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
+               if (!rast[i]) {
+                       elog(ERROR, "RASTER_dwithin: Could not deserialize the %s raster", i < 1 ? "first" : "second");
+                       for (k = 0; k <= i; k++) {
+                               if (k < i)
+                                       rt_raster_destroy(rast[k]);
+                               PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
+                       }
+                       PG_RETURN_NULL();
+               }
+
+               /* numbands */
+               numBands = rt_raster_get_num_bands(rast[i]);
+               if (numBands < 1) {
+                       elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
+                       if (i > 0) i++;
+                       for (k = 0; k < i; k++) {
+                               rt_raster_destroy(rast[k]);
+                               PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
+                       }
+                       PG_RETURN_NULL();
+               }
+
+               /* band index */
+               if (!PG_ARGISNULL(j)) {
+                       bandindex[i] = PG_GETARG_INT32(j);
+                       if (bandindex[i] < 1 || bandindex[i] > numBands) {
+                               elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
+                               if (i > 0) i++;
+                               for (k = 0; k < i; k++) {
+                                       rt_raster_destroy(rast[k]);
+                                       PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
+                               }
+                               PG_RETURN_NULL();
+                       }
+                       hasbandindex[i] = 1;
+               }
+               else
+                       hasbandindex[i] = 0;
+               POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
+               POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
+               j++;
+       }
+
+       /* distance */
+       if (PG_ARGISNULL(4)) {
+               elog(NOTICE, "Distance cannot be NULL.  Returning NULL");
+               for (k = 0; k < set_count; k++) {
+                       rt_raster_destroy(rast[k]);
+                       PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
+               }
+               PG_RETURN_NULL();
+       }
+
+       distance = PG_GETARG_FLOAT8(4);
+       if (distance < 0) {
+               elog(NOTICE, "Distance cannot be less than zero.  Returning NULL");
+               for (k = 0; k < set_count; k++) {
+                       rt_raster_destroy(rast[k]);
+                       PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
+               }
+               PG_RETURN_NULL();
+       }
+
+       /* hasbandindex must be balanced */
+       if (
+               (hasbandindex[0] && !hasbandindex[1]) ||
+               (!hasbandindex[0] && hasbandindex[1])
+       ) {
+               elog(NOTICE, "Missing band index.  Band indices must be provided for both rasters if any one is provided");
+               for (k = 0; k < set_count; k++) {
+                       rt_raster_destroy(rast[k]);
+                       PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
+               }
+               PG_RETURN_NULL();
+       }
+
+       /* SRID must match */
+       if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
+               elog(ERROR, "The two rasters provided have different SRIDs");
+               for (k = 0; k < set_count; k++) {
+                       rt_raster_destroy(rast[k]);
+                       PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
+               }
+               PG_RETURN_NULL();
+       }
+
+       rtn = rt_raster_within_distance(
+               rast[0], (hasbandindex[0] ? bandindex[0] - 1 : -1),
+               rast[1], (hasbandindex[1] ? bandindex[1] - 1 : -1),
+               distance,
+               &result
+       );
+       for (k = 0; k < set_count; k++) {
+               rt_raster_destroy(rast[k]);
+               PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
+       }
+
+       if (!rtn) {
+               elog(ERROR, "RASTER_dwithin: Unable to test that the two rasters are within the specified distance of each other");
+               PG_RETURN_NULL();
+       }
+
+       PG_RETURN_BOOL(result);
+}
+
 /**
  * See if two rasters are aligned
  */

diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c
index 811d59694155ec96b5d673d99e8af386b3e7cc0d..d2893555fb91060629a42d1283cf7f7e3b74386d 100644 (file)
--- a/raster/rt_pg/rtpostgis.sql.in.c
+++ b/raster/rt_pg/rtpostgis.sql.in.c
@@ -3444,6 +3444,28 @@ CREATE OR REPLACE FUNCTION st_within(rast1 raster, rast2 raster)
        LANGUAGE 'sql' IMMUTABLE
        COST 1000;
 
+-----------------------------------------------------------------------
+-- ST_DWithin(raster, raster)
+-----------------------------------------------------------------------
+
+CREATE OR REPLACE FUNCTION _st_dwithin(rast1 raster, nband1 integer, rast2 raster, nband2 integer, distance double precision)
+       RETURNS boolean
+       AS 'MODULE_PATHNAME', 'RASTER_dwithin'
+       LANGUAGE 'c' IMMUTABLE STRICT
+       COST 1000;
+
+CREATE OR REPLACE FUNCTION st_dwithin(rast1 raster, nband1 integer, rast2 raster, nband2 integer, distance double precision)
+       RETURNS boolean
+       AS $$ SELECT $1 && ST_Expand(ST_ConvexHull($3), $5) AND $3 && ST_Expand(ST_ConvexHull($1), $5) AND CASE WHEN $2 IS NULL OR $4 IS NULL THEN _st_dwithin(st_convexhull($1), st_convexhull($3), $5) ELSE _st_dwithin($1, $2, $3, $4, $5) END $$
+       LANGUAGE 'sql' IMMUTABLE
+       COST 1000;
+
+CREATE OR REPLACE FUNCTION st_dwithin(rast1 raster, rast2 raster, distance double precision)
+       RETURNS boolean
+       AS $$ SELECT st_dwithin($1, NULL::integer, $2, NULL::integer, $3) $$
+       LANGUAGE 'sql' IMMUTABLE
+       COST 1000;
+
 -----------------------------------------------------------------------
 -- ST_Disjoint(raster, raster)
 -----------------------------------------------------------------------

diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in
index 9e83d3f089e2e8940f808d837b7af4fd53ccf102..0917e0e0fc2f677b893b64bb7b3ce07be5ac07ce 100644 (file)
--- a/raster/test/regress/Makefile.in
+++ b/raster/test/regress/Makefile.in
@@ -133,6 +133,7 @@ TEST_SREL = \
        rt_coveredby \
        rt_within \
        rt_disjoint \
+       rt_dwithin \
        rt_samealignment
        
 TEST_BUGS = \

diff --git a/raster/test/regress/rt_dwithin.sql b/raster/test/regress/rt_dwithin.sql
new file mode 100644 (file)
index 0000000..2d83950
--- /dev/null
+++ b/raster/test/regress/rt_dwithin.sql
@@ -0,0 +1,218 @@
+--SET client_min_messages TO warning;
+
+DROP TABLE IF EXISTS raster_dwithin_rast;
+CREATE TABLE raster_dwithin_rast (
+       rid integer,
+       rast raster
+);
+CREATE OR REPLACE FUNCTION make_test_raster(rid integer, width integer DEFAULT 2, height integer DEFAULT 2, ul_x double precision DEFAULT 0, ul_y double precision DEFAULT 0, skew_x double precision DEFAULT 0, skew_y double precision DEFAULT 0)
+       RETURNS void
+       AS $$
+       DECLARE
+               x int;
+               y int;
+               rast raster;
+       BEGIN
+               rast := ST_MakeEmptyRaster(width, height, ul_x, ul_y, 1, 1, skew_x, skew_y, 0);
+               rast := ST_AddBand(rast, 1, '8BUI', 1, 0);
+
+
+               INSERT INTO raster_dwithin_rast VALUES (rid, rast);
+
+               RETURN;
+       END;
+       $$ LANGUAGE 'plpgsql';
+SELECT make_test_raster(0, 2, 2, -1, -1);
+SELECT make_test_raster(1, 2, 2);
+SELECT make_test_raster(2, 3, 3);
+SELECT make_test_raster(3, 2, 2, -1, -10);
+DROP FUNCTION make_test_raster(integer, integer, integer, double precision, double precision, double precision, double precision);
+
+INSERT INTO raster_dwithin_rast VALUES (10, (
+       SELECT
+               ST_SetValue(rast, 1, 1, 1, 0)
+       FROM raster_dwithin_rast
+       WHERE rid = 1
+));
+INSERT INTO raster_dwithin_rast VALUES (11, (
+       SELECT
+               ST_SetValue(rast, 1, 2, 1, 0)
+       FROM raster_dwithin_rast
+       WHERE rid = 1
+));
+INSERT INTO raster_dwithin_rast VALUES (12, (
+       SELECT
+               ST_SetValue(
+                       ST_SetValue(
+                               ST_SetValue(rast, 1, 1, 1, 0),
+                               1, 2, 1, 0
+                       ),
+                       1, 1, 2, 0
+               )
+       FROM raster_dwithin_rast
+       WHERE rid = 1
+));
+INSERT INTO raster_dwithin_rast VALUES (13, (
+       SELECT
+               ST_SetValue(
+                       ST_SetValue(
+                               ST_SetValue(
+                                       ST_SetValue(rast, 1, 1, 1, 0),
+                                       1, 2, 1, 0
+                               ),
+                               1, 1, 2, 0
+                       ),
+                       1, 2, 2, 0
+               )
+       FROM raster_dwithin_rast
+       WHERE rid = 1
+));
+INSERT INTO raster_dwithin_rast VALUES (14, (
+       SELECT
+               ST_SetUpperLeft(rast, 2, 0)
+       FROM raster_dwithin_rast
+       WHERE rid = 1
+));
+INSERT INTO raster_dwithin_rast VALUES (15, (
+       SELECT
+               ST_SetScale(
+                       ST_SetUpperLeft(rast, 0.1, 0.1),
+                       0.4, 0.4
+               )
+       FROM raster_dwithin_rast
+       WHERE rid = 1
+));
+INSERT INTO raster_dwithin_rast VALUES (16, (
+       SELECT
+               ST_SetScale(
+                       ST_SetUpperLeft(rast, -0.1, 0.1),
+                       0.4, 0.4
+               )
+       FROM raster_dwithin_rast
+       WHERE rid = 1
+));
+
+INSERT INTO raster_dwithin_rast VALUES (20, (
+       SELECT
+               ST_SetUpperLeft(rast, -2, -2)
+       FROM raster_dwithin_rast
+       WHERE rid = 2
+));
+INSERT INTO raster_dwithin_rast VALUES (21, (
+       SELECT
+               ST_SetValue(
+                       ST_SetValue(
+                               ST_SetValue(rast, 1, 1, 1, 0),
+                               1, 2, 2, 0
+                       ),
+                       1, 3, 3, 0
+               )
+       FROM raster_dwithin_rast
+       WHERE rid = 20
+));
+INSERT INTO raster_dwithin_rast VALUES (22, (
+       SELECT
+               ST_SetValue(
+                       ST_SetValue(
+                               rast, 1, 3, 2, 0
+                       ),
+                       1, 2, 3, 0
+               )
+       FROM raster_dwithin_rast
+       WHERE rid = 21
+));
+INSERT INTO raster_dwithin_rast VALUES (23, (
+       SELECT
+               ST_SetValue(
+                       ST_SetValue(
+                               rast, 1, 3, 1, 0
+                       ),
+                       1, 1, 3, 0
+               )
+       FROM raster_dwithin_rast
+       WHERE rid = 22
+));
+
+INSERT INTO raster_dwithin_rast VALUES (30, (
+       SELECT
+               ST_SetSkew(rast, -0.5, 0.5)
+       FROM raster_dwithin_rast
+       WHERE rid = 2
+));
+INSERT INTO raster_dwithin_rast VALUES (31, (
+       SELECT
+               ST_SetSkew(rast, -1, 1)
+       FROM raster_dwithin_rast
+       WHERE rid = 2
+));
+INSERT INTO raster_dwithin_rast VALUES (32, (
+       SELECT
+               ST_SetSkew(rast, 1, -1)
+       FROM raster_dwithin_rast
+       WHERE rid = 2
+));
+
+SELECT
+       '1.1',
+       r1.rid,
+       r2.rid,
+       ST_DWithin(r1.rast, NULL, r2.rast, NULL, NULL)
+FROM raster_dwithin_rast r1
+CROSS JOIN raster_dwithin_rast r2
+WHERE r1.rid = 0;
+
+SELECT
+       '1.2',
+       r1.rid,
+       r2.rid,
+       ST_DWithin(r1.rast, NULL, r2.rast, NULL, -1)
+FROM raster_dwithin_rast r1
+CROSS JOIN raster_dwithin_rast r2
+WHERE r1.rid = 0;
+
+SELECT
+       '1.3',
+       r1.rid,
+       r2.rid,
+       ST_DWithin(r1.rast, NULL, r2.rast, NULL, 0)
+FROM raster_dwithin_rast r1
+CROSS JOIN raster_dwithin_rast r2
+WHERE r1.rid = 0;
+
+SELECT
+       '2.1',
+       r1.rid,
+       r2.rid,
+       ST_DWithin(r1.rast, 1, r2.rast, 1, 0)
+FROM raster_dwithin_rast r1
+CROSS JOIN raster_dwithin_rast r2
+WHERE r1.rid = 0;
+
+SELECT
+       '2.2',
+       r1.rid,
+       r2.rid,
+       ST_DWithin(r1.rast, 1, r2.rast, 1, 1)
+FROM raster_dwithin_rast r1
+CROSS JOIN raster_dwithin_rast r2
+WHERE r1.rid = 0;
+
+SELECT
+       '2.3',
+       r1.rid,
+       r2.rid,
+       ST_DWithin(r1.rast, 1, r2.rast, 1, 2)
+FROM raster_dwithin_rast r1
+CROSS JOIN raster_dwithin_rast r2
+WHERE r1.rid = 0;
+
+SELECT
+       '2.3',
+       r1.rid,
+       r2.rid,
+       ST_DWithin(r1.rast, 1, r2.rast, 1, 7)
+FROM raster_dwithin_rast r1
+CROSS JOIN raster_dwithin_rast r2
+WHERE r1.rid = 0;
+
+DROP TABLE IF EXISTS raster_dwithin_rast;

diff --git a/raster/test/regress/rt_dwithin_expected b/raster/test/regress/rt_dwithin_expected
new file mode 100644 (file)
index 0000000..c7f92d5
--- /dev/null
+++ b/raster/test/regress/rt_dwithin_expected
@@ -0,0 +1,110 @@
+NOTICE:  table "raster_dwithin_rast" does not exist, skipping
+1.1|0|0|
+1.1|0|1|
+1.1|0|2|
+1.1|0|3|
+1.1|0|10|
+1.1|0|11|
+1.1|0|12|
+1.1|0|13|
+1.1|0|14|
+1.1|0|15|
+1.1|0|16|
+1.1|0|20|
+1.1|0|21|
+1.1|0|22|
+1.1|0|23|
+1.1|0|30|
+1.1|0|31|
+1.1|0|32|
+ERROR:  Tolerance cannot be less than zero
+1.3|0|0|t
+1.3|0|1|t
+1.3|0|2|t
+1.3|0|3|f
+1.3|0|10|t
+1.3|0|11|t
+1.3|0|12|t
+1.3|0|13|t
+1.3|0|14|f
+1.3|0|15|t
+1.3|0|16|t
+1.3|0|20|t
+1.3|0|21|t
+1.3|0|22|t
+1.3|0|23|t
+1.3|0|30|t
+1.3|0|31|t
+1.3|0|32|t
+2.1|0|0|t
+2.1|0|1|t
+2.1|0|2|t
+2.1|0|3|f
+2.1|0|10|t
+2.1|0|11|t
+2.1|0|12|t
+2.1|0|13|f
+2.1|0|14|f
+2.1|0|15|t
+2.1|0|16|t
+2.1|0|20|t
+2.1|0|21|t
+2.1|0|22|t
+2.1|0|23|t
+2.1|0|30|t
+2.1|0|31|t
+2.1|0|32|t
+2.2|0|0|t
+2.2|0|1|t
+2.2|0|2|t
+2.2|0|3|f
+2.2|0|10|t
+2.2|0|11|t
+2.2|0|12|t
+2.2|0|13|f
+2.2|0|14|t
+2.2|0|15|t
+2.2|0|16|t
+2.2|0|20|t
+2.2|0|21|t
+2.2|0|22|t
+2.2|0|23|t
+2.2|0|30|t
+2.2|0|31|t
+2.2|0|32|t
+2.3|0|0|t
+2.3|0|1|t
+2.3|0|2|t
+2.3|0|3|f
+2.3|0|10|t
+2.3|0|11|t
+2.3|0|12|t
+2.3|0|13|f
+2.3|0|14|t
+2.3|0|15|t
+2.3|0|16|t
+2.3|0|20|t
+2.3|0|21|t
+2.3|0|22|t
+2.3|0|23|t
+2.3|0|30|t
+2.3|0|31|t
+2.3|0|32|t
+2.3|0|0|t
+2.3|0|1|t
+2.3|0|2|t
+2.3|0|3|t
+2.3|0|10|t
+2.3|0|11|t
+2.3|0|12|t
+2.3|0|13|f
+2.3|0|14|t
+2.3|0|15|t
+2.3|0|16|t
+2.3|0|20|t
+2.3|0|21|t
+2.3|0|22|t
+2.3|0|23|t
+2.3|0|30|t
+2.3|0|31|t
+2.3|0|32|t