From: Bborie Park Date: Mon, 16 May 2011 19:36:55 +0000 (+0000) Subject: Added ST_Histogram functions. X-Git-Tag: 2.0.0alpha1~1659 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=cbaaeda73bee8d13094c69254686544a1a6d76aa;p=postgis Added ST_Histogram functions. - added function rt_band_get_histogram to rt_core/rt_api.c and rt_api.h - added test case to test/core/testapi.c - added function RASTER_histogram to rt_pg/rt_pg.c - added SQL functions for ST_Histogram - added regression tests Associated ticket is #934 git-svn-id: http://svn.osgeo.org/postgis/trunk@7152 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index cb01a5430..ac126ab3e 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -1705,6 +1705,249 @@ rt_band_get_summary_stats(rt_band band, int hasnodata, double sample, return stats; } +struct rt_histogram_t { + uint32_t count; + double proportion; + + double min; + double max; + + int inc_min; + int inc_max; +}; + +/** + * Count the distribution of data + * + * @param stats: a populated stats struct for processing + * @param bin_count: the number of bins to group the data by + * @param bin_width: the width of each bin as an array + * @param bin_width_count: number of values in bin_width + * @param right: evaluate bins by (a,b] rather than default [a,b) + * @param rtn_count: set to the number of bins being returned + * + * @return the histogram of the data + */ +rt_histogram +rt_band_get_histogram(rt_bandstats stats, + int bin_count, double *bin_width, int bin_width_count, + int right, int *rtn_count) { + rt_histogram bins = NULL; + int init_width = 0; + int i; + int j; + double tmp; + double value; + int sum = 0; + +#if POSTGIS_DEBUG_LEVEL > 0 + clock_t start, stop; + double elapsed = 0; +#endif + + RASTER_DEBUG(3, "starting"); +#if POSTGIS_DEBUG_LEVEL > 0 + start = clock(); +#endif + + assert(NULL != stats); + + if (stats->count < 1 || NULL == stats->values) { + rterror("rt_util_get_histogram: rt_bandstats object has no value"); + return NULL; + } + + /* bin width must be positive numbers and not zero */ + if (NULL != bin_width && bin_width_count > 0) { + for (i = 0; i < bin_width_count; i++) { + if (bin_width[i] <= 0.) { + rterror("rt_util_get_histogram: bin_width element is less than or equal to zero"); + return NULL; + } + } + } + + /* # of bins not provided */ + if (bin_count <= 0) { + /* + determine # of bins + http://en.wikipedia.org/wiki/Histogram + + all computed bins are assumed to have equal width + */ + /* Square-root choice for stats->count < 30 */ + if (stats->count < 30) + bin_count = ceil(sqrt(stats->count)); + /* Sturges' formula for stats->count >= 30 */ + else + bin_count = ceil(log2((double) stats->count) + 1.); + + /* bin_width_count provided and bin_width has value */ + if (bin_width_count > 0 && NULL != bin_width) { + /* user has defined something specific */ + if (bin_width_count > bin_count) + bin_count = bin_width_count; + else if (bin_width_count > 1) { + tmp = 0; + for (i = 0; i < bin_width_count; i++) tmp += bin_width[i]; + bin_count = ceil((stats->max - stats->min) / tmp) * bin_width_count; + } + else + bin_count = ceil((stats->max - stats->min) / bin_width[0]); + } + /* set bin width count to zero so that one can be calculated */ + else { + bin_width_count = 0; + } + } + + /* min and max the same */ + if (stats->min == stats->max) + bin_count = 1; + + RASTER_DEBUGF(3, "bin_count = %d", bin_count); + + /* bin count = 1, all values are in one bin */ + if (bin_count < 2) { + bins = rtalloc(sizeof(struct rt_histogram_t)); + if (NULL == bins) { + rterror("rt_util_get_histogram: Unable to allocate memory for histogram"); + return NULL; + } + + bins->count = stats->count; + bins->proportion = -1; + bins->min = stats->min; + bins->max = stats->max; + bins->inc_min = bins->inc_max = 1; + + *rtn_count = bin_count; + return bins; + } + + /* establish bin width */ + if (bin_width_count == 0) { + bin_width_count = 1; + + /* bin_width unallocated */ + if (NULL == bin_width) { + bin_width = rtalloc(sizeof(double)); + if (NULL == bin_width) { + rterror("rt_util_get_histogram: Unable to allocate memory for bin widths"); + return NULL; + } + init_width = 1; + } + + bin_width[0] = (stats->max - stats->min) / bin_count; + } + + /* initialize bins */ + bins = rtalloc(bin_count * sizeof(struct rt_histogram_t)); + if (NULL == bins) { + rterror("rt_util_get_histogram: Unable to allocate memory for histogram"); + if (init_width) rtdealloc(bin_width); + return NULL; + } + if (!right) + tmp = stats->min; + else + tmp = stats->max; + for (i = 0; i < bin_count;) { + for (j = 0; j < bin_width_count; j++) { + bins[i].count = 0; + bins->proportion = -1; + + if (!right) { + bins[i].min = tmp; + tmp += bin_width[j]; + bins[i].max = tmp; + + bins[i].inc_min = 1; + bins[i].inc_max = 0; + } + else { + bins[i].max = tmp; + tmp -= bin_width[j]; + bins[i].min = tmp; + + bins[i].inc_min = 0; + bins[i].inc_max = 1; + } + + i++; + } + } + if (!right) { + bins[bin_count - 1].inc_max = 1; + + /* align last bin to the max value */ + if (bins[bin_count - 1].max < stats->max) + bins[bin_count - 1].max = stats->max; + } + else { + bins[bin_count - 1].inc_min = 1; + + /* align first bin to the min value */ + if (bins[bin_count - 1].min > stats->min) + bins[bin_count - 1].min = stats->min; + } + + /* process the values */ + for (i = 0; i < stats->count; i++) { + value = stats->values[i]; + + /* default, [a, b) */ + if (!right) { + for (j = 0; j < bin_count; j++) { + if ( + (!bins[j].inc_max && value < bins[j].max) || + (bins[j].inc_max && value <= bins[j].max) + ) { + bins[j].count++; + sum++; + break; + } + } + } + else { + for (j = 0; j < bin_count; j++) { + if ( + (!bins[j].inc_min && value > bins[j].min) || + (bins[j].inc_min && value >= bins[j].min) + ) { + bins[j].count++; + sum++; + break; + } + } + } + } + + /* proportions */ + for (i = 0; i < bin_count; i++) { + bins[i].proportion = ((double) bins[i].count) / sum; + if (bin_width_count > 1) + bins[i].proportion /= (bins[i].max - bins[i].min); + } + +#if POSTGIS_DEBUG_LEVEL > 0 + stop = clock(); + elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC; + RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed); + + for (j = 0; j < bin_count; j++) { + RASTER_DEBUGF(5, "(min, max, inc_min, inc_max, count, sum, proportion) = (%f, %f, %d, %d, %d, %d, %f)", + bins[j].min, bins[j].max, bins[j].inc_min, bins[j].inc_max, bins[j].count, sum, bins[j].proportion); + } +#endif + + if (init_width) rtdealloc(bin_width); + *rtn_count = bin_count; + RASTER_DEBUG(3, "done"); + return bins; +} + /*- rt_raster --------------------------------------------------------*/ struct rt_raster_serialized_t { diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index 968294bca..b975ba373 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -457,6 +457,23 @@ typedef struct rt_bandstats_t* rt_bandstats; rt_bandstats rt_band_get_summary_stats(rt_band band, int hasnodata, double sample, int inc_vals); +/** + * Count the distribution of data + * + * @param stats: a populated stats struct for processing + * @param bin_count: the number of bins to group the data by + * @param bin_width: the width of each bin as an array + * @param bin_width_count: number of values in bin_width + * @param right: evaluate bins by (a,b] rather than default [a,b) + * @param rtn_count: set to the number of bins being returned + * + * @return the histogram of the data + */ +typedef struct rt_histogram_t* rt_histogram; +rt_histogram rt_band_get_histogram(rt_bandstats stats, + int bin_count, double *bin_widths, int bin_widths_count, + int right, int *rtn_count); + /*- rt_raster --------------------------------------------------------*/ diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index 4b169c8ed..3e2adb8dc 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -190,6 +190,9 @@ Datum RASTER_band(PG_FUNCTION_ARGS); /* Get summary stats */ Datum RASTER_summaryStats(PG_FUNCTION_ARGS); +/* get histogram */ +Datum RASTER_histogram(PG_FUNCTION_ARGS); + /* Replace function taken from * http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3 @@ -2961,6 +2964,310 @@ Datum RASTER_summaryStats(PG_FUNCTION_ARGS) } } +/* get histogram */ +struct rt_histogram_t { + uint32_t count; + double proportion; + + double min; + double max; + + int inc_min; + int inc_max; +}; + +/** + * Returns histogram for a band + */ +PG_FUNCTION_INFO_V1(RASTER_histogram); +Datum RASTER_histogram(PG_FUNCTION_ARGS) +{ + FuncCallContext *funcctx; + TupleDesc tupdesc; + AttInMetadata *attinmeta; + + int count; + rt_histogram hist; + rt_histogram hist2; + int call_cntr; + int max_calls; + + /* first call of function */ + if (SRF_IS_FIRSTCALL()) { + MemoryContext oldcontext; + + rt_pgraster *pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + rt_raster raster = NULL; + rt_band band = NULL; + int32_t bandindex = 0; + int num_bands = 0; + bool hasnodata = FALSE; + double sample = 0; + int bin_count = 0; + double *bin_width = NULL; + int bin_width_count = 0; + double width; + bool right = FALSE; + rt_bandstats stats = NULL; + + int i; + int j; + int n; + + 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); + + raster = rt_raster_deserialize(pgraster); + if (!raster) { + elog(ERROR, "RASTER_histogram: Could not deserialize raster"); + PG_RETURN_NULL(); + } + + /* band index is 1-based */ + bandindex = PG_GETARG_INT32(1); + num_bands = rt_raster_get_num_bands(raster); + if (bandindex < 1 || bandindex > num_bands) { + elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL"); + rt_raster_destroy(raster); + SRF_RETURN_DONE(funcctx); + } + assert(0 <= (bandindex - 1)); + + /* hasnodata flag */ + if (!PG_ARGISNULL(2)) + hasnodata = PG_GETARG_BOOL(2); + + /* sample % */ + if (!PG_ARGISNULL(3)) { + sample = PG_GETARG_FLOAT8(3); + if (sample < 0 || sample > 1) { + elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL"); + rt_raster_destroy(raster); + SRF_RETURN_DONE(funcctx); + } + else if (sample == 0) + sample = 1; + } + else + sample = 1; + + /* bin_count */ + if (!PG_ARGISNULL(4)) { + bin_count = PG_GETARG_INT32(4); + if (bin_count < 1) bin_count = 0; + } + + /* bin_width */ + if (!PG_ARGISNULL(5)) { + array = PG_GETARG_ARRAYTYPE_P(5); + etype = ARR_ELEMTYPE(array); + get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign); + + switch (etype) { + case FLOAT4OID: + case FLOAT8OID: + break; + default: + elog(ERROR, "RASTER_histogram: Invalid data type for width"); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + break; + } + + ndims = ARR_NDIM(array); + dims = ARR_DIMS(array); + lbs = ARR_LBOUND(array); + + deconstruct_array(array, etype, typlen, typbyval, typalign, &e, + &nulls, &n); + + bin_width = palloc(sizeof(double) * n); + for (i = 0, j = 0; i < n; i++) { + if (nulls[i]) continue; + + switch (etype) { + case FLOAT4OID: + width = (double) DatumGetFloat4(e[i]); + break; + case FLOAT8OID: + width = (double) DatumGetFloat8(e[i]); + break; + } + + if (width < 0) { + elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL"); + pfree(bin_width); + rt_raster_destroy(raster); + SRF_RETURN_DONE(funcctx); + } + + bin_width[j] = width; + POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]); + j++; + } + bin_width_count = j; + + if (j < 1) { + pfree(bin_width); + bin_width = NULL; + } + } + + /* right */ + if (!PG_ARGISNULL(6)) + right = PG_GETARG_BOOL(6); + + /* get band */ + band = rt_raster_get_band(raster, bandindex - 1); + if (!band) { + elog(NOTICE, "Could not find raster band of index %d. Returning NULL", bandindex); + rt_raster_destroy(raster); + SRF_RETURN_DONE(funcctx); + } + + /* get stats */ + stats = rt_band_get_summary_stats(band, (int) hasnodata, sample, 1); + rt_band_destroy(band); + rt_raster_destroy(raster); + if (NULL == stats || NULL == stats->values) { + elog(NOTICE, "Could not retrieve summary statistics of raster band of index %d", bandindex); + SRF_RETURN_DONE(funcctx); + } + + /* get histogram */ + hist = rt_band_get_histogram(stats, bin_count, bin_width, bin_width_count, right, &count); + if (bin_width_count) pfree(bin_width); + pfree(stats); + if (NULL == hist || !count) { + elog(NOTICE, "Could not retrieve histogram of raster band of index %d", bandindex); + SRF_RETURN_DONE(funcctx); + } + + POSTGIS_RT_DEBUGF(3, "%d bins returned", count); + + /* Store needed information */ + funcctx->user_fctx = hist; + + /* 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" + ) + )); + } + + /* + * generate attribute metadata needed later to produce tuples from raw + * C strings + */ + attinmeta = TupleDescGetAttInMetadata(tupdesc); + funcctx->attinmeta = attinmeta; + MemoryContextSwitchTo(oldcontext); + } + + /* stuff done on every call of the function */ + funcctx = SRF_PERCALL_SETUP(); + + call_cntr = funcctx->call_cntr; + max_calls = funcctx->max_calls; + attinmeta = funcctx->attinmeta; + hist2 = funcctx->user_fctx; + + /* do when there is more left to send */ + if (call_cntr < max_calls) { + char **values; + HeapTuple tuple; + Datum result; + + POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr); + + /* + * Prepare a values array for building the returned tuple. + * This should be an array of C strings which will + * be processed later by the type input functions. + */ + values = (char **) palloc(4 * sizeof(char *)); + + values[0] = (char *) palloc( + (snprintf(NULL, 0, "%f", hist2[call_cntr].min) + 1) * sizeof(char) + ); + values[1] = (char *) palloc( + (snprintf(NULL, 0, "%f", hist2[call_cntr].max) + 1) * sizeof(char) + ); + values[2] = (char *) palloc( + (snprintf(NULL, 0, "%d", hist2[call_cntr].count) + 1) * sizeof(char) + ); + values[3] = (char *) palloc( + (snprintf(NULL, 0, "%f", hist2[call_cntr].proportion) + 1) * sizeof(char) + ); + + snprintf( + values[0], + (snprintf(NULL, 0, "%f", hist2[call_cntr].min) + 1) * sizeof(char), + "%f", + hist2[call_cntr].min + ); + snprintf( + values[1], + (snprintf(NULL, 0, "%f", hist2[call_cntr].max) + 1) * sizeof(char), + "%f", + hist2[call_cntr].max + ); + snprintf( + values[2], + (snprintf(NULL, 0, "%d", hist2[call_cntr].count) + 1) * sizeof(char), + "%d", + hist2[call_cntr].count + ); + snprintf( + values[3], + (snprintf(NULL, 0, "%f", hist2[call_cntr].proportion) + 1) * sizeof(char), + "%f", + hist2[call_cntr].proportion + ); + + /* build a tuple */ + tuple = BuildTupleFromCStrings(attinmeta, values); + + /* make the tuple into a datum */ + result = HeapTupleGetDatum(tuple); + + /* clean up (this is not really necessary) */ + pfree(values[3]); + pfree(values[2]); + pfree(values[1]); + pfree(values[0]); + pfree(values); + + SRF_RETURN_NEXT(funcctx, result); + } + /* do when there is no more left */ + else { + pfree(hist2); + SRF_RETURN_DONE(funcctx); + } +} + /* ---------------------------------------------------------------- */ /* Memory allocation / error reporting hooks */ /* ---------------------------------------------------------------- */ diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index 2191ee85d..e72cee456 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -728,6 +728,127 @@ CREATE OR REPLACE FUNCTION st_approxminmax(rastertable text, rastercolumn text, AS $$ SELECT min, max FROM _st_summarystats($1, $2, 1, FALSE, 0.1) $$ LANGUAGE 'SQL' IMMUTABLE STRICT; +----------------------------------------------------------------------- +-- ST_Histogram and ST_ApproxHistogram +----------------------------------------------------------------------- +CREATE TYPE histogram AS ( + min double precision, + max double precision, + count integer, + proportion double precision +); + +-- Cannot be strict as "width" can be NULL +CREATE OR REPLACE FUNCTION _st_histogram(rast raster, nband int, hasnodata boolean, sample_percent double precision, bins int, width double precision[], right boolean) + RETURNS SETOF histogram + AS 'MODULE_PATHNAME','RASTER_histogram' + LANGUAGE 'C' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_histogram(rast raster, nband int, hasnodata boolean, bins int, width double precision[], right boolean) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, $3, 1, $4, $5, $6) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_histogram(rast raster, nband int, hasnodata boolean, bins int, right boolean) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, $3, 1, $4, NULL, $5) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_histogram(rast raster, nband int, hasnodata boolean, bins int) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, $3, 1, $4, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_histogram(rast raster, nband int, hasnodata boolean) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, $3, 1, 0, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_histogram(rast raster, nband int) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, 1, 0, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_histogram(rast raster) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, 1, FALSE, 1, 0, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_histogram(rast raster, nband int, bins int, width double precision[], right boolean) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, 1, $3, $4, $5) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_histogram(rast raster, nband int, bins int, right boolean) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, 1, $3, NULL, $4) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_histogram(rast raster, nband int, bins int) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, 1, $3, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, hasnodata boolean, sample_percent double precision, bins int, width double precision[], right boolean) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, $3, $4, $5, $6, $7) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, hasnodata boolean, sample_percent double precision, bins int, right boolean) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, $3, $4, $5, NULL, $6) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, hasnodata boolean, sample_percent double precision, bins int) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, $3, $4, $5, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, hasnodata boolean, sample_percent double precision) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, $3, $4, 0, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, sample_percent double precision) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, $3, 0, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, 0.1, 0, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, sample_percent double precision) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, 1, FALSE, $2, 0, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, 1, FALSE, 0.1, 0, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, sample_percent double precision, bins int, width double precision[], right boolean) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, $3, $4, $5, $6) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, sample_percent double precision, bins int, right boolean) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, $3, $4, NULL, $5) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, sample_percent double precision, bins int) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, $3, $4, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_approxhistogram(rast raster, nband int, sample_percent double precision) + RETURNS SETOF histogram + AS $$ SELECT min, max, count, proportion FROM _st_histogram($1, $2, FALSE, $3, 0, NULL, FALSE) $$ + LANGUAGE 'sql' IMMUTABLE STRICT; + ----------------------------------------------------------------------- -- MapAlgebra ----------------------------------------------------------------------- diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c index 7e2af187a..320702856 100644 --- a/raster/test/core/testapi.c +++ b/raster/test/core/testapi.c @@ -999,8 +999,21 @@ struct rt_bandstats_t { double *values; int sorted; }; +struct rt_histogram_t { + uint32_t count; + double proportion; + + double min; + double max; + + int inc_min; + int inc_max; +}; static void testBandStats() { rt_bandstats stats = NULL; + rt_histogram histogram = NULL; + double bin_width[] = {10000}; + int count = 0; rt_raster raster; rt_band band; @@ -1032,12 +1045,28 @@ static void testBandStats() { CHECK_EQUALS(stats->min, 1); CHECK_EQUALS(stats->max, 19998); + histogram = (rt_histogram) rt_band_get_histogram(stats, 0, NULL, 0, 0, &count); + CHECK(histogram); + rtdealloc(histogram); + + histogram = (rt_histogram) rt_band_get_histogram(stats, 0, NULL, 0, 1, &count); + CHECK(histogram); + rtdealloc(histogram); + + histogram = (rt_histogram) rt_band_get_histogram(stats, 0, bin_width, 1, 0, &count); + CHECK(histogram); + rtdealloc(histogram); + rtdealloc(stats->values); rtdealloc(stats); stats = (rt_bandstats) rt_band_get_summary_stats(band, 0, 0.1, 1); CHECK(stats); + histogram = (rt_histogram) rt_band_get_histogram(stats, 0, NULL, 0, 0, &count); + CHECK(histogram); + rtdealloc(histogram); + rtdealloc(stats->values); rtdealloc(stats); diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index 2712ff370..ac80ebef8 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -65,6 +65,7 @@ TEST_BANDPROPS = \ rt_mean.sql \ rt_stddev.sql \ rt_minmax.sql \ + rt_histogram.sql \ $(NULL) TEST_PIXEL = \ diff --git a/raster/test/regress/rt_histogram.sql b/raster/test/regress/rt_histogram.sql new file mode 100644 index 000000000..b6791e242 --- /dev/null +++ b/raster/test/regress/rt_histogram.sql @@ -0,0 +1,158 @@ +SELECT * FROM ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ) +); +SELECT ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ) +); +SELECT * FROM ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, TRUE, 0, ARRAY[]::double precision[], FALSE +); +SELECT * FROM ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, TRUE, 1, FALSE +); +SELECT * FROM ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, TRUE, 5 +); +SELECT * FROM ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, TRUE +); +SELECT * FROM ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1 +); +SELECT * FROM ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, 0, ARRAY[5]::double precision[], FALSE +); +SELECT * FROM ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, 3, FALSE +); +SELECT * FROM ST_Histogram( + ST_SetValue( + ST_SetValue( + ST_SetValue( + ST_AddBand( + ST_MakeEmptyRaster(10, 10, 10, 10, 2, 2, 0, 0,-1) + , 1, '64BF', 0, 0 + ) + , 1, 1, 1, -10 + ) + , 1, 5, 4, 0 + ) + , 1, 5, 5, 3.14159 + ), + 1, 5 +); diff --git a/raster/test/regress/rt_histogram_expected b/raster/test/regress/rt_histogram_expected new file mode 100644 index 000000000..1d068b835 --- /dev/null +++ b/raster/test/regress/rt_histogram_expected @@ -0,0 +1,39 @@ +-10|-3.429205|1|0.5 +-3.429205|3.14159|1|0.5 +(-10,-3.429205,1,0.5) +(-3.429205,3.14159,1,0.5) +-10|-8.357301|1|0.01 +-8.357301|-6.714602|0|0 +-6.714602|-5.071904|0|0 +-5.071904|-3.429205|0|0 +-3.429205|-1.786506|0|0 +-1.786506|-0.143807|0|0 +-0.143807|1.498891|98|0.98 +1.498891|3.14159|1|0.01 +-10|3.14159|100|-1 +-10|-7.371682|1|0.01 +-7.371682|-4.743364|0|0 +-4.743364|-2.115046|0|0 +-2.115046|0.513272|98|0.98 +0.513272|3.14159|1|0.01 +-10|-8.357301|1|0.01 +-8.357301|-6.714602|0|0 +-6.714602|-5.071904|0|0 +-5.071904|-3.429205|0|0 +-3.429205|-1.786506|0|0 +-1.786506|-0.143807|0|0 +-0.143807|1.498891|98|0.98 +1.498891|3.14159|1|0.01 +-10|-3.429205|1|0.5 +-3.429205|3.14159|1|0.5 +-10|-5|1|0.5 +-5|0|0|0 +0|5|1|0.5 +-10|-5.61947|1|0.5 +-5.61947|-1.23894|0|0 +-1.23894|3.14159|1|0.5 +-10|-7.371682|1|0.5 +-7.371682|-4.743364|0|0 +-4.743364|-2.115046|0|0 +-2.115046|0.513272|0|0 +0.513272|3.14159|1|0.5