[GRASS-SVN] r69257 - grass-addons/grass7/vector/v.surf.tps
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Aug 25 11:50:34 PDT 2016
Author: mmetz
Date: 2016-08-25 11:50:34 -0700 (Thu, 25 Aug 2016)
New Revision: 69257
Added:
grass-addons/grass7/vector/v.surf.tps/Makefile
grass-addons/grass7/vector/v.surf.tps/flag.c
grass-addons/grass7/vector/v.surf.tps/flag.h
grass-addons/grass7/vector/v.surf.tps/main.c
grass-addons/grass7/vector/v.surf.tps/rclist.c
grass-addons/grass7/vector/v.surf.tps/rclist.h
grass-addons/grass7/vector/v.surf.tps/tps.c
grass-addons/grass7/vector/v.surf.tps/tps.h
grass-addons/grass7/vector/v.surf.tps/v.surf.tps.html
Modified:
grass-addons/grass7/vector/v.surf.tps/
Log:
new addon v.surf.tps
Property changes on: grass-addons/grass7/vector/v.surf.tps
___________________________________________________________________
Added: svn:ignore
+ OBJ.* *.tmp.html
Added: grass-addons/grass7/vector/v.surf.tps/Makefile
===================================================================
--- grass-addons/grass7/vector/v.surf.tps/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.surf.tps/Makefile 2016-08-25 18:50:34 UTC (rev 69257)
@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.surf.tps
+
+LIBES = $(VECTORLIB) $(DBMILIB) $(RASTERLIB) $(BTREE2LIB) $(SEGMENTLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(VECTORDEP) $(DBMIDEP) $(RASTERDEP) $(BTREE2DEP) $(SEGMENTDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Property changes on: grass-addons/grass7/vector/v.surf.tps/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.surf.tps/flag.c
===================================================================
--- grass-addons/grass7/vector/v.surf.tps/flag.c (rev 0)
+++ grass-addons/grass7/vector/v.surf.tps/flag.c 2016-08-25 18:50:34 UTC (rev 69257)
@@ -0,0 +1,60 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "flag.h"
+
+FLAG *flag_create(int nrows, int ncols)
+{
+ unsigned char *temp;
+
+ FLAG *new_flag;
+
+ register int i;
+
+ new_flag = (FLAG *) G_malloc(sizeof(FLAG));
+ new_flag->nrows = nrows;
+ new_flag->ncols = ncols;
+ new_flag->leng = (ncols + 7) / 8;
+ new_flag->array =
+ (unsigned char **)G_malloc(nrows * sizeof(unsigned char *));
+
+ if (!new_flag->array)
+ G_fatal_error(_("Out of memory!"));
+
+ temp =
+ (unsigned char *)G_malloc(nrows * new_flag->leng *
+ sizeof(unsigned char));
+
+ if (!temp)
+ G_fatal_error(_("Out of memory!"));
+
+ for (i = 0; i < nrows; i++) {
+ new_flag->array[i] = temp;
+ temp += new_flag->leng;
+ }
+
+ flag_clear_all(new_flag);
+
+ return (new_flag);
+}
+
+int flag_destroy(FLAG * flags)
+{
+ G_free(flags->array[0]);
+ G_free(flags->array);
+ G_free(flags);
+
+ return 0;
+}
+
+int flag_clear_all(FLAG * flags)
+{
+ register int r, c;
+
+ for (r = 0; r < flags->nrows; r++) {
+ for (c = 0; c < flags->leng; c++) {
+ flags->array[r][c] = 0;
+ }
+ }
+
+ return 0;
+}
Property changes on: grass-addons/grass7/vector/v.surf.tps/flag.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.surf.tps/flag.h
===================================================================
--- grass-addons/grass7/vector/v.surf.tps/flag.h (rev 0)
+++ grass-addons/grass7/vector/v.surf.tps/flag.h 2016-08-25 18:50:34 UTC (rev 69257)
@@ -0,0 +1,64 @@
+#ifndef __FLAG_H__
+#define __FLAG_H__
+
+
+/* flag.[ch] is a set of routines which will set up an array of bits
+ ** that allow the programmer to "flag" cells in a raster map.
+ **
+ ** FLAG *
+ ** flag_create(nrows,ncols)
+ ** int nrows, ncols;
+ ** opens the structure flag.
+ ** The flag structure will be a two dimensional array of bits the
+ ** size of nrows by ncols. Will initalize flags to zero (unset).
+ **
+ ** flag_destroy(flags)
+ ** FLAG *flags;
+ ** closes flags and gives the memory back to the system.
+ **
+ ** flag_clear_all(flags)
+ ** FLAG *flags;
+ ** sets all values in flags to zero.
+ **
+ ** flag_unset(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ ** sets the value of (row, col) in flags to zero.
+ **
+ ** flag_set(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ ** will set the value of (row, col) in flags to one.
+ **
+ ** int
+ ** flag_get(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ ** returns the value in flags that is at (row, col).
+ **
+ ** idea by Michael Shapiro
+ ** code by Chuck Ehlschlaeger
+ ** April 03, 1989
+ */
+#define FLAG struct _flagsss_
+FLAG {
+ int nrows, ncols, leng;
+ unsigned char **array;
+};
+
+#define FLAG_UNSET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] &= ~(1<<((col) & 7))
+
+#define FLAG_SET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] |= (1<<((col) & 7))
+
+#define FLAG_GET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] & (1<<((col) & 7))
+
+/* flag.c */
+FLAG *flag_create(int, int);
+int flag_destroy(FLAG *);
+int flag_clear_all(FLAG *);
+
+
+#endif /* __FLAG_H__ */
Property changes on: grass-addons/grass7/vector/v.surf.tps/flag.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.surf.tps/main.c
===================================================================
--- grass-addons/grass7/vector/v.surf.tps/main.c (rev 0)
+++ grass-addons/grass7/vector/v.surf.tps/main.c 2016-08-25 18:50:34 UTC (rev 69257)
@@ -0,0 +1,512 @@
+
+/**********************************************************************
+ *
+ * MODULE: v.surf.tps
+ *
+ * AUTHOR(S): Markus Metz
+ *
+ * PURPOSE: Thin Plate Spline interpolation with covariables
+ *
+ * COPYRIGHT: (C) 2016 by by 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 <math.h>
+#include <grass/raster.h>
+#include <grass/vector.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;
+
+ struct Map_info In;
+ struct History history;
+
+ 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 Flag *c_flag;
+ struct Cell_head dstwindow;
+
+ 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;
+
+ 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;
+ double regularization, overlap;
+ double pthin;
+ int segs_mb;
+
+ /*----------------------------------------------------------------*/
+ /* Options declarations */
+ module = G_define_module();
+ G_add_keyword(_("vector"));
+ 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");
+
+ 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;
+ reg_opt->required = NO;
+ reg_opt->answer = "0";
+ reg_opt->description =
+ _("Smoothing factor");
+ reg_opt->guisection = _("Settings");
+
+ ov_opt = G_define_option();
+ ov_opt->key = "overlap";
+ ov_opt->type = TYPE_DOUBLE;
+ ov_opt->required = NO;
+ ov_opt->answer = "0.8";
+ ov_opt->label =
+ _("Overlap factor <= 1");
+ ov_opt->description =
+ _("A larger value increase the tile overlap");
+ ov_opt->guisection = _("Settings");
+
+ minpnts_opt = G_define_option();
+ minpnts_opt->key = "min";
+ minpnts_opt->type = TYPE_DOUBLE;
+ minpnts_opt->required = NO;
+ minpnts_opt->answer = "20";
+ minpnts_opt->description =
+ _("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");
+
+ 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");
+
+ out_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+ out_opt->key = "output";
+ out_opt->required = YES;
+
+ mask_opt = G_define_standard_option(G_OPT_R_INPUT);
+ mask_opt->key = "mask";
+ mask_opt->label = _("Raster map to use for masking");
+ mask_opt->description = _("Only cells that are not NULL and not zero are interpolated");
+ mask_opt->required = NO;
+
+ mem_opt = G_define_option();
+ mem_opt->key = "memory";
+ mem_opt->type = TYPE_INTEGER;
+ mem_opt->required = NO;
+ mem_opt->answer = "300";
+ mem_opt->description = _("Memory in MB");
+
+ c_flag = G_define_flag();
+ c_flag->key = 'c';
+ c_flag->description = _("Input points are dense clusters separated by empty areas");
+
+
+ /* Parsing */
+ G_gisinit(argv[0]);
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ 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);
+
+ 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);
+
+ 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"));
+ }
+ else { /* 2D */
+ if (!val_column)
+ G_fatal_error(_("Input vector map is 2D. Parameter <%s> required."), col_opt->key);
+ }
+
+ if (!use_z) {
+ val_field = Vect_get_field_number(&In, dfield_opt->answer);
+ }
+
+ n_vars = 0;
+ if (var_opt->answer) {
+ while (var_opt->answers[n_vars])
+ n_vars++;
+ }
+
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ Rast_get_window(&dstwindow);
+ nrows = dstwindow.rows;
+ ncols = dstwindow.cols;
+
+ /* count points */
+ Vect_rewind(&In);
+ n_points = 0;
+ while (1) {
+ /* register line */
+ type = Vect_read_next_line(&In, Points, Cats);
+
+ /* 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) {
+ n_points++;
+
+ if (n_vars) {
+ r = (dstwindow.north - Points->y[0]) / dstwindow.ns_res;
+ c = (Points->x[0] - dstwindow.west) / dstwindow.ew_res;
+
+ 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);
+
+ /* check memory */
+ pnts = G_malloc(n_points * sizeof(struct tps_pnt));
+ if (!pnts)
+ G_fatal_error(_("Too many points, out of memory"));
+
+ 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);
+ }
+
+ /* 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"));
+
+ driver = db_start_driver_open_database(Fi->driver, Fi->database);
+
+ 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);
+ }
+ }
+
+ G_message(_("Loading %d points..."), n_points);
+ Vect_rewind(&In);
+ i = 0;
+ j = 0;
+ while (1) {
+
+ i++;
+
+ /* register line */
+ type = Vect_read_next_line(&In, Points, Cats);
+
+ /* 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;
+ }
+ case DB_C_TYPE_DOUBLE: {
+ val = db_get_value_double(&dbval);
+ break;
+ }
+ default: {
+ Rast_set_d_null_value(&val, 1);
+ break;
+ }
+ }
+ }
+ 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);
+
+ if (driver)
+ db_close_database_shutdown_driver(driver);
+ Vect_close(&In);
+
+ if (j != n_points)
+ n_points = j;
+
+ /* sort points */
+ qsort(pnts, n_points, sizeof(struct tps_pnt), cmp_pnts);
+
+ /* 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 (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;
+
+ G_message(_("Loading %d covariables..."), n_vars);
+
+ cur_row = -1;
+
+ 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], "");
+ 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);
+ }
+ 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++;
+ }
+ }
+ }
+ G_percent(1, 1, 1);
+ n_points = ii;
+ }
+ if (n_vars) {
+ for (i = 0; i < n_vars; i++) {
+ G_free(dbuf[i]);
+ }
+ G_free(dbuf);
+ }
+
+ out_fd = Rast_open_new(outname, DCELL_TYPE);
+
+ min_points = atoi(minpnts_opt->answer);
+ if (min_points < 3 + n_vars) {
+ min_points = 3 + n_vars;
+ G_warning(_("Minimum number of points is too small, set to %d"),
+ min_points);
+ }
+
+ regularization = atof(reg_opt->answer);
+ if (regularization < 0)
+ regularization = 0;
+
+ overlap = atof(ov_opt->answer);
+ if (overlap < 0)
+ overlap = 0;
+ 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 (map_tps(out_fd, var_fd, n_vars, mask_fd, pnts, n_points,
+ regularization) != 1) {
+ G_fatal_error(_("TPS interpolation failed"));
+ }
+ }
+ else {
+ if (cell_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);
+
+ /* write map history */
+ Rast_close(out_fd);
+ Rast_short_history(outname, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(outname, &history);
+
+ G_done_msg(" ");
+
+ exit(EXIT_SUCCESS);
+} /*END MAIN */
Property changes on: grass-addons/grass7/vector/v.surf.tps/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.surf.tps/rclist.c
===================================================================
--- grass-addons/grass7/vector/v.surf.tps/rclist.c (rev 0)
+++ grass-addons/grass7/vector/v.surf.tps/rclist.c 2016-08-25 18:50:34 UTC (rev 69257)
@@ -0,0 +1,68 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "rclist.h"
+
+void rclist_init(struct rclist *list)
+{
+ list->head = list->tail = NULL;
+
+ return;
+}
+
+void rclist_add(struct rclist *list, int row, int col)
+{
+ struct rc *new = G_malloc(sizeof(struct rc));
+
+ if (!new)
+ G_fatal_error(_("rclist out of memory"));
+
+ new->next = NULL;
+ new->row = row;
+ new->col = col;
+
+ if (list->head) {
+ list->head->next = new;
+ list->head = new;
+ }
+ else {
+ list->head = list->tail = new;
+ }
+
+ return;
+}
+
+/* return 1 if an element was dropped
+ * return 0 if list is empty
+ */
+int rclist_drop(struct rclist *list, struct rc *rc)
+{
+ if (list->tail) {
+ struct rc *next = list->tail->next;
+
+ rc->row = list->tail->row;
+ rc->col = list->tail->col;
+ G_free(list->tail);
+ list->tail = next;
+ if (!list->tail)
+ list->head = NULL;
+
+ return 1;
+ }
+
+ return 0;
+}
+
+void rclist_destroy(struct rclist *list)
+{
+ struct rc *next = list->tail;
+
+ while (next) {
+ next = next->next;
+ G_free(list->tail);
+ list->tail = next;
+ }
+ list->head = NULL;
+
+ return;
+}
+
Property changes on: grass-addons/grass7/vector/v.surf.tps/rclist.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.surf.tps/rclist.h
===================================================================
--- grass-addons/grass7/vector/v.surf.tps/rclist.h (rev 0)
+++ grass-addons/grass7/vector/v.surf.tps/rclist.h 2016-08-25 18:50:34 UTC (rev 69257)
@@ -0,0 +1,20 @@
+/* row/col list */
+
+struct rc
+{
+ struct rc *next;
+ int row;
+ int col;
+};
+
+struct rclist
+{
+ struct rc *tail, *head;
+};
+
+/* rclist.c */
+void rclist_init(struct rclist *);
+void rclist_add(struct rclist *, int, int);
+int rclist_drop(struct rclist *, struct rc *);
+void rclist_destroy(struct rclist *);
+
Property changes on: grass-addons/grass7/vector/v.surf.tps/rclist.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.surf.tps/tps.c
===================================================================
--- grass-addons/grass7/vector/v.surf.tps/tps.c (rev 0)
+++ grass-addons/grass7/vector/v.surf.tps/tps.c 2016-08-25 18:50:34 UTC (rev 69257)
@@ -0,0 +1,909 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#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 "tps.h"
+#include "flag.h"
+#include "rclist.h"
+
+static int solvemat(double **m, double a[], double B[], int n)
+{
+ int i, j, i2, j2, imark;
+ double factor, temp, *tempp;
+ double pivot; /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
+
+ for (i = 0; i < n; i++) {
+ j = i;
+
+ /* find row with largest magnitude value for pivot value */
+
+ pivot = m[i][j];
+ imark = i;
+ for (i2 = i + 1; i2 < n; i2++) {
+ temp = fabs(m[i2][j]);
+ if (temp > fabs(pivot)) {
+ pivot = m[i2][j];
+ imark = i2;
+ }
+ }
+
+ /* if the pivot is very small then the points are nearly co-linear */
+ /* co-linear points result in an undefined matrix, and nearly */
+ /* co-linear points results in a solution with rounding error */
+
+ if (fabs(pivot) < GRASS_EPSILON) {
+ G_debug(4, "Matrix is unsolvable: pivot = %g", pivot);
+ return 0;
+ }
+
+ /* if row with highest pivot is not the current row, switch them */
+
+ if (imark != i) {
+ /*
+ for (j2 = 0; j2 < n; j2++) {
+ temp = m[imark][j2];
+ m[imark][j2] = m[i][j2];
+ m[i][j2] = temp;
+ }
+ */
+
+ tempp = m[imark];
+ m[imark] = m[i];
+ m[i] = tempp;
+
+ temp = a[imark];
+ a[imark] = a[i];
+ a[i] = temp;
+ }
+
+ /* compute zeros above and below the pivot, and compute
+ values for the rest of the row as well */
+
+ for (i2 = 0; i2 < n; i2++) {
+ if (i2 != i) {
+ factor = m[i2][j] / pivot;
+ for (j2 = j; j2 < n; j2++)
+ m[i2][j2] -= factor * m[i][j2];
+ a[i2] -= factor * a[i];
+ }
+ }
+ }
+
+ /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
+ COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
+
+ for (i = 0; i < n; i++) {
+ B[i] = a[i] / m[i][i];
+ }
+
+ return 1;
+}
+
+static int load_tps_pnts(struct tps_pnt *pnts, int n_points, int n_vars,
+ double regularization, double **m, double *a)
+{
+ int i, j, k;
+ double dx, dy, dist, dist2, distsum, reg;
+
+ for (i = 0; i <= n_vars; i++) {
+ a[i] = 0.0;
+ for (j = 0; j <= n_vars; j++) {
+ m[i][j] = m[j][i] = 0.0;
+ }
+ }
+
+ distsum = 0;
+ for (i = 0; i < n_points; i++) {
+ /* global */
+ a[i + 1 + n_vars] = pnts[i].val;
+ 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];
+ }
+
+ /* TPS */
+ for (k = 0; k < n_points; k++) {
+ dx = (pnts[i].c -pnts[k].c) * 2.0;
+ dy = (pnts[i].r -pnts[k].r) * 2.0;
+
+ dist2 = dx * dx + dy * dy;
+ dist = 0;
+ if (dist2 > 0)
+ dist = dist2 * log(dist2) * 0.5;
+
+ m[i + 1 + n_vars][k + 1 + n_vars] = m[k + 1 + n_vars][i + 1 + n_vars] = dist;
+ if (regularization > 0 && dist2 > 0)
+ distsum += sqrt(dist2);
+ }
+ }
+
+ if (regularization > 0) {
+ distsum /= (n_points * n_points);
+ reg = regularization * distsum * distsum;
+
+ for (i = 0; i < n_points; i++)
+ m[i + 1 + n_vars][i + 1 + n_vars] = reg;
+ }
+
+ return 1;
+}
+
+int map_tps(int out_fd, int *var_fd, int n_vars, int mask_fd,
+ struct tps_pnt *pnts, int n_points, double regularization)
+{
+ 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;
+
+ G_message(_("Global TPS interpolation with %d points..."), n_points);
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ /* 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));
+
+ 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();
+ }
+ }
+
+ maskbuf = NULL;
+ if (mask_fd >= 0)
+ maskbuf = Rast_allocate_c_buf();
+
+ /* load points to matrix */
+ load_tps_pnts(pnts, n_points, n_vars, regularization, m, a);
+
+ /* solve */
+ if (!solvemat(m, a, B, n_points + 1 + n_vars))
+ G_fatal_error(_("Matrix is unsolvable"));
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+
+ 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);
+
+ for (col = 0; col < ncols; col++) {
+
+ if (maskbuf &&
+ (Rast_is_c_null_value(&maskbuf[col]) || maskbuf[col] == 0)) {
+
+ Rast_set_d_null_value(&outbuf[col], 1);
+
+ continue;
+ }
+
+ if (n_vars) {
+ int isnull = 0;
+
+ for (j = 0; j < n_vars; j++) {
+ if (Rast_is_d_null_value(&dbuf[j][col])) {
+ isnull = 1;
+ break;
+ }
+ }
+ if (isnull)
+ continue;
+ }
+
+ 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;
+
+ dist2 = dx * dx + dy * dy;
+ dist = 0;
+ if (dist2 > 0) {
+ dist = dist2 * log(dist2) * 0.5;
+ result += B[1 + n_vars + i] * dist;
+ }
+ }
+ 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]);
+ }
+ G_free(dbuf);
+ }
+ if (maskbuf)
+ G_free(maskbuf);
+ G_free(outbuf);
+
+ return 1;
+}
+
+static int cmp_rc(const void *first, const void *second)
+{
+ struct rc *a = (struct rc *)first, *b = (struct rc *)second;
+
+ if (a->row == b->row)
+ return (a->col - b->col);
+
+ return (a->row - b->row);
+}
+
+static int bfs_search(SEGMENT *p_seg, FLAG *mask_flag, int *bfsid,
+ 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 n, found;
+ int pid;
+ struct rc next, ngbr_rc;
+ struct rclist rilist;
+ struct RB_TREE *visited;
+ double dx, dy, dist;
+
+
+ visited = rbtree_create(cmp_rc, sizeof(struct rc));
+ ngbr_rc.row = row;
+ ngbr_rc.col = col;
+ rbtree_insert(visited, &ngbr_rc);
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ /* breadth-first search */
+ next.row = row;
+ next.col = col;
+ rclist_init(&rilist);
+
+ found = 0;
+
+ 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 (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;
+
+ if (rbtree_find(visited, &ngbr_rc))
+ continue;
+
+ rbtree_insert(visited, &ngbr_rc);
+
+ Segment_get(p_seg, (void *)&pid, rown, coln);
+
+ if (pid < 0) {
+ rclist_add(&rilist, rown, coln);
+ }
+ else if (rown >= rmin && rown <= rmax && coln >= cmin && coln <= cmax) {
+ dx = coln - col;
+ dy = rown - row;
+
+ dist = dx * dx + dy * dy;
+
+ if (dist > min_dist) {
+ bfsid[found] = pid;
+ found++;
+ }
+ }
+ if (found == max_points)
+ break;
+ }
+ if (found == max_points)
+ break;
+ } while (rclist_drop(&rilist, &next)); /* while there are cells to check */
+
+
+ rclist_destroy(&rilist);
+ rbtree_destroy(visited);
+
+ return found;
+}
+
+int cell_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)
+{
+ int row, col, nrows, ncols;
+ double **m, *a, *B;
+ int i, j;
+ int kdalloc, palloc, n_cur_points;
+ struct tps_pnt *cur_pnts;
+ double dx, dy, dist, dist2;
+ DCELL **dbuf, result, *outbuf, *varbuf;
+ CELL *maskbuf;
+ int solved;
+ struct kdtree *kdt;
+ double c[2];
+ int *kduid, kdfound;
+ double *kddist;
+ int *bfsid, bfsfound, pfound;
+
+ 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;
+ double weight, dxi, dyi;
+ double wmin, wmax;
+ int rmin, rmax, cmin, cmax, rminp, rmaxp, cminp, cmaxp;
+ int irow, irow1, irow2, icol, icol1, icol2;
+ int rcnt;
+ int nskipped;
+ double pthin2;
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ kdalloc = palloc = min_points;
+ a = G_malloc((palloc + 1 + n_vars) * sizeof(double));
+ B = G_malloc((palloc + 1 + n_vars) * sizeof(double));
+ m = G_malloc((palloc + 1 + n_vars) * sizeof(double *));
+ 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));
+
+ bfsid = G_malloc(min_points * sizeof(int));
+
+ mask_flag = flag_create(nrows, ncols);
+
+ maskbuf = NULL;
+ if (mask_fd >= 0) {
+ maskbuf = Rast_allocate_c_buf();
+
+ for (row = 0; row < nrows; row++) {
+ Rast_get_c_row(mask_fd, maskbuf, row);
+ for (col = 0; col < ncols; col++) {
+ if (Rast_is_c_null_value(&maskbuf[col]) || maskbuf[col] == 0)
+ FLAG_SET(mask_flag, row, col);
+ }
+ }
+ G_free(maskbuf);
+ maskbuf = NULL;
+ }
+
+ /* create and init SEG structs */
+ varsize = n_vars * sizeof(DCELL);
+ segsize = (double)64 * 64 * (sizeof(struct tps_out) + varsize);
+ 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");
+ }
+
+ if (Segment_open(&p_seg, G_tempfile(), nrows, ncols, 64, 64,
+ sizeof(int), nsegs) != 1) {
+ G_fatal_error("Unable to create input temporary files");
+ }
+
+ 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();
+ }
+
+ if (Segment_open(&var_seg, G_tempfile(), nrows, ncols, 64, 64,
+ varsize, nsegs) != 1) {
+ G_fatal_error("Unable to create input temporary files");
+ }
+ }
+
+ if (do_bfs) {
+ 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 (i = 0; i < n_points; i++) {
+ G_percent(i, n_points, 2);
+ if (Segment_put(&p_seg, (void *)&i, pnts[i].r, pnts[i].c) != 1)
+ G_fatal_error(_("Unable to write to temporary file"));
+ }
+ }
+
+ G_message(_("Initializing output ..."));
+ tps_out.val = 0;
+ tps_out.wsum = 0;
+ 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)
+ 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);
+
+ 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);
+
+ /* G_message(_("Processing %d rows, current row:"), nrows); */
+
+ wmin = 10;
+ wmax = 0;
+
+ row = -1;
+ rcnt = 0;
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 1);
+ /* fprintf(stderr, "%6d\b\b\b\b\b\b", row); */
+ rcnt++;
+
+ for (col = 0; col < ncols; col++) {
+
+ if ((FLAG_GET(mask_flag, row, col))) {
+ continue;
+ }
+
+ 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)
+ continue;
+ }
+
+ 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) {
+
+ /* increase points if unsolvable */
+ n_cur_points += min_points;
+ if (n_cur_points > n_points)
+ n_cur_points = n_points;
+
+ /* 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));
+ }
+
+ /* collect nearest neighbors */
+ kdfound = kdtree_knn(kdt, c, kduid, kddist, n_cur_points, NULL);
+
+ /* load points to matrix */
+ rmin = nrows;
+ rmax = 0;
+ cmin = ncols;
+ 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;
+ }
+ pfound = kdfound;
+
+ if (do_bfs) {
+ /* collect points with breadth-first search
+ * min dist must be > max dist of nearest neighbors
+ * maximum min_points with dist > max dist */
+ if ((rminp < row && rmin >= row) ||
+ (cmaxp > col && cmin <= col)) {
+ bfsfound = bfs_search(&p_seg, mask_flag, bfsid,
+ row, col, n_cur_points / 4,
+ kddist[kdfound - 1],
+ 0, row - 1, col, ncols - 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]];
+ if (rmin > cur_pnts[pfound + i].r)
+ rmin = cur_pnts[pfound + i].r;
+ if (rmax < cur_pnts[pfound + i].r)
+ rmax = cur_pnts[pfound + i].r;
+ if (cmin > cur_pnts[pfound + i].c)
+ cmin = cur_pnts[pfound + i].c;
+ if (cmax < cur_pnts[pfound + i].c)
+ cmax = cur_pnts[pfound + i].c;
+ }
+ pfound += bfsfound;
+ }
+
+ if ((rminp < row && rmin >= row) ||
+ (cminp < col && cmin >= col)) {
+ bfsfound = bfs_search(&p_seg, mask_flag, bfsid,
+ row, col, n_cur_points / 4,
+ kddist[kdfound - 1],
+ 0, row - 1, 0, col - 1);
+
+ 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]];
+ if (rmin > cur_pnts[pfound + i].r)
+ rmin = cur_pnts[pfound + i].r;
+ if (rmax < cur_pnts[pfound + i].r)
+ rmax = cur_pnts[pfound + i].r;
+ if (cmin > cur_pnts[pfound + i].c)
+ cmin = cur_pnts[pfound + i].c;
+ if (cmax < cur_pnts[pfound + i].c)
+ cmax = cur_pnts[pfound + i].c;
+ }
+ pfound += bfsfound;
+ }
+
+ if ((rmaxp > row && rmax <= row) ||
+ (cminp < col && cmin >= col)) {
+ bfsfound = bfs_search(&p_seg, mask_flag, bfsid,
+ row, col, n_cur_points / 4,
+ kddist[kdfound - 1],
+ row, nrows - 1, 0, col - 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]];
+ if (rmin > cur_pnts[pfound + i].r)
+ rmin = cur_pnts[pfound + i].r;
+ if (rmax < cur_pnts[pfound + i].r)
+ rmax = cur_pnts[pfound + i].r;
+ if (cmin > cur_pnts[pfound + i].c)
+ cmin = cur_pnts[pfound + i].c;
+ if (cmax < cur_pnts[pfound + i].c)
+ cmax = cur_pnts[pfound + i].c;
+ }
+ pfound += bfsfound;
+ }
+
+ if ((rmaxp > row && rmax <= row) ||
+ (cmaxp > col && cmin <= col)) {
+ bfsfound = bfs_search(&p_seg, mask_flag, bfsid,
+ row, col, n_cur_points / 4,
+ kddist[kdfound - 1],
+ row, nrows - 1, col, ncols - 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]];
+ if (rmin > cur_pnts[pfound + i].r)
+ rmin = cur_pnts[pfound + i].r;
+ if (rmax < cur_pnts[pfound + i].r)
+ rmax = cur_pnts[pfound + i].r;
+ if (cmin > cur_pnts[pfound + i].c)
+ cmin = cur_pnts[pfound + i].c;
+ if (cmax < cur_pnts[pfound + i].c)
+ cmax = cur_pnts[pfound + i].c;
+ }
+ pfound += bfsfound;
+ }
+ }
+
+ if (palloc < pfound) {
+ G_free(a);
+ G_free(B);
+ for (i = 0; i < (palloc + 1 + n_vars); i++)
+ G_free(m[i]);
+ G_free(m);
+
+ palloc = pfound;
+ a = G_malloc((palloc + 1 + n_vars) * sizeof(double));
+ B = G_malloc((palloc + 1 + n_vars) * sizeof(double));
+ m = G_malloc((palloc + 1 + n_vars) * sizeof(double *));
+ for (i = 0; i < (palloc + 1 + n_vars); i++)
+ 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);
+
+ /* solve */
+ solved = solvemat(m, a, B, pfound + 1 + n_vars);
+
+ if (!solved && n_cur_points == n_points)
+ 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;
+ }
+
+ if (rmin == rminp)
+ rmin = 0;
+ if (rmax == rmaxp)
+ rmax = nrows - 1;
+ if (cmin == cminp)
+ cmin = 0;
+ if (cmax == cmaxp)
+ cmax = ncols - 1;
+
+ irow1 = rmin + (rmax - rmin) * overlap / 2.5;
+ irow2 = rmax - (rmax - rmin) * overlap / 2.5;
+ icol1 = cmin + (cmax - cmin) * overlap / 2.5;
+ icol2 = cmax - (cmax - cmin) * overlap / 2.5;
+
+ if (rmax == rmaxp)
+ irow2 = nrows - 1;
+ if (cmax == cmaxp)
+ icol2 = ncols - 1;
+
+ if (irow1 > row) {
+ irow2 -= irow1 - row;
+ irow1 = row;
+ }
+ if (irow2 < row) {
+ irow1 += row - irow2;
+ irow2 = row;
+ }
+ if (icol1 > col) {
+ icol2 -= icol1 - col;
+ icol1 = col;
+ }
+ if (icol2 < col) {
+ icol1 += col - icol2;
+ icol2 = col;
+ }
+
+ dx = icol2 - icol1;
+ dy = irow2 - irow1;
+
+ dxi = icol2 - icol1 + 1;
+ dyi = irow2 - irow1 + 1;
+
+ for (irow = irow1; irow <= irow2; irow++) {
+ 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)
+ continue;
+ }
+
+ Segment_get(&out_seg, (void *)&tps_out, irow, icol);
+
+ result = B[0];
+ if (n_vars) {
+ Segment_get(&var_seg, (void *)varbuf, irow, icol);
+ for (j = 0; j < n_vars; j++)
+ result += varbuf[j] * B[j + 1];
+ }
+
+ for (i = 0; i < pfound; i++) {
+ dx = (cur_pnts[i].c - icol) * 2.0;
+ dy = (cur_pnts[i].r - irow) * 2.0;
+
+ dist2 = dx * dx + dy * dy;
+ dist = 0;
+ if (dist2 > 0) {
+ dist = dist2 * log(dist2) * 0.5;
+ result += B[1 + n_vars + i] * dist;
+ }
+ }
+ weight = 1;
+
+ dx = abs(icol - col) / dxi;
+ dy = abs(irow - row) / dyi;
+ weight = (dx * dx + dy * dy) / 2;
+
+ weight = exp(-weight * weight * 4);
+
+ if (wmin > weight)
+ 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);
+ }
+ }
+ }
+ }
+ G_percent(1, 1, 1);
+
+ G_debug(1, "wmin: %g", wmin);
+ G_debug(1, "wmax: %g", wmax);
+ G_debug(1, "rcnt: %d, nrows: %d", rcnt, nrows);
+
+ if (n_vars)
+ Segment_close(&var_seg);
+
+ if (do_bfs)
+ Segment_close(&p_seg);
+
+ outbuf = Rast_allocate_d_buf();
+
+ G_message(_("Writing output..."));
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+
+ for (col = 0; col < ncols; col++) {
+
+ if ((FLAG_GET(mask_flag, row, col))) {
+ Rast_set_d_null_value(&outbuf[col], 1);
+ continue;
+ }
+
+ Segment_get(&out_seg, (void *)&tps_out, row, col);
+
+ if (tps_out.wsum == 0)
+ Rast_set_d_null_value(&outbuf[col], 1);
+ else
+ outbuf[col] = tps_out.val / tps_out.wsum;
+ }
+ Rast_put_d_row(out_fd, outbuf);
+ }
+ G_percent(1, 1, 2);
+
+ G_free(outbuf);
+ Segment_close(&out_seg);
+
+ return 1;
+}
Property changes on: grass-addons/grass7/vector/v.surf.tps/tps.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.surf.tps/tps.h
===================================================================
--- grass-addons/grass7/vector/v.surf.tps/tps.h (rev 0)
+++ grass-addons/grass7/vector/v.surf.tps/tps.h 2016-08-25 18:50:34 UTC (rev 69257)
@@ -0,0 +1,14 @@
+
+struct tps_pnt
+{
+ int r, c; /* row, col in target window */
+ double val; /* value to interpolate */
+ double *vars; /* values of covariables */
+};
+
+int map_tps(int out_fd, int *var_fd, int n_vars, int mask_fd,
+ struct tps_pnt *pnts, int n_points, double regularization);
+int cell_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);
Property changes on: grass-addons/grass7/vector/v.surf.tps/tps.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.surf.tps/v.surf.tps.html
===================================================================
--- grass-addons/grass7/vector/v.surf.tps/v.surf.tps.html (rev 0)
+++ grass-addons/grass7/vector/v.surf.tps/v.surf.tps.html 2016-08-25 18:50:34 UTC (rev 69257)
@@ -0,0 +1,104 @@
+<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.
+
+<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>
+
+<h3>Basic interpolation</h3>
+Interpolation of 30 year precipitation normals in the North Carlolina
+sample dataset:
+
+<div class="code"><pre>
+g.region -p rast=elev_state_500m
+v.surf.tps input=precip_30ynormals_3d output=precip_30ynormals_3d 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 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 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 min=20 \
+ overlap=0.1 covars=elev_state_500m smooth=0.1
+</pre></div>
+
+
+<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>
Property changes on: grass-addons/grass7/vector/v.surf.tps/v.surf.tps.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list