Compute the default set of or requested quantiles for a coverage.
A One-Pass Space-Efficient Algorithm for Finding Quantiles (1995) by Rakesh Agrawal, Arun Swami in Proc. 7th Intl. Conf. Management of Data (COMAD-95)
In the future, it may be worth exploring algorithms that don't require the size of the coverage
1019 int init_quantiles = 0;
1050 assert(NULL !=
band);
1051 assert(cov_count > 1);
1052 assert(NULL != rtn_count);
1057 rterror(
"rt_band_get_summary_stats: Cannot get band data");
1062 exclude_nodata_value = 0;
1063 RASTER_DEBUGF(3,
"exclude_nodata_value = %d", exclude_nodata_value);
1066 if (NULL == *qlls) {
1068 if (NULL == quantiles) {
1070 if (quantiles_count < 2)
1071 quantiles_count = 5;
1073 quantiles =
rtalloc(
sizeof(
double) * quantiles_count);
1075 if (NULL == quantiles) {
1076 rterror(
"rt_band_get_quantiles_stream: Could not allocate memory for quantile input");
1081 for (i = 0; i <= quantiles_count; i++)
1082 quantiles[i] = ((
double) i) / quantiles_count;
1087 for (i = 0; i < quantiles_count; i++) {
1088 if (quantiles[i] < 0. || quantiles[i] > 1.) {
1089 rterror(
"rt_band_get_quantiles_stream: Quantile value not between 0 and 1");
1090 if (init_quantiles)
rtdealloc(quantiles);
1094 quicksort(quantiles, quantiles + quantiles_count - 1);
1097 *qlls_count = quantiles_count * 2;
1100 if (NULL == *qlls) {
1101 rterror(
"rt_band_get_quantiles_stream: Could not allocate memory for quantile output");
1102 if (init_quantiles)
rtdealloc(quantiles);
1106 j = (
uint32_t) floor(MAX_VALUES / 100.) + 1;
1107 for (i = 0; i < *qlls_count; i++) {
1108 qll = &((*qlls)[i]);
1109 qll->
quantile = quantiles[(i * quantiles_count) / *qlls_count];
1118 if (NULL == qll->
index) {
1119 rterror(
"rt_band_get_quantiles_stream: Could not allocate memory for quantile output");
1120 if (init_quantiles)
rtdealloc(quantiles);
1130 if (qll->
tau < 1) qll->
tau = 1;
1135 qll->
tau = cov_count - (*qlls)[i - 1].tau + 1;
1138 RASTER_DEBUGF(4,
"qll init: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
1143 if (init_quantiles)
rtdealloc(quantiles);
1148 (sample < 0 ||
FLT_EQ(sample, 0.0)) ||
1149 (sample > 1 ||
FLT_EQ(sample, 1.0))
1160 sample_size =
band->width *
band->height;
1161 sample_per =
band->height;
1169 sample_size = round((
band->width *
band->height) * sample);
1170 sample_per = round(sample_size /
band->width);
1171 sample_int = round(
band->height / sample_per);
1174 RASTER_DEBUGF(3,
"sampling %d of %d available pixels w/ %d per set"
1175 , sample_size, (
band->width *
band->height), sample_per);
1177 for (
x = 0, j = 0, k = 0;
x <
band->width;
x++) {
1183 RASTER_DEBUG(3,
"Skipping quantile calcuation as band is NODATA");
1187 for (i = 0, z = 0; i < sample_per; i++) {
1191 offset = (rand() % sample_int) + 1;
1193 diff = sample_int - offset;
1196 if (
y >=
band->height || z > sample_per)
break;
1201 if (status ==
ES_NONE && (!exclude_nodata_value || (exclude_nodata_value && !isnodata))) {
1204 for (a = 0; a < *qlls_count; a++) {
1205 qll = &((*qlls)[a]);
1208 RASTER_DEBUGF(5,
"qll before: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
1214 if (NULL != qll->
head) {
1215 if (value < qll->head->value)
1229 if (NULL != qll->
head) {
1273 else if (qll->
count < MAX_VALUES) {
1285 if (NULL == qle)
return NULL;
1291 if (NULL == qle->
prev)
1294 if (NULL == qle->
next)
1305 RASTER_DEBUGF(5,
"qle, prev, next, head, tail = %p, %p, %p, %p, %p", qle, qle->
prev, qle->
next, qll->
head, qll->
tail);
1308 else if (qll->
algeq) {
1311 if (value < qll->head->value) {
1341 if (NULL == qle)
return NULL;
1347 if (NULL == qle->
prev)
1350 if (NULL == qle->
next)
1363 while (NULL != qle) {
1411 if (NULL == qle)
return NULL;
1417 if (NULL == qle->
prev)
1420 if (NULL == qle->
next)
1433 while (NULL != qle) {
1501 RASTER_DEBUGF(5,
"qll after: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
1516 *rtn_count = *qlls_count / 2;
1518 if (NULL == rtn)
return NULL;
1521 for (i = 0, k = 0; i < *qlls_count; i++) {
1522 qll = &((*qlls)[i]);
1525 for (
x = 0;
x < k;
x++) {
1531 if (exists)
continue;
1533 RASTER_DEBUGF(5,
"qll: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
1541 if (qll->
head == NULL || qll->
tail == NULL)
1552 for (j = i + 1; j < *qlls_count; j++) {
1555 RASTER_DEBUGF(5,
"qlls[%d]: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
1556 j, (*qlls)[j].algeq, (*qlls)[j].quantile, (*qlls)[j].
count, (*qlls)[j].tau, (*qlls)[j].sum1, (*qlls)[j].sum2);
1557 RASTER_DEBUGF(5,
"qlls[%d]: (head, tail) = (%p, %p)", j, (*qlls)[j].head, (*qlls)[j].tail);
1566 if ((*qlls)[j].algeq) {
1567 rtn[k].
value = ((qle->
value * qle->
count) + ((*qlls)[j].head->value * (*qlls)[j].head->count)) / (qle->
count + (*qlls)[j].head->count);
1568 RASTER_DEBUGF(5,
"qlls[%d].head: (value, count) = (%f, %d)", j, (*qlls)[j].head->value, (*qlls)[j].head->count);
1571 rtn[k].
value = ((qle->
value * qle->
count) + ((*qlls)[j].tail->value * (*qlls)[j].tail->count)) / (qle->
count + (*qlls)[j].tail->count);
1572 RASTER_DEBUGF(5,
"qlls[%d].tail: (value, count) = (%f, %d)", j, (*qlls)[j].tail->value, (*qlls)[j].tail->count);
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
void * rtalloc(size_t size)
Wrappers used for managing memory.
#define RASTER_DEBUG(level, msg)
#define RASTER_DEBUGF(level, msg,...)
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
void rtdealloc(void *mem)
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
static void quantile_llist_index_update(struct quantile_llist *qll, struct quantile_llist_element *qle, uint32_t idx)
static void quantile_llist_index_delete(struct quantile_llist *qll, struct quantile_llist_element *qle)
static struct quantile_llist_element * quantile_llist_index_search(struct quantile_llist *qll, double value, uint32_t *index)
static struct quantile_llist_element * quantile_llist_insert(struct quantile_llist_element *element, double value, uint32_t *idx)
static void quicksort(double *left, double *right)
static int quantile_llist_delete(struct quantile_llist_element *element)
static struct quantile_llist_element * quantile_llist_search(struct quantile_llist_element *element, double needle)
static void quantile_llist_index_reset(struct quantile_llist *qll)
struct quantile_llist_element * prev
struct quantile_llist_element * next
struct quantile_llist_element * head
struct quantile_llist_element * tail
struct quantile_llist_index * index