PostGIS  3.0.0dev-r@@SVN_REVISION@@

◆ rt_band_get_summary_stats()

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 110 of file rt_statistics.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_summaryStats_transfn(), RASTER_summaryStatsCoverage(), and test_band_stats().

114  {
115  uint32_t x = 0;
116  uint32_t y = 0;
117  uint32_t z = 0;
118  uint32_t offset = 0;
119  uint32_t diff = 0;
120  int rtn;
121  int hasnodata = FALSE;
122  double nodata = 0;
123  double *values = NULL;
124  double value;
125  int isnodata = 0;
126  rt_bandstats stats = NULL;
127 
128  uint32_t do_sample = 0;
129  uint32_t sample_size = 0;
130  uint32_t sample_per = 0;
131  uint32_t sample_int = 0;
132  uint32_t i = 0;
133  uint32_t j = 0;
134  double sum = 0;
135  uint32_t k = 0;
136  double M = 0;
137  double Q = 0;
138 
139 #if POSTGIS_DEBUG_LEVEL > 0
140  clock_t start, stop;
141  double elapsed = 0;
142 #endif
143 
144  RASTER_DEBUG(3, "starting");
145 #if POSTGIS_DEBUG_LEVEL > 0
146  start = clock();
147 #endif
148 
149  assert(NULL != band);
150 
151  /* band is empty (width < 1 || height < 1) */
152  if (band->width < 1 || band->height < 1) {
153  stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
154  if (NULL == stats) {
155  rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
156  return NULL;
157  }
158 
159  rtwarn("Band is empty as width and/or height is 0");
160 
161  stats->sample = 1;
162  stats->sorted = 0;
163  stats->values = NULL;
164  stats->count = 0;
165  stats->min = stats->max = 0;
166  stats->sum = 0;
167  stats->mean = 0;
168  stats->stddev = -1;
169 
170  return stats;
171  }
172 
173  hasnodata = rt_band_get_hasnodata_flag(band);
174  if (hasnodata != FALSE)
175  rt_band_get_nodata(band, &nodata);
176  else
177  exclude_nodata_value = 0;
178 
179  RASTER_DEBUGF(3, "nodata = %f", nodata);
180  RASTER_DEBUGF(3, "hasnodata = %d", hasnodata);
181  RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
182 
183  /* entire band is nodata */
184  if (rt_band_get_isnodata_flag(band) != FALSE) {
185  stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
186  if (NULL == stats) {
187  rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
188  return NULL;
189  }
190 
191  stats->sample = 1;
192  stats->sorted = 0;
193  stats->values = NULL;
194 
195  if (exclude_nodata_value) {
196  rtwarn("All pixels of band have the NODATA value");
197 
198  stats->count = 0;
199  stats->min = stats->max = 0;
200  stats->sum = 0;
201  stats->mean = 0;
202  stats->stddev = -1;
203  }
204  else {
205  stats->count = band->width * band->height;
206  stats->min = stats->max = nodata;
207  stats->sum = stats->count * nodata;
208  stats->mean = nodata;
209  stats->stddev = 0;
210  }
211 
212  return stats;
213  }
214 
215  /* clamp percentage */
216  if (
217  (sample < 0 || FLT_EQ(sample, 0.0)) ||
218  (sample > 1 || FLT_EQ(sample, 1.0))
219  ) {
220  do_sample = 0;
221  sample = 1;
222  }
223  else
224  do_sample = 1;
225  RASTER_DEBUGF(3, "do_sample = %d", do_sample);
226 
227  /* sample all pixels */
228  if (!do_sample) {
229  sample_size = band->width * band->height;
230  sample_per = band->height;
231  }
232  /*
233  randomly sample a percentage of available pixels
234  sampling method is known as
235  "systematic random sample without replacement"
236  */
237  else {
238  sample_size = round((band->width * band->height) * sample);
239  sample_per = round(sample_size / band->width);
240  if (sample_per < 1)
241  sample_per = 1;
242  sample_int = round(band->height / sample_per);
243  srand(time(NULL));
244  }
245 
246  RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set"
247  , sample_size, (band->width * band->height), sample_per);
248 
249  if (inc_vals) {
250  values = rtalloc(sizeof(double) * sample_size);
251  if (NULL == values) {
252  rtwarn("Could not allocate memory for values");
253  inc_vals = 0;
254  }
255  }
256 
257  /* initialize stats */
258  stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
259  if (NULL == stats) {
260  rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
261  return NULL;
262  }
263  stats->sample = sample;
264  stats->count = 0;
265  stats->sum = 0;
266  stats->mean = 0;
267  stats->stddev = -1;
268  stats->min = stats->max = 0;
269  stats->values = NULL;
270  stats->sorted = 0;
271 
272  for (x = 0, j = 0, k = 0; x < band->width; x++) {
273  y = -1;
274  diff = 0;
275 
276  for (i = 0, z = 0; i < sample_per; i++) {
277  if (!do_sample)
278  y = i;
279  else {
280  offset = (rand() % sample_int) + 1;
281  y += diff + offset;
282  diff = sample_int - offset;
283  }
284  RASTER_DEBUGF(5, "(x, y, z) = (%d, %d, %d)", x, y, z);
285  if (y >= band->height || z > sample_per) break;
286 
287  rtn = rt_band_get_pixel(band, x, y, &value, &isnodata);
288 
289  j++;
290  if (rtn == ES_NONE && (!exclude_nodata_value || (exclude_nodata_value && !isnodata))) {
291 
292  /* inc_vals set, collect pixel values */
293  if (inc_vals) values[k] = value;
294 
295  /* average */
296  k++;
297  sum += value;
298 
299  /*
300  one-pass standard deviation
301  http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
302  */
303  if (k == 1) {
304  Q = 0;
305  M = value;
306  }
307  else {
308  Q += (((k - 1) * pow(value - M, 2)) / k);
309  M += ((value - M ) / k);
310  }
311 
312  /* coverage one-pass standard deviation */
313  if (NULL != cK) {
314  (*cK)++;
315  if (*cK == 1) {
316  *cQ = 0;
317  *cM = value;
318  }
319  else {
320  *cQ += (((*cK - 1) * pow(value - *cM, 2)) / *cK);
321  *cM += ((value - *cM ) / *cK);
322  }
323  }
324 
325  /* min/max */
326  if (stats->count < 1) {
327  stats->count = 1;
328  stats->min = stats->max = value;
329  }
330  else {
331  if (value < stats->min)
332  stats->min = value;
333  if (value > stats->max)
334  stats->max = value;
335  }
336 
337  }
338 
339  z++;
340  }
341  }
342 
343  RASTER_DEBUG(3, "sampling complete");
344 
345  stats->count = k;
346  if (k > 0) {
347  if (inc_vals) {
348  /* free unused memory */
349  if (sample_size != k) {
350  values = rtrealloc(values, k * sizeof(double));
351  }
352 
353  stats->values = values;
354  }
355 
356  stats->sum = sum;
357  stats->mean = sum / k;
358 
359  /* standard deviation */
360  if (!do_sample)
361  stats->stddev = sqrt(Q / k);
362  /* sample deviation */
363  else {
364  if (k < 2)
365  stats->stddev = -1;
366  else
367  stats->stddev = sqrt(Q / (k - 1));
368  }
369  }
370  /* inc_vals thus values allocated but not used */
371  else if (inc_vals)
372  rtdealloc(values);
373 
374  /* if do_sample is one */
375  if (do_sample && k < 1)
376  rtwarn("All sampled pixels of band have the NODATA value");
377 
378 #if POSTGIS_DEBUG_LEVEL > 0
379  stop = clock();
380  elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
381  RASTER_DEBUGF(3, "(time, count, mean, stddev, min, max) = (%0.4f, %d, %f, %f, %f, %f)",
382  elapsed, stats->count, stats->mean, stats->stddev, stats->min, stats->max);
383 #endif
384 
385  RASTER_DEBUG(3, "done");
386  return stats;
387 }
uint32_t count
Definition: librtcore.h:2360
#define FLT_EQ(x, y)
Definition: librtcore.h:2234
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
void * rtrealloc(void *mem, size_t size)
Definition: rt_context.c:179
uint16_t height
Definition: librtcore.h:2317
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1730
unsigned int uint32_t
Definition: uthash.h:78
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1221
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:674
double * values
Definition: librtcore.h:2368
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
uint16_t width
Definition: librtcore.h:2316
struct rt_bandstats_t * rt_bandstats
Definition: librtcore.h:150
void rtdealloc(void *mem)
Definition: rt_context.c:186
#define FALSE
Definition: dbfopen.c:168
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:714
int value
Definition: genraster.py:61
Here is the call graph for this function:
Here is the caller graph for this function: