Add r.resamp.filter

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@41230 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.resamp.filter/Makefile

@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.resamp.filter
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

raster/r.resamp.filter/main.c

@@ -0,0 +1,566 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.resamp.filter
+ * AUTHOR(S):    Glynn Clements <glynn gclements.plus.com>
+ * PURPOSE:      
+ * COPYRIGHT:    (C) 2010 by Glynn Clements and the GRASS Development Team
+ *
+ *               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 <limits.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+static double f_box(double x)
+{
+    return (x > 1) ? 0
+	: 1;
+}
+
+static double f_bartlett(double x)
+{
+    return (x > 1) ? 0
+	: 1 - x;
+}
+
+static double f_hermite(double x)
+{
+    return (x > 1) ? 0
+	: (2 * x - 3) * x * x + 1;
+}
+
+static double f_gauss(double x)
+{
+    return exp(-2 * x * x) * sqrt(2 / M_PI);
+}
+
+static double f_normal(double x)
+{
+    return f_gauss(x/2) / 2;
+}
+
+static double f_sinc(double x)
+{
+    return (x == 0) ? 1 : sin(M_PI * x) / (M_PI * x);
+}
+
+static double lanczos(double x, int a)
+{
+    return (x > a) ? 0
+	: f_sinc(x) * f_sinc(x / a);
+}
+
+static double f_lanczos1(double x)
+{
+    return lanczos(x, 1);
+}
+
+static double f_lanczos2(double x)
+{
+    return lanczos(x, 2);
+}
+
+static double f_lanczos3(double x)
+{
+    return lanczos(x, 3);
+}
+
+static double f_hann(double x)
+{
+    return cos(M_PI * x) / 2 + 0.5;
+}
+
+static double f_hamming(double x)
+{
+    return 0.46 * cos(M_PI * x) + 0.54;
+}
+
+
+static double f_blackman(double x)
+{
+    return cos(M_PI * x) / 2 + 0.08 * cos(2 * M_PI * x) + 0.42;
+}
+
+struct filter_type {
+    const char *name;
+    double (*func)(double);
+    int radius;
+};
+
+static const struct filter_type menu[] = {
+    {"box",       f_box,       1},
+    {"bartlett",  f_bartlett,  1},
+    {"gauss",     f_gauss,     0},
+    {"normal",    f_normal,    0},
+    {"hermite",   f_hermite,   1},
+    {"sinc",      f_sinc,      0},
+    {"lanczos1",  f_lanczos1,  1},
+    {"lanczos2",  f_lanczos2,  2},
+    {"lanczos3",  f_lanczos3,  3},
+    {"hann",      f_hann,      0},
+    {"hamming",   f_hamming,   0},
+    {"blackman",  f_blackman,  0},
+    {NULL},
+};
+
+static char *build_filter_list(void)
+{
+    char *buf = G_malloc(1024);
+    char *p = buf;
+    int i;
+
+    for (i = 0; menu[i].name; i++) {
+	const char *q;
+
+	if (i)
+	    *p++ = ',';
+	for (q = menu[i].name; *q; p++, q++)
+	    *p = *q;
+    }
+    *p = '\0';
+
+    return buf;
+}
+
+static const struct filter_type *find_method(const char *name)
+{
+    int i;
+
+    for (i = 0; menu[i].name; i++)
+	if (strcmp(menu[i].name, name) == 0)
+	    return &menu[i];
+
+    G_fatal_error(_("Filter <%s> not found"), name);
+
+    return NULL;
+}
+
+struct filter {
+    double (*func)(double);
+    double x_radius;
+    double y_radius;
+};
+
+#define MAX_FILTERS 8
+
+static int infile, outfile;
+static struct filter filters[MAX_FILTERS];
+static int num_filters;
+static int nulls;
+static struct Cell_head dst_w, src_w;
+static double f_x_radius, f_y_radius;
+static int row_scale, col_scale;
+static DCELL *inbuf;
+static DCELL *outbuf;
+static DCELL **bufs;
+static double *h_weights;
+static double *v_weights;
+static int *mapcol0, *mapcol1;
+static int *maprow0, *maprow1;
+
+static void make_h_weights(void)
+{
+    int col;
+
+    h_weights = G_malloc(dst_w.cols * col_scale * sizeof(double));
+    mapcol0 = G_malloc(dst_w.cols * sizeof(int));
+    mapcol1 = G_malloc(dst_w.cols * sizeof(int));
+
+    for (col = 0; col < dst_w.cols; col++) {
+	double dx = Rast_col_to_easting(col + 0.5, &dst_w);
+	double x0 = Rast_easting_to_col(dx - f_x_radius, &src_w);
+	double x1 = Rast_easting_to_col(dx + f_x_radius, &src_w);
+	int col0 = (int)floor(x0);
+	int col1 = (int)floor(x1) + 1;
+	int cols = col1 - col0;
+	int j;
+
+	mapcol0[col] = col0;
+	mapcol1[col] = col1;
+
+	for (j = 0; j < cols; j++) {
+	    double sx = Rast_col_to_easting(col0 + j + 0.5, &src_w);
+	    double r = fabs(sx - dx);
+	    double w = 1.0;
+	    int k;
+
+	    for (k = 0; k < num_filters; k++)
+		w *= (*filters[k].func)(r / filters[k].x_radius);
+
+	    h_weights[col * col_scale + j] = w;
+	}
+
+	for (j = cols; j < col_scale; j++)
+	    h_weights[col * col_scale + j] = 0;
+    }
+}
+
+static void make_v_weights(void)
+{
+    int row;
+
+    v_weights = G_malloc(dst_w.rows * row_scale * sizeof(double));
+    maprow0 = G_malloc(dst_w.rows * sizeof(int));
+    maprow1 = G_malloc(dst_w.rows * sizeof(int));
+
+    for (row = 0; row < dst_w.rows; row++) {
+	double dy = Rast_row_to_northing(row + 0.5, &dst_w);
+	double y0 = Rast_northing_to_row(dy + f_y_radius, &src_w);
+	double y1 = Rast_northing_to_row(dy - f_y_radius, &src_w);
+	int row0 = (int)floor(y0);
+	int row1 = (int)floor(y1) + 1;
+	int rows = row1 - row0;
+	int i;
+
+	maprow0[row] = row0;
+	maprow1[row] = row1;
+
+	for (i = 0; i < rows; i++) {
+	    double sy = Rast_row_to_northing(row0 + i + 0.5, &src_w);
+	    double r = fabs(sy - dy);
+	    double w = 1.0;
+	    int k;
+
+	    for (k = 0; k < num_filters; k++)
+		w *= (*filters[k].func)(r / filters[k].y_radius);
+
+	    v_weights[row * row_scale + i] = w;
+	}
+
+	for (i = rows; i < row_scale; i++)
+	    v_weights[row * row_scale + i] = 0;
+    }
+}
+
+static void h_filter(DCELL *dst, const DCELL *src)
+{
+    int col;
+
+    for (col = 0; col < dst_w.cols; col++) {
+	int col0 = mapcol0[col];
+	int col1 = mapcol1[col];
+	int cols = col1 - col0;
+	double numer = 0.0;
+	double denom = 0.0;
+	int null = 0;
+	int j;
+
+	for (j = 0; j < cols; j++) {
+	    double w = h_weights[col * col_scale + j];
+	    const DCELL *c = &src[col0 + j];
+
+	    if (Rast_is_d_null_value(c)) {
+		if (nulls) {
+		    null = 1;
+		    break;
+		}
+	    }
+	    else {
+		numer += w * (*c);
+		denom += w;
+	    }
+	}
+
+	if (null || denom == 0)
+	    Rast_set_d_null_value(&dst[col], 1);
+	else
+	    dst[col] = numer / denom;
+    }
+}
+
+static void v_filter(DCELL *dst, DCELL **src, int row, int rows)
+{
+    int col;
+
+    for (col = 0; col < dst_w.cols; col++) {
+	double numer = 0.0;
+	double denom = 0.0;
+	int null = 0;
+	int i;
+
+	for (i = 0; i < rows; i++) {
+	    double w = v_weights[row * row_scale + i];
+	    const DCELL *c = &src[i][col];
+
+	    if (Rast_is_d_null_value(c)) {
+		if (nulls) {
+		    null = 1;
+		    break;
+		}
+	    }
+	    else {
+		numer += w * (*c);
+		denom += w;
+	    }
+	}
+
+	if (null || denom == 0)
+	    Rast_set_d_null_value(&dst[col], 1);
+	else
+	    dst[col] = numer / denom;
+    }
+}
+
+static void filter(void)
+{
+    int cur_row;
+    int num_rows = 0;
+    int row;
+
+    make_h_weights();
+    make_v_weights();
+
+    for (row = 0; row < dst_w.rows; row++) {
+	int row0 = maprow0[row];
+	int row1 = maprow1[row];
+	int rows = row1 - row0;
+	int i;
+
+	G_percent(row, dst_w.rows, 2);
+
+	if (row0 >= cur_row && row0 < cur_row + num_rows) {
+	    int m = row0 - cur_row;
+	    int n = cur_row + num_rows - row0;
+	    int i;
+
+	    for (i = 0; i < n; i++) {
+		DCELL *tmp = bufs[i];
+		bufs[i] = bufs[m + i];
+		bufs[m + i] = tmp;
+	    }
+
+	    cur_row = row0;
+	    num_rows = n;
+	}
+	else {
+	    cur_row = row0;
+	    num_rows = 0;
+	}
+
+	Rast_set_window(&src_w);
+
+	for (i = num_rows; i < rows; i++) {
+	    G_debug(5, "read: %p = %d", bufs[i], row0 + i);
+	    Rast_get_d_row(infile, inbuf, row0 + i);
+	    h_filter(bufs[i], inbuf);
+	}
+
+	num_rows = rows;
+
+	v_filter(outbuf, bufs, row, rows);
+
+	Rast_set_window(&dst_w);
+	Rast_put_d_row(outfile, outbuf);
+	G_debug(5, "write: %d", row);
+    }
+
+    G_percent(dst_w.rows, dst_w.rows, 2);
+}
+
+int main(int argc, char *argv[])
+{
+    struct GModule *module;
+    struct
+    {
+	struct Option *rastin, *rastout, *method,
+	    *radius, *x_radius, *y_radius;
+    } parm;
+    struct
+    {
+	struct Flag *nulls;
+    } flag;
+    char title[64];
+    int i;
+
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("resample"));
+    module->description =
+	_("Resamples raster map layers using an analytic kernel.");
+
+    parm.rastin = G_define_standard_option(G_OPT_R_INPUT);
+
+    parm.rastout = G_define_standard_option(G_OPT_R_OUTPUT);
+
+    parm.method = G_define_option();
+    parm.method->key = "filter";
+    parm.method->type = TYPE_STRING;
+    parm.method->required = YES;
+    parm.method->multiple = YES;
+    parm.method->description = _("Filter kernel(s)");
+    parm.method->options = build_filter_list();
+
+    parm.radius = G_define_option();
+    parm.radius->key = "radius";
+    parm.radius->type = TYPE_DOUBLE;
+    parm.radius->required = NO;
+    parm.radius->multiple = YES;
+    parm.radius->description = _("Filter radius");
+
+    parm.x_radius = G_define_option();
+    parm.x_radius->key = "x_radius";
+    parm.x_radius->type = TYPE_DOUBLE;
+    parm.x_radius->required = NO;
+    parm.x_radius->multiple = YES;
+    parm.x_radius->description = _("Filter radius (horizontal)");
+
+    parm.y_radius = G_define_option();
+    parm.y_radius->key = "y_radius";
+    parm.y_radius->type = TYPE_DOUBLE;
+    parm.y_radius->required = NO;
+    parm.y_radius->multiple = YES;
+    parm.y_radius->description = _("Filter radius (vertical)");
+
+    flag.nulls = G_define_flag();
+    flag.nulls->key = 'n';
+    flag.nulls->description = _("Propagate NULLs");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    if (parm.radius->answer) {
+	if (parm.x_radius->answer || parm.y_radius->answer)
+	    G_fatal_error(_("radius= and x_radius=/y_radius= are mutually exclusive"));
+    }
+    else {
+	if (!parm.x_radius->answer && !parm.y_radius->answer)
+	    G_fatal_error(_("Either radius= or x_radius=/y_radius= required"));
+	if (!parm.x_radius->answer || !parm.y_radius->answer)
+	    G_fatal_error(_("Both x_radius= and y_radius= required"));
+    }
+
+    nulls = flag.nulls->answer;
+
+    f_x_radius = f_y_radius = 1e100;
+
+    for (i = 0; ; i++) {
+	const char *filter_arg = parm.method->answers[i];
+	const char *x_radius_arg = parm.radius->answer
+	    ? parm.radius->answers[i]
+	    : parm.x_radius->answers[i];
+	const char *y_radius_arg = parm.radius->answer
+	    ? parm.radius->answers[i]
+	    : parm.y_radius->answers[i];
+	const struct filter_type *type;
+	struct filter *filter;
+
+	if (!filter_arg && !x_radius_arg && !y_radius_arg)
+	    break;
+
+	if (!filter_arg || !x_radius_arg || !y_radius_arg)
+	    G_fatal_error(_("Differing number of values for filter= and [xy_]radius="));
+
+	if (num_filters >= MAX_FILTERS)
+	    G_fatal_error(_("Too many filters (max: %d)"), MAX_FILTERS);
+
+	filter = &filters[num_filters++];
+	type = find_method(filter_arg);
+	filter->func = type->func;
+	filter->x_radius = fabs(atoi(x_radius_arg));
+	filter->y_radius = fabs(atoi(y_radius_arg));
+	if (type->radius) {
+	    double rx = type->radius * filter->x_radius;
+	    double ry = type->radius * filter->y_radius;
+	    if (rx < f_x_radius)
+		f_x_radius = rx;
+	    if (ry < f_y_radius)
+		f_y_radius = ry;
+	}
+    }
+
+    if (f_x_radius > 1e99 || f_y_radius > 1e99)
+	G_fatal_error(_("At least one filter must be finite"));
+
+    G_get_set_window(&dst_w);
+
+    /* set window to old map */
+    Rast_get_cellhd(parm.rastin->answer, "", &src_w);
+
+    /* enlarge source window */
+    {
+	double y0 = Rast_row_to_northing(0.5, &dst_w);
+	double y1 = Rast_row_to_northing(dst_w.rows - 0.5, &dst_w);
+	double x0 = Rast_col_to_easting(0.5, &dst_w);
+	double x1 = Rast_col_to_easting(dst_w.cols - 0.5, &dst_w);
+	int r0 = (int)floor(Rast_northing_to_row(y0 + f_y_radius, &src_w) - 0.1);
+	int r1 = (int)ceil(Rast_northing_to_row(y1 - f_y_radius, &src_w) + 0.1);
+	int c0 = (int)floor(Rast_easting_to_col(x0 - f_x_radius, &src_w) - 0.1);
+	int c1 = (int)ceil(Rast_easting_to_col(x1 + f_x_radius, &src_w) + 0.1);
+
+	src_w.south -= src_w.ns_res * (r1 - src_w.rows);
+	src_w.north += src_w.ns_res * (-r0);
+	src_w.west -= src_w.ew_res * (-c0);
+	src_w.east += src_w.ew_res * (c1 - src_w.cols);
+	src_w.rows = r1 - r0;
+	src_w.cols = c1 - c0;
+    }
+
+    row_scale = 2 + 2 * ceil(f_y_radius / src_w.ns_res);
+    col_scale = 2 + 2 * ceil(f_x_radius / src_w.ew_res);
+
+    /* allocate buffers for intermediate rows */
+    bufs = G_malloc(row_scale * sizeof(DCELL *));
+    for (i = 0; i < row_scale; i++)
+	bufs[i] = Rast_allocate_d_buf();
+
+    Rast_set_window(&src_w);
+
+    infile = Rast_open_old(parm.rastin->answer, "");
+    inbuf = Rast_allocate_d_buf();
+
+    Rast_set_window(&dst_w);
+
+    outfile = Rast_open_new(parm.rastout->answer, DCELL_TYPE);
+    outbuf = Rast_allocate_d_buf();
+
+    /* prevent complaints about window changes */
+    G_suppress_warnings(1);
+
+    filter();
+
+    Rast_close(infile);
+    Rast_close(outfile);
+
+    /* record map metadata/history info */
+    sprintf(title, "Filter resample by %s", parm.method->answer);
+    Rast_put_cell_title(parm.rastout->answer, title);
+
+    {
+	struct History history;
+	char buf_nsres[100], buf_ewres[100];
+
+	Rast_short_history(parm.rastout->answer, "raster", &history);
+	strncpy(history.datsrc_1, parm.rastin->answer, RECORD_LEN);
+	history.datsrc_1[RECORD_LEN - 1] = '\0';	/* strncpy() doesn't null terminate if maxfill */
+	G_format_resolution(src_w.ns_res, buf_nsres, src_w.proj);
+	G_format_resolution(src_w.ew_res, buf_ewres, src_w.proj);
+	sprintf(history.datsrc_2, "Source map NS res: %s   EW res: %s", buf_nsres,
+		buf_ewres);
+	Rast_command_history(&history);
+	Rast_write_history(parm.rastout->answer, &history);
+    }
+
+    /* copy color table from source map */
+    {
+	struct Colors colors;
+
+	if (Rast_read_colors(parm.rastin->answer, "", &colors) < 0)
+	    G_fatal_error(_("Unable to read color table for %s"),
+			  parm.rastin->answer);
+	Rast_mark_colors_as_fp(&colors);
+	Rast_write_colors(parm.rastout->answer, G_mapset(), &colors);
+    }
+
+    return EXIT_SUCCESS;
+}

