source: grass/trunk/raster/r.drain/main.c

Last change on this file was 72231, checked in by mmetz, 7 years ago

r.path/r.drain: accept multiple start coordinates
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id
  • Property svn:mime-type set to text/x-csrc
File size: 22.1 KB
Line 
1/****************************************************************************
2 *
3 * MODULE: r.drain
4 *
5 * AUTHOR(S): Kewan Q. Khawaja <kewan techlogix.com>
6 * Update to FP (2000): Pierre de Mouveaux <pmx@audiovu.com> <pmx@free.fr>
7 * bugfix in FCELL, DCELL: Markus Neteler 12/2000
8 * Rewritten by Roger Miller 7/2001 based on subroutines from
9 * r.fill.dir and on the original r.drain.
10 * 24 July 2004: WebValley 2004, error checking and vector points added by
11 * Matteo Franchi Liceo Leonardo Da Vinci Trento
12 * Roberto Flor ITC-irst
13 * New code added by Colin Nielsen <colin.nielsen at gmail dot com> *
14 * to use movement direction surface from r.walk and r.cost, and to
15 * output vector paths 2/2009
16 *
17 * PURPOSE: This is the main program for tracing out the path that a
18 * drop of water would take if released at a certain location
19 * on an input elevation map.
20 * COPYRIGHT: (C) 2000,2009 by the GRASS Development Team
21 *
22 * This program is free software under the GNU General Public
23 * License (>=v2). Read the file COPYING that comes with GRASS
24 * for details.
25 *
26 *****************************************************************************/
27
28#include <grass/config.h>
29#include <stdlib.h>
30#include <stdio.h>
31#include <string.h>
32#include <math.h>
33#include <float.h>
34
35/* for using the "open" statement */
36#include <sys/types.h>
37#include <sys/stat.h>
38#include <fcntl.h>
39
40/* for using the close statement */
41#include <unistd.h>
42
43#include <grass/gis.h>
44#include <grass/raster.h>
45#include <grass/glocale.h>
46#include <grass/vector.h>
47
48#define DEBUG
49#include "tinf.h"
50#include "local.h"
51
52/* should probably be updated to a pointer array & malloc/realloc as needed */
53#define POINTS_INCREMENT 1024
54
55/* define a data structure to hold the point data */
56struct point
57{
58 int row;
59 int col;
60 struct point *next;
61 double value;
62};
63
64int main(int argc, char **argv)
65{
66
67 int fe, fd, dir_fd;
68 int i, have_points = 0;
69 int new_id;
70 int nrows, ncols;
71 int *points_row = NULL, *points_col = NULL, npoints;
72 int increment_count;
73 int cell_open(), cell_open_new();
74 int map_id, dir_id;
75 char map_name[GNAME_MAX], new_map_name[GNAME_MAX], dir_name[GNAME_MAX];
76 char *tempfile1, *tempfile2, *tempfile3;
77 struct History history;
78
79 struct Cell_head window;
80 struct Option *opt1, *opt2, *coordopt, *vpointopt, *opt3, *opt4;
81 struct Flag *flag1, *flag2, *flag3, *flag4;
82 struct GModule *module;
83 int in_type;
84 void *in_buf;
85 void *dir_buf;
86 CELL *out_buf;
87 struct band3 bnd, bndC;
88 struct metrics *m = NULL;
89
90 struct point *list;
91 struct point *thispoint;
92 int ival, start_row, start_col, mode;
93 off_t bsz;
94 int costmode = 0;
95 double east, north, val;
96 struct point *drain(int, struct point *, int, int);
97 struct point *drain_cost(int, struct point *, int, int);
98 int bsort(int, struct point *);
99
100 struct line_pnts *Points;
101 struct line_cats *Cats;
102 struct Map_info vout;
103 int cat;
104 double x, y;
105
106 G_gisinit(argv[0]);
107
108 module = G_define_module();
109 G_add_keyword(_("raster"));
110 G_add_keyword(_("hydrology"));
111 G_add_keyword(_("cost surface"));
112 module->description =
113 _("Traces a flow through an elevation model or cost surface on a raster map.");
114
115 opt1 = G_define_standard_option(G_OPT_R_INPUT);
116 opt1->description = _("Name of input elevation or cost surface raster map");
117
118 opt3 = G_define_standard_option(G_OPT_R_INPUT);
119 opt3->key = "direction";
120 opt3->label =
121 _("Name of input movement direction map associated with the cost surface");
122 opt3->description =
123 _("Direction in degrees CCW from east");
124 opt3->required = NO;
125 opt3->guisection = _("Cost surface");
126
127 opt2 = G_define_standard_option(G_OPT_R_OUTPUT);
128
129 opt4 = G_define_standard_option(G_OPT_V_OUTPUT);
130 opt4->key = "drain";
131 opt4->required = NO;
132 opt4->label =
133 _("Name for output drain vector map");
134 opt4->description = _("Recommended for cost surface made using knight's move");
135
136 coordopt = G_define_standard_option(G_OPT_M_COORDS);
137 coordopt->key = "start_coordinates";
138 coordopt->multiple = YES;
139 coordopt->description = _("Coordinates of starting point(s) (E,N)");
140 coordopt->guisection = _("Start");
141
142 vpointopt = G_define_standard_option(G_OPT_V_INPUTS);
143 vpointopt->key = "start_points";
144 vpointopt->required = NO;
145 vpointopt->label = _("Name of starting vector points map(s)");
146 vpointopt->description = NULL;
147 vpointopt->guisection = _("Start");
148
149 flag1 = G_define_flag();
150 flag1->key = 'c';
151 flag1->description = _("Copy input cell values on output");
152 flag1->guisection = _("Path settings");
153
154 flag2 = G_define_flag();
155 flag2->key = 'a';
156 flag2->description = _("Accumulate input values along the path");
157 flag2->guisection = _("Path settings");
158
159 flag3 = G_define_flag();
160 flag3->key = 'n';
161 flag3->description = _("Count cell numbers along the path");
162 flag3->guisection = _("Path settings");
163
164 flag4 = G_define_flag();
165 flag4->key = 'd';
166 flag4->description =
167 _("The input raster map is a cost surface (direction surface must also be specified)");
168 flag4->guisection = _("Cost surface");
169
170 if (G_parser(argc, argv))
171 exit(EXIT_FAILURE);
172
173
174 strcpy(map_name, opt1->answer);
175 strcpy(new_map_name, opt2->answer);
176
177 if (flag4->answer) {
178 costmode = 1;
179 G_verbose_message(_("Directional drain selected... checking for direction raster map"));
180 }
181 else {
182 G_verbose_message(_("Surface/Hydrology drain selected"));
183 }
184
185 if (costmode == 1) {
186 if (!opt3->answer) {
187 G_fatal_error(_("Direction raster map <%s> not specified, if direction flag is on, "
188 "a direction raster must be given"), opt3->key);
189 }
190 strcpy(dir_name, opt3->answer);
191 }
192 if (costmode == 0) {
193 if (opt3->answer) {
194 G_fatal_error(_("Direction raster map <%s> should not be specified for Surface/Hydrology drains"),
195 opt3->answer);
196 }
197 }
198
199 if (opt4->answer) {
200 if (0 > Vect_open_new(&vout, opt4->answer, 0)) {
201 G_fatal_error(_("Unable to create vector map <%s>"),
202 opt4->answer);
203 }
204 Vect_hist_command(&vout);
205 }
206 /* allocate cell buf for the map layer */
207 in_type = Rast_map_type(map_name, "");
208
209 /* set the pointers for multi-typed functions */
210 set_func_pointers(in_type);
211
212 if ((flag1->answer + flag2->answer + flag3->answer) > 1)
213 G_fatal_error(_("Specify just one of the -c, -a and -n flags"));
214
215 mode = 0;
216 if (flag1->answer)
217 mode = 1;
218 if (flag2->answer)
219 mode = 2;
220 if (flag3->answer)
221 mode = 3;
222
223 /* get the window information */
224 G_get_window(&window);
225 nrows = Rast_window_rows();
226 ncols = Rast_window_cols();
227
228 /* calculate true cell resolution */
229 m = (struct metrics *)G_malloc(nrows * sizeof(struct metrics));
230 points_row = (int*)G_calloc(POINTS_INCREMENT, sizeof(int));
231 points_col = (int*)G_calloc(POINTS_INCREMENT, sizeof(int));
232
233 increment_count = 1;
234 npoints = 0;
235 if (coordopt->answer) {
236 for (i = 0; coordopt->answers[i] != NULL; i += 2) {
237 G_scan_easting(coordopt->answers[i], &east, G_projection());
238 G_scan_northing(coordopt->answers[i + 1], &north, G_projection());
239 start_col = (int)Rast_easting_to_col(east, &window);
240 start_row = (int)Rast_northing_to_row(north, &window);
241
242 if (start_row < 0 || start_row > nrows ||
243 start_col < 0 || start_col > ncols) {
244 G_warning(_("Starting point %d is outside the current region"),
245 i + 1);
246 continue;
247 }
248 points_row[npoints] = start_row;
249 points_col[npoints] = start_col;
250 npoints++;
251 if (npoints == POINTS_INCREMENT * increment_count)
252 {
253 increment_count++;
254 points_row = (int*)G_realloc(points_row, POINTS_INCREMENT * increment_count * sizeof(int));
255 points_col = (int*)G_realloc(points_col, POINTS_INCREMENT * increment_count * sizeof(int));
256 }
257 have_points = 1;
258 }
259 }
260 if (vpointopt->answers) {
261 for (i = 0; vpointopt->answers[i] != NULL; i++) {
262 struct Map_info In;
263 struct bound_box box;
264 int type;
265
266 Points = Vect_new_line_struct();
267 Cats = Vect_new_cats_struct();
268
269 Vect_set_open_level(1); /* topology not required */
270
271 if (1 > Vect_open_old(&In, vpointopt->answers[i], ""))
272 G_fatal_error(_("Unable to open vector map <%s>"), vpointopt->answers[i]);
273
274 G_message(_("Reading vector map <%s> with start points..."),
275 Vect_get_full_name(&In));
276
277 Vect_rewind(&In);
278
279 Vect_region_box(&window, &box);
280
281 while (1) {
282 /* register line */
283 type = Vect_read_next_line(&In, Points, Cats);
284
285 /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
286 if (type == -1) {
287 G_warning(_("Unable to read vector map"));
288 continue;
289 }
290 else if (type == -2) {
291 break;
292 }
293 if (!Vect_point_in_box(Points->x[0], Points->y[0], 0, &box))
294 continue;
295
296 start_col = (int)Rast_easting_to_col(Points->x[0], &window);
297 start_row = (int)Rast_northing_to_row(Points->y[0], &window);
298
299 /* effectively just a duplicate check to G_site_in_region() ??? */
300 if (start_row < 0 || start_row > nrows || start_col < 0 ||
301 start_col > ncols)
302 continue;
303
304 points_row[npoints] = start_row;
305 points_col[npoints] = start_col;
306 npoints++;
307 if (npoints == POINTS_INCREMENT * increment_count)
308 {
309 increment_count++;
310 points_row = (int*)G_realloc(points_row, POINTS_INCREMENT * increment_count * sizeof(int));
311 points_col = (int*)G_realloc(points_col, POINTS_INCREMENT * increment_count * sizeof(int));
312 }
313 have_points = 1;
314 }
315 Vect_close(&In);
316
317 /* only catches maps out of range until something is found, not after */
318 if (!have_points) {
319 G_warning(_("Starting vector map <%s> contains no points in the current region"),
320 vpointopt->answers[i]);
321 }
322 Vect_destroy_line_struct(Points);
323 Vect_destroy_cats_struct(Cats);
324 }
325 }
326 if (have_points == 0)
327 G_fatal_error(_("No start/stop point(s) specified"));
328
329 /* determine the drainage paths */
330
331 /* allocate storage for the first point */
332 thispoint = (struct point *)G_malloc(sizeof(struct point));
333 list = thispoint;
334 thispoint->next = NULL;
335
336 G_begin_distance_calculations();
337 {
338 double e1, n1, e2, n2;
339
340 e1 = window.east;
341 n1 = window.north;
342 e2 = e1 + window.ew_res;
343 n2 = n1 - window.ns_res;
344 for (i = 0; i < nrows; i++) {
345 m[i].ew_res = G_distance(e1, n1, e2, n1);
346 m[i].ns_res = G_distance(e1, n1, e1, n2);
347 m[i].diag_res = G_distance(e1, n1, e2, n2);
348 e2 = e1 + window.ew_res;
349 n2 = n1 - window.ns_res;
350 }
351 }
352
353 /* buffers for internal use */
354 bndC.ns = ncols;
355 bndC.sz = sizeof(CELL) * ncols;
356 bndC.b[0] = G_calloc(ncols, sizeof(CELL));
357 bndC.b[1] = G_calloc(ncols, sizeof(CELL));
358 bndC.b[2] = G_calloc(ncols, sizeof(CELL));
359
360 /* buffers for external use */
361 bnd.ns = ncols;
362 bnd.sz = ncols * bpe();
363 bnd.b[0] = G_calloc(ncols, bpe());
364 bnd.b[1] = G_calloc(ncols, bpe());
365 bnd.b[2] = G_calloc(ncols, bpe());
366
367 /* an input buffer */
368 in_buf = get_buf();
369
370 /* open the original map and get its file id */
371 map_id = Rast_open_old(map_name, "");
372
373 /* get some temp files */
374 tempfile1 = G_tempfile();
375 tempfile2 = G_tempfile();
376
377 fe = open(tempfile1, O_RDWR | O_CREAT, 0666);
378 fd = open(tempfile2, O_RDWR | O_CREAT, 0666);
379
380 /* transfer the input map to a temp file */
381 for (i = 0; i < nrows; i++) {
382 get_row(map_id, in_buf, i);
383 write(fe, in_buf, bnd.sz);
384 }
385 Rast_close(map_id);
386
387 if (costmode == 1) {
388 dir_buf = Rast_allocate_d_buf();
389 dir_id = Rast_open_old(dir_name, "");
390 tempfile3 = G_tempfile();
391 dir_fd = open(tempfile3, O_RDWR | O_CREAT, 0666);
392
393 for (i = 0; i < nrows; i++) {
394 Rast_get_d_row(dir_id, dir_buf, i);
395 write(dir_fd, dir_buf, ncols * sizeof(DCELL));
396 }
397 Rast_close(dir_id);
398 }
399
400 /* only necessary for non-dir drain */
401 if (costmode == 0) {
402 G_message(_("Calculating flow directions..."));
403
404 /* fill one-cell pits and take a first stab at flow directions */
405 filldir(fe, fd, nrows, &bnd, m);
406
407 /* determine flow directions for more ambiguous cases */
408 resolve(fd, nrows, &bndC);
409 }
410
411 /* free the buffers already used */
412 G_free(bndC.b[0]);
413 G_free(bndC.b[1]);
414 G_free(bndC.b[2]);
415
416 G_free(bnd.b[0]);
417 G_free(bnd.b[1]);
418 G_free(bnd.b[2]);
419
420 /* determine the drainage paths */
421
422 /* repeat for each starting point */
423 for (i = 0; i < npoints; i++) {
424 /* use the flow directions to determine the drainage path
425 * results are compiled as a linked list of points in downstream order */
426 thispoint->row = points_row[i];
427 thispoint->col = points_col[i];
428 thispoint->next = NULL;
429 /* drain algorithm selection (dir or non-dir) */
430 if (costmode == 0) {
431 thispoint = drain(fd, thispoint, nrows, ncols);
432 }
433 if (costmode == 1) {
434 thispoint = drain_cost(dir_fd, thispoint, nrows, ncols);
435 }
436 }
437
438 /* do the output */
439
440 if (mode == 0 || mode == 3) {
441
442 /* Output will be a cell map */
443 /* open a new file and allocate an output buffer */
444 new_id = Rast_open_c_new(new_map_name);
445 out_buf = Rast_allocate_c_buf();
446
447 /* mark each cell */
448 thispoint = list;
449 while (thispoint->next != NULL) {
450 thispoint->value = 1;
451 thispoint = thispoint->next;
452 }
453
454 if (mode == 3) {
455 /* number each cell downstream */
456 thispoint = list;
457 ival = 0;
458 while (thispoint->next != NULL) {
459 if (thispoint->row == INT_MAX) {
460 ival = 0;
461 thispoint = thispoint->next;
462 continue;
463 }
464 thispoint->value += ival;
465 ival = thispoint->value;
466 thispoint = thispoint->next;
467 }
468 }
469
470 /* build the output map */
471 G_message(_("Writing output raster map..."));
472 for (i = 0; i < nrows; i++) {
473 G_percent(i, nrows, 2);
474 Rast_set_c_null_value(out_buf, ncols);
475 thispoint = list;
476 while (thispoint->next != NULL) {
477 if (thispoint->row == i)
478 out_buf[thispoint->col] = (int)thispoint->value;
479 thispoint = thispoint->next;
480 }
481 Rast_put_c_row(new_id, out_buf);
482 }
483 G_percent(1, 1, 1);
484 }
485 else { /* mode = 1 or 2 */
486 /* Output will be of the same type as input */
487 /* open a new file and allocate an output buffer */
488 new_id = Rast_open_new(new_map_name, in_type);
489 out_buf = get_buf();
490 bsz = ncols * bpe();
491
492 /* loop through each point in the list and store the map values */
493 thispoint = list;
494 while (thispoint->next != NULL) {
495 if (thispoint->row == INT_MAX) {
496 thispoint = thispoint->next;
497 continue;
498 }
499 lseek(fe, (off_t) thispoint->row * bsz, SEEK_SET);
500 read(fe, in_buf, bsz);
501 memcpy(&thispoint->value, (char *)in_buf + bpe() * thispoint->col,
502 bpe());
503 thispoint = thispoint->next;
504 }
505
506 if (mode == 2) {
507 /* accumulate the input map values downstream */
508 thispoint = list;
509 val = 0.;
510 while (thispoint->next != NULL) {
511 if (thispoint->row == INT_MAX) {
512 val = 0.;
513 thispoint = thispoint->next;
514 continue;
515 }
516 sum(&thispoint->value, &val);
517 memcpy(&val, &thispoint->value, bpe());
518 thispoint = thispoint->next;
519 }
520 }
521
522 /* build the output map */
523 G_message(_("Writing raster map <%s>..."),
524 new_map_name);
525 for (i = 0; i < nrows; i++) {
526 G_percent(i, nrows, 2);
527 set_null_value(out_buf, ncols);
528 thispoint = list;
529 while (thispoint->next != NULL) {
530 if (thispoint->row == i)
531 memcpy((char *)out_buf + bpe() * thispoint->col,
532 &(thispoint->value), bpe());
533 thispoint = thispoint->next;
534 }
535 put_row(new_id, out_buf);
536 }
537 G_percent(1, 1, 1);
538 }
539
540 /* Output a vector path */
541 if (opt4->answer) {
542 Points = Vect_new_line_struct();
543 Cats = Vect_new_cats_struct();
544 /* Need to modify for multiple paths */
545 thispoint = list;
546 i = 1;
547 while (thispoint->next != NULL) {
548 if (thispoint->row == INT_MAX) {
549 thispoint = thispoint->next;
550 Vect_cat_set(Cats, 1, i);
551 Vect_write_line(&vout, GV_LINE, Points, Cats);
552 Vect_reset_line(Points);
553 Vect_reset_cats(Cats);
554 i++;
555 continue;
556 }
557 if (Vect_cat_get(Cats, 1, &cat) == 0) {
558 Vect_cat_set(Cats, 1, i);
559 }
560 /* Need to convert row and col to coordinates
561 * y = cell_head.north - ((double) p->row + 0.5) * cell_head.ns_res;
562 * x = cell_head.west + ((double) p->col + 0.5) * cell_head.ew_res;
563 */
564
565 x = window.west + ((double)thispoint->col + 0.5) * window.ew_res;
566 y = window.north - ((double)thispoint->row + 0.5) * window.ns_res;
567 Vect_append_point(Points, x, y, 0.0);
568 thispoint = thispoint->next;
569 } /* End while */
570 Vect_build(&vout);
571 Vect_close(&vout);
572 }
573
574 /* close files and free buffers */
575 Rast_close(new_id);
576
577 Rast_put_cell_title(new_map_name, "Surface flow trace");
578
579 Rast_short_history(new_map_name, "raster", &history);
580 Rast_command_history(&history);
581 Rast_write_history(new_map_name, &history);
582
583 close(fe);
584 close(fd);
585
586 unlink(tempfile1);
587 unlink(tempfile2);
588 G_free(in_buf);
589 G_free(out_buf);
590 if(points_row)
591 G_free(points_row);
592 if(points_col)
593 G_free(points_col);
594
595 if (costmode == 1) {
596 close(dir_fd);
597 unlink(tempfile3);
598 G_free(dir_buf);
599 }
600
601 exit(EXIT_SUCCESS);
602}
603
604struct point *drain(int fd, struct point *list, int nrow, int ncol)
605{
606 int go = 1, next_row, next_col;
607 CELL direction;
608 CELL *dir;
609
610 dir = Rast_allocate_c_buf();
611 next_row = list->row;
612 next_col = list->col;
613
614 /* begin loop */
615 while (go) {
616
617 /* find flow direction at this point */
618 lseek(fd, (off_t) list->row * ncol * sizeof(CELL), SEEK_SET);
619 read(fd, dir, ncol * sizeof(CELL));
620 direction = *(dir + list->col);
621 go = 0;
622
623 /* identify next downstream cell */
624 if (direction > 0 && direction < 256) {
625
626 if (direction == 1 || direction == 2 || direction == 4)
627 next_col += 1;
628 else if (direction == 16 || direction == 32 || direction == 64)
629 next_col -= 1;
630
631 if (direction == 64 || direction == 128 || direction == 1)
632 next_row -= 1;
633 else if (direction == 4 || direction == 8 || direction == 16)
634 next_row += 1;
635
636 if (next_col >= 0 && next_col < ncol
637 && next_row >= 0 && next_row < nrow) {
638 /* allocate and fill the next point structure */
639 list->next = (struct point *)G_malloc(sizeof(struct point));
640 list = list->next;
641 list->row = next_row;
642 list->col = next_col;
643 go = 1;
644 }
645 }
646 } /* end while */
647
648 /* allocate and fill the end-of-path flag */
649 list->next = (struct point *)G_malloc(sizeof(struct point));
650 list = list->next;
651 list->row = INT_MAX;
652
653 /* return a pointer to an empty structure */
654 list->next = (struct point *)G_malloc(sizeof(struct point));
655 list = list->next;
656 list->next = NULL;
657
658 G_free(dir);
659
660 return list;
661}
662
663struct point *drain_cost(int dir_fd, struct point *list, int nrow, int ncol)
664{
665 /*
666 * The idea is that each cell of the direction surface has a value representing
667 * the direction towards the next cell in the path. The direction is read from
668 * the input raster, and a simple case/switch is used to determine which cell to
669 * read next. This is repeated via a while loop until a null direction is found.
670 */
671
672 int neighbour, next_row, next_col, go = 1;
673 DCELL direction;
674 DCELL *dir_buf;
675
676 dir_buf = Rast_allocate_d_buf();
677
678 next_row = list->row;
679 next_col = list->col;
680
681 while (go) {
682 go = 0;
683 /* Directional algorithm
684 * 1) read cell direction
685 * 2) shift to cell in that direction
686 */
687 /* find the direction recorded at row,col */
688 lseek(dir_fd, (off_t) list->row * ncol * sizeof(DCELL), SEEK_SET);
689 read(dir_fd, dir_buf, ncol * sizeof(DCELL));
690 direction = *(dir_buf + list->col);
691 neighbour = direction * 10;
692 if (G_verbose() > 2)
693 G_message(_("direction read: %lf, neighbour found: %i"),
694 direction, neighbour);
695 switch (neighbour) {
696 case 225: /* ENE */
697 next_row = list->row - 1;
698 next_col = list->col + 2;
699 break;
700 case 450: /* NE */
701 next_row = list->row - 1;
702 next_col = list->col + 1;
703 break;
704 case 675: /* NNE */
705 next_row = list->row - 2;
706 next_col = list->col + 1;
707 break;
708 case 900: /* N */
709 next_row = list->row - 1;
710 next_col = list->col;
711 break;
712 case 1125: /* NNW */
713 next_row = list->row - 2;
714 next_col = list->col - 1;
715 break;
716 case 1350: /* NW */
717 next_col = list->col - 1;
718 next_row = list->row - 1;
719 break;
720 case 1575: /* WNW */
721 next_col = list->col - 2;
722 next_row = list->row - 1;
723 break;
724 case 1800: /* W*/
725 next_row = list->row;
726 next_col = list->col - 1;
727 break;
728 case 2025: /* WSW */
729 next_row = list->row + 1;
730 next_col = list->col - 2;
731 break;
732 case 2250: /* SW */
733 next_row = list->row + 1;
734 next_col = list->col - 1;
735 break;
736 case 2475: /* SSW */
737 next_row = list->row + 2;
738 next_col = list->col - 1;
739 break;
740 case 2700: /* S */
741 next_row = list->row + 1;
742 next_col = list->col;
743 break;
744 case 2925: /* SSE */
745 next_row = list->row + 2;
746 next_col = list->col + 1;
747 break;
748 case 3150: /* SE */
749 next_row = list->row + 1;
750 next_col = list->col + 1;
751 break;
752 case 3375: /* ESE */
753 next_row = list->row + 1;
754 next_col = list->col + 2;
755 break;
756 case 3600: /* E */
757 next_row = list->row;
758 next_col = list->col + 1;
759 break;
760 /* default:
761 break;
762 Should probably add something here for the possibility of a non-direction map
763 G_fatal_error(_("Invalid direction given (possibly not a direction map)")); */
764 } /* end switch/case */
765
766 if (next_col >= 0 && next_col < ncol && next_row >= 0 &&
767 next_row < nrow) {
768 list->next = (struct point *)G_malloc(sizeof(struct point));
769 list = list->next;
770 list->row = next_row;
771 list->col = next_col;
772 next_row = -1;
773 next_col = -1;
774 go = 1;
775 }
776
777 } /* end while */
778
779 /* allocate and fill the end-of-path flag */
780 list->next = (struct point *)G_malloc(sizeof(struct point));
781 list = list->next;
782 list->row = INT_MAX;
783
784 /* return a pointer to an empty structure */
785 list->next = (struct point *)G_malloc(sizeof(struct point));
786 list = list->next;
787 list->next = NULL;
788
789 G_free(dir_buf);
790
791 return list;
792}
Note: See TracBrowser for help on using the repository browser.