[GRASS-SVN] r69863 - grass-addons/grass7/raster/r.resamp.tps
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Nov 21 12:49:12 PST 2016
Author: mmetz
Date: 2016-11-21 12:49:11 -0800 (Mon, 21 Nov 2016)
New Revision: 69863
Added:
grass-addons/grass7/raster/r.resamp.tps/r.resamp.tps.html
Removed:
grass-addons/grass7/raster/r.resamp.tps/v.surf.tps.html
grass-addons/grass7/raster/r.resamp.tps/v_surf_tps.png
Modified:
grass-addons/grass7/raster/r.resamp.tps/Makefile
grass-addons/grass7/raster/r.resamp.tps/main.c
grass-addons/grass7/raster/r.resamp.tps/tps.c
grass-addons/grass7/raster/r.resamp.tps/tps.h
Log:
new grass-addon r.resamp.tps
Modified: grass-addons/grass7/raster/r.resamp.tps/Makefile
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/Makefile 2016-11-21 20:46:28 UTC (rev 69862)
+++ grass-addons/grass7/raster/r.resamp.tps/Makefile 2016-11-21 20:49:11 UTC (rev 69863)
@@ -1,11 +1,9 @@
MODULE_TOPDIR = ../..
-PGM = v.surf.tps
+PGM = r.resamp.tps
-LIBES = $(VECTORLIB) $(DBMILIB) $(RASTERLIB) $(BTREE2LIB) $(SEGMENTLIB) $(GISLIB) $(MATHLIB)
-DEPENDENCIES = $(VECTORDEP) $(DBMIDEP) $(RASTERDEP) $(BTREE2DEP) $(SEGMENTDEP) $(GISDEP)
-EXTRA_INC = $(VECT_INC)
-EXTRA_CFLAGS = $(VECT_CFLAGS)
+LIBES = $(RASTERLIB) $(BTREE2LIB) $(SEGMENTLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(BTREE2DEP) $(SEGMENTDEP) $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Modified: grass-addons/grass7/raster/r.resamp.tps/main.c
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/main.c 2016-11-21 20:46:28 UTC (rev 69862)
+++ grass-addons/grass7/raster/r.resamp.tps/main.c 2016-11-21 20:49:11 UTC (rev 69863)
@@ -1,7 +1,7 @@
/**********************************************************************
*
- * MODULE: v.surf.tps
+ * MODULE: r.resamp.tps
*
* AUTHOR(S): Markus Metz
*
@@ -20,86 +20,49 @@
#include <string.h>
#include <math.h>
#include <grass/raster.h>
-#include <grass/vector.h>
+#include <grass/segment.h>
#include <grass/glocale.h>
#include "tps.h"
-static int cmp_pnts(const void *a, const void *b)
-{
- struct tps_pnt *aa = (struct tps_pnt *)a;
- struct tps_pnt *bb = (struct tps_pnt *)b;
-
- if (aa->r == bb->r)
- return (aa->c - bb->c);
-
- return (aa->r - bb->r);
-}
-
int main(int argc, char *argv[])
{
- int i, j;
- const char *mapset, *outname;
- int out_fd, mask_fd;
+ int i;
+ const char *outname;
+ int in_fd, out_fd, *var_fd;
- struct Map_info In;
struct History history;
+ struct Colors colors;
struct GModule *module;
- struct Option *in_opt, *var_opt, *out_opt, *minpnts_opt, *thin_opt,
- *reg_opt, *ov_opt, *dfield_opt, *col_opt, *mask_opt,
- *mem_opt;
+ struct Option *in_opt, *ivar_opt, *ovar_opt, *out_opt, *minpnts_opt,
+ *reg_opt, *ov_opt, *mask_opt, *mem_opt;
struct Flag *c_flag;
- struct Cell_head dstwindow;
+ struct Cell_head cellhd, src, dst;
- int cur_row;
- int n_vars;
- int n_points, min_points;
- int use_z;
- int cat;
- int val_field;
- char *val_column;
- double val;
- struct field_info *Fi;
- dbDriver *driver;
- dbColumn *dbcolumn;
- dbValue dbval;
- int sqltype = 0, ctype = 0;
+ int n_ivars, n_ovars, n_vars;
+ off_t n_points;
+ int min_points;
- struct line_pnts *Points;
- struct line_cats *Cats;
- struct tps_pnt *pnts;
- int type, cnt;
- double valsum;
int r, c, nrows, ncols;
- int *var_fd;
- DCELL **dbuf;
+ DCELL **dbuf, *dval;
double regularization, overlap;
- double pthin;
- int segs_mb;
+ SEGMENT in_seg, var_seg, out_seg;
+ int insize, varsize;
+ double segsize;
+ int segs_mb, nsegs;
/*----------------------------------------------------------------*/
/* Options declarations */
module = G_define_module();
- G_add_keyword(_("vector"));
+ G_add_keyword(_("raster"));
G_add_keyword(_("surface"));
G_add_keyword(_("interpolation"));
G_add_keyword(_("TPS"));
module->description =
_("Performs thin plate spline interpolation with regularization and covariables.");
- in_opt = G_define_standard_option(G_OPT_V_INPUT);
- in_opt->label = _("Name of input vector point map");
-
- dfield_opt = G_define_standard_option(G_OPT_V_FIELD);
- dfield_opt->guisection = _("Settings");
+ in_opt = G_define_standard_option(G_OPT_R_INPUT);
- col_opt = G_define_standard_option(G_OPT_DB_COLUMN);
- col_opt->required = NO;
- col_opt->label =
- _("Name of the attribute column with values to be used for interpolation");
- col_opt->description = _("If not given, z-coordinates are used.");
- col_opt->guisection = _("Settings");
-
reg_opt = G_define_option();
reg_opt->key = "smooth";
reg_opt->type = TYPE_DOUBLE;
@@ -113,7 +76,7 @@
ov_opt->key = "overlap";
ov_opt->type = TYPE_DOUBLE;
ov_opt->required = NO;
- ov_opt->answer = "0.1";
+ ov_opt->answer = "0.8";
ov_opt->label =
_("Overlap factor <= 1");
ov_opt->description =
@@ -129,23 +92,19 @@
_("Minimum number of points to use for TPS interpolation");
minpnts_opt->guisection = _("Settings");
- var_opt = G_define_standard_option(G_OPT_R_INPUTS);
- var_opt->key = "covars";
- var_opt->required = NO;
- var_opt->label =
- _("Name of input raster map(s) to use as covariables");
- var_opt->guisection = _("Settings");
+ ivar_opt = G_define_standard_option(G_OPT_R_INPUTS);
+ ivar_opt->key = "icovars";
+ ivar_opt->required = NO;
+ ivar_opt->label =
+ _("Name of input raster map(s) to use as covariables matching the input raster");
+ ivar_opt->guisection = _("Settings");
- thin_opt = G_define_option();
- thin_opt->key = "thin";
- thin_opt->type = TYPE_DOUBLE;
- thin_opt->required = NO;
- thin_opt->answer = "1.5";
- thin_opt->label =
- _("Point cloud thinning factor in number of cells of the current region");
- thin_opt->description =
- _("Minimum distance between neighboring points for local TPS interpolation");
- minpnts_opt->guisection = _("Settings");
+ ovar_opt = G_define_standard_option(G_OPT_R_INPUTS);
+ ovar_opt->key = "ocovars";
+ ovar_opt->required = NO;
+ ovar_opt->label =
+ _("Name of input raster map(s) to use as covariables matching the current region");
+ ovar_opt->guisection = _("Settings");
out_opt = G_define_standard_option(G_OPT_R_OUTPUT);
out_opt->key = "output";
@@ -176,273 +135,185 @@
outname = out_opt->answer;
- /* open input vector without topology */
- if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
- G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
+ n_ivars = 0;
+ if (ivar_opt->answer) {
+ while (ivar_opt->answers[n_ivars])
+ n_ivars++;
+ }
- Vect_set_open_level(1); /* WITHOUT TOPOLOGY */
- if (1 > Vect_open_old(&In, in_opt->answer, mapset))
- G_fatal_error(_("Unable to open vector map <%s>"),
- in_opt->answer);
+ n_ovars = 0;
+ if (ovar_opt->answer) {
+ while (ovar_opt->answers[n_ovars])
+ n_ovars++;
+ }
- val_field = 0; /* assume 3D input */
- val_column = col_opt->answer;
-
- use_z = !val_column && Vect_is_3d(&In);
-
- if (Vect_is_3d(&In)) {
- if (!use_z)
- G_verbose_message(_("Input is 3D: using attribute values instead of z-coordinates for approximation"));
- else
- G_verbose_message(_("Input is 3D: using z-coordinates for approximation"));
+ if (n_ivars != n_ovars) {
+ G_fatal_error(_("Number of covariables matching the input raster "
+ "must be equal to the number of covariables "
+ "matching the current region"));
}
- else { /* 2D */
- if (!val_column)
- G_fatal_error(_("Input vector map is 2D. Parameter <%s> required."), col_opt->key);
- }
+ n_vars = n_ovars;
- if (!use_z) {
- val_field = Vect_get_field_number(&In, dfield_opt->answer);
+ if (!G_find_raster2(in_opt->answer, ""))
+ G_fatal_error(_("Input map <%s> not found"), in_opt->answer);
+
+ for (i = 0; i < n_vars; i++) {
+ if (!G_find_raster2(ivar_opt->answers[i], ""))
+ G_fatal_error(_("Input map <%s> not found"), ivar_opt->answers[i]);
+ if (!G_find_raster2(ovar_opt->answers[i], ""))
+ G_fatal_error(_("Input map <%s> not found"), ovar_opt->answers[i]);
}
- n_vars = 0;
- if (var_opt->answer) {
- while (var_opt->answers[n_vars])
- n_vars++;
+ if (mask_opt->answer) {
+ if (!G_find_raster2(mask_opt->answer, ""))
+ G_fatal_error(_("Mask map <%s> not found"), mask_opt->answer);
}
- Points = Vect_new_line_struct();
- Cats = Vect_new_cats_struct();
+ var_fd = NULL;
+ if (n_vars)
+ var_fd = G_malloc(n_vars * sizeof(int));
+ dbuf = G_malloc((1 + n_vars) * sizeof(DCELL *));
+ dval = G_malloc((1 + n_vars) * sizeof(DCELL));
- Rast_get_window(&dstwindow);
- nrows = dstwindow.rows;
- ncols = dstwindow.cols;
+ Rast_get_window(&dst);
- /* count points */
- Vect_rewind(&In);
- n_points = 0;
- while (1) {
- /* register line */
- type = Vect_read_next_line(&In, Points, Cats);
+ /* get cellhd of input */
+ Rast_get_cellhd(in_opt->answer, "", &cellhd);
- /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
- if (type == -1) {
- G_fatal_error(_("Unable to read vector map"));
- }
- else if (type == -2) {
- break;
- }
+ /* align dst to cellhd */
+ src = dst;
+ Rast_align_window(&src, &cellhd);
- if (type & GV_POINTS) {
- n_points++;
-
- if (n_vars) {
- r = (dstwindow.north - Points->y[0]) / dstwindow.ns_res;
- c = (Points->x[0] - dstwindow.west) / dstwindow.ew_res;
+ /* open segment structures for input and output */
- if (r < 0 || r >= nrows || c < 0 || c >= ncols)
- n_points--;
- }
- }
- }
- if (!n_points)
- G_fatal_error(_("No points in input vector <%s>"), in_opt->answer);
+ /* set input window */
+ Rast_set_window(&src);
+ nrows = src.rows;
+ ncols = src.cols;
- /* check memory */
- pnts = G_malloc(n_points * sizeof(struct tps_pnt));
- if (!pnts)
- G_fatal_error(_("Too many points, out of memory"));
+ segs_mb = atoi(mem_opt->answer);
+ if (segs_mb < 10)
+ segs_mb = 10;
- if (n_vars) {
- pnts[0].vars = G_malloc(sizeof(double) * n_points * n_vars);
- if (!pnts[0].vars)
- G_fatal_error(_("Too many points, out of memory"));
- G_free(pnts[0].vars);
- }
+ segsize = (1 + n_vars) * sizeof(DCELL) /* input */
+ + n_vars * sizeof(DCELL) /* covariables */
+ + sizeof(struct tps_out); /* output */
- /* load points */
- driver = NULL;
- Fi = NULL;
- ctype = 0;
- if (!use_z) {
- Fi = Vect_get_field(&In, val_field);
- if (Fi == NULL)
- G_fatal_error(_("Cannot read layer info"));
+ segsize = segsize * 64. * 64. / (1024. * 1024.);
+ nsegs = segs_mb / segsize;
- driver = db_start_driver_open_database(Fi->driver, Fi->database);
+ /* load input raster and corresponding covariables */
+ G_message(_("Loading input..."));
- if (db_get_column(driver, Fi->table, val_column, &dbcolumn) != DB_OK) {
- db_close_database_shutdown_driver(driver);
- G_fatal_error(_("Unable to get column <%s>"), val_column);
- }
- sqltype = db_get_column_sqltype(dbcolumn);
- ctype = db_sqltype_to_Ctype(sqltype);
-
- if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE) {
- db_close_database_shutdown_driver(driver);
- G_fatal_error(_("Column <%s> must be numeric"), val_column);
- }
+ insize = (1 + n_vars) * sizeof(DCELL);
+ if (Segment_open(&in_seg, G_tempfile(), nrows, ncols, 32, 32,
+ insize, nsegs) != 1) {
+ G_fatal_error("Unable to create input temporary files");
}
- G_message(_("Loading %d points..."), n_points);
- Vect_rewind(&In);
- i = 0;
- j = 0;
- while (1) {
-
- i++;
+ in_fd = Rast_open_old(in_opt->answer, "");
+ for (i = 0; i < n_vars; i++)
+ var_fd[i] = Rast_open_old(ivar_opt->answers[i], "");
- /* register line */
- type = Vect_read_next_line(&In, Points, Cats);
+ for (i = 0; i <= n_vars; i++)
+ dbuf[i] = Rast_allocate_d_buf();
+
+ n_points = 0;
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 5);
+ Rast_get_row(in_fd, dbuf[0], r, DCELL_TYPE);
+ for (i = 0; i < n_vars; i++)
+ Rast_get_row(var_fd[i], dbuf[i + 1], r, DCELL_TYPE);
+ for (c = 0; c < ncols; c++) {
- /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
- if (type == -1) {
- G_fatal_error(_("Unable to read vector map"));
- }
- else if (type == -2) {
- break;
- }
-
- if (!(type & GV_POINTS))
- continue;
-
- Rast_set_d_null_value(&val, 1);
- if (driver) {
- if (!Vect_cat_get(Cats, val_field, &cat))
- G_fatal_error(_("No category for point %d in layer %d"), i, val_field);
-
- db_select_value(driver, Fi->table, Fi->key, cat, val_column,
- &dbval);
- switch (ctype)
- {
- case DB_C_TYPE_INT: {
- val = db_get_value_int(&dbval);
- break;
+ dval[0] = dbuf[0][c];
+ if (!Rast_is_d_null_value(&dval[0])) {
+
+ n_points++;
+ for (i = 1; i <= n_vars; i++) {
+ dval[i] = dbuf[i][c];
+ if (Rast_is_d_null_value(&dval[i])) {
+ Rast_set_d_null_value(dval, 1 + n_vars);
+ n_points--;
+ break;
+ }
+ }
}
- case DB_C_TYPE_DOUBLE: {
- val = db_get_value_double(&dbval);
- break;
+ else if (n_vars) {
+ Rast_set_d_null_value(&dval[1], n_vars);
}
- default: {
- Rast_set_d_null_value(&val, 1);
- break;
- }
- }
+
+ if (Segment_put(&in_seg, (void *)dval, r, c) != 1)
+ G_fatal_error(_("Unable to write to temporary file"));
}
- else {
- val = Points->z[0];
- }
- r = pnts[j].r = (dstwindow.north - Points->y[0]) / dstwindow.ns_res;
- c = pnts[j].c = (Points->x[0] - dstwindow.west) / dstwindow.ew_res;
-
- if (n_vars && (r < 0 || r >= nrows || c < 0 || c >= ncols)) {
- Rast_set_d_null_value(&val, 1);
- }
-
- if (!Rast_is_d_null_value(&val)) {
- /* add point */
- G_percent(j, n_points, 4);
-
- pnts[j].val = val;
- j++;
- }
}
- G_percent(1, 1, 1);
+ Segment_flush(&in_seg);
+ G_percent(r, nrows, 5);
- if (driver)
- db_close_database_shutdown_driver(driver);
- Vect_close(&In);
+ Rast_close(in_fd);
+ for (i = 0; i < n_vars; i++)
+ Rast_close(var_fd[i]);
- if (j != n_points)
- n_points = j;
+ for (i = 0; i <= n_vars; i++)
+ G_free(dbuf[i]);
+ if (!n_points)
+ G_fatal_error(_("No valid input points"));
+
+ G_message(_("%ld input points"), n_points);
- /* sort points */
- qsort(pnts, n_points, sizeof(struct tps_pnt), cmp_pnts);
+ /* set output window */
+ Rast_set_window(&dst);
+ nrows = dst.rows;
+ ncols = dst.cols;
- /* combine duplicates */
- valsum = pnts[0].val;
- cnt = 1;
- j = 1;
- for (i = 1; i < n_points; i++) {
- if (pnts[i - 1].r == pnts[i].r && pnts[i - 1].c == pnts[i].c) {
- valsum += pnts[i].val;
- cnt++;
- }
- else {
- if (cnt > 1) {
- pnts[j - 1].val = valsum / cnt;
- valsum = pnts[i].val;
- cnt = 1;
- }
- pnts[j].r = pnts[i].r;
- pnts[j].c = pnts[i].c;
- pnts[j].val = pnts[i].val;
- j++;
- }
+ if (Segment_open(&out_seg, G_tempfile(), nrows, ncols, 64, 64,
+ sizeof(struct tps_out), nsegs) != 1) {
+ G_fatal_error("Unable to create input temporary files");
}
- if (j != n_points) {
- G_verbose_message(_("%d duplicate points pruned"), n_points - j);
- n_points = j;
- }
if (n_vars) {
- pnts[0].vars = G_malloc(sizeof(double) * n_points * n_vars);
- if (!pnts[0].vars)
- G_fatal_error(_("Too many points, out of memory"));
-
- for (i = 1; i < n_points; i++)
- pnts[i].vars = pnts[i - 1].vars + n_vars;
- }
- /* load covariables */
- var_fd = NULL;
- dbuf = NULL;
- if (n_vars) {
- int ii, isnull;
+ /* intialize output raster and load corresponding covariables */
+ G_message(_("Loading covariables for output..."));
+ varsize = (n_vars) * sizeof(DCELL);
+ if (Segment_open(&var_seg, G_tempfile(), nrows, ncols, 64, 64,
+ varsize, nsegs) != 1) {
+ G_fatal_error("Unable to create input temporary files");
+ }
- G_message(_("Loading %d covariables..."), n_vars);
+ /* open segment structure */
- cur_row = -1;
+ for (i = 0; i < n_vars; i++)
+ var_fd[i] = Rast_open_old(ovar_opt->answers[i], "");
- var_fd = G_malloc(n_vars * sizeof(int));
- dbuf = G_malloc(n_vars * sizeof(DCELL *));
- for (i = 0; var_opt->answers[i]; i++) {
- var_fd[i] = Rast_open_old(var_opt->answers[i], "");
+ for (i = 0; i < n_vars; i++)
dbuf[i] = Rast_allocate_d_buf();
- }
- ii = 0;
- for (i = 0; i < n_points; i++) {
- G_percent(i, n_points, 4);
- if (pnts[i].r >= 0 && pnts[i].r < nrows &&
- pnts[i].c >= 0 && pnts[i].c < ncols) {
- if (cur_row != pnts[i].r) {
- cur_row = pnts[i].r;
- for (j = 0; j < n_vars; j++)
- Rast_get_d_row(var_fd[j], dbuf[j], cur_row);
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 5);
+ for (i = 0; i < n_vars; i++)
+ Rast_get_row(var_fd[i], dbuf[i], r, DCELL_TYPE);
+ for (c = 0; c < ncols; c++) {
+
+ for (i = 0; i < n_vars; i++) {
+ dval[i] = dbuf[i][c];
+ if (Rast_is_d_null_value(&dval[i])) {
+ Rast_set_d_null_value(dval, n_vars);
+ break;
+ }
}
- isnull = 0;
- for (j = 0; j < n_vars; j++) {
- if (Rast_is_d_null_value(&dbuf[j][pnts[i].c]))
- isnull = 1;
- pnts[ii].vars[j] = dbuf[j][pnts[i].c];
- }
- if (!isnull) {
- pnts[ii].r = pnts[i].r;
- pnts[ii].c = pnts[i].c;
- pnts[ii].val = pnts[i].val;
- ii++;
- }
+ if (Segment_put(&var_seg, (void *)dval, r, c) != 1)
+ G_fatal_error(_("Unable to write to temporary file"));
}
}
- G_percent(1, 1, 1);
- n_points = ii;
- }
- if (n_vars) {
- for (i = 0; i < n_vars; i++) {
+ Segment_flush(&var_seg);
+ G_percent(r, nrows, 5);
+
+ for (i = 0; i < n_vars; i++)
+ Rast_close(var_fd[i]);
+
+ for (i = 0; i < n_vars; i++)
G_free(dbuf[i]);
- }
- G_free(dbuf);
}
out_fd = Rast_open_new(outname, DCELL_TYPE);
@@ -464,41 +335,17 @@
if (overlap > 1)
overlap = 1;
- mask_fd = -1;
- if (mask_opt->answer)
- mask_fd = Rast_open_old(mask_opt->answer, "");
-
- segs_mb = atoi(mem_opt->answer);
- if (segs_mb < 10)
- segs_mb = 10;
-
- pthin = atof(thin_opt->answer);
- if (pthin < 0)
- pthin = 0;
-
- if (n_points <= min_points) {
- if (global_tps(out_fd, var_fd, n_vars, mask_fd, pnts, n_points,
- regularization) != 1) {
- G_fatal_error(_("TPS interpolation failed"));
- }
+ if (local_tps(&in_seg, &var_seg, n_vars, &out_seg, out_fd,
+ mask_opt->answer, &src, &dst, n_points,
+ min_points, regularization, overlap,
+ c_flag->answer) != 1) {
+ G_fatal_error(_("TPS interpolation failed"));
}
- else {
- if (local_tps(out_fd, var_fd, n_vars, mask_fd, pnts, n_points,
- min_points, regularization, overlap, pthin,
- c_flag->answer, segs_mb) != 1) {
- G_fatal_error(_("TPS interpolation failed"));
- }
- }
- if (n_vars) {
- for (i = 0; i < n_vars; i++) {
- Rast_close(var_fd[i]);
- }
- G_free(var_fd);
- }
-
- if (mask_fd >= 0)
- Rast_close(mask_fd);
+ Segment_close(&in_seg);
+ if (n_vars)
+ Segment_close(&var_seg);
+ Segment_close(&out_seg);
/* write map history */
Rast_close(out_fd);
@@ -506,6 +353,9 @@
Rast_command_history(&history);
Rast_write_history(outname, &history);
+ if (Rast_read_colors(in_opt->answer, "", &colors) == 1)
+ Rast_write_colors(outname, G_mapset(), &colors);
+
G_done_msg(" ");
exit(EXIT_SUCCESS);
Copied: grass-addons/grass7/raster/r.resamp.tps/r.resamp.tps.html (from rev 69862, grass-addons/grass7/raster/r.resamp.tps/v.surf.tps.html)
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/r.resamp.tps.html (rev 0)
+++ grass-addons/grass7/raster/r.resamp.tps/r.resamp.tps.html 2016-11-21 20:49:11 UTC (rev 69863)
@@ -0,0 +1,76 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.resamp.tps</em> performs multivariate thin plate spline
+interpolation with regularization. The <b>input</b> is a raster
+map to be resampled to a higher resolution or where NULL cells need to
+be interpolated. Output is a raster map. Optionally, several
+raster maps can be specified to be used as covariables which will
+improve results in areas with few points. Raster maps to be used as
+covariables need to be provided separately matching the grid geometry
+of the <b>input</b> raster map with the <b>icovars</b> option and
+matching the grid geometry of the <b>output</b> raster map with the
+<b>ocovars</b> option. The module can be regarded as a combination of
+a multiple regression and spline interpolation.
+
+<p>
+The <b>min</b> options specifies the minimum number of points to be
+used for interpolation. <em>r.resamp.tps</em> always performs tiled
+local TPS interpolation. Tile sizes are variable and dependent on the
+extents of the <b>min</b> nearest neighbors when a new tile is generated.
+
+<p>
+The <b>smooth</b> option can be used to reduce the influence of the
+splines and increase the influence of the covariables. Without
+covariables, the resulting surface will be smoother. With covariables
+and a large smooting value, the resulting surface will be mainly
+determined by the multiple regression component.
+
+<p>
+The <b>overlap</b> option controls how much tiles are overlapping when
+the <b>min</b> option is smaller than the numer of input points.
+Tiling artefacts occur with low values for the <b>min</b> option and the
+<b>overlap</b> option. Increasing both options will reduce tiling
+artefacts but processing will take more time.
+
+<p>
+The module works best with evenly spaced points. In case of
+highly unevenly spaced points, e.g. remote sensing data with gaps due
+to cloud cover, the module will take a long time to finish. For data
+with large gaps, it is recommended to use first a different
+interpolation method and then optionally use <em>r.resamp.tps</em> with
+the <b>smooth</b> option to identify outliers (difference between the
+output of <em>r.resamp.tps</em> and the data interpolated with a
+different method).
+
+<p>
+The <b>memory</b> option controls only how much memory should be used
+for the covariables and the intermediate output. The data needed for
+TPS interpolation are always completely loaded to memory.
+
+
+<h2>REFERENCES</h2>
+
+<ul>
+ <li>Hutchinson MF, 1995, Interpolating mean rainfall using thin plate
+ smoothing splines. International Journal of Geographical Information
+ Systems, 9(4), pp. 385-403</li>
+ <li>Wahba G, 1990, Spline models for observational data. In CBMS-NSF
+ Regional Conference Series in Applied Mathematics. Philadelpia:
+ Society for Industrial and Applied Mathematics</li>
+</ul>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="v.surf.tps.html">v.surf.tps</a>,
+<a href="v.surf.rst.html">v.surf.rst</a>,
+<a href="v.surf.bspline.html">v.surf.rst</a>,
+<a href="v.surf.idw.html">v.surf.idw</a>
+</em>
+
+<h2>AUTHORS</h2>
+
+Markus Metz
+
+<p>
+<i>Last changed: $Date$</i>
Modified: grass-addons/grass7/raster/r.resamp.tps/tps.c
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/tps.c 2016-11-21 20:46:28 UTC (rev 69862)
+++ grass-addons/grass7/raster/r.resamp.tps/tps.c 2016-11-21 20:49:11 UTC (rev 69863)
@@ -4,14 +4,17 @@
#include <math.h>
#include <grass/gis.h>
#include <grass/raster.h>
-#include <grass/vector.h>
-#include <grass/glocale.h>
-#include <grass/kdtree.h>
#include <grass/segment.h>
+#include <grass/rbtree.h>
+#include <grass/glocale.h>
#include "tps.h"
#include "flag.h"
#include "rclist.h"
+#ifdef USE_RC
+#undef USE_RC
+#endif
+
static int solvemat(double **m, double a[], double B[], int n)
{
int i, j, i2, j2, imark;
@@ -85,34 +88,88 @@
return 1;
}
-static int load_tps_pnts(struct tps_pnt *pnts, int n_points, int n_vars,
+static int row_src2dst(int row, struct Cell_head *src, struct Cell_head *dst)
+{
+ return (dst->north - src->north + (row + 0.5) * src->ns_res) / dst->ns_res;
+}
+
+static int row_dst2src(int row, struct Cell_head *src, struct Cell_head *dst)
+{
+ return (src->north - dst->north + (row + 0.5) * dst->ns_res) / src->ns_res;
+}
+
+static int col_src2dst(int col, struct Cell_head *src, struct Cell_head *dst)
+{
+ return (src->west - dst->west + (col + 0.5) * src->ew_res) / dst->ew_res;
+}
+
+static int col_dst2src(int col, struct Cell_head *src, struct Cell_head *dst)
+{
+ return (dst->west - src->west + (col + 0.5) * dst->ew_res) / src->ew_res;
+}
+
+static int cmp_pnts(const void *first, const void *second)
+{
+ struct tps_pnt *a = (struct tps_pnt *)first;
+ struct tps_pnt *b = (struct tps_pnt *)second;
+
+ if (a->r == b->r)
+ return (a->c - b->c);
+
+ return (a->r - b->r);
+}
+
+static int load_tps_pnts(SEGMENT *in_seg, int n_vars,
+ struct tps_pnt *pnts, int n_points,
+ struct Cell_head *src, struct Cell_head *dst,
double regularization, double **m, double *a)
{
int i, j, k;
double dx, dy, dist, dist2, distsum, reg;
+ DCELL *dval;
+ dval = G_malloc((1 + n_vars) * sizeof(DCELL));
+
for (i = 0; i <= n_vars; i++) {
a[i] = 0.0;
- for (j = 0; j <= n_vars; j++) {
+ m[i][i] = 0.0;
+ for (j = 0; j < i; j++) {
m[i][j] = m[j][i] = 0.0;
}
}
distsum = 0;
for (i = 0; i < n_points; i++) {
+
+ Segment_get(in_seg, (void *)dval, pnts[i].r, pnts[i].c);
+
/* global */
- a[i + 1 + n_vars] = pnts[i].val;
+ a[i + 1 + n_vars] = dval[0];
m[0][i + 1 + n_vars] = m[i + 1 + n_vars][0] = 1.0;
for (j = 0; j < n_vars; j++) {
- m[j + 1][i + 1 + n_vars] = m[i + 1 + n_vars][j + 1] = pnts[i].vars[j];
+ m[j + 1][i + 1 + n_vars] = m[i + 1 + n_vars][j + 1] = dval[j + 1];
}
+ /* convert src r,c to dst r,c or to n,s */
+#ifdef USE_RC
+ pnts[i].r = row_src2dst(pnts[i].r, src, dst);
+ pnts[i].c = col_src2dst(pnts[i].c, src, dst);
+#else
+ pnts[i].r = src->north - (pnts[i].r + 0.5) * src->ns_res;
+ pnts[i].c = src->west + (pnts[i].c + 0.5) * src->ew_res;
+#endif
+
/* TPS */
for (k = 0; k <= i; k++) {
+
+#ifdef USE_RC
dx = (pnts[i].c -pnts[k].c) * 2.0;
dy = (pnts[i].r -pnts[k].r) * 2.0;
-
+#else
+ dx = (pnts[i].c -pnts[k].c);
+ dy = (pnts[i].r -pnts[k].r);
+#endif
dist2 = dx * dx + dy * dy;
dist = 0;
if (dist2 > 0)
@@ -131,7 +188,9 @@
}
}
}
-
+
+ G_free(dval);
+
if (regularization > 0) {
distsum /= (n_points * n_points);
reg = regularization * distsum * distsum;
@@ -143,137 +202,129 @@
return 1;
}
-int global_tps(int out_fd, int *var_fd, int n_vars, int mask_fd,
- struct tps_pnt *pnts, int n_points, double regularization)
+
+static int cmp_rc(const void *first, const void *second)
{
- int row, col, nrows, ncols;
- double **m, *a, *B;
- int i, j;
- int nalloc;
- double dx, dy, dist, dist2;
- DCELL **dbuf, result, *outbuf;
- CELL *maskbuf;
+ struct rc *a = (struct rc *)first, *b = (struct rc *)second;
- G_message(_("Global TPS interpolation with %d points..."), n_points);
+ if (a->row == b->row)
+ return (a->col - b->col);
- nrows = Rast_window_rows();
- ncols = Rast_window_cols();
+ return (a->row - b->row);
+}
- /* alloc */
- outbuf = Rast_allocate_d_buf();
- nalloc = n_points;
- a = G_malloc((nalloc + 1 + n_vars) * sizeof(double));
- B = G_malloc((nalloc + 1 + n_vars) * sizeof(double));
- m = G_malloc((nalloc + 1 + n_vars) * sizeof(double *));
- for (i = 0; i < (nalloc + 1 + n_vars); i++)
- m[i] = G_malloc((nalloc + 1 + n_vars) * sizeof(double));
+static int bfs_search_nn(FLAG *pnt_flag, struct Cell_head *src,
+ struct tps_pnt *cur_pnts, int max_points,
+ int row, int col,
+ int *rmin, int *rmax, int *cmin, int *cmax,
+ double *distmax)
+{
+ int nrows, ncols, rown, coln;
+ int nextr[8] = {0, -1, 0, 1, -1, -1, 1, 1};
+ int nextc[8] = {1, 0, -1, 0, 1, -1, -1, 1};
+ int n, found;
+ struct rc next, ngbr_rc;
+ struct rclist rilist;
+ struct RB_TREE *visited;
+ double dx, dy, dist;
- dbuf = NULL;
- if (n_vars) {
- dbuf = G_malloc(n_vars * sizeof(DCELL *));
- for (i = 0; i < n_vars; i++) {
- dbuf[i] = Rast_allocate_d_buf();
- }
- }
+ visited = rbtree_create(cmp_rc, sizeof(struct rc));
+ ngbr_rc.row = row;
+ ngbr_rc.col = col;
+ rbtree_insert(visited, &ngbr_rc);
- maskbuf = NULL;
- if (mask_fd >= 0)
- maskbuf = Rast_allocate_c_buf();
+ nrows = src->rows;
+ ncols = src->cols;
- /* load points to matrix */
- load_tps_pnts(pnts, n_points, n_vars, regularization, m, a);
+ found = 0;
+ *distmax = 0.0;
- /* solve */
- if (!solvemat(m, a, B, n_points + 1 + n_vars))
- G_fatal_error(_("Matrix is unsolvable"));
+ if (FLAG_GET(pnt_flag, row, col)) {
+ /* add to cur_pnts */
+ cur_pnts[found].r = row;
+ cur_pnts[found].c = col;
+ found++;
- for (row = 0; row < nrows; row++) {
- G_percent(row, nrows, 2);
+ if (*rmin > row)
+ *rmin = row;
+ if (*rmax < row)
+ *rmax = row;
+ if (*cmin > col)
+ *cmin = col;
+ if (*cmax < col)
+ *cmax = col;
+ }
- for (j = 0; j < n_vars; j++)
- Rast_get_d_row(var_fd[j], dbuf[j], row);
- if (maskbuf)
- Rast_get_c_row(mask_fd, maskbuf, row);
+ /* breadth-first search */
+ next.row = row;
+ next.col = col;
+ rclist_init(&rilist);
- for (col = 0; col < ncols; col++) {
+ do {
+ for (n = 0; n < 8; n++) {
+ rown = next.row + nextr[n];
+ coln = next.col + nextc[n];
+
+ if (rown < 0 || rown >= nrows || coln < 0 || coln >= ncols)
+ continue;
- if (maskbuf &&
- (Rast_is_c_null_value(&maskbuf[col]) || maskbuf[col] == 0)) {
+ ngbr_rc.row = rown;
+ ngbr_rc.col = coln;
- Rast_set_d_null_value(&outbuf[col], 1);
-
+ if (rbtree_find(visited, &ngbr_rc))
continue;
- }
- if (n_vars) {
- int isnull = 0;
+ rbtree_insert(visited, &ngbr_rc);
+ rclist_add(&rilist, rown, coln);
- for (j = 0; j < n_vars; j++) {
- if (Rast_is_d_null_value(&dbuf[j][col])) {
- isnull = 1;
- break;
- }
- }
- if (isnull) {
- Rast_set_d_null_value(&outbuf[col], 1);
- continue;
- }
- }
+ if (FLAG_GET(pnt_flag, rown, coln)) {
- result = B[0];
- for (j = 0; j < n_vars; j++)
- result += dbuf[j][col] * B[j + 1];
-
- for (i = 0; i < n_points; i++) {
- dx = (pnts[i].c - col) * 2.0;
- dy = (pnts[i].r - row) * 2.0;
+ dx = coln - col;
+ dy = rown - row;
- dist2 = dx * dx + dy * dy;
- dist = 0;
- if (dist2 > 0) {
- dist = dist2 * log(dist2) * 0.5;
- result += B[1 + n_vars + i] * dist;
- }
+ dist = dx * dx + dy * dy;
+
+ if (*distmax < dist)
+ *distmax = dist;
+
+ cur_pnts[found].r = rown;
+ cur_pnts[found].c = coln;
+ found++;
+
+ if (*rmin > rown)
+ *rmin = rown;
+ if (*rmax < rown)
+ *rmax = rown;
+ if (*cmin > coln)
+ *cmin = coln;
+ if (*cmax < coln)
+ *cmax = coln;
}
- outbuf[col] = result;
- }
- Rast_put_d_row(out_fd, outbuf);
- }
- G_percent(1, 1, 1);
- if (n_vars) {
- for (i = 0; i < n_vars; i++) {
- G_free(dbuf[i]);
+ if (found == max_points)
+ break;
}
- G_free(dbuf);
- }
- if (maskbuf)
- G_free(maskbuf);
- G_free(outbuf);
+ if (found == max_points)
+ break;
+ } while (rclist_drop(&rilist, &next)); /* while there are cells to check */
- return 1;
-}
-static int cmp_rc(const void *first, const void *second)
-{
- struct rc *a = (struct rc *)first, *b = (struct rc *)second;
+ rclist_destroy(&rilist);
+ rbtree_destroy(visited);
- if (a->row == b->row)
- return (a->col - b->col);
-
- return (a->row - b->row);
+ return found;
}
-static int bfs_search(SEGMENT *p_seg, FLAG *mask_flag, int *bfsid,
- int row, int col,
+static int bfs_search(FLAG *pnt_flag, struct Cell_head *src,
+ struct tps_pnt *cur_pnts,
+ int row, int col,
int max_points, double min_dist,
int rmin, int rmax, int cmin, int cmax)
{
int nrows, ncols, rown, coln;
- int nextr[8] = {-1, -1, -1, 0, 0, 1, 1, 1};
- int nextc[8] = {-1, 0, 1, -1, 1, -1, 0, 1};
+ int nextr[8] = {0, -1, 0, 1, -1, -1, 1, 1};
+ int nextc[8] = {1, 0, -1, 0, 1, -1, -1, 1};
int n, found;
- int pid;
struct rc next, ngbr_rc;
struct rclist rilist;
struct RB_TREE *visited;
@@ -285,8 +336,8 @@
ngbr_rc.col = col;
rbtree_insert(visited, &ngbr_rc);
- nrows = Rast_window_rows();
- ncols = Rast_window_cols();
+ nrows = src->rows;
+ ncols = src->cols;
/* breadth-first search */
next.row = row;
@@ -304,8 +355,6 @@
continue;
if (rown < rmin || rown > rmax || coln < cmin || coln > cmax)
continue;
- if ((FLAG_GET(mask_flag, rown, coln)))
- continue;
ngbr_rc.row = rown;
ngbr_rc.col = coln;
@@ -315,19 +364,18 @@
rbtree_insert(visited, &ngbr_rc);
- Segment_get(p_seg, (void *)&pid, rown, coln);
-
- if (pid < 0) {
+ if (!(FLAG_GET(pnt_flag, rown, coln))) {
rclist_add(&rilist, rown, coln);
}
- else if (rown >= rmin && rown <= rmax && coln >= cmin && coln <= cmax) {
+ else {
dx = coln - col;
dy = rown - row;
dist = dx * dx + dy * dy;
if (dist > min_dist) {
- bfsid[found] = pid;
+ cur_pnts[found].r = rown;
+ cur_pnts[found].c = coln;
found++;
}
}
@@ -338,17 +386,17 @@
break;
} while (rclist_drop(&rilist, &next)); /* while there are cells to check */
-
rclist_destroy(&rilist);
rbtree_destroy(visited);
return found;
}
-int local_tps(int out_fd, int *var_fd, int n_vars, int mask_fd,
- struct tps_pnt *pnts, int n_points, int min_points,
- double regularization, double overlap, double pthin,
- int clustered, double segs_mb)
+int local_tps(SEGMENT *in_seg, SEGMENT *var_seg, int n_vars,
+ SEGMENT *out_seg, int out_fd, char *mask_name,
+ struct Cell_head *src, struct Cell_head *dst,
+ off_t n_points, int min_points,
+ double regularization, double overlap, int clustered)
{
int ridx, cidx, row, col, nrows, ncols;
double **m, *a, *B;
@@ -356,31 +404,23 @@
int kdalloc, palloc, n_cur_points;
struct tps_pnt *cur_pnts;
double dx, dy, dist, dist2, mfactor;
- DCELL **dbuf, result, *outbuf, *varbuf;
+ DCELL *dval, result, *outbuf, *varbuf;
CELL *maskbuf;
int solved;
- struct kdtree *kdt;
- double c[2];
- int *kduid, kdfound;
- double *kddist;
- int *bfsid, bfsfound, pfound;
+ int kdfound, bfsfound, pfound;
+ double distmax;
- FLAG *mask_flag;
- SEGMENT var_seg, out_seg, p_seg;
- int nsegs;
- double segsize;
- int varsize;
- struct tps_out {
- double val;
- double wsum;
- double wmax;
- } tps_out;
+ int mask_fd;
+ FLAG *mask_flag, *pnt_flag;
+ struct tps_out tps_out;
+ double *pvar;
double weight, dxi, dyi;
double wmin, wmax;
int rmin, rmax, cmin, cmax, rminp, rmaxp, cminp, cmaxp;
int irow, irow1, irow2, icol, icol1, icol2;
- int nskipped;
- double pthin2;
+#ifndef USE_RC
+ double i_n, i_e;
+#endif
nrows = Rast_window_rows();
ncols = Rast_window_cols();
@@ -392,16 +432,26 @@
for (i = 0; i < (palloc + 1 + n_vars); i++)
m[i] = G_malloc((palloc + 1 + n_vars) * sizeof(double));
cur_pnts = G_malloc(palloc * 5 * sizeof(struct tps_pnt));
- kduid = G_malloc(palloc * sizeof(int));
- kddist = G_malloc(palloc * sizeof(double));
+
+ pvar = NULL;
+ varbuf = NULL;
+ if (n_vars) {
+ pvar = G_malloc(palloc * 5 * n_vars * sizeof(double));
+ cur_pnts[0].vars = pvar;
+ for (i = 1; i < palloc * 5; i++)
+ cur_pnts[i].vars = cur_pnts[i - 1].vars + n_vars;
- bfsid = G_malloc(min_points * sizeof(int));
+ varbuf = G_malloc(n_vars * sizeof(DCELL));
+ }
+ dval = G_malloc((1 + n_vars) * sizeof(DCELL));
+
mask_flag = flag_create(nrows, ncols);
maskbuf = NULL;
- if (mask_fd >= 0) {
+ if (mask_name) {
G_message("Loading mask map...");
+ mask_fd = Rast_open_old(mask_name, "");
maskbuf = Rast_allocate_c_buf();
for (row = 0; row < nrows; row++) {
@@ -411,60 +461,41 @@
FLAG_SET(mask_flag, row, col);
}
}
+ Rast_close(mask_fd);
G_free(maskbuf);
maskbuf = NULL;
}
- /* create and init SEG structs */
- varsize = n_vars * sizeof(DCELL);
- segsize = (double)64 * 64 * (sizeof(struct tps_out) + varsize + clustered * sizeof(int));
- nsegs = 1024. * 1024. * segs_mb / segsize;
-
- if (Segment_open(&out_seg, G_tempfile(), nrows, ncols, 64, 64,
- sizeof(struct tps_out), nsegs) != 1) {
- G_fatal_error("Unable to create input temporary files");
- }
+ G_message(_("Analyzing input ..."));
+ pnt_flag = flag_create(src->rows, src->cols);
- dbuf = NULL;
- varbuf = NULL;
- if (n_vars) {
- dbuf = G_malloc(n_vars * sizeof(DCELL *));
- varbuf = G_malloc(n_vars * sizeof(DCELL));
- for (i = 0; i < n_vars; i++) {
- dbuf[i] = Rast_allocate_d_buf();
- }
+ rminp = src->rows;
+ rmaxp = 0;
+ cminp = src->cols;
+ cmaxp = 0;
+ for (row = 0; row < src->rows; row++) {
+ G_percent(row, src->rows, 2);
- if (Segment_open(&var_seg, G_tempfile(), nrows, ncols, 64, 64,
- varsize, nsegs) != 1) {
- G_fatal_error("Unable to create input temporary files");
- }
- }
-
- if (clustered) {
- G_message("Creating temporary point map...");
-
- if (Segment_open(&p_seg, G_tempfile(), nrows, ncols, 64, 64,
- sizeof(int), nsegs) != 1) {
- G_fatal_error("Unable to create input temporary files");
- }
- i = -1;
- for (row = 0; row < nrows; row++) {
- G_percent(row, nrows, 2);
-
- for (col = 0; col < ncols; col++) {
- if (Segment_put(&p_seg, (void *)&i, row, col) != 1)
- G_fatal_error(_("Unable to write to temporary file"));
+ for (col = 0; col < src->cols; col++) {
+ Segment_get(in_seg, (void *)dval, row, col);
+ if (!Rast_is_d_null_value(dval)) {
+ FLAG_SET(pnt_flag, row, col);
+ if (rminp > row)
+ rminp = row;
+ if (rmaxp < row)
+ rmaxp = row;
+ if (cminp > col)
+ cminp = col;
+ if (cmaxp < col)
+ cmaxp = col;
}
}
- for (i = 0; i < n_points; i++) {
- G_percent(i, n_points, 2);
- if (pnts[i].r < 0 || pnts[i].r >= nrows ||
- pnts[i].c < 0 || pnts[i].c >= ncols)
- continue;
- if (Segment_put(&p_seg, (void *)&i, pnts[i].r, pnts[i].c) != 1)
- G_fatal_error(_("Unable to write to temporary file"));
- }
}
+ rminp = row_src2dst(rminp, src, dst);
+ rmaxp = row_src2dst(rmaxp, src, dst);
+ cminp = col_src2dst(cminp, src, dst);
+ cmaxp = col_src2dst(cmaxp, src, dst);
+ G_percent(1, 1, 2);
G_message(_("Initializing output ..."));
tps_out.val = 0;
@@ -472,77 +503,16 @@
tps_out.wmax = 0;
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 2);
- if (n_vars) {
- for (j = 0; j < n_vars; j++)
- Rast_get_d_row(var_fd[j], dbuf[j], row);
- }
for (col = 0; col < ncols; col++) {
- if (n_vars) {
- for (j = 0; j < n_vars; j++) {
- varbuf[j] = dbuf[j][col];
- }
- if (Segment_put(&var_seg, (void *)varbuf, row, col) != 1)
- G_fatal_error(_("Unable to write to temporary file"));
- }
- if (Segment_put(&out_seg, (void *)&tps_out, row, col) != 1)
+ if (Segment_put(out_seg, (void *)&tps_out, row, col) != 1)
G_fatal_error(_("Unable to write to temporary file"));
}
}
- if (n_vars) {
- for (j = 0; j < n_vars; j++)
- G_free(dbuf[j]);
- G_free(dbuf);
- dbuf = NULL;
- }
G_percent(1, 1, 2);
- /* create k-d tree */
- G_message(_("Creating search index ..."));
- pthin2 = pthin * pthin;
- if (pthin > 0)
- G_message(_("Thinning points with minimum distance %g"), pthin);
- kdt = kdtree_create(2, NULL);
- nskipped = 0;
- rminp = nrows;
- rmaxp = 0;
- cminp = ncols;
- cmaxp = 0;
- for (i = 0; i < n_points; i++) {
- G_percent(i, n_points, 2);
+ G_message(_("Local TPS interpolation with %ld points..."), n_points);
- c[0] = pnts[i].c;
- c[1] = pnts[i].r;
-
- if (pthin > 0 && i > 0) {
- kdfound = kdtree_knn(kdt, c, kduid, kddist, 1, NULL);
- if (kdfound && kddist[0] <= pthin2) {
- nskipped++;
- continue;
- }
- }
-
- kdtree_insert(kdt, c, i, 0);
-
- if (rminp > pnts[i].r)
- rminp = pnts[i].r;
- if (rmaxp < pnts[i].r)
- rmaxp = pnts[i].r;
- if (cminp > pnts[i].c)
- cminp = pnts[i].c;
- if (cmaxp < pnts[i].c)
- cmaxp = pnts[i].c;
- }
- kdtree_optimize(kdt, 2);
- G_percent(1, 1, 2);
-
- if (nskipped) {
- G_message(_("%d of %d points were thinned out"), nskipped, n_points);
- n_points -= nskipped;
- }
-
- G_message(_("Local TPS interpolation with %d points..."), n_points);
-
wmin = 10;
wmax = 0;
@@ -572,28 +542,20 @@
}
if (n_vars) {
- int isnull = 0;
- Segment_get(&var_seg, (void *)varbuf, row, col);
- for (j = 0; j < n_vars; j++) {
- if (Rast_is_d_null_value(&varbuf[j])) {
- isnull = 1;
- break;
- }
- }
- if (isnull)
+ Segment_get(var_seg, (void *)varbuf, row, col);
+
+ if (Rast_is_d_null_value(varbuf))
continue;
}
- Segment_get(&out_seg, (void *)&tps_out, row, col);
+ Segment_get(out_seg, (void *)&tps_out, row, col);
if (tps_out.wmax > overlap)
continue;
solved = 0;
n_cur_points = 0;
kdfound = 0;
- c[0] = col;
- c[1] = row;
while (!solved) {
@@ -604,38 +566,34 @@
/* alloc */
if (kdalloc < n_cur_points) {
- G_free(kduid);
- G_free(kddist);
- G_free(bfsid);
G_free(cur_pnts);
kdalloc = n_cur_points;
- kduid = G_malloc(kdalloc * sizeof(int));
- kddist = G_malloc(kdalloc * sizeof(double));
- bfsid = G_malloc(kdalloc * sizeof(int));
cur_pnts = G_malloc(kdalloc * 5 * sizeof(struct tps_pnt));
+ if (n_vars) {
+ G_free(pvar);
+ pvar = G_malloc(kdalloc * 5 * n_vars * sizeof(double));
+ cur_pnts[0].vars = pvar;
+ for (i = 1; i < kdalloc * 5; i++)
+ cur_pnts[i].vars = cur_pnts[i - 1].vars + n_vars;
+ }
}
/* collect nearest neighbors */
- kdfound = kdtree_knn(kdt, c, kduid, kddist, n_cur_points, NULL);
-
- /* load points to matrix */
- rmin = nrows;
+ rmin = src->rows;
rmax = 0;
- cmin = ncols;
+ cmin = src->cols;
cmax = 0;
- for (i = 0; i < kdfound; i++) {
- cur_pnts[i] = pnts[kduid[i]];
- if (rmin > cur_pnts[i].r)
- rmin = cur_pnts[i].r;
- if (rmax < cur_pnts[i].r)
- rmax = cur_pnts[i].r;
- if (cmin > cur_pnts[i].c)
- cmin = cur_pnts[i].c;
- if (cmax < cur_pnts[i].c)
- cmax = cur_pnts[i].c;
- }
+ kdfound = bfs_search_nn(pnt_flag, src, cur_pnts,
+ n_cur_points,
+ row_dst2src(row, src, dst),
+ col_dst2src(col, src, dst),
+ &rmin, &rmax, &cmin, &cmax,
+ &distmax);
+
+ /* convert src min/max to dst min/max */
+
pfound = kdfound;
if (clustered) {
@@ -644,65 +602,65 @@
/* qrt1: 0, row, col + 1, ncols - 1 */
if (rminp <= row && cmaxp > col) {
- bfsfound = bfs_search(&p_seg, mask_flag, bfsid,
- row, col, n_cur_points / 2,
- kddist[kdfound - 1],
- 0, row, col + 1, ncols - 1);
+ bfsfound = bfs_search(pnt_flag, src,
+ cur_pnts + pfound,
+ row_dst2src(row, src, dst),
+ col_dst2src(col, src, dst),
+ n_cur_points / 3, distmax,
+ 0, row_dst2src(row, src, dst),
+ col_dst2src(col, src, dst) + 1, src->cols - 1);
if (bfsfound == 0)
G_debug(4, "No BFS points for NE quadrant");
- for (i = 0; i < bfsfound; i++) {
- cur_pnts[pfound + i] = pnts[bfsid[i]];
- }
pfound += bfsfound;
}
/* qrt2: 0, row - 1, 0, col */
if (rminp < row && cminp <= col) {
- bfsfound = bfs_search(&p_seg, mask_flag, bfsid,
- row, col, n_cur_points / 2,
- kddist[kdfound - 1],
- 0, row - 1, 0, col);
+ bfsfound = bfs_search(pnt_flag, src,
+ cur_pnts + pfound,
+ row_dst2src(row, src, dst),
+ col_dst2src(col, src, dst),
+ n_cur_points / 3, distmax,
+ 0, row_dst2src(row, src, dst) - 1,
+ 0, col_dst2src(col, src, dst));
if (bfsfound == 0)
G_debug(4, "No BFS points for NW quadrant");
- for (i = 0; i < bfsfound; i++) {
- cur_pnts[pfound + i] = pnts[bfsid[i]];
- }
pfound += bfsfound;
}
/* qrt3: row, nrows - 1, 0, col - 1 */
if (rmaxp >= row && cminp < col) {
- bfsfound = bfs_search(&p_seg, mask_flag, bfsid,
- row, col, n_cur_points / 2,
- kddist[kdfound - 1],
- row, nrows - 1, 0, col - 1);
+ bfsfound = bfs_search(pnt_flag, src,
+ cur_pnts + pfound,
+ row_dst2src(row, src, dst),
+ col_dst2src(col, src, dst),
+ n_cur_points / 3, distmax,
+ row_dst2src(row, src, dst), src->rows - 1,
+ 0, col_dst2src(col, src, dst) - 1);
if (bfsfound == 0)
G_debug(4, "No BFS points for SW quadrant");
- for (i = 0; i < bfsfound; i++) {
- cur_pnts[pfound + i] = pnts[bfsid[i]];
- }
pfound += bfsfound;
}
/* qrt4: row + 1, nrows - 1, col, ncols - 1 */
if (rmaxp > row && cmaxp >= col) {
- bfsfound = bfs_search(&p_seg, mask_flag, bfsid,
- row, col, n_cur_points / 2,
- kddist[kdfound - 1],
- row + 1, nrows - 1, col, ncols - 1);
+ bfsfound = bfs_search(pnt_flag, src,
+ cur_pnts + pfound,
+ row_dst2src(row, src, dst),
+ col_dst2src(col, src, dst),
+ n_cur_points / 3, distmax,
+ row_dst2src(row, src, dst) + 1, src->rows - 1,
+ col_dst2src(col, src, dst), src->cols - 1);
if (bfsfound == 0)
G_debug(4, "No BFS points for SE quadrant");
- for (i = 0; i < bfsfound; i++) {
- cur_pnts[pfound + i] = pnts[bfsid[i]];
- }
pfound += bfsfound;
}
}
@@ -722,12 +680,14 @@
m[i] = G_malloc((palloc + 1 + n_vars) * sizeof(double));
}
-
/* sort points */
/* qsort(cur_pnts, kdfound, sizeof(struct tps_pnt), cmp_pnts); */
- load_tps_pnts(cur_pnts, pfound, n_vars, regularization, m, a);
+ qsort(cur_pnts, pfound, sizeof(struct tps_pnt), cmp_pnts);
+ load_tps_pnts(in_seg, n_vars, cur_pnts, pfound, src, dst,
+ regularization, m, a);
+
/* solve */
solved = solvemat(m, a, B, pfound + 1 + n_vars);
@@ -735,21 +695,6 @@
G_fatal_error(_("Matrix is unsolvable"));
}
- rmin = nrows;
- rmax = 0;
- cmin = ncols;
- cmax = 0;
- for (i = 0; i < kdfound; i++) {
- if (rmin > cur_pnts[i].r)
- rmin = cur_pnts[i].r;
- if (rmax < cur_pnts[i].r)
- rmax = cur_pnts[i].r;
- if (cmin > cur_pnts[i].c)
- cmin = cur_pnts[i].c;
- if (cmax < cur_pnts[i].c)
- cmax = cur_pnts[i].c;
- }
-
/* must be <= 0.5 */
mfactor = 0.0;
/* min: 0
@@ -770,12 +715,17 @@
mfactoradj = 0.0;
if (dmax > dmin) {
- mfactoradj = pow((1 - dmin / dmax), 2.0) * 0.5;
+ mfactoradj = pow(1 - dmin / dmax, 2.0) * 0.5;
G_debug(1, "adjusted mfactor: %g", mfactoradj);
}
mfactor = mfactoradj;
}
-
+
+ rmin = row_src2dst(rmin, src, dst);
+ rmax = row_src2dst(rmax, src, dst);
+ cmin = col_src2dst(cmin, src, dst);
+ cmax = col_src2dst(cmax, src, dst);
+
irow1 = rmin + (int)((rmax - rmin) * mfactor);
irow2 = rmax - (int)((rmax - rmin) * mfactor);
icol1 = cmin + (int)((cmax - cmin) * mfactor);
@@ -820,38 +770,46 @@
dyi = irow2 - irow1 + 1;
for (irow = irow1; irow <= irow2; irow++) {
+
+#ifndef USE_RC
+ i_n = dst->north - (irow + 0.5) * dst->ns_res;
+#endif
for (icol = icol1; icol <= icol2; icol++) {
if ((FLAG_GET(mask_flag, irow, icol))) {
continue;
}
if (n_vars) {
- int isnull = 0;
- Segment_get(&var_seg, (void *)varbuf, irow, icol);
- for (j = 0; j < n_vars; j++) {
- if (Rast_is_d_null_value(&varbuf[j])) {
- isnull = 1;
- break;
- }
- }
- if (isnull)
+ Segment_get(var_seg, (void *)varbuf, irow, icol);
+ if (Rast_is_d_null_value(varbuf)) {
continue;
+ }
}
- Segment_get(&out_seg, (void *)&tps_out, irow, icol);
+#ifndef USE_RC
+ i_e = dst->west + (icol + 0.5) * dst->ew_res;
+#endif
+ Segment_get(out_seg, (void *)&tps_out, irow, icol);
+ j = 0;
+
result = B[0];
if (n_vars) {
- Segment_get(&var_seg, (void *)varbuf, irow, icol);
- for (j = 0; j < n_vars; j++)
+ for (j = 0; j < n_vars; j++) {
result += varbuf[j] * B[j + 1];
+ }
}
for (i = 0; i < pfound; i++) {
+#ifdef USE_RC
dx = (cur_pnts[i].c - icol) * 2.0;
dy = (cur_pnts[i].r - irow) * 2.0;
-
+#else
+ dx = cur_pnts[i].c - i_e;
+ dy = cur_pnts[i].r - i_n;
+#endif
+
dist2 = dx * dx + dy * dy;
dist = 0;
if (dist2 > 0) {
@@ -871,28 +829,24 @@
wmin = weight;
if (wmax < weight)
wmax = weight;
-
+
if (tps_out.wmax < weight)
tps_out.wmax = weight;
tps_out.val += result * weight;
tps_out.wsum += weight;
- Segment_put(&out_seg, (void *)&tps_out, irow, icol);
+ Segment_put(out_seg, (void *)&tps_out, irow, icol);
}
}
}
}
G_percent(1, 1, 1);
- G_debug(1, "wmin: %g", wmin);
- G_debug(1, "wmax: %g", wmax);
+ G_debug(0, "wmin: %g", wmin);
+ G_debug(0, "wmax: %g", wmax);
- if (n_vars)
- Segment_close(&var_seg);
+ flag_destroy(pnt_flag);
- if (clustered)
- Segment_close(&p_seg);
-
outbuf = Rast_allocate_d_buf();
G_message(_("Writing output..."));
@@ -906,7 +860,7 @@
continue;
}
- Segment_get(&out_seg, (void *)&tps_out, row, col);
+ Segment_get(out_seg, (void *)&tps_out, row, col);
if (tps_out.wsum == 0)
Rast_set_d_null_value(&outbuf[col], 1);
@@ -918,7 +872,7 @@
G_percent(1, 1, 2);
G_free(outbuf);
- Segment_close(&out_seg);
+ flag_destroy(mask_flag);
return 1;
}
Modified: grass-addons/grass7/raster/r.resamp.tps/tps.h
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/tps.h 2016-11-21 20:46:28 UTC (rev 69862)
+++ grass-addons/grass7/raster/r.resamp.tps/tps.h 2016-11-21 20:49:11 UTC (rev 69863)
@@ -1,14 +1,19 @@
struct tps_pnt
{
- int r, c; /* row, col in target window */
+ double r, c; /* row, col in source window */
double val; /* value to interpolate */
double *vars; /* values of covariables */
};
-int global_tps(int out_fd, int *var_fd, int n_vars, int mask_fd,
- struct tps_pnt *pnts, int n_points, double regularization);
-int local_tps(int out_fd, int *var_fd, int n_vars, int mask_fd,
- struct tps_pnt *pnts, int n_points, int min_points,
- double regularization, double overlap, double pthin,
- int do_bfs, double segs_mb);
+struct tps_out {
+ double val;
+ double wsum;
+ double wmax;
+};
+
+int local_tps(SEGMENT *in_seg, SEGMENT *var_seg, int n_vars,
+ SEGMENT *out_seg, int out_fd, char *mask_name,
+ struct Cell_head *src, struct Cell_head *dst,
+ off_t n_points, int min_points,
+ double regularization, double overlap, int do_bfs);
Deleted: grass-addons/grass7/raster/r.resamp.tps/v.surf.tps.html
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/v.surf.tps.html 2016-11-21 20:46:28 UTC (rev 69862)
+++ grass-addons/grass7/raster/r.resamp.tps/v.surf.tps.html 2016-11-21 20:49:11 UTC (rev 69863)
@@ -1,141 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>v.surf.tps</em> performs multivariate thin plate spline
-interpolation with regularization. The <b>input</b> is a 2D
-or 3D vector <em>points</em> map. Values to interpolate can be the z
-values of 3D points or the values in a user-specified attribute column
-in a 2D or 3D vector map. Output is a raster map. Optionally, several
-raster maps can be specified to be used as covariables which will
-improve results in areas with few points. The module can be regarded
-as a combination of a multiple regression and spline interpolation.
-
-<p>
-The <b>min</b> options specifies the minimum number of points to be
-used for interpolation. If the number of input points is smaller than
-or equal to the minimum number of points, global TPS interpolation is
-used. If the number of input points is larger than the minimum number
-of points, tiled local TPS interpolation is used. Tile sizes are
-variable and dependent on the extents of the <b>min</b> nearest
-neighbors when a new tile is generated.
-
-<p>
-The <b>smooth</b> option can be used to reduce the influence of the
-splines and increase the influence of the covariables. Without
-covariables, the resulting surface will be smoother. With covariables
-and a large smooting value, the resulting surface will be mainly
-determined by the multiple regression component.
-
-<p>
-The <b>overlap</b> option controls how much tiles are overlapping when
-the <b>min</b> option is smaller than the numer of input points.
-Tiling artefacts occur with low values for the <b>min</b> option and the
-<b>overlap</b> option. Increasing both options will reduce tiling
-artefacts but processing will take more time. Values for the
-<b>overlap</b> option must be between 0 and 1.
-
-<p>
-The module works best with evenly spaced sparse points. In case of
-highly unevenly spaced points, e.g. remote sensing data with gaps due
-to cloud cover, the <b>thin</b> option should be used in order to avoid
-tiling artefacts, otherwise a high number of minimum points and a large
-<b>overlap</b> value are required, slowing down the module.
-
-<p>
-The <b>memory</b> option controls only how much memory should be used
-for the covariables and the intermediate output. The input points are
-always completely loaded to memory.
-
-<h2>EXAMPLES</h2>
-
-The computational region setting for the following examples:
-
-<div class="code"><pre>
-g.region -p rast=elev_state_500m
-</pre></div>
-
-<h3>Basic interpolation</h3>
-Interpolation of 30 year precipitation normals in the North Carlolina
-sample dataset:
-
-<div class="code"><pre>
-v.surf.tps input=precip_30ynormals_3d output=precip_30ynormals_3d \
- column=annual min=140
-</pre></div>
-
-<h3>Interpolation with a covariable</h3>
-
-<div class="code"><pre>
-v.surf.tps input=precip_30ynormals_3d output=precip_30ynormals_3d \
- column=annual min=140 covars=elev_state_500m
-</pre></div>
-
-<h3>Interpolation with a covariable and smoothing</h3>
-
-<div class="code"><pre>
-v.surf.tps input=precip_30ynormals_3d output=precip_30ynormals_3d \
- column=annual min=140 covars=elev_state_500m smooth=0.1
-</pre></div>
-
-<h3>Tiled interpolation with a covariable and smoothing</h3>
-
-<div class="code"><pre>
-v.surf.tps input=precip_30ynormals_3d output=precip_30ynormals_3d \
- column=annual min=20 covars=elev_state_500m smooth=0.1 \
- overlap=0.1
-</pre></div>
-
-<!--
-r.colors map=precip_30ynormals_3d rules=- <<EOF
-950 red
-1000 orange
-1200 yellow
-1400 cyan
-1600 aqua
-1800 blue
-2500 violet
-EOF
-m.nviz.image elevation_map=elev_state_500m -a \
- mode=fine resolution_fine=1 color_map=precip_30ynormals_3d \
- position=0.59,0.74 height=40959 perspective=20 twist=0 \
- zexag=10 focus=105434,77092,1073 \
- light_position=0.68,-0.68,0.80 light_brightness=80 \
- output=v_surf_tps.tif format=tif size=798,585
-d.legend -b raster=precip_30ynormals_3d label_step=300 \
- title="Annual Precipitation in Western North Carolina"
-mogrify -trim -resize 600x -format png v_surf_tps.tif
-optipng -o5 v_surf_tps.png
--->
-
-<center>
- <img src="v_surf_tps.png">
- <p><em>
- Precipitation computed based on annual normals and
- elevation as a covariable
- </em></p>
-</center>
-
-<h2>REFERENCES</h2>
-
-<ul>
- <li>Hutchinson MF, 1995, Interpolating mean rainfall using thin plate
- smoothing splines. International Journal of Geographical Information
- Systems, 9(4), pp. 385-403</li>
- <li>Wahba G, 1990, Spline models for observational data. In CBMS-NSF
- Regional Conference Series in Applied Mathematics. Philadelpia:
- Society for Industrial and Applied Mathematics</li>
-</ul>
-
-<h2>SEE ALSO</h2>
-
-<em>
-<a href="v.surf.rst.html">v.surf.rst</a>,
-<a href="v.surf.bspline.html">v.surf.rst</a>,
-<a href="v.surf.idw.html">v.surf.idw</a>
-</em>
-
-<h2>AUTHORS</h2>
-
-Markus Metz
-
-<p>
-<i>Last changed: $Date$</i>
Deleted: grass-addons/grass7/raster/r.resamp.tps/v_surf_tps.png
===================================================================
(Binary files differ)
More information about the grass-commit
mailing list