raster/r.resamp.filter/r.resamp.filter.html

@@ -0,0 +1,53 @@
+<h2>DESCRIPTION</h2>
+
+<p>
+<em>r.resamp.filter</em> resamples an input raster, filtering the
+input with an analytic kernel.
+</p>
+<p>
+All of the kernels specified by the filter= option are multiplied
+together. Typical usage will use either a single kernel or an infinite
+kernel along with a finite window.
+</p>
+
+<h2>NOTES</h2>
+
+<p>
+Resampling modules (<em>r.resample, r.resamp.stats, r.resamp.interp,
+r.resamp.rst, r.resamp.filter</em>) resample the map to match the
+current region settings.
+</p>
+ 
+<p>
+When using a kernel which can have negative values (sinc, Lanczos),
+the <em>-n</em> flag should be used. Otherwise, extreme values can
+arise due to the total weight being close (or even equal) to zero.
+</p>
+ 
+<p>
+Kernels with infinite extent (Gauss, normal, sinc, Hann, Hamming,
+Blackman) must be used in conjunction with a finite windowing function
+(box, Bartlett, Hermite, Lanczos)
+</p>
+
+<p>
+For longitude-latitude locations, the interpolation algorithm is based on
+degree fractions, not on the absolute distances between cell centers.  Any
+attempt to implement the latter would violate the integrity of the
+interpolation method.
+</p>
+
+
+<h2>SEE ALSO</h2>
+
+<em><a href="g.region.html">g.region</a></em>,
+<em><a href="r.resample.html">r.resample</a></em>,
+<em><a href="r.resamp.rst.html">r.resamp.rst</a></em>
+<em><a href="r.resamp.stats.html">r.resamp.stats</a></em>
+
+<h2>AUTHOR</h2>
+
+Glynn Clements
+
+<p>
+<i>Last changed: $Date: 2008-09-26 07:16:02 +0100 (Fri, 26 Sep 2008) $</i>