Addition of ST_PixelOfValue. Ticket is #1889.

git-svn-id: http://svn.osgeo.org/postgis/trunk@9987 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Bborie Park 2012-06-26 18:17:31 +00:00
parent 1e80f40133
commit 3908d60bc0
10 changed files with 837 additions and 8 deletions

1
NEWS
View file

@ -8,6 +8,7 @@ PostGIS 2.1.0
- #1643, Tiger Geocoder - Tiger 2011 loader (Regina Obe / Paragon Corporation)
Funded by Hunter Systems Group
- GEOMETRYCOLLECTION support for ST_MakeValid (Sandro Santilli / Vizzuality)
- ST_PixelOfValue (Bborie Park / UC Davis)
* Enhancements *
- #823, tiger geocoder: Make loader_generate_script download portion less greedy

View file

@ -3686,6 +3686,139 @@ GROUP BY (foo.geomval).val;
<para><xref linkend="RT_ST_Value" />, <xref linkend="RT_ST_DumpAsPolygons" /></para>
</refsection>
</refentry>
<refentry id="RT_ST_PixelOfValue">
<refnamediv>
<refname>ST_PixelOfValue</refname>
<refpurpose>
Get the columnx, rowy coordinates of the pixel whose value equals the search value.
</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<funcprototype>
<funcdef>setof record <function>ST_PixelOfValue</function></funcdef>
<paramdef>
<type>raster </type> <parameter>rast</parameter>
</paramdef>
<paramdef>
<type>integer </type> <parameter>nband</parameter>
</paramdef>
<paramdef>
<type>double precision[] </type> <parameter>search</parameter>
</paramdef>
<paramdef>
<type>boolean </type> <parameter>exclude_nodata_value=true</parameter>
</paramdef>
</funcprototype>
<funcprototype>
<funcdef>setof record <function>ST_PixelOfValue</function></funcdef>
<paramdef>
<type>raster </type> <parameter>rast</parameter>
</paramdef>
<paramdef>
<type>double precision[] </type> <parameter>search</parameter>
</paramdef>
<paramdef>
<type>boolean </type> <parameter>exclude_nodata_value=true</parameter>
</paramdef>
</funcprototype>
<funcprototype>
<funcdef>setof record <function>ST_PixelOfValue</function></funcdef>
<paramdef>
<type>raster </type> <parameter>rast</parameter>
</paramdef>
<paramdef>
<type>integer </type> <parameter>nband</parameter>
</paramdef>
<paramdef>
<type>double precision </type> <parameter>search</parameter>
</paramdef>
<paramdef>
<type>boolean </type> <parameter>exclude_nodata_value=true</parameter>
</paramdef>
</funcprototype>
<funcprototype>
<funcdef>setof record <function>ST_PixelOfValue</function></funcdef>
<paramdef>
<type>raster </type> <parameter>rast</parameter>
</paramdef>
<paramdef>
<type>double precision </type> <parameter>search</parameter>
</paramdef>
<paramdef>
<type>boolean </type> <parameter>exclude_nodata_value=true</parameter>
</paramdef>
</funcprototype>
</funcsynopsis>
</refsynopsisdiv>
<refsection>
<title>Description</title>
<para>
Get the columnx, rowy coordinates of the pixel whose value equals the search value. If no band is specified, then band 1 is assumed.
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting>
SELECT
(pixels).*
FROM (
SELECT
ST_PixelOfValue(
ST_SetValue(
ST_SetValue(
ST_SetValue(
ST_SetValue(
ST_SetValue(
ST_AddBand(
ST_MakeEmptyRaster(5, 5, -2, 2, 1, -1, 0, 0, 0),
'8BUI'::text, 1, 0
),
1, 1, 0
),
2, 3, 0
),
3, 5, 0
),
4, 2, 0
),
5, 4, 255
)
, 1, ARRAY[1, 255]) AS pixels
) AS foo
val | x | y
-----+---+---
1 | 1 | 2
1 | 1 | 3
1 | 1 | 4
1 | 1 | 5
1 | 2 | 1
1 | 2 | 2
1 | 2 | 4
1 | 2 | 5
1 | 3 | 1
1 | 3 | 2
1 | 3 | 3
1 | 3 | 4
1 | 4 | 1
1 | 4 | 3
1 | 4 | 4
1 | 4 | 5
1 | 5 | 1
1 | 5 | 2
1 | 5 | 3
255 | 5 | 4
1 | 5 | 5
</programlisting>
</refsection>
</refentry>
</sect1>
<sect1 id="Raster_Editors">

