PostGIS  2.5.7dev-r@@SVN_REVISION@@
rt_statistics.c
Go to the documentation of this file.
1 /*
2  *
3  * WKTRaster - Raster Types for PostGIS
4  * http://trac.osgeo.org/postgis/wiki/WKTRaster
5  *
6  * Copyright (C) 2011-2013 Regents of the University of California
7  * <bkpark@ucdavis.edu>
8  * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
9  * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
10  * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
11  * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
12  * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
13  *
14  * This program is free software; you can redistribute it and/or
15  * modify it under the terms of the GNU General Public License
16  * as published by the Free Software Foundation; either version 2
17  * of the License, or (at your option) any later version.
18  *
19  * This program is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with this program; if not, write to the Free Software Foundation,
26  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
27  *
28  */
29 
30 #include "librtcore.h"
31 #include "librtcore_internal.h"
32 
33 /******************************************************************************
34 * quicksort
35 ******************************************************************************/
36 
37 #define SWAP(x, y) { double t; t = x; x = y; y = t; }
38 #define ORDER(x, y) if (x > y) SWAP(x, y)
39 
40 static double pivot(double *left, double *right) {
41  double l, m, r, *p;
42 
43  l = *left;
44  m = *(left + (right - left) / 2);
45  r = *right;
46 
47  /* order */
48  ORDER(l, m);
49  ORDER(l, r);
50  ORDER(m, r);
51 
52  /* pivot is higher of two values */
53  if (l < m) return m;
54  if (m < r) return r;
55 
56  /* find pivot that isn't left */
57  for (p = left + 1; p <= right; ++p) {
58  if (*p != *left)
59  return (*p < *left) ? *left : *p;
60  }
61 
62  /* all values are same */
63  return -1;
64 }
65 
66 static double *partition(double *left, double *right, double pivot) {
67  while (left <= right) {
68  while (*left < pivot) ++left;
69  while (*right >= pivot) --right;
70 
71  if (left < right) {
72  SWAP(*left, *right);
73  ++left;
74  --right;
75  }
76  }
77 
78  return left;
79 }
80 
81 static void quicksort(double *left, double *right) {
82  double p = pivot(left, right);
83  double *pos;
84 
85  if (p != -1) {
86  pos = partition(left, right, p);
87  quicksort(left, pos - 1);
88  quicksort(pos, right);
89  }
90 }
91 
92 /******************************************************************************
93 * rt_band_get_summary_stats()
94 ******************************************************************************/
95 
109 rt_bandstats
110 rt_band_get_summary_stats(
111  rt_band band,
112  int exclude_nodata_value, double sample, int inc_vals,
113  uint64_t *cK, double *cM, double *cQ
114 ) {
115  uint32_t x = 0;
116  uint32_t y = 0;
117  uint32_t z = 0;
118  uint32_t offset = 0;
119  uint32_t diff = 0;
120  int rtn;
121  int hasnodata = FALSE;
122  double nodata = 0;
123  double *values = NULL;
124  double value;
125  int isnodata = 0;
126  rt_bandstats stats = NULL;
127 
128  uint32_t do_sample = 0;
129  uint32_t sample_size = 0;
130  uint32_t sample_per = 0;
131  uint32_t sample_int = 0;
132  uint32_t i = 0;
133  uint32_t j = 0;
134  double sum = 0;
135  uint32_t k = 0;
136  double M = 0;
137  double Q = 0;
138 
139 #if POSTGIS_DEBUG_LEVEL > 0
140  clock_t start, stop;
141  double elapsed = 0;
142 #endif
143 
144  RASTER_DEBUG(3, "starting");
145 #if POSTGIS_DEBUG_LEVEL > 0
146  start = clock();
147 #endif
148 
149  assert(NULL != band);
150 
151  /* band is empty (width < 1 || height < 1) */
152  if (band->width < 1 || band->height < 1) {
153  stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
154  if (NULL == stats) {
155  rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
156  return NULL;
157  }
158 
159  rtwarn("Band is empty as width and/or height is 0");
160 
161  stats->sample = 1;
162  stats->sorted = 0;
163  stats->values = NULL;
164  stats->count = 0;
165  stats->min = stats->max = 0;
166  stats->sum = 0;
167  stats->mean = 0;
168  stats->stddev = -1;
169 
170  return stats;
171  }
172 
173  hasnodata = rt_band_get_hasnodata_flag(band);
174  if (hasnodata != FALSE)
175  rt_band_get_nodata(band, &nodata);
176  else
177  exclude_nodata_value = 0;
178 
179  RASTER_DEBUGF(3, "nodata = %f", nodata);
180  RASTER_DEBUGF(3, "hasnodata = %d", hasnodata);
181  RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
182 
183  /* entire band is nodata */
184  if (rt_band_get_isnodata_flag(band) != FALSE) {
185  stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
186  if (NULL == stats) {
187  rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
188  return NULL;
189  }
190 
191  stats->sample = 1;
192  stats->sorted = 0;
193  stats->values = NULL;
194 
195  if (exclude_nodata_value) {
196  rtwarn("All pixels of band have the NODATA value");
197 
198  stats->count = 0;
199  stats->min = stats->max = 0;
200  stats->sum = 0;
201  stats->mean = 0;
202  stats->stddev = -1;
203  }
204  else {
205  stats->count = band->width * band->height;
206  stats->min = stats->max = nodata;
207  stats->sum = stats->count * nodata;
208  stats->mean = nodata;
209  stats->stddev = 0;
210  }
211 
212  return stats;
213  }
214 
215  /* clamp percentage */
216  if (
217  (sample < 0 || FLT_EQ(sample, 0.0)) ||
218  (sample > 1 || FLT_EQ(sample, 1.0))
219  ) {
220  do_sample = 0;
221  sample = 1;
222  }
223  else
224  do_sample = 1;
225  RASTER_DEBUGF(3, "do_sample = %d", do_sample);
226 
227  /* sample all pixels */
228  if (!do_sample) {
229  sample_size = band->width * band->height;
230  sample_per = band->height;
231  }
232  /*
233  randomly sample a percentage of available pixels
234  sampling method is known as
235  "systematic random sample without replacement"
236  */
237  else {
238  sample_size = round((band->width * band->height) * sample);
239  sample_per = round(sample_size / band->width);
240  if (sample_per < 1)
241  sample_per = 1;
242  sample_int = round(band->height / sample_per);
243  srand(time(NULL));
244  }
245 
246  RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set"
247  , sample_size, (band->width * band->height), sample_per);
248 
249  if (inc_vals) {
250  values = rtalloc(sizeof(double) * sample_size);
251  if (NULL == values) {
252  rtwarn("Could not allocate memory for values");
253  inc_vals = 0;
254  }
255  }
256 
257  /* initialize stats */
258  stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
259  if (NULL == stats) {
260  rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
261  return NULL;
262  }
263  stats->sample = sample;
264  stats->count = 0;
265  stats->sum = 0;
266  stats->mean = 0;
267  stats->stddev = -1;
268  stats->min = stats->max = 0;
269  stats->values = NULL;
270  stats->sorted = 0;
271 
272  for (x = 0, j = 0, k = 0; x < band->width; x++) {
273  y = -1;
274  diff = 0;
275 
276  for (i = 0, z = 0; i < sample_per; i++) {
277  if (!do_sample)
278  y = i;
279  else {
280  offset = (rand() % sample_int) + 1;
281  y += diff + offset;
282  diff = sample_int - offset;
283  }
284  RASTER_DEBUGF(5, "(x, y, z) = (%d, %d, %d)", x, y, z);
285  if (y >= band->height || z > sample_per) break;
286 
287  rtn = rt_band_get_pixel(band, x, y, &value, &isnodata);
288 
289  j++;
290  if (rtn == ES_NONE && (!exclude_nodata_value || (exclude_nodata_value && !isnodata))) {
291 
292  /* inc_vals set, collect pixel values */
293  if (inc_vals) values[k] = value;
294 
295  /* average */
296  k++;
297  sum += value;
298 
299  /*
300  one-pass standard deviation
301  http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
302  */
303  if (k == 1) {
304  Q = 0;
305  M = value;
306  }
307  else {
308  Q += (((k - 1) * pow(value - M, 2)) / k);
309  M += ((value - M ) / k);
310  }
311 
312  /* coverage one-pass standard deviation */
313  if (NULL != cK) {
314  (*cK)++;
315  if (*cK == 1) {
316  *cQ = 0;
317  *cM = value;
318  }
319  else {
320  *cQ += (((*cK - 1) * pow(value - *cM, 2)) / *cK);
321  *cM += ((value - *cM ) / *cK);
322  }
323  }
324 
325  /* min/max */
326  if (stats->count < 1) {
327  stats->count = 1;
328  stats->min = stats->max = value;
329  }
330  else {
331  if (value < stats->min)
332  stats->min = value;
333  if (value > stats->max)
334  stats->max = value;
335  }
336 
337  }
338 
339  z++;
340  }
341  }
342 
343  RASTER_DEBUG(3, "sampling complete");
344 
345  stats->count = k;
346  if (k > 0) {
347  if (inc_vals) {
348  /* free unused memory */
349  if (sample_size != k) {
350  values = rtrealloc(values, k * sizeof(double));
351  }
352 
353  stats->values = values;
354  }
355 
356  stats->sum = sum;
357  stats->mean = sum / k;
358 
359  /* standard deviation */
360  if (!do_sample)
361  stats->stddev = sqrt(Q / k);
362  /* sample deviation */
363  else {
364  if (k < 2)
365  stats->stddev = -1;
366  else
367  stats->stddev = sqrt(Q / (k - 1));
368  }
369  }
370  /* inc_vals thus values allocated but not used */
371  else if (inc_vals)
372  rtdealloc(values);
373 
374  /* if do_sample is one */
375  if (do_sample && k < 1)
376  rtwarn("All sampled pixels of band have the NODATA value");
377 
378 #if POSTGIS_DEBUG_LEVEL > 0
379  stop = clock();
380  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
381  RASTER_DEBUGF(3, "(time, count, mean, stddev, min, max) = (%0.4f, %d, %f, %f, %f, %f)",
382  elapsed, stats->count, stats->mean, stats->stddev, stats->min, stats->max);
383 #endif
384 
385  RASTER_DEBUG(3, "done");
386  return stats;
387 }
388 
389 /******************************************************************************
390 * rt_band_get_histogram()
391 ******************************************************************************/
392 
411 rt_histogram
412 rt_band_get_histogram(
413  rt_bandstats stats,
414  uint32_t bin_count, double *bin_width, uint32_t bin_width_count,
415  int right, double min, double max,
416  uint32_t *rtn_count
417 ) {
418  rt_histogram bins = NULL;
419  int init_width = 0;
420  uint32_t i;
421  uint32_t j;
422  double tmp;
423  double value;
424  int sum = 0;
425  double qmin;
426  double qmax;
427 
428 #if POSTGIS_DEBUG_LEVEL > 0
429  clock_t start, stop;
430  double elapsed = 0;
431 #endif
432 
433  RASTER_DEBUG(3, "starting");
434 #if POSTGIS_DEBUG_LEVEL > 0
435  start = clock();
436 #endif
437 
438  assert(NULL != stats);
439  assert(NULL != rtn_count);
440 
441  if (stats->count < 1 || NULL == stats->values) {
442  rterror("rt_util_get_histogram: rt_bandstats object has no value");
443  return NULL;
444  }
445 
446  /* bin width must be positive numbers and not zero */
447  if (NULL != bin_width && bin_width_count > 0) {
448  for (i = 0; i < bin_width_count; i++) {
449  if (bin_width[i] < 0 || FLT_EQ(bin_width[i], 0.0)) {
450  rterror("rt_util_get_histogram: bin_width element is less than or equal to zero");
451  return NULL;
452  }
453  }
454  }
455 
456  /* ignore min and max parameters */
457  if (FLT_EQ(max, min)) {
458  qmin = stats->min;
459  qmax = stats->max;
460  }
461  else {
462  qmin = min;
463  qmax = max;
464  if (qmin > qmax) {
465  qmin = max;
466  qmax = min;
467  }
468  }
469 
470  /* # of bins not provided */
471  if (bin_count <= 0) {
472  /*
473  determine # of bins
474  http://en.wikipedia.org/wiki/Histogram
475 
476  all computed bins are assumed to have equal width
477  */
478  /* Square-root choice for stats->count < 30 */
479  if (stats->count < 30)
480  bin_count = ceil(sqrt(stats->count));
481  /* Sturges' formula for stats->count >= 30 */
482  else
483  bin_count = ceil(log2((double) stats->count) + 1.);
484 
485  /* bin_width_count provided and bin_width has value */
486  if (bin_width_count > 0 && NULL != bin_width) {
487  /* user has defined something specific */
488  if (bin_width_count > bin_count)
489  bin_count = bin_width_count;
490  else if (bin_width_count > 1) {
491  tmp = 0;
492  for (i = 0; i < bin_width_count; i++) tmp += bin_width[i];
493  bin_count = ceil((qmax - qmin) / tmp) * bin_width_count;
494  }
495  else
496  bin_count = ceil((qmax - qmin) / bin_width[0]);
497  }
498  /* set bin width count to zero so that one can be calculated */
499  else {
500  bin_width_count = 0;
501  }
502  }
503 
504  /* min and max the same */
505  if (FLT_EQ(qmax, qmin)) bin_count = 1;
506 
507  RASTER_DEBUGF(3, "bin_count = %d", bin_count);
508 
509  /* bin count = 1, all values are in one bin */
510  if (bin_count < 2) {
511  bins = rtalloc(sizeof(struct rt_histogram_t));
512  if (NULL == bins) {
513  rterror("rt_util_get_histogram: Could not allocate memory for histogram");
514  return NULL;
515  }
516 
517  bins->count = stats->count;
518  bins->percent = -1;
519  bins->min = qmin;
520  bins->max = qmax;
521  bins->inc_min = bins->inc_max = 1;
522 
523  *rtn_count = bin_count;
524  return bins;
525  }
526 
527  /* establish bin width */
528  if (bin_width_count == 0) {
529  bin_width_count = 1;
530 
531  /* bin_width unallocated */
532  if (NULL == bin_width) {
533  bin_width = rtalloc(sizeof(double));
534  if (NULL == bin_width) {
535  rterror("rt_util_get_histogram: Could not allocate memory for bin widths");
536  return NULL;
537  }
538  init_width = 1;
539  }
540 
541  bin_width[0] = (qmax - qmin) / bin_count;
542  }
543 
544  /* initialize bins */
545  bins = rtalloc(bin_count * sizeof(struct rt_histogram_t));
546  if (NULL == bins) {
547  rterror("rt_util_get_histogram: Could not allocate memory for histogram");
548  if (init_width) rtdealloc(bin_width);
549  return NULL;
550  }
551  if (!right)
552  tmp = qmin;
553  else
554  tmp = qmax;
555  for (i = 0; i < bin_count;) {
556  for (j = 0; j < bin_width_count; j++) {
557  bins[i].count = 0;
558  bins->percent = -1;
559 
560  if (!right) {
561  bins[i].min = tmp;
562  tmp += bin_width[j];
563  bins[i].max = tmp;
564 
565  bins[i].inc_min = 1;
566  bins[i].inc_max = 0;
567  }
568  else {
569  bins[i].max = tmp;
570  tmp -= bin_width[j];
571  bins[i].min = tmp;
572 
573  bins[i].inc_min = 0;
574  bins[i].inc_max = 1;
575  }
576 
577  i++;
578  }
579  }
580  if (!right) {
581  bins[bin_count - 1].inc_max = 1;
582 
583  /* align last bin to the max value */
584  if (bins[bin_count - 1].max < qmax)
585  bins[bin_count - 1].max = qmax;
586  }
587  else {
588  bins[bin_count - 1].inc_min = 1;
589 
590  /* align first bin to the min value */
591  if (bins[bin_count - 1].min > qmin)
592  bins[bin_count - 1].min = qmin;
593  }
594 
595  /* process the values */
596  for (i = 0; i < stats->count; i++) {
597  value = stats->values[i];
598 
599  /* default, [a, b) */
600  if (!right) {
601  for (j = 0; j < bin_count; j++) {
602  if (
603  (!bins[j].inc_max && value < bins[j].max) || (
604  bins[j].inc_max && (
605  (value < bins[j].max) ||
606  FLT_EQ(value, bins[j].max)
607  )
608  )
609  ) {
610  bins[j].count++;
611  sum++;
612  break;
613  }
614  }
615  }
616  /* (a, b] */
617  else {
618  for (j = 0; j < bin_count; j++) {
619  if (
620  (!bins[j].inc_min && value > bins[j].min) || (
621  bins[j].inc_min && (
622  (value > bins[j].min) ||
623  FLT_EQ(value, bins[j].min)
624  )
625  )
626  ) {
627  bins[j].count++;
628  sum++;
629  break;
630  }
631  }
632  }
633  }
634 
635  for (i = 0; i < bin_count; i++) {
636  bins[i].percent = ((double) bins[i].count) / sum;
637  }
638 
639 #if POSTGIS_DEBUG_LEVEL > 0
640  stop = clock();
641  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
642  RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
643 
644  for (j = 0; j < bin_count; j++) {
645  RASTER_DEBUGF(5, "(min, max, inc_min, inc_max, count, sum, percent) = (%f, %f, %d, %d, %d, %d, %f)",
646  bins[j].min, bins[j].max, bins[j].inc_min, bins[j].inc_max, bins[j].count, sum, bins[j].percent);
647  }
648 #endif
649 
650  if (init_width) rtdealloc(bin_width);
651  *rtn_count = bin_count;
652  RASTER_DEBUG(3, "done");
653  return bins;
654 }
655 
656 /******************************************************************************
657 * rt_band_get_quantiles()
658 ******************************************************************************/
659 
671 rt_quantile
672 rt_band_get_quantiles(
673  rt_bandstats stats,
674  double *quantiles, int quantiles_count,
675  uint32_t *rtn_count
676 ) {
677  rt_quantile rtn;
678  int init_quantiles = 0;
679  int i = 0;
680  double h;
681  int hl;
682 
683 #if POSTGIS_DEBUG_LEVEL > 0
684  clock_t start, stop;
685  double elapsed = 0;
686 #endif
687 
688  RASTER_DEBUG(3, "starting");
689 #if POSTGIS_DEBUG_LEVEL > 0
690  start = clock();
691 #endif
692 
693  assert(NULL != stats);
694  assert(NULL != rtn_count);
695 
696  if (stats->count < 1 || NULL == stats->values) {
697  rterror("rt_band_get_quantiles: rt_bandstats object has no value");
698  return NULL;
699  }
700 
701  /* quantiles not provided */
702  if (NULL == quantiles) {
703  /* quantile count not specified, default to quartiles */
704  if (quantiles_count < 2)
705  quantiles_count = 5;
706 
707  quantiles = rtalloc(sizeof(double) * quantiles_count);
708  init_quantiles = 1;
709  if (NULL == quantiles) {
710  rterror("rt_band_get_quantiles: Could not allocate memory for quantile input");
711  return NULL;
712  }
713 
714  quantiles_count--;
715  for (i = 0; i <= quantiles_count; i++)
716  quantiles[i] = ((double) i) / quantiles_count;
717  quantiles_count++;
718  }
719 
720  /* check quantiles */
721  for (i = 0; i < quantiles_count; i++) {
722  if (quantiles[i] < 0. || quantiles[i] > 1.) {
723  rterror("rt_band_get_quantiles: Quantile value not between 0 and 1");
724  if (init_quantiles) rtdealloc(quantiles);
725  return NULL;
726  }
727  }
728  quicksort(quantiles, quantiles + quantiles_count - 1);
729 
730  /* initialize rt_quantile */
731  rtn = rtalloc(sizeof(struct rt_quantile_t) * quantiles_count);
732  if (NULL == rtn) {
733  rterror("rt_band_get_quantiles: Could not allocate memory for quantile output");
734  if (init_quantiles) rtdealloc(quantiles);
735  return NULL;
736  }
737 
738  /* sort values */
739  if (!stats->sorted) {
740  quicksort(stats->values, stats->values + stats->count - 1);
741  stats->sorted = 1;
742  }
743 
744  /*
745  make quantiles
746 
747  formula is that used in R (method 7) and Excel from
748  http://en.wikipedia.org/wiki/Quantile
749  */
750  for (i = 0; i < quantiles_count; i++) {
751  rtn[i].quantile = quantiles[i];
752 
753  h = ((stats->count - 1.) * quantiles[i]) + 1.;
754  hl = floor(h);
755 
756  /* h greater than hl, do full equation */
757  if (h > hl)
758  rtn[i].value = stats->values[hl - 1] + ((h - hl) * (stats->values[hl] - stats->values[hl - 1]));
759  /* shortcut as second part of equation is zero */
760  else
761  rtn[i].value = stats->values[hl - 1];
762  }
763 
764 #if POSTGIS_DEBUG_LEVEL > 0
765  stop = clock();
766  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
767  RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
768 #endif
769 
770  *rtn_count = quantiles_count;
771  if (init_quantiles) rtdealloc(quantiles);
772  RASTER_DEBUG(3, "done");
773  return rtn;
774 }
775 
776 /******************************************************************************
777 * rt_band_get_quantiles_stream()
778 ******************************************************************************/
779 
780 static struct quantile_llist_element *quantile_llist_search(
781  struct quantile_llist_element *element,
782  double needle
783 ) {
784  if (NULL == element)
785  return NULL;
786  else if (FLT_NEQ(needle, element->value)) {
787  if (NULL != element->next)
788  return quantile_llist_search(element->next, needle);
789  else
790  return NULL;
791  }
792  else
793  return element;
794 }
795 
796 static struct quantile_llist_element *quantile_llist_insert(
797  struct quantile_llist_element *element,
798  double value,
799  uint32_t *idx
800 ) {
801  struct quantile_llist_element *qle = NULL;
802 
803  if (NULL == element) {
804  qle = rtalloc(sizeof(struct quantile_llist_element));
805  RASTER_DEBUGF(4, "qle @ %p is only element in list", qle);
806  if (NULL == qle) return NULL;
807 
808  qle->value = value;
809  qle->count = 1;
810 
811  qle->prev = NULL;
812  qle->next = NULL;
813 
814  if (NULL != idx) *idx = 0;
815  return qle;
816  }
817  else if (value > element->value) {
818  if (NULL != idx) *idx += 1;
819  if (NULL != element->next)
820  return quantile_llist_insert(element->next, value, idx);
821  /* insert as last element in list */
822  else {
823  qle = rtalloc(sizeof(struct quantile_llist_element));
824  RASTER_DEBUGF(4, "insert qle @ %p as last element", qle);
825  if (NULL == qle) return NULL;
826 
827  qle->value = value;
828  qle->count = 1;
829 
830  qle->prev = element;
831  qle->next = NULL;
832  element->next = qle;
833 
834  return qle;
835  }
836  }
837  /* insert before current element */
838  else {
839  qle = rtalloc(sizeof(struct quantile_llist_element));
840  RASTER_DEBUGF(4, "insert qle @ %p before current element", qle);
841  if (NULL == qle) return NULL;
842 
843  qle->value = value;
844  qle->count = 1;
845 
846  if (NULL != element->prev) element->prev->next = qle;
847  qle->next = element;
848  qle->prev = element->prev;
849  element->prev = qle;
850 
851  return qle;
852  }
853 }
854 
855 static int quantile_llist_delete(struct quantile_llist_element *element) {
856  if (NULL == element) return 0;
857 
858  /* beginning of list */
859  if (NULL == element->prev && NULL != element->next) {
860  element->next->prev = NULL;
861  }
862  /* end of list */
863  else if (NULL != element->prev && NULL == element->next) {
864  element->prev->next = NULL;
865  }
866  /* within list */
867  else if (NULL != element->prev && NULL != element->next) {
868  element->prev->next = element->next;
869  element->next->prev = element->prev;
870  }
871 
872  RASTER_DEBUGF(4, "qle @ %p destroyed", element);
873  rtdealloc(element);
874 
875  return 1;
876 }
877 
878 int quantile_llist_destroy(struct quantile_llist **list, uint32_t list_count) {
879  struct quantile_llist_element *element = NULL;
880  uint32_t i;
881 
882  if (NULL == *list) return 0;
883 
884  for (i = 0; i < list_count; i++) {
885  element = (*list)[i].head;
886  while (NULL != element->next) {
887  quantile_llist_delete(element->next);
888  }
889  quantile_llist_delete(element);
890 
891  rtdealloc((*list)[i].index);
892  }
893 
894  rtdealloc(*list);
895  return 1;
896 }
897 
898 static void quantile_llist_index_update(struct quantile_llist *qll, struct quantile_llist_element *qle, uint32_t idx) {
899  uint32_t anchor = (uint32_t) floor(idx / 100);
900 
901  if (qll->tail == qle) return;
902 
903  if (
904  (anchor != 0) && (
905  NULL == qll->index[anchor].element ||
906  idx <= qll->index[anchor].index
907  )
908  ) {
909  qll->index[anchor].index = idx;
910  qll->index[anchor].element = qle;
911  }
912 
913  if (anchor != 0 && NULL == qll->index[0].element) {
914  qll->index[0].index = 0;
915  qll->index[0].element = qll->head;
916  }
917 }
918 
919 static void quantile_llist_index_delete(struct quantile_llist *qll, struct quantile_llist_element *qle) {
920  uint32_t i = 0;
921 
922  for (i = 0; i < qll->index_max; i++) {
923  if (
924  NULL == qll->index[i].element ||
925  (qll->index[i].element) != qle
926  ) {
927  continue;
928  }
929 
930  RASTER_DEBUGF(5, "deleting index: %d => %f", i, qle->value);
931  qll->index[i].index = -1;
932  qll->index[i].element = NULL;
933  }
934 }
935 
936 static struct quantile_llist_element *quantile_llist_index_search(
937  struct quantile_llist *qll,
938  double value,
939  uint32_t *index
940 ) {
941  uint32_t i = 0, j = 0;
942  RASTER_DEBUGF(5, "searching index for value %f", value);
943 
944  for (i = 0; i < qll->index_max; i++) {
945  if (NULL == qll->index[i].element) {
946  if (i < 1) break;
947  continue;
948  }
949  if (value > (qll->index[i]).element->value) continue;
950 
951  if (FLT_EQ(value, qll->index[i].element->value)) {
952  RASTER_DEBUGF(5, "using index value at %d = %f", i, qll->index[i].element->value);
953  *index = i * 100;
954  return qll->index[i].element;
955  }
956  else if (i > 0) {
957  for (j = 1; j < i; j++) {
958  if (NULL != qll->index[i - j].element) {
959  RASTER_DEBUGF(5, "using index value at %d = %f", i - j, qll->index[i - j].element->value);
960  *index = (i - j) * 100;
961  return qll->index[i - j].element;
962  }
963  }
964  }
965  }
966 
967  *index = 0;
968  return qll->head;
969 }
970 
971 static void quantile_llist_index_reset(struct quantile_llist *qll) {
972  uint32_t i = 0;
973 
974  RASTER_DEBUG(5, "resetting index");
975  for (i = 0; i < qll->index_max; i++) {
976  qll->index[i].index = -1;
977  qll->index[i].element = NULL;
978  }
979 }
980 
981 
1009 rt_quantile
1010 rt_band_get_quantiles_stream(
1011  rt_band band,
1012  int exclude_nodata_value, double sample,
1013  uint64_t cov_count,
1014  struct quantile_llist **qlls, uint32_t *qlls_count,
1015  double *quantiles, uint32_t quantiles_count,
1016  uint32_t *rtn_count
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 }
1588 
1589 /******************************************************************************
1590 * rt_band_get_value_count()
1591 ******************************************************************************/
1592 
1607 rt_valuecount
1608 rt_band_get_value_count(
1609  rt_band band, int exclude_nodata_value,
1610  double *search_values, uint32_t search_values_count, double roundto,
1611  uint32_t *rtn_total, uint32_t *rtn_count
1612 ) {
1613  rt_valuecount vcnts = NULL;
1614  rt_pixtype pixtype = PT_END;
1615  uint8_t *data = NULL;
1616  double nodata = 0;
1617 
1618  int scale = 0;
1619  int doround = 0;
1620  double tmpd = 0;
1621  uint32_t i = 0;
1622 
1623  uint32_t x = 0;
1624  uint32_t y = 0;
1625  int rtn;
1626  double pxlval;
1627  int isnodata = 0;
1628  double rpxlval;
1629  uint32_t total = 0;
1630  uint32_t vcnts_count = 0;
1631  uint32_t new_valuecount = 0;
1632 
1633 #if POSTGIS_DEBUG_LEVEL > 0
1634  clock_t start, stop;
1635  double elapsed = 0;
1636 #endif
1637 
1638  RASTER_DEBUG(3, "starting");
1639 #if POSTGIS_DEBUG_LEVEL > 0
1640  start = clock();
1641 #endif
1642 
1643  assert(NULL != band);
1644  assert(NULL != rtn_count);
1645 
1646  data = rt_band_get_data(band);
1647  if (data == NULL) {
1648  rterror("rt_band_get_summary_stats: Cannot get band data");
1649  return NULL;
1650  }
1651 
1652  pixtype = band->pixtype;
1653 
1654  if (rt_band_get_hasnodata_flag(band)) {
1655  rt_band_get_nodata(band, &nodata);
1656  RASTER_DEBUGF(3, "hasnodata, nodataval = 1, %f", nodata);
1657  }
1658  else {
1659  exclude_nodata_value = 0;
1660  RASTER_DEBUG(3, "hasnodata, nodataval = 0, 0");
1661  }
1662 
1663  RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
1664 
1665  /* process roundto */
1666  if (roundto < 0 || FLT_EQ(roundto, 0.0)) {
1667  roundto = 0;
1668  scale = 0;
1669  }
1670  /* tenths, hundredths, thousandths, etc */
1671  else if (roundto < 1) {
1672  switch (pixtype) {
1673  /* integer band types don't have digits after the decimal place */
1674  case PT_1BB:
1675  case PT_2BUI:
1676  case PT_4BUI:
1677  case PT_8BSI:
1678  case PT_8BUI:
1679  case PT_16BSI:
1680  case PT_16BUI:
1681  case PT_32BSI:
1682  case PT_32BUI:
1683  roundto = 0;
1684  break;
1685  /* floating points, check the rounding */
1686  case PT_32BF:
1687  case PT_64BF:
1688  for (scale = 0; scale <= 20; scale++) {
1689  tmpd = roundto * pow(10, scale);
1690  if (FLT_EQ((tmpd - ((int) tmpd)), 0.0)) break;
1691  }
1692  break;
1693  case PT_END:
1694  break;
1695  }
1696  }
1697  /* ones, tens, hundreds, etc */
1698  else {
1699  for (scale = 0; scale >= -20; scale--) {
1700  tmpd = roundto * pow(10, scale);
1701  if (tmpd < 1 || FLT_EQ(tmpd, 1.0)) {
1702  if (scale == 0) doround = 1;
1703  break;
1704  }
1705  }
1706  }
1707 
1708  if (scale != 0 || doround)
1709  doround = 1;
1710  else
1711  doround = 0;
1712  RASTER_DEBUGF(3, "scale = %d", scale);
1713  RASTER_DEBUGF(3, "doround = %d", doround);
1714 
1715  /* process search_values */
1716  if (search_values_count > 0 && NULL != search_values) {
1717  vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t) * search_values_count);
1718  if (NULL == vcnts) {
1719  rterror("rt_band_get_count_of_values: Could not allocate memory for value counts");
1720  *rtn_count = 0;
1721  return NULL;
1722  }
1723 
1724  for (i = 0; i < search_values_count; i++) {
1725  vcnts[i].count = 0;
1726  vcnts[i].percent = 0;
1727  if (!doround)
1728  vcnts[i].value = search_values[i];
1729  else
1730  vcnts[i].value = ROUND(search_values[i], scale);
1731  }
1732  vcnts_count = i;
1733  }
1734  else
1735  search_values_count = 0;
1736  RASTER_DEBUGF(3, "search_values_count = %d", search_values_count);
1737 
1738  /* entire band is nodata */
1739  if (rt_band_get_isnodata_flag(band) != FALSE) {
1740  if (exclude_nodata_value) {
1741  rtwarn("All pixels of band have the NODATA value");
1742  return NULL;
1743  }
1744  else {
1745  if (search_values_count > 0) {
1746  /* check for nodata match */
1747  for (i = 0; i < search_values_count; i++) {
1748  if (!doround)
1749  tmpd = nodata;
1750  else
1751  tmpd = ROUND(nodata, scale);
1752 
1753  if (FLT_NEQ(tmpd, vcnts[i].value))
1754  continue;
1755 
1756  vcnts[i].count = band->width * band->height;
1757  if (NULL != rtn_total) *rtn_total = vcnts[i].count;
1758  vcnts->percent = 1.0;
1759  }
1760 
1761  *rtn_count = vcnts_count;
1762  }
1763  /* no defined search values */
1764  else {
1765  vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t));
1766  if (NULL == vcnts) {
1767  rterror("rt_band_get_count_of_values: Could not allocate memory for value counts");
1768  *rtn_count = 0;
1769  return NULL;
1770  }
1771 
1772  vcnts->value = nodata;
1773  vcnts->count = band->width * band->height;
1774  if (NULL != rtn_total) *rtn_total = vcnts[i].count;
1775  vcnts->percent = 1.0;
1776 
1777  *rtn_count = 1;
1778  }
1779 
1780  return vcnts;
1781  }
1782  }
1783 
1784  for (x = 0; x < band->width; x++) {
1785  for (y = 0; y < band->height; y++) {
1786  rtn = rt_band_get_pixel(band, x, y, &pxlval, &isnodata);
1787 
1788  /* error getting value, continue */
1789  if (rtn != ES_NONE)
1790  continue;
1791 
1792  if (!exclude_nodata_value || (exclude_nodata_value && !isnodata)) {
1793  total++;
1794  if (doround) {
1795  rpxlval = ROUND(pxlval, scale);
1796  }
1797  else
1798  rpxlval = pxlval;
1799  RASTER_DEBUGF(5, "(pxlval, rpxlval) => (%0.6f, %0.6f)", pxlval, rpxlval);
1800 
1801  new_valuecount = 1;
1802  /* search for match in existing valuecounts */
1803  for (i = 0; i < vcnts_count; i++) {
1804  /* match found */
1805  if (FLT_EQ(vcnts[i].value, rpxlval)) {
1806  vcnts[i].count++;
1807  new_valuecount = 0;
1808  RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
1809  break;
1810  }
1811  }
1812 
1813  /*
1814  don't add new valuecount either because
1815  - no need for new one
1816  - user-defined search values
1817  */
1818  if (!new_valuecount || search_values_count > 0) continue;
1819 
1820  /* add new valuecount */
1821  vcnts = rtrealloc(vcnts, sizeof(struct rt_valuecount_t) * (vcnts_count + 1));
1822  if (NULL == vcnts) {
1823  rterror("rt_band_get_count_of_values: Could not allocate memory for value counts");
1824  *rtn_count = 0;
1825  return NULL;
1826  }
1827 
1828  vcnts[vcnts_count].value = rpxlval;
1829  vcnts[vcnts_count].count = 1;
1830  vcnts[vcnts_count].percent = 0;
1831  RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[vcnts_count].value, vcnts[vcnts_count].count);
1832  vcnts_count++;
1833  }
1834  }
1835  }
1836 
1837 #if POSTGIS_DEBUG_LEVEL > 0
1838  stop = clock();
1839  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
1840  RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
1841 #endif
1842 
1843  for (i = 0; i < vcnts_count; i++) {
1844  vcnts[i].percent = (double) vcnts[i].count / total;
1845  RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
1846  }
1847 
1848  RASTER_DEBUG(3, "done");
1849  if (NULL != rtn_total) *rtn_total = total;
1850  *rtn_count = vcnts_count;
1851  return vcnts;
1852 }
r
char * r
Definition: cu_in_wkt.c:24
FALSE
#define FALSE
Definition: dbfopen.c:168
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
FLT_NEQ
#define FLT_NEQ(x, y)
Definition: librtcore.h:2233
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:674
ROUND
#define ROUND(x, y)
Definition: librtcore.h:2241
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:1221
rt_pixtype
rt_pixtype
Definition: librtcore.h:185
PT_32BUI
@ PT_32BUI
Definition: librtcore.h:194
PT_2BUI
@ PT_2BUI
Definition: librtcore.h:187
PT_32BSI
@ PT_32BSI
Definition: librtcore.h:193
PT_END
@ PT_END
Definition: librtcore.h:197
PT_4BUI
@ PT_4BUI
Definition: librtcore.h:188
PT_32BF
@ PT_32BF
Definition: librtcore.h:195
PT_1BB
@ PT_1BB
Definition: librtcore.h:186
PT_16BUI
@ PT_16BUI
Definition: librtcore.h:192
PT_8BSI
@ PT_8BSI
Definition: librtcore.h:189
PT_16BSI
@ PT_16BSI
Definition: librtcore.h:191
PT_64BF
@ PT_64BF
Definition: librtcore.h:196
PT_8BUI
@ PT_8BUI
Definition: librtcore.h:190
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:714
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2234
rtwarn
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
rt_bandstats
struct rt_bandstats_t * rt_bandstats
Definition: librtcore.h:150
ES_NONE
@ ES_NONE
Definition: librtcore.h:180
rtrealloc
void * rtrealloc(void *mem, size_t size)
Definition: rt_context.c:179
rt_valuecount
struct rt_valuecount_t * rt_valuecount
Definition: librtcore.h:153
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1730
rtdealloc
void rtdealloc(void *mem)
Definition: rt_context.c:186
rt_band_get_data
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition: rt_band.c:400
librtcore.h
This library is the generic raster handling section of PostGIS.
librtcore_internal.h
genraster.value
int value
Definition: genraster.py:61
genraster.count
int count
Definition: genraster.py:56
ovdump.band
band
Definition: ovdump.py:57
ovdump.data
data
Definition: ovdump.py:103
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
quantile_llist_index_delete
static void quantile_llist_index_delete(struct quantile_llist *qll, struct quantile_llist_element *qle)
Definition: rt_statistics.c:919
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
SWAP
#define SWAP(x, y)
Definition: rt_statistics.c:37
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, uint32_t quantiles_count, uint32_t *rtn_count)
Compute the default set of or requested quantiles for a coverage.
Definition: rt_statistics.c:1010
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_destroy
int quantile_llist_destroy(struct quantile_llist **list, uint32_t list_count)
Definition: rt_statistics.c:878
quicksort
static void quicksort(double *left, double *right)
Definition: rt_statistics.c:81
quantile_llist_delete
static int quantile_llist_delete(struct quantile_llist_element *element)
Definition: rt_statistics.c:855
pivot
static double pivot(double *left, double *right)
Definition: rt_statistics.c:40
partition
static double * partition(double *left, double *right, double pivot)
Definition: rt_statistics.c:66
rt_band_get_histogram
rt_histogram rt_band_get_histogram(rt_bandstats stats, uint32_t bin_count, double *bin_width, uint32_t bin_width_count, int right, double min, double max, uint32_t *rtn_count)
Count the distribution of data.
Definition: rt_statistics.c:412
rt_band_get_value_count
rt_valuecount rt_band_get_value_count(rt_band band, int exclude_nodata_value, double *search_values, uint32_t search_values_count, double roundto, uint32_t *rtn_total, uint32_t *rtn_count)
Count the number of times provided value(s) occur in the band.
Definition: rt_statistics.c:1608
ORDER
#define ORDER(x, y)
Definition: rt_statistics.c:38
quantile_llist_search
static struct quantile_llist_element * quantile_llist_search(struct quantile_llist_element *element, double needle)
Definition: rt_statistics.c:780
quantile_llist_index_reset
static void quantile_llist_index_reset(struct quantile_llist *qll)
Definition: rt_statistics.c:971
rt_band_get_quantiles
rt_quantile rt_band_get_quantiles(rt_bandstats stats, double *quantiles, int quantiles_count, uint32_t *rtn_count)
Compute the default set of or requested quantiles for a set of data the quantile formula used is same...
Definition: rt_statistics.c:672
rt_band_get_summary_stats
rt_bandstats rt_band_get_summary_stats(rt_band band, int exclude_nodata_value, double sample, int inc_vals, uint64_t *cK, double *cM, double *cQ)
Compute summary statistics for a band.
Definition: rt_statistics.c:110
quantile_llist_element::count
uint32_t count
Definition: librtcore.h:2411
quantile_llist_element::value
double value
Definition: librtcore.h:2410
quantile_llist_element::prev
struct quantile_llist_element * prev
Definition: librtcore.h:2413
quantile_llist_element::next
struct quantile_llist_element * next
Definition: librtcore.h:2414
quantile_llist_element
Definition: librtcore.h:2409
quantile_llist_index::index
uint32_t index
Definition: librtcore.h:2419
quantile_llist_index::element
struct quantile_llist_element * element
Definition: librtcore.h:2418
quantile_llist_index
Definition: librtcore.h:2417
quantile_llist::count
uint32_t count
Definition: librtcore.h:2399
quantile_llist::sum1
uint64_t sum1
Definition: librtcore.h:2405
quantile_llist::tau
uint64_t tau
Definition: librtcore.h:2395
quantile_llist::sum2
uint64_t sum2
Definition: librtcore.h:2406
quantile_llist::index_max
uint32_t index_max
Definition: librtcore.h:2403
quantile_llist::head
struct quantile_llist_element * head
Definition: librtcore.h:2397
quantile_llist::quantile
double quantile
Definition: librtcore.h:2394
quantile_llist::algeq
uint8_t algeq
Definition: librtcore.h:2393
quantile_llist::tail
struct quantile_llist_element * tail
Definition: librtcore.h:2398
quantile_llist::index
struct quantile_llist_index * index
Definition: librtcore.h:2402
quantile_llist
Definition: librtcore.h:2392
rt_band_t
Definition: librtcore.h:2313
rt_bandstats_t::mean
double mean
Definition: librtcore.h:2365
rt_bandstats_t::sorted
int sorted
Definition: librtcore.h:2369
rt_bandstats_t::sample
double sample
Definition: librtcore.h:2359
rt_bandstats_t::sum
double sum
Definition: librtcore.h:2364
rt_bandstats_t::max
double max
Definition: librtcore.h:2363
rt_bandstats_t::count
uint32_t count
Definition: librtcore.h:2360
rt_bandstats_t::min
double min
Definition: librtcore.h:2362
rt_bandstats_t::values
double * values
Definition: librtcore.h:2368
rt_bandstats_t::stddev
double stddev
Definition: librtcore.h:2366
rt_bandstats_t
Definition: librtcore.h:2358
rt_histogram_t::count
uint32_t count
Definition: librtcore.h:2374
rt_histogram_t::percent
double percent
Definition: librtcore.h:2375
rt_histogram_t::min
double min
Definition: librtcore.h:2377
rt_histogram_t::inc_min
int inc_min
Definition: librtcore.h:2380
rt_histogram_t::max
double max
Definition: librtcore.h:2378
rt_histogram_t::inc_max
int inc_max
Definition: librtcore.h:2381
rt_histogram_t
Definition: librtcore.h:2373
rt_quantile_t::value
double value
Definition: librtcore.h:2387
rt_quantile_t::quantile
double quantile
Definition: librtcore.h:2386
rt_quantile_t::has_value
uint32_t has_value
Definition: librtcore.h:2388
rt_quantile_t
Definition: librtcore.h:2385
rt_valuecount_t::count
uint32_t count
Definition: librtcore.h:2425
rt_valuecount_t::percent
double percent
Definition: librtcore.h:2426
rt_valuecount_t::value
double value
Definition: librtcore.h:2424
rt_valuecount_t
Definition: librtcore.h:2423
uint32_t
unsigned int uint32_t
Definition: uthash.h:78
uint8_t
unsigned char uint8_t
Definition: uthash.h:79