PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ rt_band_get_quantiles_stream()

rt_quantile rt_band_get_quantiles_stream ( rt_band  band,
int  exclude_nodata_value,
double  sample,
uint64_t  cov_count,
struct quantile_llist **  qlls,
uint32_t *  qlls_count,
double *  quantiles,
int  quantiles_count,
uint32_t *  rtn_count 
)

Compute the default set of or requested quantiles for a coverage.

This function is based upon the algorithm described in:

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)

http://www.almaden.ibm.com/cs/projects/iis/hdb/Publications/papers/comad95.pdf

In the future, it may be worth exploring algorithms that don't require the size of the coverage

Parameters
band: the band to include in the quantile search
exclude_nodata_value: if non-zero, ignore nodata values
sample: percentage of pixels to sample
cov_count: number of values in coverage
qlls: set of quantile_llist structures
qlls_count: the number of quantile_llist structures
quantiles: the quantiles to be computed if bot qlls and quantiles provided, qlls is used
quantiles_count: the number of quantiles to be computed
rtn_count: the number of quantiles being returned
Returns
the default set of or requested quantiles for a band or NULL

Definition at line 1010 of file rt_statistics.c.

References quantile_llist::algeq, genraster::count, quantile_llist::count, quantile_llist_element::count, ovdump::data, ES_NONE, FLT_EQ, rt_quantile_t::has_value, quantile_llist::head, rt_band_t::height, quantile_llist::index, quantile_llist::index_max, quantile_llist_element::next, quantile_llist_element::prev, rt_quantile_t::quantile, quantile_llist::quantile, quantile_llist_delete(), quantile_llist_index_delete(), quantile_llist_index_reset(), quantile_llist_index_search(), quantile_llist_index_update(), quantile_llist_insert(), quantile_llist_search(), quicksort(), RASTER_DEBUG, RASTER_DEBUGF, ROUND, rt_band_get_data(), rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_band_get_pixel(), rtalloc(), rtdealloc(), rterror(), quantile_llist::sum1, quantile_llist::sum2, quantile_llist::tail, quantile_llist::tau, genraster::value, rt_quantile_t::value, quantile_llist_element::value, rt_band_t::width, pixval::x, and pixval::y.

Referenced by RASTER_quantileCoverage(), and test_band_stats().

