[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