PostGIS  2.1.10dev-r@@SVN_REVISION@@
rt_api.c
Go to the documentation of this file.
1 /*
2  * $Id: rt_api.c 14705 2016-02-26 11:37:39Z pramsey $
3  *
4  * WKTRaster - Raster Types for PostGIS
5  * http://trac.osgeo.org/postgis/wiki/WKTRaster
6  *
7  * Copyright (C) 2011-2013 Regents of the University of California
8  * <bkpark@ucdavis.edu>
9  * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
10  * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
11  * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
12  * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
13  * Copyright (C) 2008-2009 Sandro Santilli <strk@keybit.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 <math.h>
32 #include <stdio.h> /* for printf (default message handler) */
33 #include <stdarg.h> /* for va_list, va_start etc */
34 #include <string.h> /* for memcpy and strlen */
35 #include <assert.h>
36 #include <time.h> /* for time */
37 #include "rt_api.h"
38 #include "gdal_vrt.h"
39 
40 /******************************************************************************
41 * Some rules for *.(c|h) files in rt_core
42 *
43 * All functions in rt_core that receive a band index parameter
44 * must be 0-based
45 *
46 * Variables and functions internal for a public function should be prefixed
47 * with _rti_, e.g. _rti_iterator_arg.
48 ******************************************************************************/
49 
50 /*--- Utilities -------------------------------------------------*/
51 
52 static void
53 swap_char(uint8_t *a, uint8_t *b) {
54  uint8_t c = 0;
55 
56  assert(NULL != a && NULL != b);
57 
58  c = *a;
59  *a = *b;
60  *b = c;
61 }
62 
63 static void
64 flip_endian_16(uint8_t *d) {
65  assert(NULL != d);
66 
67  swap_char(d, d + 1);
68 }
69 
70 static void
71 flip_endian_32(uint8_t *d) {
72  assert(NULL != d);
73 
74  swap_char(d, d + 3);
75  swap_char(d + 1, d + 2);
76 }
77 
78 static void
79 flip_endian_64(uint8_t *d) {
80  assert(NULL != d);
81 
82  swap_char(d + 7, d);
83  swap_char(d + 6, d + 1);
84  swap_char(d + 5, d + 2);
85  swap_char(d + 4, d + 3);
86 }
87 
88 uint8_t
90  return (uint8_t)fmin(fmax((value), 0), POSTGIS_RT_1BBMAX);
91 }
92 
93 uint8_t
95  return (uint8_t)fmin(fmax((value), 0), POSTGIS_RT_2BUIMAX);
96 }
97 
98 uint8_t
100  return (uint8_t)fmin(fmax((value), 0), POSTGIS_RT_4BUIMAX);
101 }
102 
103 int8_t
105  return (int8_t)fmin(fmax((value), SCHAR_MIN), SCHAR_MAX);
106 }
107 
108 uint8_t
110  return (uint8_t)fmin(fmax((value), 0), UCHAR_MAX);
111 }
112 
113 int16_t
115  return (int16_t)fmin(fmax((value), SHRT_MIN), SHRT_MAX);
116 }
117 
118 uint16_t
120  return (uint16_t)fmin(fmax((value), 0), USHRT_MAX);
121 }
122 
123 int32_t
125  return (int32_t)fmin(fmax((value), INT_MIN), INT_MAX);
126 }
127 
128 uint32_t
130  return (uint32_t)fmin(fmax((value), 0), UINT_MAX);
131 }
132 
133 float
135  return (float)fmin(fmax((value), -FLT_MAX), FLT_MAX);
136 }
137 
138 /* quicksort */
139 #define SWAP(x, y) { double t; t = x; x = y; y = t; }
140 #define ORDER(x, y) if (x > y) SWAP(x, y)
141 
142 static double pivot(double *left, double *right) {
143  double l, m, r, *p;
144 
145  l = *left;
146  m = *(left + (right - left) / 2);
147  r = *right;
148 
149  /* order */
150  ORDER(l, m);
151  ORDER(l, r);
152  ORDER(m, r);
153 
154  /* pivot is higher of two values */
155  if (l < m) return m;
156  if (m < r) return r;
157 
158  /* find pivot that isn't left */
159  for (p = left + 1; p <= right; ++p) {
160  if (*p != *left)
161  return (*p < *left) ? *left : *p;
162  }
163 
164  /* all values are same */
165  return -1;
166 }
167 
168 static double *partition(double *left, double *right, double pivot) {
169  while (left <= right) {
170  while (*left < pivot) ++left;
171  while (*right >= pivot) --right;
172 
173  if (left < right) {
174  SWAP(*left, *right);
175  ++left;
176  --right;
177  }
178  }
179 
180  return left;
181 }
182 
183 static void quicksort(double *left, double *right) {
184  double p = pivot(left, right);
185  double *pos;
186 
187  if (p != -1) {
188  pos = partition(left, right, p);
189  quicksort(left, pos - 1);
190  quicksort(pos, right);
191  }
192 }
193 
201 GDALResampleAlg
202 rt_util_gdal_resample_alg(const char *algname) {
203  assert(algname != NULL && strlen(algname) > 0);
204 
205  if (strcmp(algname, "NEARESTNEIGHBOUR") == 0)
206  return GRA_NearestNeighbour;
207  else if (strcmp(algname, "NEARESTNEIGHBOR") == 0)
208  return GRA_NearestNeighbour;
209  else if (strcmp(algname, "BILINEAR") == 0)
210  return GRA_Bilinear;
211  else if (strcmp(algname, "CUBICSPLINE") == 0)
212  return GRA_CubicSpline;
213  else if (strcmp(algname, "CUBIC") == 0)
214  return GRA_Cubic;
215  else if (strcmp(algname, "LANCZOS") == 0)
216  return GRA_Lanczos;
217 
218  return GRA_NearestNeighbour;
219 }
220 
228 GDALDataType
230  switch (pt) {
231  case PT_1BB:
232  case PT_2BUI:
233  case PT_4BUI:
234  case PT_8BUI:
235  return GDT_Byte;
236  case PT_8BSI:
237  case PT_16BSI:
238  return GDT_Int16;
239  case PT_16BUI:
240  return GDT_UInt16;
241  case PT_32BSI:
242  return GDT_Int32;
243  case PT_32BUI:
244  return GDT_UInt32;
245  case PT_32BF:
246  return GDT_Float32;
247  case PT_64BF:
248  return GDT_Float64;
249  default:
250  return GDT_Unknown;
251  }
252 
253  return GDT_Unknown;
254 }
255 
265  switch (gdt) {
266  case GDT_Byte:
267  return PT_8BUI;
268  case GDT_UInt16:
269  return PT_16BUI;
270  case GDT_Int16:
271  return PT_16BSI;
272  case GDT_UInt32:
273  return PT_32BUI;
274  case GDT_Int32:
275  return PT_32BSI;
276  case GDT_Float32:
277  return PT_32BF;
278  case GDT_Float64:
279  return PT_64BF;
280  default:
281  return PT_END;
282  }
283 
284  return PT_END;
285 }
286 
287 /*
288  get GDAL runtime version information
289 */
290 const char*
291 rt_util_gdal_version(const char *request) {
292  if (NULL == request || !strlen(request))
293  return GDALVersionInfo("RELEASE_NAME");
294  else
295  return GDALVersionInfo(request);
296 }
297 
298 /*
299  computed extent type
300 */
302 rt_util_extent_type(const char *name) {
303  assert(name != NULL && strlen(name) > 0);
304 
305  if (strcmp(name, "UNION") == 0)
306  return ET_UNION;
307  else if (strcmp(name, "FIRST") == 0)
308  return ET_FIRST;
309  else if (strcmp(name, "SECOND") == 0)
310  return ET_SECOND;
311  else if (strcmp(name, "LAST") == 0)
312  return ET_LAST;
313  else if (strcmp(name, "CUSTOM") == 0)
314  return ET_CUSTOM;
315  else
316  return ET_INTERSECTION;
317 }
318 
319 /*
320  convert the spatial reference string from a GDAL recognized format to either WKT or Proj4
321 */
322 char*
323 rt_util_gdal_convert_sr(const char *srs, int proj4) {
324  OGRSpatialReferenceH hsrs;
325  char *rtn = NULL;
326 
327  assert(srs != NULL);
328 
329  hsrs = OSRNewSpatialReference(NULL);
330  if (OSRSetFromUserInput(hsrs, srs) == OGRERR_NONE) {
331  if (proj4)
332  OSRExportToProj4(hsrs, &rtn);
333  else
334  OSRExportToWkt(hsrs, &rtn);
335  }
336  else {
337  rterror("rt_util_gdal_convert_sr: Could not process the provided srs: %s", srs);
338  return NULL;
339  }
340 
341  OSRDestroySpatialReference(hsrs);
342  if (rtn == NULL) {
343  rterror("rt_util_gdal_convert_sr: Could not process the provided srs: %s", srs);
344  return NULL;
345  }
346 
347  return rtn;
348 }
349 
350 /*
351  is the spatial reference string supported by GDAL
352 */
353 int
354 rt_util_gdal_supported_sr(const char *srs) {
355  OGRSpatialReferenceH hsrs;
356  OGRErr rtn = OGRERR_NONE;
357 
358  assert(srs != NULL);
359 
360  hsrs = OSRNewSpatialReference(NULL);
361  rtn = OSRSetFromUserInput(hsrs, srs);
362  OSRDestroySpatialReference(hsrs);
363 
364  if (rtn == OGRERR_NONE)
365  return 1;
366  else
367  return 0;
368 }
369 
381 rt_util_gdal_sr_auth_info(GDALDatasetH hds, char **authname, char **authcode) {
382  const char *srs = NULL;
383 
384  assert(authname != NULL);
385  assert(authcode != NULL);
386 
387  *authname = NULL;
388  *authcode = NULL;
389 
390  srs = GDALGetProjectionRef(hds);
391  if (srs != NULL && srs[0] != '\0') {
392  OGRSpatialReferenceH hSRS = OSRNewSpatialReference(NULL);
393 
394  if (OSRSetFromUserInput(hSRS, srs) == OGRERR_NONE) {
395  const char* pszAuthorityName = OSRGetAuthorityName(hSRS, NULL);
396  const char* pszAuthorityCode = OSRGetAuthorityCode(hSRS, NULL);
397 
398  if (pszAuthorityName != NULL && pszAuthorityCode != NULL) {
399  *authname = rtalloc(sizeof(char) * (strlen(pszAuthorityName) + 1));
400  *authcode = rtalloc(sizeof(char) * (strlen(pszAuthorityCode) + 1));
401 
402  if (*authname == NULL || *authcode == NULL) {
403  rterror("rt_util_gdal_sr_auth_info: Could not allocate memory for auth name and code");
404  if (*authname != NULL) rtdealloc(*authname);
405  if (*authcode != NULL) rtdealloc(*authcode);
406  OSRDestroySpatialReference(hSRS);
407  return ES_ERROR;
408  }
409 
410  strncpy(*authname, pszAuthorityName, strlen(pszAuthorityName) + 1);
411  strncpy(*authcode, pszAuthorityCode, strlen(pszAuthorityCode) + 1);
412  }
413  }
414 
415  OSRDestroySpatialReference(hSRS);
416  }
417 
418  return ES_NONE;
419 }
420 
421 /*
422  is GDAL configured correctly?
423 */
425 
426  /* set of EPSG codes */
427  if (!rt_util_gdal_supported_sr("EPSG:4326"))
428  return 0;
429  if (!rt_util_gdal_supported_sr("EPSG:4269"))
430  return 0;
431  if (!rt_util_gdal_supported_sr("EPSG:4267"))
432  return 0;
433  if (!rt_util_gdal_supported_sr("EPSG:3310"))
434  return 0;
435  if (!rt_util_gdal_supported_sr("EPSG:2163"))
436  return 0;
437 
438  return 1;
439 }
440 
441 /*
442  register all GDAL drivers
443 */
444 int
445 rt_util_gdal_register_all(int force_register_all) {
446  static int registered = 0;
447 
448  if (registered && !force_register_all)
449  return 0;
450 
451  GDALAllRegister();
452  registered = 1;
453 
454  return 1;
455 }
456 
457 /*
458  is the driver registered?
459 */
460 int
462  int count = GDALGetDriverCount();
463  int i = 0;
464  GDALDriverH hdrv = NULL;
465 
466  if (drv == NULL || !strlen(drv) || count < 1)
467  return 0;
468 
469  for (i = 0; i < count; i++) {
470  hdrv = GDALGetDriver(i);
471  if (hdrv == NULL) continue;
472 
473  if (strcmp(drv, GDALGetDriverShortName(hdrv)) == 0)
474  return 1;
475  }
476 
477  return 0;
478 }
479 
480 /* variable for PostgreSQL GUC: postgis.gdal_enabled_drivers */
481 char *gdal_enabled_drivers = NULL;
482 
483 /*
484  wrapper for GDALOpen and GDALOpenShared
485 */
486 GDALDatasetH
487 rt_util_gdal_open(const char *fn, GDALAccess fn_access, int shared) {
488  assert(NULL != fn);
489 
490  if (gdal_enabled_drivers != NULL) {
491  if (strstr(gdal_enabled_drivers, GDAL_DISABLE_ALL) != NULL) {
492  rterror("rt_util_gdal_open: Cannot open file. All GDAL drivers disabled");
493  return NULL;
494  }
495  else if (strstr(gdal_enabled_drivers, GDAL_ENABLE_ALL) != NULL) {
496  /* do nothing */
497  }
498  else if (
499  (strstr(fn, "/vsicurl") != NULL) &&
500  (strstr(gdal_enabled_drivers, GDAL_VSICURL) == NULL)
501  ) {
502  rterror("rt_util_gdal_open: Cannot open VSICURL file. VSICURL disabled");
503  return NULL;
504  }
505  }
506 
507  if (shared)
508  return GDALOpenShared(fn, fn_access);
509  else
510  return GDALOpen(fn, fn_access);
511 }
512 
513 void
515  OGREnvelope env,
516  rt_envelope *ext
517 ) {
518  assert(ext != NULL);
519 
520  ext->MinX = env.MinX;
521  ext->MaxX = env.MaxX;
522  ext->MinY = env.MinY;
523  ext->MaxY = env.MaxY;
524 
525  ext->UpperLeftX = env.MinX;
526  ext->UpperLeftY = env.MaxY;
527 }
528 
529 void
531  rt_envelope ext,
532  OGREnvelope *env
533 ) {
534  assert(env != NULL);
535 
536  env->MinX = ext.MinX;
537  env->MaxX = ext.MaxX;
538  env->MinY = ext.MinY;
539  env->MaxY = ext.MaxY;
540 }
541 
542 LWPOLY *
544  rt_envelope env
545 ) {
546  LWPOLY *npoly = NULL;
547  POINTARRAY **rings = NULL;
548  POINTARRAY *pts = NULL;
549  POINT4D p4d;
550 
551  rings = (POINTARRAY **) rtalloc(sizeof (POINTARRAY*));
552  if (!rings) {
553  rterror("rt_util_envelope_to_lwpoly: Out of memory building envelope's geometry");
554  return NULL;
555  }
556  rings[0] = ptarray_construct(0, 0, 5);
557  if (!rings[0]) {
558  rterror("rt_util_envelope_to_lwpoly: Out of memory building envelope's geometry ring");
559  return NULL;
560  }
561 
562  pts = rings[0];
563 
564  /* Upper-left corner (first and last points) */
565  p4d.x = env.MinX;
566  p4d.y = env.MaxY;
567  ptarray_set_point4d(pts, 0, &p4d);
568  ptarray_set_point4d(pts, 4, &p4d);
569 
570  /* Upper-right corner (we go clockwise) */
571  p4d.x = env.MaxX;
572  p4d.y = env.MaxY;
573  ptarray_set_point4d(pts, 1, &p4d);
574 
575  /* Lower-right corner */
576  p4d.x = env.MaxX;
577  p4d.y = env.MinY;
578  ptarray_set_point4d(pts, 2, &p4d);
579 
580  /* Lower-left corner */
581  p4d.x = env.MinX;
582  p4d.y = env.MinY;
583  ptarray_set_point4d(pts, 3, &p4d);
584 
585  npoly = lwpoly_construct(SRID_UNKNOWN, 0, 1, rings);
586  if (npoly == NULL) {
587  rterror("rt_util_envelope_to_lwpoly: Could not build envelope's geometry");
588  return NULL;
589  }
590 
591  return npoly;
592 }
593 
594 int
595 rt_util_same_geotransform_matrix(double *gt1, double *gt2) {
596  int k = 0;
597 
598  if (gt1 == NULL || gt2 == NULL)
599  return FALSE;
600 
601  for (k = 0; k < 6; k++) {
602  if (FLT_NEQ(gt1[k], gt2[k]))
603  return FALSE;
604  }
605 
606  return TRUE;
607 }
608 
609 /* coordinates in RGB and HSV are floating point values between 0 and 1 */
611 rt_util_rgb_to_hsv(double rgb[3], double hsv[3]) {
612  int i;
613 
614  double minc;
615  double maxc;
616 
617  double h = 0.;
618  double s = 0.;
619  double v = 0.;
620 
621  minc = rgb[0];
622  maxc = rgb[0];
623 
624  /* get min and max values from RGB */
625  for (i = 1; i < 3; i++) {
626  if (rgb[i] > maxc)
627  maxc = rgb[i];
628  if (rgb[i] < minc)
629  minc = rgb[i];
630  }
631  v = maxc;
632 
633  if (maxc != minc) {
634  double diff = 0.;
635  double rc = 0.;
636  double gc = 0.;
637  double bc = 0.;
638  double junk = 0.;
639 
640  diff = maxc - minc;
641  s = diff / maxc;
642  rc = (maxc - rgb[0]) / diff;
643  gc = (maxc - rgb[1]) / diff;
644  bc = (maxc - rgb[2]) / diff;
645 
646  if (DBL_EQ(rgb[0], maxc))
647  h = bc - gc;
648  else if (DBL_EQ(rgb[1], maxc))
649  h = 2.0 + rc - bc;
650  else
651  h = 4.0 + gc - rc;
652 
653  h = modf((h / 6.0), &junk);
654  }
655 
656  hsv[0] = h;
657  hsv[1] = s;
658  hsv[2] = v;
659 
660  return ES_NONE;
661 }
662 
663 /* coordinates in RGB and HSV are floating point values between 0 and 1 */
665 rt_util_hsv_to_rgb(double hsv[3], double rgb[3]) {
666  double r = 0;
667  double g = 0;
668  double b = 0;
669  double v = hsv[2];
670 
671  if (DBL_EQ(hsv[1], 0.))
672  r = g = b = v;
673  else {
674  double i;
675  double f;
676  double p;
677  double q;
678  double t;
679 
680  int a;
681 
682  i = floor(hsv[0] * 6.);
683  f = (hsv[0] * 6.0) - i;
684  p = v * (1. - hsv[1]);
685  q = v * (1. - hsv[1] * f);
686  t = v * (1. - hsv[1] * (1. - f));
687 
688  a = (int) i;
689  switch (a) {
690  case 1:
691  r = q;
692  g = v;
693  b = p;
694  break;
695  case 2:
696  r = p;
697  g = v;
698  b = t;
699  break;
700  case 3:
701  r = p;
702  g = q;
703  b = v;
704  break;
705  case 4:
706  r = t;
707  g = p;
708  b = v;
709  break;
710  case 5:
711  r = v;
712  g = p;
713  b = q;
714  break;
715  case 0:
716  case 6:
717  default:
718  r = v;
719  g = t;
720  b = p;
721  break;
722  }
723  }
724 
725  rgb[0] = r;
726  rgb[1] = g;
727  rgb[2] = b;
728 
729  return ES_NONE;
730 }
731 
732 /*- rt_context -------------------------------------------------------*/
733 
734 /*
735  * Default allocators
736  *
737  * We include some default allocators that use malloc/free/realloc
738  * along with stdout/stderr since this is the most common use case
739  *
740  */
741 void *
743 {
744  void *mem = malloc(size);
745  return mem;
746 }
747 
748 void *
749 default_rt_reallocator(void *mem, size_t size)
750 {
751  void *ret = realloc(mem, size);
752  return ret;
753 }
754 
755 void
757 {
758  free(mem);
759 }
760 
761 void
762 default_rt_error_handler(const char *fmt, va_list ap) {
763 
764  static const char *label = "ERROR: ";
765  char newfmt[1024] = {0};
766  snprintf(newfmt, 1024, "%s%s\n", label, fmt);
767  newfmt[1023] = '\0';
768 
769  vprintf(newfmt, ap);
770 
771  va_end(ap);
772 }
773 
774 void
775 default_rt_warning_handler(const char *fmt, va_list ap) {
776 
777  static const char *label = "WARNING: ";
778  char newfmt[1024] = {0};
779  snprintf(newfmt, 1024, "%s%s\n", label, fmt);
780  newfmt[1023] = '\0';
781 
782  vprintf(newfmt, ap);
783 
784  va_end(ap);
785 }
786 
787 void
788 default_rt_info_handler(const char *fmt, va_list ap) {
789 
790  static const char *label = "INFO: ";
791  char newfmt[1024] = {0};
792  snprintf(newfmt, 1024, "%s%s\n", label, fmt);
793  newfmt[1023] = '\0';
794 
795  vprintf(newfmt, ap);
796 
797  va_end(ap);
798 }
799 
803 struct rt_context_t {
810 };
811 
812 /* Static variable, to be used for all rt_core functions */
813 static struct rt_context_t ctx_t = {
815  .realloc = default_rt_reallocator,
816  .dealloc = default_rt_deallocator,
820 };
821 
822 
829 void
831 {
832  ctx_t.alloc = default_rt_allocator;
838 }
839 
840 
845 void
847  rt_deallocator deallocator, rt_message_handler error_handler,
848  rt_message_handler info_handler, rt_message_handler warning_handler) {
849 
850  ctx_t.alloc = allocator;
851  ctx_t.realloc = reallocator;
852  ctx_t.dealloc = deallocator;
853 
854  ctx_t.err = error_handler;
855  ctx_t.info = info_handler;
856  ctx_t.warn = warning_handler;
857 }
858 
859 
860 
866 void *
867 rtalloc(size_t size) {
868  void * mem = ctx_t.alloc(size);
869  RASTER_DEBUGF(5, "rtalloc called: %d@%p", size, mem);
870  return mem;
871 }
872 
873 
874 void *
875 rtrealloc(void * mem, size_t size) {
876  void * result = ctx_t.realloc(mem, size);
877  RASTER_DEBUGF(5, "rtrealloc called: %d@%p", size, result);
878  return result;
879 }
880 
881 void
882 rtdealloc(void * mem) {
883  ctx_t.dealloc(mem);
884  RASTER_DEBUG(5, "rtdealloc called");
885 }
886 
894 void
895 rterror(const char *fmt, ...) {
896  va_list ap;
897 
898  va_start(ap, fmt);
899 
900  /* Call the supplied function */
901  (*ctx_t.err)(fmt, ap);
902 
903  va_end(ap);
904 }
905 
906 void
907 rtinfo(const char *fmt, ...) {
908  va_list ap;
909 
910  va_start(ap, fmt);
911 
912  /* Call the supplied function */
913  (*ctx_t.info)(fmt, ap);
914 
915  va_end(ap);
916 }
917 
918 
919 void
920 rtwarn(const char *fmt, ...) {
921  va_list ap;
922 
923  va_start(ap, fmt);
924 
925  /* Call the supplied function */
926  (*ctx_t.warn)(fmt, ap);
927 
928  va_end(ap);
929 }
930 
931 
932 
933 int
935  double initialvalue,
936  int32_t checkvalint, uint32_t checkvaluint,
937  float checkvalfloat, double checkvaldouble,
938  rt_pixtype pixtype
939 ) {
940  int result = 0;
941 
942  switch (pixtype) {
943  case PT_1BB:
944  case PT_2BUI:
945  case PT_4BUI:
946  case PT_8BSI:
947  case PT_8BUI:
948  case PT_16BSI:
949  case PT_16BUI:
950  case PT_32BSI: {
951  if (fabs(checkvalint - initialvalue) >= 1) {
952 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
953  rtwarn("Value set for %s band got clamped from %f to %d",
954  rt_pixtype_name(pixtype),
955  initialvalue, checkvalint
956  );
957 #endif
958  result = 1;
959  }
960  else if (FLT_NEQ(checkvalint, initialvalue)) {
961 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
962  rtwarn("Value set for %s band got truncated from %f to %d",
963  rt_pixtype_name(pixtype),
964  initialvalue, checkvalint
965  );
966 #endif
967  result = 1;
968  }
969  break;
970  }
971  case PT_32BUI: {
972  if (fabs(checkvaluint - initialvalue) >= 1) {
973 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
974  rtwarn("Value set for %s band got clamped from %f to %u",
975  rt_pixtype_name(pixtype),
976  initialvalue, checkvaluint
977  );
978 #endif
979  result = 1;
980  }
981  else if (FLT_NEQ(checkvaluint, initialvalue)) {
982 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
983  rtwarn("Value set for %s band got truncated from %f to %u",
984  rt_pixtype_name(pixtype),
985  initialvalue, checkvaluint
986  );
987 #endif
988  result = 1;
989  }
990  break;
991  }
992  case PT_32BF: {
993  /*
994  For float, because the initial value is a double,
995  there is very often a difference between the desired value and the obtained one
996  */
997  if (FLT_NEQ(checkvalfloat, initialvalue)) {
998 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
999  rtwarn("Value set for %s band got converted from %f to %f",
1000  rt_pixtype_name(pixtype),
1001  initialvalue, checkvalfloat
1002  );
1003 #endif
1004  result = 1;
1005  }
1006  break;
1007  }
1008  case PT_64BF: {
1009  if (FLT_NEQ(checkvaldouble, initialvalue)) {
1010 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
1011  rtwarn("Value set for %s band got converted from %f to %f",
1012  rt_pixtype_name(pixtype),
1013  initialvalue, checkvaldouble
1014  );
1015 #endif
1016  result = 1;
1017  }
1018  break;
1019  }
1020  case PT_END:
1021  break;
1022  }
1023 
1024  return result;
1025 }
1026 
1027 /*--- Debug and Testing Utilities --------------------------------------------*/
1028 
1029 #if POSTGIS_DEBUG_LEVEL > 2
1030 
1031 static char*
1032 d_binary_to_hex(const uint8_t * const raw, uint32_t size, uint32_t *hexsize) {
1033  char* hex = NULL;
1034  uint32_t i = 0;
1035 
1036 
1037  assert(NULL != raw);
1038  assert(NULL != hexsize);
1039 
1040 
1041  *hexsize = size * 2; /* hex is 2 times bytes */
1042  hex = (char*) rtalloc((*hexsize) + 1);
1043  if (!hex) {
1044  rterror("d_binary_to_hex: Out of memory hexifying raw binary");
1045  return NULL;
1046  }
1047  hex[*hexsize] = '\0'; /* Null-terminate */
1048 
1049  for (i = 0; i < size; ++i) {
1050  deparse_hex(raw[i], &(hex[2 * i]));
1051  }
1052 
1053  assert(NULL != hex);
1054  assert(0 == strlen(hex) % 2);
1055  return hex;
1056 }
1057 
1058 static void
1059 d_print_binary_hex(const char* msg, const uint8_t * const raw, uint32_t size) {
1060  char* hex = NULL;
1061  uint32_t hexsize = 0;
1062 
1063 
1064  assert(NULL != msg);
1065  assert(NULL != raw);
1066 
1067 
1068  hex = d_binary_to_hex(raw, size, &hexsize);
1069  if (NULL != hex) {
1070  rtinfo("%s\t%s", msg, hex);
1071  rtdealloc(hex);
1072  }
1073 }
1074 
1075 static size_t
1076 d_binptr_to_pos(const uint8_t * const ptr, const uint8_t * const end, size_t size) {
1077  assert(NULL != ptr && NULL != end);
1078 
1079  return (size - (end - ptr));
1080 }
1081 
1082 #define CHECK_BINPTR_POSITION(ptr, end, size, pos) \
1083  { if (pos != d_binptr_to_pos(ptr, end, size)) { \
1084  fprintf(stderr, "Check of binary stream pointer position failed on line %d\n", __LINE__); \
1085  fprintf(stderr, "\tactual = %lu, expected = %lu\n", (long unsigned)d_binptr_to_pos(ptr, end, size), (long unsigned)pos); \
1086  } }
1087 
1088 #else
1089 
1090 #define CHECK_BINPTR_POSITION(ptr, end, size, pos) ((void)0);
1091 
1092 #endif /* ifndef POSTGIS_DEBUG_LEVEL > 3 */
1093 
1094 /*- rt_pixeltype -----------------------------------------------------*/
1095 
1096 int
1098  int pixbytes = -1;
1099 
1100  switch (pixtype) {
1101  case PT_1BB:
1102  case PT_2BUI:
1103  case PT_4BUI:
1104  case PT_8BSI:
1105  case PT_8BUI:
1106  pixbytes = 1;
1107  break;
1108  case PT_16BSI:
1109  case PT_16BUI:
1110  pixbytes = 2;
1111  break;
1112  case PT_32BSI:
1113  case PT_32BUI:
1114  case PT_32BF:
1115  pixbytes = 4;
1116  break;
1117  case PT_64BF:
1118  pixbytes = 8;
1119  break;
1120  default:
1121  rterror("rt_pixtype_size: Unknown pixeltype %d", pixtype);
1122  pixbytes = -1;
1123  break;
1124  }
1125 
1126  RASTER_DEBUGF(3, "Pixel type = %s and size = %d bytes",
1127  rt_pixtype_name(pixtype), pixbytes);
1128 
1129  return pixbytes;
1130 }
1131 
1132 int
1134  return rt_pixtype_size(pixtype);
1135 }
1136 
1137 rt_pixtype
1138 rt_pixtype_index_from_name(const char* pixname) {
1139  assert(pixname && strlen(pixname) > 0);
1140 
1141  if (strcmp(pixname, "1BB") == 0)
1142  return PT_1BB;
1143  else if (strcmp(pixname, "2BUI") == 0)
1144  return PT_2BUI;
1145  else if (strcmp(pixname, "4BUI") == 0)
1146  return PT_4BUI;
1147  else if (strcmp(pixname, "8BSI") == 0)
1148  return PT_8BSI;
1149  else if (strcmp(pixname, "8BUI") == 0)
1150  return PT_8BUI;
1151  else if (strcmp(pixname, "16BSI") == 0)
1152  return PT_16BSI;
1153  else if (strcmp(pixname, "16BUI") == 0)
1154  return PT_16BUI;
1155  else if (strcmp(pixname, "32BSI") == 0)
1156  return PT_32BSI;
1157  else if (strcmp(pixname, "32BUI") == 0)
1158  return PT_32BUI;
1159  else if (strcmp(pixname, "32BF") == 0)
1160  return PT_32BF;
1161  else if (strcmp(pixname, "64BF") == 0)
1162  return PT_64BF;
1163 
1164  return PT_END;
1165 }
1166 
1167 const char*
1169  switch (pixtype) {
1170  case PT_1BB:
1171  return "1BB";
1172  case PT_2BUI:
1173  return "2BUI";
1174  case PT_4BUI:
1175  return "4BUI";
1176  case PT_8BSI:
1177  return "8BSI";
1178  case PT_8BUI:
1179  return "8BUI";
1180  case PT_16BSI:
1181  return "16BSI";
1182  case PT_16BUI:
1183  return "16BUI";
1184  case PT_32BSI:
1185  return "32BSI";
1186  case PT_32BUI:
1187  return "32BUI";
1188  case PT_32BF:
1189  return "32BF";
1190  case PT_64BF:
1191  return "64BF";
1192  default:
1193  rterror("rt_pixtype_name: Unknown pixeltype %d", pixtype);
1194  return "Unknown";
1195  }
1196 }
1197 
1205 double
1207  switch (pixtype) {
1208  case PT_1BB: {
1209  return (double) rt_util_clamp_to_1BB((double) CHAR_MIN);
1210  }
1211  case PT_2BUI: {
1212  return (double) rt_util_clamp_to_2BUI((double) CHAR_MIN);
1213  }
1214  case PT_4BUI: {
1215  return (double) rt_util_clamp_to_4BUI((double) CHAR_MIN);
1216  }
1217  case PT_8BUI: {
1218  return (double) rt_util_clamp_to_8BUI((double) CHAR_MIN);
1219  }
1220  case PT_8BSI: {
1221  return (double) rt_util_clamp_to_8BSI((double) SCHAR_MIN);
1222  }
1223  case PT_16BSI: {
1224  return (double) rt_util_clamp_to_16BSI((double) SHRT_MIN);
1225  }
1226  case PT_16BUI: {
1227  return (double) rt_util_clamp_to_16BUI((double) SHRT_MIN);
1228  }
1229  case PT_32BSI: {
1230  return (double) rt_util_clamp_to_32BSI((double) INT_MIN);
1231  }
1232  case PT_32BUI: {
1233  return (double) rt_util_clamp_to_32BUI((double) INT_MIN);
1234  }
1235  case PT_32BF: {
1236  return (double) -FLT_MAX;
1237  }
1238  case PT_64BF: {
1239  return (double) -DBL_MAX;
1240  }
1241  default: {
1242  rterror("rt_pixtype_get_min_value: Unknown pixeltype %d", pixtype);
1243  return (double) rt_util_clamp_to_8BUI((double) CHAR_MIN);
1244  }
1245  }
1246 }
1247 
1259  rt_pixtype pixtype,
1260  double val, double refval,
1261  int *isequal
1262 ) {
1263  assert(isequal != NULL);
1264  *isequal = 0;
1265 
1266  switch (pixtype) {
1267  case PT_1BB:
1268  if (rt_util_clamp_to_1BB(val) == rt_util_clamp_to_1BB(refval))
1269  *isequal = 1;
1270  break;
1271  case PT_2BUI:
1272  if (rt_util_clamp_to_2BUI(val) == rt_util_clamp_to_2BUI(refval))
1273  *isequal = 1;
1274  break;
1275  case PT_4BUI:
1276  if (rt_util_clamp_to_4BUI(val) == rt_util_clamp_to_4BUI(refval))
1277  *isequal = 1;
1278  break;
1279  case PT_8BSI:
1280  if (rt_util_clamp_to_8BSI(val) == rt_util_clamp_to_8BSI(refval))
1281  *isequal = 1;
1282  break;
1283  case PT_8BUI:
1284  if (rt_util_clamp_to_8BUI(val) == rt_util_clamp_to_8BUI(refval))
1285  *isequal = 1;
1286  break;
1287  case PT_16BSI:
1288  if (rt_util_clamp_to_16BSI(val) == rt_util_clamp_to_16BSI(refval))
1289  *isequal = 1;
1290  break;
1291  case PT_16BUI:
1292  if (rt_util_clamp_to_16BUI(val) == rt_util_clamp_to_16BUI(refval))
1293  *isequal = 1;
1294  break;
1295  case PT_32BSI:
1296  if (rt_util_clamp_to_32BSI(val) == rt_util_clamp_to_32BSI(refval))
1297  *isequal = 1;
1298  break;
1299  case PT_32BUI:
1300  if (rt_util_clamp_to_32BUI(val) == rt_util_clamp_to_32BUI(refval))
1301  *isequal = 1;
1302  break;
1303  case PT_32BF:
1305  *isequal = 1;
1306  break;
1307  case PT_64BF:
1308  if (FLT_EQ(val, refval))
1309  *isequal = 1;
1310  break;
1311  default:
1312  rterror("rt_pixtype_compare_clamped_values: Unknown pixeltype %d", pixtype);
1313  return ES_ERROR;
1314  }
1315 
1316  return ES_NONE;
1317 }
1318 
1319 /*- rt_pixel ----------------------------------------------------------*/
1320 
1321 /*
1322  * Convert an array of rt_pixel objects to two 2D arrays of value and NODATA.
1323  * The dimensions of the returned 2D array are [Y][X], going by row Y and
1324  * then column X.
1325  *
1326  * @param npixel : array of rt_pixel objects
1327  * @param count : number of elements in npixel
1328  * @param x : the column of the center pixel (0-based)
1329  * @param y : the line of the center pixel (0-based)
1330  * @param distancex : the number of pixels around the specified pixel
1331  * along the X axis
1332  * @param distancey : the number of pixels around the specified pixel
1333  * along the Y axis
1334  * @param value : pointer to pointer for 2D value array
1335  * @param nodata : pointer to pointer for 2D NODATA array
1336  * @param dimx : size of value and nodata along the X axis
1337  * @param dimy : size of value and nodata along the Y axis
1338  *
1339  * @return ES_NONE on success, ES_ERROR on error
1340  */
1342  rt_pixel npixel, int count,
1343  int x, int y,
1344  uint16_t distancex, uint16_t distancey,
1345  double ***value,
1346  int ***nodata,
1347  int *dimx, int *dimy
1348 ) {
1349  uint32_t i;
1350  uint32_t j;
1351  uint32_t dim[2] = {0};
1352  double **values = NULL;
1353  int **nodatas = NULL;
1354  int zero[2] = {0};
1355  int _x;
1356  int _y;
1357 
1358  assert(npixel != NULL && count > 0);
1359  assert(value != NULL);
1360  assert(nodata != NULL);
1361 
1362  /* dimensions */
1363  dim[0] = distancex * 2 + 1;
1364  dim[1] = distancey * 2 + 1;
1365  RASTER_DEBUGF(4, "dimensions = %d x %d", dim[0], dim[1]);
1366 
1367  /* establish 2D arrays (Y axis) */
1368  values = rtalloc(sizeof(double *) * dim[1]);
1369  nodatas = rtalloc(sizeof(int *) * dim[1]);
1370 
1371  if (values == NULL || nodatas == NULL) {
1372  rterror("rt_pixel_set_to_array: Could not allocate memory for 2D array");
1373  return ES_ERROR;
1374  }
1375 
1376  /* initialize X axis */
1377  for (i = 0; i < dim[1]; i++) {
1378  values[i] = rtalloc(sizeof(double) * dim[0]);
1379  nodatas[i] = rtalloc(sizeof(int) * dim[0]);
1380 
1381  if (values[i] == NULL || nodatas[i] == NULL) {
1382  rterror("rt_pixel_set_to_array: Could not allocate memory for dimension of 2D array");
1383 
1384  if (values[i] == NULL) {
1385  for (j = 0; j < i; j++) {
1386  rtdealloc(values[j]);
1387  rtdealloc(nodatas[j]);
1388  }
1389  }
1390  else {
1391  for (j = 0; j <= i; j++) {
1392  rtdealloc(values[j]);
1393  if (j < i)
1394  rtdealloc(nodatas[j]);
1395  }
1396  }
1397 
1398  rtdealloc(values);
1399  rtdealloc(nodatas);
1400 
1401  return ES_ERROR;
1402  }
1403 
1404  /* set values to 0 */
1405  memset(values[i], 0, sizeof(double) * dim[0]);
1406 
1407  /* set nodatas to 1 */
1408  for (j = 0; j < dim[0]; j++)
1409  nodatas[i][j] = 1;
1410  }
1411 
1412  /* get 0,0 of grid */
1413  zero[0] = x - distancex;
1414  zero[1] = y - distancey;
1415 
1416  /* populate 2D arrays */
1417  for (i = 0; i < count; i++) {
1418  if (npixel[i].nodata)
1419  continue;
1420 
1421  _x = npixel[i].x - zero[0];
1422  _y = npixel[i].y - zero[1];
1423 
1424  RASTER_DEBUGF(4, "absolute x,y: %d x %d", npixel[i].x, npixel[i].y);
1425  RASTER_DEBUGF(4, "relative x,y: %d x %d", _x, _y);
1426 
1427  values[_y][_x] = npixel[i].value;
1428  nodatas[_y][_x] = 0;
1429 
1430  RASTER_DEBUGF(4, "(x, y, nodata, value) = (%d, %d, %d, %f)", _x, _y, nodatas[_y][_x], values[_y][_x]);
1431  }
1432 
1433  *value = &(*values);
1434  *nodata = &(*nodatas);
1435  if (dimx != NULL)
1436  *dimx = dim[0];
1437  if (dimy != NULL)
1438  *dimy = dim[1];
1439 
1440  return ES_NONE;
1441 }
1442 
1443 /*- rt_band ----------------------------------------------------------*/
1444 
1465 rt_band
1467  uint16_t width, uint16_t height,
1468  rt_pixtype pixtype,
1469  uint32_t hasnodata, double nodataval,
1470  uint8_t* data
1471 ) {
1472  rt_band band = NULL;
1473 
1474  assert(NULL != data);
1475 
1476  band = rtalloc(sizeof(struct rt_band_t));
1477  if (band == NULL) {
1478  rterror("rt_band_new_inline: Out of memory allocating rt_band");
1479  return NULL;
1480  }
1481 
1482  RASTER_DEBUGF(3, "Created rt_band @ %p with pixtype %s", band, rt_pixtype_name(pixtype));
1483 
1484  band->pixtype = pixtype;
1485  band->offline = 0;
1486  band->width = width;
1487  band->height = height;
1488  band->hasnodata = hasnodata ? 1 : 0;
1489  band->isnodata = FALSE; /* we don't know what is in data, so must be FALSE */
1490  band->nodataval = 0;
1491  band->data.mem = data;
1492  band->ownsdata = 0; /* we do NOT own this data!!! */
1493  band->raster = NULL;
1494 
1495  RASTER_DEBUGF(3, "Created rt_band with dimensions %d x %d", band->width, band->height);
1496 
1497  /* properly set nodataval as it may need to be constrained to the data type */
1498  if (hasnodata && rt_band_set_nodata(band, nodataval, NULL) != ES_NONE) {
1499  rterror("rt_band_new_inline: Could not set NODATA value");
1500  rt_band_destroy(band);
1501  return NULL;
1502  }
1503 
1504  return band;
1505 }
1506 
1526 rt_band
1528  uint16_t width, uint16_t height,
1529  rt_pixtype pixtype,
1530  uint32_t hasnodata, double nodataval,
1531  uint8_t bandNum, const char* path
1532 ) {
1533  rt_band band = NULL;
1534  int pathlen = 0;
1535 
1536  assert(NULL != path);
1537 
1538  band = rtalloc(sizeof(struct rt_band_t));
1539  if (band == NULL) {
1540  rterror("rt_band_new_offline: Out of memory allocating rt_band");
1541  return NULL;
1542  }
1543 
1544  RASTER_DEBUGF(3, "Created rt_band @ %p with pixtype %s",
1545  band, rt_pixtype_name(pixtype)
1546  );
1547 
1548  band->pixtype = pixtype;
1549  band->offline = 1;
1550  band->width = width;
1551  band->height = height;
1552  band->hasnodata = hasnodata ? 1 : 0;
1553  band->nodataval = 0;
1554  band->isnodata = FALSE; /* we don't know if the offline band is NODATA */
1555  band->ownsdata = 0; /* offline, flag is useless as all offline data cache is owned internally */
1556  band->raster = NULL;
1557 
1558  /* properly set nodataval as it may need to be constrained to the data type */
1559  if (hasnodata && rt_band_set_nodata(band, nodataval, NULL) != ES_NONE) {
1560  rterror("rt_band_new_offline: Could not set NODATA value");
1561  rt_band_destroy(band);
1562  return NULL;
1563  }
1564 
1565  band->data.offline.bandNum = bandNum;
1566 
1567  /* memory for data.offline.path is managed internally */
1568  pathlen = strlen(path);
1569  band->data.offline.path = rtalloc(sizeof(char) * (pathlen + 1));
1570  if (band->data.offline.path == NULL) {
1571  rterror("rt_band_new_offline: Out of memory allocating offline path");
1572  rt_band_destroy(band);
1573  return NULL;
1574  }
1575  memcpy(band->data.offline.path, path, pathlen);
1576  band->data.offline.path[pathlen] = '\0';
1577 
1578  band->data.offline.mem = NULL;
1579 
1580  return band;
1581 }
1582 
1593 rt_band
1595  rt_band rtn = NULL;
1596 
1597  assert(band != NULL);
1598 
1599  /* offline */
1600  if (band->offline) {
1601  rtn = rt_band_new_offline(
1602  band->width, band->height,
1603  band->pixtype,
1604  band->hasnodata, band->nodataval,
1605  band->data.offline.bandNum, (const char *) band->data.offline.path
1606  );
1607  }
1608  /* online */
1609  else {
1610  uint8_t *data = NULL;
1611  data = rtalloc(rt_pixtype_size(band->pixtype) * band->width * band->height);
1612  if (data == NULL) {
1613  rterror("rt_band_duplicate: Out of memory allocating online band data");
1614  return NULL;
1615  }
1616  memcpy(data, band->data.mem, rt_pixtype_size(band->pixtype) * band->width * band->height);
1617 
1618  rtn = rt_band_new_inline(
1619  band->width, band->height,
1620  band->pixtype,
1621  band->hasnodata, band->nodataval,
1622  data
1623  );
1624  rt_band_set_ownsdata_flag(rtn, 1); /* we DO own this data!!! */
1625  }
1626 
1627  if (rtn == NULL) {
1628  rterror("rt_band_duplicate: Could not copy band");
1629  return NULL;
1630  }
1631 
1632  return rtn;
1633 }
1634 
1635 int
1637 
1638  assert(NULL != band);
1639 
1640 
1641  return band->offline ? 1 : 0;
1642 }
1643 
1649 void
1651  if (band == NULL)
1652  return;
1653 
1654  RASTER_DEBUGF(3, "Destroying rt_band @ %p", band);
1655 
1656  /* offline band */
1657  if (band->offline) {
1658  /* memory cache */
1659  if (band->data.offline.mem != NULL)
1660  rtdealloc(band->data.offline.mem);
1661  /* offline file path */
1662  if (band->data.offline.path != NULL)
1663  rtdealloc(band->data.offline.path);
1664  }
1665  /* inline band and band owns the data */
1666  else if (band->data.mem != NULL && band->ownsdata)
1667  rtdealloc(band->data.mem);
1668 
1669  rtdealloc(band);
1670 }
1671 
1672 const char*
1674 
1675  assert(NULL != band);
1676 
1677 
1678  if (!band->offline) {
1679  RASTER_DEBUG(3, "rt_band_get_ext_path: Band is not offline");
1680  return NULL;
1681  }
1682  return band->data.offline.path;
1683 }
1684 
1687  assert(NULL != band);
1688  assert(NULL != bandnum);
1689 
1690  *bandnum = 0;
1691 
1692  if (!band->offline) {
1693  RASTER_DEBUG(3, "rt_band_get_ext_band_num: Band is not offline");
1694  return ES_ERROR;
1695  }
1696 
1697  *bandnum = band->data.offline.bandNum;
1698 
1699  return ES_NONE;
1700 }
1701 
1709 void *
1711  assert(NULL != band);
1712 
1713  if (band->offline) {
1714  if (band->data.offline.mem != NULL)
1715  return band->data.offline.mem;
1716 
1717  if (rt_band_load_offline_data(band) != ES_NONE)
1718  return NULL;
1719  else
1720  return band->data.offline.mem;
1721  }
1722  else
1723  return band->data.mem;
1724 }
1725 
1726 /* variable for PostgreSQL env variable: POSTGIS_ENABLE_OUTDB_RASTERS */
1728 
1740  GDALDatasetH hdsSrc = NULL;
1741  int nband = 0;
1742  VRTDatasetH hdsDst = NULL;
1743  VRTSourcedRasterBandH hbandDst = NULL;
1744  double gt[6] = {0.};
1745  double ogt[6] = {0};
1746  double offset[2] = {0};
1747 
1748  rt_raster _rast = NULL;
1749  rt_band _band = NULL;
1750  int aligned = 0;
1751  int err = ES_NONE;
1752 
1753  assert(band != NULL);
1754  assert(band->raster != NULL);
1755 
1756  if (!band->offline) {
1757  rterror("rt_band_load_offline_data: Band is not offline");
1758  return ES_ERROR;
1759  }
1760  else if (!strlen(band->data.offline.path)) {
1761  rterror("rt_band_load_offline_data: Offline band does not a have a specified file");
1762  return ES_ERROR;
1763  }
1764 
1765  /* out-db is disabled */
1766  if (!enable_outdb_rasters) {
1767  rterror("rt_raster_load_offline_data: Access to offline bands disabled");
1768  return ES_ERROR;
1769  }
1770 
1772  /*
1773  hdsSrc = rt_util_gdal_open(band->data.offline.path, GA_ReadOnly, 1);
1774  */
1775  hdsSrc = rt_util_gdal_open(band->data.offline.path, GA_ReadOnly, 0);
1776  if (hdsSrc == NULL) {
1777  rterror("rt_band_load_offline_data: Cannot open offline raster: %s", band->data.offline.path);
1778  return ES_ERROR;
1779  }
1780 
1781  /* # of bands */
1782  nband = GDALGetRasterCount(hdsSrc);
1783  if (!nband) {
1784  rterror("rt_band_load_offline_data: No bands found in offline raster: %s", band->data.offline.path);
1785  GDALClose(hdsSrc);
1786  return ES_ERROR;
1787  }
1788  /* bandNum is 0-based */
1789  else if (band->data.offline.bandNum + 1 > nband) {
1790  rterror("rt_band_load_offline_data: Specified band %d not found in offline raster: %s", band->data.offline.bandNum, band->data.offline.path);
1791  GDALClose(hdsSrc);
1792  return ES_ERROR;
1793  }
1794 
1795  /* get raster's geotransform */
1797  RASTER_DEBUGF(3, "Raster geotransform (%f, %f, %f, %f, %f, %f)",
1798  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
1799 
1800  /* get offline raster's geotransform */
1801  if (GDALGetGeoTransform(hdsSrc, ogt) != CE_None) {
1802  RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)");
1803  ogt[0] = 0;
1804  ogt[1] = 1;
1805  ogt[2] = 0;
1806  ogt[3] = 0;
1807  ogt[4] = 0;
1808  ogt[5] = -1;
1809  }
1810  RASTER_DEBUGF(3, "Offline geotransform (%f, %f, %f, %f, %f, %f)",
1811  ogt[0], ogt[1], ogt[2], ogt[3], ogt[4], ogt[5]);
1812 
1813  /* are rasters aligned? */
1814  _rast = rt_raster_new(1, 1);
1816  rt_raster_set_srid(_rast, band->raster->srid);
1817  err = rt_raster_same_alignment(band->raster, _rast, &aligned, NULL);
1818  rt_raster_destroy(_rast);
1819 
1820  if (err != ES_NONE) {
1821  rterror("rt_band_load_offline_data: Could not test alignment of in-db representation of out-db raster");
1822  GDALClose(hdsSrc);
1823  return ES_ERROR;
1824  }
1825  else if (!aligned) {
1826  rtwarn("The in-db representation of the out-db raster is not aligned. Band data may be incorrect");
1827  }
1828 
1829  /* get offsets */
1831  band->raster,
1832  ogt[0], ogt[3],
1833  &(offset[0]), &(offset[1]),
1834  NULL
1835  );
1836 
1837  RASTER_DEBUGF(4, "offsets: (%f, %f)", offset[0], offset[1]);
1838 
1839  /* create VRT dataset */
1840  hdsDst = VRTCreate(band->width, band->height);
1841  GDALSetGeoTransform(hdsDst, gt);
1842  /*
1843  GDALSetDescription(hdsDst, "/tmp/offline.vrt");
1844  */
1845 
1846  /* add band as simple sources */
1847  GDALAddBand(hdsDst, rt_util_pixtype_to_gdal_datatype(band->pixtype), NULL);
1848  hbandDst = (VRTSourcedRasterBandH) GDALGetRasterBand(hdsDst, 1);
1849 
1850  if (band->hasnodata)
1851  GDALSetRasterNoDataValue(hbandDst, band->nodataval);
1852 
1853  VRTAddSimpleSource(
1854  hbandDst, GDALGetRasterBand(hdsSrc, band->data.offline.bandNum + 1),
1855  fabs(offset[0]), fabs(offset[1]),
1856  band->width, band->height,
1857  0, 0,
1858  band->width, band->height,
1859  "near", VRT_NODATA_UNSET
1860  );
1861 
1862  /* make sure VRT reflects all changes */
1863  VRTFlushCache(hdsDst);
1864 
1865  /* convert VRT dataset to rt_raster */
1866  _rast = rt_raster_from_gdal_dataset(hdsDst);
1867 
1868  GDALClose(hdsDst);
1869  /* XXX: need to find a way to clean up the GDALOpenShared datasets at end of transaction */
1870  /* GDALClose(hdsSrc); */
1871  GDALClose(hdsSrc);
1872 
1873  if (_rast == NULL) {
1874  rterror("rt_band_load_offline_data: Cannot load data from offline raster: %s", band->data.offline.path);
1875  return ES_ERROR;
1876  }
1877 
1878  _band = rt_raster_get_band(_rast, 0);
1879  if (_band == NULL) {
1880  rterror("rt_band_load_offline_data: Cannot load data from offline raster: %s", band->data.offline.path);
1881  rt_raster_destroy(_rast);
1882  return ES_ERROR;
1883  }
1884 
1885  /* band->data.offline.mem not NULL, free first */
1886  if (band->data.offline.mem != NULL) {
1887  rtdealloc(band->data.offline.mem);
1888  band->data.offline.mem = NULL;
1889  }
1890 
1891  band->data.offline.mem = _band->data.mem;
1892 
1893  rtdealloc(_band); /* cannot use rt_band_destroy */
1894  rt_raster_destroy(_rast);
1895 
1896  return ES_NONE;
1897 }
1898 
1899 rt_pixtype
1901 
1902  assert(NULL != band);
1903 
1904 
1905  return band->pixtype;
1906 }
1907 
1908 uint16_t
1910 
1911  assert(NULL != band);
1912 
1913 
1914  return band->width;
1915 }
1916 
1917 uint16_t
1919 
1920  assert(NULL != band);
1921 
1922 
1923  return band->height;
1924 }
1925 
1926 /* Get ownsdata flag */
1927 int
1929  assert(NULL != band);
1930 
1931  return band->ownsdata ? 1 : 0;
1932 }
1933 
1934 /* set ownsdata flag */
1935 void
1937  assert(NULL != band);
1938 
1939  band->ownsdata = flag ? 1 : 0;
1940 }
1941 
1942 #ifdef OPTIMIZE_SPACE
1943 
1944 /*
1945  * Set given number of bits of the given byte,
1946  * starting from given bitOffset (from the first)
1947  * to the given value.
1948  *
1949  * Examples:
1950  * char ch;
1951  * ch=0; setBits(&ch, 1, 1, 0) -> ch==8
1952  * ch=0; setBits(&ch, 3, 2, 1) -> ch==96 (0x60)
1953  *
1954  * Note that number of bits set must be <= 8-bitOffset
1955  *
1956  */
1957 static void
1958 setBits(char* ch, double val, int bits, int bitOffset) {
1959  char mask = 0xFF >> (8 - bits);
1960  char ival = val;
1961 
1962 
1963  assert(ch != NULL);
1964  assert(8 - bitOffset >= bits);
1965 
1966  RASTER_DEBUGF(4, "ival:%d bits:%d mask:%hhx bitoffset:%d\n",
1967  ival, bits, mask, bitOffset);
1968 
1969  /* clear all but significant bits from ival */
1970  ival &= mask;
1971 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
1972  if (ival != val) {
1973  rtwarn("Pixel value for %d-bits band got truncated"
1974  " from %g to %hhu", bits, val, ival);
1975  }
1976 #endif /* POSTGIS_RASTER_WARN_ON_TRUNCATION */
1977 
1978  RASTER_DEBUGF(4, " cleared ival:%hhx\n", ival);
1979 
1980 
1981  /* Shift ival so the significant bits start at
1982  * the first bit */
1983  ival <<= (8 - bitOffset - bits);
1984 
1985  RASTER_DEBUGF(4, " ival shifted:%hhx\n", ival);
1986  RASTER_DEBUGF(4, " ch:%hhx\n", *ch);
1987 
1988  /* clear first bits of target */
1989  *ch &= ~(mask << (8 - bits - bitOffset));
1990 
1991  RASTER_DEBUGF(4, " ch cleared:%hhx\n", *ch);
1992 
1993  /* Set the first bit of target */
1994  *ch |= ival;
1995 
1996  RASTER_DEBUGF(4, " ch ored:%hhx\n", *ch);
1997 
1998 }
1999 #endif /* OPTIMIZE_SPACE */
2000 
2001 int
2003  assert(NULL != band);
2004 
2005  return band->hasnodata ? 1 : 0;
2006 }
2007 
2008 void
2010 
2011  assert(NULL != band);
2012 
2013  band->hasnodata = (flag) ? 1 : 0;
2014 
2015  /* isnodata depends on hasnodata */
2016  if (!band->hasnodata && band->isnodata) {
2017  RASTER_DEBUG(3, "Setting isnodata to FALSE as band no longer has NODATA");
2018  band->isnodata = 0;
2019  }
2020 }
2021 
2024  assert(NULL != band);
2025 
2026  if (!band->hasnodata) {
2027  /* silently permit setting isnodata flag to FALSE */
2028  if (!flag)
2029  band->isnodata = 0;
2030  else {
2031  rterror("rt_band_set_isnodata_flag: Cannot set isnodata flag as band has no NODATA");
2032  return ES_ERROR;
2033  }
2034  }
2035  else
2036  band->isnodata = (flag) ? 1 : 0;
2037 
2038  return ES_NONE;
2039 }
2040 
2041 int
2043  assert(NULL != band);
2044 
2045  if (band->hasnodata)
2046  return band->isnodata ? 1 : 0;
2047  else
2048  return 0;
2049 }
2050 
2061 rt_band_set_nodata(rt_band band, double val, int *converted) {
2062  rt_pixtype pixtype = PT_END;
2063  int32_t checkvalint = 0;
2064  uint32_t checkvaluint = 0;
2065  float checkvalfloat = 0;
2066  double checkvaldouble = 0;
2067 
2068  assert(NULL != band);
2069 
2070  if (converted != NULL)
2071  *converted = 0;
2072 
2073  pixtype = band->pixtype;
2074 
2075  RASTER_DEBUGF(3, "rt_band_set_nodata: setting nodata value %g with band type %s", val, rt_pixtype_name(pixtype));
2076 
2077  /* return -1 on out of range */
2078  switch (pixtype) {
2079  case PT_1BB: {
2080  band->nodataval = rt_util_clamp_to_1BB(val);
2081  checkvalint = band->nodataval;
2082  break;
2083  }
2084  case PT_2BUI: {
2085  band->nodataval = rt_util_clamp_to_2BUI(val);
2086  checkvalint = band->nodataval;
2087  break;
2088  }
2089  case PT_4BUI: {
2090  band->nodataval = rt_util_clamp_to_4BUI(val);
2091  checkvalint = band->nodataval;
2092  break;
2093  }
2094  case PT_8BSI: {
2095  band->nodataval = rt_util_clamp_to_8BSI(val);
2096  checkvalint = band->nodataval;
2097  break;
2098  }
2099  case PT_8BUI: {
2100  band->nodataval = rt_util_clamp_to_8BUI(val);
2101  checkvalint = band->nodataval;
2102  break;
2103  }
2104  case PT_16BSI: {
2105  band->nodataval = rt_util_clamp_to_16BSI(val);
2106  checkvalint = band->nodataval;
2107  break;
2108  }
2109  case PT_16BUI: {
2110  band->nodataval = rt_util_clamp_to_16BUI(val);
2111  checkvalint = band->nodataval;
2112  break;
2113  }
2114  case PT_32BSI: {
2115  band->nodataval = rt_util_clamp_to_32BSI(val);
2116  checkvalint = band->nodataval;
2117  break;
2118  }
2119  case PT_32BUI: {
2120  band->nodataval = rt_util_clamp_to_32BUI(val);
2121  checkvaluint = band->nodataval;
2122  break;
2123  }
2124  case PT_32BF: {
2125  band->nodataval = rt_util_clamp_to_32F(val);
2126  checkvalfloat = band->nodataval;
2127  break;
2128  }
2129  case PT_64BF: {
2130  band->nodataval = val;
2131  checkvaldouble = band->nodataval;
2132  break;
2133  }
2134  default: {
2135  rterror("rt_band_set_nodata: Unknown pixeltype %d", pixtype);
2136  band->hasnodata = 0;
2137  return ES_ERROR;
2138  }
2139  }
2140 
2141  RASTER_DEBUGF(3, "rt_band_set_nodata: band->hasnodata = %d", band->hasnodata);
2142  RASTER_DEBUGF(3, "rt_band_set_nodata: band->nodataval = %f", band->nodataval);
2143  /* the nodata value was just set, so this band has NODATA */
2144  band->hasnodata = 1;
2145 
2146  /* also set isnodata flag to false */
2147  band->isnodata = 0;
2148 
2150  val,
2151  checkvalint, checkvaluint,
2152  checkvalfloat, checkvaldouble,
2153  pixtype
2154  ) && converted != NULL) {
2155  *converted = 1;
2156  }
2157 
2158  return ES_NONE;
2159 }
2160 
2182  rt_band band,
2183  int x, int y,
2184  void *vals, uint32_t len
2185 ) {
2186  rt_pixtype pixtype = PT_END;
2187  int size = 0;
2188  uint8_t *data = NULL;
2189  uint32_t offset = 0;
2190 
2191  assert(NULL != band);
2192  assert(vals != NULL && len > 0);
2193 
2194  RASTER_DEBUGF(3, "length of values = %d", len);
2195 
2196  if (band->offline) {
2197  rterror("rt_band_set_pixel_line not implemented yet for OFFDB bands");
2198  return ES_ERROR;
2199  }
2200 
2201  pixtype = band->pixtype;
2202  size = rt_pixtype_size(pixtype);
2203 
2204  if (
2205  x < 0 || x >= band->width ||
2206  y < 0 || y >= band->height
2207  ) {
2208  rterror("rt_band_set_pixel_line: Coordinates out of range (%d, %d) vs (%d, %d)", x, y, band->width, band->height);
2209  return ES_ERROR;
2210  }
2211 
2212  data = rt_band_get_data(band);
2213  offset = x + (y * band->width);
2214  RASTER_DEBUGF(4, "offset = %d", offset);
2215 
2216  /* make sure len of values to copy don't exceed end of data */
2217  if (len > (band->width * band->height) - offset) {
2218  rterror("rt_band_set_pixel_line: Could not apply pixels as values length exceeds end of data");
2219  return ES_ERROR;
2220  }
2221 
2222  switch (pixtype) {
2223  case PT_1BB:
2224  case PT_2BUI:
2225  case PT_4BUI:
2226  case PT_8BUI:
2227  case PT_8BSI: {
2228  uint8_t *ptr = data;
2229  ptr += offset;
2230  memcpy(ptr, vals, size * len);
2231  break;
2232  }
2233  case PT_16BUI: {
2234  uint16_t *ptr = (uint16_t *) data;
2235  ptr += offset;
2236  memcpy(ptr, vals, size * len);
2237  break;
2238  }
2239  case PT_16BSI: {
2240  int16_t *ptr = (int16_t *) data;
2241  ptr += offset;
2242  memcpy(ptr, vals, size * len);
2243  break;
2244  }
2245  case PT_32BUI: {
2246  uint32_t *ptr = (uint32_t *) data;
2247  ptr += offset;
2248  memcpy(ptr, vals, size * len);
2249  break;
2250  }
2251  case PT_32BSI: {
2252  int32_t *ptr = (int32_t *) data;
2253  ptr += offset;
2254  memcpy(ptr, vals, size * len);
2255  break;
2256  }
2257  case PT_32BF: {
2258  float *ptr = (float *) data;
2259  ptr += offset;
2260  memcpy(ptr, vals, size * len);
2261  break;
2262  }
2263  case PT_64BF: {
2264  double *ptr = (double *) data;
2265  ptr += offset;
2266  memcpy(ptr, vals, size * len);
2267  break;
2268  }
2269  default: {
2270  rterror("rt_band_set_pixel_line: Unknown pixeltype %d", pixtype);
2271  return ES_ERROR;
2272  }
2273  }
2274 
2275 #if POSTGIS_DEBUG_LEVEL > 0
2276  {
2277  double value;
2278  rt_band_get_pixel(band, x, y, &value, NULL);
2279  RASTER_DEBUGF(4, "pixel at (%d, %d) = %f", x, y, value);
2280  }
2281 #endif
2282 
2283  /* set band's isnodata flag to FALSE */
2284  if (rt_band_get_hasnodata_flag(band))
2285  rt_band_set_isnodata_flag(band, 0);
2286 
2287  return ES_NONE;
2288 }
2289 
2303  rt_band band,
2304  int x, int y,
2305  double val,
2306  int *converted
2307 ) {
2308  rt_pixtype pixtype = PT_END;
2309  unsigned char* data = NULL;
2310  uint32_t offset = 0;
2311 
2312  int32_t checkvalint = 0;
2313  uint32_t checkvaluint = 0;
2314  float checkvalfloat = 0;
2315  double checkvaldouble = 0;
2316 
2317  assert(NULL != band);
2318 
2319  if (converted != NULL)
2320  *converted = 0;
2321 
2322  if (band->offline) {
2323  rterror("rt_band_set_pixel not implemented yet for OFFDB bands");
2324  return ES_ERROR;
2325  }
2326 
2327  pixtype = band->pixtype;
2328 
2329  if (
2330  x < 0 || x >= band->width ||
2331  y < 0 || y >= band->height
2332  ) {
2333  rterror("rt_band_set_pixel: Coordinates out of range");
2334  return ES_ERROR;
2335  }
2336 
2337  /* check that clamped value isn't clamped NODATA */
2338  if (band->hasnodata && pixtype != PT_64BF) {
2339  double newval;
2340  int corrected;
2341 
2342  rt_band_corrected_clamped_value(band, val, &newval, &corrected);
2343 
2344  if (corrected) {
2345 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
2346  rtwarn("Value for pixel %d x %d has been corrected as clamped value becomes NODATA", x, y);
2347 #endif
2348  val = newval;
2349 
2350  if (converted != NULL)
2351  *converted = 1;
2352  }
2353  }
2354 
2355  data = rt_band_get_data(band);
2356  offset = x + (y * band->width);
2357 
2358  switch (pixtype) {
2359  case PT_1BB: {
2360  data[offset] = rt_util_clamp_to_1BB(val);
2361  checkvalint = data[offset];
2362  break;
2363  }
2364  case PT_2BUI: {
2365  data[offset] = rt_util_clamp_to_2BUI(val);
2366  checkvalint = data[offset];
2367  break;
2368  }
2369  case PT_4BUI: {
2370  data[offset] = rt_util_clamp_to_4BUI(val);
2371  checkvalint = data[offset];
2372  break;
2373  }
2374  case PT_8BSI: {
2375  data[offset] = rt_util_clamp_to_8BSI(val);
2376  checkvalint = (int8_t) data[offset];
2377  break;
2378  }
2379  case PT_8BUI: {
2380  data[offset] = rt_util_clamp_to_8BUI(val);
2381  checkvalint = data[offset];
2382  break;
2383  }
2384  case PT_16BSI: {
2385  int16_t *ptr = (int16_t*) data; /* we assume correct alignment */
2386  ptr[offset] = rt_util_clamp_to_16BSI(val);
2387  checkvalint = (int16_t) ptr[offset];
2388  break;
2389  }
2390  case PT_16BUI: {
2391  uint16_t *ptr = (uint16_t*) data; /* we assume correct alignment */
2392  ptr[offset] = rt_util_clamp_to_16BUI(val);
2393  checkvalint = ptr[offset];
2394  break;
2395  }
2396  case PT_32BSI: {
2397  int32_t *ptr = (int32_t*) data; /* we assume correct alignment */
2398  ptr[offset] = rt_util_clamp_to_32BSI(val);
2399  checkvalint = (int32_t) ptr[offset];
2400  break;
2401  }
2402  case PT_32BUI: {
2403  uint32_t *ptr = (uint32_t*) data; /* we assume correct alignment */
2404  ptr[offset] = rt_util_clamp_to_32BUI(val);
2405  checkvaluint = ptr[offset];
2406  break;
2407  }
2408  case PT_32BF: {
2409  float *ptr = (float*) data; /* we assume correct alignment */
2410  ptr[offset] = rt_util_clamp_to_32F(val);
2411  checkvalfloat = ptr[offset];
2412  break;
2413  }
2414  case PT_64BF: {
2415  double *ptr = (double*) data; /* we assume correct alignment */
2416  ptr[offset] = val;
2417  checkvaldouble = ptr[offset];
2418  break;
2419  }
2420  default: {
2421  rterror("rt_band_set_pixel: Unknown pixeltype %d", pixtype);
2422  return ES_ERROR;
2423  }
2424  }
2425 
2426  /* If the stored value is not NODATA, reset the isnodata flag */
2427  if (!rt_band_clamped_value_is_nodata(band, val)) {
2428  RASTER_DEBUG(3, "Band has a value that is not NODATA. Setting isnodata to FALSE");
2429  band->isnodata = FALSE;
2430  }
2431 
2432  /* Overflow checking */
2434  val,
2435  checkvalint, checkvaluint,
2436  checkvalfloat, checkvaldouble,
2437  pixtype
2438  ) && converted != NULL) {
2439  *converted = 1;
2440  }
2441 
2442  return ES_NONE;
2443 }
2444 
2466  rt_band band,
2467  int x, int y,
2468  uint16_t len,
2469  void **vals, uint16_t *nvals
2470 ) {
2471  uint8_t *_vals = NULL;
2472  int pixsize = 0;
2473  uint8_t *data = NULL;
2474  uint32_t offset = 0;
2475  uint16_t _nvals = 0;
2476  int maxlen = 0;
2477  uint8_t *ptr = NULL;
2478 
2479  assert(NULL != band);
2480  assert(vals != NULL && nvals != NULL);
2481 
2482  /* initialize to no values */
2483  *nvals = 0;
2484 
2485  if (
2486  x < 0 || x >= band->width ||
2487  y < 0 || y >= band->height
2488  ) {
2489  rtwarn("Attempting to get pixel values with out of range raster coordinates: (%d, %d)", x, y);
2490  return ES_ERROR;
2491  }
2492 
2493  if (len < 1)
2494  return ES_NONE;
2495 
2496  data = rt_band_get_data(band);
2497  if (data == NULL) {
2498  rterror("rt_band_get_pixel_line: Cannot get band data");
2499  return ES_ERROR;
2500  }
2501 
2502  /* +1 for the nodata value */
2503  offset = x + (y * band->width);
2504  RASTER_DEBUGF(4, "offset = %d", offset);
2505 
2506  pixsize = rt_pixtype_size(band->pixtype);
2507  RASTER_DEBUGF(4, "pixsize = %d", pixsize);
2508 
2509  /* cap _nvals so that it doesn't overflow */
2510  _nvals = len;
2511  maxlen = band->width * band->height;
2512 
2513  if (((int) (offset + _nvals)) > maxlen) {
2514  _nvals = maxlen - offset;
2515  rtwarn("Limiting returning number values to %d", _nvals);
2516  }
2517  RASTER_DEBUGF(4, "_nvals = %d", _nvals);
2518 
2519  ptr = data + (offset * pixsize);
2520 
2521  _vals = rtalloc(_nvals * pixsize);
2522  if (_vals == NULL) {
2523  rterror("rt_band_get_pixel_line: Could not allocate memory for pixel values");
2524  return ES_ERROR;
2525  }
2526 
2527  /* copy pixels */
2528  memcpy(_vals, ptr, _nvals * pixsize);
2529 
2530  *vals = _vals;
2531  *nvals = _nvals;
2532 
2533  return ES_NONE;
2534 }
2535 
2550  rt_band band,
2551  int x, int y,
2552  double *value,
2553  int *nodata
2554 ) {
2555  rt_pixtype pixtype = PT_END;
2556  uint8_t* data = NULL;
2557  uint32_t offset = 0;
2558 
2559  assert(NULL != band);
2560  assert(NULL != value);
2561 
2562  /* set nodata to 0 */
2563  if (nodata != NULL)
2564  *nodata = 0;
2565 
2566  if (
2567  x < 0 || x >= band->width ||
2568  y < 0 || y >= band->height
2569  ) {
2570  rtwarn("Attempting to get pixel value with out of range raster coordinates: (%d, %d)", x, y);
2571  return ES_ERROR;
2572  }
2573 
2574  /* band is NODATA */
2575  if (band->isnodata) {
2576  RASTER_DEBUG(3, "Band's isnodata flag is TRUE. Returning NODATA value");
2577  *value = band->nodataval;
2578  if (nodata != NULL) *nodata = 1;
2579  return ES_NONE;
2580  }
2581 
2582  data = rt_band_get_data(band);
2583  if (data == NULL) {
2584  rterror("rt_band_get_pixel: Cannot get band data");
2585  return ES_ERROR;
2586  }
2587 
2588  /* +1 for the nodata value */
2589  offset = x + (y * band->width);
2590 
2591  pixtype = band->pixtype;
2592 
2593  switch (pixtype) {
2594  case PT_1BB:
2595 #ifdef OPTIMIZE_SPACE
2596  {
2597  int byteOffset = offset / 8;
2598  int bitOffset = offset % 8;
2599  data += byteOffset;
2600 
2601  /* Bit to set is bitOffset into data */
2602  *value = getBits(data, val, 1, bitOffset);
2603  break;
2604  }
2605 #endif
2606  case PT_2BUI:
2607 #ifdef OPTIMIZE_SPACE
2608  {
2609  int byteOffset = offset / 4;
2610  int bitOffset = offset % 4;
2611  data += byteOffset;
2612 
2613  /* Bits to set start at bitOffset into data */
2614  *value = getBits(data, val, 2, bitOffset);
2615  break;
2616  }
2617 #endif
2618  case PT_4BUI:
2619 #ifdef OPTIMIZE_SPACE
2620  {
2621  int byteOffset = offset / 2;
2622  int bitOffset = offset % 2;
2623  data += byteOffset;
2624 
2625  /* Bits to set start at bitOffset into data */
2626  *value = getBits(data, val, 2, bitOffset);
2627  break;
2628  }
2629 #endif
2630  case PT_8BSI: {
2631  int8_t val = data[offset];
2632  *value = val;
2633  break;
2634  }
2635  case PT_8BUI: {
2636  uint8_t val = data[offset];
2637  *value = val;
2638  break;
2639  }
2640  case PT_16BSI: {
2641  int16_t *ptr = (int16_t*) data; /* we assume correct alignment */
2642  *value = ptr[offset];
2643  break;
2644  }
2645  case PT_16BUI: {
2646  uint16_t *ptr = (uint16_t*) data; /* we assume correct alignment */
2647  *value = ptr[offset];
2648  break;
2649  }
2650  case PT_32BSI: {
2651  int32_t *ptr = (int32_t*) data; /* we assume correct alignment */
2652  *value = ptr[offset];
2653  break;
2654  }
2655  case PT_32BUI: {
2656  uint32_t *ptr = (uint32_t*) data; /* we assume correct alignment */
2657  *value = ptr[offset];
2658  break;
2659  }
2660  case PT_32BF: {
2661  float *ptr = (float*) data; /* we assume correct alignment */
2662  *value = ptr[offset];
2663  break;
2664  }
2665  case PT_64BF: {
2666  double *ptr = (double*) data; /* we assume correct alignment */
2667  *value = ptr[offset];
2668  break;
2669  }
2670  default: {
2671  rterror("rt_band_get_pixel: Unknown pixeltype %d", pixtype);
2672  return ES_ERROR;
2673  }
2674  }
2675 
2676  /* set NODATA flag */
2677  if (band->hasnodata && nodata != NULL) {
2678  if (rt_band_clamped_value_is_nodata(band, *value))
2679  *nodata = 1;
2680  }
2681 
2682  return ES_NONE;
2683 }
2684 
2703  rt_band band,
2704  int x, int y,
2705  uint16_t distancex, uint16_t distancey,
2706  int exclude_nodata_value,
2707  rt_pixel *npixels
2708 ) {
2709  rt_pixel npixel = NULL;
2710  int extent[4] = {0};
2711  int max_extent[4] = {0};
2712  int d0 = 0;
2713  int distance[2] = {0};
2714  uint32_t _d[2] = {0};
2715  uint32_t i = 0;
2716  uint32_t j = 0;
2717  uint32_t k = 0;
2718  int _max = 0;
2719  int _x = 0;
2720  int _y = 0;
2721  int *_min = NULL;
2722  double pixval = 0;
2723  double minval = 0;
2724  uint32_t count = 0;
2725  int isnodata = 0;
2726 
2727  int inextent = 0;
2728 
2729  assert(NULL != band);
2730  assert(NULL != npixels);
2731 
2732  RASTER_DEBUG(3, "Starting");
2733 
2734  /* process distance */
2735  distance[0] = distancex;
2736  distance[1] = distancey;
2737 
2738  /* no distance, means get nearest pixels and return */
2739  if (!distance[0] && !distance[1])
2740  d0 = 1;
2741 
2742  RASTER_DEBUGF(4, "Selected pixel: %d x %d", x, y);
2743  RASTER_DEBUGF(4, "Distances: %d x %d", distance[0], distance[1]);
2744 
2745  /* shortcuts if outside band extent */
2746  if (
2747  exclude_nodata_value && (
2748  (x < 0 || x > band->width) ||
2749  (y < 0 || y > band->height)
2750  )
2751  ) {
2752  /* no distances specified, jump to pixel close to extent */
2753  if (d0) {
2754  if (x < 0)
2755  x = -1;
2756  else if (x > band->width)
2757  x = band->width;
2758 
2759  if (y < 0)
2760  y = -1;
2761  else if (y > band->height)
2762  y = band->height;
2763 
2764  RASTER_DEBUGF(4, "Moved selected pixel: %d x %d", x, y);
2765  }
2766  /*
2767  distances specified
2768  if distances won't capture extent of band, return 0
2769  */
2770  else if (
2771  ((x < 0 && abs(x) > distance[0]) || (x - band->width >= distance[0])) ||
2772  ((y < 0 && abs(y) > distance[1]) || (y - band->height >= distance[1]))
2773  ) {
2774  RASTER_DEBUG(4, "No nearest pixels possible for provided pixel and distances");
2775  return 0;
2776  }
2777  }
2778 
2779  /* no NODATA, exclude is FALSE */
2780  if (!band->hasnodata)
2781  exclude_nodata_value = FALSE;
2782  /* band is NODATA and excluding NODATA */
2783  else if (exclude_nodata_value && band->isnodata) {
2784  RASTER_DEBUG(4, "No nearest pixels possible as band is NODATA and excluding NODATA values");
2785  return 0;
2786  }
2787 
2788  /* determine the maximum distance to prevent an infinite loop */
2789  if (d0) {
2790  int a, b;
2791 
2792  /* X axis */
2793  a = abs(x);
2794  b = abs(x - band->width);
2795 
2796  if (a > b)
2797  distance[0] = a;
2798  else
2799  distance[0] = b;
2800 
2801  /* Y axis */
2802  a = abs(y);
2803  b = abs(y - band->height);
2804  if (a > b)
2805  distance[1] = a;
2806  else
2807  distance[1] = b;
2808 
2809  RASTER_DEBUGF(4, "Maximum distances: %d x %d", distance[0], distance[1]);
2810  }
2811 
2812  /* minimum possible value for pixel type */
2813  minval = rt_pixtype_get_min_value(band->pixtype);
2814  RASTER_DEBUGF(4, "pixtype: %s", rt_pixtype_name(band->pixtype));
2815  RASTER_DEBUGF(4, "minval: %f", minval);
2816 
2817  /* set variables */
2818  count = 0;
2819  *npixels = NULL;
2820 
2821  /* maximum extent */
2822  max_extent[0] = x - distance[0]; /* min X */
2823  max_extent[1] = y - distance[1]; /* min Y */
2824  max_extent[2] = x + distance[0]; /* max X */
2825  max_extent[3] = y + distance[1]; /* max Y */
2826  RASTER_DEBUGF(4, "Maximum Extent: (%d, %d, %d, %d)",
2827  max_extent[0], max_extent[1], max_extent[2], max_extent[3]);
2828 
2829  _d[0] = 0;
2830  _d[1] = 0;
2831  do {
2832  _d[0]++;
2833  _d[1]++;
2834 
2835  extent[0] = x - _d[0]; /* min x */
2836  extent[1] = y - _d[1]; /* min y */
2837  extent[2] = x + _d[0]; /* max x */
2838  extent[3] = y + _d[1]; /* max y */
2839 
2840  RASTER_DEBUGF(4, "Processing distances: %d x %d", _d[0], _d[1]);
2841  RASTER_DEBUGF(4, "Extent: (%d, %d, %d, %d)",
2842  extent[0], extent[1], extent[2], extent[3]);
2843 
2844  for (i = 0; i < 2; i++) {
2845 
2846  /* by row */
2847  if (i < 1)
2848  _max = extent[2] - extent[0] + 1;
2849  /* by column */
2850  else
2851  _max = extent[3] - extent[1] + 1;
2852  _max = abs(_max);
2853 
2854  for (j = 0; j < 2; j++) {
2855  /* by row */
2856  if (i < 1) {
2857  _x = extent[0];
2858  _min = &_x;
2859 
2860  /* top row */
2861  if (j < 1)
2862  _y = extent[1];
2863  /* bottom row */
2864  else
2865  _y = extent[3];
2866  }
2867  /* by column */
2868  else {
2869  _y = extent[1] + 1;
2870  _min = &_y;
2871 
2872  /* left column */
2873  if (j < 1) {
2874  _x = extent[0];
2875  _max -= 2;
2876  }
2877  /* right column */
2878  else
2879  _x = extent[2];
2880  }
2881 
2882  RASTER_DEBUGF(4, "_min, _max: %d, %d", *_min, _max);
2883  for (k = 0; k < _max; k++) {
2884  /* check that _x and _y are not outside max extent */
2885  if (
2886  _x < max_extent[0] || _x > max_extent[2] ||
2887  _y < max_extent[1] || _y > max_extent[3]
2888  ) {
2889  (*_min)++;
2890  continue;
2891  }
2892 
2893  /* outside band extent, set to NODATA */
2894  if (
2895  (_x < 0 || _x >= band->width) ||
2896  (_y < 0 || _y >= band->height)
2897  ) {
2898  /* no NODATA, set to minimum possible value */
2899  if (!band->hasnodata)
2900  pixval = minval;
2901  /* has NODATA, use NODATA */
2902  else
2903  pixval = band->nodataval;
2904  RASTER_DEBUGF(4, "NODATA pixel outside band extent: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
2905  inextent = 0;
2906  isnodata = 1;
2907  }
2908  else {
2909  if (rt_band_get_pixel(
2910  band,
2911  _x, _y,
2912  &pixval,
2913  &isnodata
2914  ) != ES_NONE) {
2915  rterror("rt_band_get_nearest_pixel: Could not get pixel value");
2916  if (count) rtdealloc(*npixels);
2917  return -1;
2918  }
2919  RASTER_DEBUGF(4, "Pixel: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
2920  inextent = 1;
2921  }
2922 
2923  /* use pixval? */
2924  if (!exclude_nodata_value || (exclude_nodata_value && !isnodata)) {
2925  /* add pixel to result set */
2926  RASTER_DEBUGF(4, "Adding pixel to set of nearest pixels: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
2927  count++;
2928 
2929  if (*npixels == NULL)
2930  *npixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count);
2931  else
2932  *npixels = (rt_pixel) rtrealloc(*npixels, sizeof(struct rt_pixel_t) * count);
2933  if (*npixels == NULL) {
2934  rterror("rt_band_get_nearest_pixel: Could not allocate memory for nearest pixel(s)");
2935  return -1;
2936  }
2937 
2938  npixel = &((*npixels)[count - 1]);
2939  npixel->x = _x;
2940  npixel->y = _y;
2941  npixel->value = pixval;
2942 
2943  /* special case for when outside band extent */
2944  if (!inextent && !band->hasnodata)
2945  npixel->nodata = 1;
2946  else
2947  npixel->nodata = 0;
2948  }
2949 
2950  (*_min)++;
2951  }
2952  }
2953  }
2954 
2955  /* distance threshholds met */
2956  if (_d[0] >= distance[0] && _d[1] >= distance[1])
2957  break;
2958  else if (d0 && count)
2959  break;
2960  }
2961  while (1);
2962 
2963  RASTER_DEBUGF(3, "Nearest pixels in return: %d", count);
2964 
2965  return count;
2966 }
2967 
2979 int
2981  rt_band band, int exclude_nodata_value,
2982  double *searchset, int searchcount,
2983  rt_pixel *pixels
2984 ) {
2985  int x;
2986  int y;
2987  int i;
2988  double pixval;
2989  int err;
2990  int count = 0;
2991  int isnodata = 0;
2992  int isequal = 0;
2993 
2994  rt_pixel pixel = NULL;
2995 
2996  assert(NULL != band);
2997  assert(NULL != pixels);
2998  assert(NULL != searchset && searchcount > 0);
2999 
3000  if (!band->hasnodata)
3001  exclude_nodata_value = FALSE;
3002  /* band is NODATA and exclude_nodata_value = TRUE, nothing to search */
3003  else if (exclude_nodata_value && band->isnodata) {
3004  RASTER_DEBUG(4, "Pixels cannot be searched as band is NODATA and excluding NODATA values");
3005  return 0;
3006  }
3007 
3008  for (x = 0; x < band->width; x++) {
3009  for (y = 0; y < band->height; y++) {
3010  err = rt_band_get_pixel(band, x, y, &pixval, &isnodata);
3011  if (err != ES_NONE) {
3012  rterror("rt_band_get_pixel_of_value: Cannot get band pixel");
3013  return -1;
3014  }
3015  else if (exclude_nodata_value && isnodata)
3016  continue;
3017 
3018  for (i = 0; i < searchcount; i++) {
3019  if (rt_pixtype_compare_clamped_values(band->pixtype, searchset[i], pixval, &isequal) != ES_NONE) {
3020  continue;
3021  }
3022 
3023  if (FLT_NEQ(pixval, searchset[i]) || !isequal)
3024  continue;
3025 
3026  /* match found */
3027  count++;
3028  if (*pixels == NULL)
3029  *pixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count);
3030  else
3031  *pixels = (rt_pixel) rtrealloc(*pixels, sizeof(struct rt_pixel_t) * count);
3032  if (*pixels == NULL) {
3033  rterror("rt_band_get_pixel_of_value: Could not allocate memory for pixel(s)");
3034  return -1;
3035  }
3036 
3037  pixel = &((*pixels)[count - 1]);
3038  pixel->x = x;
3039  pixel->y = y;
3040  pixel->nodata = 0;
3041  pixel->value = pixval;
3042  }
3043  }
3044  }
3045 
3046  return count;
3047 }
3048 
3058 rt_band_get_nodata(rt_band band, double *nodata) {
3059  assert(NULL != band);
3060  assert(NULL != nodata);
3061 
3062  *nodata = band->nodataval;
3063 
3064  if (!band->hasnodata) {
3065  rterror("rt_band_get_nodata: Band has no NODATA value");
3066  return ES_ERROR;
3067  }
3068 
3069  return ES_NONE;
3070 }
3071 
3072 double
3074  assert(NULL != band);
3075 
3076  return rt_pixtype_get_min_value(band->pixtype);
3077 }
3078 
3079 int
3081  int i, j, err;
3082  double pxValue;
3083  int isnodata = 0;
3084 
3085  assert(NULL != band);
3086 
3087  /* Check if band has nodata value */
3088  if (!band->hasnodata) {
3089  RASTER_DEBUG(3, "Band has no NODATA value");
3090  band->isnodata = FALSE;
3091  return FALSE;
3092  }
3093 
3094  pxValue = band->nodataval;
3095 
3096  /* Check all pixels */
3097  for (i = 0; i < band->width; i++) {
3098  for (j = 0; j < band->height; j++) {
3099  err = rt_band_get_pixel(band, i, j, &pxValue, &isnodata);
3100  if (err != ES_NONE) {
3101  rterror("rt_band_check_is_nodata: Cannot get band pixel");
3102  return FALSE;
3103  }
3104  else if (!isnodata) {
3105  band->isnodata = FALSE;
3106  return FALSE;
3107  }
3108  }
3109  }
3110 
3111  band->isnodata = TRUE;
3112  return TRUE;
3113 }
3114 
3125 int
3127  int isequal = 0;
3128 
3129  assert(NULL != band);
3130 
3131  /* no NODATA, so never equal */
3132  if (!band->hasnodata)
3133  return 0;
3134 
3135  /* value is exactly NODATA */
3136  if (FLT_EQ(val, band->nodataval))
3137  return 2;
3138 
3139  /* ignore error from rt_pixtype_compare_clamped_values */
3141  band->pixtype,
3142  val, band->nodataval,
3143  &isequal
3144  );
3145 
3146  return isequal ? 1 : 0;
3147 }
3148 
3163  rt_band band,
3164  double val,
3165  double *newval, int *corrected
3166 ) {
3167  double minval = 0.;
3168 
3169  assert(NULL != band);
3170  assert(NULL != newval);
3171 
3172  if (corrected != NULL)
3173  *corrected = 0;
3174 
3175  /* no need to correct if clamped values IS NOT clamped NODATA */
3176  if (rt_band_clamped_value_is_nodata(band, val) != 1) {
3177  *newval = val;
3178  return ES_NONE;
3179  }
3180 
3181  minval = rt_pixtype_get_min_value(band->pixtype);
3182  *newval = val;
3183 
3184  switch (band->pixtype) {
3185  case PT_1BB:
3186  *newval = !band->nodataval;
3187  break;
3188  case PT_2BUI:
3189  if (rt_util_clamp_to_2BUI(val) == rt_util_clamp_to_2BUI(minval))
3190  (*newval)++;
3191  else
3192  (*newval)--;
3193  break;
3194  case PT_4BUI:
3195  if (rt_util_clamp_to_4BUI(val) == rt_util_clamp_to_4BUI(minval))
3196  (*newval)++;
3197  else
3198  (*newval)--;
3199  break;
3200  case PT_8BSI:
3201  if (rt_util_clamp_to_8BSI(val) == rt_util_clamp_to_8BSI(minval))
3202  (*newval)++;
3203  else
3204  (*newval)--;
3205  break;
3206  case PT_8BUI:
3207  if (rt_util_clamp_to_8BUI(val) == rt_util_clamp_to_8BUI(minval))
3208  (*newval)++;
3209  else
3210  (*newval)--;
3211  break;
3212  case PT_16BSI:
3213  if (rt_util_clamp_to_16BSI(val) == rt_util_clamp_to_16BSI(minval))
3214  (*newval)++;
3215  else
3216  (*newval)--;
3217  break;
3218  case PT_16BUI:
3219  if (rt_util_clamp_to_16BUI(val) == rt_util_clamp_to_16BUI(minval))
3220  (*newval)++;
3221  else
3222  (*newval)--;
3223  break;
3224  case PT_32BSI:
3225  if (rt_util_clamp_to_32BSI(val) == rt_util_clamp_to_32BSI(minval))
3226  (*newval)++;
3227  else
3228  (*newval)--;
3229  break;
3230  case PT_32BUI:
3231  if (rt_util_clamp_to_32BUI(val) == rt_util_clamp_to_32BUI(minval))
3232  (*newval)++;
3233  else
3234  (*newval)--;
3235  break;
3236  case PT_32BF:
3238  *newval += FLT_EPSILON;
3239  else
3240  *newval -= FLT_EPSILON;
3241  break;
3242  case PT_64BF:
3243  break;
3244  default:
3245  rterror("rt_band_corrected_clamped_value: Unknown pixeltype %d", band->pixtype);
3246  return ES_ERROR;
3247  }
3248 
3249  if (corrected != NULL)
3250  *corrected = 1;
3251 
3252  return ES_NONE;
3253 }
3254 
3270  rt_band band,
3271  int exclude_nodata_value, double sample, int inc_vals,
3272  uint64_t *cK, double *cM, double *cQ
3273 ) {
3274  uint32_t x = 0;
3275  uint32_t y = 0;
3276  uint32_t z = 0;
3277  uint32_t offset = 0;
3278  uint32_t diff = 0;
3279  int rtn;
3280  int hasnodata = FALSE;
3281  double nodata = 0;
3282  double *values = NULL;
3283  double value;
3284  int isnodata = 0;
3285  rt_bandstats stats = NULL;
3286 
3287  uint32_t do_sample = 0;
3288  uint32_t sample_size = 0;
3289  uint32_t sample_per = 0;
3290  uint32_t sample_int = 0;
3291  uint32_t i = 0;
3292  uint32_t j = 0;
3293  double sum = 0;
3294  uint32_t k = 0;
3295  double M = 0;
3296  double Q = 0;
3297 
3298 #if POSTGIS_DEBUG_LEVEL > 0
3299  clock_t start, stop;
3300  double elapsed = 0;
3301 #endif
3302 
3303  RASTER_DEBUG(3, "starting");
3304 #if POSTGIS_DEBUG_LEVEL > 0
3305  start = clock();
3306 #endif
3307 
3308  assert(NULL != band);
3309 
3310  /* band is empty (width < 1 || height < 1) */
3311  if (band->width < 1 || band->height < 1) {
3312  stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
3313  if (NULL == stats) {
3314  rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
3315  return NULL;
3316  }
3317 
3318  rtwarn("Band is empty as width and/or height is 0");
3319 
3320  stats->sample = 1;
3321  stats->sorted = 0;
3322  stats->values = NULL;
3323  stats->count = 0;
3324  stats->min = stats->max = 0;
3325  stats->sum = 0;
3326  stats->mean = 0;
3327  stats->stddev = -1;
3328 
3329  return stats;
3330  }
3331 
3332  hasnodata = rt_band_get_hasnodata_flag(band);
3333  if (hasnodata != FALSE)
3334  rt_band_get_nodata(band, &nodata);
3335  else
3336  exclude_nodata_value = 0;
3337 
3338  RASTER_DEBUGF(3, "nodata = %f", nodata);
3339  RASTER_DEBUGF(3, "hasnodata = %d", hasnodata);
3340  RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
3341 
3342  /* entire band is nodata */
3343  if (rt_band_get_isnodata_flag(band) != FALSE) {
3344  stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
3345  if (NULL == stats) {
3346  rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
3347  return NULL;
3348  }
3349 
3350  stats->sample = 1;
3351  stats->sorted = 0;
3352  stats->values = NULL;
3353 
3354  if (exclude_nodata_value) {
3355  rtwarn("All pixels of band have the NODATA value");
3356 
3357  stats->count = 0;
3358  stats->min = stats->max = 0;
3359  stats->sum = 0;
3360  stats->mean = 0;
3361  stats->stddev = -1;
3362  }
3363  else {
3364  stats->count = band->width * band->height;
3365  stats->min = stats->max = nodata;
3366  stats->sum = stats->count * nodata;
3367  stats->mean = nodata;
3368  stats->stddev = 0;
3369  }
3370 
3371  return stats;
3372  }
3373 
3374  /* clamp percentage */
3375  if (
3376  (sample < 0 || FLT_EQ(sample, 0.0)) ||
3377  (sample > 1 || FLT_EQ(sample, 1.0))
3378  ) {
3379  do_sample = 0;
3380  sample = 1;
3381  }
3382  else
3383  do_sample = 1;
3384  RASTER_DEBUGF(3, "do_sample = %d", do_sample);
3385 
3386  /* sample all pixels */
3387  if (!do_sample) {
3388  sample_size = band->width * band->height;
3389  sample_per = band->height;
3390  }
3391  /*
3392  randomly sample a percentage of available pixels
3393  sampling method is known as
3394  "systematic random sample without replacement"
3395  */
3396  else {
3397  sample_size = round((band->width * band->height) * sample);
3398  sample_per = round(sample_size / band->width);
3399  if (sample_per < 1)
3400  sample_per = 1;
3401  sample_int = round(band->height / sample_per);
3402  srand(time(NULL));
3403  }
3404 
3405  RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set"
3406  , sample_size, (band->width * band->height), sample_per);
3407 
3408  if (inc_vals) {
3409  values = rtalloc(sizeof(double) * sample_size);
3410  if (NULL == values) {
3411  rtwarn("Could not allocate memory for values");
3412  inc_vals = 0;
3413  }
3414  }
3415 
3416  /* initialize stats */
3417  stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
3418  if (NULL == stats) {
3419  rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
3420  return NULL;
3421  }
3422  stats->sample = sample;
3423  stats->count = 0;
3424  stats->sum = 0;
3425  stats->mean = 0;
3426  stats->stddev = -1;
3427  stats->min = stats->max = 0;
3428  stats->values = NULL;
3429  stats->sorted = 0;
3430 
3431  for (x = 0, j = 0, k = 0; x < band->width; x++) {
3432  y = -1;
3433  diff = 0;
3434 
3435  for (i = 0, z = 0; i < sample_per; i++) {
3436  if (!do_sample)
3437  y = i;
3438  else {
3439  offset = (rand() % sample_int) + 1;
3440  y += diff + offset;
3441  diff = sample_int - offset;
3442  }
3443  RASTER_DEBUGF(5, "(x, y, z) = (%d, %d, %d)", x, y, z);
3444  if (y >= band->height || z > sample_per) break;
3445 
3446  rtn = rt_band_get_pixel(band, x, y, &value, &isnodata);
3447 
3448  j++;
3449  if (rtn == ES_NONE && (!exclude_nodata_value || (exclude_nodata_value && !isnodata))) {
3450 
3451  /* inc_vals set, collect pixel values */
3452  if (inc_vals) values[k] = value;
3453 
3454  /* average */
3455  k++;
3456  sum += value;
3457 
3458  /*
3459  one-pass standard deviation
3460  http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
3461  */
3462  if (k == 1) {
3463  Q = 0;
3464  M = value;
3465  }
3466  else {
3467  Q += (((k - 1) * pow(value - M, 2)) / k);
3468  M += ((value - M ) / k);
3469  }
3470 
3471  /* coverage one-pass standard deviation */
3472  if (NULL != cK) {
3473  (*cK)++;
3474  if (*cK == 1) {
3475  *cQ = 0;
3476  *cM = value;
3477  }
3478  else {
3479  *cQ += (((*cK - 1) * pow(value - *cM, 2)) / *cK);
3480  *cM += ((value - *cM ) / *cK);
3481  }
3482  }
3483 
3484  /* min/max */
3485  if (stats->count < 1) {
3486  stats->count = 1;
3487  stats->min = stats->max = value;
3488  }
3489  else {
3490  if (value < stats->min)
3491  stats->min = value;
3492  if (value > stats->max)
3493  stats->max = value;
3494  }
3495 
3496  }
3497 
3498  z++;
3499  }
3500  }
3501 
3502  RASTER_DEBUG(3, "sampling complete");
3503 
3504  stats->count = k;
3505  if (k > 0) {
3506  if (inc_vals) {
3507  /* free unused memory */
3508  if (sample_size != k) {
3509  values = rtrealloc(values, k * sizeof(double));
3510  }
3511 
3512  stats->values = values;
3513  }
3514 
3515  stats->sum = sum;
3516  stats->mean = sum / k;
3517 
3518  /* standard deviation */
3519  if (!do_sample)
3520  stats->stddev = sqrt(Q / k);
3521  /* sample deviation */
3522  else {
3523  if (k < 2)
3524  stats->stddev = -1;
3525  else
3526  stats->stddev = sqrt(Q / (k - 1));
3527  }
3528  }
3529  /* inc_vals thus values allocated but not used */
3530  else if (inc_vals)
3531  rtdealloc(values);
3532 
3533  /* if do_sample is one */
3534  if (do_sample && k < 1)
3535  rtwarn("All sampled pixels of band have the NODATA value");
3536 
3537 #if POSTGIS_DEBUG_LEVEL > 0
3538  stop = clock();
3539  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
3540  RASTER_DEBUGF(3, "(time, count, mean, stddev, min, max) = (%0.4f, %d, %f, %f, %f, %f)",
3541  elapsed, stats->count, stats->mean, stats->stddev, stats->min, stats->max);
3542 #endif
3543 
3544  RASTER_DEBUG(3, "done");
3545  return stats;
3546 }
3547 
3568  rt_bandstats stats,
3569  int bin_count, double *bin_width, int bin_width_count,
3570  int right, double min, double max,
3571  uint32_t *rtn_count
3572 ) {
3573  rt_histogram bins = NULL;
3574  int init_width = 0;
3575  int i;
3576  int j;
3577  double tmp;
3578  double value;
3579  int sum = 0;
3580  double qmin;
3581  double qmax;
3582 
3583 #if POSTGIS_DEBUG_LEVEL > 0
3584  clock_t start, stop;
3585  double elapsed = 0;
3586 #endif
3587 
3588  RASTER_DEBUG(3, "starting");
3589 #if POSTGIS_DEBUG_LEVEL > 0
3590  start = clock();
3591 #endif
3592 
3593  assert(NULL != stats);
3594  assert(NULL != rtn_count);
3595 
3596  if (stats->count < 1 || NULL == stats->values) {
3597  rterror("rt_util_get_histogram: rt_bandstats object has no value");
3598  return NULL;
3599  }
3600 
3601  /* bin width must be positive numbers and not zero */
3602  if (NULL != bin_width && bin_width_count > 0) {
3603  for (i = 0; i < bin_width_count; i++) {
3604  if (bin_width[i] < 0 || FLT_EQ(bin_width[i], 0.0)) {
3605  rterror("rt_util_get_histogram: bin_width element is less than or equal to zero");
3606  return NULL;
3607  }
3608  }
3609  }
3610 
3611  /* ignore min and max parameters */
3612  if (FLT_EQ(max, min)) {
3613  qmin = stats->min;
3614  qmax = stats->max;
3615  }
3616  else {
3617  qmin = min;
3618  qmax = max;
3619  if (qmin > qmax) {
3620  qmin = max;
3621  qmax = min;
3622  }
3623  }
3624 
3625  /* # of bins not provided */
3626  if (bin_count <= 0) {
3627  /*
3628  determine # of bins
3629  http://en.wikipedia.org/wiki/Histogram
3630 
3631  all computed bins are assumed to have equal width
3632  */
3633  /* Square-root choice for stats->count < 30 */
3634  if (stats->count < 30)
3635  bin_count = ceil(sqrt(stats->count));
3636  /* Sturges' formula for stats->count >= 30 */
3637  else
3638  bin_count = ceil(log2((double) stats->count) + 1.);
3639 
3640  /* bin_width_count provided and bin_width has value */
3641  if (bin_width_count > 0 && NULL != bin_width) {
3642  /* user has defined something specific */
3643  if (bin_width_count > bin_count)
3644  bin_count = bin_width_count;
3645  else if (bin_width_count > 1) {
3646  tmp = 0;
3647  for (i = 0; i < bin_width_count; i++) tmp += bin_width[i];
3648  bin_count = ceil((qmax - qmin) / tmp) * bin_width_count;
3649  }
3650  else
3651  bin_count = ceil((qmax - qmin) / bin_width[0]);
3652  }
3653  /* set bin width count to zero so that one can be calculated */
3654  else {
3655  bin_width_count = 0;
3656  }
3657  }
3658 
3659  /* min and max the same */
3660  if (FLT_EQ(qmax, qmin)) bin_count = 1;
3661 
3662  RASTER_DEBUGF(3, "bin_count = %d", bin_count);
3663 
3664  /* bin count = 1, all values are in one bin */
3665  if (bin_count < 2) {
3666  bins = rtalloc(sizeof(struct rt_histogram_t));
3667  if (NULL == bins) {
3668  rterror("rt_util_get_histogram: Could not allocate memory for histogram");
3669  return NULL;
3670  }
3671 
3672  bins->count = stats->count;
3673  bins->percent = -1;
3674  bins->min = qmin;
3675  bins->max = qmax;
3676  bins->inc_min = bins->inc_max = 1;
3677 
3678  *rtn_count = bin_count;
3679  return bins;
3680  }
3681 
3682  /* establish bin width */
3683  if (bin_width_count == 0) {
3684  bin_width_count = 1;
3685 
3686  /* bin_width unallocated */
3687  if (NULL == bin_width) {
3688  bin_width = rtalloc(sizeof(double));
3689  if (NULL == bin_width) {
3690  rterror("rt_util_get_histogram: Could not allocate memory for bin widths");
3691  return NULL;
3692  }
3693  init_width = 1;
3694  }
3695 
3696  bin_width[0] = (qmax - qmin) / bin_count;
3697  }
3698 
3699  /* initialize bins */
3700  bins = rtalloc(bin_count * sizeof(struct rt_histogram_t));
3701  if (NULL == bins) {
3702  rterror("rt_util_get_histogram: Could not allocate memory for histogram");
3703  if (init_width) rtdealloc(bin_width);
3704  return NULL;
3705  }
3706  if (!right)
3707  tmp = qmin;
3708  else
3709  tmp = qmax;
3710  for (i = 0; i < bin_count;) {
3711  for (j = 0; j < bin_width_count; j++) {
3712  bins[i].count = 0;
3713  bins->percent = -1;
3714 
3715  if (!right) {
3716  bins[i].min = tmp;
3717  tmp += bin_width[j];
3718  bins[i].max = tmp;
3719 
3720  bins[i].inc_min = 1;
3721  bins[i].inc_max = 0;
3722  }
3723  else {
3724  bins[i].max = tmp;
3725  tmp -= bin_width[j];
3726  bins[i].min = tmp;
3727 
3728  bins[i].inc_min = 0;
3729  bins[i].inc_max = 1;
3730  }
3731 
3732  i++;
3733  }
3734  }
3735  if (!right) {
3736  bins[bin_count - 1].inc_max = 1;
3737 
3738  /* align last bin to the max value */
3739  if (bins[bin_count - 1].max < qmax)
3740  bins[bin_count - 1].max = qmax;
3741  }
3742  else {
3743  bins[bin_count - 1].inc_min = 1;
3744 
3745  /* align first bin to the min value */
3746  if (bins[bin_count - 1].min > qmin)
3747  bins[bin_count - 1].min = qmin;
3748  }
3749 
3750  /* process the values */
3751  for (i = 0; i < stats->count; i++) {
3752  value = stats->values[i];
3753 
3754  /* default, [a, b) */
3755  if (!right) {
3756  for (j = 0; j < bin_count; j++) {
3757  if (
3758  (!bins[j].inc_max && value < bins[j].max) || (
3759  bins[j].inc_max && (
3760  (value < bins[j].max) ||
3761  FLT_EQ(value, bins[j].max)
3762  )
3763  )
3764  ) {
3765  bins[j].count++;
3766  sum++;
3767  break;
3768  }
3769  }
3770  }
3771  /* (a, b] */
3772  else {
3773  for (j = 0; j < bin_count; j++) {
3774  if (
3775  (!bins[j].inc_min && value > bins[j].min) || (
3776  bins[j].inc_min && (
3777  (value > bins[j].min) ||
3778  FLT_EQ(value, bins[j].min)
3779  )
3780  )
3781  ) {
3782  bins[j].count++;
3783  sum++;
3784  break;
3785  }
3786  }
3787  }
3788  }
3789 
3790  for (i = 0; i < bin_count; i++) {
3791  bins[i].percent = ((double) bins[i].count) / sum;
3792  }
3793 
3794 #if POSTGIS_DEBUG_LEVEL > 0
3795  stop = clock();
3796  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
3797  RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
3798 
3799  for (j = 0; j < bin_count; j++) {
3800  RASTER_DEBUGF(5, "(min, max, inc_min, inc_max, count, sum, percent) = (%f, %f, %d, %d, %d, %d, %f)",
3801  bins[j].min, bins[j].max, bins[j].inc_min, bins[j].inc_max, bins[j].count, sum, bins[j].percent);
3802  }
3803 #endif
3804 
3805  if (init_width) rtdealloc(bin_width);
3806  *rtn_count = bin_count;
3807  RASTER_DEBUG(3, "done");
3808  return bins;
3809 }
3810 
3824  rt_bandstats stats,
3825  double *quantiles, int quantiles_count,
3826  uint32_t *rtn_count
3827 ) {
3828  rt_quantile rtn;
3829  int init_quantiles = 0;
3830  int i = 0;
3831  double h;
3832  int hl;
3833 
3834 #if POSTGIS_DEBUG_LEVEL > 0
3835  clock_t start, stop;
3836  double elapsed = 0;
3837 #endif
3838 
3839  RASTER_DEBUG(3, "starting");
3840 #if POSTGIS_DEBUG_LEVEL > 0
3841  start = clock();
3842 #endif
3843 
3844  assert(NULL != stats);
3845  assert(NULL != rtn_count);
3846 
3847  if (stats->count < 1 || NULL == stats->values) {
3848  rterror("rt_band_get_quantiles: rt_bandstats object has no value");
3849  return NULL;
3850  }
3851 
3852  /* quantiles not provided */
3853  if (NULL == quantiles) {
3854  /* quantile count not specified, default to quartiles */
3855  if (quantiles_count < 2)
3856  quantiles_count = 5;
3857 
3858  quantiles = rtalloc(sizeof(double) * quantiles_count);
3859  init_quantiles = 1;
3860  if (NULL == quantiles) {
3861  rterror("rt_band_get_quantiles: Could not allocate memory for quantile input");
3862  return NULL;
3863  }
3864 
3865  quantiles_count--;
3866  for (i = 0; i <= quantiles_count; i++)
3867  quantiles[i] = ((double) i) / quantiles_count;
3868  quantiles_count++;
3869  }
3870 
3871  /* check quantiles */
3872  for (i = 0; i < quantiles_count; i++) {
3873  if (quantiles[i] < 0. || quantiles[i] > 1.) {
3874  rterror("rt_band_get_quantiles: Quantile value not between 0 and 1");
3875  if (init_quantiles) rtdealloc(quantiles);
3876  return NULL;
3877  }
3878  }
3879  quicksort(quantiles, quantiles + quantiles_count - 1);
3880 
3881  /* initialize rt_quantile */
3882  rtn = rtalloc(sizeof(struct rt_quantile_t) * quantiles_count);
3883  if (NULL == rtn) {
3884  rterror("rt_band_get_quantiles: Could not allocate memory for quantile output");
3885  if (init_quantiles) rtdealloc(quantiles);
3886  return NULL;
3887  }
3888 
3889  /* sort values */
3890  if (!stats->sorted) {
3891  quicksort(stats->values, stats->values + stats->count - 1);
3892  stats->sorted = 1;
3893  }
3894 
3895  /*
3896  make quantiles
3897 
3898  formula is that used in R (method 7) and Excel from
3899  http://en.wikipedia.org/wiki/Quantile
3900  */
3901  for (i = 0; i < quantiles_count; i++) {
3902  rtn[i].quantile = quantiles[i];
3903 
3904  h = ((stats->count - 1.) * quantiles[i]) + 1.;
3905  hl = floor(h);
3906 
3907  /* h greater than hl, do full equation */
3908  if (h > hl)
3909  rtn[i].value = stats->values[hl - 1] + ((h - hl) * (stats->values[hl] - stats->values[hl - 1]));
3910  /* shortcut as second part of equation is zero */
3911  else
3912  rtn[i].value = stats->values[hl - 1];
3913  }
3914 
3915 #if POSTGIS_DEBUG_LEVEL > 0
3916  stop = clock();
3917  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
3918  RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
3919 #endif
3920 
3921  *rtn_count = quantiles_count;
3922  if (init_quantiles) rtdealloc(quantiles);
3923  RASTER_DEBUG(3, "done");
3924  return rtn;
3925 }
3926 
3928  struct quantile_llist_element *element,
3929  double needle
3930 ) {
3931  if (NULL == element)
3932  return NULL;
3933  else if (FLT_NEQ(needle, element->value)) {
3934  if (NULL != element->next)
3935  return quantile_llist_search(element->next, needle);
3936  else
3937  return NULL;
3938  }
3939  else
3940  return element;
3941 }
3942 
3944  struct quantile_llist_element *element,
3945  double value,
3946  uint32_t *idx
3947 ) {
3948  struct quantile_llist_element *qle = NULL;
3949 
3950  if (NULL == element) {
3951  qle = rtalloc(sizeof(struct quantile_llist_element));
3952  RASTER_DEBUGF(4, "qle @ %p is only element in list", qle);
3953  if (NULL == qle) return NULL;
3954 
3955  qle->value = value;
3956  qle->count = 1;
3957 
3958  qle->prev = NULL;
3959  qle->next = NULL;
3960 
3961  if (NULL != idx) *idx = 0;
3962  return qle;
3963  }
3964  else if (value > element->value) {
3965  if (NULL != idx) *idx += 1;
3966  if (NULL != element->next)
3967  return quantile_llist_insert(element->next, value, idx);
3968  /* insert as last element in list */
3969  else {
3970  qle = rtalloc(sizeof(struct quantile_llist_element));
3971  RASTER_DEBUGF(4, "insert qle @ %p as last element", qle);
3972  if (NULL == qle) return NULL;
3973 
3974  qle->value = value;
3975  qle->count = 1;
3976 
3977  qle->prev = element;
3978  qle->next = NULL;
3979  element->next = qle;
3980 
3981  return qle;
3982  }
3983  }
3984  /* insert before current element */
3985  else {
3986  qle = rtalloc(sizeof(struct quantile_llist_element));
3987  RASTER_DEBUGF(4, "insert qle @ %p before current element", qle);
3988  if (NULL == qle) return NULL;
3989 
3990  qle->value = value;
3991  qle->count = 1;
3992 
3993  if (NULL != element->prev) element->prev->next = qle;
3994  qle->next = element;
3995  qle->prev = element->prev;
3996  element->prev = qle;
3997 
3998  return qle;
3999  }
4000 }
4001 
4002 static int quantile_llist_delete(struct quantile_llist_element *element) {
4003  if (NULL == element) return 0;
4004 
4005  /* beginning of list */
4006  if (NULL == element->prev && NULL != element->next) {
4007  element->next->prev = NULL;
4008  }
4009  /* end of list */
4010  else if (NULL != element->prev && NULL == element->next) {
4011  element->prev->next = NULL;
4012  }
4013  /* within list */
4014  else if (NULL != element->prev && NULL != element->next) {
4015  element->prev->next = element->next;
4016  element->next->prev = element->prev;
4017  }
4018 
4019  RASTER_DEBUGF(4, "qle @ %p destroyed", element);
4020  rtdealloc(element);
4021 
4022  return 1;
4023 }
4024 
4025 int quantile_llist_destroy(struct quantile_llist **list, uint32_t list_count) {
4026  struct quantile_llist_element *element = NULL;
4027  uint32_t i;
4028 
4029  if (NULL == *list) return 0;
4030 
4031  for (i = 0; i < list_count; i++) {
4032  element = (*list)[i].head;
4033  while (NULL != element->next) {
4034  quantile_llist_delete(element->next);
4035  }
4036  quantile_llist_delete(element);
4037 
4038  rtdealloc((*list)[i].index);
4039  }
4040 
4041  rtdealloc(*list);
4042  return 1;
4043 }
4044 
4045 static void quantile_llist_index_update(struct quantile_llist *qll, struct quantile_llist_element *qle, uint32_t idx) {
4046  uint32_t anchor = (uint32_t) floor(idx / 100);
4047 
4048  if (qll->tail == qle) return;
4049 
4050  if (
4051  (anchor != 0) && (
4052  NULL == qll->index[anchor].element ||
4053  idx <= qll->index[anchor].index
4054  )
4055  ) {
4056  qll->index[anchor].index = idx;
4057  qll->index[anchor].element = qle;
4058  }
4059 
4060  if (anchor != 0 && NULL == qll->index[0].element) {
4061  qll->index[0].index = 0;
4062  qll->index[0].element = qll->head;
4063  }
4064 }
4065 
4067  uint32_t i = 0;
4068 
4069  for (i = 0; i < qll->index_max; i++) {
4070  if (
4071  NULL == qll->index[i].element ||
4072  (qll->index[i].element) != qle
4073  ) {
4074  continue;
4075  }
4076 
4077  RASTER_DEBUGF(5, "deleting index: %d => %f", i, qle->value);
4078  qll->index[i].index = -1;
4079  qll->index[i].element = NULL;
4080  }
4081 }
4082 
4084  struct quantile_llist *qll,
4085  double value,
4086  uint32_t *index
4087 ) {
4088  uint32_t i = 0, j = 0;
4089  RASTER_DEBUGF(5, "searching index for value %f", value);
4090 
4091  for (i = 0; i < qll->index_max; i++) {
4092  if (NULL == qll->index[i].element) {
4093  if (i < 1) break;
4094  continue;
4095  }
4096  if (value > (qll->index[i]).element->value) continue;
4097 
4098  if (FLT_EQ(value, qll->index[i].element->value)) {
4099  RASTER_DEBUGF(5, "using index value at %d = %f", i, qll->index[i].element->value);
4100  *index = i * 100;
4101  return qll->index[i].element;
4102  }
4103  else if (i > 0) {
4104  for (j = 1; j < i; j++) {
4105  if (NULL != qll->index[i - j].element) {
4106  RASTER_DEBUGF(5, "using index value at %d = %f", i - j, qll->index[i - j].element->value);
4107  *index = (i - j) * 100;
4108  return qll->index[i - j].element;
4109  }
4110  }
4111  }
4112  }
4113 
4114  *index = 0;
4115  return qll->head;
4116 }
4117 
4118 static void quantile_llist_index_reset(struct quantile_llist *qll) {
4119  uint32_t i = 0;
4120 
4121  RASTER_DEBUG(5, "resetting index");
4122  for (i = 0; i < qll->index_max; i++) {
4123  qll->index[i].index = -1;
4124  qll->index[i].element = NULL;
4125  }
4126 }
4127 
4157  rt_band band,
4158  int exclude_nodata_value, double sample,
4159  uint64_t cov_count,
4160  struct quantile_llist **qlls, uint32_t *qlls_count,
4161  double *quantiles, int quantiles_count,
4162  uint32_t *rtn_count
4163 ) {
4164  rt_quantile rtn = NULL;
4165  int init_quantiles = 0;
4166 
4167  struct quantile_llist *qll = NULL;
4168  struct quantile_llist_element *qle = NULL;
4169  struct quantile_llist_element *qls = NULL;
4170  const uint32_t MAX_VALUES = 750;
4171 
4172  uint8_t *data = NULL;
4173  double value;
4174  int isnodata = 0;
4175 
4176  uint32_t a = 0;
4177  uint32_t i = 0;
4178  uint32_t j = 0;
4179  uint32_t k = 0;
4180  uint32_t x = 0;
4181  uint32_t y = 0;
4182  uint32_t z = 0;
4183  uint32_t idx = 0;
4184  uint32_t offset = 0;
4185  uint32_t diff = 0;
4186  uint8_t exists = 0;
4187 
4188  uint32_t do_sample = 0;
4189  uint32_t sample_size = 0;
4190  uint32_t sample_per = 0;
4191  uint32_t sample_int = 0;
4192  int status;
4193 
4194  RASTER_DEBUG(3, "starting");
4195 
4196  assert(NULL != band);
4197  assert(cov_count > 1);
4198  assert(NULL != rtn_count);
4199  RASTER_DEBUGF(3, "cov_count = %d", cov_count);
4200 
4201  data = rt_band_get_data(band);
4202  if (data == NULL) {
4203  rterror("rt_band_get_summary_stats: Cannot get band data");
4204  return NULL;
4205  }
4206 
4207  if (!rt_band_get_hasnodata_flag(band))
4208  exclude_nodata_value = 0;
4209  RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
4210 
4211  /* quantile_llist not provided */
4212  if (NULL == *qlls) {
4213  /* quantiles not provided */
4214  if (NULL == quantiles) {
4215  /* quantile count not specified, default to quartiles */
4216  if (quantiles_count < 2)
4217  quantiles_count = 5;
4218 
4219  quantiles = rtalloc(sizeof(double) * quantiles_count);
4220  init_quantiles = 1;
4221  if (NULL == quantiles) {
4222  rterror("rt_band_get_quantiles_stream: Could not allocate memory for quantile input");
4223  return NULL;
4224  }
4225 
4226  quantiles_count--;
4227  for (i = 0; i <= quantiles_count; i++)
4228  quantiles[i] = ((double) i) / quantiles_count;
4229  quantiles_count++;
4230  }
4231 
4232  /* check quantiles */
4233  for (i = 0; i < quantiles_count; i++) {
4234  if (quantiles[i] < 0. || quantiles[i] > 1.) {
4235  rterror("rt_band_get_quantiles_stream: Quantile value not between 0 and 1");
4236  if (init_quantiles) rtdealloc(quantiles);
4237  return NULL;
4238  }
4239  }
4240  quicksort(quantiles, quantiles + quantiles_count - 1);
4241 
4242  /* initialize linked-list set */
4243  *qlls_count = quantiles_count * 2;
4244  RASTER_DEBUGF(4, "qlls_count = %d", *qlls_count);
4245  *qlls = rtalloc(sizeof(struct quantile_llist) * *qlls_count);
4246  if (NULL == *qlls) {
4247  rterror("rt_band_get_quantiles_stream: Could not allocate memory for quantile output");
4248  if (init_quantiles) rtdealloc(quantiles);
4249  return NULL;
4250  }
4251 
4252  j = (uint32_t) floor(MAX_VALUES / 100.) + 1;
4253  for (i = 0; i < *qlls_count; i++) {
4254  qll = &((*qlls)[i]);
4255  qll->quantile = quantiles[(i * quantiles_count) / *qlls_count];
4256  qll->count = 0;
4257  qll->sum1 = 0;
4258  qll->sum2 = 0;
4259  qll->head = NULL;
4260  qll->tail = NULL;
4261 
4262  /* initialize index */
4263  qll->index = rtalloc(sizeof(struct quantile_llist_index) * j);
4264  if (NULL == qll->index) {
4265  rterror("rt_band_get_quantiles_stream: Could not allocate memory for quantile output");
4266  if (init_quantiles) rtdealloc(quantiles);
4267  return NULL;
4268  }
4269  qll->index_max = j;
4271 
4272  /* AL-GEQ */
4273  if (!(i % 2)) {
4274  qll->algeq = 1;
4275  qll->tau = (uint64_t) ROUND(cov_count - (cov_count * qll->quantile), 0);
4276  if (qll->tau < 1) qll->tau = 1;
4277  }
4278  /* AL-GT */
4279  else {
4280  qll->algeq = 0;
4281  qll->tau = cov_count - (*qlls)[i - 1].tau + 1;
4282  }
4283 
4284  RASTER_DEBUGF(4, "qll init: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
4285  qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
4286  RASTER_DEBUGF(4, "qll init: (head, tail) = (%p, %p)", qll->head, qll->tail);
4287  }
4288 
4289  if (init_quantiles) rtdealloc(quantiles);
4290  }
4291 
4292  /* clamp percentage */
4293  if (
4294  (sample < 0 || FLT_EQ(sample, 0.0)) ||
4295  (sample > 1 || FLT_EQ(sample, 1.0))
4296  ) {
4297  do_sample = 0;
4298  sample = 1;
4299  }
4300  else
4301  do_sample = 1;
4302  RASTER_DEBUGF(3, "do_sample = %d", do_sample);
4303 
4304  /* sample all pixels */
4305  if (!do_sample) {
4306  sample_size = band->width * band->height;
4307  sample_per = band->height;
4308  }
4309  /*
4310  randomly sample a percentage of available pixels
4311  sampling method is known as
4312  "systematic random sample without replacement"
4313  */
4314  else {
4315  sample_size = round((band->width * band->height) * sample);
4316  sample_per = round(sample_size / band->width);
4317  sample_int = round(band->height / sample_per);
4318  srand(time(NULL));
4319  }
4320  RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set"
4321  , sample_size, (band->width * band->height), sample_per);
4322 
4323  for (x = 0, j = 0, k = 0; x < band->width; x++) {
4324  y = -1;
4325  diff = 0;
4326 
4327  /* exclude_nodata_value = TRUE and band is NODATA */
4328  if (exclude_nodata_value && rt_band_get_isnodata_flag(band)) {
4329  RASTER_DEBUG(3, "Skipping quantile calcuation as band is NODATA");
4330  break;
4331  }
4332 
4333  for (i = 0, z = 0; i < sample_per; i++) {
4334  if (do_sample != 1)
4335  y = i;
4336  else {
4337  offset = (rand() % sample_int) + 1;
4338  y += diff + offset;
4339  diff = sample_int - offset;
4340  }
4341  RASTER_DEBUGF(5, "(x, y, z) = (%d, %d, %d)", x, y, z);
4342  if (y >= band->height || z > sample_per) break;
4343 
4344  status = rt_band_get_pixel(band, x, y, &value, &isnodata);
4345 
4346  j++;
4347  if (status == ES_NONE && (!exclude_nodata_value || (exclude_nodata_value && !isnodata))) {
4348 
4349  /* process each quantile */
4350  for (a = 0; a < *qlls_count; a++) {
4351  qll = &((*qlls)[a]);
4352  qls = NULL;
4353  RASTER_DEBUGF(4, "%d of %d (%f)", a + 1, *qlls_count, qll->quantile);
4354  RASTER_DEBUGF(5, "qll before: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
4355  qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
4356  RASTER_DEBUGF(5, "qll before: (head, tail) = (%p, %p)", qll->head, qll->tail);
4357 
4358  /* OPTIMIZATION: shortcuts for quantiles of zero or one */
4359  if (FLT_EQ(qll->quantile, 0.)) {
4360  if (NULL != qll->head) {
4361  if (value < qll->head->value)
4362  qll->head->value = value;
4363  }
4364  else {
4365  qle = quantile_llist_insert(qll->head, value, NULL);
4366  qll->head = qle;
4367  qll->tail = qle;
4368  qll->count = 1;
4369  }
4370 
4371  RASTER_DEBUGF(4, "quantile shortcut for %f\n\n", qll->quantile);
4372  continue;
4373  }
4374  else if (FLT_EQ(qll->quantile, 1.)) {
4375  if (NULL != qll->head) {
4376  if (value > qll->head->value)
4377  qll->head->value = value;
4378  }
4379  else {
4380  qle = quantile_llist_insert(qll->head, value, NULL);
4381  qll->head = qle;
4382  qll->tail = qle;
4383  qll->count = 1;
4384  }
4385 
4386  RASTER_DEBUGF(4, "quantile shortcut for %f\n\n", qll->quantile);
4387  continue;
4388  }
4389 
4390  /* value exists in list */
4391  /* OPTIMIZATION: check to see if value exceeds last value */
4392  if (NULL != qll->tail && value > qll->tail->value)
4393  qle = NULL;
4394  /* OPTIMIZATION: check to see if value equals last value */
4395  else if (NULL != qll->tail && FLT_EQ(value, qll->tail->value))
4396  qle = qll->tail;
4397  /* OPTIMIZATION: use index if possible */
4398  else {
4399  qls = quantile_llist_index_search(qll, value, &idx);
4400  qle = quantile_llist_search(qls, value);
4401  }
4402 
4403  /* value found */
4404  if (NULL != qle) {
4405  RASTER_DEBUGF(4, "%f found in list", value);
4406  RASTER_DEBUGF(5, "qle before: (value, count) = (%f, %d)", qle->value, qle->count);
4407 
4408  qle->count++;
4409  qll->sum1++;
4410 
4411  if (qll->algeq)
4412  qll->sum2 = qll->sum1 - qll->head->count;
4413  else
4414  qll->sum2 = qll->sum1 - qll->tail->count;
4415 
4416  RASTER_DEBUGF(4, "qle after: (value, count) = (%f, %d)", qle->value, qle->count);
4417  }
4418  /* can still add element */
4419  else if (qll->count < MAX_VALUES) {
4420  RASTER_DEBUGF(4, "Adding %f to list", value);
4421 
4422  /* insert */
4423  /* OPTIMIZATION: check to see if value exceeds last value */
4424  if (NULL != qll->tail && (value > qll->tail->value || FLT_EQ(value, qll->tail->value))) {
4425  idx = qll->count;
4426  qle = quantile_llist_insert(qll->tail, value, &idx);
4427  }
4428  /* OPTIMIZATION: use index if possible */
4429  else
4430  qle = quantile_llist_insert(qls, value, &idx);
4431  if (NULL == qle) return NULL;
4432  RASTER_DEBUGF(5, "value added at index: %d => %f", idx, value);
4433  qll->count++;
4434  qll->sum1++;
4435 
4436  /* first element */
4437  if (NULL == qle->prev)
4438  qll->head = qle;
4439  /* last element */
4440  if (NULL == qle->next)
4441  qll->tail = qle;
4442 
4443  if (qll->algeq)
4444  qll->sum2 = qll->sum1 - qll->head->count;
4445  else
4446  qll->sum2 = qll->sum1 - qll->tail->count;
4447 
4448  /* index is only needed if there are at least 100 values */
4449  quantile_llist_index_update(qll, qle, idx);
4450 
4451  RASTER_DEBUGF(5, "qle, prev, next, head, tail = %p, %p, %p, %p, %p", qle, qle->prev, qle->next, qll->head, qll->tail);
4452  }
4453  /* AL-GEQ */
4454  else if (qll->algeq) {
4455  RASTER_DEBUGF(4, "value, head->value = %f, %f", value, qll->head->value);
4456 
4457  if (value < qll->head->value) {
4458  /* ignore value if test isn't true */
4459  if (qll->sum1 >= qll->tau) {
4460  RASTER_DEBUGF(4, "skipping %f", value);
4461  }
4462  else {
4463 
4464  /* delete last element */
4465  RASTER_DEBUGF(4, "deleting %f from list", qll->tail->value);
4466  qle = qll->tail->prev;
4467  RASTER_DEBUGF(5, "to-be tail is %f with count %d", qle->value, qle->count);
4468  qle->count += qll->tail->count;
4469  quantile_llist_index_delete(qll, qll->tail);
4471  qll->tail = qle;
4472  qll->count--;
4473  RASTER_DEBUGF(5, "tail is %f with count %d", qll->tail->value, qll->tail->count);
4474 
4475  /* insert value */
4476  RASTER_DEBUGF(4, "adding %f to list", value);
4477  /* OPTIMIZATION: check to see if value exceeds last value */
4478  if (NULL != qll->tail && (value > qll->tail->value || FLT_EQ(value, qll->tail->value))) {
4479  idx = qll->count;
4480  qle = quantile_llist_insert(qll->tail, value, &idx);
4481  }
4482  /* OPTIMIZATION: use index if possible */
4483  else {
4484  qls = quantile_llist_index_search(qll, value, &idx);
4485  qle = quantile_llist_insert(qls, value, &idx);
4486  }
4487  if (NULL == qle) return NULL;
4488  RASTER_DEBUGF(5, "value added at index: %d => %f", idx, value);
4489  qll->count++;
4490  qll->sum1++;
4491 
4492  /* first element */
4493  if (NULL == qle->prev)
4494  qll->head = qle;
4495  /* last element */
4496  if (NULL == qle->next)
4497  qll->tail = qle;
4498 
4499  qll->sum2 = qll->sum1 - qll->head->count;
4500 
4501  quantile_llist_index_update(qll, qle, idx);
4502 
4503  RASTER_DEBUGF(5, "qle, head, tail = %p, %p, %p", qle, qll->head, qll->tail);
4504 
4505  }
4506  }
4507  else {
4508  qle = qll->tail;
4509  while (NULL != qle) {
4510  if (qle->value < value) {
4511  qle->count++;
4512  qll->sum1++;
4513  qll->sum2 = qll->sum1 - qll->head->count;
4514  RASTER_DEBUGF(4, "incremented count of %f by 1 to %d", qle->value, qle->count);
4515  break;
4516  }
4517 
4518  qle = qle->prev;
4519  }
4520  }
4521  }
4522  /* AL-GT */
4523  else {
4524  RASTER_DEBUGF(4, "value, tail->value = %f, %f", value, qll->tail->value);
4525 
4526  if (value > qll->tail->value) {
4527  /* ignore value if test isn't true */
4528  if (qll->sum1 >= qll->tau) {
4529  RASTER_DEBUGF(4, "skipping %f", value);
4530  }
4531  else {
4532 
4533  /* delete last element */
4534  RASTER_DEBUGF(4, "deleting %f from list", qll->head->value);
4535  qle = qll->head->next;
4536  RASTER_DEBUGF(5, "to-be tail is %f with count %d", qle->value, qle->count);
4537  qle->count += qll->head->count;
4538  quantile_llist_index_delete(qll, qll->head);
4540  qll->head = qle;
4541  qll->count--;
4542  quantile_llist_index_update(qll, NULL, 0);
4543  RASTER_DEBUGF(5, "tail is %f with count %d", qll->head->value, qll->head->count);
4544 
4545  /* insert value */
4546  RASTER_DEBUGF(4, "adding %f to list", value);
4547  /* OPTIMIZATION: check to see if value exceeds last value */
4548  if (NULL != qll->tail && (value > qll->tail->value || FLT_EQ(value, qll->tail->value))) {
4549  idx = qll->count;
4550  qle = quantile_llist_insert(qll->tail, value, &idx);
4551  }
4552  /* OPTIMIZATION: use index if possible */
4553  else {
4554  qls = quantile_llist_index_search(qll, value, &idx);
4555  qle = quantile_llist_insert(qls, value, &idx);
4556  }
4557  if (NULL == qle) return NULL;
4558  RASTER_DEBUGF(5, "value added at index: %d => %f", idx, value);
4559  qll->count++;
4560  qll->sum1++;
4561 
4562  /* first element */
4563  if (NULL == qle->prev)
4564  qll->head = qle;
4565  /* last element */
4566  if (NULL == qle->next)
4567  qll->tail = qle;
4568 
4569  qll->sum2 = qll->sum1 - qll->tail->count;
4570 
4571  quantile_llist_index_update(qll, qle, idx);
4572 
4573  RASTER_DEBUGF(5, "qle, head, tail = %p, %p, %p", qle, qll->head, qll->tail);
4574 
4575  }
4576  }
4577  else {
4578  qle = qll->head;
4579  while (NULL != qle) {
4580  if (qle->value > value) {
4581  qle->count++;
4582  qll->sum1++;
4583  qll->sum2 = qll->sum1 - qll->tail->count;
4584  RASTER_DEBUGF(4, "incremented count of %f by 1 to %d", qle->value, qle->count);
4585  break;
4586  }
4587 
4588  qle = qle->next;
4589  }
4590  }
4591  }
4592 
4593  RASTER_DEBUGF(5, "sum2, tau = %d, %d", qll->sum2, qll->tau);
4594  if (qll->sum2 >= qll->tau) {
4595  /* AL-GEQ */
4596  if (qll->algeq) {
4597  RASTER_DEBUGF(4, "deleting first element %f from list", qll->head->value);
4598 
4599  if (NULL != qll->head->next) {
4600  qle = qll->head->next;
4601  qll->sum1 -= qll->head->count;
4602  qll->sum2 = qll->sum1 - qle->count;
4603  quantile_llist_index_delete(qll, qll->head);
4605  qll->head = qle;
4606  qll->count--;
4607 
4608  quantile_llist_index_update(qll, NULL, 0);
4609  }
4610  else {
4612  qll->head = NULL;
4613  qll->tail = NULL;
4614  qll->sum1 = 0;
4615  qll->sum2 = 0;
4616  qll->count = 0;
4617 
4619  }
4620  }
4621  /* AL-GT */
4622  else {
4623  RASTER_DEBUGF(4, "deleting first element %f from list", qll->tail->value);
4624 
4625  if (NULL != qll->tail->prev) {
4626  qle = qll->tail->prev;
4627  qll->sum1 -= qll->tail->count;
4628  qll->sum2 = qll->sum1 - qle->count;
4629  quantile_llist_index_delete(qll, qll->tail);
4631  qll->tail = qle;
4632  qll->count--;
4633  }
4634  else {
4636  qll->head = NULL;
4637  qll->tail = NULL;
4638  qll->sum1 = 0;
4639  qll->sum2 = 0;
4640  qll->count = 0;
4641 
4643  }
4644  }
4645  }
4646 
4647  RASTER_DEBUGF(5, "qll after: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
4648  qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
4649  RASTER_DEBUGF(5, "qll after: (head, tail) = (%p, %p)\n\n", qll->head, qll->tail);
4650  }
4651 
4652  }
4653  else {
4654  RASTER_DEBUGF(5, "skipping value at (x, y) = (%d, %d)", x, y);
4655  }
4656 
4657  z++;
4658  }
4659  }
4660 
4661  /* process quantiles */
4662  *rtn_count = *qlls_count / 2;
4663  rtn = rtalloc(sizeof(struct rt_quantile_t) * *rtn_count);
4664  if (NULL == rtn) return NULL;
4665 
4666  RASTER_DEBUGF(3, "returning %d quantiles", *rtn_count);
4667  for (i = 0, k = 0; i < *qlls_count; i++) {
4668  qll = &((*qlls)[i]);
4669 
4670  exists = 0;
4671  for (x = 0; x < k; x++) {
4672  if (FLT_EQ(qll->quantile, rtn[x].quantile)) {
4673  exists = 1;
4674  break;
4675  }
4676  }
4677  if (exists) continue;
4678 
4679  RASTER_DEBUGF(5, "qll: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
4680  qll->algeq, qll->quantile, qll->count, qll->tau, qll->sum1, qll->sum2);
4681  RASTER_DEBUGF(5, "qll: (head, tail) = (%p, %p)", qll->head, qll->tail);
4682 
4683  rtn[k].quantile = qll->quantile;
4684  rtn[k].has_value = 0;
4685 
4686  /* check that qll->head and qll->tail have value */
4687  if (qll->head == NULL || qll->tail == NULL)
4688  continue;
4689 
4690  /* AL-GEQ */
4691  if (qll->algeq)
4692  qle = qll->head;
4693  /* AM-GT */
4694  else
4695  qle = qll->tail;
4696 
4697  exists = 0;
4698  for (j = i + 1; j < *qlls_count; j++) {
4699  if (FLT_EQ((*qlls)[j].quantile, qll->quantile)) {
4700 
4701  RASTER_DEBUGF(5, "qlls[%d]: (algeq, quantile, count, tau, sum1, sum2) = (%d, %f, %d, %d, %d, %d)",
4702  j, (*qlls)[j].algeq, (*qlls)[j].quantile, (*qlls)[j].count, (*qlls)[j].tau, (*qlls)[j].sum1, (*qlls)[j].sum2);
4703  RASTER_DEBUGF(5, "qlls[%d]: (head, tail) = (%p, %p)", j, (*qlls)[j].head, (*qlls)[j].tail);
4704 
4705  exists = 1;
4706  break;
4707  }
4708  }
4709 
4710  /* weighted average for quantile */
4711  if (exists) {
4712  if ((*qlls)[j].algeq) {
4713  rtn[k].value = ((qle->value * qle->count) + ((*qlls)[j].head->value * (*qlls)[j].head->count)) / (qle->count + (*qlls)[j].head->count);
4714  RASTER_DEBUGF(5, "qlls[%d].head: (value, count) = (%f, %d)", j, (*qlls)[j].head->value, (*qlls)[j].head->count);
4715  }
4716  else {
4717  rtn[k].value = ((qle->value * qle->count) + ((*qlls)[j].tail->value * (*qlls)[j].tail->count)) / (qle->count + (*qlls)[j].tail->count);
4718  RASTER_DEBUGF(5, "qlls[%d].tail: (value, count) = (%f, %d)", j, (*qlls)[j].tail->value, (*qlls)[j].tail->count);
4719  }
4720  }
4721  /* straight value for quantile */
4722  else {
4723  rtn[k].value = qle->value;
4724  }
4725  rtn[k].has_value = 1;
4726  RASTER_DEBUGF(3, "(quantile, value) = (%f, %f)\n\n", rtn[k].quantile, rtn[k].value);
4727 
4728  k++;
4729  }
4730 
4731  RASTER_DEBUG(3, "done");
4732  return rtn;
4733 }
4734 
4751  rt_band band, int exclude_nodata_value,
4752  double *search_values, uint32_t search_values_count, double roundto,
4753  uint32_t *rtn_total, uint32_t *rtn_count
4754 ) {
4755  rt_valuecount vcnts = NULL;
4756  rt_pixtype pixtype = PT_END;
4757  uint8_t *data = NULL;
4758  double nodata = 0;
4759 
4760  int scale = 0;
4761  int doround = 0;
4762  double tmpd = 0;
4763  int i = 0;
4764 
4765  uint32_t x = 0;
4766  uint32_t y = 0;
4767  int rtn;
4768  double pxlval;
4769  int isnodata = 0;
4770  double rpxlval;
4771  uint32_t total = 0;
4772  int vcnts_count = 0;
4773  int new_valuecount = 0;
4774 
4775 #if POSTGIS_DEBUG_LEVEL > 0
4776  clock_t start, stop;
4777  double elapsed = 0;
4778 #endif
4779 
4780  RASTER_DEBUG(3, "starting");
4781 #if POSTGIS_DEBUG_LEVEL > 0
4782  start = clock();
4783 #endif
4784 
4785  assert(NULL != band);
4786  assert(NULL != rtn_count);
4787 
4788  data = rt_band_get_data(band);
4789  if (data == NULL) {
4790  rterror("rt_band_get_summary_stats: Cannot get band data");
4791  return NULL;
4792  }
4793 
4794  pixtype = band->pixtype;
4795 
4796  if (rt_band_get_hasnodata_flag(band)) {
4797  rt_band_get_nodata(band, &nodata);
4798  RASTER_DEBUGF(3, "hasnodata, nodataval = 1, %f", nodata);
4799  }
4800  else {
4801  exclude_nodata_value = 0;
4802  RASTER_DEBUG(3, "hasnodata, nodataval = 0, 0");
4803  }
4804 
4805  RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
4806 
4807  /* process roundto */
4808  if (roundto < 0 || FLT_EQ(roundto, 0.0)) {
4809  roundto = 0;
4810  scale = 0;
4811  }
4812  /* tenths, hundredths, thousandths, etc */
4813  else if (roundto < 1) {
4814  switch (pixtype) {
4815  /* integer band types don't have digits after the decimal place */
4816  case PT_1BB:
4817  case PT_2BUI:
4818  case PT_4BUI:
4819  case PT_8BSI:
4820  case PT_8BUI:
4821  case PT_16BSI:
4822  case PT_16BUI:
4823  case PT_32BSI:
4824  case PT_32BUI:
4825  roundto = 0;
4826  break;
4827  /* floating points, check the rounding */
4828  case PT_32BF:
4829  case PT_64BF:
4830  for (scale = 0; scale <= 20; scale++) {
4831  tmpd = roundto * pow(10, scale);
4832  if (FLT_EQ((tmpd - ((int) tmpd)), 0.0)) break;
4833  }
4834  break;
4835  case PT_END:
4836  break;
4837  }
4838  }
4839  /* ones, tens, hundreds, etc */
4840  else {
4841  for (scale = 0; scale >= -20; scale--) {
4842  tmpd = roundto * pow(10, scale);
4843  if (tmpd < 1 || FLT_EQ(tmpd, 1.0)) {
4844  if (scale == 0) doround = 1;
4845  break;
4846  }
4847  }
4848  }
4849 
4850  if (scale != 0 || doround)
4851  doround = 1;
4852  else
4853  doround = 0;
4854  RASTER_DEBUGF(3, "scale = %d", scale);
4855  RASTER_DEBUGF(3, "doround = %d", doround);
4856 
4857  /* process search_values */
4858  if (search_values_count > 0 && NULL != search_values) {
4859  vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t) * search_values_count);
4860  if (NULL == vcnts) {
4861  rterror("rt_band_get_count_of_values: Could not allocate memory for value counts");
4862  *rtn_count = 0;
4863  return NULL;
4864  }
4865 
4866  for (i = 0; i < search_values_count; i++) {
4867  vcnts[i].count = 0;
4868  vcnts[i].percent = 0;
4869  if (!doround)
4870  vcnts[i].value = search_values[i];
4871  else
4872  vcnts[i].value = ROUND(search_values[i], scale);
4873  }
4874  vcnts_count = i;
4875  }
4876  else
4877  search_values_count = 0;
4878  RASTER_DEBUGF(3, "search_values_count = %d", search_values_count);
4879 
4880  /* entire band is nodata */
4881  if (rt_band_get_isnodata_flag(band) != FALSE) {
4882  if (exclude_nodata_value) {
4883  rtwarn("All pixels of band have the NODATA value");
4884  return NULL;
4885  }
4886  else {
4887  if (search_values_count > 0) {
4888  /* check for nodata match */
4889  for (i = 0; i < search_values_count; i++) {
4890  if (!doround)
4891  tmpd = nodata;
4892  else
4893  tmpd = ROUND(nodata, scale);
4894 
4895  if (FLT_NEQ(tmpd, vcnts[i].value))
4896  continue;
4897 
4898  vcnts[i].count = band->width * band->height;
4899  if (NULL != rtn_total) *rtn_total = vcnts[i].count;
4900  vcnts->percent = 1.0;
4901  }
4902 
4903  *rtn_count = vcnts_count;
4904  }
4905  /* no defined search values */
4906  else {
4907  vcnts = (rt_valuecount) rtalloc(sizeof(struct rt_valuecount_t));
4908  if (NULL == vcnts) {
4909  rterror("rt_band_get_count_of_values: Could not allocate memory for value counts");
4910  *rtn_count = 0;
4911  return NULL;
4912  }
4913 
4914  vcnts->value = nodata;
4915  vcnts->count = band->width * band->height;
4916  if (NULL != rtn_total) *rtn_total = vcnts[i].count;
4917  vcnts->percent = 1.0;
4918 
4919  *rtn_count = 1;
4920  }
4921 
4922  return vcnts;
4923  }
4924  }
4925 
4926  for (x = 0; x < band->width; x++) {
4927  for (y = 0; y < band->height; y++) {
4928  rtn = rt_band_get_pixel(band, x, y, &pxlval, &isnodata);
4929 
4930  /* error getting value, continue */
4931  if (rtn != ES_NONE)
4932  continue;
4933 
4934  if (!exclude_nodata_value || (exclude_nodata_value && !isnodata)) {
4935  total++;
4936  if (doround) {
4937  rpxlval = ROUND(pxlval, scale);
4938  }
4939  else
4940  rpxlval = pxlval;
4941  RASTER_DEBUGF(5, "(pxlval, rpxlval) => (%0.6f, %0.6f)", pxlval, rpxlval);
4942 
4943  new_valuecount = 1;
4944  /* search for match in existing valuecounts */
4945  for (i = 0; i < vcnts_count; i++) {
4946  /* match found */
4947  if (FLT_EQ(vcnts[i].value, rpxlval)) {
4948  vcnts[i].count++;
4949  new_valuecount = 0;
4950  RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
4951  break;
4952  }
4953  }
4954 
4955  /*
4956  don't add new valuecount either because
4957  - no need for new one
4958  - user-defined search values
4959  */
4960  if (!new_valuecount || search_values_count > 0) continue;
4961 
4962  /* add new valuecount */
4963  vcnts = rtrealloc(vcnts, sizeof(struct rt_valuecount_t) * (vcnts_count + 1));
4964  if (NULL == vcnts) {
4965  rterror("rt_band_get_count_of_values: Could not allocate memory for value counts");
4966  *rtn_count = 0;
4967  return NULL;
4968  }
4969 
4970  vcnts[vcnts_count].value = rpxlval;
4971  vcnts[vcnts_count].count = 1;
4972  vcnts[vcnts_count].percent = 0;
4973  RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[vcnts_count].value, vcnts[vcnts_count].count);
4974  vcnts_count++;
4975  }
4976  }
4977  }
4978 
4979 #if POSTGIS_DEBUG_LEVEL > 0
4980  stop = clock();
4981  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
4982  RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
4983 #endif
4984 
4985  for (i = 0; i < vcnts_count; i++) {
4986  vcnts[i].percent = (double) vcnts[i].count / total;
4987  RASTER_DEBUGF(5, "(value, count) => (%0.6f, %d)", vcnts[i].value, vcnts[i].count);
4988  }
4989 
4990  RASTER_DEBUG(3, "done");
4991  if (NULL != rtn_total) *rtn_total = total;
4992  *rtn_count = vcnts_count;
4993  return vcnts;
4994 }
4995 
5008 rt_band
5010  rt_band srcband, rt_pixtype pixtype,
5011  uint32_t hasnodata, double nodataval,
5012  rt_reclassexpr *exprset, int exprcount
5013 ) {
5014  rt_band band = NULL;
5015  uint32_t width = 0;
5016  uint32_t height = 0;
5017  int numval = 0;
5018  int memsize = 0;
5019  void *mem = NULL;
5020  uint32_t src_hasnodata = 0;
5021  double src_nodataval = 0.0;
5022  int isnodata = 0;
5023 
5024  int rtn;
5025  uint32_t x;
5026  uint32_t y;
5027  int i;
5028  double or = 0;
5029  double ov = 0;
5030  double nr = 0;
5031  double nv = 0;
5032  int do_nv = 0;
5033  rt_reclassexpr expr = NULL;
5034 
5035  assert(NULL != srcband);
5036  assert(NULL != exprset && exprcount > 0);
5037  RASTER_DEBUGF(4, "exprcount = %d", exprcount);
5038  RASTER_DEBUGF(4, "exprset @ %p", exprset);
5039 
5040  /* source nodata */
5041  src_hasnodata = rt_band_get_hasnodata_flag(srcband);
5042  if (src_hasnodata)
5043  rt_band_get_nodata(srcband, &src_nodataval);
5044 
5045  /* size of memory block to allocate */
5046  width = rt_band_get_width(srcband);
5047  height = rt_band_get_height(srcband);
5048  numval = width * height;
5049  memsize = rt_pixtype_size(pixtype) * numval;
5050  mem = (int *) rtalloc(memsize);
5051  if (!mem) {
5052  rterror("rt_band_reclass: Could not allocate memory for band");
5053  return 0;
5054  }
5055 
5056  /* initialize to zero */
5057  if (!hasnodata) {
5058  memset(mem, 0, memsize);
5059  }
5060  /* initialize to nodataval */
5061  else {
5062  int32_t checkvalint = 0;
5063  uint32_t checkvaluint = 0;
5064  double checkvaldouble = 0;
5065  float checkvalfloat = 0;
5066 
5067  switch (pixtype) {
5068  case PT_1BB:
5069  {
5070  uint8_t *ptr = mem;
5071  uint8_t clamped_initval = rt_util_clamp_to_1BB(nodataval);
5072  for (i = 0; i < numval; i++)
5073  ptr[i] = clamped_initval;
5074  checkvalint = ptr[0];
5075  break;
5076  }
5077  case PT_2BUI:
5078  {
5079  uint8_t *ptr = mem;
5080  uint8_t clamped_initval = rt_util_clamp_to_2BUI(nodataval);
5081  for (i = 0; i < numval; i++)
5082  ptr[i] = clamped_initval;
5083  checkvalint = ptr[0];
5084  break;
5085  }
5086  case PT_4BUI:
5087  {
5088  uint8_t *ptr = mem;
5089  uint8_t clamped_initval = rt_util_clamp_to_4BUI(nodataval);
5090  for (i = 0; i < numval; i++)
5091  ptr[i] = clamped_initval;
5092  checkvalint = ptr[0];
5093  break;
5094  }
5095  case PT_8BSI:
5096  {
5097  int8_t *ptr = mem;
5098  int8_t clamped_initval = rt_util_clamp_to_8BSI(nodataval);
5099  for (i = 0; i < numval; i++)
5100  ptr[i] = clamped_initval;
5101  checkvalint = ptr[0];
5102  break;
5103  }
5104  case PT_8BUI:
5105  {
5106  uint8_t *ptr = mem;
5107  uint8_t clamped_initval = rt_util_clamp_to_8BUI(nodataval);
5108  for (i = 0; i < numval; i++)
5109  ptr[i] = clamped_initval;
5110  checkvalint = ptr[0];
5111  break;
5112  }
5113  case PT_16BSI:
5114  {
5115  int16_t *ptr = mem;
5116  int16_t clamped_initval = rt_util_clamp_to_16BSI(nodataval);
5117  for (i = 0; i < numval; i++)
5118  ptr[i] = clamped_initval;
5119  checkvalint = ptr[0];
5120  break;
5121  }
5122  case PT_16BUI:
5123  {
5124  uint16_t *ptr = mem;
5125  uint16_t clamped_initval = rt_util_clamp_to_16BUI(nodataval);
5126  for (i = 0; i < numval; i++)
5127  ptr[i] = clamped_initval;
5128  checkvalint = ptr[0];
5129  break;
5130  }
5131  case PT_32BSI:
5132  {
5133  int32_t *ptr = mem;
5134  int32_t clamped_initval = rt_util_clamp_to_32BSI(nodataval);
5135  for (i = 0; i < numval; i++)
5136  ptr[i] = clamped_initval;
5137  checkvalint = ptr[0];
5138  break;
5139  }
5140  case PT_32BUI:
5141  {
5142  uint32_t *ptr = mem;
5143  uint32_t clamped_initval = rt_util_clamp_to_32BUI(nodataval);
5144  for (i = 0; i < numval; i++)
5145  ptr[i] = clamped_initval;
5146  checkvaluint = ptr[0];
5147  break;
5148  }
5149  case PT_32BF:
5150  {
5151  float *ptr = mem;
5152  float clamped_initval = rt_util_clamp_to_32F(nodataval);
5153  for (i = 0; i < numval; i++)
5154  ptr[i] = clamped_initval;
5155  checkvalfloat = ptr[0];
5156  break;
5157  }
5158  case PT_64BF:
5159  {
5160  double *ptr = mem;
5161  for (i = 0; i < numval; i++)
5162  ptr[i] = nodataval;
5163  checkvaldouble = ptr[0];
5164  break;
5165  }
5166  default:
5167  {
5168  rterror("rt_band_reclass: Unknown pixeltype %d", pixtype);
5169  rtdealloc(mem);
5170  return 0;
5171  }
5172  }
5173 
5174  /* Overflow checking */
5176  nodataval,
5177  checkvalint, checkvaluint,
5178  checkvalfloat, checkvaldouble,
5179  pixtype
5180  );
5181  }
5182  RASTER_DEBUGF(3, "rt_band_reclass: width = %d height = %d", width, height);
5183 
5184  band = rt_band_new_inline(width, height, pixtype, hasnodata, nodataval, mem);
5185  if (!band) {
5186  rterror("rt_band_reclass: Could not create new band");
5187  rtdealloc(mem);
5188  return 0;
5189  }
5190  rt_band_set_ownsdata_flag(band, 1); /* we DO own this data!!! */
5191  RASTER_DEBUGF(3, "rt_band_reclass: new band @ %p", band);
5192 
5193  for (x = 0; x < width; x++) {
5194  for (y = 0; y < height; y++) {
5195  rtn = rt_band_get_pixel(srcband, x, y, &ov, &isnodata);
5196 
5197  /* error getting value, skip */
5198  if (rtn != ES_NONE) {
5199  RASTER_DEBUGF(3, "Cannot get value at %d, %d", x, y);
5200  continue;
5201  }
5202 
5203  do {
5204  do_nv = 0;
5205 
5206  /* no data*/
5207  if (hasnodata && isnodata) {
5208  do_nv = 1;
5209  break;
5210  }
5211 
5212  for (i = 0; i < exprcount; i++) {
5213  expr = exprset[i];
5214 
5215  /* ov matches min and max*/
5216  if (
5217  FLT_EQ(expr->src.min, ov) &&
5218  FLT_EQ(expr->src.max, ov)
5219  ) {
5220  do_nv = 1;
5221  break;
5222  }
5223 
5224  /* process min */
5225  if ((
5226  expr->src.exc_min && (
5227  expr->src.min > ov ||
5228  FLT_EQ(expr->src.min, ov)
5229  )) || (
5230  expr->src.inc_min && (
5231  expr->src.min < ov ||
5232  FLT_EQ(expr->src.min, ov)
5233  )) || (
5234  expr->src.min < ov
5235  )) {
5236  /* process max */
5237  if ((
5238  expr->src.exc_max && (
5239  ov > expr->src.max ||
5240  FLT_EQ(expr->src.max, ov)
5241  )) || (
5242  expr->src.inc_max && (
5243  ov < expr->src.max ||
5244  FLT_EQ(expr->src.max, ov)
5245  )) || (
5246  ov < expr->src.max
5247  )) {
5248  do_nv = 1;
5249  break;
5250  }
5251  }
5252  }
5253  }
5254  while (0);
5255 
5256  /* no expression matched, do not continue */
5257  if (!do_nv) continue;
5258  RASTER_DEBUGF(3, "Using exprset[%d] unless NODATA", i);
5259 
5260  /* converting a value from one range to another range
5261  OldRange = (OldMax - OldMin)
5262  NewRange = (NewMax - NewMin)
5263  NewValue = (((OldValue - OldMin) * NewRange) / OldRange) + NewMin
5264  */
5265 
5266  /* NODATA */
5267  if (hasnodata && isnodata) {
5268  nv = nodataval;
5269  }
5270  /*
5271  "src" min and max is the same, prevent division by zero
5272  set nv to "dst" min, which should be the same as "dst" max
5273  */
5274  else if (FLT_EQ(expr->src.max, expr->src.min)) {
5275  nv = expr->dst.min;
5276  }
5277  else {
5278  or = expr->src.max - expr->src.min;
5279  nr = expr->dst.max - expr->dst.min;
5280  nv = (((ov - expr->src.min) * nr) / or) + expr->dst.min;
5281 
5282  /* if dst range is from high to low */
5283  if (expr->dst.min > expr->dst.max) {
5284  if (nv > expr->dst.min)
5285  nv = expr->dst.min;
5286  else if (nv < expr->dst.max)
5287  nv = expr->dst.max;
5288  }
5289  /* if dst range is from low to high */
5290  else {
5291  if (nv < expr->dst.min)
5292  nv = expr->dst.min;
5293  else if (nv > expr->dst.max)
5294  nv = expr->dst.max;
5295  }
5296  }
5297 
5298  /* round the value for integers */
5299  switch (pixtype) {
5300  case PT_1BB:
5301  case PT_2BUI:
5302  case PT_4BUI:
5303  case PT_8BSI:
5304  case PT_8BUI:
5305  case PT_16BSI:
5306  case PT_16BUI:
5307  case PT_32BSI:
5308  case PT_32BUI:
5309  nv = round(nv);
5310  break;
5311  default:
5312  break;
5313  }
5314 
5315  RASTER_DEBUGF(4, "(%d, %d) ov: %f or: %f - %f nr: %f - %f nv: %f"
5316  , x
5317  , y
5318  , ov
5319  , (NULL != expr) ? expr->src.min : 0
5320  , (NULL != expr) ? expr->src.max : 0
5321  , (NULL != expr) ? expr->dst.min : 0
5322  , (NULL != expr) ? expr->dst.max : 0
5323  , nv
5324  );
5325  if (rt_band_set_pixel(band, x, y, nv, NULL) != ES_NONE) {
5326  rterror("rt_band_reclass: Could not assign value to new band");
5327  rt_band_destroy(band);
5328  rtdealloc(mem);
5329  return 0;
5330  }
5331 
5332  expr = NULL;
5333  }
5334  }
5335 
5336  return band;
5337 }
5338 
5339 /*- rt_raster --------------------------------------------------------*/
5340 
5352 rt_raster
5353 rt_raster_new(uint32_t width, uint32_t height) {
5354  rt_raster ret = NULL;
5355 
5356  ret = (rt_raster) rtalloc(sizeof (struct rt_raster_t));
5357  if (!ret) {
5358  rterror("rt_raster_new: Out of virtual memory creating an rt_raster");
5359  return NULL;
5360  }
5361 
5362  RASTER_DEBUGF(3, "Created rt_raster @ %p", ret);
5363 
5364  if (width > 65535 || height > 65535) {
5365  rterror("rt_raster_new: Dimensions requested exceed the maximum (65535 x 65535) permitted for a raster");
5366  rt_raster_destroy(ret);
5367  return NULL;
5368  }
5369 
5370  ret->width = width;
5371  ret->height = height;
5372  ret->scaleX = 1;
5373  ret->scaleY = -1;
5374  ret->ipX = 0.0;
5375  ret->ipY = 0.0;
5376  ret->skewX = 0.0;
5377  ret->skewY = 0.0;
5378  ret->srid = SRID_UNKNOWN;
5379 
5380  ret->numBands = 0;
5381  ret->bands = NULL;
5382 
5383  return ret;
5384 }
5385 
5386 void
5388  if (raster == NULL)
5389  return;
5390 
5391  RASTER_DEBUGF(3, "Destroying rt_raster @ %p", raster);
5392 
5393  if (raster->bands)
5394  rtdealloc(raster->bands);
5395 
5396  rtdealloc(raster);
5397 }
5398 
5399 static void
5401  int numband = 0;
5402  int i = 0;
5403  rt_band band = NULL;
5404 
5405  if (raster == NULL)
5406  return;
5407 
5408  numband = rt_raster_get_num_bands(raster);
5409  if (numband < 1)
5410  return;
5411 
5412  for (i = 0; i < numband; i++) {
5413  band = rt_raster_get_band(raster, i);
5414  if (NULL == band)
5415  continue;
5416 
5417  if (!rt_band_is_offline(band))
5418  continue;
5419 
5420  rtwarn("Changes made to raster geotransform matrix may affect out-db band data. Returned band data may be incorrect");
5421  break;
5422  }
5423 }
5424 
5425 uint16_t
5427 
5428  assert(NULL != raster);
5429 
5430  return raster->width;
5431 }
5432 
5433 uint16_t
5435 
5436  assert(NULL != raster);
5437 
5438  return raster->height;
5439 }
5440 
5441 void
5443  rt_raster raster,
5444  double scaleX, double scaleY
5445 ) {
5446  assert(NULL != raster);
5447 
5448  raster->scaleX = scaleX;
5449  raster->scaleY = scaleY;
5450 
5452 }
5453 
5454 double
5456 
5457 
5458  assert(NULL != raster);
5459 
5460  return raster->scaleX;
5461 }
5462 
5463 double
5465 
5466 
5467  assert(NULL != raster);
5468 
5469  return raster->scaleY;
5470 }
5471 
5472 void
5474  rt_raster raster,
5475  double skewX, double skewY
5476 ) {
5477  assert(NULL != raster);
5478 
5479  raster->skewX = skewX;
5480  raster->skewY = skewY;
5481 
5483 }
5484 
5485 double
5487 
5488 
5489  assert(NULL != raster);
5490 
5491  return raster->skewX;
5492 }
5493 
5494 double
5496 
5497 
5498  assert(NULL != raster);
5499 
5500  return raster->skewY;
5501 }
5502 
5503 void
5505  rt_raster raster,
5506  double x, double y
5507 ) {
5508 
5509  assert(NULL != raster);
5510 
5511  raster->ipX = x;
5512  raster->ipY = y;
5513 
5515 }
5516 
5517 double
5519 
5520 
5521  assert(NULL != raster);
5522 
5523  return raster->ipX;
5524 }
5525 
5526 double
5528 
5529 
5530  assert(NULL != raster);
5531 
5532  return raster->ipY;
5533 }
5534 
5535 void
5537  double *i_mag, double *j_mag, double *theta_i, double *theta_ij)
5538 {
5539  double o11, o12, o21, o22 ; /* geotransform coefficients */
5540 
5541  if (rast == NULL) return ;
5542  if ( (i_mag==NULL) || (j_mag==NULL) || (theta_i==NULL) || (theta_ij==NULL))
5543  return ;
5544 
5545  /* retrieve coefficients from raster */
5546  o11 = rt_raster_get_x_scale(rast) ;
5547  o12 = rt_raster_get_x_skew(rast) ;
5548  o21 = rt_raster_get_y_skew(rast) ;
5549  o22 = rt_raster_get_y_scale(rast) ;
5550 
5551  rt_raster_calc_phys_params(o11, o12, o21, o22, i_mag, j_mag, theta_i, theta_ij);
5552 }
5553 
5554 void
5555 rt_raster_calc_phys_params(double xscale, double xskew, double yskew, double yscale,
5556  double *i_mag, double *j_mag, double *theta_i, double *theta_ij)
5557 
5558 {
5559  double theta_test ;
5560 
5561  if ( (i_mag==NULL) || (j_mag==NULL) || (theta_i==NULL) || (theta_ij==NULL))
5562  return ;
5563 
5564  /* pixel size in the i direction */
5565  *i_mag = sqrt(xscale*xscale + yskew*yskew) ;
5566 
5567  /* pixel size in the j direction */
5568  *j_mag = sqrt(xskew*xskew + yscale*yscale) ;
5569 
5570  /* Rotation
5571  * ========
5572  * Two steps:
5573  * 1] calculate the magnitude of the angle between the x axis and
5574  * the i basis vector.
5575  * 2] Calculate the sign of theta_i based on the angle between the y axis
5576  * and the i basis vector.
5577  */
5578  *theta_i = acos(xscale/(*i_mag)) ; /* magnitude */
5579  theta_test = acos(yskew/(*i_mag)) ; /* sign */
5580  if (theta_test < M_PI_2){
5581  *theta_i = -(*theta_i) ;
5582  }
5583 
5584 
5585  /* Angular separation of basis vectors
5586  * ===================================
5587  * Two steps:
5588  * 1] calculate the magnitude of the angle between the j basis vector and
5589  * the i basis vector.
5590  * 2] Calculate the sign of theta_ij based on the angle between the
5591  * perpendicular of the i basis vector and the j basis vector.
5592  */
5593  *theta_ij = acos(((xscale*xskew) + (yskew*yscale))/((*i_mag)*(*j_mag))) ;
5594  theta_test = acos( ((-yskew*xskew)+(xscale*yscale)) /
5595  ((*i_mag)*(*j_mag)));
5596  if (theta_test > M_PI_2) {
5597  *theta_ij = -(*theta_ij) ;
5598  }
5599 }
5600 
5601 void
5602 rt_raster_set_phys_params(rt_raster rast,double i_mag, double j_mag, double theta_i, double theta_ij)
5603 {
5604  double o11, o12, o21, o22 ; /* calculated geotransform coefficients */
5605  int success ;
5606 
5607  if (rast == NULL) return ;
5608 
5609  success = rt_raster_calc_gt_coeff(i_mag, j_mag, theta_i, theta_ij,
5610  &o11, &o12, &o21, &o22) ;
5611 
5612  if (success) {
5613  rt_raster_set_scale(rast, o11, o22) ;
5614  rt_raster_set_skews(rast, o12, o21) ;
5615  }
5616 }
5617 
5618 int
5619 rt_raster_calc_gt_coeff(double i_mag, double j_mag, double theta_i, double theta_ij,
5620  double *xscale, doubl