From: Bborie Park <bkpark at ucdavis.edu>
Date: Thu, 26 May 2011 04:23:54 +0000 (+0000)
Subject: Adds ST_ValueCount to count the number of times a user-provided value or all values... 
X-Git-Tag: 2.0.0alpha1~1562
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=43f42cd6b8b419fbe3595e3073ec71c179d9ecfb;p=postgis

Adds ST_ValueCount to count the number of times a user-provided value or all values occurs in a raster's band.

Associated ticket is #953


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

diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c
index b1fc8114f..d3a9ec19d 100644
--- a/raster/rt_core/rt_api.c
+++ b/raster/rt_core/rt_api.c
@@ -130,8 +130,8 @@ rt_util_clamp_to_32F(double value) {
 }
 
 /* quicksort */
-#define   swap(x, y)    { double t; t = x; x = y; y = t; }
-#define   order(x, y)   if (x > y) swap(x, y)
+#define SWAP(x, y) { double t; t = x; x = y; y = t; }
+#define ORDER(x, y) if (x > y) SWAP(x, y)
 
 static double pivot(double *left, double *right) {
 	double l, m, r, *p;
@@ -141,9 +141,9 @@ static double pivot(double *left, double *right) {
 	r = *right;
 
 	/* order */
-	order(l, m);
-	order(l, r);
-	order(m, r);
+	ORDER(l, m);
+	ORDER(l, r);
+	ORDER(m, r);
 
 	/* pivot is higher of two values */
 	if (l < m) return m;
@@ -165,7 +165,7 @@ static double *partition(double *left, double *right, double pivot) {
 		while (*right >= pivot) --right;
 
 		if (left < right) {
-			swap(*left, *right);
+			SWAP(*left, *right);
 			++left;
 			--right;
 		}
@@ -1500,7 +1500,7 @@ struct rt_bandstats_t {
 /**
  * Compute summary statistics for a band
  *
- * @param band: the band to query for minimum and maximum pixel values
+ * @param band: the band to query for summary stats
  * @param hasnodata: if non-zero, ignore nodata values
  * @param sample: percentage of pixels to sample
  * @param inc_vals: flag to include values in return struct
@@ -1510,7 +1510,6 @@ struct rt_bandstats_t {
 rt_bandstats
 rt_band_get_summary_stats(rt_band band, int hasnodata, double sample,
 	int inc_vals) {
-	rt_pixtype pixtype = PT_END;
 	uint8_t *data = NULL;
 	uint32_t x = 0;
 	uint32_t y = 0;
@@ -1556,15 +1555,12 @@ rt_band_get_summary_stats(rt_band band, int hasnodata, double sample,
 	}
 
 	data = rt_band_get_data(band);
-	pixtype = band->pixtype;
 
 	hasnodata_flag = rt_band_get_hasnodata_flag(band);
-	if (hasnodata_flag != FALSE) {
+	if (hasnodata_flag != FALSE)
 		nodata = rt_band_get_nodata(band);
-	}
-	else {
+	else
 		hasnodata = 0;
-	}
 
 	RASTER_DEBUGF(3, "nodata = %f", nodata);
 	RASTER_DEBUGF(3, "hasnodata_flag = %d", hasnodata_flag);
@@ -2143,6 +2139,273 @@ rt_band_get_quantiles(rt_bandstats stats,
 	return rtn;
 }
 
+/* symmetrical rounding */
+#define ROUND(x, y) (((x > 0.0) ? floor((x * pow(10, y) + 0.5)) : ceil((x * pow(10, y) - 0.5))) / pow(10, y));
+
+struct rt_valuecount_t {
+	double value;
+	uint32_t count;
+	double proportion;
+};
+
+/**
+ * Count the number of times provided value(s) occur in
+ * the band
+ *
+ * @param band: the band to query for minimum and maximum pixel values
+ * @param hasnodata: if non-zero, ignore nodata values
+ * @param search_values: array of values to count
+ * @param search_values_count: the number of search values
+ * @param roundto: the decimal place to round the values to
+ * @param rtn_count: the number of value counts being returned
+ *
+ * @return the default set of or requested quantiles for a band
+ */
+rt_valuecount
+rt_band_get_value_count(rt_band band, int hasnodata,
+	double *search_values, uint32_t search_values_count, double roundto,
+	int *rtn_count) {
+	rt_valuecount vcnts = NULL;
+	rt_pixtype pixtype = PT_END;
+	uint8_t *data = NULL;
+	int hasnodata_flag = FALSE;
+	double nodata = 0;
+
+	int scale = 0;
+	int doround = 0;
+	double tmpd = 0;
+	int i = 0;
+
+	uint32_t x = 0;
+	uint32_t y = 0;
+	int rtn;
+	double pxlval;
+	double rpxlval;
+	uint32_t total = 0;
+	int vcnts_count = 0;
+	int new_valuecount = 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 != band);
+
+	if (band->offline) {
+		rterror("rt_band_get_count_of_values not implemented yet for OFFDB bands");
+		return NULL;
+	}
+
+	data = rt_band_get_data(band);
+	pixtype = band->pixtype;
+
+	hasnodata_flag = rt_band_get_hasnodata_flag(band);
+	if (hasnodata_flag != FALSE)
+		nodata = rt_band_get_nodata(band);
+	else
+		hasnodata = 0;
+
+	RASTER_DEBUGF(3, "nodata = %f", nodata);
+	RASTER_DEBUGF(3, "hasnodata_flag = %d", hasnodata_flag);
+	RASTER_DEBUGF(3, "user hasnodata = %d", hasnodata);
+
+	/* process roundto */
+	if (roundto < 0 || (fabs(roundto - 0.0) < FLT_EPSILON)) {
+		roundto = 0;
+		scale = 0;
+	}
+	/* tenths, hundredths, thousandths, etc */
+	else if (roundto < 1) {
+    switch (pixtype) {
+			/* integer band types don't have digits after the decimal place */
+			case PT_1BB:
+			case PT_2BUI:
+			case PT_4BUI:
+			case PT_8BSI:
+			case PT_8BUI:
+			case PT_16BSI:
+			case PT_16BUI:
+			case PT_32BSI:
+			case PT_32BUI:
+				roundto = 0;
+				break;
+			/* floating points, check the rounding */
+			case PT_32BF:
+			case PT_64BF:
+				for (scale = 0; scale <= 20; scale++) {
+					tmpd = roundto * pow(10, scale);
+					if (fabs((tmpd - ((int) tmpd)) - 0.0) < FLT_EPSILON) break;
+				}
+				break;
+			case PT_END:
+				break;
+		}
+	}
+	/* ones, tens, hundreds, etc */
+	else {
+		for (scale = 0; scale >= -20; scale--) {
+			tmpd = roundto * pow(10, scale);
+			if (tmpd < 1 || fabs(tmpd - 1.0) < FLT_EPSILON) {
+				if (scale == 0) doround = 1;
+				break;
+			}
+		}
+	}
+
+	if (scale != 0 || doround)
+		doround = 1;
+	else
+		doround = 0;
+	RASTER_DEBUGF(3, "scale = %d", scale);
+	RASTER_DEBUGF(3, "doround = %d", doround);
+
+	/* process search_values */
+	if (search_values_count > 0 && NULL != search_values) {
+		vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t) * search_values_count);
+		if (NULL == vcnts) {
+			rterror("rt_band_get_count_of_values: Unable to allocate memory for value counts");
+			*rtn_count = 0;
+			return NULL;
+		}
+
+		for (i = 0; i < search_values_count; i++) {
+			vcnts[i].count = 0;
+			vcnts[i].proportion = 0;
+			if (!doround)
+				vcnts[i].value = search_values[i];
+			else
+				vcnts[i].value = ROUND(search_values[i], scale);
+		}
+		vcnts_count = i;
+	}
+	else
+		search_values_count = 0;
+	RASTER_DEBUGF(3, "search_values_count = %d", search_values_count);
+
+	/* entire band is nodata */
+	if (rt_band_get_isnodata_flag(band) != FALSE) {
+		if (hasnodata) {
+			rtwarn("All pixels of band have the NODATA value");
+			return NULL;
+		}
+		else {
+			if (search_values_count > 0) {
+				/* check for nodata match */
+				for (i = 0; i < search_values_count; i++) {
+					if (!doround)
+						tmpd = nodata;
+					else
+						tmpd = ROUND(nodata, scale);
+
+					if (fabs(tmpd - vcnts[i].value) > FLT_EPSILON)
+						continue;
+
+					vcnts[i].count = band->width * band->height;
+					vcnts->proportion = 1.0;
+				}
+
+				*rtn_count = vcnts_count;
+			}
+			/* no defined search values */
+			else {
+				vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t));
+				if (NULL == vcnts) {
+					rterror("rt_band_get_count_of_values: Unable to allocate memory for value counts");
+					*rtn_count = 0;
+					return NULL;
+				}
+
+				vcnts->value = nodata;
+				vcnts->count = band->width * band->height;
+				vcnts->proportion = 1.0;
+
+				*rtn_count = 1;
+			}
+
+			return vcnts;
+		}
+	}
+
+	for (x = 0; x < band->width; x++) {
+		for (y = 0; y < band->height; y++) {
+			rtn = rt_band_get_pixel(band, x, y, &pxlval);
+
+			/* error getting value, continue */
+			if (rtn == -1) continue;
+
+			if (
+				!hasnodata || (
+					hasnodata &&
+					(hasnodata_flag != FALSE) &&
+					(fabs(pxlval - nodata) > FLT_EPSILON)
+				)
+			) {
+				total++;
+				if (doround) {
+					rpxlval = ROUND(pxlval, scale);
+				}
+				else
+					rpxlval = pxlval;
+				RASTER_DEBUGF(5, "(pxlval, rpxlval) => (%0.6f, %0.6f)", pxlval, rpxlval);
+
+				new_valuecount = 1;
+				/* search for match in existing valuecounts */
+				for (i = 0; i < vcnts_count; i++) {
+					/* match found */
+					if (fabs(vcnts[i].value - rpxlval) < FLT_EPSILON) {
+						vcnts[i].count++;
+						new_valuecount = 0;
+						RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
+						break;
+					}
+				}
+
+				/*
+					don't add new valuecount either because
+						- no need for new one
+						- user-defined search values
+				*/
+				if (!new_valuecount || search_values_count > 0) continue;
+
+				/* add new valuecount */
+				vcnts = rtrealloc(vcnts, sizeof(struct rt_valuecount_t) * (vcnts_count + 1));
+				if (NULL == vcnts) {
+					rterror("rt_band_get_count_of_values: Unable to allocate memory for value counts");
+					*rtn_count = 0;
+					return NULL;
+				}
+
+				vcnts[vcnts_count].value = rpxlval;
+				vcnts[vcnts_count].count = 1;
+				vcnts[vcnts_count].proportion = 0;
+				RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[vcnts_count].value, vcnts[vcnts_count].count);
+				vcnts_count++;
+			}
+		}
+	}
+
+#if POSTGIS_DEBUG_LEVEL > 0
+	stop = clock();
+	elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
+	RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
+#endif
+
+	for (i = 0; i < vcnts_count; i++) {
+		vcnts[i].proportion = (double) vcnts[i].count / total;
+		RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
+	}
+
+	RASTER_DEBUG(3, "done");
+	*rtn_count = vcnts_count;
+	return vcnts;
+}
+
 struct rt_reclassexpr_t {
 	struct rt_reclassrange {
 		double min;
@@ -2171,8 +2434,8 @@ rt_band_reclass(rt_band srcband, rt_pixtype pixtype,
 	uint32_t hasnodata, double nodataval, rt_reclassexpr *exprset,
 	int exprcount) {
 	rt_band band = NULL;
-	int width = 0;
-	int height = 0;
+	uint32_t width = 0;
+	uint32_t height = 0;
 	int numval = 0;
 	int memsize = 0;
 	void *mem = NULL;
@@ -2180,8 +2443,8 @@ rt_band_reclass(rt_band srcband, rt_pixtype pixtype,
 	double src_nodataval = 0.0;
 
 	int rtn;
-	int x;
-	int y;
+	uint32_t x;
+	uint32_t y;
 	int i;
 	double or = 0;
 	double ov = 0;
diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h
index 86145a157..b07926281 100644
--- a/raster/rt_core/rt_api.h
+++ b/raster/rt_core/rt_api.h
@@ -450,7 +450,7 @@ int rt_band_check_is_nodata(rt_band band);
 /**
  * Compute summary statistics for a band
  *
- * @param band: the band to query for minimum and maximum pixel values
+ * @param band: the band to query for summary stats 
  * @param hasnodata: if non-zero, ignore nodata values
  * @param sample: percentage of pixels to sample
  * @param inc_vals: flag to include values in return struct
@@ -493,6 +493,24 @@ typedef struct rt_quantile_t* rt_quantile;
 rt_quantile rt_band_get_quantiles(rt_bandstats stats,
 	double *quantiles, int quantiles_count, int *rtn_count);
 
+/**
+ * Count the number of times provided value(s) occur in
+ * the band
+ *
+ * @param band: the band to query for minimum and maximum pixel values
+ * @param hasnodata: if non-zero, ignore nodata values
+ * @param search_values: array of values to count
+ * @param search_values_count: the number of search values
+ * @param roundto: the decimal place to round the values to
+ * @param rtn_count: the number of value counts being returned
+ *
+ * @return the default set of or requested quantiles for a band
+ */
+typedef struct rt_valuecount_t* rt_valuecount;
+rt_valuecount rt_band_get_value_count(rt_band band, int hasnodata,
+	double *search_values, uint32_t search_values_count,
+	double roundto, int *rtn_count);
+
 /**
  * Returns new band with values reclassified
  * 
diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c
index 9c0adf148..aa1c4eb6a 100644
--- a/raster/rt_pg/rt_pg.c
+++ b/raster/rt_pg/rt_pg.c
@@ -209,6 +209,9 @@ Datum RASTER_histogram(PG_FUNCTION_ARGS);
 /* get quantiles */
 Datum RASTER_quantile(PG_FUNCTION_ARGS);
 
+/* get counts of values */
+Datum RASTER_valueCount(PG_FUNCTION_ARGS);
+
 /* reclassify specified bands of a raster */
 Datum RASTER_reclass(PG_FUNCTION_ARGS);
 
@@ -3656,6 +3659,250 @@ Datum RASTER_quantile(PG_FUNCTION_ARGS)
 	}
 }
 
+struct rt_valuecount_t {
+	double value;
+	uint32_t count;
+	double proportion;
+};
+
+/* get counts of values */
+PG_FUNCTION_INFO_V1(RASTER_valueCount);
+Datum RASTER_valueCount(PG_FUNCTION_ARGS) {
+	FuncCallContext *funcctx;
+	TupleDesc tupdesc;
+	AttInMetadata *attinmeta;
+
+	int count;
+	rt_valuecount vcnts;
+	rt_valuecount vcnts2;
+	int call_cntr;
+	int max_calls;
+
+	/* first call of function */
+	if (SRF_IS_FIRSTCALL()) {
+		MemoryContext oldcontext;
+
+		rt_pgraster *pgraster = NULL;
+		rt_raster raster = NULL;
+		rt_band band = NULL;
+		int32_t bandindex = 0;
+		int num_bands = 0;
+		bool hasnodata = TRUE;
+		double *search_values = NULL;
+		int search_values_count = 0;
+		double roundto = 0;
+
+		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);
+
+		/* pgraster is null, return nothing */
+		if (PG_ARGISNULL(0)) SRF_RETURN_DONE(funcctx);
+		pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+
+		raster = rt_raster_deserialize(pgraster);
+		if (!raster) {
+			elog(ERROR, "RASTER_valueCount: 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);
+
+		/* search values */
+		if (!PG_ARGISNULL(3)) {
+			array = PG_GETARG_ARRAYTYPE_P(3);
+			etype = ARR_ELEMTYPE(array);
+			get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
+
+			switch (etype) {
+				case FLOAT4OID:
+				case FLOAT8OID:
+					break;
+				default:
+					elog(ERROR, "RASTER_valueCount: Invalid data type for values");
+					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);
+
+			search_values = palloc(sizeof(double) * n);
+			for (i = 0, j = 0; i < n; i++) {
+				if (nulls[i]) continue;
+
+				switch (etype) {
+					case FLOAT4OID:
+						search_values[j] = (double) DatumGetFloat4(e[i]);
+						break;
+					case FLOAT8OID:
+						search_values[j] = (double) DatumGetFloat8(e[i]);
+						break;
+				}
+
+				POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
+				j++;
+			}
+			search_values_count = j;
+
+			if (j < 1) {
+				pfree(search_values);
+				search_values = NULL;
+			}
+		}
+
+		/* roundto */
+		if (!PG_ARGISNULL(4)) {
+			roundto = PG_GETARG_FLOAT8(4);
+			if (roundto < 0.) roundto = 0;
+		}
+
+		/* 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 counts of values */
+		vcnts = rt_band_get_value_count(band, (int) hasnodata, search_values, search_values_count, roundto, &count);
+		rt_band_destroy(band);
+		rt_raster_destroy(raster);
+		if (NULL == vcnts || !count) {
+			elog(NOTICE, "Could not count the values of raster band of index %d", bandindex);
+			SRF_RETURN_DONE(funcctx);
+		}
+
+		POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
+
+		/* Store needed information */
+		funcctx->user_fctx = vcnts;
+
+		/* 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;
+	vcnts2 = 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(3 * sizeof(char *));
+
+		values[0] = (char *) palloc(sizeof(char) * (MAX_DBL_CHARLEN + 1));
+		values[1] = (char *) palloc(sizeof(char) * (MAX_INT_CHARLEN + 1));
+		values[2] = (char *) palloc(sizeof(char) * (MAX_DBL_CHARLEN + 1));
+
+		snprintf(
+			values[0],
+			sizeof(char) * (MAX_DBL_CHARLEN + 1),
+			"%f",
+			vcnts2[call_cntr].value
+		);
+		snprintf(
+			values[1],
+			sizeof(char) * (MAX_INT_CHARLEN + 1),
+			"%d",
+			vcnts2[call_cntr].count
+		);
+		snprintf(
+			values[2],
+			sizeof(char) * (MAX_DBL_CHARLEN + 1),
+			"%f",
+			vcnts2[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[2]);
+		pfree(values[1]);
+		pfree(values[0]);
+		pfree(values);
+
+		SRF_RETURN_NEXT(funcctx, result);
+	}
+	/* do when there is no more left */
+	else {
+		pfree(vcnts2);
+		SRF_RETURN_DONE(funcctx);
+	}
+}
+
 struct rt_reclassexpr_t {
 	struct rt_reclassrange {
 		double min;
diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c
index 9d0ca23e6..19bf95225 100644
--- a/raster/rt_pg/rtpostgis.sql.in.c
+++ b/raster/rt_pg/rtpostgis.sql.in.c
@@ -1163,6 +1163,188 @@ CREATE OR REPLACE FUNCTION st_approxquantile(rast raster, hasnodata boolean, qua
 	AS $$ SELECT quantile, value FROM _st_quantile($1, 1, $2, 1, ARRAY[$3]::double precision[]) $$
 	LANGUAGE 'sql' IMMUTABLE STRICT;
 
+-----------------------------------------------------------------------
+-- ST_ValueCount
+-----------------------------------------------------------------------
+CREATE TYPE valuecount AS (
+	value double precision,
+	count integer,
+	proportion double precision
+);
+
+-- None of the functions can be strict as "searchvalues" and "roundto" can be NULL
+
+CREATE OR REPLACE FUNCTION _st_valuecount(rast raster, nband integer, hasnodata boolean, searchvalues double precision[], roundto double precision)
+	RETURNS SETOF valuecount
+	AS 'MODULE_PATHNAME', 'RASTER_valueCount'
+	LANGUAGE 'C' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, hasnodata boolean, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count integer)
+	RETURNS SETOF record
+	AS $$ SELECT value, count FROM _st_valuecount($1, $2, $3, $4, $5) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count integer)
+	RETURNS SETOF record
+	AS $$ SELECT value, count FROM _st_valuecount($1, $2, TRUE, $3, $4) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, searchvalues double precision[], OUT value double precision, OUT count integer)
+	RETURNS SETOF record
+	AS $$ SELECT value, count FROM _st_valuecount($1, $2, TRUE, $3, 0) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count integer)
+	RETURNS SETOF record
+	AS $$ SELECT value, count FROM _st_valuecount($1, 1, TRUE, $2, $3) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, searchvalues double precision[], OUT value double precision, OUT count integer)
+	RETURNS SETOF record
+	AS $$ SELECT value, count FROM _st_valuecount($1, 1, TRUE, $2, 0) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, hasnodata boolean, searchvalue double precision, roundto double precision, OUT value double precision, OUT count integer)
+	RETURNS record
+	AS $$ SELECT value, count FROM _st_valuecount($1, $2, $3, ARRAY[$4]::double precision[], $5) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, searchvalue double precision, roundto double precision, OUT value double precision, OUT count integer)
+	RETURNS record
+	AS $$ SELECT value, count FROM _st_valuecount($1, $2, TRUE, ARRAY[$3]::double precision[], $4) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, nband integer, searchvalue double precision, OUT value double precision, OUT count integer)
+	RETURNS record
+	AS $$ SELECT value, count FROM _st_valuecount($1, $2, TRUE, ARRAY[$3]::double precision[], 0) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, searchvalue double precision, roundto double precision, OUT value double precision, OUT count integer)
+	RETURNS record
+	AS $$ SELECT value, count FROM _st_valuecount($1, 1, TRUE, ARRAY[$2]::double precision[], $3) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rast raster, searchvalue double precision, OUT value double precision, OUT count integer)
+	RETURNS record
+	AS $$ SELECT value, count FROM _st_valuecount($1, 1, TRUE, ARRAY[$2]::double precision[], 0) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION _st_valuecount(rastertable text, rastercolumn text, nband integer, hasnodata boolean, searchvalues double precision[], roundto double precision)
+	RETURNS SETOF valuecount
+	AS $$
+	DECLARE
+		curs refcursor;
+
+		ctable text;
+		ccolumn text;
+		rast raster;
+
+		vcnts valuecount;
+	BEGIN		
+		-- nband
+		IF nband < 1 THEN
+			RAISE WARNING 'Invalid band index (must use 1-based). Returning NULL';
+			RETURN;
+		END IF;
+
+		-- rastertable and rastercolumn
+		IF rastertable IS NULL THEN
+			RAISE WARNING 'rastertable cannot be NULL. Returning NULL';
+			RETURN;
+		ELSEIF rastercolumn IS NULL THEN
+			RAISE WARNING 'rastercolumn cannot be NULL. Returning NULL';
+			RETURN;
+		END IF;
+
+		-- clean rastertable and rastercolumn
+		ctable := quote_ident(rastertable);
+		ccolumn := quote_ident(rastercolumn);
+
+		BEGIN
+			OPEN curs FOR EXECUTE 'SELECT '
+					|| ccolumn
+					|| ' FROM '
+					|| ctable
+					|| ' WHERE '
+					|| ccolumn
+					|| ' IS NOT NULL';
+		EXCEPTION
+			WHEN OTHERS THEN
+				RAISE WARNING 'Invalid table or column name. Returning NULL';
+				RETURN;
+		END;
+
+		LOOP
+			FETCH curs INTO rast;
+			EXIT WHEN NOT FOUND;
+
+			FOR vcnts IN SELECT * FROM _st_valuecount(rast, nband, hasnodata, searchvalues, roundto) LOOP
+				IF vcnts IS NULL THEN
+					CONTINUE;
+				END IF;
+				--RAISE NOTICE 'vcnts = %', vcnts;
+				
+				RETURN NEXT vcnts;
+			END LOOP;
+
+		END LOOP;
+
+		CLOSE curs;
+
+		RETURN;
+	END;
+	$$ LANGUAGE 'plpgsql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, hasnodata boolean, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count bigint)
+	RETURNS SETOF record
+	AS $$ SELECT value, sum(count) AS count FROM _st_valuecount($1, $2, $3, $4, $5, $6) GROUP BY 1 ORDER BY 1 $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count bigint)
+	RETURNS SETOF record
+	AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, TRUE, $4, $5) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, searchvalues double precision[], OUT value double precision, OUT count bigint)
+	RETURNS SETOF record
+	AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, TRUE, $4, 0) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, searchvalues double precision[], roundto double precision, OUT value double precision, OUT count bigint)
+	RETURNS SETOF record
+	AS $$ SELECT value, count FROM st_valuecount($1, $2, 1, TRUE, $3, $4) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, searchvalues double precision[], OUT value double precision, OUT count bigint)
+	RETURNS SETOF record
+	AS $$ SELECT value, count FROM st_valuecount($1, $2, 1, TRUE, $3, 0) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, hasnodata boolean, searchvalue double precision, roundto double precision, OUT value double precision, OUT count bigint)
+	RETURNS record
+	AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, $4, ARRAY[$5]::double precision[], $6) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, searchvalue double precision, roundto double precision, OUT value double precision, OUT count bigint)
+	RETURNS record
+	AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, TRUE, ARRAY[$4]::double precision[], $5) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, nband integer, searchvalue double precision, OUT value double precision, OUT count bigint)
+	RETURNS record
+	AS $$ SELECT value, count FROM st_valuecount($1, $2, $3, TRUE, ARRAY[$4]::double precision[], 0) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, searchvalue double precision, roundto double precision, OUT value double precision, OUT count bigint)
+	RETURNS record
+	AS $$ SELECT value, count FROM st_valuecount($1, $2, 1, TRUE, ARRAY[$3]::double precision[], $4) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_valuecount(rastertable text, rastercolumn text, searchvalue double precision, OUT value double precision, OUT count bigint)
+	RETURNS record
+	AS $$ SELECT value, count FROM st_valuecount($1, $2, 1, TRUE, ARRAY[$3]::double precision[], 0) $$
+	LANGUAGE 'sql' IMMUTABLE;
+
 -----------------------------------------------------------------------
 -- ST_Reclass
 -----------------------------------------------------------------------
diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c
index 0d9319e6d..8cff63f95 100644
--- a/raster/test/core/testapi.c
+++ b/raster/test/core/testapi.c
@@ -1124,8 +1124,7 @@ static void testBandStats() {
 	rtdealloc(stats->values);
 	rtdealloc(stats);
 
-	rt_band_destroy(band);
-	rt_raster_destroy(raster);
+	deepRelease(raster);
 }
 
 static void testRasterReplaceBand() {
@@ -1257,9 +1256,8 @@ static void testBandReclass() {
 	for (i = cnt - 1; i >= 0; i--)
 		free(exprset[i]);
 	free(exprset);
-	rt_band_destroy(band);
 	rt_band_destroy(newband);
-	rt_raster_destroy(raster);
+	deepRelease(raster);
 }
 
 struct rt_gdaldriver_t {
@@ -1307,6 +1305,69 @@ static void testRasterToGDAL(rt_raster raster) {
 	rt_raster_destroy(rast);
 }
 
+struct rt_valuecount_t {
+	double value;
+	uint32_t count;
+	double proportion;
+};
+static void testValueCount() {
+	rt_valuecount vcnts = NULL;
+
+	rt_raster raster;
+	rt_band band;
+	uint32_t x;
+	uint32_t xmax = 100;
+	uint32_t y;
+	uint32_t ymax = 100;
+	int rtn = 0;
+
+	double count[] = {3, 4, 5};
+
+	raster = rt_raster_new(xmax, ymax);
+	assert(raster); /* or we're out of virtual memory */
+	band = addBand(raster, PT_32BF, 0, 0);
+	CHECK(band);
+	rt_band_set_nodata(band, 0);
+
+	for (x = 0; x < xmax; x++) {
+		for (y = 0; y < ymax; y++) {
+			rtn = rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1));
+			CHECK((rtn != -1));
+		}
+	}
+	vcnts = rt_band_get_value_count(band, 1, NULL, 0, 0, &rtn);
+	CHECK(vcnts);
+	CHECK((rtn > 0));
+	rtdealloc(vcnts);
+
+	vcnts = rt_band_get_value_count(band, 1, NULL, 0, 0.01, &rtn);
+	CHECK(vcnts);
+	CHECK((rtn > 0));
+	rtdealloc(vcnts);
+
+	vcnts = rt_band_get_value_count(band, 1, NULL, 0, 0.1, &rtn);
+	CHECK(vcnts);
+	CHECK((rtn > 0));
+	rtdealloc(vcnts);
+
+	vcnts = rt_band_get_value_count(band, 1, NULL, 0, 1, &rtn);
+	CHECK(vcnts);
+	CHECK((rtn > 0));
+	rtdealloc(vcnts);
+
+	vcnts = rt_band_get_value_count(band, 1, NULL, 0, 10, &rtn);
+	CHECK(vcnts);
+	CHECK((rtn > 0));
+	rtdealloc(vcnts);
+
+	vcnts = rt_band_get_value_count(band, 1, count, 3, 1, &rtn);
+	CHECK(vcnts);
+	CHECK((rtn > 0));
+	rtdealloc(vcnts);
+
+	deepRelease(raster);
+}
+
 int
 main()
 {
@@ -1637,6 +1698,10 @@ main()
 		testGDALDrivers();
 		printf("Successfully tested rt_raster_gdal_drivers\n");
 
+		printf("Testing rt_band_get_value_count\n");
+		testValueCount();
+		printf("Successfully tested rt_band_get_value_count\n");
+
     deepRelease(raster);
 
     return EXIT_SUCCESS;
diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in
index 5829db928..ce9a7322f 100644
--- a/raster/test/regress/Makefile.in
+++ b/raster/test/regress/Makefile.in
@@ -75,6 +75,7 @@ TEST_BANDPROPS = \
 	rt_minmax.sql \
 	rt_histogram.sql \
 	rt_quantile.sql \
+	rt_valuecount.sql \
 	$(NULL)
 
 TEST_PIXEL = \
diff --git a/raster/test/regress/rt_valuecount.sql b/raster/test/regress/rt_valuecount.sql
new file mode 100644
index 000000000..85018f12e
--- /dev/null
+++ b/raster/test/regress/rt_valuecount.sql
@@ -0,0 +1,246 @@
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, ARRAY[]::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, NULL::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, FALSE, NULL::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, NULL::double precision[], 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, ARRAY[3.1], 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, ARRAY[-10]);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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
+	)
+, ARRAY[-10], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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
+	)
+, ARRAY[-10, 3]);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, 3.14, 0.19);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, FALSE, 3.14, 0.01);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, -10, 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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, -10);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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
+	)
+, -10., 10);
+SELECT round(value::numeric, 3), count FROM st_valuecount(
+	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
+	)
+, 3.14159);
+BEGIN;
+CREATE TEMP TABLE test
+	ON COMMIT DROP AS
+	SELECT
+		rast.rast
+	FROM (
+		SELECT 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
+		) AS rast
+	) AS rast
+	FULL JOIN (
+		SELECT generate_series(1, 10) AS id
+	) AS id
+		ON 1 = 1;
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, FALSE, NULL::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, ARRAY[-10]::double precision[], 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, ARRAY[1]::double precision[]);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', NULL::double precision[], 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', ARRAY[-1, 3.1]::double precision[], 0.1);
+
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, TRUE, NULL::double precision, 0);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, 3.14, 1);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 1, -1);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', 3.1, 0.1);
+SELECT round(value::numeric, 3), count FROM st_valuecount('test', 'rast', -9);
+ROLLBACK;
diff --git a/raster/test/regress/rt_valuecount_expected b/raster/test/regress/rt_valuecount_expected
new file mode 100644
index 000000000..48751d8ab
--- /dev/null
+++ b/raster/test/regress/rt_valuecount_expected
@@ -0,0 +1,36 @@
+-10.000|1
+3.142|1
+-10.000|1
+3.142|1
+-10.000|1
+0.000|98
+3.142|1
+-10.000|1
+3.100|1
+3.100|1
+-10.000|1
+-10.000|1
+-10.000|1
+3.000|0
+3.140|1
+3.140|1
+-10.000|1
+-10.000|1
+-10.000|1
+3.142|1
+BEGIN
+-10.000|10
+0.000|980
+3.142|10
+-10.000|10
+1.000|0
+-10.000|10
+3.100|10
+-1.000|0
+3.100|10
+-10.000|10
+3.000|10
+-1.000|0
+3.100|10
+-9.000|0
+COMMIT