1017  {
1018  rt_quantile rtn = NULL;
1019  int init_quantiles = 0;
1020 
1021  struct quantile_llist *qll = NULL;
1022  struct quantile_llist_element *qle = NULL;
1023  struct quantile_llist_element *qls = NULL;
1024  const uint32_t MAX_VALUES = 750;
1025 
1026  uint8_t *data = NULL;
1027  double value;
1028  int isnodata = 0;
1029 
1030  uint32_t a = 0;
1031  uint32_t i = 0;
1032  uint32_t j = 0;
1033  uint32_t k = 0;
1034  uint32_t x = 0;
1035  uint32_t y = 0;
1036  uint32_t z = 0;
1037  uint32_t idx = 0;
1038  uint32_t offset = 0;
1039  uint32_t diff = 0;
1040  uint8_t exists = 0;
1041 
1042  uint32_t do_sample = 0;
1043  uint32_t sample_size = 0;
1044  uint32_t sample_per = 0;
1045  uint32_t sample_int = 0;
1046  int status;
1047 
1048  RASTER_DEBUG(3, "starting");
1049 
1050  assert(NULL != band);
1051  assert(cov_count > 1);
1052  assert(NULL != rtn_count);
1053  RASTER_DEBUGF(3, "cov_count = %d", cov_count);
1054 
1055  data = rt_band_get_data(band);
1056  if (data == NULL) {
1057  rterror("rt_band_get_summary_stats: Cannot get band data");
1058  return NULL;
1059  }
1060 
1061  if (!rt_band_get_hasnodata_flag(band))
1062  exclude_nodata_value = 0;
1063  RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
1064 
1065  /* quantile_llist not provided */
1066  if (NULL == *qlls) {
1067  /* quantiles not provided */
1068  if (NULL == quantiles) {
1069  /* quantile count not specified, default to quartiles */
1070  if (quantiles_count < 2)
1071  quantiles_count = 5;
1072 
1073  quantiles = rtalloc(sizeof(double) * quantiles_count);
1074  init_quantiles = 1;
1075  if (NULL == quantiles) {
1076  rterror("rt_band_get_quantiles_stream: Could not allocate memory for quantile input");
1077  return NULL;
1078  }
1079 
1080  quantiles_count--;
1081  for (i = 0; i <= quantiles_count; i++)
1082  quantiles[i] = ((double) i) / quantiles_count;
1083  quantiles_count++;
1084  }
1085 
1086  /* check quantiles */
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);
1091  return NULL;
1092  }
1093  }
1094  quicksort(quantiles, quantiles + quantiles_count - 1);
1095 
1096  /* initialize linked-list set */
1097  *qlls_count = quantiles_count * 2;
1098  RASTER_DEBUGF(4, "qlls_count = %d", *qlls_count);
1099  *qlls = rtalloc(sizeof(struct quantile_llist) * *qlls_count);
1100  if (NULL == *qlls) {
1101  rterror("rt_band_get_quantiles_stream: Could not allocate memory for quantile output");
1102  if (init_quantiles) rtdealloc(quantiles);
1103  return NULL;
1104  }
1105 
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];
1110  qll->count = 0;
1111  qll->sum1 = 0;
1112  qll->sum2 = 0;
1113  qll->head = NULL;
1114  qll->tail = NULL;
1115 
1116  /* initialize index */
1117  qll->index = rtalloc(sizeof(struct quantile_llist_index) * j);
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);
1121  return NULL;
1122  }
1123  qll->index_max = j;
1124  quantile_llist_index_reset(qll);
1125 
1126  /* AL-GEQ */
1127  if (!(i % 2)) {
1128  qll->algeq = 1;
1129  qll->tau = (uint64_t) ROUND(cov_count - (cov_count * qll->quantile), 0);
1130  if (qll->tau < 1) qll->tau = 1;
1131  }
1132  /* AL-GT */
1133  else {
1134  qll->algeq = 0;
1135  qll->tau = cov_count - (*qlls)[i - 1].tau + 1;
1136  }
1137 
1138  RASTER_DEBUGF(4, "qll init: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
1139  qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
1140  RASTER_DEBUGF(4, "qll init: (head, tail) = (%p, %p)", qll->head, qll->tail);
1141  }
1142 
1143  if (init_quantiles) rtdealloc(quantiles);
1144  }
1145 
1146  /* clamp percentage */
1147  if (
1148  (sample < 0 || FLT_EQ(sample, 0.0)) ||
1149  (sample > 1 || FLT_EQ(sample, 1.0))
1150  ) {
1151  do_sample = 0;
1152  sample = 1;
1153  }
1154  else
1155  do_sample = 1;
1156  RASTER_DEBUGF(3, "do_sample = %d", do_sample);
1157 
1158  /* sample all pixels */
1159  if (!do_sample) {
1160  sample_size = band->width * band->height;
1161  sample_per = band->height;
1162  }
1163  /*
1164  randomly sample a percentage of available pixels
1165  sampling method is known as
1166  "systematic random sample without replacement"
1167  */
1168  else {
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);
1172  srand(time(NULL));
1173  }
1174  RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set"
1175  , sample_size, (band->width * band->height), sample_per);
1176 
1177  for (x = 0, j = 0, k = 0; x < band->width; x++) {
1178  y = -1;
1179  diff = 0;
1180 
1181  /* exclude_nodata_value = TRUE and band is NODATA */
1182  if (exclude_nodata_value && rt_band_get_isnodata_flag(band)) {
1183  RASTER_DEBUG(3, "Skipping quantile calcuation as band is NODATA");
1184  break;
1185  }
1186 
1187  for (i = 0, z = 0; i < sample_per; i++) {
1188  if (do_sample != 1)
1189  y = i;
1190  else {
1191  offset = (rand() % sample_int) + 1;
1192  y += diff + offset;
1193  diff = sample_int - offset;
1194  }
1195  RASTER_DEBUGF(5, "(x, y, z) = (%d, %d, %d)", x, y, z);
1196  if (y >= band->height || z > sample_per) break;
1197 
1198  status = rt_band_get_pixel(band, x, y, &value, &isnodata);
1199 
1200  j++;
1201  if (status == ES_NONE && (!exclude_nodata_value || (exclude_nodata_value && !isnodata))) {
1202 
1203  /* process each quantile */
1204  for (a = 0; a < *qlls_count; a++) {
1205  qll = &((*qlls)[a]);
1206  qls = NULL;
1207  RASTER_DEBUGF(4, "%d of %d (%f)", a + 1, *qlls_count, qll->quantile);
1208  RASTER_DEBUGF(5, "qll before: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
1209  qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
1210  RASTER_DEBUGF(5, "qll before: (head, tail) = (%p, %p)", qll->head, qll->tail);
1211 
1212  /* OPTIMIZATION: shortcuts for quantiles of zero or one */
1213  if (FLT_EQ(qll->quantile, 0.)) {
1214  if (NULL != qll->head) {
1215  if (value < qll->head->value)
1216  qll->head->value = value;
1217  }
1218  else {
1219  qle = quantile_llist_insert(qll->head, value, NULL);
1220  qll->head = qle;
1221  qll->tail = qle;
1222  qll->count = 1;
1223  }
1224 
1225  RASTER_DEBUGF(4, "quantile shortcut for %f\n\n", qll->quantile);
1226  continue;
1227  }
1228  else if (FLT_EQ(qll->quantile, 1.)) {
1229  if (NULL != qll->head) {
1230  if (value > qll->head->value)
1231  qll->head->value = value;
1232  }
1233  else {
1234  qle = quantile_llist_insert(qll->head, value, NULL);
1235  qll->head = qle;
1236  qll->tail = qle;
1237  qll->count = 1;
1238  }
1239 
1240  RASTER_DEBUGF(4, "quantile shortcut for %f\n\n", qll->quantile);
1241  continue;
1242  }
1243 
1244  /* value exists in list */
1245  /* OPTIMIZATION: check to see if value exceeds last value */
1246  if (NULL != qll->tail && value > qll->tail->value)
1247  qle = NULL;
1248  /* OPTIMIZATION: check to see if value equals last value */
1249  else if (NULL != qll->tail && FLT_EQ(value, qll->tail->value))
1250  qle = qll->tail;
1251  /* OPTIMIZATION: use index if possible */
1252  else {
1253  qls = quantile_llist_index_search(qll, value, &idx);
1254  qle = quantile_llist_search(qls, value);
1255  }
1256 
1257  /* value found */
1258  if (NULL != qle) {
1259  RASTER_DEBUGF(4, "%f found in list", value);
1260  RASTER_DEBUGF(5, "qle before: (value, count) = (%f, %d)", qle->value, qle->count);
1261 
1262  qle->count++;
1263  qll->sum1++;
1264 
1265  if (qll->algeq)
1266  qll->sum2 = qll->sum1 - qll->head->count;
1267  else
1268  qll->sum2 = qll->sum1 - qll->tail->count;
1269 
1270  RASTER_DEBUGF(4, "qle after: (value, count) = (%f, %d)", qle->value, qle->count);
1271  }
1272  /* can still add element */
1273  else if (qll->count < MAX_VALUES) {
1274  RASTER_DEBUGF(4, "Adding %f to list", value);
1275 
1276  /* insert */
1277  /* OPTIMIZATION: check to see if value exceeds last value */
1278  if (NULL != qll->tail && (value > qll->tail->value || FLT_EQ(value, qll->tail->value))) {
1279  idx = qll->count;
1280  qle = quantile_llist_insert(qll->tail, value, &idx);
1281  }
1282  /* OPTIMIZATION: use index if possible */
1283  else
1284  qle = quantile_llist_insert(qls, value, &idx);
1285  if (NULL == qle) return NULL;
1286  RASTER_DEBUGF(5, "value added at index: %d => %f", idx, value);
1287  qll->count++;
1288  qll->sum1++;
1289 
1290  /* first element */
1291  if (NULL == qle->prev)
1292  qll->head = qle;
1293  /* last element */
1294  if (NULL == qle->next)
1295  qll->tail = qle;
1296 
1297  if (qll->algeq)
1298  qll->sum2 = qll->sum1 - qll->head->count;
1299  else
1300  qll->sum2 = qll->sum1 - qll->tail->count;
1301 
1302  /* index is only needed if there are at least 100 values */
1303  quantile_llist_index_update(qll, qle, idx);
1304 
1305  RASTER_DEBUGF(5, "qle, prev, next, head, tail = %p, %p, %p, %p, %p", qle, qle->prev, qle->next, qll->head, qll->tail);
1306  }
1307  /* AL-GEQ */
1308  else if (qll->algeq) {
1309  RASTER_DEBUGF(4, "value, head->value = %f, %f", value, qll->head->value);
1310 
1311  if (value < qll->head->value) {
1312  /* ignore value if test isn't true */
1313  if (qll->sum1 >= qll->tau) {
1314  RASTER_DEBUGF(4, "skipping %f", value);
1315  }
1316  else {
1317 
1318  /* delete last element */
1319  RASTER_DEBUGF(4, "deleting %f from list", qll->tail->value);
1320  qle = qll->tail->prev;
1321  RASTER_DEBUGF(5, "to-be tail is %f with count %d", qle->value, qle->count);
1322  qle->count += qll->tail->count;
1323  quantile_llist_index_delete(qll, qll->tail);
1324  quantile_llist_delete(qll->tail);
1325  qll->tail = qle;
1326  qll->count--;
1327  RASTER_DEBUGF(5, "tail is %f with count %d", qll->tail->value, qll->tail->count);
1328 
1329  /* insert value */
1330  RASTER_DEBUGF(4, "adding %f to list", value);
1331  /* OPTIMIZATION: check to see if value exceeds last value */
1332  if (NULL != qll->tail && (value > qll->tail->value || FLT_EQ(value, qll->tail->value))) {
1333  idx = qll->count;
1334  qle = quantile_llist_insert(qll->tail, value, &idx);
1335  }
1336  /* OPTIMIZATION: use index if possible */
1337  else {
1338  qls = quantile_llist_index_search(qll, value, &idx);
1339  qle = quantile_llist_insert(qls, value, &idx);
1340  }
1341  if (NULL == qle) return NULL;
1342  RASTER_DEBUGF(5, "value added at index: %d => %f", idx, value);
1343  qll->count++;
1344  qll->sum1++;
1345 
1346  /* first element */
1347  if (NULL == qle->prev)
1348  qll->head = qle;
1349  /* last element */
1350  if (NULL == qle->next)
1351  qll->tail = qle;
1352 
1353  qll->sum2 = qll->sum1 - qll->head->count;
1354 
1355  quantile_llist_index_update(qll, qle, idx);
1356 
1357  RASTER_DEBUGF(5, "qle, head, tail = %p, %p, %p", qle, qll->head, qll->tail);
1358 
1359  }
1360  }
1361  else {
1362  qle = qll->tail;
1363  while (NULL != qle) {
1364  if (qle->value < value) {
1365  qle->count++;
1366  qll->sum1++;
1367  qll->sum2 = qll->sum1 - qll->head->count;
1368  RASTER_DEBUGF(4, "incremented count of %f by 1 to %d", qle->value, qle->count);
1369  break;
1370  }
1371 
1372  qle = qle->prev;
1373  }
1374  }
1375  }
1376  /* AL-GT */
1377  else {
1378  RASTER_DEBUGF(4, "value, tail->value = %f, %f", value, qll->tail->value);
1379 
1380  if (value > qll->tail->value) {
1381  /* ignore value if test isn't true */
1382  if (qll->sum1 >= qll->tau) {
1383  RASTER_DEBUGF(4, "skipping %f", value);
1384  }
1385  else {
1386 
1387  /* delete last element */
1388  RASTER_DEBUGF(4, "deleting %f from list", qll->head->value);
1389  qle = qll->head->next;
1390  RASTER_DEBUGF(5, "to-be tail is %f with count %d", qle->value, qle->count);
1391  qle->count += qll->head->count;
1392  quantile_llist_index_delete(qll, qll->head);
1393  quantile_llist_delete(qll->head);
1394  qll->head = qle;
1395  qll->count--;
1396  quantile_llist_index_update(qll, NULL, 0);
1397  RASTER_DEBUGF(5, "tail is %f with count %d", qll->head->value, qll->head->count);
1398 
1399  /* insert value */
1400  RASTER_DEBUGF(4, "adding %f to list", value);
1401  /* OPTIMIZATION: check to see if value exceeds last value */
1402  if (NULL != qll->tail && (value > qll->tail->value || FLT_EQ(value, qll->tail->value))) {
1403  idx = qll->count;
1404  qle = quantile_llist_insert(qll->tail, value, &idx);
1405  }
1406  /* OPTIMIZATION: use index if possible */
1407  else {
1408  qls = quantile_llist_index_search(qll, value, &idx);
1409  qle = quantile_llist_insert(qls, value, &idx);
1410  }
1411  if (NULL == qle) return NULL;
1412  RASTER_DEBUGF(5, "value added at index: %d => %f", idx, value);
1413  qll->count++;
1414  qll->sum1++;
1415 
1416  /* first element */
1417  if (NULL == qle->prev)
1418  qll->head = qle;
1419  /* last element */
1420  if (NULL == qle->next)
1421  qll->tail = qle;
1422 
1423  qll->sum2 = qll->sum1 - qll->tail->count;
1424 
1425  quantile_llist_index_update(qll, qle, idx);
1426 
1427  RASTER_DEBUGF(5, "qle, head, tail = %p, %p, %p", qle, qll->head, qll->tail);
1428 
1429  }
1430  }
1431  else {
1432  qle = qll->head;
1433  while (NULL != qle) {
1434  if (qle->value > value) {
1435  qle->count++;
1436  qll->sum1++;
1437  qll->sum2 = qll->sum1 - qll->tail->count;
1438  RASTER_DEBUGF(4, "incremented count of %f by 1 to %d", qle->value, qle->count);
1439  break;
1440  }
1441 
1442  qle = qle->next;
1443  }
1444  }
1445  }
1446 
1447  RASTER_DEBUGF(5, "sum2, tau = %d, %d", qll->sum2, qll->tau);
1448  if (qll->sum2 >= qll->tau) {
1449  /* AL-GEQ */
1450  if (qll->algeq) {
1451  RASTER_DEBUGF(4, "deleting first element %f from list", qll->head->value);
1452 
1453  if (NULL != qll->head->next) {
1454  qle = qll->head->next;
1455  qll->sum1 -= qll->head->count;
1456  qll->sum2 = qll->sum1 - qle->count;
1457  quantile_llist_index_delete(qll, qll->head);
1458  quantile_llist_delete(qll->head);
1459  qll->head = qle;
1460  qll->count--;
1461 
1462  quantile_llist_index_update(qll, NULL, 0);
1463  }
1464  else {
1465  quantile_llist_delete(qll->head);
1466  qll->head = NULL;
1467  qll->tail = NULL;
1468  qll->sum1 = 0;
1469  qll->sum2 = 0;
1470  qll->count = 0;
1471 
1472  quantile_llist_index_reset(qll);
1473  }
1474  }
1475  /* AL-GT */
1476  else {
1477  RASTER_DEBUGF(4, "deleting first element %f from list", qll->tail->value);
1478 
1479  if (NULL != qll->tail->prev) {
1480  qle = qll->tail->prev;
1481  qll->sum1 -= qll->tail->count;
1482  qll->sum2 = qll->sum1 - qle->count;
1483  quantile_llist_index_delete(qll, qll->tail);
1484  quantile_llist_delete(qll->tail);
1485  qll->tail = qle;
1486  qll->count--;
1487  }
1488  else {
1489  quantile_llist_delete(qll->tail);
1490  qll->head = NULL;
1491  qll->tail = NULL;
1492  qll->sum1 = 0;
1493  qll->sum2 = 0;
1494  qll->count = 0;
1495 
1496  quantile_llist_index_reset(qll);
1497  }
1498  }
1499  }
1500 
1501  RASTER_DEBUGF(5, "qll after: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
1502  qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
1503  RASTER_DEBUGF(5, "qll after: (head, tail) = (%p, %p)\n\n", qll->head, qll->tail);
1504  }
1505 
1506  }
1507  else {
1508  RASTER_DEBUGF(5, "skipping value at (x, y) = (%d, %d)", x, y);
1509  }
1510 
1511  z++;
1512  }
1513  }
1514 
1515  /* process quantiles */
1516  *rtn_count = *qlls_count / 2;
1517  rtn = rtalloc(sizeof(struct rt_quantile_t) * *rtn_count);
1518  if (NULL == rtn) return NULL;
1519 
1520  RASTER_DEBUGF(3, "returning %d quantiles", *rtn_count);
1521  for (i = 0, k = 0; i < *qlls_count; i++) {
1522  qll = &((*qlls)[i]);
1523 
1524  exists = 0;
1525  for (x = 0; x < k; x++) {
1526  if (FLT_EQ(qll->quantile, rtn[x].quantile)) {
1527  exists = 1;
1528  break;
1529  }
1530  }
1531  if (exists) continue;
1532 
1533  RASTER_DEBUGF(5, "qll: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
1534  qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
1535  RASTER_DEBUGF(5, "qll: (head, tail) = (%p, %p)", qll->head, qll->tail);
1536 
1537  rtn[k].quantile = qll->quantile;
1538  rtn[k].has_value = 0;
1539 
1540  /* check that qll->head and qll->tail have value */
1541  if (qll->head == NULL || qll->tail == NULL)
1542  continue;
1543 
1544  /* AL-GEQ */
1545  if (qll->algeq)
1546  qle = qll->head;
1547  /* AM-GT */
1548  else
1549  qle = qll->tail;
1550 
1551  exists = 0;
1552  for (j = i + 1; j < *qlls_count; j++) {
1553  if (FLT_EQ((*qlls)[j].quantile, qll->quantile)) {
1554 
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);
1558 
1559  exists = 1;
1560  break;
1561  }
1562  }
1563 
1564  /* weighted average for quantile */
1565  if (exists) {
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);
1569  }
1570  else {
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);
1573  }
1574  }
1575  /* straight value for quantile */
1576  else {
1577  rtn[k].value = qle->value;
1578  }
1579  rtn[k].has_value = 1;
1580  RASTER_DEBUGF(3, "(quantile, value) = (%f, %f)\n\n", rtn[k].quantile, rtn[k].value);
1581 
1582  k++;
1583  }
1584 
1585  RASTER_DEBUG(3, "done");
1586  return rtn;
1587 }
quantile_llist::sum1
uint64_t sum1
Definition: librtcore.h:2356
quantile_llist_delete
static int quantile_llist_delete(struct quantile_llist_element *element)
Definition: rt_statistics.c:855
quantile_llist::quantile
double quantile
Definition: librtcore.h:2345
quantile_llist_element::prev
struct quantile_llist_element * prev
Definition: librtcore.h:2364
rt_quantile_t::quantile
double quantile
Definition: librtcore.h:2337
quantile_llist_element::value
double value
Definition: librtcore.h:2361
quantile_llist_index_reset
static void quantile_llist_index_reset(struct quantile_llist *qll)
Definition: rt_statistics.c:971
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2185
rterror
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
quantile_llist_index_search
static struct quantile_llist_element * quantile_llist_index_search(struct quantile_llist *qll, double value, uint32_t *index)
Definition: rt_statistics.c:936
ovdump.data
data
Definition: ovdump.py:103
quantile_llist_element::count
uint32_t count
Definition: librtcore.h:2362
rt_band_t::height
uint16_t height
Definition: librtcore.h:2268
rt_quantile_t::value
double value
Definition: librtcore.h:2338
quantile_llist_element::next
struct quantile_llist_element * next
Definition: librtcore.h:2365
uint32_t
unsigned int uint32_t
Definition: uthash.h:78
quantile_llist_element
Definition: librtcore.h:2360
quantile_llist::tail
struct quantile_llist_element * tail
Definition: librtcore.h:2349
quicksort
static void quicksort(double *left, double *right)
Definition: rt_statistics.c:81
quantile_llist
Definition: librtcore.h:2343
rt_band_get_pixel
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1088
quantile_llist::tau
uint64_t tau
Definition: librtcore.h:2346
quantile_llist_index_delete
static void quantile_llist_index_delete(struct quantile_llist *qll, struct quantile_llist_element *qle)
Definition: rt_statistics.c:919
rt_quantile_t
Definition: librtcore.h:2336
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
genraster.count
int count
Definition: genraster.py:56
rt_band_get_data
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition: rt_band.c:302
quantile_llist::index_max
uint32_t index_max
Definition: librtcore.h:2354
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
rt_band_t::width
uint16_t width
Definition: librtcore.h:2267
quantile_llist_index_update
static void quantile_llist_index_update(struct quantile_llist *qll, struct quantile_llist_element *qle, uint32_t idx)
Definition: rt_statistics.c:898
rtdealloc
void rtdealloc(void *mem)
Definition: rt_context.c:186
ROUND
#define ROUND(x, y)
Definition: librtcore.h:2192
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
quantile_llist_insert
static struct quantile_llist_element * quantile_llist_insert(struct quantile_llist_element *element, double value, uint32_t *idx)
Definition: rt_statistics.c:796
quantile_llist::index
struct quantile_llist_index * index
Definition: librtcore.h:2353
rt_quantile_t::has_value
uint32_t has_value
Definition: librtcore.h:2339
quantile_llist::algeq
uint8_t algeq
Definition: librtcore.h:2344
quantile_llist::count
uint32_t count
Definition: librtcore.h:2350
quantile_llist_search
static struct quantile_llist_element * quantile_llist_search(struct quantile_llist_element *element, double needle)
Definition: rt_statistics.c:780
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:581
genraster.value
int value
Definition: genraster.py:61
quantile_llist::head
struct quantile_llist_element * head
Definition: librtcore.h:2348
quantile_llist_index
Definition: librtcore.h:2368
uint8_t
unsigned char uint8_t
Definition: uthash.h:79
quantile_llist::sum2
uint64_t sum2
Definition: librtcore.h:2357
Here is the call graph for this function:
Here is the caller graph for this function: