[GRASS-SVN] r70362 - grass-addons/grass7/raster/r.resamp.tps
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jan 13 06:10:15 PST 2017
Author: mmetz
Date: 2017-01-13 06:10:15 -0800 (Fri, 13 Jan 2017)
New Revision: 70362
Modified:
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:
r.resamp.tps: more safety checks, Rast_align_window() not working in latlon
Modified: grass-addons/grass7/raster/r.resamp.tps/main.c
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/main.c 2017-01-13 05:56:55 UTC (rev 70361)
+++ grass-addons/grass7/raster/r.resamp.tps/main.c 2017-01-13 14:10:15 UTC (rev 70362)
@@ -35,17 +35,18 @@
struct GModule *module;
struct Option *in_opt, *ivar_opt, *ovar_opt, *out_opt, *minpnts_opt,
- *radius_opt, *reg_opt, *ov_opt, *mask_opt, *mem_opt;
+ *maxpnts_opt, *radius_opt, *reg_opt, *ov_opt,
+ *lm_opt, *mask_opt, *mem_opt;
struct Flag *c_flag;
struct Cell_head cellhd, src, dst;
int n_ivars, n_ovars, n_vars;
off_t n_points;
- int min_points, radius;
+ int min_points, max_points, radius;
int r, c, nrows, ncols;
DCELL **dbuf, *dval;
- double regularization, overlap;
+ double regularization, overlap, lm_thresh;
SEGMENT in_seg, var_seg, out_seg;
int insize, varsize;
double segsize;
@@ -76,7 +77,7 @@
ov_opt->key = "overlap";
ov_opt->type = TYPE_DOUBLE;
ov_opt->required = NO;
- ov_opt->answer = "0.8";
+ ov_opt->answer = "0.2";
ov_opt->label =
_("Overlap factor <= 1");
ov_opt->description =
@@ -87,11 +88,19 @@
minpnts_opt->key = "min";
minpnts_opt->type = TYPE_DOUBLE;
minpnts_opt->required = NO;
- minpnts_opt->answer = "20";
+ minpnts_opt->answer = "100";
minpnts_opt->description =
_("Minimum number of points to use for TPS interpolation");
minpnts_opt->guisection = _("Settings");
+ maxpnts_opt = G_define_option();
+ maxpnts_opt->key = "max";
+ maxpnts_opt->type = TYPE_DOUBLE;
+ maxpnts_opt->required = NO;
+ maxpnts_opt->description =
+ _("Maximum number of points to use for TPS interpolation");
+ maxpnts_opt->guisection = _("Settings");
+
radius_opt = G_define_option();
radius_opt->key = "radius";
radius_opt->type = TYPE_INTEGER;
@@ -117,6 +126,17 @@
_("Name of input raster map(s) to use as covariables matching the current region");
ovar_opt->guisection = _("Settings");
+ lm_opt = G_define_option();
+ lm_opt->key = "lmfilter";
+ lm_opt->type = TYPE_DOUBLE;
+ lm_opt->required = NO;
+ lm_opt->answer = "0";
+ lm_opt->label =
+ _("Threshold to avoid interpolation outliers when using covariables");
+ lm_opt->description =
+ _("Disabled when set to zero");
+ lm_opt->guisection = _("Settings");
+
out_opt = G_define_standard_option(G_OPT_R_OUTPUT);
out_opt->key = "output";
out_opt->required = YES;
@@ -138,12 +158,15 @@
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);
+ if (!minpnts_opt->answer && !radius_opt->answer)
+ G_fatal_error(_("Either <%s> or <%s> must be given"),
+ minpnts_opt->key, radius_opt->key);
+
outname = out_opt->answer;
n_ivars = 0;
@@ -193,8 +216,22 @@
/* align dst to cellhd */
src = dst;
+ /*
Rast_align_window(&src, &cellhd);
+ */
+ src.ns_res = cellhd.ns_res;
+ src.ew_res = cellhd.ew_res;
+ src.south =
+ cellhd.north - ceil((cellhd.north - src.south) / cellhd.ns_res) * cellhd.ns_res;
+ src.north =
+ cellhd.north - floor((cellhd.north - src.north) / cellhd.ns_res) * cellhd.ns_res;
+ src.east =
+ cellhd.west + ceil((src.east - cellhd.west) / cellhd.ew_res) * cellhd.ew_res;
+ src.west =
+ cellhd.west + floor((src.west - cellhd.west) / cellhd.ew_res) * cellhd.ew_res;
+
+
/* open segment structures for input and output */
/* set input window */
@@ -329,16 +366,30 @@
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);
+ min_points = 0;
+ if (minpnts_opt->answer) {
+ 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);
+ }
}
+ max_points = 0;
+ if (maxpnts_opt->answer) {
+ max_points = atoi(maxpnts_opt->answer);
+ if (max_points < min_points) {
+ G_warning(_("Maximum number of points must be equal to or larger than minimum number of points, disabling"));
+ max_points = 0;
+ }
+ }
- radius = atoi(radius_opt->answer);
- if (radius < 0)
- radius = 0;
+ radius = 0;
+ if (radius_opt->answer) {
+ radius = atoi(radius_opt->answer);
+ if (radius < 0)
+ radius = 0;
+ }
regularization = atof(reg_opt->answer);
if (regularization < 0)
@@ -350,19 +401,26 @@
if (overlap > 1)
overlap = 1;
+ lm_thresh = 0;
+ if (lm_opt->answer) {
+ lm_thresh = atof(lm_opt->answer);
+ if (lm_thresh < 0)
+ lm_thresh = 0;
+ }
+
if (radius) {
if (tps_window(&in_seg, &var_seg, n_vars, &out_seg, out_fd,
mask_opt->answer, &src, &dst, n_points,
regularization, overlap,
- radius) != 1) {
+ radius, lm_thresh) != 1) {
G_fatal_error(_("TPS interpolation failed"));
}
}
else {
if (tps_nn(&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) {
+ min_points, max_points, regularization, overlap,
+ c_flag->answer, lm_thresh) != 1) {
G_fatal_error(_("TPS interpolation failed"));
}
}
Modified: grass-addons/grass7/raster/r.resamp.tps/tps.c
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/tps.c 2017-01-13 05:56:55 UTC (rev 70361)
+++ grass-addons/grass7/raster/r.resamp.tps/tps.c 2017-01-13 14:10:15 UTC (rev 70362)
@@ -11,10 +11,6 @@
#include "rclist.h"
#include "pavl.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;
@@ -40,7 +36,7 @@
/* 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) {
+ if (fabs(pivot) < 1.0e-7) {
G_debug(4, "Matrix is unsolvable: pivot = %g", pivot);
return 0;
}
@@ -109,17 +105,15 @@
return (a->r - b->r);
}
-static int load_tps_pnts(SEGMENT *in_seg, int n_vars,
+static int load_tps_pnts(SEGMENT *in_seg, DCELL *dval, int n_vars,
struct tps_pnt *pnts, int n_points,
- struct Cell_head *src, struct Cell_head *dst,
- double regularization, double **m, double *a)
+ struct Cell_head *src, struct Cell_head *dst,
+ double regularization, double **m, double *a,
+ double **mvars, double *avars)
{
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;
m[i][i] = 0.0;
@@ -128,6 +122,16 @@
}
}
+ if (n_vars) {
+ for (i = 0; i <= n_vars; i++) {
+ avars[i] = 0.0;
+ mvars[i][i] = 0.0;
+ for (j = 0; j < i; j++) {
+ mvars[i][j] = mvars[j][i] = 0.0;
+ }
+ }
+ }
+
distsum = 0;
for (i = 0; i < n_points; i++) {
@@ -141,25 +145,30 @@
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
+ if (n_vars) {
+ avars[0] += dval[0];
+ mvars[0][0] += 1;
+ for (j = 1; j <= n_vars; j++) {
+ mvars[j][0] += dval[j];
+ mvars[0][j] = mvars[j][0];
+ mvars[j][j] += dval[j] * dval[j];
+ for (k = 1; k < j; k++) {
+ mvars[j][k] += dval[j] * dval[k];
+ mvars[k][j] = mvars[j][k];
+ }
+ avars[j] += dval[0] * dval[j];
+ }
+ }
+
+ /* convert src r,c to n,s */
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)
@@ -179,8 +188,6 @@
}
}
- G_free(dval);
-
if (regularization > 0) {
distsum /= (n_points * n_points);
reg = regularization * distsum * distsum;
@@ -469,21 +476,117 @@
return found;
}
+static double interp_wa(SEGMENT *in_seg, DCELL *dval,
+ struct tps_pnt *cur_pnts, int pfound,
+ int row, int col,
+ struct Cell_head *src, struct Cell_head *dst,
+ double distmax, double *tps_weight)
+{
+ int i;
+ int src_row, src_col;
+ double wsum, maxweight, weight, result;
+ double i_n, i_e, dxyw;
+ double dx, dy, dist2;
+
+ G_fatal_error("Weighted average");
+
+ i_n = dst->north - (row + 0.5) * dst->ns_res;
+ i_e = dst->west + (col + 0.5) * dst->ew_res;
+ dxyw = distmax * (src->ns_res + src->ew_res) / 2.0;
+
+ result = 0;
+ wsum = 0;
+ maxweight = 0;
+ for (i = 0; i < pfound; i++) {
+
+ src_row = (src->north - cur_pnts[i].r) / src->ns_res;
+ src_col = (cur_pnts[i].c - src->west) / src->ew_res;
+
+ Segment_get(in_seg, (void *)dval, src_row, src_col);
+
+ dx = cur_pnts[i].c - i_e;
+ dy = cur_pnts[i].r - i_n;
+
+ /* weight for average */
+ dx = fabs(dx) / dxyw;
+ dy = fabs(dy) / dxyw;
+ dist2 = dx * dx + dy * dy;
+ weight = exp(-dist2 * 4.0);
+
+ /* weight for tps */
+ if (maxweight < weight)
+ maxweight = weight;
+
+ wsum += weight;
+ result += dval[0] * weight;
+ }
+ result /= wsum;
+ *tps_weight = maxweight;
+
+ return result;
+}
+
+static double lm_rsqr(SEGMENT *in_seg, int n_vars, struct Cell_head *src,
+ struct tps_pnt *cur_pnts, int pfound, double *B)
+{
+ int i, j, r, c;
+ double *obs, *est;
+ DCELL *dval;
+ double sumX, sumY, sumsqX, sumsqY, sumXY, R2;
+
+ obs = G_malloc(pfound * sizeof(double));
+ est = G_malloc(pfound * sizeof(double));
+ dval = G_malloc((1 + n_vars) * sizeof(DCELL));
+
+ sumX = sumY = sumsqX = sumsqY = sumXY = 0.0;
+ for (i = 0; i < pfound; i++) {
+
+ r = (int)((src->north - cur_pnts[i].r) / src->ns_res);
+ c = (int)((cur_pnts[i].c - src->west) / src->ew_res);
+
+ Segment_get(in_seg, (void *)dval, r, c);
+ obs[i] = dval[0];
+ est[i] = B[0];
+ for (j = 1; j <= n_vars; j++) {
+ est[i] += dval[j] * B[j];
+ }
+
+ sumX += obs[i];
+ sumY += est[i];
+ sumsqX += obs[i] * obs[i];
+ sumsqY += est[i] * est[i];
+ sumXY += obs[i] * est[i];
+ }
+
+ R2 = (sumXY * pfound - sumX * sumY) * (sumXY * pfound - sumX * sumY) /
+ ((sumsqX * pfound - sumX * sumX) * (sumsqY * pfound - sumY * sumY));
+
+ G_free(obs);
+ G_free(est);
+ G_free(dval);
+
+ return R2;
+}
+
+
int tps_nn(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)
+ off_t n_points, int min_points, int max_points,
+ double regularization, double overlap, int clustered,
+ double lm_thresh)
{
int ridx, cidx, row, col, nrows, ncols, src_row, src_col;
double **m, *a, *B;
+ double **mvars, *avars, *Bvars;
int i, j;
int kdalloc, palloc, n_cur_points;
struct tps_pnt *cur_pnts;
double dx, dy, dist, dist2, mfactor;
DCELL *dval, result, *outbuf, *varbuf;
CELL *maskbuf;
- int solved;
+ int solved, solved_tps_lm, solved_tps, solved_lm, n_vars_i;
+ double rsqr;
int kdfound, bfsfound, pfound;
double distmax, mindist;
int do_clustered;
@@ -495,9 +598,8 @@
double wmin, wmax;
int rmin, rmax, cmin, cmax, rminp, rmaxp, cminp, cmaxp;
int irow, irow1, irow2, icol, icol1, icol2;
-#ifndef USE_RC
+ unsigned int cnt_wa, cnt_tps_lm, cnt_tps;
double i_n, i_e;
-#endif
nrows = Rast_window_rows();
ncols = Rast_window_cols();
@@ -512,6 +614,9 @@
pvar = NULL;
varbuf = NULL;
+ mvars = NULL;
+ avars = NULL;
+ Bvars = NULL;
if (n_vars) {
pvar = G_malloc(palloc * 5 * n_vars * sizeof(double));
cur_pnts[0].vars = pvar;
@@ -519,6 +624,12 @@
cur_pnts[i].vars = cur_pnts[i - 1].vars + n_vars;
varbuf = G_malloc(n_vars * sizeof(DCELL));
+
+ avars = G_malloc((1 + n_vars) * sizeof(double));
+ Bvars = G_malloc((1 + n_vars) * sizeof(double));
+ mvars = G_malloc((1 + n_vars) * sizeof(double *));
+ for (i = 0; i < (1 + n_vars); i++)
+ mvars[i] = G_malloc((1 + n_vars) * sizeof(double));
}
dval = G_malloc((1 + n_vars) * sizeof(DCELL));
@@ -592,6 +703,9 @@
wmin = 10;
wmax = 0;
+ cnt_wa = 0;
+ cnt_tps_lm = 0;
+ cnt_tps = 0;
if (overlap > 1.0)
overlap = 1.0;
@@ -632,15 +746,29 @@
if (tps_out.wmax > overlap)
continue;
+ rmin = src->rows;
+ rmax = 0;
+ cmin = src->cols;
+ cmax = 0;
+
solved = 0;
+ solved_tps_lm = solved_tps = 0;
n_cur_points = 0;
- kdfound = 0;
+ pfound = kdfound = 0;
src_col = col_src2dst(col, dst, src);
-
- while (!solved) {
+ distmax = 0;
+ n_vars_i = n_vars;
+
+ while (!solved &&
+ (max_points == 0 || n_cur_points < max_points)) {
/* increase points if unsolvable */
n_cur_points += min_points;
+
+ if (n_cur_points > min_points)
+ G_verbose_message(_("Increasing number of points to %d"),
+ n_cur_points);
+
if (n_cur_points > n_points)
n_cur_points = n_points;
@@ -769,16 +897,103 @@
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);
+ load_tps_pnts(in_seg, dval, n_vars, cur_pnts, pfound,
+ src, dst, regularization, m, a, mvars, avars);
/* solve */
- solved = solvemat(m, a, B, pfound + 1 + n_vars);
+ /* it can happen that a simple linear model with the
+ * given covariables can be solved but TPS with
+ * covariables can not be solved
+ */
+ solved_tps_lm = solved_tps = 0;
+ n_vars_i = n_vars;
+ if (n_vars) {
+ solved_tps_lm = solvemat(m, a, B, pfound + 1 + n_vars);
+
+ if (solved_tps_lm && lm_thresh > 0) {
+
+ solved_lm = solvemat(mvars, avars, Bvars, 1 + n_vars);
+ if (!solved_lm) {
+ G_debug(1, "LM with covariables not working at row %d, col %d",
+ row, col);
+
+ solved_tps_lm = 0;
+ }
+ else {
+ rsqr = lm_rsqr(in_seg, n_vars, src, cur_pnts, pfound, B);
+
+ if (rsqr < lm_thresh) {
+ solved_tps_lm = 0;
+ }
+ else {
+ for (i = 1; i <= n_vars; i++) {
+ if (fabs(B[i]) > fabs(5 * Bvars[i])) {
+ G_debug(1, "LM B%d is %g but TPS B%d is %g",
+ i, Bvars[i], i, B[i]);
+ solved_tps_lm = 0;
+ }
+ }
+ }
+ }
+ }
+
+ if (!solved_tps_lm) {
+ n_vars_i = 0;
+ for (i = 0; i < pfound; i++) {
+ cur_pnts[i].r = (int)((src->north - cur_pnts[i].r) / src->ns_res);
+ cur_pnts[i].c = (int)((cur_pnts[i].c - src->west) / src->ew_res);
+ }
+ load_tps_pnts(in_seg, dval, 0, cur_pnts, pfound,
+ src, dst, regularization, m, a, NULL, NULL);
+ }
+ }
+
+ if (!solved_tps_lm) {
+ solved_tps = solvemat(m, a, B, pfound + 1);
+ }
+
+ solved = (solved_tps_lm || solved_tps);
+
if (!solved && n_cur_points == n_points)
G_fatal_error(_("Matrix is unsolvable"));
}
+ /* priorities for interpolation
+ * TPS with covariables
+ * TPS without covariables
+ * weighted average
+ */
+
+ if (solved_tps_lm)
+ cnt_tps_lm++;
+ else if (solved_tps)
+ cnt_tps++;
+
+ if (!solved) {
+ if (pfound > 0 && distmax > 0) {
+
+ G_debug(1, "Weighted average");
+
+ result = interp_wa(in_seg, dval, cur_pnts, pfound,
+ row, col, src, dst, distmax,
+ &weight);
+
+ Segment_get(out_seg, (void *)&tps_out, irow, icol);
+
+ /* weight according to distance to nearest point */
+ 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, row, col);
+
+ cnt_wa++;
+ }
+ continue;
+ }
+
/* must be <= 0.5 */
mfactor = 0.0;
/* min: 0
@@ -855,50 +1070,54 @@
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) {
+ if (n_vars_i) {
Segment_get(var_seg, (void *)varbuf, irow, icol);
if (Rast_is_d_null_value(varbuf)) {
continue;
}
+
+ /* TODO: check current vars against range
+ * of observed vars
+ */
+ /*
+ use_tps = 0;
+ for (i = 0; i < n_vars_i; i++) {
+ if (varbuf[i] < min[i] - factor * (max[i] - min[i]) ||
+ varbuf[i] > max[i] + factor * (max[i] - min[i]))
+ use_tps = 1;
+ }
+ */
+
}
-#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) {
- for (j = 0; j < n_vars; j++) {
+ if (n_vars_i) {
+ for (j = 0; j < n_vars_i; 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) {
dist = dist2 * log(dist2) * 0.5;
- result += B[1 + n_vars + i] * dist;
+ result += B[1 + n_vars_i + i] * dist;
}
}
@@ -926,9 +1145,22 @@
}
G_percent(1, 1, 1);
- G_debug(1, "wmin: %g", wmin);
- G_debug(1, "wmax: %g", wmax);
+ G_debug(1, "min weight: %g", wmin);
+ G_debug(1, "max weight: %g", wmax);
+ G_debug(1, "Weighted average count: %u", cnt_wa);
+ if (n_vars > 0 && cnt_tps) {
+ double perc_no_vars;
+
+ G_verbose_message(_("Number of TPS interpolations with covariables: %u"),
+ cnt_tps_lm);
+ G_verbose_message(_("Number of TPS interpolations without covariables: %u"),
+ cnt_tps);
+ perc_no_vars = (double) 100.0 * cnt_tps / (cnt_tps + cnt_tps_lm);
+ G_verbose_message(_("Percentage of TPS interpolations without covariables: %.2f"),
+ perc_no_vars);
+ }
+
flag_destroy(pnt_flag);
outbuf = Rast_allocate_d_buf();
@@ -964,18 +1196,20 @@
int tps_window(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,
- double regularization, double overlap, int radius)
+ off_t n_points, double regularization, double overlap,
+ int radius, double lm_thresh)
{
int ridx, cidx, row, col, nrows, ncols, src_row, src_col;
double **m, *a, *B;
+ double **mvars, *avars, *Bvars;
int i, j;
int palloc;
struct tps_pnt *cur_pnts;
double dx, dy, dist, dist2, mfactor;
DCELL *dval, result, *outbuf, *varbuf;
CELL *maskbuf;
- int solved;
+ int solved, solved_tps_lm, solved_tps, solved_lm, n_vars_i;
+ double rsqr;
int pfound;
double distmax;
int mask_fd;
@@ -988,16 +1222,12 @@
int irow, irow1, irow2, icol, icol1, icol2;
int wsize;
unsigned int wacnt;
- double dxyw;
-#ifndef USE_RC
double i_n, i_e;
-#endif
nrows = Rast_window_rows();
ncols = Rast_window_cols();
wsize = (radius * 2 + 1) * (radius * 2 + 1);
- dxyw = (src->ew_res + src->ns_res) * 0.5 * (radius + 1);
palloc = wsize;
a = G_malloc((palloc + 1 + n_vars) * sizeof(double));
@@ -1009,6 +1239,9 @@
pvar = NULL;
varbuf = NULL;
+ mvars = NULL;
+ avars = NULL;
+ Bvars = NULL;
if (n_vars) {
pvar = G_malloc(palloc * n_vars * sizeof(double));
cur_pnts[0].vars = pvar;
@@ -1016,6 +1249,12 @@
cur_pnts[i].vars = cur_pnts[i - 1].vars + n_vars;
varbuf = G_malloc(n_vars * sizeof(DCELL));
+
+ avars = G_malloc((1 + n_vars) * sizeof(double));
+ Bvars = G_malloc((1 + n_vars) * sizeof(double));
+ mvars = G_malloc((1 + n_vars) * sizeof(double *));
+ for (i = 0; i < (1 + n_vars); i++)
+ mvars[i] = G_malloc((1 + n_vars) * sizeof(double));
}
dval = G_malloc((1 + n_vars) * sizeof(DCELL));
@@ -1145,60 +1384,73 @@
&rmin, &rmax, &cmin, &cmax,
&distmax);
+ /* sort points */
+ qsort(cur_pnts, pfound, sizeof(struct tps_pnt), cmp_pnts);
+
+ load_tps_pnts(in_seg, dval, n_vars, cur_pnts, pfound,
+ src, dst, regularization, m, a, mvars, avars);
+
+ n_vars_i = n_vars;
if (pfound > 2) {
- /* sort points */
- qsort(cur_pnts, pfound, sizeof(struct tps_pnt), cmp_pnts);
+ /* solve */
+ solved_tps_lm = solved_tps = 0;
- load_tps_pnts(in_seg, n_vars, cur_pnts, pfound, src, dst,
- regularization, m, a);
+ if (n_vars) {
- /* solve */
- solved = solvemat(m, a, B, pfound + 1 + n_vars);
- if (!solved) {
- for (i = 0; i < pfound; i++) {
- cur_pnts[i].r = (src->north - cur_pnts[i].r) / src->ns_res;
- cur_pnts[i].c = (cur_pnts[i].c - src->west) / src->ew_res;
+ solved_tps_lm = solvemat(m, a, B, pfound + 1 + n_vars);
+
+ if (solved_tps_lm && lm_thresh > 0) {
+
+ solved_lm = solvemat(mvars, avars, Bvars, 1 + n_vars);
+ if (!solved_lm) {
+ G_debug(1, "LM with covariables not working at row %d, col %d",
+ row, col);
+
+ solved_tps_lm = 0;
+ }
+ else {
+ rsqr = lm_rsqr(in_seg, n_vars, src, cur_pnts, pfound, B);
+
+ if (rsqr < lm_thresh) {
+ for (i = 1; i <= n_vars; i++) {
+ if (fabs(B[i]) > fabs(5 * Bvars[i])) {
+ G_debug(0, "LM B%d is %g but TPS B%d is %g",
+ i, Bvars[i], i, B[i]);
+ solved_tps_lm = 0;
+ }
+ }
+ }
+ }
}
+
+ if (!solved_tps_lm) {
+ n_vars_i = 0;
+ for (i = 0; i < pfound; i++) {
+ cur_pnts[i].r = (int)((src->north - cur_pnts[i].r) / src->ns_res);
+ cur_pnts[i].c = (int)((cur_pnts[i].c - src->west) / src->ew_res);
+ }
+ load_tps_pnts(in_seg, dval, 0, cur_pnts, pfound,
+ src, dst, regularization, m, a, NULL, NULL);
+ }
}
+
+ if (!solved_tps_lm) {
+ solved_tps = solvemat(m, a, B, pfound + 1);
+ }
+
+ solved = (solved_tps_lm | solved_tps);
}
if (!solved) {
if (pfound > 0) {
/* weighted average */
- double wsum, maxweight;
- i_n = dst->north - (row + 0.5) * dst->ns_res;
- i_e = dst->west + (col + 0.5) * dst->ew_res;
+ G_debug(1, "Weighted average");
+
+ result = interp_wa(in_seg, dval, cur_pnts, pfound,
+ row, col, src, dst, distmax,
+ &weight);
- result = 0;
- wsum = 0;
- maxweight = 0;
- for (i = 0; i < pfound; i++) {
-
- Segment_get(in_seg, (void *)dval, cur_pnts[i].r, cur_pnts[i].c);
-
- cur_pnts[i].r = src->north - (cur_pnts[i].r + 0.5) * src->ns_res;
- cur_pnts[i].c = src->west + (cur_pnts[i].c + 0.5) * src->ew_res;
-
- dx = cur_pnts[i].c - i_e;
- dy = cur_pnts[i].r - i_n;
-
- /* weight for average */
- dx = fabs(dx) / dxyw;
- dy = fabs(dx) / dxyw;
- dist2 = dx * dx + dy * dy;
- weight = exp(-dist2 * 4.0);
-
- /* weight for tps */
- if (maxweight < weight)
- maxweight = weight;
-
- wsum += weight;
- result += dval[0] * weight;
- }
- result /= wsum;
- weight = maxweight;
-
/* weight according to distance to nearest point */
if (tps_out.wmax < weight)
tps_out.wmax = weight;
@@ -1270,15 +1522,14 @@
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) {
+ if (n_vars_i) {
Segment_get(var_seg, (void *)varbuf, irow, icol);
if (Rast_is_d_null_value(varbuf)) {
@@ -1286,28 +1537,23 @@
}
}
-#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) {
+ if (n_vars_i) {
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;
Modified: grass-addons/grass7/raster/r.resamp.tps/tps.h
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/tps.h 2017-01-13 05:56:55 UTC (rev 70361)
+++ grass-addons/grass7/raster/r.resamp.tps/tps.h 2017-01-13 14:10:15 UTC (rev 70362)
@@ -15,11 +15,12 @@
int tps_nn(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);
+ off_t n_points, int min_points, int max_points,
+ double regularization, double overlap, int do_bfs,
+ double lm_thresh);
int tps_window(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,
- double regularization, double overlap, int radius);
+ off_t n_points, double regularization, double overlap,
+ int radius, double lm_thresh);
More information about the grass-commit
mailing list