r.cost/r.walk: add option to control which direction is used in case of multiple directions with equal costs

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

raster/r.cost/main.c

@@ -105,7 +105,8 @@ int main(int argc, char *argv[])
     const char *cost_layer;
     const char *cost_mapset, *search_mapset;
     void *cell, *cell2, *dir_cell, *nearest_cell;
-    SEGMENT cost_seg, dir_seg;
+    SEGMENT cost_seg, dir_seg, solve_seg;
+    int have_solver;
     double *value;
     extern struct Cell_head window;
     double NS_fac, EW_fac, DIAG_fac, H_DIAG_fac, V_DIAG_fac;
@@ -115,7 +116,7 @@ int main(int argc, char *argv[])
     double zero = 0.0;
     int col, row, nrows, ncols;
     int maxcost;
-    int nseg;
+    int nseg, nbytes;
     int maxmem;
     int segments_in_memory;
     int cost_fd, cum_fd, dir_fd, nearest_fd;
@@ -132,7 +133,7 @@ int main(int argc, char *argv[])
     struct GModule *module;
     struct Flag *flag2, *flag3, *flag4, *flag5, *flag6;
     struct Option *opt1, *opt2, *opt3, *opt4, *opt5, *opt6, *opt7, *opt8;
-    struct Option *opt9, *opt10, *opt11, *opt12;
+    struct Option *opt9, *opt10, *opt11, *opt12, *opt_solve;
     struct cost *pres_cell;
     struct start_pt *head_start_pt = NULL;
     struct start_pt *next_start_pt;
@@ -150,6 +151,7 @@ int main(int argc, char *argv[])
     int dsize, nearest_size;
     double disk_mb, mem_mb, pq_mb;
     int dir_bin;
+    DCELL mysolvedir[2], solvedir[2];
 
     G_gisinit(argv[0]);
 
@@ -164,11 +166,19 @@ int main(int argc, char *argv[])
 	  "geographic locations on an input raster map "
 	  "whose cell category values represent cost.");
 
+    opt1 = G_define_standard_option(G_OPT_R_OUTPUT);
+
     opt2 = G_define_standard_option(G_OPT_R_INPUT);
     opt2->description =
 	_("Name of input raster map containing grid cell cost information");
 
-    opt1 = G_define_standard_option(G_OPT_R_OUTPUT);
+    opt_solve = G_define_standard_option(G_OPT_R_INPUT);
+    opt_solve->key = "solver";
+    opt_solve->required = NO;
+    opt_solve->label =
+	_("Name of input raster map solving equal costs");
+    opt_solve->description =
+	_("Helper variable to pick a direction if two directions have equal cumulative costs (smaller is better)");
 
     opt12 = G_define_standard_option(G_OPT_R_OUTPUT);
     opt12->key = "nearest";
@@ -305,7 +315,9 @@ int main(int argc, char *argv[])
 
     start_with_raster_vals = flag4->answer;
 
-    dir_bin = flag6->answer;
+    dir_bin = 0;
+    if (dir)
+	dir_bin = flag6->answer;
 
     {
 	int count = 0;
@@ -352,6 +364,14 @@ int main(int argc, char *argv[])
 	    G_fatal_error(_("Vector map <%s> not found"), opt7->answer);
     }
 
+    have_solver = 0;
+    if (dir && opt_solve->answer) {
+	search_mapset = G_find_raster2(opt_solve->answer, "");
+	if (search_mapset == NULL)
+	    G_fatal_error(_("Raster map <%s> not found"), opt_solve->answer);
+	have_solver = 1;
+    }
+
     if (!Rast_is_d_null_value(&null_cost)) {
 	if (null_cost < 0.0) {
 	    G_warning(_("Assigning negative cost to null cell. Null cells excluded."));
@@ -414,26 +434,21 @@ int main(int argc, char *argv[])
     maxmem -= pq_mb;
     if (maxmem < 10)
 	maxmem = 10;
-    if (dir == TRUE) {
-	disk_mb = (double) nrows * ncols * 28. / 1048576.;
-	segments_in_memory = maxmem / 
-	                     ((double) srows * scols * (28. / 1048576.));
-	if (segments_in_memory < 4)
-	    segments_in_memory = 4;
-	if (segments_in_memory > nseg)
-	    segments_in_memory = nseg;
-	mem_mb = (double) srows * scols * (28. / 1048576.) * segments_in_memory;
-    }
-    else {
-	disk_mb = (double) nrows * ncols * 24. / 1048576.;
-	segments_in_memory = maxmem / 
-	                     ((double) srows * scols * (24. / 1048576.));
-	if (segments_in_memory < 4)
-	    segments_in_memory = 4;
-	if (segments_in_memory > nseg)
-	    segments_in_memory = nseg;
-	mem_mb = (double) srows * scols * (24. / 1048576.) * segments_in_memory;
-    }
+
+    nbytes = 24;
+    if (dir == TRUE)
+	nbytes += 4;
+    if (have_solver)
+	nbytes += 16;
+
+    disk_mb = (double) nrows * ncols * nbytes / 1048576.;
+    segments_in_memory = maxmem / 
+			 ((double) srows * scols * (nbytes / 1048576.));
+    if (segments_in_memory < 4)
+	segments_in_memory = 4;
+    if (segments_in_memory > nseg)
+	segments_in_memory = nseg;
+    mem_mb = (double) srows * scols * (nbytes / 1048576.) * segments_in_memory;
 
     if (flag5->answer) {
 	fprintf(stdout, _("Will need at least %.2f MB of disk space"), disk_mb);
@@ -467,6 +482,31 @@ int main(int argc, char *argv[])
 	    G_fatal_error(_("Can not create temporary file"));
     }
 
+    if (have_solver) {
+	int sfd;
+
+	if (Segment_open(&solve_seg, G_tempfile(), nrows, ncols, srows, scols,
+		         sizeof(DCELL) * 2, segments_in_memory) != 1)
+	    G_fatal_error(_("Can not create temporary file"));
+
+	sfd = Rast_open_old(opt_solve->answer, "");
+	cell = Rast_allocate_buf(DCELL_TYPE);
+	Rast_set_d_null_value(&solvedir[1], 1); 
+	dsize = Rast_cell_size(DCELL_TYPE);
+	for (row = 0; row < nrows; row++) {
+	    G_percent(row, nrows, 2);
+	    Rast_get_d_row(sfd, cell, row);
+	    ptr2 = cell;
+	    for (col = 0; col < ncols; col++) {
+		solvedir[0] = *(DCELL *)ptr2;
+		Segment_put(&solve_seg, solvedir, row, col);
+		ptr2 = G_incr_void_ptr(ptr2, dsize);
+	    }
+	}
+	Rast_close(sfd);
+	G_free(cell);
+    }
+
     /* Write the cost layer in the segmented file */
     G_message(_("Reading raster map <%s>, initializing output..."),
 	      G_fully_qualified_name(cost_layer, cost_mapset));
@@ -791,6 +831,9 @@ int main(int argc, char *argv[])
 	    }
 	}
 
+	if (have_solver)
+	    Segment_get(&solve_seg, mysolvedir, pres_cell->row, pres_cell->col);
+
 	my_cost = costs.cost_in;
 	nearest = costs.nearest;
 
@@ -1067,6 +1110,11 @@ int main(int argc, char *argv[])
 			cur_dir = (1 << (int)cur_dir);
 		    Segment_put(&dir_seg, &cur_dir, row, col);
 		}
+		if (have_solver) {
+		    Segment_get(&solve_seg, solvedir, row, col);
+		    solvedir[1] = mysolvedir[0];
+		    Segment_put(&solve_seg, solvedir, row, col);
+		}
 	    }
 	    /* update with lower costs */
 	    else if (old_min_cost > min_cost) {
@@ -1079,22 +1127,48 @@ int main(int argc, char *argv[])
 			cur_dir = (1 << (int)cur_dir);
 		    Segment_put(&dir_seg, &cur_dir, row, col);
 		}
+		if (have_solver) {
+		    Segment_get(&solve_seg, solvedir, row, col);
+		    solvedir[1] = mysolvedir[0];
+		    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)) {
 		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;
 
-		/* this can create circular paths:
-		 * set only if current cell does not point to neighbor
-                 * does not avoid longer circular paths */
-		Segment_get(&dir_seg, &old_dir, pres_cell->row, pres_cell->col);
-		dir_fwd = (1 << dir_inv[(int)cur_dir]);
-		if (!((int)old_dir & dir_fwd)) {
-		    Segment_get(&dir_seg, &old_dir, row, col);
-		    cur_dir = ((1 << (int)cur_dir) | (int)old_dir);
-		    Segment_put(&dir_seg, &cur_dir, row, col);
+		if (have_solver) {
+		    Segment_get(&solve_seg, solvedir, row, col);
+		    equal = (solvedir[1] == mysolvedir[0]);
+		    if (solvedir[1] > mysolvedir[0]) {
+			solvedir[1] = mysolvedir[0];
+			Segment_put(&solve_seg, solvedir, row, col);
+
+			costs.nearest = nearest;
+			Segment_put(&cost_seg, &costs, row, col);
+
+			if (dir == 1) {
+			    if (dir_bin)
+				cur_dir = (1 << (int)cur_dir);
+			    Segment_put(&dir_seg, &cur_dir, row, col);
+			}
+		    }
+		}
+
+		if (equal) {
+		    /* this can create circular paths:
+		     * set only if current cell does not point to neighbor
+		     * does not avoid longer circular paths */
+		    Segment_get(&dir_seg, &old_dir, pres_cell->row, pres_cell->col);
+		    dir_fwd = (1 << dir_inv[(int)cur_dir]);
+		    if (!((int)old_dir & dir_fwd)) {
+			Segment_get(&dir_seg, &old_dir, row, col);
+			cur_dir = ((1 << (int)cur_dir) | (int)old_dir);
+			Segment_put(&dir_seg, &cur_dir, row, col);
+		    }
 		}
 	    }
 	}
@@ -1113,7 +1187,11 @@ int main(int argc, char *argv[])
 
     /* free heap */
     free_heap();
-    
+
+    if (have_solver) {
+	Segment_close(&solve_seg);
+    }
+
     /* Open cumulative cost layer for writing */
     cum_fd = Rast_open_new(cum_cost_layer, cum_data_type);
     cell = Rast_allocate_buf(cum_data_type);

raster/r.cost/r.cost.html

@@ -76,7 +76,7 @@ also considered.
 Knight's move example:
 <center>
 <img src="rcost_knightsmove.png" border="1"><br>
-<table border="0" width="590">
+<table border="0" width="600" height="221">
 <tr><td><center>
 <i>Flat cost surface without (left pane) and with the knight's move (right pane).
 The default is to grow the cost outwards in 8 directions.
@@ -86,10 +86,32 @@ Using the knight's move grows it outwards in 16 directions.</i>
 </center>
 
 <p>
-If the <em>nearest</em> output parameter is specified, the module will calculate 
+If the <b>nearest</b> output parameter is specified, the module will calculate 
 for each cell its nearest starting point based on the minimized accumulative cost
 while moving over the cost map.
 
+<p>
+The <b>solver</b> option helps to select a particular direction in case 
+of multiple directions with equal costs. Sometimes fields with equal 
+cumulative costs exist and multiple paths with equal costs would lead 
+from a start point to a stop point. By default, a path along the edge 
+of such a field would be produced or multiple paths of equal costs with 
+the <b>-b</b> flag. An additional variable can be supplied with the 
+<b>solver</b> option to help the algorithm pick a particular direction.
+<p>
+Example for solving multiple directions:
+<center>
+<img src="rcost_solvedir.png" border="1"><br>
+<table border="0" width="400" height="206">
+<tr><td><center>
+<i>A field of equal cumulative costs with multiple paths (black). By 
+default a path along the edge will be selected (red). Path selection 
+can be controlled with the solver option (blue).</i>
+</center></td></tr>
+</table>
+</center>
+
+
 <h2>NULL CELLS</h2>
 
 By default null cells in the input raster map are excluded from
@@ -114,24 +136,23 @@ follows:
 The user generates a raster map indicating the cost of
 traversing each cell in the north-south and east-west directions.
 This map, along with a set of starting points are submitted to
-<em>r.cost</em>. The starting points are put into a list cells from which
+<em>r.cost</em>. The starting points are put into a heap of cells from which
 costs to the adjacent cells are to be calculated. The cell on the
-list with the lowest cumulative cost is selected for computing costs to
+heap with the lowest cumulative cost is selected for computing costs to
 the neighboring cells. Costs are computed and those cells are put
-on the list and the originating cell is finalized. This process
+on the heap and the originating cell is finalized. This process
 of selecting the lowest cumulative cost cell, computing costs to the
-neighbors, putting the neighbors on the list and removing the
-originating cell from the list continues until the list is empty.
+neighbors, putting the neighbors on the heap and removing the
+originating cell from the heap continues until the heap is empty.
 <p>
 The most time consuming aspect of this algorithm is the management of
-the list of cells for which cumulative costs have been at least
-initially computed. <em>r.cost</em> uses a binary tree with an linked list
-at each node in the tree for efficiently holding cells with identical
-cumulative costs.
+the heap of cells for which cumulative costs have been at least
+initially computed. <em>r.cost</em> uses a minimum heap for efficiently 
+tracking the next cell with the lowest cumulative costs.
 <p>
 <em>r.cost</em>, like most all GRASS raster programs, is also made to 
 be run on maps larger that can fit in available computer memory. As the 
-algorithm works through the dynamic list of cells it can move almost 
+algorithm works through the dynamic heap of cells it can move almost 
 randomly around the entire area. <em>r.cost</em> divides the entire 
 area into a number of pieces and swaps these pieces in and out of 
 memory (to and from disk) as needed. This provides a virtual memory 

raster/r.walk/main.c

@@ -139,7 +139,8 @@ int main(int argc, char *argv[])
     const char *cost_layer, *dtm_layer;
     const char *dtm_mapset, *cost_mapset, *search_mapset;
     void *dtm_cell, *cost_cell, *cum_cell, *dir_cell, *cell2 = NULL;
-    SEGMENT cost_seg, dir_seg;
+    SEGMENT cost_seg, dir_seg, solve_seg;
+    int have_solver;
     double *value;
     char buf[400];
     extern struct Cell_head window;
@@ -150,7 +151,7 @@ int main(int argc, char *argv[])
     double zero = 0.0;
     int col, row, nrows, ncols;
     int maxcost, par_number;
-    int nseg;
+    int nseg, nbytes;
     int maxmem;
     int segments_in_memory;
     int cost_fd, cum_fd, dtm_fd, dir_fd;
@@ -169,6 +170,7 @@ int main(int argc, char *argv[])
     struct Flag *flag2, *flag3, *flag4, *flag5, *flag6;
     struct Option *opt1, *opt2, *opt3, *opt4, *opt5, *opt6, *opt7, *opt8;
     struct Option *opt9, *opt10, *opt11, *opt12, *opt13, *opt14, *opt15;
+    struct Option *opt_solve;
     struct cost *pres_cell;
     struct start_pt *head_start_pt = NULL;
     struct start_pt *next_start_pt;
@@ -186,6 +188,7 @@ int main(int argc, char *argv[])
     int dtm_dsize, cost_dsize;
     double disk_mb, mem_mb, pq_mb;
     int dir_bin;
+    DCELL mysolvedir[2], solvedir[2];
 
     /* Definition for dimension and region check */
     struct Cell_head dtm_cellhd, cost_cellhd;
@@ -213,6 +216,14 @@ int main(int argc, char *argv[])
     opt1 = G_define_standard_option(G_OPT_R_OUTPUT);
     opt1->description = _("Name for output raster map to contain walking costs");
 
+    opt_solve = G_define_standard_option(G_OPT_R_INPUT);
+    opt_solve->key = "solver";
+    opt_solve->required = NO;
+    opt_solve->label =
+	_("Name of input raster map solving equal costs");
+    opt_solve->description =
+	_("Helper variable to pick a direction if two directions have equal cumulative costs (smaller is better)");
+
     opt11 = G_define_standard_option(G_OPT_R_OUTPUT);
     opt11->key = "outdir";
     opt11->required = NO;
@@ -438,6 +449,14 @@ int main(int argc, char *argv[])
 	    G_fatal_error(_("Vector map <%s> not found"), opt7->answer);
     }
 
+    have_solver = 0;
+    if (dir && opt_solve->answer) {
+	search_mapset = G_find_raster2(opt_solve->answer, "");
+	if (search_mapset == NULL)
+	    G_fatal_error(_("Raster map <%s> not found"), opt_solve->answer);
+	have_solver = 1;
+    }
+
     if (!Rast_is_d_null_value(&null_cost)) {
 	if (null_cost < 0.0) {
 	    G_warning(_("Assigning negative cost to null cell. Null cells excluded."));
@@ -528,26 +547,21 @@ int main(int argc, char *argv[])
     maxmem -= pq_mb;
     if (maxmem < 10)
 	maxmem = 10;
-    if (dir == TRUE) {
-	disk_mb = (double) nrows * ncols * 28. / 1048576.;
-	segments_in_memory = maxmem / 
-	                     ((double) srows * scols * (28. / 1048576.));
-	if (segments_in_memory < 4)
-	    segments_in_memory = 4;
-	if (segments_in_memory > nseg)
-	    segments_in_memory = nseg;
-	mem_mb = (double) srows * scols * (28. / 1048576.) * segments_in_memory;
-    }
-    else {
-	disk_mb = (double) nrows * ncols * 24. / 1048576.;
-	segments_in_memory = maxmem / 
-	                     ((double) srows * scols * (24. / 1048576.));
-	if (segments_in_memory < 4)
-	    segments_in_memory = 4;
-	if (segments_in_memory > nseg)
-	    segments_in_memory = nseg;
-	mem_mb = (double) srows * scols * (24. / 1048576.) * segments_in_memory;
-    }
+
+    nbytes = 24;
+    if (dir == TRUE)
+	nbytes += 4;
+    if (have_solver)
+	nbytes += 16;
+
+    disk_mb = (double) nrows * ncols * nbytes / 1048576.;
+    segments_in_memory = maxmem / 
+			 ((double) srows * scols * (nbytes / 1048576.));
+    if (segments_in_memory < 4)
+	segments_in_memory = 4;
+    if (segments_in_memory > nseg)
+	segments_in_memory = nseg;
+    mem_mb = (double) srows * scols * (nbytes / 1048576.) * segments_in_memory;
 
     if (flag5->answer) {
 	fprintf(stdout, _("Will need at least %.2f MB of disk space"), disk_mb);
@@ -582,6 +596,32 @@ int main(int argc, char *argv[])
 	    G_fatal_error(_("Can not create temporary file"));
     }
 
+    if (have_solver) {
+	int sfd, dsize;
+	void *cell;
+
+	if (Segment_open(&solve_seg, G_tempfile(), nrows, ncols, srows, scols,
+		         sizeof(DCELL) * 2, segments_in_memory) != 1)
+	    G_fatal_error(_("Can not create temporary file"));
+
+	sfd = Rast_open_old(opt_solve->answer, "");
+	cell = Rast_allocate_buf(DCELL_TYPE);
+	Rast_set_d_null_value(&solvedir[1], 1); 
+	dsize = Rast_cell_size(DCELL_TYPE);
+	for (row = 0; row < nrows; row++) {
+	    G_percent(row, nrows, 2);
+	    Rast_get_d_row(sfd, cell, row);
+	    ptr2 = cell;
+	    for (col = 0; col < ncols; col++) {
+		solvedir[0] = *(DCELL *)ptr2;
+		Segment_put(&solve_seg, solvedir, row, col);
+		ptr2 = G_incr_void_ptr(ptr2, dsize);
+	    }
+	}
+	Rast_close(sfd);
+	G_free(cell);
+    }
+
     /* Write the dtm and cost layers in the segmented file */
     G_message(_("Reading raster maps <%s> and <%s>, initializing output..."),
 	      G_fully_qualified_name(dtm_layer, dtm_mapset),
@@ -947,6 +987,9 @@ int main(int argc, char *argv[])
 	    continue;
 	}
 
+	if (have_solver)
+	    Segment_get(&solve_seg, mysolvedir, pres_cell->row, pres_cell->col);
+
 	row = pres_cell->row;
 	col = pres_cell->col;
 
@@ -1419,6 +1462,11 @@ int main(int argc, char *argv[])
 			cur_dir = (1 << (int)cur_dir);
 		    Segment_put(&dir_seg, &cur_dir, row, col);
 		}
+		if (have_solver) {
+		    Segment_get(&solve_seg, solvedir, row, col);
+		    solvedir[1] = mysolvedir[0];
+		    Segment_put(&solve_seg, solvedir, row, col);
+		}
 	    }
 	    /* update with lower costs */
 	    else if (old_min_cost > min_cost) {
@@ -1430,22 +1478,47 @@ int main(int argc, char *argv[])
 			cur_dir = (1 << (int)cur_dir);
 		    Segment_put(&dir_seg, &cur_dir, row, col);
 		}
+		if (have_solver) {
+		    Segment_get(&solve_seg, solvedir, row, col);
+		    solvedir[1] = mysolvedir[0];
+		    Segment_put(&solve_seg, solvedir, row, col);
+		}
 	    }
 	    else if (dir && dir_bin && old_min_cost == min_cost) {
 		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;
 
-		/* this can create circular paths:
-		 * set only if current cell does not point to neighbor
-                 * does not avoid longer circular paths */
-		Segment_get(&dir_seg, &old_dir, pres_cell->row, pres_cell->col);
-		dir_fwd = (1 << dir_inv[(int)cur_dir]);
-		if (!((int)old_dir & dir_fwd)) {
-		    Segment_get(&dir_seg, &old_dir, row, col);
-		    cur_dir = ((1 << (int)cur_dir) | (int)old_dir);
-		    Segment_put(&dir_seg, &cur_dir, row, col);
+		if (have_solver) {
+		    Segment_get(&solve_seg, solvedir, row, col);
+		    equal = (solvedir[1] == mysolvedir[0]);
+		    if (solvedir[1] > mysolvedir[0]) {
+			solvedir[1] = mysolvedir[0];
+			Segment_put(&solve_seg, solvedir, row, col);
+
+			Segment_put(&cost_seg, &costs, row, col);
+
+			if (dir == 1) {
+			    if (dir_bin)
+				cur_dir = (1 << (int)cur_dir);
+			    Segment_put(&dir_seg, &cur_dir, row, col);
+			}
+		    }
+		}
+
+		if (equal) {
+		    /* this can create circular paths:
+		     * set only if current cell does not point to neighbor
+		     * does not avoid longer circular paths */
+		    Segment_get(&dir_seg, &old_dir, pres_cell->row, pres_cell->col);
+		    dir_fwd = (1 << dir_inv[(int)cur_dir]);
+		    if (!((int)old_dir & dir_fwd)) {
+			Segment_get(&dir_seg, &old_dir, row, col);
+			cur_dir = ((1 << (int)cur_dir) | (int)old_dir);
+			Segment_put(&dir_seg, &cur_dir, row, col);
+		    }
 		}
 	    }
 	}
@@ -1464,6 +1537,10 @@ int main(int argc, char *argv[])
 
     /* free heap */
     free_heap();
+
+    if (have_solver) {
+	Segment_close(&solve_seg);
+    }
     
     /* Open cumulative cost layer for writing */
     cum_fd = Rast_open_new(cum_cost_layer, cum_data_type);