Index: main.c =================================================================== --- main.c (revision 34032) +++ main.c (working copy) @@ -7,6 +7,9 @@ * James Westervelt - CERL * Pierre de Mouveaux * Eric G. Miller + * + * Updated for calculation errors and directional surface generation + * Colin Nielsen * * PURPOSE: Outputs a raster map layer showing the cumulative cost * of moving between different geographic locations on an @@ -72,25 +75,26 @@ int main(int argc, char *argv[]) { - void *cell, *cell2; - SEGMENT in_seg, out_seg; - char *cost_mapset; - char *in_file, *out_file; + void *cell, *cell2, *dir_cell; + SEGMENT in_seg, out_seg, out_seg2; + char *cost_mapset, *move_dir_mapset; + char *in_file, *out_file, *dir_out_file; char *search_mapset; double *value; extern struct Cell_head window; double NS_fac, EW_fac, DIAG_fac, H_DIAG_fac, V_DIAG_fac; double fcost; double min_cost, old_min_cost; + double cur_dir, old_cur_dir; double zero = 0.0; int at_percent = 0; int col, row, nrows, ncols; int maxcost; int maxmem; double cost; - int cost_fd, cum_fd; + int cost_fd, cum_fd, dir_fd; int have_stop_points = 0; - int in_fd, out_fd; + int in_fd, out_fd, dir_out_fd; double my_cost; double null_cost; int srows, scols; @@ -104,13 +108,13 @@ struct GModule *module; struct Flag *flag2, *flag3, *flag4; struct Option *opt1, *opt2, *opt3, *opt4, *opt5, *opt6, *opt7, *opt8; - struct Option *opt9, *opt10; + struct Option *opt9, *opt10, *opt11; struct cost *pres_cell, *new_cell; struct start_pt *pres_start_pt = NULL; struct start_pt *pres_stop_pt = NULL; void *ptr2; - RASTER_MAP_TYPE data_type; + RASTER_MAP_TYPE data_type, dir_data_type = DCELL_TYPE; struct History history; double peak = 0.0; int dsize; @@ -198,6 +202,14 @@ opt10->answer = "100"; opt10->description = _("Percent of map to keep in memory"); + opt11 = G_define_option(); + opt11->key = "outdir"; + opt11->type = TYPE_STRING; + opt11->required = NO; + opt11->gisprompt = "new,cell,raster"; + opt11->description = + _("Name of output raster map to contain movement directions"); + flag2 = G_define_flag(); flag2->key = 'k'; flag2->description = @@ -231,6 +243,7 @@ in_file = G_tempfile(); out_file = G_tempfile(); + dir_out_file = G_tempfile(); /* Get database window parameters */ @@ -321,11 +334,19 @@ if (cost_mapset == NULL) G_fatal_error(_("Raster map <%s> not found"), cost_layer); + strcpy(move_dir_layer, opt11->answer); + + search_mapset = ""; + move_dir_mapset = G_find_cell2(move_dir_layer, search_mapset); + /* Check if specified output layer name is legal */ if (G_legal_filename(cum_cost_layer) < 0) G_fatal_error(_("<%s> is an illegal file name"), cum_cost_layer); + if (G_legal_filename(move_dir_layer) < 0) + G_fatal_error(_("<%s> is an illegal file name"), move_dir_layer); + /* Find number of rows and columns in window */ nrows = G_window_rows(); @@ -379,6 +400,10 @@ segment_format(out_fd, nrows, ncols, srows, scols, sizeof(double)); close(out_fd); + dir_out_fd = creat(dir_out_file, 0600); + segment_format(dir_out_fd, nrows, ncols, srows, scols, sizeof(double)); + close(dir_out_fd); + /* Open initialize and segment all files */ in_fd = open(in_file, 2); @@ -388,6 +413,9 @@ out_fd = open(out_file, 2); segment_init(&out_seg, out_fd, segments_in_memory); + dir_out_fd = open(dir_out_file, 2); + segment_init(&out_seg2, dir_out_fd, segments_in_memory); + /* Write the cost layer in the segmented file */ G_message(_("Reading %s"), cost_layer); @@ -479,6 +507,32 @@ G_free(fbuff); } + G_message(_("Initializing directional output ")); + { + double *fbuff; + int i; + + fbuff = (double *)G_malloc((unsigned int)(ncols * sizeof(double))); + + if (fbuff == NULL) + G_fatal_error(_ + ("Unable to allocate memory for segment fbuff == NULL")); + + G_set_d_null_value(fbuff, ncols); + + for (row = 0; row < nrows; row++) { + { + G_percent(row, nrows, 2); + } + for (i = 0; i < ncols; i++) { + segment_put(&out_seg2, &fbuff[i], row, i); + } + } + segment_flush(&out_seg2); + G_percent(row, nrows, 2); + G_free(fbuff); + } + /* Scan the start_points layer searching for starting points. * Create a btree of starting points ordered by increasing costs. */ @@ -706,55 +760,71 @@ case 1: row = pres_cell->row; col = pres_cell->col - 1; + cur_dir = 180.0; break; case 2: col = pres_cell->col + 1; + cur_dir = 0.0; break; case 3: row = pres_cell->row - 1; col = pres_cell->col; + cur_dir = 90.0; break; case 4: row = pres_cell->row + 1; + cur_dir = 270.0; break; case 5: row = pres_cell->row - 1; col = pres_cell->col - 1; + cur_dir = 135.0; break; case 6: col = pres_cell->col + 1; + cur_dir = 45.0; break; case 7: row = pres_cell->row + 1; + cur_dir = 315.0; break; case 8: col = pres_cell->col - 1; + cur_dir = 225.0; break; case 9: row = pres_cell->row - 2; col = pres_cell->col - 1; + cur_dir = 112.5; break; case 10: col = pres_cell->col + 1; + cur_dir = 67.5; break; case 11: row = pres_cell->row + 2; + cur_dir = 292.5; break; case 12: col = pres_cell->col - 1; + cur_dir = 247.5; break; case 13: row = pres_cell->row - 1; col = pres_cell->col - 2; + cur_dir = 157.5; break; case 14: col = pres_cell->col + 2; + cur_dir = 22.5; break; case 15: row = pres_cell->row + 1; + cur_dir = 337.5; break; case 16: col = pres_cell->col - 2; + cur_dir = 202.5; break; } @@ -769,97 +839,113 @@ case 1: value = &W; segment_get(&in_seg, value, row, col); - fcost = (double)(W + my_cost) / 2.0; + fcost = (double)((W / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost * EW_fac; break; case 2: value = &E; segment_get(&in_seg, value, row, col); - fcost = (double)(E + my_cost) / 2.0; + fcost = (double)((E / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost * EW_fac; break; case 3: value = &N; segment_get(&in_seg, value, row, col); - fcost = (double)(N + my_cost) / 2.0; + fcost = (double)((N / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost * NS_fac; break; case 4: value = &S; segment_get(&in_seg, value, row, col); - fcost = (double)(S + my_cost) / 2.0; + fcost = (double)((S / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost * NS_fac; break; case 5: value = &NW; segment_get(&in_seg, value, row, col); - fcost = (double)(NW + my_cost) / 2.0; + fcost = (double)((NW / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost * DIAG_fac; break; case 6: value = &NE; segment_get(&in_seg, value, row, col); - fcost = (double)(NE + my_cost) / 2.0; + fcost = (double)((NE / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost * DIAG_fac; break; case 7: value = &SE; segment_get(&in_seg, value, row, col); - fcost = (double)(SE + my_cost) / 2.0; + fcost = (double)((SE / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost * DIAG_fac; break; case 8: value = &SW; segment_get(&in_seg, value, row, col); - fcost = (double)(SW + my_cost) / 2.0; + fcost = (double)((SW / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost * DIAG_fac; break; case 9: value = &NNW; segment_get(&in_seg, value, row, col); - fcost = (double)(N + NW + NNW + my_cost) / 4.0; + fcost = + (double)((N / 4.0) + (NW / 4.0) + (NNW / 4.0) + + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost * V_DIAG_fac; break; case 10: value = &NNE; segment_get(&in_seg, value, row, col); - fcost = (double)(N + NE + NNE + my_cost) / 4.0; + fcost = + (double)((N / 4.0) + (NE / 4.0) + (NNE / 4.0) + + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost * V_DIAG_fac; break; case 11: value = &SSE; segment_get(&in_seg, value, row, col); - fcost = (double)(S + SE + SSE + my_cost) / 4.0; + fcost = + (double)((S / 4.0) + (SE / 4.0) + (SSE / 4.0) + + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost * V_DIAG_fac; break; case 12: value = &SSW; segment_get(&in_seg, value, row, col); - fcost = (double)(S + SW + SSW + my_cost) / 4.0; + fcost = + (double)((S / 4.0) + (SW / 4.0) + (SSW / 4.0) + + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost * V_DIAG_fac; break; case 13: value = &WNW; segment_get(&in_seg, value, row, col); - fcost = (double)(W + NW + WNW + my_cost) / 4.0; + fcost = + (double)((W / 4.0) + (NW / 4.0) + (WNW / 4.0) + + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost * H_DIAG_fac; break; case 14: value = &ENE; segment_get(&in_seg, value, row, col); - fcost = (double)(E + NE + ENE + my_cost) / 4.0; + fcost = + (double)((E / 4.0) + (NE / 4.0) + (ENE / 4.0) + + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost * H_DIAG_fac; break; case 15: value = &ESE; segment_get(&in_seg, value, row, col); - fcost = (double)(E + SE + ESE + my_cost) / 4.0; + fcost = + (double)((E / 4.0) + (SE / 4.0) + (ESE / 4.0) + + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost * H_DIAG_fac; break; case 16: value = &WSW; segment_get(&in_seg, value, row, col); - fcost = (double)(W + SW + WSW + my_cost) / 4.0; + fcost = + (double)((W / 4.0) + (SW / 4.0) + (WSW / 4.0) + + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost * H_DIAG_fac; break; } @@ -868,15 +954,18 @@ continue; segment_get(&out_seg, &old_min_cost, row, col); + segment_get(&out_seg2, &old_cur_dir, row, col); if (G_is_d_null_value(&old_min_cost)) { segment_put(&out_seg, &min_cost, row, col); new_cell = insert(min_cost, row, col); + segment_put(&out_seg2, &cur_dir, row, col); } else { if (old_min_cost > min_cost) { segment_put(&out_seg, &min_cost, row, col); new_cell = insert(min_cost, row, col); + segment_put(&out_seg2, &cur_dir, row, col); } else { } @@ -902,9 +991,13 @@ cum_fd = G_open_raster_new(cum_cost_layer, data_type); + dir_fd = G_open_raster_new(move_dir_layer, dir_data_type); + dir_cell = G_allocate_raster_buf(dir_data_type); + /* Write pending updates by segment_put() to output map */ segment_flush(&out_seg); + segment_flush(&out_seg2); /* Copy segmented map to output map */ G_message(_("Writing %s"), cum_cost_layer); @@ -1014,23 +1107,44 @@ } } + double *p; + + G_message(_("Writing movement direction file %s..."), move_dir_layer); + for (row = 0; row < nrows; row++) { + p = dir_cell; + for (col = 0; col < ncols; col++) { + segment_get(&out_seg2, &cur_dir, row, col); + *(p + col) = cur_dir; + } + G_put_raster_row(dir_fd, dir_cell, dir_data_type); + G_percent(row, nrows, 2); + } + G_percent(row, nrows, 2); G_message(_("Peak cost value: %f"), peak); segment_release(&in_seg); /* release memory */ segment_release(&out_seg); + segment_release(&out_seg2); G_close_cell(cost_fd); G_close_cell(cum_fd); + G_close_cell(dir_fd); close(in_fd); /* close all files */ close(out_fd); + close(dir_out_fd); unlink(in_file); /* remove submatrix files */ unlink(out_file); + unlink(dir_out_file); G_short_history(cum_cost_layer, "raster", &history); G_command_history(&history); G_write_history(cum_cost_layer, &history); + G_short_history(move_dir_layer, "raster", &history); + G_command_history(&history); + G_write_history(move_dir_layer, &history); + /* Create colours for output map */ /* Index: description.html =================================================================== --- description.html (revision 34032) +++ description.html (working copy) @@ -6,11 +6,13 @@ other user-specified cell(s) whose locations are specified by their geographic coordinate(s). Each cell in the original cost surface map will contain a category value which represents the cost of traversing -that cell. r.cost will produce an output raster map in +that cell. r.cost will produce 1) an output raster map in which each cell contains the lowest total cost of traversing the -space between each cell and the user-specified points. (Diagonal +space between each cell and the user-specified points (diagonal costs are multiplied by a factor that depends on the dimensions of -the cell.) This program uses the current geographic region settings. +the cell) and 2) a second raster map layer showing the movement +direction to the next cell on the path back to the start point (see +Movement Direction). This program uses the current geographic region settings. The output map will be of the same data format as the input map, integer or floating point.

@@ -18,7 +20,8 @@ The input name is the name of a raster map whose category values represent the surface cost. The output name is the name of the -resultant raster map of cumulative cost. +resultant raster map of cumulative cost. The outdir name is the +name of the resultant raster map of movement directions (see Movement Direction).

r.cost can be run with three different methods of identifying the @@ -101,6 +104,7 @@ propagate the adjacent costs. These cells can be retained as null cells in the output map by using the -n flag. +

NOTES

Sometimes, when the differences among integer cell category values in the @@ -203,11 +207,11 @@ The output map can be viewed, for example, as an elevation model in which the starting location(s) is/are the lowest point(s). Outputs from r.cost -can be used as inputs to r.drain, in order -to trace the least-cost path given by this model between any given cell -and the r.cost starting location(s). The two programs, when -used together, generate least-cost paths or corridors between any two -map locations (cells). +can be used as inputs to r.drain with +the direction flag -d, in order to trace the least-cost path given by this +model between any given cell and the r.cost starting location(s). The +two programs, when used together, generate least-cost paths or corridors between any +two map locations (cells).

Shortest distance surfaces

The r.cost module allows for computing the shortest distance @@ -228,6 +232,26 @@ d.rast dist_meters + +

Movement Direction

+

+The movement direction surface is created to record the sequence of +movements that created the cost accumulation surface. Without it +r.drain would not correctly create a path from an end point +back to the start point. The direction shown in each cell points away +from the cell that came before it. The directions are recorded as +GRASS standard directions:

+       112.5 90  67.5         i.e. a cell with the value 135 
+157.5  135   0   45   22.5    means the cell before it is 
+       180   x   0            to the south-east.
+202.5  225  270  315  337.5
+       247.5     292.5
+
+

+Once r.cost computes the cumulative cost map, r.drain +can be used to find the minimum cost path. Make sure to use the -d flag +and the movement direction raster map when running r.drain to ensure +the path is computed according to the proper movement directions.

BUGS

Index: stash.h =================================================================== --- stash.h (revision 34032) +++ stash.h (working copy) @@ -29,6 +29,7 @@ #define CUM_COST_LAYER 1 #define COST_LAYER 2 #define START_PT 3 +#define MOVE_DIR_LAYER 4 struct start_pt { @@ -49,17 +50,18 @@ { "output", CUM_COST_LAYER}, { "input", COST_LAYER}, { - "coor", START_PT} + "coor", START_PT}, { + "outdir", MOVE_DIR_LAYER} }; -char cum_cost_layer[64]; +char cum_cost_layer[64], move_dir_layer[64]; char cost_layer[64]; struct start_pt *head_start_pt = NULL; struct start_pt *head_end_pt = NULL; #else -extern char cum_cost_layer[]; +extern char cum_cost_layer[], move_dir_layer[]; extern char cost_layer[]; extern struct start_pt *head_start_pt; extern struct start_pt *head_end_pt;