View file

@ -2474,6 +2474,78 @@ int rt_band_get_nearest_pixel(
return count;
}
/**
* Search band for pixel(s) with search values
*
* @param band: the band to query for minimum and maximum pixel values
* @param exclude_nodata_value: if non-zero, ignore nodata values
* @param search_values: array of values to count
* @param search_values_count: the number of search values
*
* @return -1 on error, otherwise number of pixels
*/
int
rt_band_get_pixel_of_value(
rt_band band, int exclude_nodata_value,
double *searchset, int searchcount,
rt_pixel *pixels
) {
int x;
int y;
int i;
double pixval;
int err;
int count = 0;
rt_pixel pixel = NULL;
assert(NULL != band);
assert(NULL != pixels);
for (x = 0; x < band->width; x++) {
for (y = 0; y < band->height; y++) {
err = rt_band_get_pixel(band, x, y, &pixval);
if (err != 0) {
rterror("rt_band_get_pixel_of_value: Cannot get band pixel");
return -1;
}
else if (
exclude_nodata_value &&
(band->hasnodata != FALSE) && (
FLT_EQ(pixval, band->nodataval) ||
rt_band_clamped_value_is_nodata(band, pixval) == 1
)
) {
continue;
}
for (i = 0; i < searchcount; i++) {
if (FLT_NEQ(pixval, searchset[i]))
continue;
/* match found */
count++;
if (*pixels == NULL)
*pixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count);
else
*pixels = (rt_pixel) rtrealloc(*pixels, sizeof(struct rt_pixel_t) * count);
if (*pixels == NULL) {
rterror("rt_band_get_pixel_of_value: Unable to allocate memory for pixel(s)");
return -1;
}
pixel = &((*pixels)[count - 1]);
pixel->x = x;
pixel->y = y;
pixel->nodata = 0;
pixel->value = pixval;
}
}
}
return count;
}
double
rt_band_get_nodata(rt_band band) {
@ -2717,7 +2789,6 @@ rt_band_get_summary_stats(
int exclude_nodata_value, double sample, int inc_vals,
uint64_t *cK, double *cM, double *cQ
) {
uint8_t *data = NULL;
uint32_t x = 0;
uint32_t y = 0;
uint32_t z = 0;
@ -2775,12 +2846,6 @@ rt_band_get_summary_stats(
return stats;
}
data = rt_band_get_data(band);
if (data == NULL) {
rterror("rt_band_get_summary_stats: Cannot get band data");
return NULL;
}
hasnodata = rt_band_get_hasnodata_flag(band);
if (hasnodata != FALSE)
nodata = rt_band_get_nodata(band);

View file

@ -582,6 +582,22 @@ int rt_band_get_nearest_pixel(
rt_pixel *npixels
);
/**
* Search band for pixel(s) with search values
*
* @param band: the band to query for minimum and maximum pixel values
* @param exclude_nodata_value: if non-zero, ignore nodata values
* @param search_values: array of values to count
* @param search_values_count: the number of search values
*
* @return -1 on error, otherwise number of pixels
*/
int rt_band_get_pixel_of_value(
rt_band band, int exclude_nodata_value,
double *searchset, int searchcount,
rt_pixel *pixels
);
/**
* Returns the minimal possible value for the band according to the pixel type.
* @param band: the band to get info from

View file

@ -202,6 +202,9 @@ Datum RASTER_setPixelValue(PG_FUNCTION_ARGS);
/* Get pixel geographical shape */
Datum RASTER_getPixelPolygon(PG_FUNCTION_ARGS);
/* Get pixels of value */
Datum RASTER_pixelOfValue(PG_FUNCTION_ARGS);
/* Get nearest value to a point */
Datum RASTER_nearestValue(PG_FUNCTION_ARGS);
@ -2497,6 +2500,241 @@ Datum RASTER_getPixelPolygon(PG_FUNCTION_ARGS)
PG_RETURN_POINTER(gser);
}
/**
* Get pixels of value
*/
PG_FUNCTION_INFO_V1(RASTER_pixelOfValue);
Datum RASTER_pixelOfValue(PG_FUNCTION_ARGS)
{
FuncCallContext *funcctx;
TupleDesc tupdesc;
rt_pixel pixels = NULL;
rt_pixel pixels2 = NULL;
int count = 0;
int i = 0;
int n = 0;
int call_cntr;
int max_calls;
if (SRF_IS_FIRSTCALL()) {
MemoryContext oldcontext;
rt_pgraster *pgraster = NULL;
rt_raster raster = NULL;
rt_band band = NULL;
int nband = 1;
int num_bands = 0;
double *search = NULL;
int nsearch = 0;
double val;
bool exclude_nodata_value = TRUE;
ArrayType *array;
Oid etype;
Datum *e;
bool *nulls;
int16 typlen;
bool typbyval;
char typalign;
int ndims = 1;
int *dims;
int *lbs;
/* create a function context for cross-call persistence */
funcctx = SRF_FIRSTCALL_INIT();
/* switch to memory context appropriate for multiple function calls */
oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
if (PG_ARGISNULL(0)) {
MemoryContextSwitchTo(oldcontext);
SRF_RETURN_DONE(funcctx);
}
pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
raster = rt_raster_deserialize(pgraster, FALSE);
if (!raster) {
elog(ERROR, "RASTER_pixelOfValue: Could not deserialize raster");
PG_FREE_IF_COPY(pgraster, 0);
MemoryContextSwitchTo(oldcontext);
SRF_RETURN_DONE(funcctx);
}
/* num_bands */
num_bands = rt_raster_get_num_bands(raster);
if (num_bands < 1) {
elog(NOTICE, "Raster provided has no bands");
rt_raster_destroy(raster);
PG_FREE_IF_COPY(pgraster, 0);
MemoryContextSwitchTo(oldcontext);
SRF_RETURN_DONE(funcctx);
}
/* band index is 1-based */
if (!PG_ARGISNULL(1))
nband = PG_GETARG_INT32(1);
if (nband < 1 || nband > num_bands) {
elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
rt_raster_destroy(raster);
PG_FREE_IF_COPY(pgraster, 0);
MemoryContextSwitchTo(oldcontext);
SRF_RETURN_DONE(funcctx);
}
/* search values */
array = PG_GETARG_ARRAYTYPE_P(2);
etype = ARR_ELEMTYPE(array);
get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
switch (etype) {
case FLOAT4OID:
case FLOAT8OID:
break;
default:
elog(ERROR, "RASTER_pixelOfValue: Invalid data type for pixel values");
rt_raster_destroy(raster);
PG_FREE_IF_COPY(pgraster, 0);
MemoryContextSwitchTo(oldcontext);
SRF_RETURN_DONE(funcctx);
break;
}
ndims = ARR_NDIM(array);
dims = ARR_DIMS(array);
lbs = ARR_LBOUND(array);
deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
&nulls, &n);
search = palloc(sizeof(double) * n);
for (i = 0, nsearch = 0; i < n; i++) {
if (nulls[i]) continue;
switch (etype) {
case FLOAT4OID:
val = (double) DatumGetFloat4(e[i]);
break;
case FLOAT8OID:
val = (double) DatumGetFloat8(e[i]);
break;
}
search[nsearch] = val;
POSTGIS_RT_DEBUGF(3, "search[%d] = %d", nsearch, search[nsearch]);
nsearch++;
}
/* not searching for anything */
if (nsearch < 1) {
elog(NOTICE, "No search values provided. Returning NULL");
pfree(search);
rt_raster_destroy(raster);
PG_FREE_IF_COPY(pgraster, 0);
MemoryContextSwitchTo(oldcontext);
SRF_RETURN_DONE(funcctx);
}
else if (nsearch < n)
search = repalloc(search, sizeof(double) * nsearch);
/* exclude_nodata_value flag */
if (!PG_ARGISNULL(3))
exclude_nodata_value = PG_GETARG_BOOL(3);
/* get band */
band = rt_raster_get_band(raster, nband - 1);
if (!band) {
elog(NOTICE, "Could not find band at index %d. Returning NULL", nband);
rt_raster_destroy(raster);
PG_FREE_IF_COPY(pgraster, 0);
MemoryContextSwitchTo(oldcontext);
SRF_RETURN_DONE(funcctx);
}
/* get pixels of values */
count = rt_band_get_pixel_of_value(
band, exclude_nodata_value,
search, nsearch,
&pixels
);
pfree(search);
rt_band_destroy(band);
rt_raster_destroy(raster);
PG_FREE_IF_COPY(pgraster, 0);
if (count < 1) {
/* error */
if (count < 0)
elog(NOTICE, "Unable to get the pixels of search values for band at index %d", nband);
/* no nearest pixel */
else
elog(NOTICE, "No pixels of search values found for band at index %d", nband);
MemoryContextSwitchTo(oldcontext);
SRF_RETURN_DONE(funcctx);
}
/* Store needed information */
funcctx->user_fctx = pixels;
/* total number of tuples to be returned */
funcctx->max_calls = count;
/* Build a tuple descriptor for our result type */
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
ereport(ERROR, (
errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg(
"function returning record called in context "
"that cannot accept type record"
)
));
}
BlessTupleDesc(tupdesc);
funcctx->tuple_desc = tupdesc;
MemoryContextSwitchTo(oldcontext);
}
/* stuff done on every call of the function */
funcctx = SRF_PERCALL_SETUP();
call_cntr = funcctx->call_cntr;
max_calls = funcctx->max_calls;
tupdesc = funcctx->tuple_desc;
pixels2 = funcctx->user_fctx;
/* do when there is more left to send */
if (call_cntr < max_calls) {
int values_length = 3;
Datum values[values_length];
bool nulls[values_length];
HeapTuple tuple;
Datum result;
memset(nulls, FALSE, values_length);
/* 0-based to 1-based */
pixels2[call_cntr].x += 1;
pixels2[call_cntr].y += 1;
values[0] = Float8GetDatum(pixels2[call_cntr].value);
values[1] = Int64GetDatum(pixels2[call_cntr].x);
values[2] = Int64GetDatum(pixels2[call_cntr].y);
/* build a tuple */
tuple = heap_form_tuple(tupdesc, values, nulls);
/* make the tuple into a datum */
result = HeapTupleGetDatum(tuple);
SRF_RETURN_NEXT(funcctx, result);
}
else {
pfree(pixels2);
SRF_RETURN_DONE(funcctx);
}
}
/**
* Return nearest value to a point
*/

