[GRASS-SVN] r72720 - in grass/trunk/raster: r.cost r.walk
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat May 19 08:50:47 PDT 2018
Author: mmetz
Date: 2018-05-19 08:50:47 -0700 (Sat, 19 May 2018)
New Revision: 72720
Added:
grass/trunk/raster/r.cost/rcost_solvedir.png
Modified:
grass/trunk/raster/r.cost/main.c
grass/trunk/raster/r.cost/r.cost.html
grass/trunk/raster/r.walk/main.c
Log:
r.cost/r.walk: add option to control which direction is used in case of multiple directions with equal costs
Modified: grass/trunk/raster/r.cost/main.c
===================================================================
--- grass/trunk/raster/r.cost/main.c 2018-05-19 15:11:58 UTC (rev 72719)
+++ grass/trunk/raster/r.cost/main.c 2018-05-19 15:50:47 UTC (rev 72720)
@@ -105,7 +105,8 @@
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 @@
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 @@
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 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 @@
"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 @@
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 @@
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,27 +434,22 @@
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);
fprintf(stdout, "\n");
@@ -467,6 +482,31 @@
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 @@
}
}
+ 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 @@
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,23 +1127,49 @@
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 @@
/* 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);
Modified: grass/trunk/raster/r.cost/r.cost.html
===================================================================
--- grass/trunk/raster/r.cost/r.cost.html 2018-05-19 15:11:58 UTC (rev 72719)
+++ grass/trunk/raster/r.cost/r.cost.html 2018-05-19 15:50:47 UTC (rev 72720)
@@ -76,7 +76,7 @@
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 @@
</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 @@
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
Added: grass/trunk/raster/r.cost/rcost_solvedir.png
===================================================================
(Binary files differ)
Index: grass/trunk/raster/r.cost/rcost_solvedir.png
===================================================================
--- grass/trunk/raster/r.cost/rcost_solvedir.png 2018-05-19 15:11:58 UTC (rev 72719)
+++ grass/trunk/raster/r.cost/rcost_solvedir.png 2018-05-19 15:50:47 UTC (rev 72720)
Property changes on: grass/trunk/raster/r.cost/rcost_solvedir.png
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+image/png
\ No newline at end of property
Modified: grass/trunk/raster/r.walk/main.c
===================================================================
--- grass/trunk/raster/r.walk/main.c 2018-05-19 15:11:58 UTC (rev 72719)
+++ grass/trunk/raster/r.walk/main.c 2018-05-19 15:50:47 UTC (rev 72720)
@@ -139,7 +139,8 @@
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 @@
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 @@
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 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 @@
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 @@
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,27 +547,22 @@
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);
fprintf(stdout, "\n");
@@ -582,6 +596,32 @@
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 @@
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 @@
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,6 +1478,11 @@
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;
@@ -1436,17 +1489,37 @@
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 @@
/* 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);
More information about the grass-commit
mailing list