[GRASS-SVN] r41230 - in grass/trunk/raster: . r.resamp.filter
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Feb 28 13:42:35 EST 2010
Author: glynn
Date: 2010-02-28 13:42:34 -0500 (Sun, 28 Feb 2010)
New Revision: 41230
Added:
grass/trunk/raster/r.resamp.filter/
grass/trunk/raster/r.resamp.filter/Makefile
grass/trunk/raster/r.resamp.filter/main.c
grass/trunk/raster/r.resamp.filter/r.resamp.filter.html
Log:
Add r.resamp.filter
Property changes on: grass/trunk/raster/r.resamp.filter
___________________________________________________________________
Added: svn:ignore
+ OBJ.*
Added: grass/trunk/raster/r.resamp.filter/Makefile
===================================================================
--- grass/trunk/raster/r.resamp.filter/Makefile (rev 0)
+++ grass/trunk/raster/r.resamp.filter/Makefile 2010-02-28 18:42:34 UTC (rev 41230)
@@ -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
Added: grass/trunk/raster/r.resamp.filter/main.c
===================================================================
--- grass/trunk/raster/r.resamp.filter/main.c (rev 0)
+++ grass/trunk/raster/r.resamp.filter/main.c 2010-02-28 18:42:34 UTC (rev 41230)
@@ -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;
+}
Added: grass/trunk/raster/r.resamp.filter/r.resamp.filter.html
===================================================================
--- grass/trunk/raster/r.resamp.filter/r.resamp.filter.html (rev 0)
+++ grass/trunk/raster/r.resamp.filter/r.resamp.filter.html 2010-02-28 18:42:34 UTC (rev 41230)
@@ -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>
More information about the grass-commit
mailing list