PostGIS  2.1.10dev-r@@SVN_REVISION@@
rt_bandstats rt_band_get_summary_stats ( rt_band  band,
int  exclude_nodata_value,
double  sample,
int  inc_vals,
uint64_t *  cK,
double *  cM,
double *  cQ 
)

Compute summary statistics for a band.

Parameters
band: the band to query for summary stats
exclude_nodata_value: if non-zero, ignore nodata values
sample: percentage of pixels to sample
inc_vals: flag to include values in return struct
cK: number of pixels counted thus far in coverage
cM: M component of 1-pass stddev for coverage
cQ: Q component of 1-pass stddev for coverage
Returns
the summary statistics for a band or NULL

Definition at line 3269 of file rt_api.c.

References rt_bandstats_t::count, ES_NONE, FALSE, FLT_EQ, rt_band_t::height, rt_bandstats_t::max, rt_bandstats_t::mean, rt_bandstats_t::min, RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_band_get_nodata(), rt_band_get_pixel(), rtalloc(), rtdealloc(), rterror(), rtrealloc(), rtwarn(), rt_bandstats_t::sample, rt_bandstats_t::sorted, rt_bandstats_t::stddev, rt_bandstats_t::sum, genraster::value, rt_bandstats_t::values, rt_band_t::width, pixval::x, and pixval::y.

Referenced by RASTER_colorMap(), RASTER_histogram(), RASTER_histogramCoverage(), RASTER_quantile(), RASTER_summaryStats(), RASTER_summaryStatsCoverage(), and test_band_stats().

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 }
double mean
Definition: rt_api.h:2282
void rtdealloc(void *mem)
Definition: rt_api.c:882
uint32_t count
Definition: rt_api.h:2277
struct rt_bandstats_t * rt_bandstats
Definition: rt_api.h:137
double sum
Definition: rt_api.h:2281
uint16_t height
Definition: rt_api.h:2242
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_api.c:2042
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
double min
Definition: rt_api.h:2279
void rtwarn(const char *fmt,...)
Definition: rt_api.c:920
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_api.c:3058
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_api.c:2549
double * values
Definition: rt_api.h:2285
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
uint16_t width
Definition: rt_api.h:2241
tuple x
Definition: pixval.py:53
double sample
Definition: rt_api.h:2276
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
#define FALSE
Definition: dbfopen.c:169
double max
Definition: rt_api.h:2280
double stddev
Definition: rt_api.h:2283
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_api.c:2002
tuple y
Definition: pixval.py:54
void * rtrealloc(void *mem, size_t size)
Definition: rt_api.c:875

Here is the call graph for this function:

Here is the caller graph for this function: