PostGIS  2.1.10dev-r@@SVN_REVISION@@
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: the number of quantiles being returned
Returns
the default set of or requested quantiles for a band or NULL
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 3823 of file rt_api.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().

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 }
static void quicksort(double *left, double *right)
Definition: rt_api.c:183
void rtdealloc(void *mem)
Definition: rt_api.c:882
uint32_t count
Definition: rt_api.h:2277
double quantile
Definition: rt_api.h:2303
double value
Definition: rt_api.h:2304
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
double * values
Definition: rt_api.h:2285
void * rtalloc(size_t size)
Raster core memory management functions.
Definition: rt_api.c:867
void rterror(const char *fmt,...)
Raster core error and info handlers.
Definition: rt_api.c:895

Here is the call graph for this function:

Here is the caller graph for this function: