[GRASS-SVN] r69916 - grass-addons/grass7/raster/r.resamp.tps
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Nov 25 14:19:42 PST 2016
Author: mmetz
Date: 2016-11-25 14:19:42 -0800 (Fri, 25 Nov 2016)
New Revision: 69916
Modified:
grass-addons/grass7/raster/r.resamp.tps/main.c
grass-addons/grass7/raster/r.resamp.tps/pavl.c
grass-addons/grass7/raster/r.resamp.tps/pavl.h
grass-addons/grass7/raster/r.resamp.tps/tps.c
grass-addons/grass7/raster/r.resamp.tps/tps.h
Log:
r.resamp.tps: +moving window method, optimized AVL tree
Modified: grass-addons/grass7/raster/r.resamp.tps/main.c
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/main.c 2016-11-25 12:21:53 UTC (rev 69915)
+++ grass-addons/grass7/raster/r.resamp.tps/main.c 2016-11-25 22:19:42 UTC (rev 69916)
@@ -35,13 +35,13 @@
struct GModule *module;
struct Option *in_opt, *ivar_opt, *ovar_opt, *out_opt, *minpnts_opt,
- *reg_opt, *ov_opt, *mask_opt, *mem_opt;
+ *radius_opt, *reg_opt, *ov_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;
+ int min_points, radius;
int r, c, nrows, ncols;
DCELL **dbuf, *dval;
@@ -92,6 +92,17 @@
_("Minimum number of points to use for TPS interpolation");
minpnts_opt->guisection = _("Settings");
+ radius_opt = G_define_option();
+ radius_opt->key = "radius";
+ radius_opt->type = TYPE_INTEGER;
+ radius_opt->required = NO;
+ radius_opt->answer = "0";
+ radius_opt->label =
+ _("Radius for moving window interpolation");
+ radius_opt->description =
+ _("If radius is > 0, moving window interpolation will be used instead of nearest neighbor search");
+ radius_opt->guisection = _("Settings");
+
ivar_opt = G_define_standard_option(G_OPT_R_INPUTS);
ivar_opt->key = "icovars";
ivar_opt->required = NO;
@@ -325,6 +336,10 @@
min_points);
}
+ radius = atoi(radius_opt->answer);
+ if (radius < 0)
+ radius = 0;
+
regularization = atof(reg_opt->answer);
if (regularization < 0)
regularization = 0;
@@ -335,12 +350,22 @@
if (overlap > 1)
overlap = 1;
- 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"));
+ 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) {
+ 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) {
+ G_fatal_error(_("TPS interpolation failed"));
+ }
+ }
Segment_close(&in_seg);
if (n_vars)
Modified: grass-addons/grass7/raster/r.resamp.tps/pavl.c
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/pavl.c 2016-11-25 12:21:53 UTC (rev 69915)
+++ grass-addons/grass7/raster/r.resamp.tps/pavl.c 2016-11-25 22:19:42 UTC (rev 69916)
@@ -20,8 +20,9 @@
02110-1301 USA.
*/
-/* from libavl-2.0.3
- some speed optimizations by Markus Metz
+/* Nov 2016, Markus Metz
+ * from libavl-2.0.3
+ * added safety checks and speed optimizations
*/
#include <assert.h>
@@ -33,7 +34,7 @@
with comparison function |compare| using parameter |param|
and memory allocator |allocator|.
Returns |NULL| if memory allocation failed. */
-struct pavl_table *pavl_create(pavl_comparison_func * compare, void *param,
+struct pavl_table *pavl_create(pavl_comparison_func * compare,
struct libavl_allocator *allocator)
{
struct pavl_table *tree;
@@ -43,13 +44,12 @@
if (allocator == NULL)
allocator = &pavl_allocator_default;
- tree = allocator->libavl_malloc(allocator, sizeof *tree);
+ tree = allocator->libavl_malloc(sizeof *tree);
if (tree == NULL)
return NULL;
tree->pavl_root = NULL;
tree->pavl_compare = compare;
- tree->pavl_param = param;
tree->pavl_alloc = allocator;
tree->pavl_count = 0;
@@ -66,7 +66,7 @@
p = tree->pavl_root;
while (p != NULL) {
- int cmp = tree->pavl_compare(item, p->pavl_data, tree->pavl_param);
+ int cmp = tree->pavl_compare(item, p->pavl_data);
if (cmp == 0)
return p->pavl_data;
@@ -95,7 +95,7 @@
q = NULL;
dir = 0;
while (p != NULL) {
- int cmp = tree->pavl_compare(item, p->pavl_data, tree->pavl_param);
+ int cmp = tree->pavl_compare(item, p->pavl_data);
if (cmp == 0)
return &p->pavl_data;
@@ -108,7 +108,7 @@
q = p, p = p->pavl_link[dir];
}
- n = tree->pavl_alloc->libavl_malloc(tree->pavl_alloc, sizeof *p);
+ n = tree->pavl_alloc->libavl_malloc(sizeof *p);
if (n == NULL)
return NULL;
@@ -117,13 +117,12 @@
n->pavl_parent = q;
n->pavl_data = item;
n->pavl_balance = 0;
- if (q != NULL)
- q->pavl_link[dir] = n;
- else {
+ if (q == NULL) {
tree->pavl_root = n;
return &n->pavl_data;
}
+ q->pavl_link[dir] = n;
p = n;
while (p != y) {
@@ -270,14 +269,14 @@
p = tree->pavl_root;
dir = 0;
- cmp = tree->pavl_compare(item, p->pavl_data, tree->pavl_param);
+ cmp = tree->pavl_compare(item, p->pavl_data);
while (cmp != 0) {
dir = cmp > 0;
p = p->pavl_link[dir];
if (p == NULL)
return NULL;
- cmp = tree->pavl_compare(item, p->pavl_data, tree->pavl_param);
+ cmp = tree->pavl_compare(item, p->pavl_data);
}
item = p->pavl_data;
@@ -326,7 +325,7 @@
dir = 0;
}
}
- tree->pavl_alloc->libavl_free(tree->pavl_alloc, p);
+ tree->pavl_alloc->libavl_free(p);
while (q != (struct pavl_node *)&tree->pavl_root) {
struct pavl_node *y = q;
@@ -501,13 +500,14 @@
void *item)
{
struct pavl_node *p;
- int dir;
assert(trav != NULL && tree != NULL && item != NULL);
trav->pavl_table = tree;
- for (p = tree->pavl_root; p != NULL; p = p->pavl_link[dir]) {
- int cmp = tree->pavl_compare(item, p->pavl_data, tree->pavl_param);
+
+ p = tree->pavl_root;
+ while (p != NULL) {
+ int cmp = tree->pavl_compare(item, p->pavl_data);
if (cmp == 0) {
trav->pavl_node = p;
@@ -515,10 +515,11 @@
return p->pavl_data;
}
- dir = cmp > 0;
+ p = p->pavl_link[cmp > 0];
}
trav->pavl_node = NULL;
+
return NULL;
}
@@ -539,8 +540,9 @@
p = pavl_probe(tree, item);
if (p != NULL) {
trav->pavl_table = tree;
- trav->pavl_node = ((struct pavl_node *) ((char *)p -
- offsetof(struct pavl_node, pavl_data)));
+ trav->pavl_node = ((struct pavl_node *)((char *)p -
+ offsetof(struct pavl_node,
+ pavl_data)));
return *p;
}
@@ -580,8 +582,8 @@
if (q == NULL || p == q->pavl_link[0]) {
trav->pavl_node = q;
- return trav->pavl_node != NULL ?
- trav->pavl_node->pavl_data : NULL;
+ return trav->pavl_node != NULL ?
+ trav->pavl_node->pavl_data : NULL;
}
}
else {
@@ -610,8 +612,8 @@
if (q == NULL || p == q->pavl_link[1]) {
trav->pavl_node = q;
- return trav->pavl_node != NULL ?
- trav->pavl_node->pavl_data : NULL;
+ return trav->pavl_node != NULL ?
+ trav->pavl_node->pavl_data : NULL;
}
}
else {
@@ -686,7 +688,7 @@
struct pavl_node *y;
assert(org != NULL);
- new = pavl_create(org->pavl_compare, org->pavl_param,
+ new = pavl_create(org->pavl_compare,
allocator != NULL ? allocator : org->pavl_alloc);
if (new == NULL)
return NULL;
@@ -697,11 +699,10 @@
x = (const struct pavl_node *)&org->pavl_root;
y = (struct pavl_node *)&new->pavl_root;
- for (;;) {
+ while (x != NULL) {
while (x->pavl_link[0] != NULL) {
y->pavl_link[0] =
- new->pavl_alloc->libavl_malloc(new->pavl_alloc,
- sizeof *y->pavl_link[0]);
+ new->pavl_alloc->libavl_malloc(sizeof *y->pavl_link[0]);
if (y->pavl_link[0] == NULL) {
if (y != (struct pavl_node *)&new->pavl_root) {
y->pavl_data = NULL;
@@ -724,7 +725,7 @@
if (copy == NULL)
y->pavl_data = x->pavl_data;
else {
- y->pavl_data = copy(x->pavl_data, org->pavl_param);
+ y->pavl_data = copy(x->pavl_data);
if (y->pavl_data == NULL) {
y->pavl_link[1] = NULL;
copy_error_recovery(y, new, destroy);
@@ -735,8 +736,7 @@
if (x->pavl_link[1] != NULL) {
y->pavl_link[1] =
- new->pavl_alloc->libavl_malloc(new->pavl_alloc,
- sizeof *y->pavl_link[1]);
+ new->pavl_alloc->libavl_malloc(sizeof *y->pavl_link[1]);
if (y->pavl_link[1] == NULL) {
copy_error_recovery(y, new, destroy);
@@ -783,8 +783,8 @@
if (p->pavl_link[0] == NULL) {
q = p->pavl_link[1];
if (destroy != NULL && p->pavl_data != NULL)
- destroy(p->pavl_data, tree->pavl_param);
- tree->pavl_alloc->libavl_free(tree->pavl_alloc, p);
+ destroy(p->pavl_data);
+ tree->pavl_alloc->libavl_free(p);
}
else {
q = p->pavl_link[0];
@@ -794,24 +794,24 @@
p = q;
}
- tree->pavl_alloc->libavl_free(tree->pavl_alloc, tree);
+ tree->pavl_alloc->libavl_free(tree);
}
/* Allocates |size| bytes of space using |malloc()|.
Returns a null pointer if allocation fails. */
-void *pavl_malloc(struct libavl_allocator *allocator, size_t size)
+void *pavl_malloc(size_t size)
{
- assert(allocator != NULL && size > 0);
+ if (size > 0)
+ return malloc(size);
- return malloc(size);
+ return NULL;
}
/* Frees |block|. */
-void pavl_free(struct libavl_allocator *allocator, void *block)
+void pavl_free(void *block)
{
- assert(allocator != NULL && block != NULL);
-
- free(block);
+ if (block)
+ free(block);
}
/* Default memory allocator that uses |malloc()| and |free()|. */
Modified: grass-addons/grass7/raster/r.resamp.tps/pavl.h
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/pavl.h 2016-11-25 12:21:53 UTC (rev 69915)
+++ grass-addons/grass7/raster/r.resamp.tps/pavl.h 2016-11-25 22:19:42 UTC (rev 69916)
@@ -18,7 +18,7 @@
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301 USA.
-*/
+ */
#ifndef PAVL_H
#define PAVL_H 1
@@ -26,25 +26,24 @@
#include <stddef.h>
/* Function types. */
-typedef int pavl_comparison_func (const void *pavl_a, const void *pavl_b,
- void *pavl_param);
-typedef void pavl_item_func (void *pavl_item, void *pavl_param);
-typedef void *pavl_copy_func (void *pavl_item, void *pavl_param);
+typedef int pavl_comparison_func(const void *pavl_a, const void *pavl_b);
+typedef void pavl_item_func(void *pavl_item);
+typedef void *pavl_copy_func(void *pavl_item);
#ifndef LIBAVL_ALLOCATOR
#define LIBAVL_ALLOCATOR
/* Memory allocator. */
struct libavl_allocator
{
- void *(*libavl_malloc) (struct libavl_allocator *, size_t libavl_size);
- void (*libavl_free) (struct libavl_allocator *, void *libavl_block);
+ void *(*libavl_malloc) (size_t libavl_size);
+ void (*libavl_free) (void *libavl_block);
};
#endif
/* Default memory allocator. */
extern struct libavl_allocator pavl_allocator_default;
-void *pavl_malloc (struct libavl_allocator *, size_t);
-void pavl_free (struct libavl_allocator *, void *);
+void *pavl_malloc(size_t);
+void pavl_free(void *);
/* Maximum PAVL height, unused. */
#ifndef PAVL_MAX_HEIGHT
@@ -54,55 +53,54 @@
/* Tree data structure. */
struct pavl_table
{
- struct pavl_node *pavl_root; /* Tree's root. */
- pavl_comparison_func *pavl_compare; /* Comparison function. */
- void *pavl_param; /* Extra argument to |pavl_compare|. */
- struct libavl_allocator *pavl_alloc; /* Memory allocator. */
- size_t pavl_count; /* Number of items in tree. */
+ struct pavl_node *pavl_root; /* Tree's root. */
+ pavl_comparison_func *pavl_compare; /* Comparison function. */
+ struct libavl_allocator *pavl_alloc; /* Memory allocator. */
+ size_t pavl_count; /* Number of items in tree. */
};
/* An PAVL tree node. */
struct pavl_node
{
- struct pavl_node *pavl_link[2]; /* Subtrees. */
- struct pavl_node *pavl_parent; /* Parent node. */
- void *pavl_data; /* Pointer to data. */
- signed char pavl_balance; /* Balance factor. */
+ struct pavl_node *pavl_link[2]; /* Subtrees. */
+ struct pavl_node *pavl_parent; /* Parent node. */
+ void *pavl_data; /* Pointer to data. */
+ signed char pavl_balance; /* Balance factor. */
};
/* PAVL traverser structure. */
struct pavl_traverser
{
- struct pavl_table *pavl_table; /* Tree being traversed. */
- struct pavl_node *pavl_node; /* Current node in tree. */
+ struct pavl_table *pavl_table; /* Tree being traversed. */
+ struct pavl_node *pavl_node; /* Current node in tree. */
};
/* Table functions. */
-struct pavl_table *pavl_create (pavl_comparison_func *, void *,
- struct libavl_allocator *);
-struct pavl_table *pavl_copy (const struct pavl_table *, pavl_copy_func *,
- pavl_item_func *, struct libavl_allocator *);
-void pavl_destroy (struct pavl_table *, pavl_item_func *);
-void **pavl_probe (struct pavl_table *, void *);
-void *pavl_insert (struct pavl_table *, void *);
-void *pavl_replace (struct pavl_table *, void *);
-void *pavl_delete (struct pavl_table *, const void *);
-void *pavl_find (const struct pavl_table *, const void *);
-void pavl_assert_insert (struct pavl_table *, void *);
-void *pavl_assert_delete (struct pavl_table *, void *);
+struct pavl_table *pavl_create(pavl_comparison_func *,
+ struct libavl_allocator *);
+struct pavl_table *pavl_copy(const struct pavl_table *, pavl_copy_func *,
+ pavl_item_func *, struct libavl_allocator *);
+void pavl_destroy(struct pavl_table *, pavl_item_func *);
+void **pavl_probe(struct pavl_table *, void *);
+void *pavl_insert(struct pavl_table *, void *);
+void *pavl_replace(struct pavl_table *, void *);
+void *pavl_delete(struct pavl_table *, const void *);
+void *pavl_find(const struct pavl_table *, const void *);
+void pavl_assert_insert(struct pavl_table *, void *);
+void *pavl_assert_delete(struct pavl_table *, void *);
#define pavl_count(table) ((size_t) (table)->pavl_count)
/* Table traverser functions. */
-void pavl_t_init (struct pavl_traverser *, struct pavl_table *);
-void *pavl_t_first (struct pavl_traverser *, struct pavl_table *);
-void *pavl_t_last (struct pavl_traverser *, struct pavl_table *);
-void *pavl_t_find (struct pavl_traverser *, struct pavl_table *, void *);
-void *pavl_t_insert (struct pavl_traverser *, struct pavl_table *, void *);
-void *pavl_t_copy (struct pavl_traverser *, const struct pavl_traverser *);
-void *pavl_t_next (struct pavl_traverser *);
-void *pavl_t_prev (struct pavl_traverser *);
-void *pavl_t_cur (struct pavl_traverser *);
-void *pavl_t_replace (struct pavl_traverser *, void *);
+void pavl_t_init(struct pavl_traverser *, struct pavl_table *);
+void *pavl_t_first(struct pavl_traverser *, struct pavl_table *);
+void *pavl_t_last(struct pavl_traverser *, struct pavl_table *);
+void *pavl_t_find(struct pavl_traverser *, struct pavl_table *, void *);
+void *pavl_t_insert(struct pavl_traverser *, struct pavl_table *, void *);
+void *pavl_t_copy(struct pavl_traverser *, const struct pavl_traverser *);
+void *pavl_t_next(struct pavl_traverser *);
+void *pavl_t_prev(struct pavl_traverser *);
+void *pavl_t_cur(struct pavl_traverser *);
+void *pavl_t_replace(struct pavl_traverser *, void *);
#endif /* pavl.h */
Modified: grass-addons/grass7/raster/r.resamp.tps/tps.c
===================================================================
--- grass-addons/grass7/raster/r.resamp.tps/tps.c 2016-11-25 12:21:53 UTC (rev 69915)
+++ grass-addons/grass7/raster/r.resamp.tps/tps.c 2016-11-25 22:19:42 UTC (rev 69916)
@@ -193,7 +193,7 @@
}
-static int cmp_rc(const void *first, const void *second, void *avl_param)
+static int cmp_rc(const void *first, const void *second)
{
struct rc *a = (struct rc *)first, *b = (struct rc *)second;
@@ -203,7 +203,7 @@
return (a->row - b->row);
}
-static void avl_free_item(void *avl_item, void *avl_param)
+static void avl_free_item(void *avl_item)
{
G_free(avl_item);
}
@@ -223,7 +223,7 @@
struct pavl_table *visited;
double dx, dy, dist;
- visited = pavl_create(cmp_rc, NULL, NULL);
+ visited = pavl_create(cmp_rc, NULL);
ngbr_rc.row = row;
ngbr_rc.col = col;
@@ -338,7 +338,7 @@
struct pavl_table *visited;
double dx, dy, dist;
- visited = pavl_create(cmp_rc, NULL, NULL);
+ visited = pavl_create(cmp_rc, NULL);
ngbr_rc.row = row;
ngbr_rc.col = col;
@@ -411,12 +411,70 @@
return found;
}
-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)
+static int window_pnts(FLAG *pnt_flag, struct Cell_head *src,
+ struct tps_pnt *cur_pnts, int radius,
+ int row, int col,
+ int *rmin, int *rmax, int *cmin, int *cmax,
+ double *distmax)
{
+ int rown, coln, row1, row2, col1, col2;
+ int found;
+ double dx, dy, dist;
+
+ found = 0;
+ *distmax = 0.0;
+
+ row1 = row - radius;
+ if (row1 < 0)
+ row1 = 0;
+ row2 = row + radius;
+ if (row2 > src->rows - 1)
+ row2 = src->rows - 1;
+
+ col1 = col - radius;
+ if (col1 < 0)
+ col1 = 0;
+ col2 = col + radius;
+ if (col2 > src->cols - 1)
+ col2 = src->cols - 1;
+
+ for (rown = row1; rown <= row2; rown++) {
+ for (coln = col1; coln <= col2; coln++) {
+ if (FLAG_GET(pnt_flag, rown, coln)) {
+ /* add to cur_pnts */
+ 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;
+
+ dx = coln - col;
+ dy = rown - row;
+
+ dist = dx * dx + dy * dy;
+
+ if (*distmax < dist)
+ *distmax = dist;
+ }
+ }
+ }
+
+ return found;
+}
+
+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)
+{
int ridx, cidx, row, col, nrows, ncols, src_row, src_col;
double **m, *a, *B;
int i, j;
@@ -530,7 +588,7 @@
}
G_percent(1, 1, 2);
- G_message(_("Local TPS interpolation with %ld points..."), n_points);
+ G_message(_("Nearest neighbor TPS interpolation with %ld points..."), n_points);
wmin = 10;
wmax = 0;
@@ -868,8 +926,8 @@
}
G_percent(1, 1, 1);
- G_debug(0, "wmin: %g", wmin);
- G_debug(0, "wmax: %g", wmax);
+ G_debug(1, "wmin: %g", wmin);
+ G_debug(1, "wmax: %g", wmax);
flag_destroy(pnt_flag);
@@ -902,3 +960,415 @@
return 1;
}
+
+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)
+{
+ int ridx, cidx, row, col, nrows, ncols, src_row, src_col;
+ double **m, *a, *B;
+ 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 pfound;
+ double distmax, mindist;
+ 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 wsize;
+ 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));
+ 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 * sizeof(struct tps_pnt));
+
+ pvar = NULL;
+ varbuf = NULL;
+ if (n_vars) {
+ pvar = G_malloc(palloc * n_vars * sizeof(double));
+ cur_pnts[0].vars = pvar;
+ for (i = 1; i < palloc; i++)
+ cur_pnts[i].vars = cur_pnts[i - 1].vars + n_vars;
+
+ 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_name) {
+ G_message("Loading mask map...");
+ mask_fd = Rast_open_old(mask_name, "");
+ 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);
+ }
+ }
+ Rast_close(mask_fd);
+ G_free(maskbuf);
+ maskbuf = NULL;
+ }
+
+ G_message(_("Analyzing input ..."));
+ pnt_flag = flag_create(src->rows, src->cols);
+
+ rminp = src->rows;
+ rmaxp = 0;
+ cminp = src->cols;
+ cmaxp = 0;
+ for (row = 0; row < src->rows; row++) {
+ G_percent(row, src->rows, 2);
+
+ 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;
+ }
+ }
+ }
+ 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;
+ tps_out.wsum = 0;
+ tps_out.wmax = 0;
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+
+ for (col = 0; col < ncols; col++) {
+ if (Segment_put(out_seg, (void *)&tps_out, row, col) != 1)
+ G_fatal_error(_("Unable to write to temporary file"));
+ }
+ }
+ G_percent(1, 1, 2);
+
+ G_message(_("Moving window TPS interpolation with %ld points..."), n_points);
+
+ wmin = 10;
+ wmax = 0;
+
+ if (overlap > 1.0)
+ overlap = 1.0;
+ if (overlap < 0.0)
+ overlap = 0.0;
+ /* keep in sync with weight calculation below */
+ overlap = exp((overlap - 1.0) * 8.0);
+
+ for (ridx = 0; ridx < nrows; ridx++) {
+
+ G_percent(ridx, nrows, 1);
+
+ row = (ridx >> 1);
+ if (ridx & 1)
+ row = nrows - 1 - row;
+
+ src_row = row_src2dst(row, dst, src);
+
+ for (cidx = 0; cidx < ncols; cidx++) {
+
+ col = (cidx >> 1);
+ if (cidx & 1)
+ col = ncols - 1 - col;
+
+ if ((FLAG_GET(mask_flag, row, col))) {
+ continue;
+ }
+
+ if (n_vars) {
+
+ 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);
+ if (tps_out.wmax > overlap)
+ continue;
+
+ solved = 0;
+ pfound = 0;
+ src_col = col_src2dst(col, dst, src);
+
+ /* collect points within moving window */
+ rmin = src->rows;
+ rmax = 0;
+ cmin = src->cols;
+ cmax = 0;
+ pfound = window_pnts(pnt_flag, src, cur_pnts,
+ radius,
+ src_row, src_col,
+ &rmin, &rmax, &cmin, &cmax,
+ &distmax);
+
+ if (pfound > 2) {
+ /* sort points */
+ 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);
+ 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;
+ }
+ }
+ }
+
+ 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;
+
+ mindist = -1;
+ 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 (mindist > dist2 || mindist > 0)
+ 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;
+
+ tps_out.val += result * weight;
+ tps_out.wsum += weight;
+ Segment_put(out_seg, (void *)&tps_out, row, col);
+ }
+ continue;
+ }
+
+ /* must be <= 0.5 */
+ mfactor = 0.0;
+ /* min: 0
+ * max: 0.5
+ * 1 - 1 / d: -> 0 for dense spacing
+ * 1 for sparse points
+ */
+
+ 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);
+ icol2 = cmax - (int)((cmax - cmin) * mfactor);
+
+ 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;
+ }
+
+ if (rmin == rminp)
+ irow1 = 0;
+ if (rmax == rmaxp)
+ irow2 = nrows - 1;
+ if (cmin == cminp)
+ icol1 = 0;
+ if (cmax == cmaxp)
+ icol2 = ncols - 1;
+
+ if (irow1 < 0)
+ irow1 = 0;
+ if (irow2 > nrows - 1)
+ irow2 = nrows - 1;
+ if (icol1 < 0)
+ icol1 = 0;
+ if (icol2 > ncols - 1)
+ icol2 = ncols - 1;
+
+ dxi = icol2 - icol1 + 1;
+ 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) {
+
+ Segment_get(var_seg, (void *)varbuf, irow, icol);
+ if (Rast_is_d_null_value(varbuf)) {
+ continue;
+ }
+ }
+
+#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++) {
+ 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;
+ }
+ }
+
+ dx = fabs(2.0 * icol - (icol2 + icol1)) / dxi;
+ dy = fabs(2.0 * irow - (irow2 + irow1)) / dyi;
+
+ dist2 = (dx * dx + dy * dy);
+
+ weight = exp(-dist2 * 4.0);
+
+ 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);
+
+ flag_destroy(pnt_flag);
+
+ 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);
+ 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-25 12:21:53 UTC (rev 69915)
+++ grass-addons/grass7/raster/r.resamp.tps/tps.h 2016-11-25 22:19:42 UTC (rev 69916)
@@ -12,8 +12,14 @@
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);
+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);
+
+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);
More information about the grass-commit
mailing list