[GRASS-SVN] r71962 - in grass/trunk/raster: . r.path r.path/tests
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Dec 20 14:06:06 PST 2017
Author: mmetz
Date: 2017-12-20 14:06:06 -0800 (Wed, 20 Dec 2017)
New Revision: 71962
Added:
grass/trunk/raster/r.path/
grass/trunk/raster/r.path/pavlrc.c
grass/trunk/raster/r.path/pavlrc.h
grass/trunk/raster/r.path/r.path.html
grass/trunk/raster/r.path/r_path.png
grass/trunk/raster/r.path/r_path_with_r_watershed_direction.png
Removed:
grass/trunk/raster/r.path/README
grass/trunk/raster/r.path/filldir.c
grass/trunk/raster/r.path/r.drain.html
grass/trunk/raster/r.path/r_drain.png
grass/trunk/raster/r.path/r_drain_with_r_watershed_direction.png
grass/trunk/raster/r.path/resolve.c
grass/trunk/raster/r.path/tests/test.r.drain.sh
grass/trunk/raster/r.path/tests/testascii_nc.asc
grass/trunk/raster/r.path/tinf.c
grass/trunk/raster/r.path/tinf.h
Modified:
grass/trunk/raster/Makefile
grass/trunk/raster/r.path/Makefile
grass/trunk/raster/r.path/local.h
grass/trunk/raster/r.path/main.c
Log:
add r.path
Modified: grass/trunk/raster/Makefile
===================================================================
--- grass/trunk/raster/Makefile 2017-12-19 22:55:14 UTC (rev 71961)
+++ grass/trunk/raster/Makefile 2017-12-20 22:06:06 UTC (rev 71962)
@@ -62,6 +62,7 @@
r.out.vtk \
r.param.scale \
r.patch \
+ r.path \
r.profile \
r.proj \
r.quant \
Modified: grass/trunk/raster/r.path/Makefile
===================================================================
--- grass/trunk/raster/r.drain/Makefile 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/Makefile 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,6 +1,6 @@
MODULE_TOPDIR = ../..
-PGM = r.drain
+PGM = r.path
LIBES = $(VECTORLIB) $(RASTERLIB) $(GISLIB)
DEPENDENCIES = $(VECTORDEP) $(RASTERDEP) $(GISDEP)
Deleted: grass/trunk/raster/r.path/README
===================================================================
--- grass/trunk/raster/r.drain/README 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/README 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,42 +0,0 @@
-Date: Sun, 15 Jul 2001 18:26:51 -0600
-From: "Roger S. Miller" <rgrmill at rt66.com>
-X-Mailer: Mozilla 4.76 [en] (X11; U; Linux 2.0.27 i586)
-X-Accept-Language: en
-To: Markus Neteler <neteler at geog.uni-hannover.de>
-Subject: r.drain -- new version
-
-Markus,
-
-I worked on the original r.drain for a while and wasn't really able to
-understand what it was doing, much less where the problem was. I
-figured that I might be able to spend the rest of the weekend figuring
-it out and in the end I would still have a program with unusual behavior
-in flat spots. I decided instead to finish the version that I started
-based on r.fill.dir.
-
-My version of r.drain (which I called r.d2) is in the attached zip file.
-
-r.d2 reproduces most of the behavior of r.drain. It does not use the
-"null_value" option and it won't read input from a sites file. It uses
-all of the flags defined for r.drain. My first draft of the model
-didn't correctly handle the ew,ns and diagonal resolutions (as
-r.fill.dir does not). I believe that I fixed that problem. There are
-also some other minor behaviors that differ a little from r.drain.
-
-r.d2 runs pretty quickly on a small, clean dem, but it can take quite a
-while to run on a large cell map - particularly one with very many
-areas where the flow directions are ambiguous. I didn't write anything
-into the program to warn the user about the execution time.
-
-I tested it on the fcell map that you sent me when this issue first came
-up and on a lat-long cell map. It worked well in both cases. The
-version of r.fill.dir that I started with was the version I had archived
-here, which was not updated for all the changes needed to make
-r.fill.dir run on all of our platforms. I built in the changes that I
-could remember, but there may be a few more that come up when you try to
-compile it on the Cray.
-
-I hope this meets your needs.
-
-
-Roger
Deleted: grass/trunk/raster/r.path/filldir.c
===================================================================
--- grass/trunk/raster/r.drain/filldir.c 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/filldir.c 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,129 +0,0 @@
-#include <grass/config.h>
-#include <sys/types.h>
-#include <stdlib.h>
-#include <math.h>
-#include <limits.h>
-#include <float.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include "tinf.h"
-#include "local.h"
-
-/* get the slope between two cells and return a slope direction */
-void check(CELL newdir, CELL * dir, void *center, void *edge, double cnst,
- double *oldslope)
-{
- double newslope;
-
- /* never discharge to a null boundary */
- if (is_null(edge))
- return;
-
- newslope = slope(center, edge, cnst);
- if (newslope == *oldslope) {
- *dir += newdir;
- }
- else if (newslope > *oldslope) {
- *oldslope = newslope;
- *dir = newdir;
- }
-
- return;
-
-}
-
-/* determine the flow direction at each cell on one row */
-void build_one_row(int i, int nl, int ns, struct band3 *bnd, CELL * dir,
- struct metrics m)
-{
- int j, inc;
- off_t offset;
- CELL sdir;
- double slope;
- char *center;
- char *edge;
-
- inc = bpe();
-
- for (j = 0; j < ns; j += 1) {
- offset = j * bpe();
- center = bnd->b[1] + offset;
- if (is_null(center)) {
- Rast_set_c_null_value(dir + j, 1);
- continue;
- }
-
- sdir = 0;
- /* slope=HUGE; */
- slope = HUGE_VAL;
- if (i == 0) {
- sdir = 128;
- }
- else if (i == nl - 1) {
- sdir = 8;
- }
- else if (j == 0) {
- sdir = 32;
- }
- else if (j == ns - 1) {
- sdir = 2;
- }
- else {
- /* slope=-HUGE; */
- slope = -HUGE_VAL;
-
- /* check one row back */
- edge = bnd->b[0] + offset;
- check(64, &sdir, center, edge - inc, m.diag_res, &slope);
- check(128, &sdir, center, edge, m.ns_res, &slope);
- check(1, &sdir, center, edge + inc, m.diag_res, &slope);
-
- /* check this row */
- check(32, &sdir, center, center - inc, m.ew_res, &slope);
- check(2, &sdir, center, center + inc, m.ew_res, &slope);
-
- /* check one row forward */
- edge = bnd->b[2] + offset;
- check(16, &sdir, center, edge - inc, m.diag_res, &slope);
- check(8, &sdir, center, edge, m.ns_res, &slope);
- check(4, &sdir, center, edge + inc, m.diag_res, &slope);
- }
-
- if (slope == 0.)
- sdir = -sdir;
- else if (slope < 0.)
- sdir = -256;
- dir[j] = sdir;
- }
- return;
-}
-
-void filldir(int fe, int fd, int nl, struct band3 *bnd, struct metrics *m)
-{
- int i, bufsz; /* use off_t bufsz for large files ? MM */
- CELL *dir;
-
- /* determine the flow direction in each cell. On outer rows and columns
- * the flow direction is always directly out of the map */
-
- dir = G_calloc(bnd->ns, sizeof(CELL));
- bufsz = bnd->ns * sizeof(CELL);
-
- lseek(fe, 0, SEEK_SET);
- lseek(fd, 0, SEEK_SET);
- advance_band3(fe, bnd);
- for (i = 0; i < nl; i++) {
- G_percent(i, nl, 5);
- advance_band3(fe, bnd);
- build_one_row(i, nl, bnd->ns, bnd, dir, m[i]);
- write(fd, dir, bufsz);
- }
- G_percent(1, 1, 1);
- advance_band3(fe, bnd);
- build_one_row(nl - 1, nl, bnd->ns, bnd, dir, m[nl - 1]);
- write(fd, dir, bufsz);
-
- G_free(dir);
-
- return;
-}
Modified: grass/trunk/raster/r.path/local.h
===================================================================
--- grass/trunk/raster/r.drain/local.h 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/local.h 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,3 +1,9 @@
+
+#define OUT_PID 1
+#define OUT_CNT 2
+#define OUT_CPY 3
+#define OUT_ACC 4
+
struct metrics
{
double ew_res, ns_res, diag_res;
@@ -2,7 +8 @@
};
-
-void filldir(int, int, int, struct band3 *, struct metrics *);
-void resolve(int, int, struct band3 *);
-int dopolys(int, int, int, int);
-void wtrshed(int, int, int, int, int);
-void ppupdate(int, int, int, int, struct band3 *, struct band3 *);
Modified: grass/trunk/raster/r.path/main.c
===================================================================
--- grass/trunk/raster/r.drain/main.c 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/main.c 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,23 +1,13 @@
/****************************************************************************
*
- * MODULE: r.drain
+ * MODULE: r.path
*
- * AUTHOR(S): Kewan Q. Khawaja <kewan techlogix.com>
- * Update to FP (2000): Pierre de Mouveaux <pmx at audiovu.com> <pmx at free.fr>
- * bugfix in FCELL, DCELL: Markus Neteler 12/2000
- * Rewritten by Roger Miller 7/2001 based on subroutines from
- * r.fill.dir and on the original r.drain.
- * 24 July 2004: WebValley 2004, error checking and vector points added by
- * Matteo Franchi Liceo Leonardo Da Vinci Trento
- * Roberto Flor ITC-irst
- * New code added by Colin Nielsen <colin.nielsen at gmail dot com> *
- * to use movement direction surface from r.walk and r.cost, and to
- * output vector paths 2/2009
+ * AUTHOR(S): based on r.drain
+ * Markus Metz
*
- * PURPOSE: This is the main program for tracing out the path that a
- * drop of water would take if released at a certain location
- * on an input elevation map.
- * COPYRIGHT: (C) 2000,2009 by the GRASS Development Team
+ * PURPOSE: Tracing paths from starting points following
+ * input directions.
+ * COPYRIGHT: (C) 2017 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -42,66 +32,100 @@
#include <grass/gis.h>
#include <grass/raster.h>
+#include <grass/vector.h>
#include <grass/glocale.h>
-#include <grass/vector.h>
-#define DEBUG
-#include "tinf.h"
#include "local.h"
+#include "pavlrc.h"
/* should probably be updated to a pointer array & malloc/realloc as needed */
#define POINTS_INCREMENT 1024
-/* define a data structure to hold the point data */
+/* start points */
struct point
{
int row;
int col;
+ double value;
struct point *next;
+};
+
+/* stack points for bitmask directions */
+struct spoint
+{
+ int row;
+ int col;
+ int dir;
double value;
+ struct spoint *next;
};
+/* output path points */
+struct ppoint
+{
+ int row;
+ int col;
+ double value;
+};
+
+/* managed list of output path points */
+struct point_list
+{
+ struct ppoint *p;
+ int n;
+ int nalloc;
+};
+
+int dir_degree(int dir_fd, int val_fd, struct point *startp, struct Cell_head *window,
+ struct Map_info *Out, struct point_list *pl, int out_mode);
+int dir_bitmask(int dir_fd, int val_fd, struct point *startp, struct Cell_head *window,
+ struct Map_info *Out, struct point_list *pl, int out_mode);
+void pl_add(struct point_list *, struct ppoint *);
+
+/* comparing path points with qsort */
+int cmp_pp(const void *a, const void *b)
+{
+ const struct ppoint *ap = (const struct ppoint *)a;
+ const struct ppoint *bp = (const struct ppoint *)b;
+
+ /* 1. ascending by row */
+ if (ap->row != bp->row)
+ return (ap->row - bp->row);
+ /* 2. ascending by col */
+ if (ap->col != bp->col)
+ return (ap->col - bp->col);
+ /* 3. descending by value */
+ if (ap->value > bp->value)
+ return -1;
+ return (ap->value < bp->value);
+}
+
int main(int argc, char **argv)
{
-
- int fe, fd, dir_fd;
- int i, have_points = 0;
- int new_id;
+ int fd, dir_fd, val_fd;
+ int i, j, have_points = 0;
int nrows, ncols;
- int *points_row = NULL, *points_col = NULL, npoints;
- int increment_count;
- int cell_open(), cell_open_new();
- int map_id, dir_id;
- char map_name[GNAME_MAX], new_map_name[GNAME_MAX], dir_name[GNAME_MAX];
- char *tempfile1, *tempfile2, *tempfile3;
+ int npoints;
+ int out_id, dir_id;
+ char map_name[GNAME_MAX], out_name[GNAME_MAX], dir_name[GNAME_MAX];
+ char *tempfile1, *tempfile2;
struct History history;
struct Cell_head window;
struct Option *opt1, *opt2, *coordopt, *vpointopt, *opt3, *opt4;
struct Flag *flag1, *flag2, *flag3, *flag4;
struct GModule *module;
- int in_type;
- void *in_buf;
void *dir_buf;
- CELL *out_buf;
- struct band3 bnd, bndC;
- struct metrics *m = NULL;
- struct point *list;
- struct point *thispoint;
- int ival, start_row, start_col, mode;
- off_t bsz;
- int costmode = 0;
- double east, north, val;
- struct point *drain(int, struct point *, int, int);
- struct point *drain_cost(int, struct point *, int, int);
- int bsort(int, struct point *);
+ struct point *head_start_pt = NULL;
+ struct point *next_start_pt;
+ struct point_list pl, *ppl;
+ int start_row, start_col, out_mode;
+ double east, north;
struct line_pnts *Points;
struct line_cats *Cats;
- struct Map_info vout;
- int cat;
- double x, y;
+ struct Map_info vout, *pvout;
G_gisinit(argv[0]);
@@ -110,28 +134,28 @@
G_add_keyword(_("hydrology"));
G_add_keyword(_("cost surface"));
module->description =
- _("Traces a flow through an elevation model or cost surface on a raster map.");
+ _("Traces paths from starting points following input directions.");
opt1 = G_define_standard_option(G_OPT_R_INPUT);
- opt1->description = _("Name of input elevation or cost surface raster map");
+ opt1->description = _("Name of input direction");
+ opt1->description =
+ _("Direction in degrees CCW from east, or bitmask encoded (-b flag)");
- opt3 = G_define_standard_option(G_OPT_R_INPUT);
- opt3->key = "direction";
- opt3->label =
- _("Name of input movement direction map associated with the cost surface");
- opt3->description =
- _("Direction in degrees CCW from east");
+ opt2 = G_define_standard_option(G_OPT_R_INPUT);
+ opt2->key = "values";
+ opt2->label =
+ _("Name of input raster values to be used for output");
+ opt2->required = NO;
+
+ opt3 = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt3->key = "raster_path";
opt3->required = NO;
- opt3->guisection = _("Cost surface");
-
- opt2 = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt3->label = _("Name for output raster path map");
opt4 = G_define_standard_option(G_OPT_V_OUTPUT);
- opt4->key = "drain";
+ opt4->key = "vector_path";
opt4->required = NO;
- opt4->label =
- _("Name for output drain vector map");
- opt4->description = _("Recommended for cost surface made using knight's move");
+ opt4->label = _("Name for output vector path map");
coordopt = G_define_standard_option(G_OPT_M_COORDS);
coordopt->key = "start_coordinates";
@@ -161,76 +185,57 @@
flag3->guisection = _("Path settings");
flag4 = G_define_flag();
- flag4->key = 'd';
+ flag4->key = 'b';
flag4->description =
- _("The input raster map is a cost surface (direction surface must also be specified)");
- flag4->guisection = _("Cost surface");
+ _("The input direction map is bitmask encoded");
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
-
- strcpy(map_name, opt1->answer);
- strcpy(new_map_name, opt2->answer);
-
- if (flag4->answer) {
- costmode = 1;
- G_verbose_message(_("Directional drain selected... checking for direction raster map"));
+ /* check args */
+ if (!opt3->answer && !opt4->answer) {
+ G_fatal_error(_("No output requested"));
}
- else {
- G_verbose_message(_("Surface/Hydrology drain selected"));
- }
- if (costmode == 1) {
- if (!opt3->answer) {
- G_fatal_error(_("Direction raster map <%s> not specified, if direction flag is on, "
- "a direction raster must be given"), opt3->key);
- }
- strcpy(dir_name, opt3->answer);
+ strcpy(dir_name, opt1->answer);
+ *map_name = '\0';
+ *out_name = '\0';
+ if (opt3->answer) {
+ strcpy(out_name, opt3->answer);
+ if (opt2->answer)
+ strcpy(map_name, opt2->answer);
}
- if (costmode == 0) {
- if (opt3->answer) {
- G_fatal_error(_("Direction raster map <%s> should not be specified for Surface/Hydrology drains"),
- opt3->answer);
- }
- }
+ pvout = NULL;
if (opt4->answer) {
if (0 > Vect_open_new(&vout, opt4->answer, 0)) {
G_fatal_error(_("Unable to create vector map <%s>"),
opt4->answer);
}
Vect_hist_command(&vout);
+ pvout = &vout;
}
- /* allocate cell buf for the map layer */
- in_type = Rast_map_type(map_name, "");
- /* set the pointers for multi-typed functions */
- set_func_pointers(in_type);
-
if ((flag1->answer + flag2->answer + flag3->answer) > 1)
G_fatal_error(_("Specify just one of the -c, -a and -n flags"));
- mode = 0;
- if (flag1->answer)
- mode = 1;
- if (flag2->answer)
- mode = 2;
- if (flag3->answer)
- mode = 3;
+ out_mode = OUT_PID;
+ if (opt3->answer && opt2->answer) {
+ if (flag1->answer)
+ out_mode = OUT_CPY;
+ if (flag2->answer)
+ out_mode = OUT_ACC;
+ if (flag3->answer)
+ out_mode = OUT_CNT;
+ }
/* get the window information */
G_get_window(&window);
nrows = Rast_window_rows();
ncols = Rast_window_cols();
- /* calculate true cell resolution */
- m = (struct metrics *)G_malloc(nrows * sizeof(struct metrics));
- points_row = (int*)G_calloc(POINTS_INCREMENT, sizeof(int));
- points_col = (int*)G_calloc(POINTS_INCREMENT, sizeof(int));
-
- increment_count = 1;
npoints = 0;
+ /* TODO: use r.cost method to create a list of start points */
if (coordopt->answer) {
for (i = 0; coordopt->answers[i] != NULL; i += 2) {
G_scan_easting(coordopt->answers[i], &east, G_projection());
@@ -244,23 +249,24 @@
i + 1);
continue;
}
- points_row[npoints] = start_row;
- points_col[npoints] = start_col;
npoints++;
- if (npoints == POINTS_INCREMENT * increment_count)
- {
- increment_count++;
- points_row = (int*)G_realloc(points_row, POINTS_INCREMENT * increment_count * sizeof(int));
- points_col = (int*)G_realloc(points_col, POINTS_INCREMENT * increment_count * sizeof(int));
- }
have_points = 1;
+
+ next_start_pt =
+ (struct point *)(G_malloc(sizeof(struct point)));
+
+ next_start_pt->row = start_row;
+ next_start_pt->col = start_col;
+ next_start_pt->value = npoints;
+ next_start_pt->next = head_start_pt;
+ head_start_pt = next_start_pt;
}
}
if (vpointopt->answers) {
for (i = 0; vpointopt->answers[i] != NULL; i++) {
struct Map_info In;
struct bound_box box;
- int type;
+ int cat, type;
Points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
@@ -276,6 +282,7 @@
Vect_rewind(&In);
Vect_region_box(&window, &box);
+ box.T = box.B = 0;
while (1) {
/* register line */
@@ -300,16 +307,18 @@
start_col > ncols)
continue;
- points_row[npoints] = start_row;
- points_col[npoints] = start_col;
npoints++;
- if (npoints == POINTS_INCREMENT * increment_count)
- {
- increment_count++;
- points_row = (int*)G_realloc(points_row, POINTS_INCREMENT * increment_count * sizeof(int));
- points_col = (int*)G_realloc(points_col, POINTS_INCREMENT * increment_count * sizeof(int));
- }
have_points = 1;
+
+ next_start_pt =
+ (struct point *)(G_malloc(sizeof(struct point)));
+
+ next_start_pt->row = start_row;
+ next_start_pt->col = start_col;
+ Vect_cat_get(Cats, 1, &cat);
+ next_start_pt->value = cat;
+ next_start_pt->next = head_start_pt;
+ head_start_pt = next_start_pt;
}
Vect_close(&In);
@@ -323,360 +332,604 @@
}
}
if (have_points == 0)
- G_fatal_error(_("No start/stop point(s) specified"));
+ G_fatal_error(_("No start point(s) specified"));
/* determine the drainage paths */
- /* allocate storage for the first point */
- thispoint = (struct point *)G_malloc(sizeof(struct point));
- list = thispoint;
- thispoint->next = NULL;
+ /* get some temp files */
+ val_fd = -1;
+ tempfile1 = NULL;
+ if (opt2->answer && opt3->answer) {
+ DCELL *map_buf;
- G_begin_distance_calculations();
- {
- double e1, n1, e2, n2;
+ tempfile1 = G_tempfile();
+ val_fd = open(tempfile1, O_RDWR | O_CREAT, 0666);
- e1 = window.east;
- n1 = window.north;
- e2 = e1 + window.ew_res;
- n2 = n1 - window.ns_res;
+ map_buf = Rast_allocate_d_buf();
+
+ fd = Rast_open_old(map_name, "");
+ /* transfer the values map to a temp file */
for (i = 0; i < nrows; i++) {
- m[i].ew_res = G_distance(e1, n1, e2, n1);
- m[i].ns_res = G_distance(e1, n1, e1, n2);
- m[i].diag_res = G_distance(e1, n1, e2, n2);
- e2 = e1 + window.ew_res;
- n2 = n1 - window.ns_res;
+ Rast_get_d_row(fd, map_buf, i);
+ if (write(val_fd, map_buf, ncols * sizeof(DCELL)) !=
+ ncols * sizeof(DCELL)) {
+ G_fatal_error(_("Unable to write to tempfile"));
+ }
}
+ Rast_close(fd);
+ G_free(map_buf);
}
- /* buffers for internal use */
- bndC.ns = ncols;
- bndC.sz = sizeof(CELL) * ncols;
- bndC.b[0] = G_calloc(ncols, sizeof(CELL));
- bndC.b[1] = G_calloc(ncols, sizeof(CELL));
- bndC.b[2] = G_calloc(ncols, sizeof(CELL));
-
- /* buffers for external use */
- bnd.ns = ncols;
- bnd.sz = ncols * bpe();
- bnd.b[0] = G_calloc(ncols, bpe());
- bnd.b[1] = G_calloc(ncols, bpe());
- bnd.b[2] = G_calloc(ncols, bpe());
-
- /* an input buffer */
- in_buf = get_buf();
-
- /* open the original map and get its file id */
- map_id = Rast_open_old(map_name, "");
-
- /* get some temp files */
- tempfile1 = G_tempfile();
+ dir_id = Rast_open_old(dir_name, "");
tempfile2 = G_tempfile();
+ dir_fd = open(tempfile2, O_RDWR | O_CREAT, 0666);
- fe = open(tempfile1, O_RDWR | O_CREAT, 0666);
- fd = open(tempfile2, O_RDWR | O_CREAT, 0666);
-
- /* transfer the input map to a temp file */
- for (i = 0; i < nrows; i++) {
- get_row(map_id, in_buf, i);
- write(fe, in_buf, bnd.sz);
+ if (flag4->answer) {
+ dir_buf = Rast_allocate_c_buf();
+ for (i = 0; i < nrows; i++) {
+ Rast_get_c_row(dir_id, dir_buf, i);
+ if (write(dir_fd, dir_buf, ncols * sizeof(CELL)) !=
+ ncols * sizeof(CELL)) {
+ G_fatal_error(_("Unable to write to tempfile"));
+ }
+ }
}
- Rast_close(map_id);
-
- if (costmode == 1) {
+ else {
dir_buf = Rast_allocate_d_buf();
- dir_id = Rast_open_old(dir_name, "");
- tempfile3 = G_tempfile();
- dir_fd = open(tempfile3, O_RDWR | O_CREAT, 0666);
-
for (i = 0; i < nrows; i++) {
Rast_get_d_row(dir_id, dir_buf, i);
- write(dir_fd, dir_buf, ncols * sizeof(DCELL));
+ if (write(dir_fd, dir_buf, ncols * sizeof(DCELL)) !=
+ ncols * sizeof(DCELL)) {
+ G_fatal_error(_("Unable to write to tempfile"));
+ }
}
- Rast_close(dir_id);
}
+ Rast_close(dir_id);
+ G_free(dir_buf);
- /* only necessary for non-dir drain */
- if (costmode == 0) {
- G_message(_("Calculating flow directions..."));
+ ppl = NULL;
+ if (opt3->answer) {
+ /* raster output */
+ pl.p = NULL;
+ pl.n = 0;
+ pl.nalloc = 0;
+ ppl = &pl;
+ }
- /* fill one-cell pits and take a first stab at flow directions */
- filldir(fe, fd, nrows, &bnd, m);
+ /* determine the drainage paths */
- /* determine flow directions for more ambiguous cases */
- resolve(fd, nrows, &bndC);
- }
+ /* repeat for each starting point */
- /* free the buffers already used */
- G_free(bndC.b[0]);
- G_free(bndC.b[1]);
- G_free(bndC.b[2]);
+ next_start_pt = head_start_pt;
+ while (next_start_pt) {
+ /* follow directions from start points to determine paths */
+ /* path tracing algorithm selection */
+ if (flag4->answer) {
+ struct Map_info Tmp;
- G_free(bnd.b[0]);
- G_free(bnd.b[1]);
- G_free(bnd.b[2]);
+ if (pvout) {
+ if (Vect_open_tmp_new(&Tmp, NULL, 0) < 0)
+ G_fatal_error(_("Unable to create temporary vector map"));
+ pvout = &Tmp;
+ }
- /* determine the drainage paths */
-
- /* repeat for each starting point */
- for (i = 0; i < npoints; i++) {
- /* use the flow directions to determine the drainage path
- * results are compiled as a linked list of points in downstream order */
- thispoint->row = points_row[i];
- thispoint->col = points_col[i];
- thispoint->next = NULL;
- /* drain algorithm selection (dir or non-dir) */
- if (costmode == 0) {
- thispoint = drain(fd, thispoint, nrows, ncols);
+ if (!dir_bitmask(dir_fd, val_fd, next_start_pt, &window,
+ pvout, ppl, out_mode)) {
+ G_warning(_("No path at row %d, col %d"),
+ next_start_pt->row, next_start_pt->col);
+ }
+ if (pvout) {
+ Vect_build_partial(&Tmp, GV_BUILD_BASE);
+ G_message(_("Breaking lines..."));
+ Vect_break_lines(&Tmp, GV_LINE, NULL);
+ Vect_copy_map_lines(&Tmp, &vout);
+ Vect_set_release_support(&Tmp);
+ Vect_close(&Tmp); /* temporary map is deleted automatically */
+ pvout = &vout;
+ }
}
- if (costmode == 1) {
- thispoint = drain_cost(dir_fd, thispoint, nrows, ncols);
+ else {
+ if (!dir_degree(dir_fd, val_fd, next_start_pt, &window,
+ pvout, ppl, out_mode)) {
+ G_warning(_("No path at row %d, col %d"),
+ next_start_pt->row, next_start_pt->col);
+ }
}
+ next_start_pt = next_start_pt->next;
}
- /* do the output */
+ /* raster output */
+ if (opt3->answer) {
+ int row;
- if (mode == 0 || mode == 3) {
+ if (pl.n > 1) {
+ /* sort points */
+ qsort(pl.p, pl.n, sizeof(struct ppoint), cmp_pp);
- /* Output will be a cell map */
- /* open a new file and allocate an output buffer */
- new_id = Rast_open_c_new(new_map_name);
- out_buf = Rast_allocate_c_buf();
- /* mark each cell */
- thispoint = list;
- while (thispoint->next != NULL) {
- thispoint->value = 1;
- thispoint = thispoint->next;
- }
-
- if (mode == 3) {
- /* number each cell downstream */
- thispoint = list;
- ival = 0;
- while (thispoint->next != NULL) {
- if (thispoint->row == INT_MAX) {
- ival = 0;
- thispoint = thispoint->next;
- continue;
+ /* remove duplicates */
+ j = 1;
+ for (i = 1; i < pl.n; i++) {
+ if (pl.p[j - 1].row != pl.p[i].row ||
+ pl.p[j - 1].col != pl.p[i].col) {
+
+ pl.p[j].row = pl.p[i].row;
+ pl.p[j].col = pl.p[i].col;
+ pl.p[j].value = pl.p[i].value;
+ j++;
}
- thispoint->value += ival;
- ival = thispoint->value;
- thispoint = thispoint->next;
}
+ pl.n = j;
}
- /* build the output map */
- G_message(_("Writing output raster map..."));
- for (i = 0; i < nrows; i++) {
- G_percent(i, nrows, 2);
- Rast_set_c_null_value(out_buf, ncols);
- thispoint = list;
- while (thispoint->next != NULL) {
- if (thispoint->row == i)
- out_buf[thispoint->col] = (int)thispoint->value;
- thispoint = thispoint->next;
- }
- Rast_put_c_row(new_id, out_buf);
- }
- G_percent(1, 1, 1);
- }
- else { /* mode = 1 or 2 */
- /* Output will be of the same type as input */
- /* open a new file and allocate an output buffer */
- new_id = Rast_open_new(new_map_name, in_type);
- out_buf = get_buf();
- bsz = ncols * bpe();
+ if (out_mode == OUT_PID || out_mode == OUT_CNT) {
+ CELL *out_buf;
- /* loop through each point in the list and store the map values */
- thispoint = list;
- while (thispoint->next != NULL) {
- if (thispoint->row == INT_MAX) {
- thispoint = thispoint->next;
- continue;
+ /* Output will be a cell map */
+ /* open a new file and allocate an output buffer */
+ out_id = Rast_open_c_new(out_name);
+ out_buf = Rast_allocate_c_buf();
+ Rast_set_c_null_value(out_buf, window.cols);
+ row = 0;
+
+ /* build the output map */
+ G_message(_("Writing output raster map..."));
+ for (i = 0; i < pl.n; i++) {
+ while (row < pl.p[i].row) {
+ G_percent(row, nrows, 2);
+ Rast_put_c_row(out_id, out_buf);
+ Rast_set_c_null_value(out_buf, window.cols);
+ row++;
+ }
+ out_buf[pl.p[i].col] = pl.p[i].value;
}
- lseek(fe, (off_t) thispoint->row * bsz, SEEK_SET);
- read(fe, in_buf, bsz);
- memcpy(&thispoint->value, (char *)in_buf + bpe() * thispoint->col,
- bpe());
- thispoint = thispoint->next;
+ while (row < window.rows) {
+ G_percent(row, nrows, 2);
+ Rast_put_c_row(out_id, out_buf);
+ Rast_set_c_null_value(out_buf, window.cols);
+ row++;
+ }
+ G_percent(1, 1, 1);
+ G_free(out_buf);
}
+ else { /* mode = OUT_CPY or OUT_ACC */
+ DCELL *out_buf;
- if (mode == 2) {
- /* accumulate the input map values downstream */
- thispoint = list;
- val = 0.;
- while (thispoint->next != NULL) {
- if (thispoint->row == INT_MAX) {
- val = 0.;
- thispoint = thispoint->next;
- continue;
+ /* Output could be of the same type as input */
+ /* open a new file and allocate an output buffer */
+ out_id = Rast_open_new(out_name, DCELL_TYPE);
+ out_buf = Rast_allocate_d_buf();
+ Rast_set_d_null_value(out_buf, window.cols);
+ row = 0;
+
+ /* build the output map */
+ G_message(_("Writing output raster map..."));
+ for (i = 0; i < pl.n; i++) {
+ while (row < pl.p[i].row) {
+ G_percent(row, nrows, 2);
+ Rast_put_d_row(out_id, out_buf);
+ Rast_set_d_null_value(out_buf, window.cols);
+ row++;
}
- sum(&thispoint->value, &val);
- memcpy(&val, &thispoint->value, bpe());
- thispoint = thispoint->next;
+ out_buf[pl.p[i].col] = pl.p[i].value;
}
+ while (row < window.rows) {
+ G_percent(row, nrows, 2);
+ Rast_put_d_row(out_id, out_buf);
+ Rast_set_d_null_value(out_buf, window.cols);
+ row++;
+ }
+ G_percent(1, 1, 1);
+ G_free(out_buf);
}
- /* build the output map */
- G_message(_("Writing raster map <%s>..."),
- new_map_name);
- for (i = 0; i < nrows; i++) {
- G_percent(i, nrows, 2);
- set_null_value(out_buf, ncols);
- thispoint = list;
- while (thispoint->next != NULL) {
- if (thispoint->row == i)
- memcpy((char *)out_buf + bpe() * thispoint->col,
- &(thispoint->value), bpe());
- thispoint = thispoint->next;
- }
- put_row(new_id, out_buf);
- }
- G_percent(1, 1, 1);
+ Rast_close(out_id);
+
+ Rast_put_cell_title(out_name, "Path trace");
+
+ Rast_short_history(out_name, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(out_name, &history);
}
- /* Output a vector path */
+ /* vector output */
if (opt4->answer) {
- Points = Vect_new_line_struct();
- Cats = Vect_new_cats_struct();
- /* Need to modify for multiple paths */
- thispoint = list;
- i = 1;
- while (thispoint->next != NULL) {
- if (thispoint->row == INT_MAX) {
- thispoint = thispoint->next;
- Vect_cat_set(Cats, 1, i);
- Vect_write_line(&vout, GV_LINE, Points, Cats);
- Vect_reset_line(Points);
- Vect_reset_cats(Cats);
- i++;
- continue;
- }
- if (Vect_cat_get(Cats, 1, &cat) == 0) {
- Vect_cat_set(Cats, 1, i);
- }
- /* Need to convert row and col to coordinates
- * y = cell_head.north - ((double) p->row + 0.5) * cell_head.ns_res;
- * x = cell_head.west + ((double) p->col + 0.5) * cell_head.ew_res;
- */
-
- x = window.west + ((double)thispoint->col + 0.5) * window.ew_res;
- y = window.north - ((double)thispoint->row + 0.5) * window.ns_res;
- Vect_append_point(Points, x, y, 0.0);
- thispoint = thispoint->next;
- } /* End while */
Vect_build(&vout);
Vect_close(&vout);
}
/* close files and free buffers */
- Rast_close(new_id);
-
- Rast_put_cell_title(new_map_name, "Surface flow trace");
-
- Rast_short_history(new_map_name, "raster", &history);
- Rast_command_history(&history);
- Rast_write_history(new_map_name, &history);
-
- close(fe);
- close(fd);
-
- unlink(tempfile1);
+ close(dir_fd);
unlink(tempfile2);
- G_free(in_buf);
- G_free(out_buf);
- if(points_row)
- G_free(points_row);
- if(points_col)
- G_free(points_col);
- if (costmode == 1) {
- close(dir_fd);
- unlink(tempfile3);
- G_free(dir_buf);
+ if (opt2->answer && opt3->answer) {
+ close(val_fd);
+ unlink(tempfile1);
}
-
+
exit(EXIT_SUCCESS);
}
-struct point *drain(int fd, struct point *list, int nrow, int ncol)
+void pl_add(struct point_list *pl, struct ppoint *p)
{
- int go = 1, next_row, next_col;
+ if (pl->n == pl->nalloc) {
+ pl->nalloc += POINTS_INCREMENT;
+ pl->p = (struct ppoint *)G_realloc(pl->p, pl->nalloc * sizeof(struct ppoint));
+ if (!pl->p)
+ G_fatal_error(_("Unable to increase point list"));
+ }
+ pl->p[pl->n] = *p;
+ pl->n++;
+}
+
+
+/* Jenson Domingue 1988, r.fill.dir:
+ * bitmask encoded directions CW from North
+ * value log2(value) direction
+ * 1 0 NE
+ * 2 1 E
+ * 4 2 SE
+ * 8 3 S
+ * 16 4 SW
+ * 32 5 W
+ * 64 6 NW
+ * 128 7 N
+ * 256 8 NEN
+ * 512 9 ENE
+ * 1024 10 ESE
+ * 2048 11 SES
+ * 4096 12 SWS
+ * 8192 13 WSW
+ * 16384 14 WNW
+ * 32768 15 NWN
+ *
+ *
+ * 8 neighbours
+ * 64 128 1
+ * 32 x 2
+ * 16 8 4
+ *
+ * log2(value) CW from North starting at NE
+ * expanded for Knight's move
+ *
+ * 15 8
+ * 14 6 7 0 9
+ * 5 X 1
+ * 13 4 3 2 10
+ * 12 11
+ */
+int dir_bitmask(int dir_fd, int val_fd, struct point *startp, struct Cell_head *window,
+ struct Map_info *Out, struct point_list *pl, int out_mode)
+{
+ int go = 1, next_row, next_col, next_dir;
+ int npoints;
+ int dir_row = -1, val_row = -1;
CELL direction;
- CELL *dir;
+ CELL *dir_buf;
+ DCELL *val_buf;
+ struct spoint *stackp, *newp;
+ struct pavlrc_table *visited;
+ struct pavlrc ngbr_rc;
+ struct ppoint pp;
+ int is_stack;
+ int cur_dir, i, npaths;
+ struct line_pnts *Points;
+ struct line_cats *Cats;
+ double x, y;
+ double value;
+ int col_offset[16] = { 1, 1, 1, 0, -1, -1, -1, 0,
+ 1, 2, 2, 1, -1, -2, -2, -1 };
+ int row_offset[16] = { -1, 0, 1, 1, 1, 0, -1, -1,
+ -2, -1, 1, 2, 2, 1, -1, -2 };
- dir = Rast_allocate_c_buf();
- next_row = list->row;
- next_col = list->col;
+ dir_buf = Rast_allocate_c_buf();
+ val_buf = NULL;
- /* begin loop */
- while (go) {
+ stackp = (struct spoint *)G_malloc(sizeof(struct spoint));
- /* find flow direction at this point */
- lseek(fd, (off_t) list->row * ncol * sizeof(CELL), SEEK_SET);
- read(fd, dir, ncol * sizeof(CELL));
- direction = *(dir + list->col);
- go = 0;
+ stackp->row = startp->row;
+ stackp->col = startp->col;
+ stackp->value = startp->value;
+ value = startp->value;
+ stackp->dir = -1;
+ stackp->next = NULL;
- /* identify next downstream cell */
- if (direction > 0 && direction < 256) {
+ visited = pavlrc_create(NULL);
+ ngbr_rc.row = stackp->row;
+ ngbr_rc.col = stackp->col;
+ pavlrc_insert(visited, &ngbr_rc);
- if (direction == 1 || direction == 2 || direction == 4)
- next_col += 1;
- else if (direction == 16 || direction == 32 || direction == 64)
- next_col -= 1;
+ if (Out) {
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+ Vect_cat_set(Cats, 1, (int)stackp->value);
+ }
- if (direction == 64 || direction == 128 || direction == 1)
- next_row -= 1;
- else if (direction == 4 || direction == 8 || direction == 16)
- next_row += 1;
-
- if (next_col >= 0 && next_col < ncol
- && next_row >= 0 && next_row < nrow) {
- /* allocate and fill the next point structure */
- list->next = (struct point *)G_malloc(sizeof(struct point));
- list = list->next;
- list->row = next_row;
- list->col = next_col;
- go = 1;
+ if (pl) {
+ if (out_mode == OUT_CNT)
+ value = 1;
+ else if (out_mode == OUT_CPY || out_mode == OUT_ACC) {
+ /* read input raster */
+ val_buf = Rast_allocate_d_buf();
+ if (val_row != stackp->row) {
+ lseek(val_fd, (off_t) stackp->row * window->cols * sizeof(DCELL),
+ SEEK_SET);
+ if (read(val_fd, val_buf, window->cols * sizeof(DCELL)) !=
+ window->cols * sizeof(DCELL)) {
+ G_fatal_error(_("Unable to read from temp file"));
+ }
+ val_row = stackp->row;
}
+ value = val_buf[stackp->col];
}
- } /* end while */
+ stackp->value = value;
+ }
- /* allocate and fill the end-of-path flag */
- list->next = (struct point *)G_malloc(sizeof(struct point));
- list = list->next;
- list->row = INT_MAX;
+ /* multi-path tracing:
+ * multiple paths from the same starting point */
- /* return a pointer to an empty structure */
- list->next = (struct point *)G_malloc(sizeof(struct point));
- list = list->next;
- list->next = NULL;
+ npoints = 0;
+ while (stackp) {
+ is_stack = 1;
+ stackp->dir++;
+ next_row = stackp->row;
+ next_col = stackp->col;
+ value = stackp->value;
- G_free(dir);
+ /* begin loop */
+ go = 1;
+ while (go) {
+ go = 0;
- return list;
+ /* find direction from this point */
+ if (dir_row != next_row) {
+ lseek(dir_fd, (off_t) next_row * window->cols * sizeof(CELL),
+ SEEK_SET);
+ if (read(dir_fd, dir_buf, window->cols * sizeof(CELL)) !=
+ window->cols * sizeof(CELL)) {
+ G_fatal_error(_("Unable to read from temp file"));
+ }
+ dir_row = next_row;
+ }
+ direction = *(dir_buf + next_col);
+
+ if (direction <= 0 || Rast_is_c_null_value(&direction)) {
+ /* end of current path: write out the result */
+ if (Out && Points->n_points > 1) {
+ Vect_write_line(Out, GV_LINE, Points, Cats);
+ }
+
+ /* leave this loop */
+ break;
+ }
+
+ /* start direction */
+ cur_dir = 0;
+ if (is_stack)
+ cur_dir = stackp->dir;
+
+ /* count paths going from current point and
+ * get next direction as log2(direction) */
+ next_dir = -1;
+ npaths = 0;
+ for (i = cur_dir; i < 16; i++) {
+ if ((direction & (1 << i)) != 0) {
+ npaths++;
+ if (next_dir < 0)
+ next_dir = i;
+ }
+ }
+
+ if (is_stack) {
+ if (npaths > 0) {
+ stackp->dir = next_dir;
+ /* start a new path */
+ if (Out) {
+ Vect_reset_line(Points);
+ x = window->west + (next_col + 0.5) * window->ew_res;
+ y = window->north - (next_row + 0.5) * window->ns_res;
+ Vect_append_point(Points, x, y, 0.0);
+ }
+ if (pl) {
+ value = stackp->value;
+ pp.row = next_row;
+ pp.col = next_col;
+ pp.value = value;
+ pl_add(pl, &pp);
+ }
+ npoints++;
+ }
+ else {
+ /* stack point without path ? */
+ if (stackp->dir == 0)
+ G_warning(_("No path from row %d, col %d"),
+ stackp->row, stackp->col);
+
+ /* drop this point from the stack */
+ G_debug(1, "drop point from stack");
+ newp = stackp->next;
+ G_free(stackp);
+ stackp = newp;
+ break;
+ }
+ }
+ else { /* not a stack point */
+ if (npaths == 0) {
+ /* should not happen */
+ G_fatal_error(_("Invalid direction %d"), direction);
+ }
+ if (npaths > 1) {
+ /* write out the result */
+ if (Out && Points->n_points > 1) {
+ Vect_write_line(Out, GV_LINE, Points, Cats);
+ }
+
+ ngbr_rc.row = next_row;
+ ngbr_rc.col = next_col;
+
+ /* add it to the stack and leave this loop */
+ G_debug(1, "add point to stack: row %d, col %d, dir %d",
+ next_row, next_col, next_dir);
+ newp = (struct spoint *)G_malloc(sizeof(struct spoint));
+ newp->row = next_row;
+ newp->col = next_col;
+ newp->dir = next_dir - 1;
+ newp->value = value;
+ newp->next = stackp;
+ stackp = newp;
+ break;
+ }
+ }
+
+ is_stack = 0;
+
+ /* identify next downstream cell */
+ next_row += row_offset[next_dir];
+ next_col += col_offset[next_dir];
+
+ G_debug(1, "next cell at row %d, col %d", next_row, next_col);
+
+ if (next_col >= 0 && next_col < window->cols
+ && next_row >= 0 && next_row < window->rows) {
+
+ /* add point */
+ if (Out) {
+ x = window->west + (next_col + 0.5) * window->ew_res;
+ y = window->north - (next_row + 0.5) * window->ns_res;
+ Vect_append_point(Points, x, y, 0.0);
+ }
+ if (pl) {
+ if (out_mode == OUT_CNT)
+ value += 1;
+ else if (out_mode == OUT_CPY || out_mode == OUT_ACC) {
+ /* read input raster */
+ if (val_row != next_row) {
+ lseek(val_fd, (off_t) next_row * window->cols * sizeof(DCELL),
+ SEEK_SET);
+ if (read(val_fd, val_buf, window->cols * sizeof(DCELL)) !=
+ window->cols * sizeof(DCELL)) {
+ G_fatal_error(_("Unable to read from temp file"));
+ }
+ val_row = next_row;
+ }
+
+ if (out_mode == OUT_CPY)
+ value = val_buf[next_col];
+ else
+ value += val_buf[next_col];
+ }
+ pp.row = next_row;
+ pp.col = next_col;
+ pp.value = value;
+ pl_add(pl, &pp);
+ }
+
+ /* avoid circular paths
+ * avoid tracing the same path segment several times:
+ * a path can split and later on merge again,
+ * then split again, merge again
+ * path segments are identical from every merge point
+ * to the next split point */
+ ngbr_rc.row = next_row;
+ ngbr_rc.col = next_col;
+ if (!pavlrc_insert(visited, &ngbr_rc)) {
+ go = 1;
+ npoints++;
+ }
+ else {
+ /* reached a merge point */
+ /* write out the result */
+ if (Out && Points->n_points > 1) {
+ Vect_write_line(Out, GV_LINE, Points, Cats);
+ }
+ }
+ }
+ else {
+ G_warning(_("Path is leaving the current region"));
+ }
+ } /* end while */
+ }
+
+ pavlrc_destroy(visited);
+
+ G_free(dir_buf);
+ if (val_buf)
+ G_free(val_buf);
+
+ if (Out) {
+ Vect_destroy_line_struct(Points);
+ Vect_destroy_cats_struct(Cats);
+ }
+
+ return (npoints > 1);
}
-struct point *drain_cost(int dir_fd, struct point *list, int nrow, int ncol)
+int dir_degree(int dir_fd, int val_fd, struct point *startp, struct Cell_head *window,
+ struct Map_info *Out, struct point_list *pl, int out_mode)
{
/*
* The idea is that each cell of the direction surface has a value representing
* the direction towards the next cell in the path. The direction is read from
* the input raster, and a simple case/switch is used to determine which cell to
* read next. This is repeated via a while loop until a null direction is found.
+ *
+ * directions are degrees CCW from East with East = 360
*/
int neighbour, next_row, next_col, go = 1;
+ int npoints;
+ int dir_row = -1, val_row = -1;
DCELL direction;
DCELL *dir_buf;
+ DCELL *val_buf;
+ double x, y;
+ double value;
+ struct ppoint pp;
+ struct line_pnts *Points;
+ struct line_cats *Cats;
dir_buf = Rast_allocate_d_buf();
+ val_buf = NULL;
- next_row = list->row;
- next_col = list->col;
+ next_row = startp->row;
+ next_col = startp->col;
+ value = startp->value;
+ if (Out) {
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+ Vect_cat_set(Cats, 1, (int)value);
+
+ x = window->west + (next_col + 0.5) * window->ew_res;
+ y = window->north - (next_row + 0.5) * window->ns_res;
+ Vect_append_point(Points, x, y, 0.0);
+ }
+
+ if (pl) {
+ if (out_mode == OUT_CNT)
+ value = 1;
+ else if (out_mode == OUT_CPY || out_mode == OUT_ACC) {
+ /* read input raster */
+ val_buf = Rast_allocate_d_buf();
+ if (val_row != next_row) {
+ lseek(val_fd, (off_t) next_row * window->cols * sizeof(DCELL),
+ SEEK_SET);
+ if (read(val_fd, val_buf, window->cols * sizeof(DCELL)) !=
+ window->cols * sizeof(DCELL)) {
+ G_fatal_error(_("Unable to read from temp file"));
+ }
+ val_row = next_row;
+ }
+
+ value = val_buf[next_col];
+ }
+ pp.row = next_row;
+ pp.col = next_col;
+ pp.value = value;
+ pl_add(pl, &pp);
+ }
+
+ npoints = 1;
while (go) {
go = 0;
/* Directional algorithm
@@ -684,108 +937,142 @@
* 2) shift to cell in that direction
*/
/* find the direction recorded at row,col */
- lseek(dir_fd, (off_t) list->row * ncol * sizeof(DCELL), SEEK_SET);
- read(dir_fd, dir_buf, ncol * sizeof(DCELL));
- direction = *(dir_buf + list->col);
- neighbour = direction * 10;
- if (G_verbose() > 2)
- G_message(_("direction read: %lf, neighbour found: %i"),
- direction, neighbour);
+ if (dir_row != next_row) {
+ lseek(dir_fd, (off_t) next_row * window->cols * sizeof(DCELL),
+ SEEK_SET);
+ if (read(dir_fd, dir_buf, window->cols * sizeof(DCELL)) !=
+ window->cols * sizeof(DCELL)) {
+ G_fatal_error(_("Unable to read from temp file"));
+ }
+ dir_row = next_row;
+ }
+ direction = *(dir_buf + next_col);
+ neighbour = 0;
+ if (!Rast_is_d_null_value(&direction)) {
+ neighbour = direction * 10;
+ G_debug(2, "direction read: %lf, neighbour found: %i",
+ direction, neighbour);
+ }
switch (neighbour) {
case 225: /* ENE */
- next_row = list->row - 1;
- next_col = list->col + 2;
+ next_row -= 1;
+ next_col += 2;
break;
case 450: /* NE */
- next_row = list->row - 1;
- next_col = list->col + 1;
+ next_row -= 1;
+ next_col += 1;
break;
case 675: /* NNE */
- next_row = list->row - 2;
- next_col = list->col + 1;
+ next_row -= 2;
+ next_col += 1;
break;
case 900: /* N */
- next_row = list->row - 1;
- next_col = list->col;
+ next_row -= 1;
break;
case 1125: /* NNW */
- next_row = list->row - 2;
- next_col = list->col - 1;
+ next_row -= 2;
+ next_col -= 1;
break;
case 1350: /* NW */
- next_col = list->col - 1;
- next_row = list->row - 1;
+ next_col -= 1;
+ next_row -= 1;
break;
case 1575: /* WNW */
- next_col = list->col - 2;
- next_row = list->row - 1;
+ next_col -= 2;
+ next_row -= 1;
break;
case 1800: /* W*/
- next_row = list->row;
- next_col = list->col - 1;
+ next_col -= 1;
break;
case 2025: /* WSW */
- next_row = list->row + 1;
- next_col = list->col - 2;
+ next_row += 1;
+ next_col -= 2;
break;
case 2250: /* SW */
- next_row = list->row + 1;
- next_col = list->col - 1;
+ next_row += 1;
+ next_col -= 1;
break;
case 2475: /* SSW */
- next_row = list->row + 2;
- next_col = list->col - 1;
+ next_row += 2;
+ next_col -= 1;
break;
case 2700: /* S */
- next_row = list->row + 1;
- next_col = list->col;
+ next_row += 1;
break;
case 2925: /* SSE */
- next_row = list->row + 2;
- next_col = list->col + 1;
+ next_row += 2;
+ next_col += 1;
break;
case 3150: /* SE */
- next_row = list->row + 1;
- next_col = list->col + 1;
+ next_row += 1;
+ next_col += 1;
break;
case 3375: /* ESE */
- next_row = list->row + 1;
- next_col = list->col + 2;
+ next_row += 1;
+ next_col += 2;
break;
case 3600: /* E */
- next_row = list->row;
- next_col = list->col + 1;
+ next_col += 1;
break;
- /* default:
- break;
- Should probably add something here for the possibility of a non-direction map
- G_fatal_error(_("Invalid direction given (possibly not a direction map)")); */
+ default:
+ /* end of path */
+ next_row = -1;
+ next_col = -1;
+ break;
} /* end switch/case */
- if (next_col >= 0 && next_col < ncol && next_row >= 0 &&
- next_row < nrow) {
- list->next = (struct point *)G_malloc(sizeof(struct point));
- list = list->next;
- list->row = next_row;
- list->col = next_col;
- next_row = -1;
- next_col = -1;
+ if (next_col >= 0 && next_col < window->cols && next_row >= 0 &&
+ next_row < window->rows) {
+
+ if (Out) {
+ x = window->west + (next_col + 0.5) * window->ew_res;
+ y = window->north - (next_row + 0.5) * window->ns_res;
+ Vect_append_point(Points, x, y, 0.0);
+ }
+ if (pl) {
+ if (out_mode == OUT_CNT)
+ value += 1;
+ else if (out_mode == OUT_CPY || out_mode == OUT_ACC) {
+ /* read input raster */
+ if (val_row != next_row) {
+ lseek(val_fd, (off_t) next_row * window->cols * sizeof(DCELL),
+ SEEK_SET);
+ if (read(val_fd, val_buf, window->cols * sizeof(DCELL)) !=
+ window->cols * sizeof(DCELL)) {
+ G_fatal_error(_("Unable to read from temp file"));
+ }
+ val_row = next_row;
+ }
+
+ if (out_mode == OUT_CPY)
+ value = val_buf[next_col];
+ else
+ value += val_buf[next_col];
+ }
+ pp.row = next_row;
+ pp.col = next_col;
+ pp.value = value;
+ pl_add(pl, &pp);
+ }
+
go = 1;
+ npoints++;
}
} /* end while */
- /* allocate and fill the end-of-path flag */
- list->next = (struct point *)G_malloc(sizeof(struct point));
- list = list->next;
- list->row = INT_MAX;
+ if (Out && Points->n_points > 1) {
+ Vect_write_line(Out, GV_LINE, Points, Cats);
+ }
- /* return a pointer to an empty structure */
- list->next = (struct point *)G_malloc(sizeof(struct point));
- list = list->next;
- list->next = NULL;
-
G_free(dir_buf);
+ if (val_buf)
+ G_free(val_buf);
- return list;
+ if (Out) {
+ Vect_destroy_line_struct(Points);
+ Vect_destroy_cats_struct(Cats);
+ }
+
+ return (npoints > 1);
}
Added: grass/trunk/raster/r.path/pavlrc.c
===================================================================
--- grass/trunk/raster/r.path/pavlrc.c (rev 0)
+++ grass/trunk/raster/r.path/pavlrc.c 2017-12-20 22:06:06 UTC (rev 71962)
@@ -0,0 +1,843 @@
+/* Produced by texiweb from libavl.w. */
+
+/* libavl - library for manipulation of binary trees.
+ Copyright (C) 1998, 1999, 2000, 2001, 2002, 2004 Free Software
+ Foundation, Inc.
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 3 of the License, or (at your option) any later version.
+
+ This library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with this library; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ 02110-1301 USA.
+ */
+
+/* Nov 2016, Markus Metz
+ * from libavl-2.0.3
+ * added safety checks and speed optimizations
+ */
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "pavlrc.h"
+
+#define CMP_RC(a, b) ((a)->row == (b)->row ? (a)->col - (b)->col : (a)->row - (b)->row)
+
+/* Creates and returns a new table
+ with comparison function |compare| using parameter |param|
+ and memory allocator |allocator|.
+ Returns |NULL| if memory allocation failed. */
+struct pavlrc_table *pavlrc_create(struct libavl_allocator *allocator)
+{
+ struct pavlrc_table *tree;
+
+ if (allocator == NULL)
+ allocator = &pavlrc_allocator_default;
+
+ tree = allocator->libavl_malloc(allocator, sizeof *tree);
+ if (tree == NULL)
+ return NULL;
+
+ tree->pavl_root = NULL;
+ tree->pavl_alloc = allocator;
+ tree->pavl_count = 0;
+
+ return tree;
+}
+
+/* Search |tree| for an item matching |item|, and return it if found.
+ Otherwise return |NULL|. */
+struct pavlrc *pavlrc_find(const struct pavlrc_table *tree, const struct pavlrc *item)
+{
+ struct pavlrc_node *p;
+
+ assert(tree != NULL && item != NULL);
+
+ p = tree->pavl_root;
+ while (p != NULL) {
+ int cmp = CMP_RC(item, &p->pavl_data);
+
+ if (cmp == 0)
+ return &p->pavl_data;
+
+ p = p->pavl_link[cmp > 0];
+ }
+
+ return NULL;
+}
+
+/* Inserts |item| into |tree| and returns a pointer to |item|.
+ If a duplicate item is found in the tree,
+ returns a pointer to the duplicate without inserting |item|.
+ Returns |NULL| in case of memory allocation failure. */
+struct pavlrc *pavlrc_probe(struct pavlrc_table *tree, struct pavlrc *item)
+{
+ struct pavlrc_node *y; /* Top node to update balance factor, and parent. */
+ struct pavlrc_node *p, *q; /* Iterator, and parent. */
+ struct pavlrc_node *n; /* Newly inserted node. */
+ struct pavlrc_node *w; /* New root of rebalanced subtree. */
+ int dir; /* Direction to descend. */
+
+ assert(tree != NULL && item != NULL);
+
+ y = p = tree->pavl_root;
+ q = NULL;
+ dir = 0;
+ while (p != NULL) {
+ int cmp = CMP_RC(item, &p->pavl_data);
+
+ if (cmp == 0)
+ return &p->pavl_data;
+
+ dir = cmp > 0;
+
+ if (p->pavl_balance != 0)
+ y = p;
+
+ q = p, p = p->pavl_link[dir];
+ }
+
+ n = tree->pavl_alloc->libavl_malloc(tree->pavl_alloc, sizeof *p);
+ if (n == NULL)
+ return NULL;
+
+ tree->pavl_count++;
+ n->pavl_link[0] = n->pavl_link[1] = NULL;
+ n->pavl_parent = q;
+ n->pavl_data = *item;
+ n->pavl_balance = 0;
+ if (q == NULL) {
+ tree->pavl_root = n;
+
+ return item /* &n->pavl_data */;
+ }
+ q->pavl_link[dir] = n;
+
+ p = n;
+ while (p != y) {
+ q = p->pavl_parent;
+ /*
+ dir = q->pavl_link[0] != p;
+ if (dir == 0)
+ q->pavl_balance--;
+ else
+ q->pavl_balance++;
+ */
+ if (q->pavl_link[0] != p)
+ q->pavl_balance++;
+ else
+ q->pavl_balance--;
+
+ p = q;
+ }
+
+ if (y->pavl_balance == -2) {
+ struct pavlrc_node *x = y->pavl_link[0];
+
+ if (x->pavl_balance == -1) {
+ w = x;
+ y->pavl_link[0] = x->pavl_link[1];
+ x->pavl_link[1] = y;
+ x->pavl_balance = y->pavl_balance = 0;
+ x->pavl_parent = y->pavl_parent;
+ y->pavl_parent = x;
+ if (y->pavl_link[0] != NULL)
+ y->pavl_link[0]->pavl_parent = y;
+ }
+ else {
+ assert(x->pavl_balance == +1);
+ w = x->pavl_link[1];
+ x->pavl_link[1] = w->pavl_link[0];
+ w->pavl_link[0] = x;
+ y->pavl_link[0] = w->pavl_link[1];
+ w->pavl_link[1] = y;
+ if (w->pavl_balance == -1)
+ x->pavl_balance = 0, y->pavl_balance = +1;
+ else if (w->pavl_balance == 0)
+ x->pavl_balance = y->pavl_balance = 0;
+ else /* |w->pavl_balance == +1| */
+ x->pavl_balance = -1, y->pavl_balance = 0;
+ w->pavl_balance = 0;
+ w->pavl_parent = y->pavl_parent;
+ x->pavl_parent = y->pavl_parent = w;
+ if (x->pavl_link[1] != NULL)
+ x->pavl_link[1]->pavl_parent = x;
+ if (y->pavl_link[0] != NULL)
+ y->pavl_link[0]->pavl_parent = y;
+ }
+ }
+ else if (y->pavl_balance == +2) {
+ struct pavlrc_node *x = y->pavl_link[1];
+
+ if (x->pavl_balance == +1) {
+ w = x;
+ y->pavl_link[1] = x->pavl_link[0];
+ x->pavl_link[0] = y;
+ x->pavl_balance = y->pavl_balance = 0;
+ x->pavl_parent = y->pavl_parent;
+ y->pavl_parent = x;
+ if (y->pavl_link[1] != NULL)
+ y->pavl_link[1]->pavl_parent = y;
+ }
+ else {
+ assert(x->pavl_balance == -1);
+ w = x->pavl_link[0];
+ x->pavl_link[0] = w->pavl_link[1];
+ w->pavl_link[1] = x;
+ y->pavl_link[1] = w->pavl_link[0];
+ w->pavl_link[0] = y;
+ if (w->pavl_balance == +1)
+ x->pavl_balance = 0, y->pavl_balance = -1;
+ else if (w->pavl_balance == 0)
+ x->pavl_balance = y->pavl_balance = 0;
+ else /* |w->pavl_balance == -1| */
+ x->pavl_balance = +1, y->pavl_balance = 0;
+ w->pavl_balance = 0;
+ w->pavl_parent = y->pavl_parent;
+ x->pavl_parent = y->pavl_parent = w;
+ if (x->pavl_link[0] != NULL)
+ x->pavl_link[0]->pavl_parent = x;
+ if (y->pavl_link[1] != NULL)
+ y->pavl_link[1]->pavl_parent = y;
+ }
+ }
+ else
+ return item /* &n->pavl_data */;
+
+ if (w->pavl_parent != NULL)
+ w->pavl_parent->pavl_link[y != w->pavl_parent->pavl_link[0]] = w;
+ else
+ tree->pavl_root = w;
+
+ return item /* &n->pavl_data */;
+}
+
+/* Inserts |item| into |table|.
+ Returns |NULL| if |item| was successfully inserted
+ or if a memory allocation error occurred.
+ Otherwise, returns the duplicate item. */
+struct pavlrc *pavlrc_insert(struct pavlrc_table *table, struct pavlrc *item)
+{
+ struct pavlrc *p = pavlrc_probe(table, item);
+
+ return p == NULL || p == item ? NULL : p;
+}
+
+/* Inserts |item| into |table|, replacing any duplicate item.
+ Returns |NULL| if |item| was inserted without replacing a duplicate,
+ or if a memory allocation error occurred.
+ Otherwise, returns the item that was replaced. */
+struct pavlrc *pavlrc_replace(struct pavlrc_table *table, struct pavlrc *item)
+{
+ struct pavlrc *p = pavlrc_probe(table, item);
+
+ if (p == NULL || p == item)
+ return NULL;
+ else {
+ struct pavlrc *r = p;
+
+ *p = *item;
+
+ return r;
+ }
+}
+
+/* Deletes from |tree| and returns an item matching |item|.
+ Returns a null pointer if no matching item found. */
+struct pavlrc *pavlrc_delete(struct pavlrc_table *tree, struct pavlrc *item)
+{
+ struct pavlrc_node *p; /* Traverses tree to find node to delete. */
+ struct pavlrc_node *q; /* Parent of |p|. */
+ int dir; /* Side of |q| on which |p| is linked. */
+ int cmp; /* Result of comparison between |item| and |p|. */
+
+ assert(tree != NULL && item != NULL);
+
+ p = tree->pavl_root;
+ dir = 0;
+ while (p != NULL) {
+ cmp = CMP_RC(item, &p->pavl_data);
+
+ if (cmp == 0)
+ break;
+
+ dir = cmp > 0;
+ p = p->pavl_link[dir];
+ }
+ if (p == NULL)
+ return NULL;
+
+ /* item = p->pavl_data; */
+
+ q = p->pavl_parent;
+ if (q == NULL) {
+ q = (struct pavlrc_node *)&tree->pavl_root;
+ dir = 0;
+ }
+
+ if (p->pavl_link[1] == NULL) {
+ q->pavl_link[dir] = p->pavl_link[0];
+ if (q->pavl_link[dir] != NULL)
+ q->pavl_link[dir]->pavl_parent = p->pavl_parent;
+ }
+ else {
+ struct pavlrc_node *r = p->pavl_link[1];
+
+ if (r->pavl_link[0] == NULL) {
+ r->pavl_link[0] = p->pavl_link[0];
+ q->pavl_link[dir] = r;
+ r->pavl_parent = p->pavl_parent;
+ if (r->pavl_link[0] != NULL)
+ r->pavl_link[0]->pavl_parent = r;
+ r->pavl_balance = p->pavl_balance;
+ q = r;
+ dir = 1;
+ }
+ else {
+ struct pavlrc_node *s = r->pavl_link[0];
+
+ while (s->pavl_link[0] != NULL)
+ s = s->pavl_link[0];
+ r = s->pavl_parent;
+ r->pavl_link[0] = s->pavl_link[1];
+ s->pavl_link[0] = p->pavl_link[0];
+ s->pavl_link[1] = p->pavl_link[1];
+ q->pavl_link[dir] = s;
+ if (s->pavl_link[0] != NULL)
+ s->pavl_link[0]->pavl_parent = s;
+ s->pavl_link[1]->pavl_parent = s;
+ s->pavl_parent = p->pavl_parent;
+ if (r->pavl_link[0] != NULL)
+ r->pavl_link[0]->pavl_parent = r;
+ s->pavl_balance = p->pavl_balance;
+ q = r;
+ dir = 0;
+ }
+ }
+ tree->pavl_alloc->libavl_free(tree->pavl_alloc, p);
+
+ while (q != (struct pavlrc_node *)&tree->pavl_root) {
+ struct pavlrc_node *y = q;
+
+ if (y->pavl_parent != NULL)
+ q = y->pavl_parent;
+ else
+ q = (struct pavlrc_node *)&tree->pavl_root;
+
+ if (dir == 0) {
+ dir = q->pavl_link[0] != y;
+ y->pavl_balance++;
+ if (y->pavl_balance == +1)
+ break;
+ else if (y->pavl_balance == +2) {
+ struct pavlrc_node *x = y->pavl_link[1];
+
+ if (x->pavl_balance == -1) {
+ struct pavlrc_node *w;
+
+ assert(x->pavl_balance == -1);
+ w = x->pavl_link[0];
+ x->pavl_link[0] = w->pavl_link[1];
+ w->pavl_link[1] = x;
+ y->pavl_link[1] = w->pavl_link[0];
+ w->pavl_link[0] = y;
+ if (w->pavl_balance == +1)
+ x->pavl_balance = 0, y->pavl_balance = -1;
+ else if (w->pavl_balance == 0)
+ x->pavl_balance = y->pavl_balance = 0;
+ else /* |w->pavl_balance == -1| */
+ x->pavl_balance = +1, y->pavl_balance = 0;
+ w->pavl_balance = 0;
+ w->pavl_parent = y->pavl_parent;
+ x->pavl_parent = y->pavl_parent = w;
+ if (x->pavl_link[0] != NULL)
+ x->pavl_link[0]->pavl_parent = x;
+ if (y->pavl_link[1] != NULL)
+ y->pavl_link[1]->pavl_parent = y;
+ q->pavl_link[dir] = w;
+ }
+ else {
+ y->pavl_link[1] = x->pavl_link[0];
+ x->pavl_link[0] = y;
+ x->pavl_parent = y->pavl_parent;
+ y->pavl_parent = x;
+ if (y->pavl_link[1] != NULL)
+ y->pavl_link[1]->pavl_parent = y;
+ q->pavl_link[dir] = x;
+ if (x->pavl_balance == 0) {
+ x->pavl_balance = -1;
+ y->pavl_balance = +1;
+ break;
+ }
+ else {
+ x->pavl_balance = y->pavl_balance = 0;
+ y = x;
+ }
+ }
+ }
+ }
+ else {
+ dir = q->pavl_link[0] != y;
+ y->pavl_balance--;
+ if (y->pavl_balance == -1)
+ break;
+ else if (y->pavl_balance == -2) {
+ struct pavlrc_node *x = y->pavl_link[0];
+
+ if (x->pavl_balance == +1) {
+ struct pavlrc_node *w;
+
+ assert(x->pavl_balance == +1);
+ w = x->pavl_link[1];
+ x->pavl_link[1] = w->pavl_link[0];
+ w->pavl_link[0] = x;
+ y->pavl_link[0] = w->pavl_link[1];
+ w->pavl_link[1] = y;
+ if (w->pavl_balance == -1)
+ x->pavl_balance = 0, y->pavl_balance = +1;
+ else if (w->pavl_balance == 0)
+ x->pavl_balance = y->pavl_balance = 0;
+ else /* |w->pavl_balance == +1| */
+ x->pavl_balance = -1, y->pavl_balance = 0;
+ w->pavl_balance = 0;
+ w->pavl_parent = y->pavl_parent;
+ x->pavl_parent = y->pavl_parent = w;
+ if (x->pavl_link[1] != NULL)
+ x->pavl_link[1]->pavl_parent = x;
+ if (y->pavl_link[0] != NULL)
+ y->pavl_link[0]->pavl_parent = y;
+ q->pavl_link[dir] = w;
+ }
+ else {
+ y->pavl_link[0] = x->pavl_link[1];
+ x->pavl_link[1] = y;
+ x->pavl_parent = y->pavl_parent;
+ y->pavl_parent = x;
+ if (y->pavl_link[0] != NULL)
+ y->pavl_link[0]->pavl_parent = y;
+ q->pavl_link[dir] = x;
+ if (x->pavl_balance == 0) {
+ x->pavl_balance = +1;
+ y->pavl_balance = -1;
+ break;
+ }
+ else {
+ x->pavl_balance = y->pavl_balance = 0;
+ y = x;
+ }
+ }
+ }
+ }
+ }
+
+ tree->pavl_count--;
+
+ return item;
+}
+
+/* Initializes |trav| for use with |tree|
+ and selects the null node. */
+void pavlrc_t_init(struct pavlrc_traverser *trav, struct pavlrc_table *tree)
+{
+ trav->pavl_table = tree;
+ trav->pavl_node = NULL;
+}
+
+/* Initializes |trav| for |tree|.
+ Returns data item in |tree| with the least value,
+ or |NULL| if |tree| is empty. */
+struct pavlrc *pavlrc_t_first(struct pavlrc_traverser *trav, struct pavlrc_table *tree)
+{
+ assert(tree != NULL && trav != NULL);
+
+ trav->pavl_table = tree;
+ trav->pavl_node = tree->pavl_root;
+ if (trav->pavl_node != NULL) {
+ while (trav->pavl_node->pavl_link[0] != NULL)
+ trav->pavl_node = trav->pavl_node->pavl_link[0];
+
+ return &trav->pavl_node->pavl_data;
+ }
+ else
+ return NULL;
+}
+
+/* Initializes |trav| for |tree|.
+ Returns data item in |tree| with the greatest value,
+ or |NULL| if |tree| is empty. */
+struct pavlrc *pavlrc_t_last(struct pavlrc_traverser *trav, struct pavlrc_table *tree)
+{
+ assert(tree != NULL && trav != NULL);
+
+ trav->pavl_table = tree;
+ trav->pavl_node = tree->pavl_root;
+ if (trav->pavl_node != NULL) {
+ while (trav->pavl_node->pavl_link[1] != NULL)
+ trav->pavl_node = trav->pavl_node->pavl_link[1];
+
+ return &trav->pavl_node->pavl_data;
+ }
+ else
+ return NULL;
+}
+
+/* Searches for |item| in |tree|.
+ If found, initializes |trav| to the item found and returns the item
+ as well.
+ If there is no matching item, initializes |trav| to the null item
+ and returns |NULL|. */
+struct pavlrc *pavlrc_t_find(struct pavlrc_traverser *trav, struct pavlrc_table *tree,
+ struct pavlrc *item)
+{
+ struct pavlrc_node *p;
+
+ assert(trav != NULL && tree != NULL && item != NULL);
+
+ trav->pavl_table = tree;
+
+ p = tree->pavl_root;
+ while (p != NULL) {
+ int cmp = CMP_RC(item, &p->pavl_data);
+
+ if (cmp == 0) {
+ trav->pavl_node = p;
+
+ return &p->pavl_data;
+ }
+
+ p = p->pavl_link[cmp > 0];
+ }
+
+ trav->pavl_node = NULL;
+
+ return NULL;
+}
+
+/* Attempts to insert |item| into |tree|.
+ If |item| is inserted successfully, it is returned and |trav| is
+ initialized to its location.
+ If a duplicate is found, it is returned and |trav| is initialized to
+ its location. No replacement of the item occurs.
+ If a memory allocation failure occurs, |NULL| is returned and |trav|
+ is initialized to the null item. */
+struct pavlrc *pavlrc_t_insert(struct pavlrc_traverser *trav,
+ struct pavlrc_table *tree, struct pavlrc *item)
+{
+ struct pavlrc *p;
+
+ assert(trav != NULL && tree != NULL && item != NULL);
+
+ p = pavlrc_probe(tree, item);
+ if (p != NULL) {
+ trav->pavl_table = tree;
+ trav->pavl_node = ((struct pavlrc_node *)((char *)p -
+ offsetof(struct pavlrc_node,
+ pavl_data)));
+
+ return p;
+ }
+ else {
+ pavlrc_t_init(trav, tree);
+
+ return NULL;
+ }
+}
+
+/* Initializes |trav| to have the same current node as |src|. */
+struct pavlrc *pavlrc_t_copy(struct pavlrc_traverser *trav,
+ const struct pavlrc_traverser *src)
+{
+ assert(trav != NULL && src != NULL);
+
+ trav->pavl_table = src->pavl_table;
+ trav->pavl_node = src->pavl_node;
+
+ return trav->pavl_node != NULL ? &trav->pavl_node->pavl_data : NULL;
+}
+
+/* Returns the next data item in inorder
+ within the tree being traversed with |trav|,
+ or if there are no more data items returns |NULL|. */
+struct pavlrc *pavlrc_t_next(struct pavlrc_traverser *trav)
+{
+ assert(trav != NULL);
+
+ if (trav->pavl_node == NULL)
+ return pavlrc_t_first(trav, trav->pavl_table);
+ else if (trav->pavl_node->pavl_link[1] == NULL) {
+ struct pavlrc_node *q, *p; /* Current node and its child. */
+
+ for (p = trav->pavl_node, q = p->pavl_parent;;
+ p = q, q = q->pavl_parent)
+ if (q == NULL || p == q->pavl_link[0]) {
+ trav->pavl_node = q;
+
+ return trav->pavl_node != NULL ?
+ &trav->pavl_node->pavl_data : NULL;
+ }
+ }
+ else {
+ trav->pavl_node = trav->pavl_node->pavl_link[1];
+ while (trav->pavl_node->pavl_link[0] != NULL)
+ trav->pavl_node = trav->pavl_node->pavl_link[0];
+
+ return &trav->pavl_node->pavl_data;
+ }
+}
+
+/* Returns the previous data item in inorder
+ within the tree being traversed with |trav|,
+ or if there are no more data items returns |NULL|. */
+struct pavlrc *pavlrc_t_prev(struct pavlrc_traverser *trav)
+{
+ assert(trav != NULL);
+
+ if (trav->pavl_node == NULL)
+ return pavlrc_t_last(trav, trav->pavl_table);
+ else if (trav->pavl_node->pavl_link[0] == NULL) {
+ struct pavlrc_node *q, *p; /* Current node and its child. */
+
+ for (p = trav->pavl_node, q = p->pavl_parent;;
+ p = q, q = q->pavl_parent)
+ if (q == NULL || p == q->pavl_link[1]) {
+ trav->pavl_node = q;
+
+ return trav->pavl_node != NULL ?
+ &trav->pavl_node->pavl_data : NULL;
+ }
+ }
+ else {
+ trav->pavl_node = trav->pavl_node->pavl_link[0];
+ while (trav->pavl_node->pavl_link[1] != NULL)
+ trav->pavl_node = trav->pavl_node->pavl_link[1];
+
+ return &trav->pavl_node->pavl_data;
+ }
+}
+
+/* Returns |trav|'s current item. */
+struct pavlrc *pavlrc_t_cur(struct pavlrc_traverser *trav)
+{
+ assert(trav != NULL);
+
+ return trav->pavl_node != NULL ? &trav->pavl_node->pavl_data : NULL;
+}
+
+/* Replaces the current item in |trav| by |new| and returns the item replaced.
+ |trav| must not have the null item selected.
+ The new item must not upset the ordering of the tree. */
+struct pavlrc *pavlrc_t_replace(struct pavlrc_traverser *trav, struct pavlrc *new)
+{
+ struct pavlrc *old;
+
+ assert(trav != NULL && trav->pavl_node != NULL && new != NULL);
+ old = &trav->pavl_node->pavl_data;
+ trav->pavl_node->pavl_data = *new;
+
+ return old;
+}
+
+/* Destroys |new| with |pavl_destroy (new, destroy)|,
+ first initializing right links in |new| that have
+ not yet been initialized at time of call. */
+static void
+copy_error_recovery(struct pavlrc_node *q,
+ struct pavlrc_table *new)
+{
+ assert(q != NULL && new != NULL);
+
+ for (;;) {
+ struct pavlrc_node *p = q;
+
+ q = q->pavl_parent;
+ if (q == NULL)
+ break;
+
+ if (p == q->pavl_link[0])
+ q->pavl_link[1] = NULL;
+ }
+
+ pavlrc_destroy(new);
+}
+
+/* Copies |org| to a newly created tree, which is returned.
+ If |copy != NULL|, each data item in |org| is first passed to |copy|,
+ and the return values are inserted into the tree;
+ |NULL| return values are taken as indications of failure.
+ On failure, destroys the partially created new tree,
+ applying |destroy|, if non-null, to each item in the new tree so far,
+ and returns |NULL|.
+ If |allocator != NULL|, it is used for allocation in the new tree.
+ Otherwise, the same allocator used for |org| is used. */
+struct pavlrc_table *pavlrc_copy(const struct pavlrc_table *org,
+ pavlrc_copy_func * copy,
+ struct libavl_allocator *allocator)
+{
+ struct pavlrc_table *new;
+ struct pavlrc_node *x;
+ struct pavlrc_node *y;
+
+ assert(org != NULL);
+ new = pavlrc_create(allocator != NULL ? allocator : org->pavl_alloc);
+ if (new == NULL)
+ return NULL;
+
+ new->pavl_count = org->pavl_count;
+ if (new->pavl_count == 0)
+ return new;
+
+ x = (struct pavlrc_node *)&org->pavl_root;
+ y = (struct pavlrc_node *)&new->pavl_root;
+ while (x != NULL) {
+ while (x->pavl_link[0] != NULL) {
+ y->pavl_link[0] =
+ new->pavl_alloc->libavl_malloc(new->pavl_alloc,
+ sizeof *y->pavl_link[0]);
+ if (y->pavl_link[0] == NULL) {
+ if (y != (struct pavlrc_node *)&new->pavl_root) {
+ /* y->pavl_data = NULL; */
+ y->pavl_link[1] = NULL;
+ }
+
+ copy_error_recovery(y, new);
+
+ return NULL;
+ }
+ y->pavl_link[0]->pavl_parent = y;
+
+ x = x->pavl_link[0];
+ y = y->pavl_link[0];
+ }
+ y->pavl_link[0] = NULL;
+
+ for (;;) {
+ y->pavl_balance = x->pavl_balance;
+ if (copy == NULL)
+ y->pavl_data = x->pavl_data;
+ else {
+ y->pavl_data = *(struct pavlrc *)copy(&x->pavl_data);
+ /*
+ if (y->pavl_data == NULL) {
+ y->pavl_link[1] = NULL;
+ copy_error_recovery(y, new);
+
+ return NULL;
+ }
+ */
+ }
+
+ if (x->pavl_link[1] != NULL) {
+ y->pavl_link[1] =
+ new->pavl_alloc->libavl_malloc(new->pavl_alloc,
+ sizeof *y->pavl_link[1]);
+ if (y->pavl_link[1] == NULL) {
+ copy_error_recovery(y, new);
+
+ return NULL;
+ }
+ y->pavl_link[1]->pavl_parent = y;
+
+ x = x->pavl_link[1];
+ y = y->pavl_link[1];
+ break;
+ }
+ else
+ y->pavl_link[1] = NULL;
+
+ for (;;) {
+ const struct pavlrc_node *w = x;
+
+ x = x->pavl_parent;
+ if (x == NULL) {
+ new->pavl_root->pavl_parent = NULL;
+ return new;
+ }
+ y = y->pavl_parent;
+
+ if (w == x->pavl_link[0])
+ break;
+ }
+ }
+ }
+
+ return new;
+}
+
+/* Frees storage allocated for |tree|.
+ If |destroy != NULL|, applies it to each data item in inorder. */
+void pavlrc_destroy(struct pavlrc_table *tree)
+{
+ struct pavlrc_node *p, *q;
+
+ assert(tree != NULL);
+
+ p = tree->pavl_root;
+ while (p != NULL) {
+ if (p->pavl_link[0] == NULL) {
+ q = p->pavl_link[1];
+ tree->pavl_alloc->libavl_free(tree->pavl_alloc, p);
+ }
+ else {
+ q = p->pavl_link[0];
+ p->pavl_link[0] = q->pavl_link[1];
+ q->pavl_link[1] = p;
+ }
+ p = q;
+ }
+
+ tree->pavl_alloc->libavl_free(tree->pavl_alloc, tree);
+}
+
+/* Allocates |size| bytes of space using |malloc()|.
+ Returns a null pointer if allocation fails. */
+void *pavlrc_malloc(struct libavl_allocator *allocator, size_t size)
+{
+ assert(allocator != NULL && size > 0);
+
+ return malloc(size);
+}
+
+/* Frees |block|. */
+void pavlrc_free(struct libavl_allocator *allocator, void *block)
+{
+ assert(allocator != NULL && block != NULL);
+
+ free(block);
+}
+
+/* Default memory allocator that uses |malloc()| and |free()|. */
+struct libavl_allocator pavlrc_allocator_default = {
+ pavlrc_malloc,
+ pavlrc_free
+};
+
+#undef NDEBUG
+#include <assert.h>
+
+/* Asserts that |pavl_insert()| succeeds at inserting |item| into |table|. */
+void (pavlrc_assert_insert) (struct pavlrc_table * table, struct pavlrc *item)
+{
+ struct pavlrc *p = pavlrc_probe(table, item);
+
+ assert(p != NULL && p == item);
+}
+
+/* Asserts that |pavl_delete()| really removes |item| from |table|,
+ and returns the removed item. */
+struct pavlrc *(pavlrc_assert_delete) (struct pavlrc_table * table, struct pavlrc *item)
+{
+ struct pavlrc *p = pavlrc_delete(table, item);
+
+ assert(p != NULL);
+
+ return p;
+}
Property changes on: grass/trunk/raster/r.path/pavlrc.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass/trunk/raster/r.path/pavlrc.h
===================================================================
--- grass/trunk/raster/r.path/pavlrc.h (rev 0)
+++ grass/trunk/raster/r.path/pavlrc.h 2017-12-20 22:06:06 UTC (rev 71962)
@@ -0,0 +1,107 @@
+/* Produced by texiweb from libavl.w. */
+
+/* libavl - library for manipulation of binary trees.
+ Copyright (C) 1998, 1999, 2000, 2001, 2002, 2004 Free Software
+ Foundation, Inc.
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 3 of the License, or (at your option) any later version.
+
+ This library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with this library; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ 02110-1301 USA.
+ */
+
+#ifndef PAVLRC_H
+#define PAVLRC_H 1
+
+#include <stddef.h>
+
+struct pavlrc
+{
+ int row, col;
+};
+
+/* Function types. */
+typedef void *pavlrc_copy_func(void *pavl_item);
+
+#ifndef LIBAVL_ALLOCATOR
+#define LIBAVL_ALLOCATOR
+/* Memory allocator. */
+struct libavl_allocator
+{
+ void *(*libavl_malloc) (struct libavl_allocator *, size_t libavl_size);
+ void (*libavl_free) (struct libavl_allocator *, void *libavl_block);
+};
+#endif
+
+/* Default memory allocator. */
+extern struct libavl_allocator pavlrc_allocator_default;
+void *pavlrc_malloc(struct libavl_allocator *, size_t);
+void pavlrc_free(struct libavl_allocator *, void *);
+
+/* Maximum PAVL height, unused. */
+#ifndef PAVL_MAX_HEIGHT
+#define PAVL_MAX_HEIGHT 32
+#endif
+
+/* Tree data structure. */
+struct pavlrc_table
+{
+ struct pavlrc_node *pavl_root; /* Tree's root. */
+ struct libavl_allocator *pavl_alloc; /* Memory allocator. */
+ size_t pavl_count; /* Number of items in tree. */
+};
+
+/* An PAVL tree node. */
+struct pavlrc_node
+{
+ struct pavlrc_node *pavl_link[2]; /* Subtrees. */
+ struct pavlrc_node *pavl_parent; /* Parent node. */
+ struct pavlrc pavl_data; /* data. */
+ signed char pavl_balance; /* Balance factor. */
+};
+
+/* PAVL traverser structure. */
+struct pavlrc_traverser
+{
+ struct pavlrc_table *pavl_table; /* Tree being traversed. */
+ struct pavlrc_node *pavl_node; /* Current node in tree. */
+};
+
+/* Table functions. */
+struct pavlrc_table *pavlrc_create(struct libavl_allocator *);
+struct pavlrc_table *pavlrc_copy(const struct pavlrc_table *, pavlrc_copy_func *,
+ struct libavl_allocator *);
+void pavlrc_destroy(struct pavlrc_table *);
+struct pavlrc *pavlrc_probe(struct pavlrc_table *, struct pavlrc *);
+struct pavlrc *pavlrc_insert(struct pavlrc_table *, struct pavlrc *);
+struct pavlrc *pavlrc_replace(struct pavlrc_table *, struct pavlrc *);
+struct pavlrc *pavlrc_delete(struct pavlrc_table *, struct pavlrc *);
+struct pavlrc *pavlrc_find(const struct pavlrc_table *, const struct pavlrc *);
+void pavlrc_assert_insert(struct pavlrc_table *, struct pavlrc *);
+struct pavlrc *pavlrc_assert_delete(struct pavlrc_table *, struct pavlrc *);
+
+#define pavlrc_count(table) ((size_t) (table)->pavlrc_count)
+
+/* Table traverser functions. */
+void pavlrc_t_init(struct pavlrc_traverser *, struct pavlrc_table *);
+struct pavlrc *pavlrc_t_first(struct pavlrc_traverser *, struct pavlrc_table *);
+struct pavlrc *pavlrc_t_last(struct pavlrc_traverser *, struct pavlrc_table *);
+struct pavlrc *pavlrc_t_find(struct pavlrc_traverser *, struct pavlrc_table *, struct pavlrc *);
+struct pavlrc *pavlrc_t_insert(struct pavlrc_traverser *, struct pavlrc_table *, struct pavlrc *);
+struct pavlrc *pavlrc_t_copy(struct pavlrc_traverser *, const struct pavlrc_traverser *);
+struct pavlrc *pavlrc_t_next(struct pavlrc_traverser *);
+struct pavlrc *pavlrc_t_prev(struct pavlrc_traverser *);
+struct pavlrc *pavlrc_t_cur(struct pavlrc_traverser *);
+struct pavlrc *pavlrc_t_replace(struct pavlrc_traverser *, struct pavlrc *);
+
+#endif /* pavlrc.h */
Property changes on: grass/trunk/raster/r.path/pavlrc.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Deleted: grass/trunk/raster/r.path/r.drain.html
===================================================================
--- grass/trunk/raster/r.drain/r.drain.html 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/r.drain.html 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,316 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>r.drain</em> traces a flow through a least-cost path in an
-elevation model or cost surface. For cost surfaces, a movement
-direction map must be specified with the
-<b>direction</b> option and the <b>-d</b> flag to trace a flow path following
-the given directions. Such a movement direction map can be generated with
-<em><a href="r.walk.html">r.walk</a></em>,
-<em><a href="r.cost.html">r.cost</a></em>,
-<em><a href="r.slope.aspect.html">r.slope.aspect</a></em> or
-<em><a href="r.watershed.html">r.watershed</a></em>
-provided that the direction is in degrees, measured counterclockwise
-from east.
-
-<p>
-The <b>output</b> raster map will show one or more least-cost paths
-between each user-provided location(s) and the minima (low category
-values) in the raster <b>input</b> map. If the <b>-d</b> flag is used
-the output least-cost paths will be found using the direction raster
-map. By default, the <b>output</b> will be an integer CELL map with
-category 1 along the least cost path, and null cells elsewhere.
-
-<p>With the <b>-c</b> (<em>copy</em>) flag, the input raster map cell values are
-copied verbatim along the path. With the <b>-a</b> (<em>accumulate</em>)
-flag, the accumulated cell value from the starting point up to the current
-cell is written on output. With either the <b>-c</b> or the <b>-a</b> flags, the
-<b>output</b> map is created with the same cell type as
-the <b>input</b> raster map (integer, float or double). With
-the <b>-n</b> (<em>number</em>) flag, the cells are numbered
-consecutively from the starting point to the final point.
-The <b>-c</b>, <b>-a</b>, and <b>-n</b> flags are mutually
-incompatible.
-
-<p>For an elevation surface, the path is calculated by choosing the
-steeper "slope" between adjacent cells. The slope calculation
-accurately acounts for the variable scale in lat-lon projections. For
-a cost surface, the path is calculated by following the movement
-direction surface back to the start point given
-in <em><a href="r.walk.html">r.walk</a></em> or
-<em><a href="r.cost.html">r.cost</a></em>. The path search stops
-as soon as a region border or a neighboring NULL cell is encountered,
-because in these cases the direction can not be determined (the path
-could continue outside the current region).
-
-<p>The <b>start_coordinates</b> parameter consists of map E and N grid coordinates of
-a starting point. Each x,y pair is the easting and northing (respectively) of
-a starting point from which a least-cost corridor will be developed.
-The <b>start_points</b> parameter can take multiple vector maps containing
-additional starting points.
-Up to 1024 starting points can be input from a combination of the
-<b>start_coordinates</b> and <b>start_points</b> parameters.
-
-<h3>Explanation of output values</h3>
-
-Consider the following example:
-
-<div class="code"><pre>
-Input: Output:
- ELEVATION SURFACE LEAST COST PATH
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 19. 20. 18. 19. 16. 15. 15. . . . . . . . .
-. . --- . . . . . . . . . . . . . . . . . . . . . . . . .
-. 20| 19| 17. 16. 17. 16. 16. . . 1 . 1 . 1 . . . .
-. . --- . . . . . . . . . . . . . . . . . . . . . . . . .
-. 18. 18. 24. 18. 15. 12. 11. . . . . . 1 . . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 22. 16. 16. 18. 10. 10. 10. . . . . . 1 . . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 17. 15. 15. 15. 10. 8 . 8 . . . . . . . 1 . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 24. 16. 8 . 7 . 8 . 0 . 12. . . . . . . 1 . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 17. 9 . 8 . 7 . 8 . 6 . 12. . . . . . . . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-</pre></div>
-
-<p>
-The user-provided starting location in the above example is
-the boxed <b>19</b> in the left-hand map. The path in the output
-shows the least-cost corridor for moving from the starting
-box to the lowest (smallest) possible point. This is the path a raindrop
-would take in this landscape.
-<p>
-With the <b>-c</b> <em>(copy)</em> flag, you get the following result:
-
-<div class="code"><pre>
-Input: Output:
- ELEVATION SURFACE LEAST COST PATH
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 19. 20. 18. 19. 16. 15. 15. . . . . . . . .
-. . --- . . . . . . . . . . . . . . . . . . . . . . . . .
-. 20| 19| 17. 16. 17. 16. 16. . . 19. 17. 16. . . .
-. . --- . . . . . . . . . . . . . . . . . . . . . . . . .
-. 18. 18. 24. 18. 15. 12. 11. . . . . . 15. . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 22. 16. 16. 18. 10. 10. 10. . . . . . 10. . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 17. 15. 15. 15. 10. 8 . 8 . . . . . . . 8 . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 24. 16. 8 . 7 . 8 . 0 .12 . . . . . . . 0 . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 17. 9 . 8 . 7 . 8 . 6 .12 . . . . . . . . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-
-Note that the last <em>0</em> will not be put in the null values map.
-</pre></div>
-
-<p>
-With the <b>-a</b> <em>(accumulate)</em> flag, you get the following result:
-
-<div class="code"><pre>
-Input: Output:
- ELEVATION SURFACE LEAST COST PATH
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 19. 20. 18. 19. 16. 15. 15. . . . . . . . .
-. . --- . . . . . . . . . . . . . . . . . . . . . . . . .
-. 20| 19| 17. 16. 17. 16. 16. . . 19. 36. 52. . . .
-. . --- . . . . . . . . . . . . . . . . . . . . . . . . .
-. 18. 18. 24. 18. 15. 12. 11. . . . . . 67. . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 22. 16. 16. 18. 10. 10. 10. . . . . . 77. . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 17. 15. 15. 15. 10. 8 . 8 . . . . . . . 85. .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 24. 16. 8 . 7 . 8 . 0 .12 . . . . . . . 85. .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 17. 9 . 8 . 7 . 8 . 6 .12 . . . . . . . . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-</pre></div>
-
-<p>
-With the <b>-n</b> <em>(number)</em> flag, you get the following result:
-
-<div class="code"><pre>
-Input: Output:
- ELEVATION SURFACE LEAST COST PATH
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 19. 20. 18. 19. 16. 15. 15. . . . . . . . .
-. . --- . . . . . . . . . . . . . . . . . . . . . . . . .
-. 20| 19| 17. 16. 17. 16. 16. . . 1 . 2 . 3 . . . .
-. . --- . . . . . . . . . . . . . . . . . . . . . . . . .
-. 18. 18. 24. 18. 15. 12. 11. . . . . . 4 . . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 22. 16. 16. 18. 10. 10. 10. . . . . . 5 . . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 17. 15. 15. 15. 10. 8 . 8 . . . . . . . 6 . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 24. 16. 8 . 7 . 8 . 0 .12 . . . . . . . 7 . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-. 17. 9 . 8 . 7 . 8 . 6 .12 . . . . . . . . .
-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-</pre></div>
-
-With the <b>-d</b> <em>(direction)</em> flag, the direction raster is used
-for the input, rather than the elevation surface. The output is then created
-according to one of the <b>-can</b> flags.
-<div class="code"><pre>
-The directions are recorded as degrees CCW from East:
- 112.5 67.5 i.e. a cell with the value 135
-157.5 135 90 45 22.5 means the next cell is to the North-West
- 180 x 0
-202.5 225 270 315 337.5
- 247.5 292.5
-</pre></div>
-
-<h2>NOTES</h2>
-
-If no direction input map is given, <em>r.drain</em> currently finds
-only the lowest point (the cell having the smallest category value) in
-the input file that can be reached through directly adjacent cells
-that are less than or equal in value to the cell reached immediately
-prior to it; therefore, it will not necessarily reach the lowest point
-in the input file. It currently finds <em>pits</em> in the data,
-rather than the lowest point in the entire input
-map. The <em><a href="r.fill.dir.html">r.fill.dir</a></em>,
-<em><a href="r.terraflow.html">r.terraflow</a></em>,
-and <em><a href="r.basins.fill.html">r.basins.fill</a></em> modules
-can be used to fill in subbasins prior to processing
-with <em>r.drain</em>.
-
-<p>
-<em>r.drain</em> will not give sane results at the region boundary. On outer rows
-and columns bordering the edge of the region, the flow direction is always directly out
-of the map. In this case, the user could try adjusting the region extents slightly with
-<em>g.region</em> to allow additional outlet paths for <em>r.drain</em>.
-
-<h2>EXAMPLES</h2>
-
-<h3>Path to the lowest point</h3>
-
-In this example we compute drainage paths from two given points
-following decreasing elevation values to the lowest point.
-We are using the full North Carolina sample dataset.
-First we create the two points from a text file using
-<em><a href="v.in.ascii.html">v.in.ascii</a></em> module
-(here the text file is CSV and we are using unix here-file syntax
-with EOF, in GUI just enter the values directly for the parameter input):
-
-<div class="code"><pre>
-v.in.ascii input=- output=start format=point separator=comma <<EOF
-638667.15686275,220610.29411765
-638610.78431373,220223.03921569
-EOF
-</pre></div>
-
-Now we compute the drainage path:
-
-<div class="code"><pre>
-r.drain input=elev_lid792_1m output=drain_path drain=drain start_points=start
-</pre></div>
-
-Before we visualize the result, we set a color table for the elevation
-we are using and we create a shaded relief map:
-
-<div class="code"><pre>
-r.colors map=elev_lid792_1m color=elevation
-r.relief input=elev_lid792_1m output=relief
-</pre></div>
-
-Finally we visualize all the input and output data:
-
-<div class="code"><pre>
-d.shade shade=relief color=elev_lid792_1m
-d.vect map=drain_path color=0:0:61 width=4 legend_label="drainage paths"
-d.vect map=start color=none fill_color=224:0:0 icon=basic/circle size=15 legend_label=origins
-d.legend.vect -b
-</pre></div>
-
-<div align="center">
-<a href="r_drain.png"><img src="r_drain.png" alt="drainage using r.watershed" width="300" height="280"></a>
-<br>
-<i>Figure: Drainage paths from two points flowing into the points with
-lowest values</i>
-</div>
-
-<h3>Path following directions</h3>
-
-To continue flow even after it hits a depression, we need to supply a
-direction raster map which will tell the <em>r.drain</em> module how to
-continue from the depression. To get these directions, we use the
-<em><a href="r.watershed.html">r.watershed</a></em> module:
-
-<div class="code"><pre>
-r.watershed elevation=elev_lid792_1m accumulation=accum drainage=drain_dir
-</pre></div>
-
-The directions are categorical and we convert them to degrees using
-raster algebra:
-
-<div class="code"><pre>
-r.mapcalc "drain_deg = if(drain_dir != 0, 45. * abs(drain_dir), null())"
-</pre></div>
-
-Together with directions, we need to provide the <em>r.drain</em> module
-with cost values. We don't have any cost to assign to specific cells,
-so we create a constant surface:
-
-<div class="code"><pre>
-r.mapcalc "const1 = 1"
-</pre></div>
-
-Now we are ready to compute the drainage paths.
-We are using the two points from the previous example.
-
-<div class="code"><pre>
-r.drain -d input=const1 direction=drain_deg output=drain_path_2 drain=drain_2 start_points=start
-</pre></div>
-
-We visualize the result in the same way as in the previous example.
-
-<div align="center">
-<a href="r_drain_with_r_watershed_direction.png"><img src="r_drain_with_r_watershed_direction.png" alt="drainage using r.watershed" width="300" height="280"></a>
-<br>
-<i>Figure: Drainage paths from two points where directions from
-r.watershed were used</i>
-</div>
-
-<h2>KNOWN ISSUES</h2>
-
-<p>Sometimes, when the differences among integer cell category values in the
-<em><a href="r.cost.html">r.cost</a></em> cumulative cost surface output are
-small, this cumulative cost output cannot accurately be used as input to
-<em>r.drain</em> (<em>r.drain</em> will output bad results).
-This problem can be circumvented by making the differences
-between cell category values in the cumulative cost output bigger. It
-is recommended that if the output from <em>r.cost</em> is to be used
-as input to <em>r.drain</em>, the user multiply the <em>r.cost</em>
-input cost surface map by the value of the map's cell resolution,
-before running <em>r.cost</em>. This can be done using
-<em><a href="r.mapcalc.html">r.mapcalc</a></em>. The map resolution can be
-found using <em><a href="g.region.html">g.region</a></em>.
-This problem doesn't arise with floating point maps.
-
-
-<h2>SEE ALSO</h2>
-
-<em>
-<a href="g.region.html">g.region</a>,
-<a href="r.cost.html">r.cost</a>,
-<a href="r.fill.dir.html">r.fill.dir</a>,
-<a href="r.basins.fill.html">r.basins.fill</a>,
-<a href="r.watershed.html">r.watershed</a>,
-<a href="r.terraflow.html">r.terraflow</a>,
-<a href="r.mapcalc.html">r.mapcalc</a>,
-<a href="r.walk.html">r.walk</a>
-</em>
-
-<h2>AUTHORS</h2>
-
-Completely rewritten by Roger S. Miller, 2001<br>
-July 2004 at WebValley 2004, error checking and vector points added by
-Matteo Franchi (Liceo Leonardo Da Vinci, Trento) and
-Roberto Flor (ITC-irst, Trento, Italy)
-
-<p>
-<i>Last changed: $Date$</i>
Copied: grass/trunk/raster/r.path/r.path.html (from rev 71888, grass/trunk/raster/r.drain/r.drain.html)
===================================================================
--- grass/trunk/raster/r.path/r.path.html (rev 0)
+++ grass/trunk/raster/r.path/r.path.html 2017-12-20 22:06:06 UTC (rev 71962)
@@ -0,0 +1,192 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.path</em> traces a path from starting points following input
+directions. Such a movement direction map can be generated with
+<em><a href="r.walk.html">r.walk</a></em>,
+<em><a href="r.cost.html">r.cost</a></em>,
+<em><a href="r.slope.aspect.html">r.slope.aspect</a></em>,
+<em><a href="r.watershed.html">r.watershed</a></em>, or
+<em><a href="r.fill.dir.html">r.fill.dir</a></em>,
+provided that the direction is in degrees, measured counterclockwise
+from east.
+
+<p>
+Alternatively, bitmask-encoded directions can be provided where each
+bit position corresponds to a specific neighbour. A path will continue
+to all neighbours with their bit set. This means a path can split and
+merge. Such bitmasked directions can be created with the <b>-b</b>
+flag of <em><a href="r.cost.html">r.cost</a></em> and
+<em><a href="r.walk.html">r.walk</a></em>. The corresponding <b>-b</b>
+flag must then be set for <em>r.path</em>.
+
+<div class="code"><pre>
+Direction encoding for neighbors of x
+
+ 135 90 45 7 8 1
+ 180 x 360 6 x 2
+ 225 270 315 5 4 3
+
+ degrees bit positions
+ CCW from East
+</pre></div>
+
+A path stops when the direction is zero or negative, indicating a stop
+point or outlet.
+<p>
+The <b>output</b> raster map will show one or more least-cost paths
+between each user-provided location(s) and the target points (direction
+≤ 0). By default, the <b>output</b> will be an integer CELL map with
+the id of the start points along the least cost path, and null cells elsewhere.
+
+<p>
+With the <b>-c</b> (<em>copy</em>) flag, the values raster map cell values are
+copied verbatim along the path. With the <b>-a</b> (<em>accumulate</em>)
+flag, the accumulated cell value from the starting point up to the current
+cell is written on output. With either the <b>-c</b> or the <b>-a</b> flags, the
+<b>raster_path</b> map is created with the same cell type as
+the <b>values</b> raster map (integer, float or double). With
+the <b>-n</b> (<em>number</em>) flag, the cells are numbered
+consecutively from the starting point to the final point.
+The <b>-c</b>, <b>-a</b>, and <b>-n</b> flags are mutually
+incompatible.
+
+
+<p>
+The <b>start_coordinates</b> parameter consists of map E and N grid
+coordinates of a starting point. Each x,y pair is the easting and
+northing (respectively) of a starting point from which a path will be
+traced following <b>input</b> directions. The <b>start_points</b>
+parameter can take multiple vector maps containing additional starting
+points.
+
+<h2>NOTES</h2>
+The directions are recorded as degrees CCW from East, the Knight's move
+of r.cost and r.walk is considered:
+<div class="code"><pre>
+ 112.5 67.5
+157.5 135 90 45 22.5
+ 180 x 0
+202.5 225 270 315 337.5
+ 247.5 292.5
+</pre></div>
+i.e. a cell with the value 135 means the next cell is to the North-West,
+and a cell with the value 157.5 means that the next cell is to the
+West-North-West-
+
+
+<h2>EXAMPLES</h2>
+
+<h3>Hydrological path</h3>
+
+We are using the full North Carolina sample dataset.
+First we create the two points from a text file using
+<em><a href="v.in.ascii.html">v.in.ascii</a></em> module
+(here the text file is CSV and we are using unix here-file syntax
+with EOF, in GUI just enter the values directly for the parameter input):
+
+<div class="code"><pre>
+v.in.ascii input=- output=start format=point separator=comma <<EOF
+638667.15686275,220610.29411765
+638610.78431373,220223.03921569
+EOF
+</pre></div>
+
+We need to supply a direction raster map to the <em>r.path</em> module.
+To get these directions, we use the
+<em><a href="r.watershed.html">r.watershed</a></em> module:
+
+<div class="code"><pre>
+r.watershed elevation=elev_lid792_1m accumulation=accum drainage=drain_dir
+</pre></div>
+
+The directions are categorical and we convert them to degrees using
+raster algebra:
+
+<div class="code"><pre>
+r.mapcalc "drain_deg = if(drain_dir != 0, 45. * abs(drain_dir), null())"
+</pre></div>
+
+Now we are ready to extract the drainage paths starting at the two points.
+
+<div class="code"><pre>
+r.path input=drain_deg raster_path=drain_path vector_path=drain_path start_points=start
+</pre></div>
+
+Before we visualize the result, we set a color table for the elevation
+we are using and create a shaded relief map:
+
+<div class="code"><pre>
+r.colors map=elev_lid792_1m color=elevation
+r.relief input=elev_lid792_1m output=relief
+</pre></div>
+
+We visualize the input and output data:
+
+<div class="code"><pre>
+d.shade shade=relief color=elev_lid792_1m
+d.vect map=drain_path color=0:0:61 width=4 legend_label="drainage paths"
+d.vect map=start color=none fill_color=224:0:0 icon=basic/circle size=15 legend_label=origins
+d.legend.vect -b
+</pre></div>
+
+<div align="center">
+<a href="r_drain_with_r_watershed_direction.png"><img src="r_drain_with_r_watershed_direction.png" alt="drainage using r.watershed" width="300" height="280"></a>
+<br>
+<i>Figure: Drainage paths from two points where directions from
+r.watershed were used</i>
+</div>
+
+<h3>Least-cost path</h3>
+
+We compute bitmask encoded movement directions using <em>r.walk:</em>
+<div class="code"><pre>
+g.region swwake_30m -p
+
+# create friction map based on land cover
+r.recode landclass96 out=friction << EOF
+1:3:0.1:0.1
+4:5:10.:10.
+6:6:1000.0:1000.0
+7:7:0.3:0.3
+EOF
+
+# without Knight's move
+r.walk -b elevation=elev_ned_30m friction=friction output=walkcost \
+ outdir=walkdir start_coordinates=635576,216485
+
+r.path -b input=walkdir start_coordinates=640206,222795 \
+ raster_path=walkpath vector_path=walkpath
+
+# with Knight's move
+r.walk -b -k elevation=elev_ned_30m friction=friction output=walkcost_k \
+ outdir=walkdir_k start_coordinates=635576,216485
+
+r.path -b input=walkdir_k start_coordinates=640206,222795 \
+ raster_path=walkpath_k vector_path=walkpath_k
+</pre></div>
+
+The extracted least-cost path splits and merges on the way from
+the start point to the stop point (start point for r.walk). Note the
+gaps in the raster path when using the Knight's move.
+<div class="code"><pre>
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="g.region.html">g.region</a>,
+<a href="r.cost.html">r.cost</a>,
+<a href="r.walk.html">r.walk</a>,
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.fill.dir.html">r.fill.dir</a>,
+<a href="r.basins.fill.html">r.basins.fill</a>,
+<a href="r.terraflow.html">r.terraflow</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+</em>
+
+<h2>AUTHORS</h2>
+
+Markus Metz
+
+<p>
+<i>Last changed: $Date$</i>
Deleted: grass/trunk/raster/r.path/r_drain.png
===================================================================
(Binary files differ)
Deleted: grass/trunk/raster/r.path/r_drain_with_r_watershed_direction.png
===================================================================
(Binary files differ)
Copied: grass/trunk/raster/r.path/r_path.png (from rev 71888, grass/trunk/raster/r.drain/r_drain.png)
===================================================================
(Binary files differ)
Copied: grass/trunk/raster/r.path/r_path_with_r_watershed_direction.png (from rev 71888, grass/trunk/raster/r.drain/r_drain_with_r_watershed_direction.png)
===================================================================
(Binary files differ)
Deleted: grass/trunk/raster/r.path/resolve.c
===================================================================
--- grass/trunk/raster/r.drain/resolve.c 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/resolve.c 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,209 +0,0 @@
-#include <grass/config.h>
-#include <stdlib.h>
-#include <string.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include "tinf.h"
-
-CELL select_dir(CELL i)
-{
- CELL dir[256] = { 0, 1, 2, 2, 4, 1, 2, 2, 8, 1,
- 8, 2, 8, 4, 4, 2, 16, 16, 16, 2, 16, 4, 4,
- 2, 8, 8, 8, 8, 8, 8, 8, 4, 32, 1, 2, 2,
- 4, 4, 2, 2, 32, 8, 8, 2, 8, 8, 4, 4, 32,
- 32, 32, 32, 16, 32, 4, 2, 16, 16, 16, 16, 8, 16,
- 8, 8, 64, 64, 64, 1, 64, 1, 2, 2, 64, 64, 8,
- 2, 8, 8, 4, 2, 16, 64, 64, 2, 16, 64, 2, 2,
- 16, 8, 8, 8, 8, 8, 8, 4, 32, 64, 32, 1, 32,
- 32, 32, 2, 32, 32, 32, 2, 32, 8, 4, 4, 32, 32,
- 32, 32, 32, 32, 32, 32, 32, 32, 16, 16, 16, 16, 8,
- 8, 128, 128, 128, 1, 4, 1, 2, 2, 128, 128, 2, 1,
- 8, 4, 4, 2, 16, 128, 2, 1, 4, 128, 2, 1, 8,
- 128, 8, 1, 8, 8, 4, 2, 32, 128, 1, 1, 128, 128,
- 2, 1, 32, 128, 32, 1, 8, 128, 4, 2, 32, 32, 32,
- 1, 32, 128, 32, 1, 16, 16, 16, 1, 16, 16, 8, 4,
- 128, 128, 128, 128, 128, 128, 2, 1, 128, 128, 128, 1, 128,
- 128, 4, 2, 64, 128, 128, 1, 128, 128, 128, 1, 8, 128,
- 8, 1, 8, 8, 8, 2, 64, 128, 64, 128, 64, 128, 64,
- 128, 32, 64, 64, 128, 64, 64, 64, 1, 32, 64, 64, 128,
- 64, 64, 64, 128, 32, 32, 32, 64, 32, 32, 16, 128
- };
-
- return dir[i];
-}
-
-void flink(int i, int j, int nl, int ns, CELL * p1, CELL * p2, CELL * p3,
- int *active, int *goagain)
-{
- CELL bitmask[8] = { 1, 2, 4, 8, 16, 32, 64, 128 };
- CELL outflow, cwork, c[8];
- int k;
-
- cwork = p2[j];
- if (Rast_is_c_null_value(p2 + j) || cwork >= 0 || cwork == -256)
- return;
- cwork = -cwork;
-
- for (k = 7; k >= 0; k--) {
- c[k] = 0;
- if (cwork >= bitmask[k]) {
- c[k] = 1;
- cwork -= bitmask[k];
- }
- }
-
- outflow = 0;
-
- /* There's no need to resolve directions at cells adjacent to null cells,
- * as those directions will be resolved before we get here */
- /* look one row back */
- cwork = p1[j - 1];
- if (cwork > 0 && cwork != 4 && c[6])
- outflow += 64;
- cwork = p1[j];
- if (cwork > 0 && cwork != 8 && c[7])
- outflow += 128;
- cwork = p1[j + 1];
- if (cwork > 0 && cwork != 16 && c[0])
- outflow += 1;
-
- /* look at this row */
- cwork = p2[j - 1];
- if (cwork > 0 && cwork != 2 && c[5])
- outflow += 32;
- cwork = p2[j + 1];
- if (cwork > 0 && cwork != 32 && c[1])
- outflow += 2;
-
- /* look one row forward */
- cwork = p3[j - 1];
- if (cwork > 0 && cwork != 1 && c[4])
- outflow += 16;
- cwork = p3[j];
- if (cwork > 0 && cwork != 128 && c[3])
- outflow += 8;
- cwork = p3[j + 1];
- if (cwork > 0 && cwork != 64 && c[2])
- outflow += 4;
-
- if (outflow == 0) {
- *active = 1;
- }
- else {
- *goagain = 1;
- p2[j] = select_dir(outflow);
- }
- return;
-}
-
-void resolve(int fd, int nl, struct band3 *bnd)
-{
- CELL cvalue;
- int *active;
- int offset, isz, i, j, pass, activity, goagain, done;
-
- active = (int *)G_calloc(nl, sizeof(int));
-
- isz = sizeof(CELL);
-
- /* select a direction when there are multiple non-flat links */
- lseek(fd, bnd->sz, SEEK_SET);
- for (i = 1; i < nl - 1; i += 1) {
- read(fd, bnd->b[0], bnd->sz);
- for (j = 1; j < bnd->ns - 1; j += 1) {
- offset = j * isz;
- if (Rast_is_c_null_value((void *)(bnd->b[0] + offset)))
- continue;
- memcpy(&cvalue, bnd->b[0] + offset, isz);
- if (cvalue > 0)
- cvalue = select_dir(cvalue);
- memcpy(bnd->b[0] + offset, &cvalue, isz);
- }
- lseek(fd, -bnd->sz, SEEK_CUR);
- write(fd, bnd->b[0], bnd->sz);
- }
-
- pass = 0;
- for (i = 1; i < nl - 1; i += 1)
- active[i] = 1;
-
- /* select a direction when there are multiple flat links */
- do {
- done = 1;
- pass += 1;
-
- activity = 0;
-
- lseek(fd, 0, SEEK_SET);
- advance_band3(fd, bnd);
- advance_band3(fd, bnd);
- for (i = 1; i < nl - 1; i++) {
- lseek(fd, (off_t) (i + 1) * bnd->sz, SEEK_SET);
- advance_band3(fd, bnd);
-
- if (!active[i])
- continue;
-
- done = 0;
- active[i] = 0;
- do {
- goagain = 0;
- for (j = 1; j < bnd->ns - 1; j += 1) {
- flink(i, j, nl, bnd->ns,
- (void *)bnd->b[0], (void *)bnd->b[1],
- (void *)bnd->b[2], &active[i], &goagain);
- if (goagain)
- activity = 1;
- }
- } while (goagain);
-
- lseek(fd, (off_t) i * bnd->sz, SEEK_SET);
- write(fd, bnd->b[1], bnd->sz);
-
- }
-
- if (!activity) {
- done = 1;
- continue;
- }
-
- activity = 0;
-
- lseek(fd, (off_t) (nl - 1) * bnd->sz, SEEK_SET);
- retreat_band3(fd, bnd);
- retreat_band3(fd, bnd);
- for (i = nl - 2; i >= 1; i -= 1) {
- lseek(fd, (off_t) (i - 1) * bnd->sz, SEEK_SET);
- retreat_band3(fd, bnd);
-
- if (!active[i])
- continue;
-
- done = 0;
- active[i] = 0;
- do {
- goagain = 0;
- for (j = 1; j < bnd->ns - 1; j++) {
- flink(i, j, nl, bnd->ns,
- (void *)bnd->b[0], (void *)bnd->b[1],
- (void *)bnd->b[2], &active[i], &goagain);
- if (goagain)
- activity = 1;
- }
- } while (goagain);
-
- lseek(fd, (off_t) i * bnd->sz, SEEK_SET);
- write(fd, bnd->b[1], bnd->sz);
- }
-
- if (!activity) {
- done = 1;
- }
-
- } while (!done);
-
- G_free(active);
-
- return;
-
-}
Deleted: grass/trunk/raster/r.path/tests/test.r.drain.sh
===================================================================
--- grass/trunk/raster/r.drain/tests/test.r.drain.sh 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/tests/test.r.drain.sh 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,24 +0,0 @@
-#!/bin/sh
-
-# to be run in the North Carolina location
-
-export GRASS_OVERWRITE=1
-
-r.in.ascii testascii_nc.asc out=testascii
-g.region raster=testascii -p
-
-d.mon wx0
-sleep 2
-d.rast.num testascii
-
-r.drain input=testascii output=drain coord=638442.5,220548.5
-d.rast drain
-
-d.erase
-r.mapcalc "testascii_float = float(testascii)"
-r.drain input=testascii_float output=drain coord=638442.5,220548.5
-d.rast drain
-
-# for statistical test results:
-# - verify number of resulting pixels (r.univar -g)
-# - verify minimum and maximum of result (r.univar -g)
Deleted: grass/trunk/raster/r.path/tests/testascii_nc.asc
===================================================================
--- grass/trunk/raster/r.drain/tests/testascii_nc.asc 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/tests/testascii_nc.asc 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,13 +0,0 @@
-north: 220600
-south: 220000
-west: 638300
-east: 639000
-rows: 6
-cols: 7
-null: -9999
-20 19 17 16 17 16 16
-18 18 24 18 15 12 11
-22 16 16 18 10 10 10
-17 15 15 15 10 8 8
-24 16 8 7 8 0 12
-17 9 8 7 8 6 12
Deleted: grass/trunk/raster/r.path/tinf.c
===================================================================
--- grass/trunk/raster/r.drain/tinf.c 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/tinf.c 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,409 +0,0 @@
-#include <grass/config.h>
-/* #include <limits.h> */
-#include <float.h>
-#include <math.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include "tinf.h"
-
-int (*is_null) (void *);
-void (*set_null_value) (void *, int);
-int (*bpe) ();
-void *(*get_max) (void *, void *);
-void *(*get_min) (void *, void *);
-void (*get_row) (int, void *, int);
-void *(*get_buf) ();
-void (*put_row) (int, void *);
-double (*slope) (void *, void *, double);
-void (*set_min) (void *);
-void (*set_max) (void *);
-void (*diff) (void *, void *);
-void (*sum) (void *, void *);
-void (*quot) (void *, void *);
-void (*prod) (void *, void *);
-
-/* To add a new multitype function, use the function below to initialize
- * the function pointer to each of the three typed functions. The function
- * pointers and the function prototypes are defined in a header file.
- * The actual functions follow. */
-
-void set_func_pointers(int in_type)
-{
- switch (in_type) {
- case CELL_TYPE:
- is_null = is_null_c;
- bpe = bpe_c;
- get_max = get_max_c;
- get_min = get_min_c;
- get_row = get_row_c;
- get_buf = get_buf_c;
- put_row = put_row_c;
- slope = slope_c;
- set_min = set_min_c;
- set_max = set_max_c;
- diff = diff_c;
- sum = sum_c;
- quot = quot_c;
- prod = prod_c;
- set_null_value = set_null_value_c;
-
- break;
-
- case FCELL_TYPE:
- is_null = is_null_f;
- bpe = bpe_f;
- get_max = get_max_f;
- get_min = get_min_f;
- get_row = get_row_f;
- get_buf = get_buf_f;
- put_row = put_row_f;
- slope = slope_f;
- set_min = set_min_f;
- set_max = set_max_f;
- diff = diff_f;
- sum = sum_f;
- quot = quot_f;
- prod = prod_f;
- set_null_value = set_null_value_f;
-
- break;
-
- case DCELL_TYPE:
- is_null = is_null_d;
- bpe = bpe_d;
- get_max = get_max_d;
- get_min = get_min_d;
- get_row = get_row_d;
- get_buf = get_buf_d;
- put_row = put_row_d;
- slope = slope_d;
- set_min = set_min_d;
- set_max = set_max_d;
- diff = diff_d;
- sum = sum_d;
- quot = quot_d;
- prod = prod_d;
- set_null_value = set_null_value_d;
-
- }
-
- return;
-
-}
-
-/* check for null values */
-int is_null_c(void *value)
-{
- return Rast_is_c_null_value((CELL *) value);
-}
-int is_null_f(void *value)
-{
- return Rast_is_f_null_value((FCELL *) value);
-}
-int is_null_d(void *value)
-{
- return Rast_is_d_null_value((DCELL *) value);
-}
-
-/* set null values in buffer */
-void set_null_value_c(void *value, int num)
-{
- Rast_set_c_null_value((CELL *) value, num);
-}
-void set_null_value_f(void *value, int num)
-{
- Rast_set_f_null_value((FCELL *) value, num);
-}
-void set_null_value_d(void *value, int num)
-{
- Rast_set_d_null_value((DCELL *) value, num);
-}
-
-/* return the size of the current type */
-int bpe_c()
-{
- return sizeof(CELL);
-}
-
-int bpe_f()
-{
- return sizeof(FCELL);
-}
-
-int bpe_d()
-{
- return sizeof(DCELL);
-}
-
-/* return the pointer that points to the smaller of two value */
-void *get_min_c(void *v1, void *v2)
-{
- void *rc;
-
- rc = v2;
- if (*(CELL *) v1 < *(CELL *) v2)
- rc = v1;
- return rc;
-}
-
-void *get_min_f(void *v1, void *v2)
-{
- void *rc;
-
- rc = v2;
- if (*(FCELL *) v1 < *(FCELL *) v2)
- rc = v1;
- return rc;
-}
-
-void *get_min_d(void *v1, void *v2)
-{
- void *rc;
-
- rc = v2;
- if (*(DCELL *) v1 < *(DCELL *) v2)
- rc = v1;
- return rc;
-}
-
-/* return the pointer that points to the larger value */
-void *get_max_c(void *v1, void *v2)
-{
- void *rc;
-
- rc = v2;
- if (*(CELL *) v1 > *(CELL *) v2)
- rc = v1;
- return rc;
-}
-
-void *get_max_f(void *v1, void *v2)
-{
- void *rc;
-
- rc = v2;
- if (*(FCELL *) v1 > *(FCELL *) v2)
- rc = v1;
- return rc;
-}
-
-void *get_max_d(void *v1, void *v2)
-{
- void *rc;
-
- rc = v2;
- if (*(DCELL *) v1 > *(DCELL *) v2)
- rc = v1;
- return rc;
-}
-
-/* Read one line from a raster map */
-void get_row_c(int fd, void *row, int n)
-{
- Rast_get_c_row(fd, (CELL *) row, n);
-}
-
-void get_row_f(int fd, void *row, int n)
-{
- Rast_get_f_row(fd, (FCELL *) row, n);
-}
-
-void get_row_d(int fd, void *row, int n)
-{
- Rast_get_d_row(fd, (DCELL *) row, n);
-}
-
-/* Write one row to a raster map */
-void put_row_c(int fd, void *row)
-{
- Rast_put_c_row(fd, (CELL *) row);
-}
-
-void put_row_f(int fd, void *row)
-{
- Rast_put_f_row(fd, (FCELL *) row);
-}
-
-void put_row_d(int fd, void *row)
-{
- Rast_put_d_row(fd, (DCELL *) row);
-}
-
-/* Allocate memory for one line of data */
-void *get_buf_c(void)
-{
- return (void *)Rast_allocate_c_buf();
-}
-
-void *get_buf_f(void)
-{
- return (void *)Rast_allocate_f_buf();
-}
-
-void *get_buf_d(void)
-{
- return (void *)Rast_allocate_d_buf();
-}
-
-/* initialize memory to a minimum value */
-void set_min_c(void *v)
-{
- *(CELL *) v = INT_MIN;
-}
-void set_min_f(void *v)
-{
- *(FCELL *) v = FLT_MIN;
-}
-void set_min_d(void *v)
-{
- *(DCELL *) v = DBL_MIN;
-}
-
-/* initialize memory to a maximum value */
-void set_max_c(void *v)
-{
- *(CELL *) v = INT_MAX;
-}
-void set_max_f(void *v)
-{
- *(FCELL *) v = FLT_MAX;
-}
-void set_max_d(void *v)
-{
- *(DCELL *) v = DBL_MAX;
-}
-
-/* get the difference between two values, returned in the first pointer */
-void diff_c(void *v1, void *v2)
-{
- *(CELL *) v1 -= *(CELL *) v2;
-}
-void diff_f(void *v1, void *v2)
-{
- *(FCELL *) v1 -= *(FCELL *) v2;
-}
-void diff_d(void *v1, void *v2)
-{
- *(DCELL *) v1 -= *(DCELL *) v2;
-}
-
-/* get the sum of two values, returned in the first pointer */
-void sum_c(void *v1, void *v2)
-{
- *(CELL *) v1 += *(CELL *) v2;
-}
-void sum_f(void *v1, void *v2)
-{
- *(FCELL *) v1 += *(FCELL *) v2;
-}
-void sum_d(void *v1, void *v2)
-{
- *(DCELL *) v1 += *(DCELL *) v2;
-}
-
-/* get the quotient of two values, returned in the first pointer */
-void quot_c(void *v1, void *v2)
-{
- *(CELL *) v1 /= *(CELL *) v2;
-}
-void quot_f(void *v1, void *v2)
-{
- *(FCELL *) v1 /= *(FCELL *) v2;
-}
-void quot_d(void *v1, void *v2)
-{
- *(DCELL *) v1 /= *(DCELL *) v2;
-}
-
-/* get the product of two values, returned in the first pointer */
-void prod_c(void *v1, void *v2)
-{
- *(CELL *) v1 *= *(CELL *) v2;
-}
-void prod_f(void *v1, void *v2)
-{
- *(FCELL *) v1 *= *(FCELL *) v2;
-}
-void prod_d(void *v1, void *v2)
-{
- *(DCELL *) v1 *= *(DCELL *) v2;
-}
-
-/* probably not a function of general interest */
-/* calculate the slope between two cells, returned as a double */
-double slope_c(void *line1, void *line2, double cnst)
-{
- double rc;
- CELL *pedge;
-
- rc = -HUGE_VAL;
- pedge = (CELL *) line2;
- if (!Rast_is_c_null_value(pedge)) {
- rc = (*(CELL *) line1 - *pedge) / cnst;
- }
- return rc;
-}
-
-double slope_f(void *line1, void *line2, double cnst)
-{
- double rc;
- FCELL *pedge;
-
- rc = -HUGE_VAL;
- pedge = (FCELL *) line2;
- if (!Rast_is_f_null_value(pedge)) {
- rc = (*(FCELL *) line1 - *pedge) / cnst;
- }
- return rc;
-}
-
-double slope_d(void *line1, void *line2, double cnst)
-{
- double rc;
- DCELL *pedge;
-
- rc = -HUGE_VAL;
- pedge = (DCELL *) line2;
- if (!Rast_is_d_null_value(pedge)) {
- rc = (*(DCELL *) line1 - *pedge) / cnst;
- }
- return rc;
-}
-
-/* read a line and update a three-line buffer */
-/* moving forward through a file */
-int advance_band3(int fh, struct band3 *bnd)
-{
- int rc;
- void *hold;
-
- hold = bnd->b[0];
- bnd->b[0] = bnd->b[1];
- bnd->b[1] = bnd->b[2];
- bnd->b[2] = hold;
- if (fh == 0)
- rc = 0;
- else
- rc = read(fh, bnd->b[2], bnd->sz);
- return rc;
-}
-
-/* read a line and update a three-line buffer */
-/* moving backward through a file */
-int retreat_band3(int fh, struct band3 *bnd)
-{
- int rc;
- void *hold;
-
- hold = bnd->b[2];
- bnd->b[2] = bnd->b[1];
- bnd->b[1] = bnd->b[0];
- bnd->b[0] = hold;
- if (fh == 0)
- rc = 0;
- else {
- rc = read(fh, bnd->b[0], bnd->sz);
- lseek(fh, (off_t) - 2 * bnd->sz, SEEK_CUR);
- }
- return rc;
-}
Deleted: grass/trunk/raster/r.path/tinf.h
===================================================================
--- grass/trunk/raster/r.drain/tinf.h 2017-12-04 15:37:44 UTC (rev 71888)
+++ grass/trunk/raster/r.path/tinf.h 2017-12-20 22:06:06 UTC (rev 71962)
@@ -1,102 +0,0 @@
-#include <limits.h>
-#include <math.h>
-#include <unistd.h>
-#include <sys/types.h>
-/* #include <values.h> */
-
-/* to add a new multiple-type function first add three prototypes
- * (one for each type). The functions themselves must be defined
- * elsewhere */
-
-void set_func_pointers(int);
-
-int is_null_c(void *);
-int is_null_f(void *);
-int is_null_d(void *);
-
-void set_null_value_c(void *, int);
-void set_null_value_f(void *, int);
-void set_null_value_d(void *, int);
-
-int bpe_c();
-int bpe_f();
-int bpe_d();
-
-void *get_min_c(void *, void *);
-void *get_min_f(void *, void *);
-void *get_min_d(void *, void *);
-
-void *get_max_c(void *, void *);
-void *get_max_f(void *, void *);
-void *get_max_d(void *, void *);
-
-void get_row_c(int, void *, int);
-void get_row_f(int, void *, int);
-void get_row_d(int, void *, int);
-
-void put_row_c(int, void *);
-void put_row_f(int, void *);
-void put_row_d(int, void *);
-
-void *get_buf_c(void);
-void *get_buf_f(void);
-void *get_buf_d(void);
-
-void set_min_c(void *);
-void set_min_f(void *);
-void set_min_d(void *);
-
-void set_max_c(void *);
-void set_max_f(void *);
-void set_max_d(void *);
-
-void diff_c(void *, void *);
-void diff_f(void *, void *);
-void diff_d(void *, void *);
-
-void sum_c(void *, void *);
-void sum_f(void *, void *);
-void sum_d(void *, void *);
-
-void quot_c(void *, void *);
-void quot_f(void *, void *);
-void quot_d(void *, void *);
-
-void prod_c(void *, void *);
-void prod_f(void *, void *);
-void prod_d(void *, void *);
-
-/* to add a new multitype function, add a pointer for the function and
- * its argument list to the list below */
-
-extern int (*is_null) (void *);
-extern void (*set_null_value) (void *, int);
-extern int (*bpe) ();
-extern void *(*get_max) (void *, void *);
-extern void *(*get_min) (void *, void *);
-extern void (*get_row) (int, void *, int);
-extern void *(*get_buf) ();
-extern void (*put_row) (int, void *);
-extern double (*slope) (void *, void *, double);
-extern void (*set_min) (void *);
-extern void (*set_max) (void *);
-extern void (*diff) (void *, void *);
-extern void (*sum) (void *, void *);
-extern void (*quot) (void *, void *);
-extern void (*prod) (void *, void *);
-
-/* probably not something of general interest */
-
-double slope_c(void *, void *, double);
-double slope_f(void *, void *, double);
-double slope_d(void *, void *, double);
-
-struct band3
-{
- int ns; /* samples per line */
- off_t sz; /* bytes per line */
- char *b[3]; /* pointers to start of each line */
-};
-
-int advance_band3(int, struct band3 *);
-int retreat_band3(int, struct band3 *);
More information about the grass-commit
mailing list