]> granicus.if.org Git - postgis/commitdiff
Added ST_Histogram functions.
authorBborie Park <bkpark at ucdavis.edu>
Mon, 16 May 2011 19:36:55 +0000 (19:36 +0000)
committerBborie Park <bkpark at ucdavis.edu>
Mon, 16 May 2011 19:36:55 +0000 (19:36 +0000)
- 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

raster/rt_core/rt_api.c
raster/rt_core/rt_api.h
raster/rt_pg/rt_pg.c
raster/rt_pg/rtpostgis.sql.in.c
raster/test/core/testapi.c
raster/test/regress/Makefile.in
raster/test/regress/rt_histogram.sql [new file with mode: 0644]
raster/test/regress/rt_histogram_expected [new file with mode: 0644]

index cb01a5430e8e9615d8da6c14fd30aace55f71525..ac126ab3e92bf020041466e160d312ee8d54ce21 100644 (file)
@@ -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 {
index 968294bca737ccddb3b4bc74556ae2a578775f73..b975ba373fa0793c7cd0618793dccc0d9d50df06 100644 (file)
@@ -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 --------------------------------------------------------*/
 
index 4b169c8ed018f431c7490771f62301529246b373..3e2adb8dc44234bcb995330ff81613b08634e351 100644 (file)
@@ -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                       */
 /* ---------------------------------------------------------------- */
index 2191ee85d76c37a274805d55c6df3c946ea50fe3..e72cee4566511a992051cc0b26a0fbe27c7b883d 100644 (file)
@@ -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
 -----------------------------------------------------------------------
index 7e2af187a614c06c17f12a05660cd0fc491d6469..320702856a6d25c09e520b1be69f9d8adc6f4ba9 100644 (file)
@@ -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);
 
index 2712ff370c354f11f3d51ee5b371b71c40b53f4b..ac80ebef85553cf50f544b62cc722bd019cb9fc1 100644 (file)
@@ -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 (file)
index 0000000..b6791e2
--- /dev/null
@@ -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 (file)
index 0000000..1d068b8
--- /dev/null
@@ -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