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
pivot(double * left,double * right)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
partition(double * left,double * right,double pivot)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
quicksort(double * left,double * right)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
96 /**
97 * Compute summary statistics for a band
98 *
99 * @param band : the band to query for summary stats
100 * @param exclude_nodata_value : if non-zero, ignore nodata values
101 * @param sample : percentage of pixels to sample
102 * @param inc_vals : flag to include values in return struct
103 * @param cK : number of pixels counted thus far in coverage
104 * @param cM : M component of 1-pass stddev for coverage
105 * @param cQ : Q component of 1-pass stddev for coverage
106 *
107 * @return the summary statistics for a band or NULL
108 */
109 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)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 int64_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
393 /**
394 * Count the distribution of data
395 *
396 * @param stats : a populated stats struct for processing
397 * @param bin_count : the number of bins to group the data by
398 * @param bin_width : the width of each bin as an array
399 * @param bin_width_count : number of values in bin_width
400 * @param right : evaluate bins by (a,b] rather than default [a,b)
401 * @param min : user-defined minimum value of the histogram
402 * a value less than the minimum value is not counted in any bins
403 * if min = max, min and max are not used
404 * @param max : user-defined maximum value of the histogram
405 * a value greater than the max value is not counted in any bins
406 * if min = max, min and max are not used
407 * @param rtn_count : set to the number of bins being returned
408 *
409 * @return the histogram of the data or NULL
410 */
411 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)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
660 /**
661 * Compute the default set of or requested quantiles for a set of data
662 * the quantile formula used is same as Excel and R default method
663 *
664 * @param stats : a populated stats struct for processing
665 * @param quantiles : the quantiles to be computed
666 * @param quantiles_count : the number of quantiles to be computed
667 * @param rtn_count : set to the number of quantiles being returned
668 *
669 * @return the default set of or requested quantiles for a band or NULL
670 */
671 rt_quantile
rt_band_get_quantiles(rt_bandstats stats,double * quantiles,int quantiles_count,uint32_t * rtn_count)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
quantile_llist_search(struct quantile_llist_element * element,double needle)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
quantile_llist_insert(struct quantile_llist_element * element,double value,uint32_t * idx)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
quantile_llist_delete(struct quantile_llist_element * element)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
quantile_llist_destroy(struct quantile_llist ** list,uint32_t list_count)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
quantile_llist_index_update(struct quantile_llist * qll,struct quantile_llist_element * qle,uint32_t idx)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
quantile_llist_index_delete(struct quantile_llist * qll,struct quantile_llist_element * qle)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 = UINT32_MAX;
932 qll->index[i].element = NULL;
933 }
934 }
935
quantile_llist_index_search(struct quantile_llist * qll,double value,uint32_t * index)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
quantile_llist_index_reset(struct quantile_llist * qll)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 = UINT32_MAX;
977 qll->index[i].element = NULL;
978 }
979 }
980
981
982 /**
983 * Compute the default set of or requested quantiles for a coverage
984 *
985 * This function is based upon the algorithm described in:
986 *
987 * A One-Pass Space-Efficient Algorithm for Finding Quantiles (1995)
988 * by Rakesh Agrawal, Arun Swami
989 * in Proc. 7th Intl. Conf. Management of Data (COMAD-95)
990 *
991 * http://www.almaden.ibm.com/cs/projects/iis/hdb/Publications/papers/comad95.pdf
992 *
993 * In the future, it may be worth exploring algorithms that don't
994 * require the size of the coverage
995 *
996 * @param band : the band to include in the quantile search
997 * @param exclude_nodata_value : if non-zero, ignore nodata values
998 * @param sample : percentage of pixels to sample
999 * @param cov_count : number of values in coverage
1000 * @param qlls : set of quantile_llist structures
1001 * @param qlls_count : the number of quantile_llist structures
1002 * @param quantiles : the quantiles to be computed
1003 * if bot qlls and quantiles provided, qlls is used
1004 * @param quantiles_count : the number of quantiles to be computed
1005 * @param rtn_count : the number of quantiles being returned
1006 *
1007 * @return the default set of or requested quantiles for a band or NULL
1008 */
1009 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)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 int64_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
1593 /**
1594 * Count the number of times provided value(s) occur in
1595 * the band
1596 *
1597 * @param band : the band to query for minimum and maximum pixel values
1598 * @param exclude_nodata_value : if non-zero, ignore nodata values
1599 * @param search_values : array of values to count
1600 * @param search_values_count : the number of search values
1601 * @param roundto : the decimal place to round the values to
1602 * @param rtn_total : the number of pixels examined in the band
1603 * @param rtn_count : the number of value counts being returned
1604 *
1605 * @return the number of times the provide value(s) occur or NULL
1606 */
1607 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)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 }
1853