PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches
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 * Copyright (C) 2025 Darafei Praliaskouski <me@komzpa.net>
14 *
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * This program is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with this program; if not, write to the Free Software Foundation,
27 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28 *
29 */
30
31#include "librtcore.h"
32#include "librtcore_internal.h"
33
34/******************************************************************************
35* quicksort
36******************************************************************************/
37
38#define SWAP(x, y) { double t; t = x; x = y; y = t; }
39#define ORDER(x, y) if (x > y) SWAP(x, y)
40
41static double pivot(double *left, double *right) {
42 double l, m, r, *p;
43
44 l = *left;
45 m = *(left + (right - left) / 2);
46 r = *right;
47
48 /* order */
49 ORDER(l, m);
50 ORDER(l, r);
51 ORDER(m, r);
52
53 /* pivot is higher of two values */
54 if (l < m) return m;
55 if (m < r) return r;
56
57 /* find pivot that isn't left */
58 for (p = left + 1; p <= right; ++p) {
59 if (*p != *left)
60 return (*p < *left) ? *left : *p;
61 }
62
63 /* all values are same */
64 return -1;
65}
66
67static double *partition(double *left, double *right, double pivot) {
68 while (left <= right) {
69 while (*left < pivot) ++left;
70 while (*right >= pivot) --right;
71
72 if (left < right) {
73 SWAP(*left, *right);
74 ++left;
75 --right;
76 }
77 }
78
79 return left;
80}
81
82static void quicksort(double *left, double *right) {
83 double p = pivot(left, right);
84 double *pos;
85
86 if (p != -1) {
87 pos = partition(left, right, p);
88 quicksort(left, pos - 1);
89 quicksort(pos, right);
90 }
91}
92
93/******************************************************************************
94* rt_band_get_summary_stats()
95******************************************************************************/
96
110rt_bandstats
111rt_band_get_summary_stats(
112 rt_band band,
113 int exclude_nodata_value, double sample, int inc_vals,
114 uint64_t *cK, double *cM, double *cQ
115) {
116 uint32_t x = 0;
117 int64_t y = 0;
118 uint32_t z = 0;
119 uint32_t offset = 0;
120 uint32_t diff = 0;
121 int rtn;
122 int hasnodata = FALSE;
123 double nodata = 0;
124 double *values = NULL;
125 double value;
126 int isnodata = 0;
127 rt_bandstats stats = NULL;
128
129 uint32_t do_sample = 0;
130 uint32_t sample_size = 0;
131 uint32_t sample_per = 0;
132 uint32_t sample_int = 0;
133 uint32_t i = 0;
134 uint32_t j = 0;
135 double sum = 0;
136 uint32_t k = 0;
137 double M = 0;
138 double Q = 0;
139
140#if POSTGIS_DEBUG_LEVEL > 0
141 clock_t start, stop;
142 double elapsed = 0;
143#endif
144
145 RASTER_DEBUG(3, "starting");
146#if POSTGIS_DEBUG_LEVEL > 0
147 start = clock();
148#endif
149
150 assert(NULL != band);
151
152 /* band is empty (width < 1 || height < 1) */
153 if (band->width < 1 || band->height < 1) {
154 stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
155 if (NULL == stats) {
156 rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
157 return NULL;
158 }
159
160 rtwarn("Band is empty as width and/or height is 0");
161
162 stats->sample = 1;
163 stats->sorted = 0;
164 stats->values = NULL;
165 stats->count = 0;
166 stats->min = stats->max = 0;
167 stats->sum = 0;
168 stats->mean = 0;
169 stats->stddev = -1;
170
171 return stats;
172 }
173
174 hasnodata = rt_band_get_hasnodata_flag(band);
175 if (hasnodata != FALSE)
176 rt_band_get_nodata(band, &nodata);
177 else
178 exclude_nodata_value = 0;
179
180 RASTER_DEBUGF(3, "nodata = %f", nodata);
181 RASTER_DEBUGF(3, "hasnodata = %d", hasnodata);
182 RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
183
184 /* entire band is nodata */
185 if (rt_band_get_isnodata_flag(band) != FALSE) {
186 stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
187 if (NULL == stats) {
188 rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
189 return NULL;
190 }
191
192 stats->sample = 1;
193 stats->sorted = 0;
194 stats->values = NULL;
195
196 if (exclude_nodata_value) {
197 rtwarn("All pixels of band have the NODATA value");
198
199 stats->count = 0;
200 stats->min = stats->max = 0;
201 stats->sum = 0;
202 stats->mean = 0;
203 stats->stddev = -1;
204 }
205 else {
206 stats->count = band->width * band->height;
207 stats->min = stats->max = nodata;
208 stats->sum = stats->count * nodata;
209 stats->mean = nodata;
210 stats->stddev = 0;
211 }
212
213 return stats;
214 }
215
216 /* clamp percentage */
217 if (
218 (sample < 0 || FLT_EQ(sample, 0.0)) ||
219 (sample > 1 || FLT_EQ(sample, 1.0))
220 ) {
221 do_sample = 0;
222 sample = 1;
223 }
224 else
225 do_sample = 1;
226 RASTER_DEBUGF(3, "do_sample = %d", do_sample);
227
228 /* sample all pixels */
229 if (!do_sample) {
230 sample_size = band->width * band->height;
231 sample_per = band->height;
232 }
233 /*
234 randomly sample a percentage of available pixels
235 sampling method is known as
236 "systematic random sample without replacement"
237 */
238 else {
239 sample_size = round((band->width * band->height) * sample);
240 sample_per = round(sample_size / band->width);
241 if (sample_per < 1)
242 sample_per = 1;
243 sample_int = round(band->height / sample_per);
244 srand(time(NULL));
245 }
246
247 RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set"
248 , sample_size, (band->width * band->height), sample_per);
249
250 if (inc_vals) {
251 values = rtalloc(sizeof(double) * sample_size);
252 if (NULL == values) {
253 rtwarn("Could not allocate memory for values");
254 inc_vals = 0;
255 }
256 }
257
258 /* initialize stats */
259 stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
260 if (NULL == stats) {
261 rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
262 return NULL;
263 }
264 stats->sample = sample;
265 stats->count = 0;
266 stats->sum = 0;
267 stats->mean = 0;
268 stats->stddev = -1;
269 stats->min = stats->max = 0;
270 stats->values = NULL;
271 stats->sorted = 0;
272
273 for (x = 0, j = 0, k = 0; x < band->width; x++) {
274 y = -1;
275 diff = 0;
276
277 for (i = 0, z = 0; i < sample_per; i++) {
278 if (!do_sample)
279 y = i;
280 else {
281 offset = (rand() % sample_int) + 1;
282 y += diff + offset;
283 diff = sample_int - offset;
284 }
285 RASTER_DEBUGF(5, "(x, y, z) = (%u, %lld, %u)", x, y, z);
286 if (y >= band->height || z > sample_per) break;
287
288 rtn = rt_band_get_pixel(band, x, y, &value, &isnodata);
289
290 j++;
291 if (rtn == ES_NONE && (!exclude_nodata_value || (exclude_nodata_value && !isnodata))) {
292
293 /* inc_vals set, collect pixel values */
294 if (inc_vals) values[k] = value;
295
296 /* average */
297 k++;
298 sum += value;
299
300 /*
301 one-pass standard deviation
302 http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
303 */
304 if (k == 1) {
305 Q = 0;
306 M = value;
307 }
308 else {
309 Q += (((k - 1) * pow(value - M, 2)) / k);
310 M += ((value - M ) / k);
311 }
312
313 /* coverage one-pass standard deviation */
314 if (NULL != cK) {
315 (*cK)++;
316 if (*cK == 1) {
317 *cQ = 0;
318 *cM = value;
319 }
320 else {
321 *cQ += (((*cK - 1) * pow(value - *cM, 2)) / *cK);
322 *cM += ((value - *cM ) / *cK);
323 }
324 }
325
326 /* min/max */
327 if (stats->count < 1) {
328 stats->count = 1;
329 stats->min = stats->max = value;
330 }
331 else {
332 if (value < stats->min)
333 stats->min = value;
334 if (value > stats->max)
335 stats->max = value;
336 }
337
338 }
339
340 z++;
341 }
342 }
343
344 RASTER_DEBUG(3, "sampling complete");
345
346 stats->count = k;
347 if (k > 0) {
348 if (inc_vals) {
349 /* free unused memory */
350 if (sample_size != k) {
351 values = rtrealloc(values, k * sizeof(double));
352 }
353
354 stats->values = values;
355 }
356
357 stats->sum = sum;
358 stats->mean = sum / k;
359
360 /* standard deviation */
361 if (!do_sample)
362 stats->stddev = sqrt(Q / k);
363 /* sample deviation */
364 else {
365 if (k < 2)
366 stats->stddev = -1;
367 else
368 stats->stddev = sqrt(Q / (k - 1));
369 }
370 }
371 /* inc_vals thus values allocated but not used */
372 else if (inc_vals)
373 rtdealloc(values);
374
375 /* if do_sample is one */
376 if (do_sample && k < 1)
377 rtwarn("All sampled pixels of band have the NODATA value");
378
379#if POSTGIS_DEBUG_LEVEL > 0
380 stop = clock();
381 elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
382 RASTER_DEBUGF(3, "(time, count, mean, stddev, min, max) = (%0.4f, %d, %f, %f, %f, %f)",
383 elapsed, stats->count, stats->mean, stats->stddev, stats->min, stats->max);
384#endif
385
386 RASTER_DEBUG(3, "done");
387 return stats;
388}
389
390/******************************************************************************
391* rt_band_get_histogram()
392******************************************************************************/
393
412rt_histogram
413rt_band_get_histogram(
414 rt_bandstats stats,
415 uint32_t bin_count, double *bin_width, uint32_t bin_width_count,
416 int right, double min, double max,
417 uint32_t *rtn_count
418) {
419 rt_histogram bins = NULL;
420 int init_width = 0;
421 uint32_t i;
422 uint32_t j;
423 double tmp;
424 double value;
425 int sum = 0;
426 double qmin;
427 double qmax;
428
429#if POSTGIS_DEBUG_LEVEL > 0
430 clock_t start, stop;
431 double elapsed = 0;
432#endif
433
434 RASTER_DEBUG(3, "starting");
435#if POSTGIS_DEBUG_LEVEL > 0
436 start = clock();
437#endif
438
439 assert(NULL != stats);
440 assert(NULL != rtn_count);
441
442 if (stats->count < 1 || NULL == stats->values) {
443 rterror("rt_util_get_histogram: rt_bandstats object has no value");
444 return NULL;
445 }
446
447 /* bin width must be positive numbers and not zero */
448 if (NULL != bin_width && bin_width_count > 0) {
449 for (i = 0; i < bin_width_count; i++) {
450 if (bin_width[i] < 0 || FLT_EQ(bin_width[i], 0.0)) {
451 rterror("rt_util_get_histogram: bin_width element is less than or equal to zero");
452 return NULL;
453 }
454 }
455 }
456
457 /* ignore min and max parameters */
458 if (FLT_EQ(max, min)) {
459 qmin = stats->min;
460 qmax = stats->max;
461 }
462 else {
463 qmin = min;
464 qmax = max;
465 if (qmin > qmax) {
466 qmin = max;
467 qmax = min;
468 }
469 }
470
471 /* # of bins not provided */
472 if (bin_count <= 0) {
473 /*
474 determine # of bins
475 http://en.wikipedia.org/wiki/Histogram
476
477 all computed bins are assumed to have equal width
478 */
479 /* Square-root choice for stats->count < 30 */
480 if (stats->count < 30)
481 bin_count = ceil(sqrt(stats->count));
482 /* Sturges' formula for stats->count >= 30 */
483 else
484 bin_count = ceil(log2((double) stats->count) + 1.);
485
486 /* bin_width_count provided and bin_width has value */
487 if (bin_width_count > 0 && NULL != bin_width) {
488 /* user has defined something specific */
489 if (bin_width_count > bin_count)
490 bin_count = bin_width_count;
491 else if (bin_width_count > 1) {
492 tmp = 0;
493 for (i = 0; i < bin_width_count; i++) tmp += bin_width[i];
494 bin_count = ceil((qmax - qmin) / tmp) * bin_width_count;
495 }
496 else
497 bin_count = ceil((qmax - qmin) / bin_width[0]);
498 }
499 /* set bin width count to zero so that one can be calculated */
500 else {
501 bin_width_count = 0;
502 }
503 }
504
505 /* min and max the same */
506 if (FLT_EQ(qmax, qmin)) bin_count = 1;
507
508 RASTER_DEBUGF(3, "bin_count = %d", bin_count);
509
510 /* bin count = 1, all values are in one bin */
511 if (bin_count < 2) {
512 bins = rtalloc(sizeof(struct rt_histogram_t));
513 if (NULL == bins) {
514 rterror("rt_util_get_histogram: Could not allocate memory for histogram");
515 return NULL;
516 }
517
518 bins->count = stats->count;
519 bins->percent = -1;
520 bins->min = qmin;
521 bins->max = qmax;
522 bins->inc_min = bins->inc_max = 1;
523
524 *rtn_count = bin_count;
525 return bins;
526 }
527
528 /* establish bin width */
529 if (bin_width_count == 0) {
530 bin_width_count = 1;
531
532 /* bin_width unallocated */
533 if (NULL == bin_width) {
534 bin_width = rtalloc(sizeof(double));
535 if (NULL == bin_width) {
536 rterror("rt_util_get_histogram: Could not allocate memory for bin widths");
537 return NULL;
538 }
539 init_width = 1;
540 }
541
542 bin_width[0] = (qmax - qmin) / bin_count;
543 }
544
545 /* initialize bins */
546 bins = rtalloc(bin_count * sizeof(struct rt_histogram_t));
547 if (NULL == bins) {
548 rterror("rt_util_get_histogram: Could not allocate memory for histogram");
549 if (init_width) rtdealloc(bin_width);
550 return NULL;
551 }
552 if (!right)
553 tmp = qmin;
554 else
555 tmp = qmax;
556 for (i = 0; i < bin_count;) {
557 for (j = 0; j < bin_width_count; j++) {
558 bins[i].count = 0;
559 bins->percent = -1;
560
561 if (!right) {
562 bins[i].min = tmp;
563 tmp += bin_width[j];
564 bins[i].max = tmp;
565
566 bins[i].inc_min = 1;
567 bins[i].inc_max = 0;
568 }
569 else {
570 bins[i].max = tmp;
571 tmp -= bin_width[j];
572 bins[i].min = tmp;
573
574 bins[i].inc_min = 0;
575 bins[i].inc_max = 1;
576 }
577
578 i++;
579 }
580 }
581 if (!right) {
582 bins[bin_count - 1].inc_max = 1;
583
584 /* align last bin to the max value */
585 if (bins[bin_count - 1].max < qmax)
586 bins[bin_count - 1].max = qmax;
587 }
588 else {
589 bins[bin_count - 1].inc_min = 1;
590
591 /* align first bin to the min value */
592 if (bins[bin_count - 1].min > qmin)
593 bins[bin_count - 1].min = qmin;
594 }
595
596 /* process the values */
597 for (i = 0; i < stats->count; i++) {
598 value = stats->values[i];
599
600 /* default, [a, b) */
601 if (!right) {
602 for (j = 0; j < bin_count; j++) {
603 if (
604 (!bins[j].inc_max && value < bins[j].max) || (
605 bins[j].inc_max && (
606 (value < bins[j].max) ||
607 FLT_EQ(value, bins[j].max)
608 )
609 )
610 ) {
611 bins[j].count++;
612 sum++;
613 break;
614 }
615 }
616 }
617 /* (a, b] */
618 else {
619 for (j = 0; j < bin_count; j++) {
620 if (
621 (!bins[j].inc_min && value > bins[j].min) || (
622 bins[j].inc_min && (
623 (value > bins[j].min) ||
624 FLT_EQ(value, bins[j].min)
625 )
626 )
627 ) {
628 bins[j].count++;
629 sum++;
630 break;
631 }
632 }
633 }
634 }
635
636 for (i = 0; i < bin_count; i++) {
637 bins[i].percent = ((double) bins[i].count) / sum;
638 }
639
640#if POSTGIS_DEBUG_LEVEL > 0
641 stop = clock();
642 elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
643 RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
644
645 for (j = 0; j < bin_count; j++) {
646 RASTER_DEBUGF(5, "(min, max, inc_min, inc_max, count, sum, percent) = (%f, %f, %d, %d, %d, %d, %f)",
647 bins[j].min, bins[j].max, bins[j].inc_min, bins[j].inc_max, bins[j].count, sum, bins[j].percent);
648 }
649#endif
650
651 if (init_width) rtdealloc(bin_width);
652 *rtn_count = bin_count;
653 RASTER_DEBUG(3, "done");
654 return bins;
655}
656
657/******************************************************************************
658* rt_band_get_quantiles()
659******************************************************************************/
660
672rt_quantile
673rt_band_get_quantiles(
674 rt_bandstats stats,
675 double *quantiles, int quantiles_count,
676 uint32_t *rtn_count
677) {
678 rt_quantile rtn;
679 int init_quantiles = 0;
680 int i = 0;
681 double h;
682 int hl;
683
684#if POSTGIS_DEBUG_LEVEL > 0
685 clock_t start, stop;
686 double elapsed = 0;
687#endif
688
689 RASTER_DEBUG(3, "starting");
690#if POSTGIS_DEBUG_LEVEL > 0
691 start = clock();
692#endif
693
694 assert(NULL != stats);
695 assert(NULL != rtn_count);
696
697 if (stats->count < 1 || NULL == stats->values) {
698 rterror("rt_band_get_quantiles: rt_bandstats object has no value");
699 return NULL;
700 }
701
702 /* quantiles not provided */
703 if (NULL == quantiles) {
704 /* quantile count not specified, default to quartiles */
705 if (quantiles_count < 2)
706 quantiles_count = 5;
707
708 quantiles = rtalloc(sizeof(double) * quantiles_count);
709 init_quantiles = 1;
710 if (NULL == quantiles) {
711 rterror("rt_band_get_quantiles: Could not allocate memory for quantile input");
712 return NULL;
713 }
714
715 quantiles_count--;
716 for (i = 0; i <= quantiles_count; i++)
717 quantiles[i] = ((double) i) / quantiles_count;
718 quantiles_count++;
719 }
720
721 /* check quantiles */
722 for (i = 0; i < quantiles_count; i++) {
723 if (quantiles[i] < 0. || quantiles[i] > 1.) {
724 rterror("rt_band_get_quantiles: Quantile value not between 0 and 1");
725 if (init_quantiles) rtdealloc(quantiles);
726 return NULL;
727 }
728 }
729 quicksort(quantiles, quantiles + quantiles_count - 1);
730
731 /* initialize rt_quantile */
732 rtn = rtalloc(sizeof(struct rt_quantile_t) * quantiles_count);
733 if (NULL == rtn) {
734 rterror("rt_band_get_quantiles: Could not allocate memory for quantile output");
735 if (init_quantiles) rtdealloc(quantiles);
736 return NULL;
737 }
738
739 /* sort values */
740 if (!stats->sorted) {
741 quicksort(stats->values, stats->values + stats->count - 1);
742 stats->sorted = 1;
743 }
744
745 /*
746 make quantiles
747
748 formula is that used in R (method 7) and Excel from
749 http://en.wikipedia.org/wiki/Quantile
750 */
751 for (i = 0; i < quantiles_count; i++) {
752 rtn[i].quantile = quantiles[i];
753
754 h = ((stats->count - 1.) * quantiles[i]) + 1.;
755 hl = floor(h);
756
757 /* h greater than hl, do full equation */
758 if (h > hl)
759 rtn[i].value = stats->values[hl - 1] + ((h - hl) * (stats->values[hl] - stats->values[hl - 1]));
760 /* shortcut as second part of equation is zero */
761 else
762 rtn[i].value = stats->values[hl - 1];
763 }
764
765#if POSTGIS_DEBUG_LEVEL > 0
766 stop = clock();
767 elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
768 RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
769#endif
770
771 *rtn_count = quantiles_count;
772 if (init_quantiles) rtdealloc(quantiles);
773 RASTER_DEBUG(3, "done");
774 return rtn;
775}
776
777/******************************************************************************
778* rt_band_get_quantiles_stream()
779******************************************************************************/
780
781static struct quantile_llist_element *quantile_llist_search(
782 struct quantile_llist_element *element,
783 double needle
784) {
785 if (NULL == element)
786 return NULL;
787 else if (FLT_NEQ(needle, element->value)) {
788 if (NULL != element->next)
789 return quantile_llist_search(element->next, needle);
790 else
791 return NULL;
792 }
793 else
794 return element;
795}
796
797static struct quantile_llist_element *quantile_llist_insert(
798 struct quantile_llist_element *element,
799 double value,
800 uint32_t *idx
801) {
802 struct quantile_llist_element *qle = NULL;
803
804 if (NULL == element) {
805 qle = rtalloc(sizeof(struct quantile_llist_element));
806 RASTER_DEBUGF(4, "qle @ %p is only element in list", qle);
807 if (NULL == qle) return NULL;
808
809 qle->value = value;
810 qle->count = 1;
811
812 qle->prev = NULL;
813 qle->next = NULL;
814
815 if (NULL != idx) *idx = 0;
816 return qle;
817 }
818 else if (value > element->value) {
819 if (NULL != idx) *idx += 1;
820 if (NULL != element->next)
821 return quantile_llist_insert(element->next, value, idx);
822 /* insert as last element in list */
823 else {
824 qle = rtalloc(sizeof(struct quantile_llist_element));
825 RASTER_DEBUGF(4, "insert qle @ %p as last element", qle);
826 if (NULL == qle) return NULL;
827
828 qle->value = value;
829 qle->count = 1;
830
831 qle->prev = element;
832 qle->next = NULL;
833 element->next = qle;
834
835 return qle;
836 }
837 }
838 /* insert before current element */
839 else {
840 qle = rtalloc(sizeof(struct quantile_llist_element));
841 RASTER_DEBUGF(4, "insert qle @ %p before current element", qle);
842 if (NULL == qle) return NULL;
843
844 qle->value = value;
845 qle->count = 1;
846
847 if (NULL != element->prev) element->prev->next = qle;
848 qle->next = element;
849 qle->prev = element->prev;
850 element->prev = qle;
851
852 return qle;
853 }
854}
855
856static int quantile_llist_delete(struct quantile_llist_element *element) {
857 if (NULL == element) return 0;
858
859 /* beginning of list */
860 if (NULL == element->prev && NULL != element->next) {
861 element->next->prev = NULL;
862 }
863 /* end of list */
864 else if (NULL != element->prev && NULL == element->next) {
865 element->prev->next = NULL;
866 }
867 /* within list */
868 else if (NULL != element->prev && NULL != element->next) {
869 element->prev->next = element->next;
870 element->next->prev = element->prev;
871 }
872
873 RASTER_DEBUGF(4, "qle @ %p destroyed", element);
874 rtdealloc(element);
875
876 return 1;
877}
878
879int quantile_llist_destroy(struct quantile_llist **list, uint32_t list_count) {
880 struct quantile_llist_element *element = NULL;
881 uint32_t i;
882
883 if (NULL == *list) return 0;
884
885 for (i = 0; i < list_count; i++) {
886 element = (*list)[i].head;
887 while (NULL != element->next) {
888 quantile_llist_delete(element->next);
889 }
890 quantile_llist_delete(element);
891
892 rtdealloc((*list)[i].index);
893 }
894
895 rtdealloc(*list);
896 return 1;
897}
898
899static void quantile_llist_index_update(struct quantile_llist *qll, struct quantile_llist_element *qle, uint32_t idx) {
900 uint32_t anchor = (uint32_t) floor(idx / 100);
901
902 if (qll->tail == qle) return;
903
904 if (
905 (anchor != 0) && (
906 NULL == qll->index[anchor].element ||
907 idx <= qll->index[anchor].index
908 )
909 ) {
910 qll->index[anchor].index = idx;
911 qll->index[anchor].element = qle;
912 }
913
914 if (anchor != 0 && NULL == qll->index[0].element) {
915 qll->index[0].index = 0;
916 qll->index[0].element = qll->head;
917 }
918}
919
920static void quantile_llist_index_delete(struct quantile_llist *qll, struct quantile_llist_element *qle) {
921 uint32_t i = 0;
922
923 for (i = 0; i < qll->index_max; i++) {
924 if (
925 NULL == qll->index[i].element ||
926 (qll->index[i].element) != qle
927 ) {
928 continue;
929 }
930
931 RASTER_DEBUGF(5, "deleting index: %d => %f", i, qle->value);
932 qll->index[i].index = UINT32_MAX;
933 qll->index[i].element = NULL;
934 }
935}
936
937static struct quantile_llist_element *quantile_llist_index_search(
938 struct quantile_llist *qll,
939 double value,
940 uint32_t *index
941) {
942 uint32_t i = 0, j = 0;
943 RASTER_DEBUGF(5, "searching index for value %f", value);
944
945 for (i = 0; i < qll->index_max; i++) {
946 if (NULL == qll->index[i].element) {
947 if (i < 1) break;
948 continue;
949 }
950 if (value > (qll->index[i]).element->value) continue;
951
952 if (FLT_EQ(value, qll->index[i].element->value)) {
953 RASTER_DEBUGF(5, "using index value at %d = %f", i, qll->index[i].element->value);
954 *index = i * 100;
955 return qll->index[i].element;
956 }
957 else if (i > 0) {
958 for (j = 1; j < i; j++) {
959 if (NULL != qll->index[i - j].element) {
960 RASTER_DEBUGF(5, "using index value at %d = %f", i - j, qll->index[i - j].element->value);
961 *index = (i - j) * 100;
962 return qll->index[i - j].element;
963 }
964 }
965 }
966 }
967
968 *index = 0;
969 return qll->head;
970}
971
972static void quantile_llist_index_reset(struct quantile_llist *qll) {
973 uint32_t i = 0;
974
975 RASTER_DEBUG(5, "resetting index");
976 for (i = 0; i < qll->index_max; i++) {
977 qll->index[i].index = UINT32_MAX;
978 qll->index[i].element = NULL;
979 }
980}
981
982
1010rt_quantile
1011rt_band_get_quantiles_stream(
1012 rt_band band,
1013 int exclude_nodata_value, double sample,
1014 uint64_t cov_count,
1015 struct quantile_llist **qlls, uint32_t *qlls_count,
1016 double *quantiles, uint32_t quantiles_count,
1017 uint32_t *rtn_count
1018) {
1019 rt_quantile rtn = NULL;
1020 int init_quantiles = 0;
1021
1022 struct quantile_llist *qll = NULL;
1023 struct quantile_llist_element *qle = NULL;
1024 struct quantile_llist_element *qls = NULL;
1025 const uint32_t MAX_VALUES = 750;
1026
1027 uint8_t *data = NULL;
1028 double value;
1029 int isnodata = 0;
1030
1031 uint32_t a = 0;
1032 uint32_t i = 0;
1033 uint32_t j = 0;
1034 uint32_t k = 0;
1035 uint32_t x = 0;
1036 int64_t y = 0;
1037 uint32_t z = 0;
1038 uint32_t idx = 0;
1039 uint32_t offset = 0;
1040 uint32_t diff = 0;
1041 uint8_t exists = 0;
1042
1043 uint32_t do_sample = 0;
1044 uint32_t sample_size = 0;
1045 uint32_t sample_per = 0;
1046 uint32_t sample_int = 0;
1047 int status;
1048
1049 RASTER_DEBUG(3, "starting");
1050
1051 assert(NULL != band);
1052 assert(cov_count > 1);
1053 assert(NULL != rtn_count);
1054 RASTER_DEBUGF(3, "cov_count = %llu", cov_count);
1055
1056 data = rt_band_get_data(band);
1057 if (data == NULL) {
1058 rterror("rt_band_get_summary_stats: Cannot get band data");
1059 return NULL;
1060 }
1061
1062 if (!rt_band_get_hasnodata_flag(band))
1063 exclude_nodata_value = 0;
1064 RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
1065
1066 /* quantile_llist not provided */
1067 if (NULL == *qlls) {
1068 /* quantiles not provided */
1069 if (NULL == quantiles) {
1070 /* quantile count not specified, default to quartiles */
1071 if (quantiles_count < 2)
1072 quantiles_count = 5;
1073
1074 quantiles = rtalloc(sizeof(double) * quantiles_count);
1075 init_quantiles = 1;
1076 if (NULL == quantiles) {
1077 rterror("rt_band_get_quantiles_stream: Could not allocate memory for quantile input");
1078 return NULL;
1079 }
1080
1081 quantiles_count--;
1082 for (i = 0; i <= quantiles_count; i++)
1083 quantiles[i] = ((double) i) / quantiles_count;
1084 quantiles_count++;
1085 }
1086
1087 /* check quantiles */
1088 for (i = 0; i < quantiles_count; i++) {
1089 if (quantiles[i] < 0. || quantiles[i] > 1.) {
1090 rterror("rt_band_get_quantiles_stream: Quantile value not between 0 and 1");
1091 if (init_quantiles) rtdealloc(quantiles);
1092 return NULL;
1093 }
1094 }
1095 quicksort(quantiles, quantiles + quantiles_count - 1);
1096
1097 /* initialize linked-list set */
1098 *qlls_count = quantiles_count * 2;
1099 RASTER_DEBUGF(4, "qlls_count = %d", *qlls_count);
1100 *qlls = rtalloc(sizeof(struct quantile_llist) * *qlls_count);
1101 if (NULL == *qlls) {
1102 rterror("rt_band_get_quantiles_stream: Could not allocate memory for quantile output");
1103 if (init_quantiles) rtdealloc(quantiles);
1104 return NULL;
1105 }
1106
1107 j = (uint32_t) floor(MAX_VALUES / 100.) + 1;
1108 for (i = 0; i < *qlls_count; i++) {
1109 qll = &((*qlls)[i]);
1110 qll->quantile = quantiles[(i * quantiles_count) / *qlls_count];
1111 qll->count = 0;
1112 qll->sum1 = 0;
1113 qll->sum2 = 0;
1114 qll->head = NULL;
1115 qll->tail = NULL;
1116
1117 /* initialize index */
1118 qll->index = rtalloc(sizeof(struct quantile_llist_index) * j);
1119 if (NULL == qll->index) {
1120 rterror("rt_band_get_quantiles_stream: Could not allocate memory for quantile output");
1121 if (init_quantiles) rtdealloc(quantiles);
1122 return NULL;
1123 }
1124 qll->index_max = j;
1125 quantile_llist_index_reset(qll);
1126
1127 /* AL-GEQ */
1128 if (!(i % 2)) {
1129 qll->algeq = 1;
1130 qll->tau = (uint64_t) ROUND(cov_count - (cov_count * qll->quantile), 0);
1131 if (qll->tau < 1) qll->tau = 1;
1132 }
1133 /* AL-GT */
1134 else {
1135 qll->algeq = 0;
1136 qll->tau = cov_count - (*qlls)[i - 1].tau + 1;
1137 }
1138
1139 RASTER_DEBUGF(4, "qll init: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %llu, %llu, %llu)",
1140 qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
1141 RASTER_DEBUGF(4, "qll init: (head, tail) = (%p, %p)", qll->head, qll->tail);
1142 }
1143
1144 if (init_quantiles) rtdealloc(quantiles);
1145 }
1146
1147 /* clamp percentage */
1148 if (
1149 (sample < 0 || FLT_EQ(sample, 0.0)) ||
1150 (sample > 1 || FLT_EQ(sample, 1.0))
1151 ) {
1152 do_sample = 0;
1153 sample = 1;
1154 }
1155 else
1156 do_sample = 1;
1157 RASTER_DEBUGF(3, "do_sample = %d", do_sample);
1158
1159 /* sample all pixels */
1160 if (!do_sample) {
1161 sample_size = band->width * band->height;
1162 sample_per = band->height;
1163 }
1164 /*
1165 randomly sample a percentage of available pixels
1166 sampling method is known as
1167 "systematic random sample without replacement"
1168 */
1169 else {
1170 sample_size = round((band->width * band->height) * sample);
1171 sample_per = round(sample_size / band->width);
1172 sample_int = round(band->height / sample_per);
1173 srand(time(NULL));
1174 }
1175 RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set"
1176 , sample_size, (band->width * band->height), sample_per);
1177
1178 for (x = 0, j = 0, k = 0; x < band->width; x++) {
1179 y = -1;
1180 diff = 0;
1181
1182 /* exclude_nodata_value = TRUE and band is NODATA */
1183 if (exclude_nodata_value && rt_band_get_isnodata_flag(band)) {
1184 RASTER_DEBUG(3, "Skipping quantile calculation as band is NODATA");
1185 break;
1186 }
1187
1188 for (i = 0, z = 0; i < sample_per; i++) {
1189 if (do_sample != 1)
1190 y = i;
1191 else {
1192 offset = (rand() % sample_int) + 1;
1193 y += diff + offset;
1194 diff = sample_int - offset;
1195 }
1196 RASTER_DEBUGF(5, "(x, y, z) = (%u, %lld, %u)", x, y, z);
1197 if (y >= band->height || z > sample_per) break;
1198
1199 status = rt_band_get_pixel(band, x, y, &value, &isnodata);
1200
1201 j++;
1202 if (status == ES_NONE && (!exclude_nodata_value || (exclude_nodata_value && !isnodata))) {
1203
1204 /* process each quantile */
1205 for (a = 0; a < *qlls_count; a++) {
1206 qll = &((*qlls)[a]);
1207 qls = NULL;
1208 RASTER_DEBUGF(4, "%d of %d (%f)", a + 1, *qlls_count, qll->quantile);
1209 RASTER_DEBUGF(5, "qll before: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %lld, %lld, %lld)",
1210 qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
1211 RASTER_DEBUGF(5, "qll before: (head, tail) = (%p, %p)", qll->head, qll->tail);
1212
1213 /* OPTIMIZATION: shortcuts for quantiles of zero or one */
1214 if (FLT_EQ(qll->quantile, 0.)) {
1215 if (NULL != qll->head) {
1216 if (value < qll->head->value)
1217 qll->head->value = value;
1218 }
1219 else {
1220 qle = quantile_llist_insert(qll->head, value, NULL);
1221 qll->head = qle;
1222 qll->tail = qle;
1223 qll->count = 1;
1224 }
1225
1226 RASTER_DEBUGF(4, "quantile shortcut for %f\n\n", qll->quantile);
1227 continue;
1228 }
1229 else if (FLT_EQ(qll->quantile, 1.)) {
1230 if (NULL != qll->head) {
1231 if (value > qll->head->value)
1232 qll->head->value = value;
1233 }
1234 else {
1235 qle = quantile_llist_insert(qll->head, value, NULL);
1236 qll->head = qle;
1237 qll->tail = qle;
1238 qll->count = 1;
1239 }
1240
1241 RASTER_DEBUGF(4, "quantile shortcut for %f\n\n", qll->quantile);
1242 continue;
1243 }
1244
1245 /* value exists in list */
1246 /* OPTIMIZATION: check to see if value exceeds last value */
1247 if (NULL != qll->tail && value > qll->tail->value)
1248 qle = NULL;
1249 /* OPTIMIZATION: check to see if value equals last value */
1250 else if (NULL != qll->tail && FLT_EQ(value, qll->tail->value))
1251 qle = qll->tail;
1252 /* OPTIMIZATION: use index if possible */
1253 else {
1254 qls = quantile_llist_index_search(qll, value, &idx);
1255 qle = quantile_llist_search(qls, value);
1256 }
1257
1258 /* value found */
1259 if (NULL != qle) {
1260 RASTER_DEBUGF(4, "%f found in list", value);
1261 RASTER_DEBUGF(5, "qle before: (value, count) = (%f, %d)", qle->value, qle->count);
1262
1263 qle->count++;
1264 qll->sum1++;
1265
1266 if (qll->algeq)
1267 qll->sum2 = qll->sum1 - qll->head->count;
1268 else
1269 qll->sum2 = qll->sum1 - qll->tail->count;
1270
1271 RASTER_DEBUGF(4, "qle after: (value, count) = (%f, %d)", qle->value, qle->count);
1272 }
1273 /* can still add element */
1274 else if (qll->count < MAX_VALUES) {
1275 RASTER_DEBUGF(4, "Adding %f to list", value);
1276
1277 /* insert */
1278 /* OPTIMIZATION: check to see if value exceeds last value */
1279 if (NULL != qll->tail && (value > qll->tail->value || FLT_EQ(value, qll->tail->value))) {
1280 idx = qll->count;
1281 qle = quantile_llist_insert(qll->tail, value, &idx);
1282 }
1283 /* OPTIMIZATION: use index if possible */
1284 else
1285 qle = quantile_llist_insert(qls, value, &idx);
1286 if (NULL == qle) return NULL;
1287 RASTER_DEBUGF(5, "value added at index: %d => %f", idx, value);
1288 qll->count++;
1289 qll->sum1++;
1290
1291 /* first element */
1292 if (NULL == qle->prev)
1293 qll->head = qle;
1294 /* last element */
1295 if (NULL == qle->next)
1296 qll->tail = qle;
1297
1298 if (qll->algeq)
1299 qll->sum2 = qll->sum1 - qll->head->count;
1300 else
1301 qll->sum2 = qll->sum1 - qll->tail->count;
1302
1303 /* index is only needed if there are at least 100 values */
1304 quantile_llist_index_update(qll, qle, idx);
1305
1306 RASTER_DEBUGF(5, "qle, prev, next, head, tail = %p, %p, %p, %p, %p", qle, qle->prev, qle->next, qll->head, qll->tail);
1307 }
1308 /* AL-GEQ */
1309 else if (qll->algeq) {
1310 RASTER_DEBUGF(4, "value, head->value = %f, %f", value, qll->head->value);
1311
1312 if (value < qll->head->value) {
1313 /* ignore value if test isn't true */
1314 if (qll->sum1 >= qll->tau) {
1315 RASTER_DEBUGF(4, "skipping %f", value);
1316 }
1317 else {
1318
1319 /* delete last element */
1320 RASTER_DEBUGF(4, "deleting %f from list", qll->tail->value);
1321 qle = qll->tail->prev;
1322 RASTER_DEBUGF(5, "to-be tail is %f with count %d", qle->value, qle->count);
1323 qle->count += qll->tail->count;
1324 quantile_llist_index_delete(qll, qll->tail);
1325 quantile_llist_delete(qll->tail);
1326 qll->tail = qle;
1327 qll->count--;
1328 RASTER_DEBUGF(5, "tail is %f with count %d", qll->tail->value, qll->tail->count);
1329
1330 /* insert value */
1331 RASTER_DEBUGF(4, "adding %f to list", value);
1332 /* OPTIMIZATION: check to see if value exceeds last value */
1333 if (NULL != qll->tail && (value > qll->tail->value || FLT_EQ(value, qll->tail->value))) {
1334 idx = qll->count;
1335 qle = quantile_llist_insert(qll->tail, value, &idx);
1336 }
1337 /* OPTIMIZATION: use index if possible */
1338 else {
1339 qls = quantile_llist_index_search(qll, value, &idx);
1340 qle = quantile_llist_insert(qls, value, &idx);
1341 }
1342 if (NULL == qle) return NULL;
1343 RASTER_DEBUGF(5, "value added at index: %d => %f", idx, value);
1344 qll->count++;
1345 qll->sum1++;
1346
1347 /* first element */
1348 if (NULL == qle->prev)
1349 qll->head = qle;
1350 /* last element */
1351 if (NULL == qle->next)
1352 qll->tail = qle;
1353
1354 qll->sum2 = qll->sum1 - qll->head->count;
1355
1356 quantile_llist_index_update(qll, qle, idx);
1357
1358 RASTER_DEBUGF(5, "qle, head, tail = %p, %p, %p", qle, qll->head, qll->tail);
1359
1360 }
1361 }
1362 else {
1363 qle = qll->tail;
1364 while (NULL != qle) {
1365 if (qle->value < value) {
1366 qle->count++;
1367 qll->sum1++;
1368 qll->sum2 = qll->sum1 - qll->head->count;
1369 RASTER_DEBUGF(4, "incremented count of %f by 1 to %d", qle->value, qle->count);
1370 break;
1371 }
1372
1373 qle = qle->prev;
1374 }
1375 }
1376 }
1377 /* AL-GT */
1378 else {
1379 RASTER_DEBUGF(4, "value, tail->value = %f, %f", value, qll->tail->value);
1380
1381 if (value > qll->tail->value) {
1382 /* ignore value if test isn't true */
1383 if (qll->sum1 >= qll->tau) {
1384 RASTER_DEBUGF(4, "skipping %f", value);
1385 }
1386 else {
1387
1388 /* delete last element */
1389 RASTER_DEBUGF(4, "deleting %f from list", qll->head->value);
1390 qle = qll->head->next;
1391 RASTER_DEBUGF(5, "to-be tail is %f with count %d", qle->value, qle->count);
1392 qle->count += qll->head->count;
1393 quantile_llist_index_delete(qll, qll->head);
1394 quantile_llist_delete(qll->head);
1395 qll->head = qle;
1396 qll->count--;
1397 quantile_llist_index_update(qll, NULL, 0);
1398 RASTER_DEBUGF(5, "tail is %f with count %d", qll->head->value, qll->head->count);
1399
1400 /* insert value */
1401 RASTER_DEBUGF(4, "adding %f to list", value);
1402 /* OPTIMIZATION: check to see if value exceeds last value */
1403 if (NULL != qll->tail && (value > qll->tail->value || FLT_EQ(value, qll->tail->value))) {
1404 idx = qll->count;
1405 qle = quantile_llist_insert(qll->tail, value, &idx);
1406 }
1407 /* OPTIMIZATION: use index if possible */
1408 else {
1409 qls = quantile_llist_index_search(qll, value, &idx);
1410 qle = quantile_llist_insert(qls, value, &idx);
1411 }
1412 if (NULL == qle) return NULL;
1413 RASTER_DEBUGF(5, "value added at index: %d => %f", idx, value);
1414 qll->count++;
1415 qll->sum1++;
1416
1417 /* first element */
1418 if (NULL == qle->prev)
1419 qll->head = qle;
1420 /* last element */
1421 if (NULL == qle->next)
1422 qll->tail = qle;
1423
1424 qll->sum2 = qll->sum1 - qll->tail->count;
1425
1426 quantile_llist_index_update(qll, qle, idx);
1427
1428 RASTER_DEBUGF(5, "qle, head, tail = %p, %p, %p", qle, qll->head, qll->tail);
1429
1430 }
1431 }
1432 else {
1433 qle = qll->head;
1434 while (NULL != qle) {
1435 if (qle->value > value) {
1436 qle->count++;
1437 qll->sum1++;
1438 qll->sum2 = qll->sum1 - qll->tail->count;
1439 RASTER_DEBUGF(4, "incremented count of %f by 1 to %d", qle->value, qle->count);
1440 break;
1441 }
1442
1443 qle = qle->next;
1444 }
1445 }
1446 }
1447
1448 RASTER_DEBUGF(5, "sum2, tau = %lld, %lld", qll->sum2, qll->tau);
1449 if (qll->sum2 >= qll->tau) {
1450 /* AL-GEQ */
1451 if (qll->algeq) {
1452 RASTER_DEBUGF(4, "deleting first element %f from list", qll->head->value);
1453
1454 if (NULL != qll->head->next) {
1455 qle = qll->head->next;
1456 qll->sum1 -= qll->head->count;
1457 qll->sum2 = qll->sum1 - qle->count;
1458 quantile_llist_index_delete(qll, qll->head);
1459 quantile_llist_delete(qll->head);
1460 qll->head = qle;
1461 qll->count--;
1462
1463 quantile_llist_index_update(qll, NULL, 0);
1464 }
1465 else {
1466 quantile_llist_delete(qll->head);
1467 qll->head = NULL;
1468 qll->tail = NULL;
1469 qll->sum1 = 0;
1470 qll->sum2 = 0;
1471 qll->count = 0;
1472
1473 quantile_llist_index_reset(qll);
1474 }
1475 }
1476 /* AL-GT */
1477 else {
1478 RASTER_DEBUGF(4, "deleting first element %f from list", qll->tail->value);
1479
1480 if (NULL != qll->tail->prev) {
1481 qle = qll->tail->prev;
1482 qll->sum1 -= qll->tail->count;
1483 qll->sum2 = qll->sum1 - qle->count;
1484 quantile_llist_index_delete(qll, qll->tail);
1485 quantile_llist_delete(qll->tail);
1486 qll->tail = qle;
1487 qll->count--;
1488 }
1489 else {
1490 quantile_llist_delete(qll->tail);
1491 qll->head = NULL;
1492 qll->tail = NULL;
1493 qll->sum1 = 0;
1494 qll->sum2 = 0;
1495 qll->count = 0;
1496
1497 quantile_llist_index_reset(qll);
1498 }
1499 }
1500 }
1501
1502 RASTER_DEBUGF(5, "qll after: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %lld, %lld, %lld)",
1503 qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
1504 RASTER_DEBUGF(5, "qll after: (head, tail) = (%p, %p)\n\n", qll->head, qll->tail);
1505 }
1506
1507 }
1508 else {
1509 RASTER_DEBUGF(5, "skipping value at (x, y) = (%u, %lld)", x, y);
1510 }
1511
1512 z++;
1513 }
1514 }
1515
1516 /* process quantiles */
1517 *rtn_count = *qlls_count / 2;
1518 rtn = rtalloc(sizeof(struct rt_quantile_t) * *rtn_count);
1519 if (NULL == rtn) return NULL;
1520
1521 RASTER_DEBUGF(3, "returning %d quantiles", *rtn_count);
1522 for (i = 0, k = 0; i < *qlls_count; i++) {
1523 qll = &((*qlls)[i]);
1524
1525 exists = 0;
1526 for (x = 0; x < k; x++) {
1527 if (FLT_EQ(qll->quantile, rtn[x].quantile)) {
1528 exists = 1;
1529 break;
1530 }
1531 }
1532 if (exists) continue;
1533
1534 RASTER_DEBUGF(5, "qll: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %lld, %lld, %lld)",
1535 qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
1536 RASTER_DEBUGF(5, "qll: (head, tail) = (%p, %p)", qll->head, qll->tail);
1537
1538 rtn[k].quantile = qll->quantile;
1539 rtn[k].has_value = 0;
1540
1541 /* check that qll->head and qll->tail have value */
1542 if (qll->head == NULL || qll->tail == NULL)
1543 continue;
1544
1545 /* AL-GEQ */
1546 if (qll->algeq)
1547 qle = qll->head;
1548 /* AM-GT */
1549 else
1550 qle = qll->tail;
1551
1552 exists = 0;
1553 for (j = i + 1; j < *qlls_count; j++) {
1554 if (FLT_EQ((*qlls)[j].quantile, qll->quantile)) {
1555
1556 RASTER_DEBUGF(5, "qlls[%d]: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %lld, %lld, %lld)",
1557 j, (*qlls)[j].algeq, (*qlls)[j].quantile, (*qlls)[j].count, (*qlls)[j].tau, (*qlls)[j].sum1, (*qlls)[j].sum2);
1558 RASTER_DEBUGF(5, "qlls[%d]: (head, tail) = (%p, %p)", j, (*qlls)[j].head, (*qlls)[j].tail);
1559
1560 exists = 1;
1561 break;
1562 }
1563 }
1564
1565 /* weighted average for quantile */
1566 if (exists) {
1567 if ((*qlls)[j].algeq) {
1568 rtn[k].value = ((qle->value * qle->count) + ((*qlls)[j].head->value * (*qlls)[j].head->count)) / (qle->count + (*qlls)[j].head->count);
1569 RASTER_DEBUGF(5, "qlls[%d].head: (value, count) = (%f, %d)", j, (*qlls)[j].head->value, (*qlls)[j].head->count);
1570 }
1571 else {
1572 rtn[k].value = ((qle->value * qle->count) + ((*qlls)[j].tail->value * (*qlls)[j].tail->count)) / (qle->count + (*qlls)[j].tail->count);
1573 RASTER_DEBUGF(5, "qlls[%d].tail: (value, count) = (%f, %d)", j, (*qlls)[j].tail->value, (*qlls)[j].tail->count);
1574 }
1575 }
1576 /* straight value for quantile */
1577 else {
1578 rtn[k].value = qle->value;
1579 }
1580 rtn[k].has_value = 1;
1581 RASTER_DEBUGF(3, "(quantile, value) = (%f, %f)\n\n", rtn[k].quantile, rtn[k].value);
1582
1583 k++;
1584 }
1585
1586 RASTER_DEBUG(3, "done");
1587 return rtn;
1588}
1589
1590/******************************************************************************
1591* rt_band_get_value_count()
1592******************************************************************************/
1593
1608rt_valuecount
1609rt_band_get_value_count(
1610 rt_band band, int exclude_nodata_value,
1611 double *search_values, uint32_t search_values_count, double roundto,
1612 uint32_t *rtn_total, uint32_t *rtn_count
1613) {
1614 rt_valuecount vcnts = NULL;
1615 rt_pixtype pixtype = PT_END;
1616 uint8_t *data = NULL;
1617 double nodata = 0;
1618
1619 int scale = 0;
1620 int doround = 0;
1621 double tmpd = 0;
1622 uint32_t i = 0;
1623
1624 uint32_t x = 0;
1625 uint32_t y = 0;
1626 int rtn;
1627 double pxlval;
1628 int isnodata = 0;
1629 double rpxlval;
1630 uint32_t total = 0;
1631 uint32_t vcnts_count = 0;
1632 uint32_t new_valuecount = 0;
1633
1634#if POSTGIS_DEBUG_LEVEL > 0
1635 clock_t start, stop;
1636 double elapsed = 0;
1637#endif
1638
1639 RASTER_DEBUG(3, "starting");
1640#if POSTGIS_DEBUG_LEVEL > 0
1641 start = clock();
1642#endif
1643
1644 assert(NULL != band);
1645 assert(NULL != rtn_count);
1646
1647 data = rt_band_get_data(band);
1648 if (data == NULL) {
1649 rterror("rt_band_get_summary_stats: Cannot get band data");
1650 return NULL;
1651 }
1652
1653 pixtype = band->pixtype;
1654
1655 if (rt_band_get_hasnodata_flag(band)) {
1656 rt_band_get_nodata(band, &nodata);
1657 RASTER_DEBUGF(3, "hasnodata, nodataval = 1, %f", nodata);
1658 }
1659 else {
1660 exclude_nodata_value = 0;
1661 RASTER_DEBUG(3, "hasnodata, nodataval = 0, 0");
1662 }
1663
1664 RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
1665
1666 /* process roundto */
1667 if (roundto < 0 || FLT_EQ(roundto, 0.0)) {
1668 roundto = 0;
1669 scale = 0;
1670 }
1671 /* tenths, hundredths, thousandths, etc */
1672 else if (roundto < 1) {
1673 switch (pixtype) {
1674 /* integer band types don't have digits after the decimal place */
1675 case PT_1BB:
1676 case PT_2BUI:
1677 case PT_4BUI:
1678 case PT_8BSI:
1679 case PT_8BUI:
1680 case PT_16BSI:
1681 case PT_16BUI:
1682 case PT_32BSI:
1683 case PT_32BUI:
1684 roundto = 0;
1685 break;
1686 /* floating points, check the rounding */
1687 case PT_16BF:
1688 case PT_32BF:
1689 case PT_64BF:
1690 for (scale = 0; scale <= 20; scale++)
1691 {
1692 tmpd = roundto * pow(10, scale);
1693 if (FLT_EQ((tmpd - ((int)tmpd)), 0.0))
1694 break;
1695 }
1696 break;
1697 case PT_END:
1698 break;
1699 }
1700 }
1701 /* ones, tens, hundreds, etc */
1702 else {
1703 for (scale = 0; scale >= -20; scale--) {
1704 tmpd = roundto * pow(10, scale);
1705 if (tmpd < 1 || FLT_EQ(tmpd, 1.0)) {
1706 if (scale == 0) doround = 1;
1707 break;
1708 }
1709 }
1710 }
1711
1712 if (scale != 0 || doround)
1713 doround = 1;
1714 else
1715 doround = 0;
1716 RASTER_DEBUGF(3, "scale = %d", scale);
1717 RASTER_DEBUGF(3, "doround = %d", doround);
1718
1719 /* process search_values */
1720 if (search_values_count > 0 && NULL != search_values) {
1721 vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t) * search_values_count);
1722 if (NULL == vcnts) {
1723 rterror("rt_band_get_count_of_values: Could not allocate memory for value counts");
1724 *rtn_count = 0;
1725 return NULL;
1726 }
1727
1728 for (i = 0; i < search_values_count; i++) {
1729 vcnts[i].count = 0;
1730 vcnts[i].percent = 0;
1731 if (!doround)
1732 vcnts[i].value = search_values[i];
1733 else
1734 vcnts[i].value = ROUND(search_values[i], scale);
1735 }
1736 vcnts_count = i;
1737 }
1738 else
1739 search_values_count = 0;
1740 RASTER_DEBUGF(3, "search_values_count = %d", search_values_count);
1741
1742 /* entire band is nodata */
1743 if (rt_band_get_isnodata_flag(band) != FALSE) {
1744 if (exclude_nodata_value) {
1745 rtwarn("All pixels of band have the NODATA value");
1746 return NULL;
1747 }
1748 else {
1749 if (search_values_count > 0) {
1750 /* check for nodata match */
1751 for (i = 0; i < search_values_count; i++) {
1752 if (!doround)
1753 tmpd = nodata;
1754 else
1755 tmpd = ROUND(nodata, scale);
1756
1757 if (FLT_NEQ(tmpd, vcnts[i].value))
1758 continue;
1759
1760 vcnts[i].count = band->width * band->height;
1761 if (NULL != rtn_total) *rtn_total = vcnts[i].count;
1762 vcnts->percent = 1.0;
1763 }
1764
1765 *rtn_count = vcnts_count;
1766 }
1767 /* no defined search values */
1768 else {
1769 vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t));
1770 if (NULL == vcnts) {
1771 rterror("rt_band_get_count_of_values: Could not allocate memory for value counts");
1772 *rtn_count = 0;
1773 return NULL;
1774 }
1775
1776 vcnts->value = nodata;
1777 vcnts->count = band->width * band->height;
1778 if (NULL != rtn_total) *rtn_total = vcnts[i].count;
1779 vcnts->percent = 1.0;
1780
1781 *rtn_count = 1;
1782 }
1783
1784 return vcnts;
1785 }
1786 }
1787
1788 for (x = 0; x < band->width; x++) {
1789 for (y = 0; y < band->height; y++) {
1790 rtn = rt_band_get_pixel(band, x, y, &pxlval, &isnodata);
1791
1792 /* error getting value, continue */
1793 if (rtn != ES_NONE)
1794 continue;
1795
1796 if (!exclude_nodata_value || (exclude_nodata_value && !isnodata)) {
1797 total++;
1798 if (doround) {
1799 rpxlval = ROUND(pxlval, scale);
1800 }
1801 else
1802 rpxlval = pxlval;
1803 RASTER_DEBUGF(5, "(pxlval, rpxlval) => (%0.6f, %0.6f)", pxlval, rpxlval);
1804
1805 new_valuecount = 1;
1806 /* search for match in existing valuecounts */
1807 for (i = 0; i < vcnts_count; i++) {
1808 /* match found */
1809 if (FLT_EQ(vcnts[i].value, rpxlval)) {
1810 vcnts[i].count++;
1811 new_valuecount = 0;
1812 RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
1813 break;
1814 }
1815 }
1816
1817 /*
1818 don't add new valuecount either because
1819 - no need for new one
1820 - user-defined search values
1821 */
1822 if (!new_valuecount || search_values_count > 0) continue;
1823
1824 /* add new valuecount */
1825 vcnts = rtrealloc(vcnts, sizeof(struct rt_valuecount_t) * (vcnts_count + 1));
1826 if (NULL == vcnts) {
1827 rterror("rt_band_get_count_of_values: Could not allocate memory for value counts");
1828 *rtn_count = 0;
1829 return NULL;
1830 }
1831
1832 vcnts[vcnts_count].value = rpxlval;
1833 vcnts[vcnts_count].count = 1;
1834 vcnts[vcnts_count].percent = 0;
1835 RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[vcnts_count].value, vcnts[vcnts_count].count);
1836 vcnts_count++;
1837 }
1838 }
1839 }
1840
1841#if POSTGIS_DEBUG_LEVEL > 0
1842 stop = clock();
1843 elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
1844 RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
1845#endif
1846
1847 for (i = 0; i < vcnts_count; i++) {
1848 vcnts[i].percent = (double) vcnts[i].count / total;
1849 RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
1850 }
1851
1852 RASTER_DEBUG(3, "done");
1853 if (NULL != rtn_total) *rtn_total = total;
1854 *rtn_count = vcnts_count;
1855 return vcnts;
1856}
r
char * r
Definition cu_in_wkt.c:24
FALSE
#define FALSE
Definition dbfopen.c:72
rterror
void rterror(const char *fmt,...) __attribute__((format(printf
Wrappers used for reporting errors and info.
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition rt_context.c:191
FLT_NEQ
#define FLT_NEQ(x, y)
Definition librtcore.h:2435
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition librtcore.h:304
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:308
rtwarn
void void void rtwarn(const char *fmt,...) __attribute__((format(printf
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:833
rt_band_get_data
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition rt_band.c:559
ROUND
#define ROUND(x, y)
Definition librtcore.h:2443
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:1551
rt_pixtype
rt_pixtype
Definition librtcore.h:188
PT_32BUI
@ PT_32BUI
Definition librtcore.h:197
PT_16BF
@ PT_16BF
Definition librtcore.h:198
PT_2BUI
@ PT_2BUI
Definition librtcore.h:190
PT_32BSI
@ PT_32BSI
Definition librtcore.h:196
PT_END
@ PT_END
Definition librtcore.h:201
PT_4BUI
@ PT_4BUI
Definition librtcore.h:191
PT_32BF
@ PT_32BF
Definition librtcore.h:199
PT_1BB
@ PT_1BB
Definition librtcore.h:189
PT_16BUI
@ PT_16BUI
Definition librtcore.h:195
PT_8BSI
@ PT_8BSI
Definition librtcore.h:192
PT_16BSI
@ PT_16BSI
Definition librtcore.h:194
PT_64BF
@ PT_64BF
Definition librtcore.h:200
PT_8BUI
@ PT_8BUI
Definition librtcore.h:193
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition rt_band.c:873
FLT_EQ
#define FLT_EQ(x, y)
Definition librtcore.h:2436
rt_bandstats
struct rt_bandstats_t * rt_bandstats
Definition librtcore.h:151
ES_NONE
@ ES_NONE
Definition librtcore.h:182
rtrealloc
void * rtrealloc(void *mem, size_t size)
Definition rt_context.c:199
rt_valuecount
struct rt_valuecount_t * rt_valuecount
Definition librtcore.h:154
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition rt_band.c:2067
rtdealloc
void rtdealloc(void *mem)
Definition rt_context.c:206
librtcore.h
This library is the generic raster handling section of PostGIS.
librtcore_internal.h
UINT32_MAX
#define UINT32_MAX
Definition lwin_wkt_lex.c:343
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:937
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:899
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:797
quantile_llist_index_delete
static void quantile_llist_index_delete(struct quantile_llist *qll, struct quantile_llist_element *qle)
Definition rt_statistics.c:920
quantile_llist_search
static struct quantile_llist_element * quantile_llist_search(struct quantile_llist_element *element, double needle)
Definition rt_statistics.c:781
SWAP
#define SWAP(x, y)
Definition rt_statistics.c:38
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:1011
quantile_llist_destroy
int quantile_llist_destroy(struct quantile_llist **list, uint32_t list_count)
Definition rt_statistics.c:879
quicksort
static void quicksort(double *left, double *right)
Definition rt_statistics.c:82
quantile_llist_delete
static int quantile_llist_delete(struct quantile_llist_element *element)
Definition rt_statistics.c:856
partition
static double * partition(double *left, double *right, double pivot)
Definition rt_statistics.c:67
pivot
static double pivot(double *left, double *right)
Definition rt_statistics.c:41
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:413
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:1609
ORDER
#define ORDER(x, y)
Definition rt_statistics.c:39
quantile_llist_index_reset
static void quantile_llist_index_reset(struct quantile_llist *qll)
Definition rt_statistics.c:972
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:673
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:111
quantile_llist_element::count
uint32_t count
Definition librtcore.h:2613
quantile_llist_element::value
double value
Definition librtcore.h:2612
quantile_llist_element::prev
struct quantile_llist_element * prev
Definition librtcore.h:2615
quantile_llist_element::next
struct quantile_llist_element * next
Definition librtcore.h:2616
quantile_llist_element
Definition librtcore.h:2611
quantile_llist_index::index
uint32_t index
Definition librtcore.h:2621
quantile_llist_index::element
struct quantile_llist_element * element
Definition librtcore.h:2620
quantile_llist_index
Definition librtcore.h:2619
quantile_llist::count
uint32_t count
Definition librtcore.h:2601
quantile_llist::sum1
uint64_t sum1
Definition librtcore.h:2607
quantile_llist::tau
uint64_t tau
Definition librtcore.h:2597
quantile_llist::sum2
uint64_t sum2
Definition librtcore.h:2608
quantile_llist::index_max
uint32_t index_max
Definition librtcore.h:2605
quantile_llist::head
struct quantile_llist_element * head
Definition librtcore.h:2599
quantile_llist::quantile
double quantile
Definition librtcore.h:2596
quantile_llist::algeq
uint8_t algeq
Definition librtcore.h:2595
quantile_llist::tail
struct quantile_llist_element * tail
Definition librtcore.h:2600
quantile_llist::index
struct quantile_llist_index * index
Definition librtcore.h:2604
quantile_llist
Definition librtcore.h:2594
rt_band_t
Definition librtcore.h:2515
rt_bandstats_t::mean
double mean
Definition librtcore.h:2567
rt_bandstats_t::sorted
int sorted
Definition librtcore.h:2571
rt_bandstats_t::sample
double sample
Definition librtcore.h:2561
rt_bandstats_t::sum
double sum
Definition librtcore.h:2566
rt_bandstats_t::max
double max
Definition librtcore.h:2565
rt_bandstats_t::count
uint32_t count
Definition librtcore.h:2562
rt_bandstats_t::min
double min
Definition librtcore.h:2564
rt_bandstats_t::values
double * values
Definition librtcore.h:2570
rt_bandstats_t::stddev
double stddev
Definition librtcore.h:2568
rt_bandstats_t
Definition librtcore.h:2560
rt_histogram_t::count
uint32_t count
Definition librtcore.h:2576
rt_histogram_t::percent
double percent
Definition librtcore.h:2577
rt_histogram_t::min
double min
Definition librtcore.h:2579
rt_histogram_t::inc_min
int inc_min
Definition librtcore.h:2582
rt_histogram_t::max
double max
Definition librtcore.h:2580
rt_histogram_t::inc_max
int inc_max
Definition librtcore.h:2583
rt_histogram_t
Definition librtcore.h:2575
rt_quantile_t::value
double value
Definition librtcore.h:2589
rt_quantile_t::quantile
double quantile
Definition librtcore.h:2588
rt_quantile_t::has_value
uint32_t has_value
Definition librtcore.h:2590
rt_quantile_t
Definition librtcore.h:2587
rt_valuecount_t::count
uint32_t count
Definition librtcore.h:2627
rt_valuecount_t::percent
double percent
Definition librtcore.h:2628
rt_valuecount_t::value
double value
Definition librtcore.h:2626
rt_valuecount_t
Definition librtcore.h:2625