PostGIS  2.3.8dev-r@@SVN_REVISION@@

◆ rt_band_get_quantiles()

rt_quantile rt_band_get_quantiles ( rt_bandstats  stats,
double *  quantiles,
int  quantiles_count,
uint32_t *  rtn_count 
)

Compute the default set of or requested quantiles for a set of data the quantile formula used is same as Excel and R default method.

Parameters
stats: a populated stats struct for processing
quantiles: the quantiles to be computed
quantiles_count: the number of quantiles to be computed
rtn_count: set to the number of quantiles being returned
Returns
the default set of or requested quantiles for a band or NULL

Definition at line 672 of file rt_statistics.c.

References rt_bandstats_t::count, rt_quantile_t::quantile, quicksort(), RASTER_DEBUG, RASTER_DEBUGF, rtalloc(), rtdealloc(), rterror(), rt_bandstats_t::sorted, rt_quantile_t::value, and rt_bandstats_t::values.

Referenced by RASTER_quantile(), and test_band_stats().

676  {
677  rt_quantile rtn;
678  int init_quantiles = 0;
679  int i = 0;
680  double h;
681  int hl;
682 
683 #if POSTGIS_DEBUG_LEVEL > 0
684  clock_t start, stop;
685  double elapsed = 0;
686 #endif
687 
688  RASTER_DEBUG(3, "starting");
689 #if POSTGIS_DEBUG_LEVEL > 0
690  start = clock();
691 #endif
692 
693  assert(NULL != stats);
694  assert(NULL != rtn_count);
695 
696  if (stats->count < 1 || NULL == stats->values) {
697  rterror("rt_band_get_quantiles: rt_bandstats object has no value");
698  return NULL;
699  }
700 
701  /* quantiles not provided */
702  if (NULL == quantiles) {
703  /* quantile count not specified, default to quartiles */
704  if (quantiles_count < 2)
705  quantiles_count = 5;
706 
707  quantiles = rtalloc(sizeof(double) * quantiles_count);
708  init_quantiles = 1;
709  if (NULL == quantiles) {
710  rterror("rt_band_get_quantiles: Could not allocate memory for quantile input");
711  return NULL;
712  }
713 
714  quantiles_count--;
715  for (i = 0; i <= quantiles_count; i++)
716  quantiles[i] = ((double) i) / quantiles_count;
717  quantiles_count++;
718  }
719 
720  /* check quantiles */
721  for (i = 0; i < quantiles_count; i++) {
722  if (quantiles[i] < 0. || quantiles[i] > 1.) {
723  rterror("rt_band_get_quantiles: Quantile value not between 0 and 1");
724  if (init_quantiles) rtdealloc(quantiles);
725  return NULL;
726  }
727  }
728  quicksort(quantiles, quantiles + quantiles_count - 1);
729 
730  /* initialize rt_quantile */
731  rtn = rtalloc(sizeof(struct rt_quantile_t) * quantiles_count);
732  if (NULL == rtn) {
733  rterror("rt_band_get_quantiles: Could not allocate memory for quantile output");
734  if (init_quantiles) rtdealloc(quantiles);
735  return NULL;
736  }
737 
738  /* sort values */
739  if (!stats->sorted) {
740  quicksort(stats->values, stats->values + stats->count - 1);
741  stats->sorted = 1;
742  }
743 
744  /*
745  make quantiles
746 
747  formula is that used in R (method 7) and Excel from
748  http://en.wikipedia.org/wiki/Quantile
749  */
750  for (i = 0; i < quantiles_count; i++) {
751  rtn[i].quantile = quantiles[i];
752 
753  h = ((stats->count - 1.) * quantiles[i]) + 1.;
754  hl = floor(h);
755 
756  /* h greater than hl, do full equation */
757  if (h > hl)
758  rtn[i].value = stats->values[hl - 1] + ((h - hl) * (stats->values[hl] - stats->values[hl - 1]));
759  /* shortcut as second part of equation is zero */
760  else
761  rtn[i].value = stats->values[hl - 1];
762  }
763 
764 #if POSTGIS_DEBUG_LEVEL > 0
765  stop = clock();
766  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
767  RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
768 #endif
769 
770  *rtn_count = quantiles_count;
771  if (init_quantiles) rtdealloc(quantiles);
772  RASTER_DEBUG(3, "done");
773  return rtn;
774 }
uint32_t count
Definition: librtcore.h:2323
double quantile
Definition: librtcore.h:2349
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
double value
Definition: librtcore.h:2350
static void quicksort(double *left, double *right)
Definition: rt_statistics.c:81
double * values
Definition: librtcore.h:2331
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:311
void rtdealloc(void *mem)
Definition: rt_context.c:186
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:307
Here is the call graph for this function:
Here is the caller graph for this function: