[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