r.cost,r.walk: avoid circular paths when using solver

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@73001 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.cost/flag.c

@@ -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;
+}

raster/r.cost/flag.h

@@ -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__ */

raster/r.cost/main.c

@@ -11,6 +11,8 @@
  *               Updated for calculation errors and directional surface generation
  *                 Colin Nielsen <colin.nielsen gmail com>
  *               Use min heap instead of btree (faster, less memory)
+ *               multiple directions with bitmask encoding
+ *               avoid circular paths
  *                 Markus Metz
  *
  * PURPOSE:      Outputs a raster map layer showing the cumulative cost 
@@ -74,6 +76,7 @@
 #include <grass/glocale.h>
 #include "cost.h"
 #include "stash.h"
+#include "flag.h"
 
 #define SEGCOLSIZE 	64
 
@@ -140,6 +143,7 @@ int main(int argc, char *argv[])
     struct cc {
 	double cost_in, cost_out, nearest;
     } costs;
+    FLAG *visited;
 
     void *ptr2;
     RASTER_MAP_TYPE data_type,			 	/* input cost type */
@@ -806,6 +810,7 @@ int main(int argc, char *argv[])
     G_debug(1, "nrows x ncols: %d", nrows * ncols);
     G_message(_("Finding cost path..."));
     n_processed = 0;
+    visited = flag_create(nrows, ncols);
 
     pres_cell = get_lowest();
     while (pres_cell != NULL) {
@@ -830,6 +835,12 @@ int main(int argc, char *argv[])
 		continue;
 	    }
 	}
+	if (FLAG_GET(visited, pres_cell->row, pres_cell->col)) {
+	    delete(pres_cell);
+	    pres_cell = get_lowest();
+	    continue;
+	}
+	FLAG_SET(visited, pres_cell->row, pres_cell->col);
 
 	if (have_solver)
 	    Segment_get(&solve_seg, mysolvedir, pres_cell->row, pres_cell->col);
@@ -1006,6 +1017,8 @@ int main(int argc, char *argv[])
 	    if (col < 0 || col >= ncols)
 		continue;
 
+	    /* skip already processed neighbors here ? */
+
 	    min_cost = dnullval;
 	    Segment_get(&cost_seg, &costs, row, col);
 
@@ -1133,12 +1146,17 @@ int main(int argc, char *argv[])
 		    Segment_put(&solve_seg, solvedir, row, col);
 		}
 	    }
-	    else if (old_min_cost == min_cost && (dir_bin || have_solver)) {
+	    else if (old_min_cost == min_cost &&
+	             (dir_bin || have_solver) && 
+		     !(FLAG_GET(visited, row, col))) {
 		FCELL old_dir;
 		int dir_inv[16] = { 4, 5, 6, 7, 0, 1, 2, 3,
 		                   12, 13, 14, 15, 8, 9, 10, 11 };
 		int dir_fwd;
 		int equal = 1;
+		
+		/* only update neighbors that have not yet been processed,
+		 * otherwise we might get circular paths */
 
 		if (have_solver) {
 		    Segment_get(&solve_seg, solvedir, row, col);
@@ -1187,6 +1205,7 @@ int main(int argc, char *argv[])
 
     /* free heap */
     free_heap();
+    flag_destroy(visited);
 
     if (have_solver) {
 	Segment_close(&solve_seg);

raster/r.walk/flag.c

@@ -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;
+}

raster/r.walk/flag.h

@@ -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__ */

raster/r.walk/main.c

@@ -28,6 +28,8 @@
  *               Updated for calculation errors and directional surface generation
  *                 Colin Nielsen <colin.nielsen gmail com>
  *               Use min heap instead of btree (faster, less memory)
+ *               multiple directions with bitmask encoding
+ *               avoid circular paths
  *                 Markus Metz
  * PURPOSE:      anisotropic movements on cost surfaces
  * COPYRIGHT:    (C) 1999-2015 by the GRASS Development Team
@@ -108,6 +110,7 @@
 #include <grass/glocale.h>
 #include "cost.h"
 #include "stash.h"
+#include "flag.h"
 
 #define SEGCOLSIZE 	64
 
@@ -179,6 +182,7 @@ int main(int argc, char *argv[])
 	double cost_in;		/* friction costs */
 	double cost_out;	/* cumulative costs */
     } costs;
+    FLAG *visited;
 
     void *ptr1, *ptr2;
     RASTER_MAP_TYPE dtm_data_type, cost_data_type, cum_data_type =
@@ -941,6 +945,7 @@ int main(int argc, char *argv[])
     G_debug(1, "nrows x ncols: %d", nrows * ncols);
     G_message(_("Finding cost path..."));
     n_processed = 0;
+    visited = flag_create(nrows, ncols);
 
     pres_cell = get_lowest();
     while (pres_cell != NULL) {
@@ -973,7 +978,6 @@ int main(int argc, char *argv[])
 		continue;
 	    }
 	}
-
 	my_dtm = costs.dtm;
 	if (Rast_is_d_null_value(&my_dtm)) {
 	    delete(pres_cell);
@@ -986,6 +990,12 @@ int main(int argc, char *argv[])
 	    pres_cell = get_lowest();
 	    continue;
 	}
+	if (FLAG_GET(visited, pres_cell->row, pres_cell->col)) {
+	    delete(pres_cell);
+	    pres_cell = get_lowest();
+	    continue;
+	}
+	FLAG_SET(visited, pres_cell->row, pres_cell->col);
 
 	if (have_solver)
 	    Segment_get(&solve_seg, mysolvedir, pres_cell->row, pres_cell->col);
@@ -1159,6 +1169,8 @@ int main(int argc, char *argv[])
 	    if (col < 0 || col >= ncols)
 		continue;
 
+	    /* skip already processed neighbors here ? */
+
 	    min_cost = dnullval;
 	    Segment_get(&cost_seg, &costs, row, col);
 
@@ -1484,12 +1496,17 @@ int main(int argc, char *argv[])
 		    Segment_put(&solve_seg, solvedir, row, col);
 		}
 	    }
-	    else if (dir && dir_bin && old_min_cost == min_cost) {
+	    else if (old_min_cost == min_cost &&
+	             (dir_bin || have_solver) && 
+		     !(FLAG_GET(visited, row, col))) {
 		FCELL old_dir;
 		int dir_inv[16] = { 4, 5, 6, 7, 0, 1, 2, 3,
 		                   12, 13, 14, 15, 8, 9, 10, 11 };
 		int dir_fwd;
 		int equal = 1;
+		
+		/* only update neighbors that have not yet been processed,
+		 * otherwise we might get circular paths */
 
 		if (have_solver) {
 		    Segment_get(&solve_seg, solvedir, row, col);
@@ -1508,7 +1525,7 @@ int main(int argc, char *argv[])
 		    }
 		}
 
-		if (equal) {
+		if (dir_bin && equal) {
 		    /* this can create circular paths:
 		     * set only if current cell does not point to neighbor
 		     * does not avoid longer circular paths */
@@ -1537,6 +1554,7 @@ int main(int argc, char *argv[])
 
     /* free heap */
     free_heap();
+    flag_destroy(visited);
 
     if (have_solver) {
 	Segment_close(&solve_seg);