From f3875ab86aca42125b79d552d2fafecbdde33205 Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Wed, 25 Jul 2012 03:35:45 +0000 Subject: [PATCH] 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 --- raster/rt_core/rt_api.c | 2 + raster/rt_pg/rt_pg.c | 144 ++++++++++++++++ raster/rt_pg/rtpostgis.sql.in.c | 22 +++ raster/test/regress/Makefile.in | 1 + raster/test/regress/rt_dwithin.sql | 218 ++++++++++++++++++++++++ raster/test/regress/rt_dwithin_expected | 110 ++++++++++++ 6 files changed, 497 insertions(+) create mode 100644 raster/test/regress/rt_dwithin.sql create mode 100644 raster/test/regress/rt_dwithin_expected diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index c0d3def85..e418f3950 100644 --- 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 eaa40f436..87933b528 100644 --- 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 811d59694..d2893555f 100644 --- 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 9e83d3f08..0917e0e0f 100644 --- 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 index 000000000..2d83950b4 --- /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 index 000000000..c7f92d524 --- /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 -- 2.40.0