PostGIS  2.1.10dev-r@@SVN_REVISION@@
rt_histogram rt_band_get_histogram ( rt_bandstats  stats,
int  bin_count,
double *  bin_width,
int  bin_width_count,
int  right,
double  min,
double  max,
uint32_t *  rtn_count 
)

Count the distribution of data.

Parameters
stats: a populated stats struct for processing
bin_count: the number of bins to group the data by
bin_width: the width of each bin as an array
bin_width_count: number of values in bin_width
right: evaluate bins by (a,b] rather than default [a,b)
min: user-defined minimum value of the histogram a value less than the minimum value is not counted in any bins if min = max, min and max are not used
max: user-defined maximum value of the histogram a value greater than the max value is not counted in any bins if min = max, min and max are not used
rtn_count: set to the number of bins being returned
Returns
the histogram of the data or NULL

Definition at line 3567 of file rt_api.c.

References genraster::count, rt_bandstats_t::count, rt_histogram_t::count, FLT_EQ, rt_histogram_t::inc_max, rt_histogram_t::inc_min, rt_bandstats_t::max, rt_histogram_t::max, rt_bandstats_t::min, rt_histogram_t::min, rt_histogram_t::percent, RASTER_DEBUG, RASTER_DEBUGF, rtalloc(), rtdealloc(), rterror(), genraster::value, and rt_bandstats_t::values.

Referenced by RASTER_histogram(), RASTER_histogramCoverage(), and test_band_stats().

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 }
void rtdealloc(void *mem)
Definition: rt_api.c:882
uint32_t count
Definition: rt_api.h:2277
double max
Definition: rt_api.h:2295
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
double min
Definition: rt_api.h:2279
double percent
Definition: rt_api.h:2292
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
int count
Definition: genraster.py:57
double * values
Definition: rt_api.h:2285
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
double min
Definition: rt_api.h:2294
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
double max
Definition: rt_api.h:2280
uint32_t count
Definition: rt_api.h:2291

Here is the call graph for this function:

Here is the caller graph for this function: