1
2 /****************************************************************************
3 *
4 * MODULE: r.walk
5 * AUTHOR(S): Based on r.cost written by :
6 * Antony Awaida,
7 * Intelligent Engineering
8 * Systems Laboratory,
9 * M.I.T.
10 * James Westervelt,
11 * U.S.Army Construction Engineering Research Laboratory
12 *
13 * Updated for Grass 5
14 * Pierre de Mouveaux (pmx@audiovu.com)
15 *
16 * Initial version of r.walk:
17 * Steno Fontanari, 2002, ITC-irst
18 *
19 * GRASS 6.0 version of r.walk:
20 * Franceschetti Simone, Sorrentino Diego, Mussi Fabiano and Pasolli Mattia
21 * Correction by: Fontanari Steno, Napolitano Maurizio and Flor Roberto
22 * In collaboration with: Franchi Matteo, Vaglia Beatrice, Bartucca Luisa,
23 * Fava Valentina and Tolotti Mathias, 2004
24 *
25 * Updated for GRASS 6.1
26 * Roberto Flor and Markus Neteler
27 * Glynn Clements <glynn gclements.plus.com>, Soeren Gebbert <soeren.gebbert gmx.de>
28 * Updated for calculation errors and directional surface generation
29 * Colin Nielsen <colin.nielsen gmail com>
30 * Use min heap instead of btree (faster, less memory)
31 * multiple directions with bitmask encoding
32 * avoid circular paths
33 * Markus Metz
34 * PURPOSE: anisotropic movements on cost surfaces
35 * COPYRIGHT: (C) 1999-2015 by the GRASS Development Team
36 *
37 * This program is free software under the GNU General Public
38 * License (>=v2). Read the file COPYING that comes with GRASS
39 * for details.
40 *
41 ***************************************************************************/
42
43 /*********************************************************************
44 *
45 * This is the main program for the minimum path cost analysis.
46 * It generates a cumulative cost map (output) from an elevation (inputdtm
47 * and cost map (inputcost) with respect to starting locations (coor).
48 *
49 * It takes as input the following:
50 * 1) Cost of traversing each grid cell as given by an elevation map and
51 * a cost map cell (inputcost).
52 * 2) If starting points are not specified on the command line
53 * then the output map must exist and contain the starting locations
54 *
55 * Otherwise the output map need not exist and the coor points
56 * from the command line are used.
57 *
58 *********************************************************************/
59
60 /*********************************************************************
61 * The walking energy is computed for the human walk, based on Aitken,
62 * 1977, Langmuir, 1984:
63 *
64 * {T= [(a)x(Delta S)] + [(b)x(Delta H Climb)]
65 * +[(c)*(Delta H moderate downhill)]+[(d)*(Delta H steep downhill]}
66 *
67 * where T is time in seconds, Delta S distance in meter, Delta H the heigth difference
68 *
69 * The default a,b,c,d parameters used below have been measured using oxygen consumption in biomechanical
70 * experiments.
71 * Refs:
72 * * Aitken, R. 1977. Wilderness areas in Scotland. Unpublished Ph.D. thesis. University of Aberdeen.
73 * * Steno Fontanari, University of Trento, Italy, Ingegneria per l'Ambiente e
74 * il Territorio, 2000-2001. Svilluppo di metodologie GIS per la determinazione dell'accessibilita'
75 * territoriale come supporto alle decisioni nella gestione ambientale.
76 * * Langmuir, E. 1984. Mountaincraft and leadership. The Scottish Sports Council/MLTB. Cordee,
77 * Leicester.
78 *
79 * The total cost is computed as a linear combination of walking energy and a given friction cost map:
80 *
81 * TOTAL COST = [(WALKING ENERGY ) + (LAMBDA*FRICTION)]
82 *
83 * TODO: generalize formula to other species
84 *************/
85 /*
86 *
87 * 20 july 2004 - Pierre de Mouveaux. pmx@audiovu.com
88 * Updated to use the Grass 5.0 floating point raster cell format.
89 * Convert floats to double. Done ;)
90 * 2001: original r.walk by Steno Fontanari, ITC-irst
91 * 24 July 2004: WebValley 2004, fixed and enhanced by
92 * Matteo Franchi Liceo Leonardo Da Vinci Trento
93 * Roberto Flor ITC-irst
94 * 7 December 2005: Grass 6.1 cleanup
95 * Roberto Flor ITC-irst
96 * Markus Neteler CEA
97 */
98
99 #include <stdlib.h>
100 #include <unistd.h>
101 #include <string.h>
102 #include <math.h>
103 #include <sys/types.h>
104 #include <sys/stat.h>
105 #include <fcntl.h>
106 #include <grass/gis.h>
107 #include <grass/raster.h>
108 #include <grass/vector.h>
109 #include <grass/segment.h>
110 #include <grass/glocale.h>
111 #include "cost.h"
112 #include "stash.h"
113 #include "flag.h"
114
115 #define SEGCOLSIZE 64
116
117 struct Cell_head window;
118
119 struct rc
120 {
121 int r;
122 int c;
123 };
124
125 static struct rc *stop_pnts = NULL;
126 static int n_stop_pnts = 0;
127 static int stop_pnts_alloc = 0;
128
cmp_rc(struct rc * a,struct rc * b)129 int cmp_rc(struct rc *a, struct rc *b)
130 {
131 if (a->r == b->r)
132 return (a->c - b->c);
133
134 return (a->r - b->r);
135 }
136
137 void add_stop_pnt(int r, int c);
138
main(int argc,char * argv[])139 int main(int argc, char *argv[])
140 {
141 const char *cum_cost_layer, *move_dir_layer, *nearest_layer;
142 const char *cost_layer, *dtm_layer;
143 const char *dtm_mapset, *cost_mapset, *search_mapset;
144 void *dtm_cell, *cost_cell, *cum_cell, *dir_cell, *cell2 = NULL, *nearest_cell;
145 SEGMENT cost_seg, dir_seg, solve_seg;
146 int have_solver;
147 double *value;
148 char buf[400];
149 extern struct Cell_head window;
150 double NS_fac, EW_fac, DIAG_fac, H_DIAG_fac, V_DIAG_fac;
151 double fcost_dtm, fcost_cost;
152 double min_cost, old_min_cost;
153 FCELL cur_dir;
154 double zero = 0.0;
155 int col, row, nrows, ncols;
156 int maxcost, par_number;
157 int nseg, nbytes;
158 int maxmem;
159 int segments_in_memory;
160 int cost_fd, cum_fd, dtm_fd, dir_fd, nearest_fd;
161 int dir = 0;
162 double my_dtm, my_cost, check_dtm, nearest;
163 double null_cost, dnullval;
164 double a, b, c, d, lambda, slope_factor;
165 int srows, scols;
166 int total_reviewed;
167 int keep_nulls = 1;
168 int start_with_raster_vals = 1;
169 int neighbor;
170 long n_processed = 0;
171 long total_cells;
172 struct GModule *module;
173 struct Flag *flag2, *flag3, *flag4, *flag5, *flag6;
174 struct Option *opt1, *opt2, *opt3, *opt4, *opt5, *opt6, *opt7, *opt8;
175 struct Option *opt9, *opt10, *opt11, *opt12, *opt13, *opt14, *opt15, *opt16;
176 struct Option *opt_solve;
177 struct cost *pres_cell;
178 struct start_pt *head_start_pt = NULL;
179 struct start_pt *next_start_pt;
180 struct cc {
181 double dtm; /* elevation model */
182 double cost_in; /* friction costs */
183 double cost_out; /* cumulative costs */
184 double nearest; /* nearest start point */
185 } costs;
186 FLAG *visited;
187
188 void *ptr1, *ptr2;
189 RASTER_MAP_TYPE dtm_data_type, cost_data_type, cum_data_type =
190 DCELL_TYPE, dir_data_type = FCELL_TYPE,
191 nearest_data_type = CELL_TYPE; /* output nearest type */
192 struct History history;
193 double peak = 0.0;
194 int dtm_dsize, cost_dsize, nearest_size;
195 double disk_mb, mem_mb, pq_mb;
196
197 int dir_bin;
198 DCELL mysolvedir[2], solvedir[2];
199
200 /* Definition for dimension and region check */
201 struct Cell_head dtm_cellhd, cost_cellhd;
202
203 G_gisinit(argv[0]);
204
205 module = G_define_module();
206 G_add_keyword(_("raster"));
207 G_add_keyword(_("cost surface"));
208 G_add_keyword(_("cumulative costs"));
209 G_add_keyword(_("cost allocation"));
210 module->description =
211 _("Creates a raster map showing the "
212 "anisotropic cumulative cost of moving between different "
213 "geographic locations on an input raster map "
214 "whose cell category values represent cost.");
215
216 opt12 = G_define_standard_option(G_OPT_R_ELEV);
217
218 opt2 = G_define_standard_option(G_OPT_R_INPUT);
219 opt2->key = "friction";
220 opt2->description =
221 _("Name of input raster map containing friction costs");
222
223 opt1 = G_define_standard_option(G_OPT_R_OUTPUT);
224 opt1->description = _("Name for output raster map to contain walking costs");
225
226 opt_solve = G_define_standard_option(G_OPT_R_INPUT);
227 opt_solve->key = "solver";
228 opt_solve->required = NO;
229 opt_solve->label =
230 _("Name of input raster map solving equal costs");
231 opt_solve->description =
232 _("Helper variable to pick a direction if two directions have equal cumulative costs (smaller is better)");
233
234 opt16 = G_define_standard_option(G_OPT_R_OUTPUT);
235 opt16->key = "nearest";
236 opt16->required = NO;
237 opt16->description =
238 _("Name for output raster map with nearest start point");
239 opt16->guisection = _("Optional outputs");
240
241 opt11 = G_define_standard_option(G_OPT_R_OUTPUT);
242 opt11->key = "outdir";
243 opt11->required = NO;
244 opt11->description =
245 _("Name for output raster map to contain movement directions");
246 opt11->guisection = _("Optional outputs");
247
248 opt7 = G_define_standard_option(G_OPT_V_INPUT);
249 opt7->key = "start_points";
250 opt7->required = NO;
251 opt7->label = _("Name of starting vector points map");
252 opt7->guisection = _("Start");
253
254 opt8 = G_define_standard_option(G_OPT_V_INPUT);
255 opt8->key = "stop_points";
256 opt8->required = NO;
257 opt8->label = _("Name of stopping vector points map");
258 opt8->guisection = _("Stop");
259
260 opt9 = G_define_standard_option(G_OPT_R_INPUT);
261 opt9->key = "start_raster";
262 opt9->required = NO;
263 opt9->description = _("Name of starting raster points map");
264 opt9->guisection = _("Start");
265
266 opt3 = G_define_standard_option(G_OPT_M_COORDS);
267 opt3->key = "start_coordinates";
268 opt3->multiple = YES;
269 opt3->description =
270 _("Coordinates of starting point(s) (E,N)");
271 opt3->guisection = _("Start");
272
273 opt4 = G_define_standard_option(G_OPT_M_COORDS);
274 opt4->key = "stop_coordinates";
275 opt4->multiple = YES;
276 opt4->description =
277 _("Coordinates of stopping point(s) (E,N)");
278 opt4->guisection = _("Stop");
279
280 opt5 = G_define_option();
281 opt5->key = "max_cost";
282 opt5->type = TYPE_INTEGER;
283 opt5->key_desc = "value";
284 opt5->required = NO;
285 opt5->multiple = NO;
286 opt5->answer = "0";
287 opt5->description = _("Maximum cumulative cost");
288
289 opt6 = G_define_option();
290 opt6->key = "null_cost";
291 opt6->type = TYPE_DOUBLE;
292 opt6->key_desc = "value";
293 opt6->required = NO;
294 opt6->multiple = NO;
295 opt6->description =
296 _("Cost assigned to null cells. By default, null cells are excluded");
297 opt6->guisection = _("NULL cells");
298
299 opt10 = G_define_standard_option(G_OPT_MEMORYMB);
300
301 opt15 = G_define_option();
302 opt15->key = "walk_coeff";
303 opt15->type = TYPE_STRING;
304 opt15->key_desc = "a,b,c,d";
305 opt15->required = NO;
306 opt15->multiple = NO;
307 opt15->answer = "0.72,6.0,1.9998,-1.9998";
308 opt15->description =
309 _("Coefficients for walking energy formula parameters a,b,c,d");
310 opt15->guisection = _("Settings");
311
312 opt14 = G_define_option();
313 opt14->key = "lambda";
314 opt14->type = TYPE_DOUBLE;
315 opt14->required = NO;
316 opt14->multiple = NO;
317 opt14->answer = "1.0";
318 opt14->description =
319 _("Lambda coefficients for combining walking energy and friction cost");
320 opt14->guisection = _("Settings");
321
322 opt13 = G_define_option();
323 opt13->key = "slope_factor";
324 opt13->type = TYPE_DOUBLE;
325 opt13->required = NO;
326 opt13->multiple = NO;
327 opt13->answer = "-0.2125";
328 opt13->description =
329 _("Slope factor determines travel energy cost per height step");
330 opt13->guisection = _("Settings");
331
332 flag2 = G_define_flag();
333 flag2->key = 'k';
334 flag2->description =
335 _("Use the 'Knight's move'; slower, but more accurate");
336
337 flag3 = G_define_flag();
338 flag3->key = 'n';
339 flag3->description = _("Keep null values in output raster map");
340 flag3->guisection = _("NULL cells");
341
342 flag4 = G_define_flag();
343 flag4->key = 'r';
344 flag4->description = _("Start with values in raster map");
345 flag4->guisection = _("Start");
346
347 flag5 = G_define_flag();
348 flag5->key = 'i';
349 flag5->description = _("Print info about disk space and memory requirements and exit");
350
351 flag6 = G_define_flag();
352 flag6->key = 'b';
353 flag6->description = _("Create bitmask encoded directions");
354 flag6->guisection = _("Optional outputs");
355
356 /* Parse options */
357 if (G_parser(argc, argv))
358 exit(EXIT_FAILURE);
359
360 /* If no outdir is specified, set flag to skip all dir */
361 if (opt11->answer != NULL)
362 dir = 1;
363
364 /* Get database window parameters */
365 Rast_get_window(&window);
366
367 /* Find north-south, east_west and diagonal factors */
368 EW_fac = window.ew_res; /* Must be the physical distance */
369 NS_fac = window.ns_res;
370 DIAG_fac = (double)sqrt((double)(NS_fac * NS_fac + EW_fac * EW_fac));
371 V_DIAG_fac =
372 (double)sqrt((double)(4 * NS_fac * NS_fac + EW_fac * EW_fac));
373 H_DIAG_fac =
374 (double)sqrt((double)(NS_fac * NS_fac + 4 * EW_fac * EW_fac));
375
376 Rast_set_d_null_value(&null_cost, 1);
377
378 if (flag2->answer)
379 total_reviewed = 16;
380 else
381 total_reviewed = 8;
382
383 keep_nulls = flag3->answer;
384
385 start_with_raster_vals = flag4->answer;
386
387 dir_bin = flag6->answer;
388
389 {
390 int count = 0;
391
392 if (opt3->answers)
393 count++;
394 if (opt7->answers)
395 count++;
396 if (opt9->answers)
397 count++;
398
399 if (count != 1)
400 G_fatal_error(_("Must specify exactly one of start_points, start_rast or coordinate"));
401 }
402
403 if (opt3->answers) {
404 head_start_pt = process_start_coords(opt3->answers, head_start_pt);
405 if (!head_start_pt)
406 G_fatal_error(_("No start points"));
407 }
408
409 if (opt4->answers) {
410 if (!process_stop_coords(opt4->answers))
411 G_fatal_error(_("No stop points"));
412 }
413
414 if (sscanf(opt5->answer, "%d", &maxcost) != 1 || maxcost < 0)
415 G_fatal_error(_("Inappropriate maximum cost: %d"), maxcost);
416
417 if (sscanf(opt10->answer, "%d", &maxmem) != 1 || maxmem <= 0)
418 G_fatal_error(_("Inappropriate amount of memory: %d"), maxmem);
419
420 /* Getting walking energy formula parameters */
421 if ((par_number =
422 sscanf(opt15->answer, "%lf,%lf,%lf,%lf", &a, &b, &c, &d)) != 4)
423 G_fatal_error(_("Missing required value: got %d instead of 4"),
424 par_number);
425 else {
426 G_message(_("Walking costs are a=%g b=%g c=%g d=%g"), a, b, c, d);
427 }
428
429 /* Getting lambda */
430 if ((par_number = sscanf(opt14->answer, "%lf", &lambda)) != 1)
431 G_fatal_error(_("Missing required value: %d"), par_number);
432 else {
433
434 G_message(_("Lambda is %g"), lambda);
435 }
436
437 /* Getting slope_factor */
438 if ((par_number = sscanf(opt13->answer, "%lf", &slope_factor)) != 1)
439 G_fatal_error(_("Missing required value: %d"), par_number);
440 else {
441
442 G_message(_("Slope_factor is %g"), slope_factor);
443 }
444
445 if ((opt6->answer == NULL) ||
446 (sscanf(opt6->answer, "%lf", &null_cost) != 1)) {
447 G_debug(1, "Null cells excluded from cost evaluation");
448 Rast_set_d_null_value(&null_cost, 1);
449 }
450 else if (keep_nulls)
451 G_debug(1,"Input null cell will be retained into output map");
452
453 if (opt7->answer) {
454 search_mapset = G_find_vector2(opt7->answer, "");
455 if (search_mapset == NULL)
456 G_fatal_error(_("Vector map <%s> not found"), opt7->answer);
457 }
458
459 have_solver = 0;
460 if (dir && opt_solve->answer) {
461 search_mapset = G_find_raster2(opt_solve->answer, "");
462 if (search_mapset == NULL)
463 G_fatal_error(_("Raster map <%s> not found"), opt_solve->answer);
464 have_solver = 1;
465 }
466
467 if (!Rast_is_d_null_value(&null_cost)) {
468 if (null_cost < 0.0) {
469 G_warning(_("Assigning negative cost to null cell. Null cells excluded."));
470 Rast_set_d_null_value(&null_cost, 1);
471 }
472 }
473 else {
474 keep_nulls = 0; /* handled automagically... */
475 }
476
477 cum_cost_layer = opt1->answer;
478 cost_layer = opt2->answer;
479 move_dir_layer = opt11->answer;
480 dtm_layer = opt12->answer;
481 nearest_layer = opt16->answer;
482
483 /* Find number of rows and columns in window */
484 nrows = Rast_window_rows();
485 ncols = Rast_window_cols();
486
487 /* Open cost cell layer for reading */
488 dtm_mapset = G_find_raster2(dtm_layer, "");
489 if (dtm_mapset == NULL)
490 G_fatal_error(_("Raster map <%s> not found"), dtm_layer);
491 dtm_fd = Rast_open_old(dtm_layer, "");
492
493 cost_mapset = G_find_raster2(cost_layer, "");
494 if (cost_mapset == NULL)
495 G_fatal_error(_("Raster map <%s> not found"), cost_layer);
496 cost_fd = Rast_open_old(cost_layer, cost_mapset);
497
498 Rast_get_cellhd(dtm_layer, "", &dtm_cellhd);
499 Rast_get_cellhd(cost_layer, "", &cost_cellhd);
500
501 dtm_data_type = Rast_get_map_type(dtm_fd);
502 cost_data_type = Rast_get_map_type(cost_fd);
503
504 /* Parameters for map submatrices */
505 switch (dtm_data_type) {
506 case (CELL_TYPE):
507 G_debug(1, "DTM_Source map is: Integer cell type");
508 break;
509 case (FCELL_TYPE):
510 G_debug(1, "DTM_Source map is: Floating point (float) cell type");
511 break;
512 case (DCELL_TYPE):
513 G_debug(1, "DTM_Source map is: Floating point (double) cell type");
514 break;
515 }
516 G_debug(1, "DTM %d rows, %d cols", dtm_cellhd.rows, dtm_cellhd.cols);
517
518 switch (cost_data_type) {
519 case (CELL_TYPE):
520 G_debug(1, "COST_Source map is: Integer cell type");
521 break;
522 case (FCELL_TYPE):
523 G_debug(1, "COST_Source map is: Floating point (float) cell type");
524 break;
525 case (DCELL_TYPE):
526 G_debug(1, "COST_Source map is: Floating point (double) cell type");
527 break;
528 }
529 G_debug(1, "COST %d rows, %d cols", cost_cellhd.rows, cost_cellhd.cols);
530
531 G_debug(1, " %d rows, %d cols", nrows, ncols);
532 G_format_resolution(window.ew_res, buf, window.proj);
533 G_debug(1, " EW resolution %s (%g)", buf, window.ew_res);
534 G_format_resolution(window.ns_res, buf, window.proj);
535 G_debug(1, " NS resolution %s (%g)", buf, window.ns_res);
536
537 /* this is most probably the limitation of r.walk for large datasets
538 * segment size needs to be reduced to avoid unnecessary disk IO
539 * but it doesn't make sense to go down to 1
540 * so use 64 segment rows and cols for <= 200 million cells
541 * for larger regions, 32 segment rows and cols
542 * maybe go down to 16 for > 500 million cells ? */
543 if ((double) nrows * ncols > 200000000)
544 srows = scols = SEGCOLSIZE / 2;
545 else
546 srows = scols = SEGCOLSIZE;
547
548 /* calculate total number of segments */
549 nseg = ((nrows + srows - 1) / srows) * ((ncols + scols - 1) / scols);
550
551 /* calculate disk space and memory requirements */
552 /* (nrows + ncols) * 8. * 20.0 / 1048576. for Dijkstra search */
553 pq_mb = ((double)nrows + ncols) * 8. * 20.0 / 1048576.;
554 G_debug(1, "pq MB: %g", pq_mb);
555 maxmem -= pq_mb;
556 if (maxmem < 10)
557 maxmem = 10;
558
559 nbytes = 24;
560 if (dir == TRUE)
561 nbytes += 4;
562 if (have_solver)
563 nbytes += 16;
564
565 disk_mb = (double) nrows * ncols * nbytes / 1048576.;
566 segments_in_memory = maxmem /
567 ((double) srows * scols * (nbytes / 1048576.));
568 if (segments_in_memory < 4)
569 segments_in_memory = 4;
570 if (segments_in_memory > nseg)
571 segments_in_memory = nseg;
572 mem_mb = (double) srows * scols * (nbytes / 1048576.) * segments_in_memory;
573
574 if (flag5->answer) {
575 fprintf(stdout, _("Will need at least %.2f MB of disk space"), disk_mb);
576 fprintf(stdout, "\n");
577 fprintf(stdout, _("Will need at least %.2f MB of memory"), mem_mb);
578 fprintf(stdout, "\n");
579 fprintf(stdout, _("%d of %d segments are kept in memory"),
580 segments_in_memory, nseg);
581 fprintf(stdout, "\n");
582 Rast_close(cost_fd);
583 Rast_close(dtm_fd);
584 exit(EXIT_SUCCESS);
585 }
586
587 G_verbose_message("--------------------------------------------");
588 G_verbose_message(_("Will need at least %.2f MB of disk space"), disk_mb);
589 G_verbose_message(_("Will need at least %.2f MB of memory"), mem_mb);
590 G_verbose_message(_("%d of %d segments are kept in memory"),
591 segments_in_memory, nseg);
592 G_verbose_message("--------------------------------------------");
593
594 /* Create segmented format files for cost layer and output layer */
595 G_verbose_message(_("Creating some temporary files..."));
596
597 if (Segment_open(&cost_seg, G_tempfile(), nrows, ncols, srows, scols,
598 sizeof(struct cc), segments_in_memory) != 1)
599 G_fatal_error(_("Can not create temporary file"));
600
601 if (dir == 1) {
602 if (Segment_open(&dir_seg, G_tempfile(), nrows, ncols, srows, scols,
603 sizeof(FCELL), segments_in_memory) != 1)
604 G_fatal_error(_("Can not create temporary file"));
605 }
606
607 if (have_solver) {
608 int sfd, dsize;
609 void *cell;
610
611 if (Segment_open(&solve_seg, G_tempfile(), nrows, ncols, srows, scols,
612 sizeof(DCELL) * 2, segments_in_memory) != 1)
613 G_fatal_error(_("Can not create temporary file"));
614
615 sfd = Rast_open_old(opt_solve->answer, "");
616 cell = Rast_allocate_buf(DCELL_TYPE);
617 Rast_set_d_null_value(&solvedir[1], 1);
618 dsize = Rast_cell_size(DCELL_TYPE);
619 for (row = 0; row < nrows; row++) {
620 G_percent(row, nrows, 2);
621 Rast_get_d_row(sfd, cell, row);
622 ptr2 = cell;
623 for (col = 0; col < ncols; col++) {
624 solvedir[0] = *(DCELL *)ptr2;
625 if (Segment_put(&solve_seg, solvedir, row, col) < 0)
626 G_fatal_error(_("Can not write to temporary file"));
627 ptr2 = G_incr_void_ptr(ptr2, dsize);
628 }
629 }
630 Rast_close(sfd);
631 G_free(cell);
632 }
633
634 /* Write the dtm and cost layers in the segmented file */
635 G_message(_("Reading raster maps <%s> and <%s>, initializing output..."),
636 G_fully_qualified_name(dtm_layer, dtm_mapset),
637 G_fully_qualified_name(cost_layer, cost_mapset));
638
639 /* read required maps cost and dtm */
640 {
641 int skip_nulls;
642 double p_dtm, p_cost;
643
644 Rast_set_d_null_value(&dnullval, 1);
645 costs.cost_out = dnullval;
646 costs.nearest = 0;
647
648 total_cells = nrows * ncols;
649
650 skip_nulls = Rast_is_d_null_value(&null_cost);
651
652 dtm_dsize = Rast_cell_size(dtm_data_type);
653 cost_dsize = Rast_cell_size(cost_data_type);
654 dtm_cell = Rast_allocate_buf(dtm_data_type);
655 cost_cell = Rast_allocate_buf(cost_data_type);
656 p_dtm = 0.0;
657 p_cost = 0.0;
658
659 for (row = 0; row < nrows; row++) {
660 G_percent(row, nrows, 2);
661 Rast_get_row(dtm_fd, dtm_cell, row, dtm_data_type);
662 Rast_get_row(cost_fd, cost_cell, row, cost_data_type);
663 /* INPUT NULL VALUES: ??? */
664 ptr1 = cost_cell;
665 ptr2 = dtm_cell;
666
667 for (col = 0; col < ncols; col++) {
668 if (Rast_is_null_value(ptr1, cost_data_type)) {
669 p_cost = null_cost;
670 if (skip_nulls) {
671 total_cells--;
672 }
673 }
674 else {
675 switch (cost_data_type) {
676 case CELL_TYPE:
677 p_cost = *(CELL *)ptr1;
678 break;
679 case FCELL_TYPE:
680 p_cost = *(FCELL *)ptr1;
681 break;
682 case DCELL_TYPE:
683 p_cost = *(DCELL *)ptr1;
684 break;
685 }
686 }
687 costs.cost_in = p_cost;
688
689 if (Rast_is_null_value(ptr2, dtm_data_type)) {
690 p_dtm = null_cost;
691 if (skip_nulls && !Rast_is_null_value(ptr1, cost_data_type)) {
692 total_cells--;
693 }
694 }
695 else {
696 switch (dtm_data_type) {
697 case CELL_TYPE:
698 p_dtm = *(CELL *)ptr2;
699 break;
700 case FCELL_TYPE:
701 p_dtm = *(FCELL *)ptr2;
702 break;
703 case DCELL_TYPE:
704 p_dtm = *(DCELL *)ptr2;
705 break;
706 }
707 }
708
709 costs.dtm = p_dtm;
710 if (Segment_put(&cost_seg, &costs, row, col) < 0)
711 G_fatal_error(_("Can not write to temporary file"));
712 ptr1 = G_incr_void_ptr(ptr1, cost_dsize);
713 ptr2 = G_incr_void_ptr(ptr2, dtm_dsize);
714 }
715 }
716 G_free(dtm_cell);
717 G_free(cost_cell);
718 G_percent(1, 1, 1);
719 }
720
721 if (dir == 1) {
722 FCELL fnullval;
723
724 G_message(_("Initializing directional output..."));
725 Rast_set_f_null_value(&fnullval, 1);
726 for (row = 0; row < nrows; row++) {
727 G_percent(row, nrows, 2);
728 for (col = 0; col < ncols; col++) {
729 if (Segment_put(&dir_seg, &fnullval, row, col) < 0)
730 G_fatal_error(_("Can not write to temporary file"));
731 }
732 }
733 G_percent(1, 1, 1);
734 }
735 /* Scan the start_points layer searching for starting points.
736 * Create a heap of starting points ordered by increasing costs.
737 */
738 init_heap();
739
740 /* read vector with start points */
741 if (opt7->answer) {
742 struct Map_info In;
743 struct line_pnts *Points;
744 struct line_cats *Cats;
745 struct bound_box box;
746 int cat, type, npoints = 0;
747
748 Points = Vect_new_line_struct();
749 Cats = Vect_new_cats_struct();
750
751 Vect_set_open_level(1); /* topology not required */
752
753 if (1 > Vect_open_old(&In, opt7->answer, ""))
754 G_fatal_error(_("Unable to open vector map <%s>"), opt7->answer);
755
756 G_message(_("Reading vector map <%s> with start points..."),
757 Vect_get_full_name(&In));
758
759 Vect_rewind(&In);
760
761 Vect_region_box(&window, &box);
762
763 while (1) {
764 /* register line */
765 type = Vect_read_next_line(&In, Points, Cats);
766
767 /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
768 if (type == -1) {
769 G_warning(_("Unable to read vector map"));
770 continue;
771 }
772 else if (type == -2) {
773 break;
774 }
775 if (!Vect_point_in_box(Points->x[0], Points->y[0], 0, &box))
776 continue;
777 npoints++;
778
779 col = (int)Rast_easting_to_col(Points->x[0], &window);
780 row = (int)Rast_northing_to_row(Points->y[0], &window);
781
782 next_start_pt =
783 (struct start_pt *)(G_malloc(sizeof(struct start_pt)));
784
785 next_start_pt->row = row;
786 next_start_pt->col = col;
787 Vect_cat_get(Cats, 1, &cat);
788 next_start_pt->value = cat;
789 next_start_pt->next = head_start_pt;
790 head_start_pt = next_start_pt;
791 }
792
793 if (npoints < 1)
794 G_fatal_error(_("No start points found in vector map <%s>"), Vect_get_full_name(&In));
795 else
796 G_verbose_message(n_("%d point found", "%d points found", npoints), npoints);
797
798 Vect_close(&In);
799 }
800
801 /* read vector with stop points */
802 if (opt8->answer) {
803 struct Map_info In;
804 struct line_pnts *Points;
805 struct line_cats *Cats;
806 struct bound_box box;
807 int type;
808
809 G_message(_("Reading vector map <%s> with stop points..."), opt8->answer);
810
811 Points = Vect_new_line_struct();
812 Cats = Vect_new_cats_struct();
813
814 Vect_set_open_level(1); /* topology not required */
815
816 if (1 > Vect_open_old(&In, opt8->answer, ""))
817 G_fatal_error(_("Unable to open vector map <%s>"), opt8->answer);
818
819 Vect_rewind(&In);
820
821 Vect_region_box(&window, &box);
822
823 while (1) {
824 /* register line */
825 type = Vect_read_next_line(&In, Points, Cats);
826
827 /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
828 if (type == -1) {
829 G_warning(_("Unable to read vector map"));
830 continue;
831 }
832 else if (type == -2) {
833 break;
834 }
835 if (!Vect_point_in_box(Points->x[0], Points->y[0], 0, &box))
836 continue;
837
838 col = (int)Rast_easting_to_col(Points->x[0], &window);
839 row = (int)Rast_northing_to_row(Points->y[0], &window);
840
841 add_stop_pnt(row, col);
842 }
843
844 Vect_close(&In);
845
846 if (!stop_pnts)
847 G_fatal_error(_("No stop points found in vector <%s>"), opt8->answer);
848 }
849
850 /* read raster with start points */
851 if (opt9->answer) {
852 int dsize2;
853 int fd;
854 RASTER_MAP_TYPE data_type2;
855 int got_one = 0;
856
857 search_mapset = G_find_raster(opt9->answer, "");
858
859 if (search_mapset == NULL)
860 G_fatal_error(_("Raster map <%s> not found"), opt9->answer);
861
862 fd = Rast_open_old(opt9->answer, search_mapset);
863 data_type2 = Rast_get_map_type(fd);
864 nearest_data_type = data_type2;
865 dsize2 = Rast_cell_size(data_type2);
866 cell2 = Rast_allocate_buf(data_type2);
867 if (!cell2)
868 G_fatal_error(_("Unable to allocate memory"));
869
870 G_message(_("Reading raster map <%s> with start points..."), opt9->answer);
871 for (row = 0; row < nrows; row++) {
872 G_percent(row, nrows, 2);
873 Rast_get_row(fd, cell2, row, data_type2);
874 ptr2 = cell2;
875 for (col = 0; col < ncols; col++) {
876 /* Did I understand that concept of cumulative cost map? - (pmx) 12 april 2000 */
877 if (!Rast_is_null_value(ptr2, data_type2)) {
878 double cellval;
879
880 if (Segment_get(&cost_seg, &costs, row, col) < 0)
881 G_fatal_error(_("Can not read from temporary file"));
882
883 cellval = Rast_get_d_value(ptr2, data_type2);
884 if (start_with_raster_vals == 1) {
885 insert(cellval, row, col);
886 costs.cost_out = cellval;
887 costs.nearest = cellval;
888 if (Segment_put(&cost_seg, &costs, row, col) < 0)
889 G_fatal_error(_("Can not write to temporary file"));
890 }
891 else {
892 value = &zero;
893 insert(zero, row, col);
894 costs.cost_out = *value;
895 costs.nearest = cellval;
896 if (Segment_put(&cost_seg, &costs, row, col) < 0)
897 G_fatal_error(_("Can not write to temporary file"));
898 }
899 got_one = 1;
900 }
901 ptr2 = G_incr_void_ptr(ptr2, dsize2);
902 }
903 }
904 G_percent(1, 1, 1);
905
906 Rast_close(fd);
907 G_free(cell2);
908
909 if (!got_one)
910 G_fatal_error(_("No start points"));
911 }
912
913 /* Insert start points into min heap */
914 if (head_start_pt) {
915
916 next_start_pt = head_start_pt;
917 while (next_start_pt != NULL) {
918 value = &zero;
919 if (next_start_pt->row < 0 || next_start_pt->row >= nrows
920 || next_start_pt->col < 0 || next_start_pt->col >= ncols)
921 G_fatal_error(_("Specified starting location outside database window"));
922 insert(zero, next_start_pt->row, next_start_pt->col);
923 if (Segment_get(&cost_seg, &costs, next_start_pt->row,
924 next_start_pt->col) < 0)
925 G_fatal_error(_("Can not read from temporary file"));
926 costs.cost_out = *value;
927 costs.nearest = next_start_pt->value;
928
929 if (Segment_put(&cost_seg, &costs, next_start_pt->row,
930 next_start_pt->col) < 0)
931 G_fatal_error(_("Can not write to temporary file"));
932 next_start_pt = next_start_pt->next;
933 }
934 }
935
936 if (n_stop_pnts > 1) {
937 int i, j;
938
939 /* prune stop points */
940 j = 1;
941 for (i = 1; i < n_stop_pnts; i++) {
942 if (stop_pnts[i].r != stop_pnts[j - 1].r ||
943 stop_pnts[i].c != stop_pnts[j - 1].c) {
944 stop_pnts[j].r = stop_pnts[i].r;
945 stop_pnts[j].c = stop_pnts[i].c;
946 j++;
947 }
948 }
949 if (n_stop_pnts > j) {
950 G_message(_("Number of duplicate stop points: %d"), n_stop_pnts - j);
951 n_stop_pnts = j;
952 }
953 }
954
955 /* Loop through the heap and perform at each cell the following:
956 * 1) If an adjacent cell has not already been assigned a value compute
957 * the min cost and assign it.
958 * 2) Insert the adjacent cell in the heap.
959 * 3) Free the memory allocated to the present cell.
960 */
961
962 G_debug(1, "total cells: %ld", total_cells);
963 G_debug(1, "nrows x ncols: %d", nrows * ncols);
964 G_message(_("Finding cost path..."));
965 n_processed = 0;
966 visited = flag_create(nrows, ncols);
967
968 pres_cell = get_lowest();
969 while (pres_cell != NULL) {
970 struct cost *ct;
971 double N_dtm, NE_dtm, E_dtm, SE_dtm, S_dtm, SW_dtm, W_dtm, NW_dtm;
972 double NNE_dtm, ENE_dtm, ESE_dtm, SSE_dtm, SSW_dtm, WSW_dtm, WNW_dtm,
973 NNW_dtm;
974 double N_cost, NE_cost, E_cost, SE_cost, S_cost, SW_cost, W_cost,
975 NW_cost;
976 double NNE_cost, ENE_cost, ESE_cost, SSE_cost, SSW_cost, WSW_cost,
977 WNW_cost, NNW_cost;
978
979 N_dtm = NE_dtm = E_dtm = SE_dtm = S_dtm = SW_dtm = W_dtm = NW_dtm = dnullval;
980 NNE_dtm = ENE_dtm = ESE_dtm = SSE_dtm = SSW_dtm = WSW_dtm = WNW_dtm = NNW_dtm = dnullval;
981
982 N_cost = NE_cost = E_cost = SE_cost = S_cost = SW_cost = W_cost = NW_cost = dnullval;
983 NNE_cost = ENE_cost = ESE_cost = SSE_cost = SSW_cost = WSW_cost = WNW_cost = NNW_cost = dnullval;
984
985 /* If we have surpassed the user specified maximum cost, then quit */
986 if (maxcost && ((double)maxcost < pres_cell->min_cost))
987 break;
988
989 /* If I've already been updated, delete me */
990 if (Segment_get(&cost_seg, &costs, pres_cell->row, pres_cell->col) < 0)
991 G_fatal_error(_("Can not read from temporary file"));
992 old_min_cost = costs.cost_out;
993 if (!Rast_is_d_null_value(&old_min_cost)) {
994 if (pres_cell->min_cost > old_min_cost) {
995 delete(pres_cell);
996 pres_cell = get_lowest();
997 continue;
998 }
999 }
1000 my_dtm = costs.dtm;
1001 if (Rast_is_d_null_value(&my_dtm)) {
1002 delete(pres_cell);
1003 pres_cell = get_lowest();
1004 continue;
1005 }
1006 my_cost = costs.cost_in;
1007 if (Rast_is_d_null_value(&my_cost)) {
1008 delete(pres_cell);
1009 pres_cell = get_lowest();
1010 continue;
1011 }
1012 if (FLAG_GET(visited, pres_cell->row, pres_cell->col)) {
1013 delete(pres_cell);
1014 pres_cell = get_lowest();
1015 continue;
1016 }
1017 FLAG_SET(visited, pres_cell->row, pres_cell->col);
1018
1019 if (have_solver) {
1020 if (Segment_get(&solve_seg, mysolvedir, pres_cell->row, pres_cell->col) < 0)
1021 G_fatal_error(_("Can not read from temporary file"));
1022 }
1023
1024 nearest = costs.nearest;
1025
1026 row = pres_cell->row;
1027 col = pres_cell->col;
1028
1029 G_percent(n_processed++, total_cells, 1);
1030
1031 /* 9 10 Order in which neighbors
1032 * 13 5 3 6 14 are visited (Knight move).
1033 * 1 2
1034 * 16 8 4 7 15
1035 * 12 11
1036 */
1037
1038 /* drainage directions in degrees CCW from East
1039 * drainage directions are set for each neighbor and must be
1040 * read as from neighbor to current cell
1041 *
1042 * X = neighbor:
1043 *
1044 * 112.5 67.5
1045 * 157.5 135 90 45 22.5
1046 * 180 X 360
1047 * 202.5 225 270 315 337.5
1048 * 247.5 292.5
1049 *
1050 * X = current cell:
1051 *
1052 * 292.5 247.5
1053 * 337.5 315 270 225 202.5
1054 * 360 X 180
1055 * 22.5 45 90 135 157.5
1056 * 67.5 112.5
1057 */
1058
1059 /* drainage directions bitmask encoded CW from North
1060 * drainage directions are set for each neighbor and must be
1061 * read as from neighbor to current cell
1062 *
1063 * bit positions, zero-based, from neighbor to current cell
1064 *
1065 * X = neighbor X = current cell
1066 *
1067 * 15 8 11 12
1068 * 14 6 7 0 9 10 2 3 4 13
1069 * 5 X 1 1 X 5
1070 * 13 4 3 2 10 9 0 7 6 14
1071 * 12 11 8 15
1072 */
1073
1074 for (neighbor = 1; neighbor <= total_reviewed; neighbor++) {
1075 switch (neighbor) {
1076 case 1:
1077 row = pres_cell->row;
1078 col = pres_cell->col - 1;
1079 cur_dir = 360.0;
1080 if (dir_bin)
1081 cur_dir = 1;
1082 break;
1083 case 2:
1084 row = pres_cell->row;
1085 col = pres_cell->col + 1;
1086 cur_dir = 180.0;
1087 if (dir_bin)
1088 cur_dir = 5;
1089 break;
1090 case 3:
1091 row = pres_cell->row - 1;
1092 col = pres_cell->col;
1093 cur_dir = 270.0;
1094 if (dir_bin)
1095 cur_dir = 3;
1096 break;
1097 case 4:
1098 row = pres_cell->row + 1;
1099 col = pres_cell->col;
1100 cur_dir = 90.0;
1101 if (dir_bin)
1102 cur_dir = 7;
1103 break;
1104 case 5:
1105 row = pres_cell->row - 1;
1106 col = pres_cell->col - 1;
1107 cur_dir = 315.0;
1108 if (dir_bin)
1109 cur_dir = 2;
1110 break;
1111 case 6:
1112 row = pres_cell->row - 1;
1113 col = pres_cell->col + 1;
1114 cur_dir = 225.0;
1115 if (dir_bin)
1116 cur_dir = 4;
1117 break;
1118 case 7:
1119 row = pres_cell->row + 1;
1120 col = pres_cell->col + 1;
1121 cur_dir = 135.0;
1122 if (dir_bin)
1123 cur_dir = 6;
1124 break;
1125 case 8:
1126 row = pres_cell->row + 1;
1127 col = pres_cell->col - 1;
1128 cur_dir = 45.0;
1129 if (dir_bin)
1130 cur_dir = 0;
1131 break;
1132 case 9:
1133 row = pres_cell->row - 2;
1134 col = pres_cell->col - 1;
1135 cur_dir = 292.5;
1136 if (dir_bin)
1137 cur_dir = 11;
1138 break;
1139 case 10:
1140 row = pres_cell->row - 2;
1141 col = pres_cell->col + 1;
1142 cur_dir = 247.5;
1143 if (dir_bin)
1144 cur_dir = 12;
1145 break;
1146 case 11:
1147 row = pres_cell->row + 2;
1148 col = pres_cell->col + 1;
1149 cur_dir = 112.5;
1150 if (dir_bin)
1151 cur_dir = 15;
1152 break;
1153 case 12:
1154 row = pres_cell->row + 2;
1155 col = pres_cell->col - 1;
1156 cur_dir = 67.5;
1157 if (dir_bin)
1158 cur_dir = 8;
1159 break;
1160 case 13:
1161 row = pres_cell->row - 1;
1162 col = pres_cell->col - 2;
1163 cur_dir = 337.5;
1164 if (dir_bin)
1165 cur_dir = 10;
1166 break;
1167 case 14:
1168 row = pres_cell->row - 1;
1169 col = pres_cell->col + 2;
1170 cur_dir = 202.5;
1171 if (dir_bin)
1172 cur_dir = 13;
1173 break;
1174 case 15:
1175 row = pres_cell->row + 1;
1176 col = pres_cell->col + 2;
1177 cur_dir = 157.5;
1178 if (dir_bin)
1179 cur_dir = 14;
1180 break;
1181 case 16:
1182 row = pres_cell->row + 1;
1183 col = pres_cell->col - 2;
1184 cur_dir = 22.5;
1185 if (dir_bin)
1186 cur_dir = 9;
1187 break;
1188 }
1189
1190 if (row < 0 || row >= nrows)
1191 continue;
1192 if (col < 0 || col >= ncols)
1193 continue;
1194
1195 /* skip already processed neighbors here ? */
1196
1197 min_cost = dnullval;
1198 if (Segment_get(&cost_seg, &costs, row, col) < 0)
1199 G_fatal_error(_("Can not read from temporary file"));
1200
1201 switch (neighbor) {
1202 case 1:
1203 W_dtm = costs.dtm;
1204 W_cost = costs.cost_in;
1205 if (Rast_is_d_null_value(&W_cost))
1206 continue;
1207 check_dtm = (W_dtm - my_dtm) / EW_fac;
1208 if (check_dtm >= 0)
1209 fcost_dtm = (double)(W_dtm - my_dtm) * b;
1210 else if (check_dtm < (slope_factor))
1211 fcost_dtm = (double)(W_dtm - my_dtm) * d;
1212 else
1213 fcost_dtm = (double)(W_dtm - my_dtm) * c;
1214 fcost_cost = (double)(W_cost + my_cost) / 2.0;
1215 min_cost =
1216 pres_cell->min_cost + fcost_dtm + (EW_fac * a) +
1217 lambda * fcost_cost * EW_fac;
1218 break;
1219 case 2:
1220 E_dtm = costs.dtm;
1221 E_cost = costs.cost_in;
1222 if (Rast_is_d_null_value(&E_cost))
1223 continue;
1224 check_dtm = (E_dtm - my_dtm) / EW_fac;
1225 if (check_dtm >= 0)
1226 fcost_dtm = (double)(E_dtm - my_dtm) * b;
1227 else if (check_dtm < (slope_factor))
1228 fcost_dtm = (double)(E_dtm - my_dtm) * d;
1229 else
1230 fcost_dtm = (double)(E_dtm - my_dtm) * c;
1231 fcost_cost = (double)(E_cost + my_cost) / 2.0;
1232 min_cost =
1233 pres_cell->min_cost + fcost_dtm + (EW_fac * a) +
1234 lambda * fcost_cost * EW_fac;
1235 break;
1236 case 3:
1237 N_dtm = costs.dtm;
1238 N_cost = costs.cost_in;
1239 if (Rast_is_d_null_value(&N_cost))
1240 continue;
1241 check_dtm = (N_dtm - my_dtm) / NS_fac;
1242 if (check_dtm >= 0)
1243 fcost_dtm = (double)(N_dtm - my_dtm) * b;
1244 else if (check_dtm < (slope_factor))
1245 fcost_dtm = (double)(N_dtm - my_dtm) * d;
1246 else
1247 fcost_dtm = (double)(N_dtm - my_dtm) * c;
1248 fcost_cost = (double)(N_cost + my_cost) / 2.0;
1249 min_cost =
1250 pres_cell->min_cost + fcost_dtm + (NS_fac * a) +
1251 lambda * fcost_cost * NS_fac;
1252 break;
1253 case 4:
1254 S_dtm = costs.dtm;
1255 S_cost = costs.cost_in;
1256 if (Rast_is_d_null_value(&S_cost))
1257 continue;
1258 check_dtm = (S_dtm - my_dtm) / NS_fac;
1259 if (check_dtm >= 0)
1260 fcost_dtm = (double)(S_dtm - my_dtm) * b;
1261 else if (check_dtm < (slope_factor))
1262 fcost_dtm = (double)(S_dtm - my_dtm) * d;
1263 else
1264 fcost_dtm = (double)(S_dtm - my_dtm) * c;
1265 fcost_cost = (double)(S_cost + my_cost) / 2.0;
1266 min_cost =
1267 pres_cell->min_cost + fcost_dtm + (NS_fac * a) +
1268 lambda * fcost_cost * NS_fac;
1269 break;
1270 case 5:
1271 NW_dtm = costs.dtm;
1272 NW_cost = costs.cost_in;
1273 if (Rast_is_d_null_value(&NW_cost))
1274 continue;
1275 check_dtm = (NW_dtm - my_dtm) / DIAG_fac;
1276 if (check_dtm >= 0)
1277 fcost_dtm = (double)(NW_dtm - my_dtm) * b;
1278 else if (check_dtm < (slope_factor))
1279 fcost_dtm = (double)(NW_dtm - my_dtm) * d;
1280 else
1281 fcost_dtm = (double)(NW_dtm - my_dtm) * c;
1282 fcost_cost = (double)(NW_cost + my_cost) / 2.0;
1283 min_cost =
1284 pres_cell->min_cost + fcost_dtm + (DIAG_fac * a) +
1285 lambda * fcost_cost * DIAG_fac;
1286 break;
1287 case 6:
1288 NE_dtm = costs.dtm;
1289 NE_cost = costs.cost_in;
1290 if (Rast_is_d_null_value(&NE_cost))
1291 continue;
1292 check_dtm = (NE_dtm - my_dtm) / DIAG_fac;
1293 if (check_dtm >= 0)
1294 fcost_dtm = (double)(NE_dtm - my_dtm) * b;
1295 else if (check_dtm < (slope_factor))
1296 fcost_dtm = (double)(NE_dtm - my_dtm) * d;
1297 else
1298 fcost_dtm = (double)(NE_dtm - my_dtm) * c;
1299 fcost_cost = (double)(NE_cost + my_cost) / 2.0;
1300 min_cost =
1301 pres_cell->min_cost + fcost_dtm + (DIAG_fac * a) +
1302 lambda * fcost_cost * DIAG_fac;
1303 break;
1304 case 7:
1305 SE_dtm = costs.dtm;
1306 SE_cost = costs.cost_in;
1307 if (Rast_is_d_null_value(&SE_cost))
1308 continue;
1309 check_dtm = (SE_dtm - my_dtm) / DIAG_fac;
1310 if (check_dtm >= 0)
1311 fcost_dtm = (double)(SE_dtm - my_dtm) * b;
1312 else if (check_dtm < (slope_factor))
1313 fcost_dtm = (double)(SE_dtm - my_dtm) * d;
1314 else
1315 fcost_dtm = (double)(SE_dtm - my_dtm) * c;
1316 fcost_cost = (double)(SE_cost + my_cost) / 2.0;
1317 min_cost =
1318 pres_cell->min_cost + fcost_dtm + (DIAG_fac * a) +
1319 lambda * fcost_cost * DIAG_fac;
1320 break;
1321 case 8:
1322 SW_dtm = costs.dtm;
1323 SW_cost = costs.cost_in;
1324 if (Rast_is_d_null_value(&SW_cost))
1325 continue;
1326 check_dtm = (SW_dtm - my_dtm) / DIAG_fac;
1327 if (check_dtm >= 0)
1328 fcost_dtm = (double)(SW_dtm - my_dtm) * b;
1329 else if (check_dtm < (slope_factor))
1330 fcost_dtm = (double)(SW_dtm - my_dtm) * d;
1331 else
1332 fcost_dtm = (double)(SW_dtm - my_dtm) * c;
1333 fcost_cost = (double)(SW_cost + my_cost) / 2.0;
1334 min_cost =
1335 pres_cell->min_cost + fcost_dtm + (DIAG_fac * a) +
1336 lambda * fcost_cost * DIAG_fac;
1337 break;
1338 case 9:
1339 NNW_dtm = costs.dtm;
1340 NNW_cost = costs.cost_in;
1341 if (Rast_is_d_null_value(&NNW_cost))
1342 continue;
1343 check_dtm = (NNW_dtm - my_dtm) / V_DIAG_fac;
1344 if (check_dtm >= 0)
1345 fcost_dtm = (double)(NNW_dtm - my_dtm) * b;
1346 else if (check_dtm < (slope_factor))
1347 fcost_dtm = (double)(NNW_dtm - my_dtm) * d;
1348 else
1349 fcost_dtm = (double)(NNW_dtm - my_dtm) * c;
1350 fcost_cost =
1351 (double)(N_cost + NW_cost + NNW_cost + my_cost) / 4.0;
1352 min_cost =
1353 pres_cell->min_cost + fcost_dtm + (V_DIAG_fac * a) +
1354 lambda * fcost_cost * V_DIAG_fac;
1355 break;
1356 case 10:
1357 NNE_dtm = costs.dtm;
1358 NNE_cost = costs.cost_in;
1359 if (Rast_is_d_null_value(&NNE_cost))
1360 continue;
1361 check_dtm = ((NNE_dtm - my_dtm) / V_DIAG_fac);
1362 if (check_dtm >= 0)
1363 fcost_dtm = (double)(NNE_dtm - my_dtm) * b;
1364 else if (check_dtm < (slope_factor))
1365 fcost_dtm = (double)(NNE_dtm - my_dtm) * d;
1366 else
1367 fcost_dtm = (double)(NNE_dtm - my_dtm) * c;
1368 fcost_cost =
1369 (double)(N_cost + NE_cost + NNE_cost + my_cost) / 4.0;
1370 min_cost =
1371 pres_cell->min_cost + fcost_dtm + (V_DIAG_fac * a) +
1372 lambda * fcost_cost * V_DIAG_fac;
1373 break;
1374 case 11:
1375 SSE_dtm = costs.dtm;
1376 SSE_cost = costs.cost_in;
1377 if (Rast_is_d_null_value(&SSE_cost))
1378 continue;
1379 check_dtm = (SSE_dtm - my_dtm) / V_DIAG_fac;
1380 if (check_dtm >= 0)
1381 fcost_dtm = (double)(SSE_dtm - my_dtm) * b;
1382 else if (check_dtm < (slope_factor))
1383 fcost_dtm = (double)(SSE_dtm - my_dtm) * d;
1384 else
1385 fcost_dtm = (double)(SSE_dtm - my_dtm) * c;
1386 fcost_cost =
1387 (double)(S_cost + SE_cost + SSE_cost + my_cost) / 4.0;
1388 min_cost =
1389 pres_cell->min_cost + fcost_dtm + (V_DIAG_fac * a) +
1390 lambda * fcost_cost * V_DIAG_fac;
1391 break;
1392 case 12:
1393 SSW_dtm = costs.dtm;
1394 SSW_cost = costs.cost_in;
1395 if (Rast_is_d_null_value(&SSW_cost))
1396 continue;
1397 check_dtm = (SSW_dtm - my_dtm) / V_DIAG_fac;
1398 if (check_dtm >= 0)
1399 fcost_dtm = (double)(SSW_dtm - my_dtm) * b;
1400 else if (check_dtm < (slope_factor))
1401 fcost_dtm = (double)(SSW_dtm - my_dtm) * d;
1402 else
1403 fcost_dtm = (double)(SSW_dtm - my_dtm) * c;
1404 fcost_cost =
1405 (double)(S_cost + SW_cost + SSW_cost + my_cost) / 4.0;
1406 min_cost =
1407 pres_cell->min_cost + fcost_dtm + (V_DIAG_fac * a) +
1408 lambda * fcost_cost * V_DIAG_fac;
1409 break;
1410 case 13:
1411 WNW_dtm = costs.dtm;
1412 WNW_cost = costs.cost_in;
1413 if (Rast_is_d_null_value(&WNW_cost))
1414 continue;
1415 check_dtm = (WNW_dtm - my_dtm) / H_DIAG_fac;
1416 if (check_dtm >= 0)
1417 fcost_dtm = (double)(WNW_dtm - my_dtm) * b;
1418 else if (check_dtm < (slope_factor))
1419 fcost_dtm = (double)(WNW_dtm - my_dtm) * d;
1420 else
1421 fcost_dtm = (double)(WNW_dtm - my_dtm) * c;
1422 fcost_cost =
1423 (double)(W_cost + NW_cost + WNW_cost + my_cost) / 4.0;
1424 min_cost =
1425 pres_cell->min_cost + fcost_dtm + (H_DIAG_fac * a) +
1426 lambda * fcost_cost * H_DIAG_fac;
1427 break;
1428 case 14:
1429 ENE_dtm = costs.dtm;
1430 ENE_cost = costs.cost_in;
1431 if (Rast_is_d_null_value(&ENE_cost))
1432 continue;
1433 check_dtm = (ENE_dtm - my_dtm) / H_DIAG_fac;
1434 if (check_dtm >= 0)
1435 fcost_dtm = (double)(ENE_dtm - my_dtm) * b;
1436 else if (check_dtm < (slope_factor))
1437 fcost_dtm = (double)(ENE_dtm - my_dtm) * d;
1438 else
1439 fcost_dtm = (double)(ENE_dtm - my_dtm) * c;
1440 fcost_cost =
1441 (double)(E_cost + NE_cost + ENE_cost + my_cost) / 4.0;
1442 min_cost =
1443 pres_cell->min_cost + fcost_dtm + (H_DIAG_fac * a) +
1444 lambda * fcost_cost * H_DIAG_fac;
1445 break;
1446 case 15:
1447 ESE_dtm = costs.dtm;
1448 ESE_cost = costs.cost_in;
1449 if (Rast_is_d_null_value(&ESE_cost))
1450 continue;
1451 check_dtm = (ESE_dtm - my_dtm) / H_DIAG_fac;
1452 if (check_dtm >= 0)
1453 fcost_dtm = (double)(ESE_dtm - my_dtm) * b;
1454 else if (check_dtm < (slope_factor))
1455 fcost_dtm = (double)(ESE_dtm - my_dtm) * d;
1456 else
1457 fcost_dtm = (double)(ESE_dtm - my_dtm) * c;
1458 fcost_cost =
1459 (double)(E_cost + SE_cost + ESE_cost + my_cost) / 4.0;
1460 min_cost =
1461 pres_cell->min_cost + fcost_dtm + (H_DIAG_fac * a) +
1462 lambda * fcost_cost * H_DIAG_fac;
1463 break;
1464 case 16:
1465 WSW_dtm = costs.dtm;
1466 WSW_cost = costs.cost_in;
1467 if (Rast_is_d_null_value(&WSW_cost))
1468 continue;
1469 check_dtm = (WSW_dtm - my_dtm) / H_DIAG_fac;
1470 if (check_dtm >= 0)
1471 fcost_dtm = (double)(WSW_dtm - my_dtm) * b;
1472 else if (check_dtm < (slope_factor))
1473 fcost_dtm = (double)(WSW_dtm - my_dtm) * d;
1474 else
1475 fcost_dtm = (double)(WSW_dtm - my_dtm) * c;
1476 fcost_cost =
1477 (double)(W_cost + SW_cost + WSW_cost + my_cost) / 4.0;
1478 min_cost =
1479 pres_cell->min_cost + fcost_dtm + (H_DIAG_fac * a) +
1480 lambda * fcost_cost * H_DIAG_fac;
1481 break;
1482 }
1483
1484 /* skip if costs could not be calculated */
1485 if (Rast_is_d_null_value(&min_cost))
1486 continue;
1487
1488 if (Segment_get(&cost_seg, &costs, row, col) < 0)
1489 G_fatal_error(_("Can not read from temporary file"));
1490 old_min_cost = costs.cost_out;
1491
1492 /* add to list */
1493 if (Rast_is_d_null_value(&old_min_cost)) {
1494 costs.cost_out = min_cost;
1495 costs.nearest = nearest;
1496 if (Segment_put(&cost_seg, &costs, row, col) < 0)
1497 G_fatal_error(_("Can not write to temporary file"));
1498 insert(min_cost, row, col);
1499 if (dir == 1) {
1500 if (dir_bin)
1501 cur_dir = (1 << (int)cur_dir);
1502 if (Segment_put(&dir_seg, &cur_dir, row, col) < 0)
1503 G_fatal_error(_("Can not write to temporary file"));
1504 }
1505 if (have_solver) {
1506 if (Segment_get(&solve_seg, solvedir, row, col) < 0)
1507 G_fatal_error(_("Can not read from temporary file"));
1508 solvedir[1] = mysolvedir[0];
1509 if (Segment_put(&solve_seg, solvedir, row, col) < 0)
1510 G_fatal_error(_("Can not write to temporary file"));
1511 }
1512 }
1513 /* update with lower costs */
1514 else if (old_min_cost > min_cost) {
1515 costs.cost_out = min_cost;
1516 costs.nearest = nearest;
1517 if (Segment_put(&cost_seg, &costs, row, col) < 0)
1518 G_fatal_error(_("Can not write to temporary file"));
1519 insert(min_cost, row, col);
1520 if (dir == 1) {
1521 if (dir_bin)
1522 cur_dir = (1 << (int)cur_dir);
1523 if (Segment_put(&dir_seg, &cur_dir, row, col) < 0)
1524 G_fatal_error(_("Can not write to temporary file"));
1525 }
1526 if (have_solver) {
1527 if (Segment_get(&solve_seg, solvedir, row, col) < 0)
1528 G_fatal_error(_("Can not read from temporary file"));
1529 solvedir[1] = mysolvedir[0];
1530 if (Segment_put(&solve_seg, solvedir, row, col) < 0)
1531 G_fatal_error(_("Can not write to temporary file"));
1532 }
1533 }
1534 else if (old_min_cost == min_cost &&
1535 (dir_bin || have_solver) &&
1536 !(FLAG_GET(visited, row, col))) {
1537 FCELL old_dir;
1538 int dir_inv[16] = { 4, 5, 6, 7, 0, 1, 2, 3,
1539 12, 13, 14, 15, 8, 9, 10, 11 };
1540 int dir_fwd;
1541 int equal = 1;
1542
1543 /* only update neighbors that have not yet been processed,
1544 * otherwise we might get circular paths */
1545
1546 if (have_solver) {
1547 if (Segment_get(&solve_seg, solvedir, row, col) < 0)
1548 G_fatal_error(_("Can not read from temporary file"));
1549 equal = (solvedir[1] == mysolvedir[0]);
1550 if (solvedir[1] > mysolvedir[0]) {
1551 solvedir[1] = mysolvedir[0];
1552 if (Segment_put(&solve_seg, solvedir, row, col) < 0)
1553 G_fatal_error(_("Can not write to temporary file"));
1554
1555 costs.nearest = nearest;
1556 if (Segment_put(&cost_seg, &costs, row, col) < 0)
1557 G_fatal_error(_("Can not write to temporary file"));
1558
1559 if (dir == 1) {
1560 if (dir_bin)
1561 cur_dir = (1 << (int)cur_dir);
1562 if (Segment_put(&dir_seg, &cur_dir, row, col) < 0)
1563 G_fatal_error(_("Can not write to temporary file"));
1564 }
1565 }
1566 }
1567
1568 if (dir_bin && equal) {
1569 /* this can create circular paths:
1570 * set only if current cell does not point to neighbor
1571 * does not avoid longer circular paths */
1572 if (Segment_get(&dir_seg, &old_dir, pres_cell->row, pres_cell->col) < 0)
1573 G_fatal_error(_("Can not read from temporary file"));
1574 dir_fwd = (1 << dir_inv[(int)cur_dir]);
1575 if (!((int)old_dir & dir_fwd)) {
1576 if (Segment_get(&dir_seg, &old_dir, row, col) < 0)
1577 G_fatal_error(_("Can not read from temporary file"));
1578 cur_dir = ((1 << (int)cur_dir) | (int)old_dir);
1579 if (Segment_put(&dir_seg, &cur_dir, row, col) < 0)
1580 G_fatal_error(_("Can not write to temporary file"));
1581 }
1582 }
1583 }
1584 }
1585
1586 if (stop_pnts && time_to_stop(pres_cell->row, pres_cell->col))
1587 break;
1588
1589 ct = pres_cell;
1590 delete(pres_cell);
1591 pres_cell = get_lowest();
1592
1593 if (ct == pres_cell)
1594 G_warning(_("Error, ct == pres_cell"));
1595 }
1596 G_percent(1, 1, 1);
1597
1598 /* free heap */
1599 free_heap();
1600 flag_destroy(visited);
1601
1602 if (have_solver) {
1603 Segment_close(&solve_seg);
1604 }
1605
1606 /* Open cumulative cost layer for writing */
1607 cum_fd = Rast_open_new(cum_cost_layer, cum_data_type);
1608 cum_cell = Rast_allocate_buf(cum_data_type);
1609
1610 /* Open nearest start point layer */
1611 if (nearest_layer) {
1612 nearest_fd = Rast_open_new(nearest_layer, nearest_data_type);
1613 nearest_cell = Rast_allocate_buf(nearest_data_type);
1614 }
1615 else {
1616 nearest_fd = -1;
1617 nearest_cell = NULL;
1618 }
1619 nearest_size = Rast_cell_size(nearest_data_type);
1620
1621 /* Copy segmented map to output map */
1622 G_message(_("Writing output raster map <%s>... "), cum_cost_layer);
1623 if (nearest_layer) {
1624 G_message(_("Writing raster map with nearest start point <%s>..."), nearest_layer);
1625 }
1626
1627 cell2 = Rast_allocate_buf(dtm_data_type);
1628 {
1629 void *p;
1630 void *p2;
1631 void *p3;
1632 int cum_dsize = Rast_cell_size(cum_data_type);
1633
1634 Rast_set_null_value(cell2, ncols, dtm_data_type);
1635
1636 for (row = 0; row < nrows; row++) {
1637 G_percent(row, nrows, 2);
1638 if (keep_nulls)
1639 Rast_get_row(dtm_fd, cell2, row, dtm_data_type);
1640
1641 p = cum_cell;
1642 p2 = cell2;
1643 p3 = nearest_cell;
1644 for (col = 0; col < ncols; col++) {
1645 if (keep_nulls) {
1646 if (Rast_is_null_value(p2, dtm_data_type)) {
1647 Rast_set_null_value(p, 1, cum_data_type);
1648 p = G_incr_void_ptr(p, cum_dsize);
1649 p2 = G_incr_void_ptr(p2, dtm_dsize);
1650 if (nearest_layer) {
1651 Rast_set_null_value(p3, 1, nearest_data_type);
1652 p3 = G_incr_void_ptr(p3, nearest_size);
1653 }
1654
1655 continue;
1656 }
1657 }
1658 if (Segment_get(&cost_seg, &costs, row, col) < 0)
1659 G_fatal_error(_("Can not read from temporary file"));
1660 min_cost = costs.cost_out;
1661 nearest = costs.nearest;
1662 if (Rast_is_d_null_value(&min_cost)) {
1663 Rast_set_null_value((p), 1, cum_data_type);
1664 if (nearest_layer)
1665 Rast_set_null_value(p3, 1, nearest_data_type);
1666 }
1667 else {
1668 if (min_cost > peak)
1669 peak = min_cost;
1670
1671 switch (cum_data_type) {
1672 case CELL_TYPE:
1673 *(CELL *)p = (CELL)(min_cost + .5);
1674 break;
1675 case FCELL_TYPE:
1676 *(FCELL *)p = (FCELL)(min_cost);
1677 break;
1678 case DCELL_TYPE:
1679 *(DCELL *)p = (DCELL)(min_cost);
1680 break;
1681 }
1682
1683 if (nearest_layer) {
1684 switch (nearest_data_type) {
1685 case CELL_TYPE:
1686 *(CELL *)p3 = (CELL)(nearest);
1687 break;
1688 case FCELL_TYPE:
1689 *(FCELL *)p3 = (FCELL)(nearest);
1690 break;
1691 case DCELL_TYPE:
1692 *(DCELL *)p3 = (DCELL)(nearest);
1693 break;
1694 }
1695 }
1696 }
1697 p = G_incr_void_ptr(p, cum_dsize);
1698 p2 = G_incr_void_ptr(p2, dtm_dsize);
1699 if (nearest_layer)
1700 p3 = G_incr_void_ptr(p3, nearest_size);
1701 }
1702 Rast_put_row(cum_fd, cum_cell, cum_data_type);
1703 if (nearest_layer)
1704 Rast_put_row(nearest_fd, nearest_cell, nearest_data_type);
1705 }
1706 G_percent(1, 1, 1);
1707 G_free(cum_cell);
1708 G_free(cell2);
1709 if (nearest_layer)
1710 G_free(nearest_cell);
1711 }
1712
1713 if (dir == 1) {
1714 void *p;
1715 size_t dir_size = Rast_cell_size(dir_data_type);
1716
1717 dir_fd = Rast_open_new(move_dir_layer, dir_data_type);
1718 dir_cell = Rast_allocate_buf(dir_data_type);
1719
1720 G_message(_("Writing output movement direction raster map <%s>..."), move_dir_layer);
1721 for (row = 0; row < nrows; row++) {
1722 p = dir_cell;
1723 for (col = 0; col < ncols; col++) {
1724 if (Segment_get(&dir_seg, &cur_dir, row, col) < 0)
1725 G_fatal_error(_("Can not read from temporary file"));
1726 *((FCELL *) p) = cur_dir;
1727 p = G_incr_void_ptr(p, dir_size);
1728 }
1729 Rast_put_row(dir_fd, dir_cell, dir_data_type);
1730 G_percent(row, nrows, 2);
1731 }
1732 G_percent(1, 1, 1);
1733 G_free(dir_cell);
1734 }
1735
1736 Segment_close(&cost_seg); /* release memory */
1737 if (dir == 1)
1738 Segment_close(&dir_seg);
1739
1740 Rast_close(dtm_fd);
1741 Rast_close(cost_fd);
1742 Rast_close(cum_fd);
1743 if (dir == 1)
1744 Rast_close(dir_fd);
1745 if (nearest_layer)
1746 Rast_close(nearest_fd);
1747
1748 /* writing history file */
1749 Rast_short_history(cum_cost_layer, "raster", &history);
1750 Rast_command_history(&history);
1751 Rast_write_history(cum_cost_layer, &history);
1752
1753 if (dir == 1) {
1754 Rast_short_history(move_dir_layer, "raster", &history);
1755 Rast_command_history(&history);
1756 Rast_write_history(move_dir_layer, &history);
1757 }
1758
1759 if (nearest_layer) {
1760 Rast_short_history(nearest_layer, "raster", &history);
1761 Rast_command_history(&history);
1762 Rast_write_history(nearest_layer, &history);
1763 if (opt9->answer) {
1764 struct Colors colors;
1765 Rast_read_colors(opt9->answer, "", &colors);
1766 Rast_write_colors(nearest_layer, G_mapset(), &colors);
1767 }
1768 else {
1769 struct Colors colors;
1770 struct Range range;
1771 CELL min, max;
1772
1773 Rast_read_range(nearest_layer, G_mapset(), &range);
1774 Rast_get_range_min_max(&range, &min, &max);
1775 Rast_make_random_colors(&colors, min, max);
1776 Rast_write_colors(nearest_layer, G_mapset(), &colors);
1777 }
1778 }
1779
1780 /* Create colours for output map */
1781
1782 /*
1783 * Rast_read_range (cum_cost_layer, "", &range);
1784 * Rast_get_range_min_max(&range, &min, &max);
1785 * G_make_color_wave(&colors,min, max);
1786 * Rast_write_colors (cum_cost_layer,"",&colors);
1787 */
1788
1789 G_done_msg(_("Peak cost value: %g"), peak);
1790
1791 exit(EXIT_SUCCESS);
1792 }
1793
1794 struct start_pt *
process_start_coords(char ** answers,struct start_pt * top_start_pt)1795 process_start_coords(char **answers, struct start_pt *top_start_pt)
1796 {
1797 int col, row;
1798 double east, north;
1799 struct start_pt *new_start_pt;
1800 int point_no = 0;
1801
1802 if (!answers)
1803 return (0);
1804
1805 for (; *answers != NULL; answers += 2) {
1806 if (!G_scan_easting(*answers, &east, G_projection()))
1807 G_fatal_error(_("Illegal x coordinate <%s>"), *answers);
1808 if (!G_scan_northing(*(answers + 1), &north, G_projection()))
1809 G_fatal_error(_("Illegal y coordinate <%s>"), *(answers + 1));
1810
1811 if (east < window.west || east > window.east ||
1812 north < window.south || north > window.north) {
1813 G_warning(_("Warning, ignoring point outside window: %g, %g"),
1814 east, north);
1815 continue;
1816 }
1817
1818 row = (window.north - north) / window.ns_res;
1819 col = (east - window.west) / window.ew_res;
1820
1821 new_start_pt = (struct start_pt *)(G_malloc(sizeof(struct start_pt)));
1822
1823 new_start_pt->row = row;
1824 new_start_pt->col = col;
1825 new_start_pt->value = ++point_no;
1826 new_start_pt->next = top_start_pt;
1827 top_start_pt = new_start_pt;
1828 }
1829
1830 return top_start_pt;
1831 }
1832
process_stop_coords(char ** answers)1833 int process_stop_coords(char **answers)
1834 {
1835 int col, row;
1836 double east, north;
1837
1838 if (!answers)
1839 return 0;
1840
1841 for (; *answers != NULL; answers += 2) {
1842 if (!G_scan_easting(*answers, &east, G_projection()))
1843 G_fatal_error(_("Illegal x coordinate <%s>"), *answers);
1844 if (!G_scan_northing(*(answers + 1), &north, G_projection()))
1845 G_fatal_error(_("Illegal y coordinate <%s>"), *(answers + 1));
1846
1847 if (east < window.west || east > window.east ||
1848 north < window.south || north > window.north) {
1849 G_warning(_("Warning, ignoring point outside window: %g, %g"),
1850 east, north);
1851 continue;
1852 }
1853
1854 row = (window.north - north) / window.ns_res;
1855 col = (east - window.west) / window.ew_res;
1856
1857 add_stop_pnt(row, col);
1858 }
1859
1860 return (stop_pnts != NULL);
1861 }
1862
add_stop_pnt(int r,int c)1863 void add_stop_pnt(int r, int c)
1864 {
1865 int i;
1866 struct rc sp;
1867
1868 if (n_stop_pnts == stop_pnts_alloc) {
1869 stop_pnts_alloc += 100;
1870 stop_pnts = (struct rc *)G_realloc(stop_pnts, stop_pnts_alloc * sizeof(struct rc));
1871 }
1872
1873 sp.r = r;
1874 sp.c = c;
1875 i = n_stop_pnts;
1876 while (i > 0 && cmp_rc(stop_pnts + i - 1, &sp) > 0) {
1877 stop_pnts[i] = stop_pnts[i - 1];
1878 i--;
1879 }
1880 stop_pnts[i] = sp;
1881
1882 n_stop_pnts++;
1883 }
1884
time_to_stop(int row,int col)1885 int time_to_stop(int row, int col)
1886 {
1887 int lo, mid, hi;
1888 struct rc sp;
1889 static int hits = 0;
1890
1891 sp.r = row;
1892 sp.c = col;
1893
1894 lo = 0;
1895 hi = n_stop_pnts - 1;
1896
1897 /* bsearch with deferred test for equality
1898 * slightly more efficient for worst case: no match */
1899 while (lo < hi) {
1900 mid = lo + ((hi - lo) >> 1);
1901 if (cmp_rc(stop_pnts + mid, &sp) < 0)
1902 lo = mid + 1;
1903 else
1904 hi = mid;
1905 }
1906 if (cmp_rc(stop_pnts + lo, &sp) == 0) {
1907 return (++hits == n_stop_pnts);
1908 }
1909
1910 return 0;
1911 }
1912