View file

@ -2208,6 +2208,58 @@ CREATE OR REPLACE FUNCTION st_value(rast raster, pt geometry, hasnodata boolean
AS $$ SELECT st_value($1, 1, $2, $3) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
-----------------------------------------------------------------------
-- ST_PixelOfValue()
-----------------------------------------------------------------------
CREATE OR REPLACE FUNCTION st_pixelofvalue(
rast raster,
nband integer,
search double precision[],
exclude_nodata_value boolean DEFAULT TRUE,
OUT val double precision,
OUT x integer,
OUT y integer
)
RETURNS SETOF record
AS 'MODULE_PATHNAME', 'RASTER_pixelOfValue'
LANGUAGE 'c' IMMUTABLE STRICT;
CREATE OR REPLACE FUNCTION st_pixelofvalue(
rast raster,
search double precision[],
exclude_nodata_value boolean DEFAULT TRUE,
OUT val double precision,
OUT x integer,
OUT y integer
)
RETURNS SETOF record
AS $$ SELECT val, x, y FROM st_pixelofvalue($1, 1, $2, $3) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
CREATE OR REPLACE FUNCTION st_pixelofvalue(
rast raster,
nband integer,
search double precision,
exclude_nodata_value boolean DEFAULT TRUE,
OUT x integer,
OUT y integer
)
RETURNS SETOF record
AS $$ SELECT x, y FROM st_pixelofvalue($1, $2, ARRAY[$3], $4) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
CREATE OR REPLACE FUNCTION st_pixelofvalue(
rast raster,
search double precision,
exclude_nodata_value boolean DEFAULT TRUE,
OUT x integer,
OUT y integer
)
RETURNS SETOF record
AS $$ SELECT x, y FROM st_pixelofvalue($1, 1, ARRAY[$2], $3) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
-----------------------------------------------------------------------
-- Raster Accessors ST_Georeference()
-----------------------------------------------------------------------

View file

@ -2813,6 +2813,87 @@ static void testNearestPixel() {
deepRelease(rast);
}
static void testPixelOfValue() {
rt_raster rast;
rt_band band;
uint32_t x, y;
int rtn;
const int maxX = 10;
const int maxY = 10;
rt_pixel pixels = NULL;
double search0[1] = {0};
double search1[1] = {1};
double search2[2] = {3, 5};
rast = rt_raster_new(maxX, maxY);
assert(rast);
band = addBand(rast, PT_32BUI, 1, 0);
CHECK(band);
for (x = 0; x < maxX; x++) {
for (y = 0; y < maxY; y++) {
rtn = rt_band_set_pixel(band, x, y, 1);
CHECK((rtn != -1));
}
}
rt_band_set_pixel(band, 0, 0, 0);
rt_band_set_pixel(band, 3, 0, 0);
rt_band_set_pixel(band, 6, 0, 0);
rt_band_set_pixel(band, 9, 0, 0);
rt_band_set_pixel(band, 1, 2, 0);
rt_band_set_pixel(band, 4, 2, 0);
rt_band_set_pixel(band, 7, 2, 0);
rt_band_set_pixel(band, 2, 4, 0);
rt_band_set_pixel(band, 5, 4, 0);
rt_band_set_pixel(band, 8, 4, 0);
rt_band_set_pixel(band, 0, 6, 0);
rt_band_set_pixel(band, 3, 6, 0);
rt_band_set_pixel(band, 6, 6, 0);
rt_band_set_pixel(band, 9, 6, 0);
rt_band_set_pixel(band, 1, 8, 0);
rt_band_set_pixel(band, 4, 8, 0);
rt_band_set_pixel(band, 7, 8, 0);
pixels = NULL;
rtn = rt_band_get_pixel_of_value(
band, TRUE,
search1, 1,
&pixels
);
CHECK((rtn == 83));
if (rtn)
rtdealloc(pixels);
pixels = NULL;
rtn = rt_band_get_pixel_of_value(
band, FALSE,
search0, 1,
&pixels
);
CHECK((rtn == 17));
if (rtn)
rtdealloc(pixels);
rt_band_set_pixel(band, 4, 2, 3);
rt_band_set_pixel(band, 7, 2, 5);
rt_band_set_pixel(band, 1, 8, 3);
pixels = NULL;
rtn = rt_band_get_pixel_of_value(
band, TRUE,
search2, 2,
&pixels
);
CHECK((rtn == 3));
if (rtn)
rtdealloc(pixels);
deepRelease(rast);
}
int
main()
{
@ -3083,6 +3164,10 @@ main()
testNearestPixel();
printf("OK\n");
printf("Testing rt_band_get_pixel_of_value... ");
testPixelOfValue();
printf("OK\n");
deepRelease(raster);
return EXIT_SUCCESS;

View file

@ -81,7 +81,8 @@ TEST_BANDPROPS = \
rt_pixelvalue \
drop_rt_band_properties_test \
rt_neighborhood \
rt_nearestvalue
rt_nearestvalue \
rt_pixelofvalue
TEST_UTILITY = \
rt_utility \

View file

@ -0,0 +1,113 @@
DROP TABLE IF EXISTS raster_pixelofvalue;
CREATE TABLE raster_pixelofvalue (
rast raster
);
CREATE OR REPLACE FUNCTION make_test_raster()
RETURNS void
AS $$
DECLARE
width int := 10;
height int := 10;
x int;
y int;
rast raster;
BEGIN
rast := ST_MakeEmptyRaster(width, height, 0, 0, 1, -1, 0, 0, 0);
rast := ST_AddBand(rast, 1, '32BUI', 0, 0);
FOR x IN 1..width LOOP
FOR y IN 1..height LOOP
IF (x + y) % 2 = 1 THEN
rast := ST_SetValue(rast, 1, x, y, x + y);
END IF;
END LOOP;
END LOOP;
INSERT INTO raster_pixelofvalue VALUES (rast);
RETURN;
END;
$$ LANGUAGE 'plpgsql';
SELECT make_test_raster();
DROP FUNCTION make_test_raster();
SELECT
(pixval).*
FROM (
SELECT
ST_PixelOfValue(
rast, 1,
ARRAY[3, 11]
) AS pixval
FROM raster_pixelofvalue
) foo;
SELECT
(pixval).*
FROM (
SELECT
ST_PixelOfValue(
rast, 1,
ARRAY[5]
) AS pixval
FROM raster_pixelofvalue
) foo;
SELECT
(pixval).*
FROM (
SELECT
ST_PixelOfValue(
rast,
ARRAY[0]
) AS pixval
FROM raster_pixelofvalue
) foo;
SELECT
(pixval).*
FROM (
SELECT
ST_PixelOfValue(
rast,
ARRAY[0],
FALSE
) AS pixval
FROM raster_pixelofvalue
) foo;
SELECT
(pixval).*
FROM (
SELECT
ST_PixelOfValue(
rast, 1,
7
) AS pixval
FROM raster_pixelofvalue
) foo;
SELECT
(pixval).*
FROM (
SELECT
ST_PixelOfValue(
rast,
2
) AS pixval
FROM raster_pixelofvalue
) foo;
SELECT
(pixval).*
FROM (
SELECT
ST_PixelOfValue(
rast,
0,
FALSE
) AS pixval
FROM raster_pixelofvalue
) foo;
DROP TABLE IF EXISTS raster_pixelofvalue;

View file

@ -0,0 +1,125 @@
NOTICE: table "raster_pixelofvalue" does not exist, skipping
3|1|2
11|1|10
3|2|1
11|2|9
11|3|8
11|4|7
11|5|6
11|6|5
11|7|4
11|8|3
11|9|2
11|10|1
5|1|4
5|2|3
5|3|2
5|4|1
NOTICE: No pixels of search values found for band at index 1
0|1|1
0|1|3
0|1|5
0|1|7
0|1|9
0|2|2
0|2|4
0|2|6
0|2|8
0|2|10
0|3|1
0|3|3
0|3|5
0|3|7
0|3|9
0|4|2
0|4|4
0|4|6
0|4|8
0|4|10
0|5|1
0|5|3
0|5|5
0|5|7
0|5|9
0|6|2
0|6|4
0|6|6
0|6|8
0|6|10
0|7|1
0|7|3
0|7|5
0|7|7
0|7|9
0|8|2
0|8|4
0|8|6
0|8|8
0|8|10
0|9|1
0|9|3
0|9|5
0|9|7
0|9|9
0|10|2
0|10|4
0|10|6
0|10|8
0|10|10
1|6
2|5
3|4
4|3
5|2
6|1
NOTICE: No pixels of search values found for band at index 1
1|1
1|3
1|5
1|7
1|9
2|2
2|4
2|6
2|8
2|10
3|1
3|3
3|5
3|7
3|9
4|2
4|4
4|6
4|8
4|10
5|1
5|3
5|5
5|7
5|9
6|2
6|4
6|6
6|8
6|10
7|1
7|3
7|5
7|7
7|9
8|2
8|4
8|6
8|8
8|10
9|1
9|3
9|5
9|7
9|9
10|2
10|4
10|6
10|8
10|10