[GRASS-SVN] r30301 - in grass/trunk/raster: . r.quantile
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Feb 23 05:12:26 EST 2008
Author: glynn
Date: 2008-02-23 05:12:26 -0500 (Sat, 23 Feb 2008)
New Revision: 30301
Added:
grass/trunk/raster/r.quantile/
grass/trunk/raster/r.quantile/Makefile
grass/trunk/raster/r.quantile/main.c
Log:
New module: r.quantile
computes quantiles in a manner suitable for use with large amounts of data
Added: grass/trunk/raster/r.quantile/Makefile
===================================================================
--- grass/trunk/raster/r.quantile/Makefile (rev 0)
+++ grass/trunk/raster/r.quantile/Makefile 2008-02-23 10:12:26 UTC (rev 30301)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.quantile
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass/trunk/raster/r.quantile/main.c
===================================================================
--- grass/trunk/raster/r.quantile/main.c (rev 0)
+++ grass/trunk/raster/r.quantile/main.c 2008-02-23 10:12:26 UTC (rev 30301)
@@ -0,0 +1,290 @@
+/****************************************************************************
+ *
+ * MODULE: r.quantile
+ * AUTHOR(S): Glynn Clements <glynn gclements.plus.com> (original contributor),
+ * PURPOSE: Compute quantiles using two passes
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+struct bin
+{
+ unsigned long origin;
+ DCELL min, max;
+ int base, count;
+};
+
+static int rows, cols;
+
+static DCELL min, max;
+static int num_quant;
+static int num_slots;
+static unsigned int *slots;
+static DCELL slot_size;
+static unsigned long total;
+static int num_values;
+static unsigned char *slot_bins;
+static int num_bins;
+static struct bin *bins;
+static DCELL *values;
+
+static inline int get_slot(DCELL c)
+{
+ int i = (int) floor((c - min) / slot_size);
+ if (i < 0)
+ i = 0;
+ if (i > num_slots - 1)
+ i = num_slots - 1;
+ return i;
+}
+
+static inline double get_quantile(int n)
+{
+ return (double) total * n / num_quant;
+}
+
+static void get_slot_counts(int infile)
+{
+ DCELL *inbuf = G_allocate_d_raster_buf();
+ int row, col;
+
+ G_message(_("Computing histogram"));
+
+ total = 0;
+
+ for (row = 0; row < rows; row++)
+ {
+ G_get_d_raster_row(infile, inbuf, row);
+
+ for (col = 0; col < cols; col++)
+ {
+ int i;
+
+ if (G_is_d_null_value(&inbuf[col]))
+ continue;
+
+ i = get_slot(inbuf[col]);
+
+ slots[i]++;
+ total++;
+ }
+
+ G_percent(row, rows, 2);
+ }
+
+ G_percent(rows, rows, 2);
+ G_free(inbuf);
+}
+
+static void initialize_bins(void)
+{
+ int slot;
+ double next;
+ int bin = 0;
+ unsigned long accum = 0;
+ int quant = 1;
+
+ G_message(_("Computing bins"));
+
+ num_values = 0;
+ next = get_quantile(quant);
+
+ for (slot = 0; slot < num_slots; slot++)
+ {
+ unsigned int count = slots[slot];
+ unsigned long accum2 = accum + count;
+
+ if (accum2 > next)
+ {
+ struct bin *b = &bins[bin];
+
+ slot_bins[slot] = ++bin;
+
+ b->origin = accum;
+ b->base = num_values;
+ b->count = 0;
+ b->min = min + slot_size * slot;
+ b->max = min + slot_size * (slot + 1);
+
+ while (accum2 > next)
+ next = get_quantile(++quant);
+
+ num_values += count;
+ }
+
+ accum = accum2;
+ }
+
+ num_bins = bin;
+}
+
+static void fill_bins(int infile)
+{
+ DCELL *inbuf = G_allocate_d_raster_buf();
+ int row, col;
+
+ G_message(_("Binning data"));
+
+ for (row = 0; row < rows; row++)
+ {
+ G_get_d_raster_row(infile, inbuf, row);
+
+ for (col = 0; col < cols; col++)
+ {
+ int i, bin;
+ struct bin *b;
+
+ if (G_is_d_null_value(&inbuf[col]))
+ continue;
+
+ i = get_slot(inbuf[col]);
+ if (!slot_bins[i])
+ continue;
+
+ bin = slot_bins[i] - 1;
+ b = &bins[bin];
+
+ values[b->base + b->count++] = inbuf[col];
+ }
+
+ G_percent(row, rows, 2);
+ }
+
+ G_percent(rows, rows, 2);
+ G_free(inbuf);
+}
+
+static int compare_dcell(const void *aa, const void *bb)
+{
+ DCELL a = *(const DCELL *) aa;
+ DCELL b = *(const DCELL *) bb;
+
+ if (a < b)
+ return -1;
+ if (a > b)
+ return 1;
+ return 0;
+}
+
+static void sort_bins(void)
+{
+ int bin;
+
+ G_message(_("Sorting bins"));
+
+ for (bin = 0; bin < num_bins; bin++)
+ {
+ struct bin *b = &bins[bin];
+ qsort(&values[b->base], b->count, sizeof(DCELL), compare_dcell);
+ G_percent(bin, num_bins, 2);
+ }
+ G_percent(bin, num_bins, 2);
+
+ G_message(_("Computing quantiles"));
+}
+
+static void compute_quantiles(void)
+{
+ struct bin *b = &bins[0];
+ int quant;
+
+ for (quant = 1; quant < num_quant; quant++)
+ {
+ double next = get_quantile(quant);
+ double k, v;
+ int i0, i1;
+
+ while (b->origin + b->count < next)
+ b++;
+
+ k = next - b->origin;
+ i0 = (int) floor(k);
+ i1 = (int) ceil(k);
+
+ v = (i0 == i1)
+ ? values[b->base + i0]
+ : values[b->base + i0] * (i1 - k) + values[b->base + i1] * (k - i0);
+
+ printf("%2d: %f\n", quant, v);
+ }
+}
+
+int main( int argc, char *argv[])
+{
+ struct GModule *module;
+ struct {
+ struct Option *input, *quant, *slots;
+ } opt;
+ int infile;
+ struct FPRange range;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("raster, statistics");
+ module->description =
+ _("Compute quantiles using two passes.");
+
+ opt.input = G_define_standard_option(G_OPT_R_INPUT);
+
+ opt.quant = G_define_option();
+ opt.quant->key = "quantiles";
+ opt.quant->type = TYPE_INTEGER;
+ opt.quant->required = NO;
+ opt.quant->description = _("Number of quantiles");
+ opt.quant->answer = "4";
+
+ opt.slots = G_define_option();
+ opt.slots->key = "bins";
+ opt.slots->type = TYPE_INTEGER;
+ opt.slots->required = NO;
+ opt.slots->description = _("Number of bins to use");
+ opt.slots->answer = "1000000";
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ num_quant = atoi(opt.quant->answer);
+ num_slots = atoi(opt.slots->answer);
+
+ infile = G_open_cell_old(opt.input->answer, "");
+ if (infile < 0)
+ G_fatal_error(_("Unable to open raster map <%s>"), opt.input->answer);
+
+ G_read_fp_range(opt.input->answer, "", &range);
+ G_get_fp_range_min_max(&range, &min, &max);
+
+ slots = G_calloc(num_slots, sizeof(unsigned int));
+ slot_bins = G_calloc(num_slots, sizeof(unsigned char));
+
+ slot_size = (max - min) / num_slots;
+
+ rows = G_window_rows();
+ cols = G_window_cols();
+
+ get_slot_counts(infile);
+
+ bins = G_calloc(num_quant, sizeof(struct bin));
+ initialize_bins();
+ G_free(slots);
+
+ values = G_calloc(num_values, sizeof(DCELL));
+ fill_bins(infile);
+
+ G_close_cell(infile);
+ G_free(slot_bins);
+
+ sort_bins();
+ compute_quantiles();
+
+ return(EXIT_SUCCESS);
+}
+
More information about the grass-commit
mailing list