move r.stream.extract to trunk

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@57164 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.stream.extract/Makefile

@@ -0,0 +1,12 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.stream.extract
+
+LIBES     = $(SEGMENTLIB) $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(SEGMENTDEP) $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

raster/r.stream.extract/bseg.c

@@ -0,0 +1,159 @@
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/glocale.h>
+#include "seg.h"
+
+int bseg_open(BSEG *bseg, int srows, int scols, int nsegs_in_memory)
+{
+    char *filename;
+    int errflag;
+    int fd;
+
+    bseg->filename = NULL;
+    bseg->fd = -1;
+    bseg->name = NULL;
+    bseg->mapset = NULL;
+
+    filename = G_tempfile();
+    if (-1 == (fd = creat(filename, 0666))) {
+	G_warning(_("bseg_open(): unable to create segment file"));
+	return -2;
+    }
+    if (0 > (errflag = segment_format(fd, Rast_window_rows(),
+				      Rast_window_cols(), srows, scols,
+				      sizeof(char)))) {
+	close(fd);
+	unlink(filename);
+	if (errflag == -1) {
+	    G_warning(_("bseg_open(): could not write segment file"));
+	    return -1;
+	}
+	else {
+	    G_warning(_("bseg_open(): illegal configuration parameter(s)"));
+	    return -3;
+	}
+    }
+    close(fd);
+    if (-1 == (fd = open(filename, 2))) {
+	unlink(filename);
+	G_warning(_("bseg_open(): unable to re-open segment file"));
+	return -4;
+    }
+    if (0 > (errflag = segment_init(&(bseg->seg), fd, nsegs_in_memory))) {
+	close(fd);
+	unlink(filename);
+	if (errflag == -1) {
+	    G_warning(_("bseg_open(): could not read segment file"));
+	    return -5;
+	}
+	else {
+	    G_warning(_("bseg_open(): out of memory"));
+	    return -6;
+	}
+    }
+    bseg->filename = filename;
+    bseg->fd = fd;
+    return 0;
+}
+
+int bseg_close(BSEG *bseg)
+{
+    segment_release(&(bseg->seg));
+    close(bseg->fd);
+    unlink(bseg->filename);
+    if (bseg->name) {
+	G_free(bseg->name);
+	bseg->name = NULL;
+    }
+    if (bseg->mapset) {
+	G_free(bseg->mapset);
+	bseg->mapset = NULL;
+    }
+    return 0;
+}
+
+int bseg_put(BSEG *bseg, char *value, int row, int col)
+{
+    if (segment_put(&(bseg->seg), value, row, col) < 0) {
+	G_warning(_("bseg_put(): could not write segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int bseg_put_row(BSEG *bseg, char *value, int row)
+{
+    if (segment_put_row(&(bseg->seg), value, row) < 0) {
+	G_warning(_("bseg_put_row(): could not write segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int bseg_get(BSEG *bseg, char *value, int row, int col)
+{
+    if (segment_get(&(bseg->seg), value, row, col) < 0) {
+	G_warning(_("bseg_get(): could not read segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+
+int bseg_read_raster(BSEG *bseg, char *map_name, char *mapset)
+{
+    int row, nrows;
+    int col, ncols;
+    int map_fd;
+    CELL *buffer;
+    char cbuf;
+
+    bseg->name = NULL;
+    bseg->mapset = NULL;
+
+    map_fd = Rast_open_old(map_name, mapset);
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+    buffer = Rast_allocate_c_buf();
+    for (row = 0; row < nrows; row++) {
+	Rast_get_c_row(map_fd, buffer, row);
+	for (col = ncols; col >= 0; col--) {
+	    cbuf = (char) buffer[col];
+	    bseg_put(bseg, &cbuf, row, col);
+	}
+    }
+
+    Rast_close(map_fd);
+    G_free(buffer);
+
+    bseg->name = G_store(map_name);
+    bseg->mapset = G_store(mapset);
+
+    return 0;
+}
+
+int bseg_write_raster(BSEG *bseg, char *map_name)
+{
+    int map_fd;
+    int row, nrows;
+    int col, ncols;
+    CELL *buffer;
+    char value;
+
+    map_fd = Rast_open_c_new(map_name);
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+    buffer = Rast_allocate_c_buf();
+    for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 1);
+	for (col = 0; col < ncols; col++) {
+	    bseg_get(bseg, &value, row, col);
+	    buffer[col] = value;
+	}
+	Rast_put_row(map_fd, buffer, CELL_TYPE);
+    }
+    G_percent(row, nrows, 1);    /* finish it */
+    G_free(buffer);
+    Rast_close(map_fd);
+    return 0;
+}

raster/r.stream.extract/close.c

@@ -0,0 +1,328 @@
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include "local_proto.h"
+
+int close_streamvect(char *stream_vect)
+{
+    int r, c, r_nbr, c_nbr, done;
+    GW_LARGE_INT i;
+    CELL stream_id, stream_nbr;
+    ASP_FLAG af;
+    int next_node;
+    struct sstack
+    {
+	int stream_id;
+	int next_trib;
+    } *nodestack;
+    int top = 0, stack_step = 1000;
+    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    struct Map_info Out;
+    static struct line_pnts *Points;
+    struct line_cats *Cats;
+    dbDriver *driver;
+    dbHandle handle;
+    dbString table_name, dbsql, valstr;
+    struct field_info *Fi;
+    char *cat_col_name = "cat", buf[2000];
+    struct Cell_head window;
+    double north_offset, west_offset, ns_res, ew_res;
+    int next_cat;
+
+    G_message(_("Writing vector map <%s>..."), stream_vect);
+
+    if (0 > Vect_open_new(&Out, stream_vect, 0)) {
+	G_fatal_error(_("Unable to create vector map <%s>"), stream_vect);
+    }
+
+    nodestack = (struct sstack *)G_malloc(stack_step * sizeof(struct sstack));
+
+    Points = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+
+    G_get_set_window(&window);
+    ns_res = window.ns_res;
+    ew_res = window.ew_res;
+    north_offset = window.north - 0.5 * ns_res;
+    west_offset = window.west + 0.5 * ew_res;
+
+    next_cat = n_stream_nodes + 1;
+
+    for (i = 0; i < n_outlets; i++, next_cat++) {
+	G_percent(i, n_outlets, 2);
+	r = outlets[i].r;
+	c = outlets[i].c;
+	cseg_get(&stream, &stream_id, r, c);
+
+	if (!stream_id)
+	    continue;
+
+	Vect_reset_line(Points);
+	Vect_reset_cats(Cats);
+
+	/* outlet */
+	Vect_cat_set(Cats, 1, stream_id);
+	Vect_cat_set(Cats, 2, 2);
+	Vect_append_point(Points, west_offset + c * ew_res,
+			  north_offset - r * ns_res, 0);
+	Vect_write_line(&Out, GV_POINT, Points, Cats);
+
+	/* add root node to stack */
+	G_debug(3, "add root node");
+	top = 0;
+	nodestack[top].stream_id = stream_id;
+	nodestack[top].next_trib = 0;
+
+	/* depth first post order traversal */
+	G_debug(3, "traverse");
+	while (top >= 0) {
+
+	    done = 1;
+	    stream_id = nodestack[top].stream_id;
+	    G_debug(3, "stream_id %d", stream_id);
+	    if (nodestack[top].next_trib < stream_node[stream_id].n_trib) {
+		/* add to stack */
+		next_node =
+		    stream_node[stream_id].trib[nodestack[top].next_trib];
+		G_debug(3, "add to stack: next %d, trib %d, n trib %d",
+			next_node, nodestack[top].next_trib,
+			stream_node[stream_id].n_trib);
+		nodestack[top].next_trib++;
+		top++;
+		if (top >= stack_step) {
+		    /* need more space */
+		    stack_step += 1000;
+		    nodestack =
+			(struct sstack *)G_realloc(nodestack,
+						   stack_step *
+						   sizeof(struct sstack));
+		}
+		nodestack[top].next_trib = 0;
+		nodestack[top].stream_id = next_node;
+		done = 0;
+		G_debug(3, "go further down");
+	    }
+	    if (done) {
+		G_debug(3, "write stream segment");
+
+		Vect_reset_line(Points);
+		Vect_reset_cats(Cats);
+
+		r_nbr = stream_node[stream_id].r;
+		c_nbr = stream_node[stream_id].c;
+
+		cseg_get(&stream, &stream_nbr, r_nbr, c_nbr);
+		if (stream_nbr <= 0)
+		    G_fatal_error("stream id %d not set, top is %d, parent is %d", stream_id, top, nodestack[top - 1].stream_id);
+
+		Vect_cat_set(Cats, 1, stream_id);
+		if (stream_node[stream_id].n_trib == 0)
+		    Vect_cat_set(Cats, 2, 0);
+		else
+		    Vect_cat_set(Cats, 2, 1);
+
+		Vect_append_point(Points, west_offset + c_nbr * ew_res,
+				  north_offset - r_nbr * ns_res, 0);
+
+		Vect_write_line(&Out, GV_POINT, Points, Cats);
+
+		seg_get(&aspflag, (char *)&af, r_nbr, c_nbr);
+		while (af.asp > 0) {
+		    r_nbr = r_nbr + asp_r[(int)af.asp];
+		    c_nbr = c_nbr + asp_c[(int)af.asp];
+		    
+		    cseg_get(&stream, &stream_nbr, r_nbr, c_nbr);
+		    if (stream_nbr <= 0)
+			G_fatal_error("stream id not set while tracing");
+
+		    Vect_append_point(Points, west_offset + c_nbr * ew_res,
+				      north_offset - r_nbr * ns_res, 0);
+		    if (stream_nbr != stream_id) {
+			/* first point of parent stream */
+			break;
+		    }
+		    seg_get(&aspflag, (char *)&af, r_nbr, c_nbr);
+		}
+
+		Vect_write_line(&Out, GV_LINE, Points, Cats);
+
+		top--;
+	    }
+	}
+    }
+    G_percent(n_outlets, n_outlets, 1);	/* finish it */
+
+    G_message(_("Write vector attribute table"));
+
+    /* Prepeare strings for use in db_* calls */
+    db_init_string(&dbsql);
+    db_init_string(&valstr);
+    db_init_string(&table_name);
+    db_init_handle(&handle);
+
+    /* Preparing database for use */
+    /* Create database for new vector map */
+    Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
+    driver = db_start_driver_open_database(Fi->driver,
+					   Vect_subst_var(Fi->database,
+							          &Out));
+    if (driver == NULL) {
+	G_fatal_error(_("Unable to start driver <%s>"), Fi->driver);
+    }
+
+    G_debug(1, "table: %s", Fi->table);
+    G_debug(1, "driver: %s", Fi->driver);
+    G_debug(1, "database: %s", Fi->database);
+
+    sprintf(buf,
+	    "create table %s (%s integer, stream_type varchar(20), type_code integer)",
+	    Fi->table, cat_col_name);
+    db_set_string(&dbsql, buf);
+
+    if (db_execute_immediate(driver, &dbsql) != DB_OK) {
+	db_close_database(driver);
+	db_shutdown_driver(driver);
+	G_fatal_error(_("Cannot create table: %s"), db_get_string(&dbsql));
+    }
+
+    if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
+	G_warning(_("Cannot create index"));
+
+    if (db_grant_on_table(driver, Fi->table, DB_PRIV_SELECT,
+			  DB_GROUP | DB_PUBLIC) != DB_OK)
+	G_fatal_error(_("Cannot grant privileges on table %s"), Fi->table);
+
+    db_begin_transaction(driver);
+
+    /* stream nodes */
+    for (i = 1; i <= n_stream_nodes; i++) {
+
+	sprintf(buf, "insert into %s values ( %lld, \'%s\', %d )",
+		Fi->table, i,
+		(stream_node[i].n_trib > 0 ? "intermediate" : "start"),
+		(stream_node[i].n_trib > 0));
+
+	db_set_string(&dbsql, buf);
+
+	if (db_execute_immediate(driver, &dbsql) != DB_OK) {
+	    db_close_database(driver);
+	    db_shutdown_driver(driver);
+	    G_fatal_error(_("Cannot insert new row: %s"),
+			  db_get_string(&dbsql));
+	}
+    }
+
+    db_commit_transaction(driver);
+    db_close_database_shutdown_driver(driver);
+
+    Vect_map_add_dblink(&Out, 1, NULL, Fi->table,
+			cat_col_name, Fi->database, Fi->driver);
+
+    G_debug(1, "close vector");
+
+    Vect_hist_command(&Out);
+    Vect_build(&Out);
+    Vect_close(&Out);
+
+    G_free(nodestack);
+
+    return 1;
+}
+
+
+int close_maps(char *stream_rast, char *stream_vect, char *dir_rast)
+{
+    int stream_fd, dir_fd, r, c, i;
+    CELL *cell_buf1, *cell_buf2;
+    struct History history;
+    CELL stream_id;
+    ASP_FLAG af;
+
+    /* cheating... */
+    stream_fd = dir_fd = -1;
+    cell_buf1 = cell_buf2 = NULL;
+
+    G_message(_("Writing raster %s"),
+              (stream_rast != NULL) + (dir_rast != NULL) > 1 ? "maps" : "map");
+
+    /* write requested output rasters */
+    if (stream_rast) {
+	stream_fd = Rast_open_new(stream_rast, CELL_TYPE);
+	cell_buf1 = Rast_allocate_c_buf();
+    }
+    if (dir_rast) {
+	dir_fd = Rast_open_new(dir_rast, CELL_TYPE);
+	cell_buf2 = Rast_allocate_c_buf();
+    }
+
+    for (r = 0; r < nrows; r++) {
+	G_percent(r, nrows, 2);
+	if (stream_rast)
+	    Rast_set_c_null_value(cell_buf1, ncols);	/* reset row to all NULL */
+	if (dir_rast)
+	    Rast_set_c_null_value(cell_buf2, ncols);	/* reset row to all NULL */
+
+	for (c = 0; c < ncols; c++) {
+	    if (stream_rast) {
+		cseg_get(&stream, &stream_id, r, c);
+		if (stream_id)
+		    cell_buf1[c] = stream_id;
+	    }
+	    if (dir_rast) {
+		seg_get(&aspflag, (char *)&af, r, c);
+		if (!FLAG_GET(af.flag, NULLFLAG)) {
+		    cell_buf2[c] = af.asp;
+		}
+	    }
+	    
+	}
+	if (stream_rast)
+	    Rast_put_row(stream_fd, cell_buf1, CELL_TYPE);
+	if (dir_rast)
+	    Rast_put_row(dir_fd, cell_buf2, CELL_TYPE);
+    }
+    G_percent(nrows, nrows, 2);	/* finish it */
+
+    if (stream_rast) {
+	Rast_close(stream_fd);
+	G_free(cell_buf1);
+	Rast_short_history(stream_rast, "raster", &history);
+	Rast_command_history(&history);
+	Rast_write_history(stream_rast, &history);
+    }
+    if (dir_rast) {
+	struct Colors colors;
+
+	Rast_close(dir_fd);
+	G_free(cell_buf2);
+
+	Rast_short_history(dir_rast, "raster", &history);
+	Rast_command_history(&history);
+	Rast_write_history(dir_rast, &history);
+
+	Rast_init_colors(&colors);
+	Rast_make_aspect_colors(&colors, -8, 8);
+	Rast_write_colors(dir_rast, G_mapset(), &colors);
+    }
+
+    /* close stream vector */
+    if (stream_vect) {
+	if (close_streamvect(stream_vect) < 0)
+	    G_fatal_error(_("Unable to write vector map <%s>"), stream_vect);
+    }
+
+    /* rearranging desk chairs on the Titanic... */
+    G_free(outlets);
+
+    /* free stream nodes */
+    for (i = 1; i <= n_stream_nodes; i++) {
+	if (stream_node[i].n_alloc > 0) {
+	    G_free(stream_node[i].trib);
+	}
+    }
+    G_free(stream_node);
+
+    return 1;
+}

raster/r.stream.extract/cseg.c

@@ -0,0 +1,154 @@
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/glocale.h>
+#include "seg.h"
+
+int cseg_open(CSEG *cseg, int srows, int scols, int nsegs_in_memory)
+{
+    char *filename;
+    int errflag;
+    int fd;
+
+    cseg->filename = NULL;
+    cseg->fd = -1;
+    cseg->name = NULL;
+    cseg->mapset = NULL;
+
+    filename = G_tempfile();
+    if (-1 == (fd = creat(filename, 0666))) {
+	G_warning(_("cseg_open(): unable to create segment file"));
+	return -2;
+    }
+    if (0 >
+	(errflag =
+	 segment_format(fd, Rast_window_rows(), Rast_window_cols(), srows, scols,
+			sizeof(CELL)))) {
+	close(fd);
+	unlink(filename);
+	if (errflag == -1) {
+	    G_warning(_("cseg_open(): could not write segment file"));
+	    return -1;
+	}
+	else {
+	    G_warning(_("cseg_open(): illegal configuration parameter(s)"));
+	    return -3;
+	}
+    }
+    close(fd);
+    if (-1 == (fd = open(filename, 2))) {
+	unlink(filename);
+	G_warning(_("cseg_open(): unable to re-open segment file"));
+	return -4;
+    }
+    if (0 > (errflag = segment_init(&(cseg->seg), fd, nsegs_in_memory))) {
+	close(fd);
+	unlink(filename);
+	if (errflag == -1) {
+	    G_warning(_("cseg_open(): could not read segment file"));
+	    return -5;
+	}
+	else {
+	    G_warning(_("cseg_open(): out of memory"));
+	    return -6;
+	}
+    }
+    cseg->filename = filename;
+    cseg->fd = fd;
+    return 0;
+}
+
+int cseg_close(CSEG *cseg)
+{
+    segment_release(&(cseg->seg));
+    close(cseg->fd);
+    unlink(cseg->filename);
+    if (cseg->name) {
+	G_free(cseg->name);
+	cseg->name = NULL;
+    }
+    if (cseg->mapset) {
+	G_free(cseg->mapset);
+	cseg->mapset = NULL;
+    }
+    return 0;
+}
+
+int cseg_put(CSEG *cseg, CELL *value, int row, int col)
+{
+    if (segment_put(&(cseg->seg), value, row, col) < 0) {
+	G_warning(_("cseg_put(): could not write segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int cseg_put_row(CSEG *cseg, CELL *value, int row)
+{
+    if (segment_put_row(&(cseg->seg), value, row) < 0) {
+	G_warning(_("cseg_put_row(): could not write segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int cseg_get(CSEG *cseg, CELL *value, int row, int col)
+{
+    if (segment_get(&(cseg->seg), value, row, col) < 0) {
+	G_warning(_("cseg_get(): could not read segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int cseg_read_raster(CSEG *cseg, char *map_name, char *mapset)
+{
+    int row, nrows;
+    int map_fd;
+    CELL *buffer;
+
+    cseg->name = NULL;
+    cseg->mapset = NULL;
+
+    map_fd = Rast_open_old(map_name, mapset);
+    nrows = Rast_window_rows();
+    buffer = Rast_allocate_c_buf();
+    for (row = 0; row < nrows; row++) {
+	Rast_get_c_row(map_fd, buffer, row);
+	if (segment_put_row(&(cseg->seg), buffer, row) < 0) {
+	    G_free(buffer);
+	    Rast_close(map_fd);
+	    G_warning(_("cseg_read_cell(): unable to segment put row for <%s> in <%s>"),
+		    map_name, mapset);
+	    return (-1);
+	}
+    }
+
+    Rast_close(map_fd);
+    G_free(buffer);
+
+    cseg->name = G_store(map_name);
+    cseg->mapset = G_store(mapset);
+
+    return 0;
+}
+
+int cseg_write_raster(CSEG *cseg, char *map_name)
+{
+    int map_fd;
+    int row, nrows;
+    CELL *buffer;
+
+    map_fd = Rast_open_c_new(map_name);
+    nrows = Rast_window_rows();
+    buffer = Rast_allocate_c_buf();
+    segment_flush(&(cseg->seg));
+    for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 1);
+	segment_get_row(&(cseg->seg), buffer, row);
+	Rast_put_row(map_fd, buffer, CELL_TYPE);
+    }
+    G_percent(row, nrows, 1);    /* finish it */
+    G_free(buffer);
+    Rast_close(map_fd);
+    return 0;
+}

raster/r.stream.extract/del_streams.c

@@ -0,0 +1,194 @@
+#include <stdlib.h>
+#include <math.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+/* get stream segment length */
+int seg_length(CELL stream_id, CELL *next_stream_id)
+{
+    int r, c, r_nbr, c_nbr;
+    int slength = 1;
+    CELL curr_stream;
+    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    ASP_FLAG af;
+
+    r = stream_node[stream_id].r;
+    c = stream_node[stream_id].c;
+    if (next_stream_id)
+	*next_stream_id = stream_id;
+
+    /* get next downstream point */
+    seg_get(&aspflag, (char *)&af, r, c);
+    while (af.asp > 0) {
+	r_nbr = r + asp_r[(int)af.asp];
+	c_nbr = c + asp_c[(int)af.asp];
+
+	/* user-defined depression */
+	if (r_nbr == r && c_nbr == c)
+	    break;
+	/* outside region */
+	if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
+	    break;
+	/* next stream */
+	cseg_get(&stream, &curr_stream, r_nbr, c_nbr);
+	if (next_stream_id)
+	    *next_stream_id = curr_stream;
+	if (curr_stream != stream_id)
+	    break;
+	slength++;
+	r = r_nbr;
+	c = c_nbr;
+	seg_get(&aspflag, (char *)&af, r, c);
+    }
+
+    return slength;
+}
+
+/* change downstream id: update or remove */
+int update_stream_id(CELL stream_id, CELL new_stream_id)
+{
+    int i, r, c, r_nbr, c_nbr;
+    CELL new_stream = new_stream_id;
+    CELL curr_stream;
+    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    ASP_FLAG af;
+
+    r = stream_node[stream_id].r;
+    c = stream_node[stream_id].c;
+    cseg_get(&stream, &curr_stream, r, c);
+    if (curr_stream != stream_id)
+	G_fatal_error("update downstream id: curr_stream != stream_id");
+    cseg_put(&stream, &new_stream, r, c);
+    curr_stream = stream_id;
+
+    /* get next downstream point */
+    seg_get(&aspflag, (char *)&af, r, c);
+    while (af.asp > 0) {
+	r_nbr = r + asp_r[(int)af.asp];
+	c_nbr = c + asp_c[(int)af.asp];
+
+	/* user-defined depression */
+	if (r_nbr == r && c_nbr == c)
+	    break;
+	/* outside region */
+	if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
+	    break;
+	/* next stream */
+	cseg_get(&stream, &curr_stream, r_nbr, c_nbr);
+	if (curr_stream != stream_id)
+	    break;
+	r = r_nbr;
+	c = c_nbr;
+	cseg_put(&stream, &new_stream, r, c);
+	seg_get(&aspflag, (char *)&af, r, c);
+    }
+    
+    if (curr_stream <= 0)
+	return curr_stream;
+
+    /* update tributaries */
+    if (curr_stream != stream_id) {
+	int last_i = -1;
+	
+	for (i = 0; i < stream_node[curr_stream].n_trib; i++) {
+	    if (stream_node[curr_stream].trib[i] == stream_id) {
+		if (new_stream_id) {
+		    stream_node[curr_stream].trib[i] = new_stream_id;
+		}
+		else {
+		    stream_node[curr_stream].n_trib--;
+		    stream_node[curr_stream].trib[i] = stream_node[curr_stream].trib[stream_node[curr_stream].n_trib];
+		    stream_node[curr_stream].trib[stream_node[curr_stream].n_trib] = 0;
+		}
+		last_i = i;
+		break;
+	    }
+	}
+	for (i = 0; i < stream_node[curr_stream].n_trib; i++) {
+	    if (stream_node[curr_stream].trib[i] == stream_id) {
+		G_warning("last_i %d, i %d", last_i, i);
+		G_warning("failed updating stream %d at node %d", stream_id, curr_stream);
+	    }
+	}
+    }
+
+    return curr_stream;
+}
+
+/* delete stream segments shorter than threshold */
+int del_streams(int min_length)
+{
+    int i;
+    int n_deleted = 0;
+    CELL curr_stream, stream_id;
+    int other_trib, tmp_trib;
+    int slength;
+
+    G_message(_("Deleting stream segments shorter than %d cells..."), min_length);
+
+    /* TODO: proceed from stream heads to outlets
+     *       -> use depth first post order traversal */
+
+    /* go through all nodes */
+    for (i = 1; i <= n_stream_nodes; i++) {
+	G_percent(i, n_stream_nodes, 2);
+
+	/* not a stream head */
+	if (stream_node[i].n_trib > 0)
+	    continue;
+
+	/* already deleted */
+	cseg_get(&stream, &curr_stream, stream_node[i].r, stream_node[i].c);
+	if (curr_stream == 0)
+	    continue;
+
+	/* get length counted as n cells */
+	if ((slength = seg_length(i, &curr_stream)) >= min_length)
+	    continue;
+
+	stream_id = i;
+	
+	/* check n sibling tributaries */
+	if (curr_stream != stream_id) {
+	    /* only one sibling tributary */
+	    if (stream_node[curr_stream].n_trib == 2) {
+		if (stream_node[curr_stream].trib[0] != stream_id)
+		    other_trib = stream_node[curr_stream].trib[0];
+		else
+		    other_trib = stream_node[curr_stream].trib[1];
+
+		/* other trib is also stream head */
+		if (stream_node[other_trib].n_trib == 0) {
+		    /* use shorter one */
+		    if (seg_length(other_trib, NULL) < slength) {
+			tmp_trib = stream_id;
+			stream_id = other_trib;
+			other_trib = tmp_trib;
+		    }
+		}
+		update_stream_id(stream_id, 0);
+		n_deleted++;
+		
+		/* update downstream IDs */
+		update_stream_id(curr_stream, other_trib);
+	    }
+	    /* more than one other tributary */
+	    else {
+		update_stream_id(stream_id, 0);
+		n_deleted++;
+	    }
+	}
+	/* stream head is also outlet */
+	else {
+	    update_stream_id(stream_id, 0);
+	    n_deleted++;
+	}
+    }
+
+    G_verbose_message(_("%d stream segments deleted"), n_deleted);
+
+    return n_deleted;
+}

raster/r.stream.extract/do_astar.c

@@ -0,0 +1,288 @@
+#include <stdlib.h>
+#include <math.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+#define GET_PARENT(c) ((((c) - 2) >> 3) + 1)
+#define GET_CHILD(p) (((p) << 3) - 6)
+
+HEAP_PNT heap_drop(void);
+static double get_slope(CELL, CELL, double);
+
+int do_astar(void)
+{
+    int r, c, r_nbr, c_nbr, ct_dir;
+    GW_LARGE_INT first_cum, count;
+    int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+    int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+    CELL ele_val, ele_up, ele_nbr[8];
+    WAT_ALT wa;
+    ASP_FLAG af;
+    char is_in_list, is_worked;
+    HEAP_PNT heap_p;
+    /* sides
+     * |7|1|4|
+     * |2| |3|
+     * |5|0|6|
+     */
+    int nbr_ew[8] = { 0, 1, 2, 3, 1, 0, 0, 1 };
+    int nbr_ns[8] = { 0, 1, 2, 3, 3, 2, 3, 2 };
+    double dx, dy, dist_to_nbr[8], ew_res, ns_res;
+    double slope[8];
+    struct Cell_head window;
+    int skip_diag;
+
+    count = 0;
+
+    first_cum = n_points;
+
+    G_message(_("A* Search..."));
+
+    Rast_get_window(&window);
+
+    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+	/* get r, c (r_nbr, c_nbr) for neighbours */
+	r_nbr = nextdr[ct_dir];
+	c_nbr = nextdc[ct_dir];
+	/* account for rare cases when ns_res != ew_res */
+	dy = abs(r_nbr) * window.ns_res;
+	dx = abs(c_nbr) * window.ew_res;
+	if (ct_dir < 4)
+	    dist_to_nbr[ct_dir] = dx + dy;
+	else
+	    dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+    }
+    ew_res = window.ew_res;
+    ns_res = window.ns_res;
+    
+    while (heap_size > 0) {
+	G_percent(count++, n_points, 1);
+	if (count > n_points)
+	    G_fatal_error(_("BUG in A* Search: %lld surplus points"),
+	                  heap_size);
+
+	if (heap_size > n_points)
+	    G_fatal_error
+		(_("BUG in A* Search: too many points in heap %lld, should be %lld"),
+		 heap_size, n_points);
+
+	heap_p = heap_drop();
+
+	r = heap_p.pnt.r;
+	c = heap_p.pnt.c;
+
+	ele_val = heap_p.ele;
+
+	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+	    /* get r, c (r_nbr, c_nbr) for neighbours */
+	    r_nbr = r + nextdr[ct_dir];
+	    c_nbr = c + nextdc[ct_dir];
+	    slope[ct_dir] = ele_nbr[ct_dir] = 0;
+	    skip_diag = 0;
+
+	    /* check that neighbour is within region */
+	    if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
+		continue;
+
+	    seg_get(&aspflag, (char *)&af, r_nbr, c_nbr);
+	    is_in_list = FLAG_GET(af.flag, INLISTFLAG);
+	    is_worked = FLAG_GET(af.flag, WORKEDFLAG);
+	    if (!is_worked) {
+		seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
+		ele_nbr[ct_dir] = wa.ele;
+		slope[ct_dir] = get_slope(ele_val, ele_nbr[ct_dir],
+			                  dist_to_nbr[ct_dir]);
+	    }
+	    /* avoid diagonal flow direction bias */
+	    if (!is_in_list) {
+		if (ct_dir > 3 && slope[ct_dir] > 0) {
+		    if (slope[nbr_ew[ct_dir]] > 0) {
+			/* slope to ew nbr > slope to center */
+			if (slope[ct_dir] <
+			    get_slope(ele_nbr[nbr_ew[ct_dir]],
+				       ele_nbr[ct_dir], ew_res))
+			    skip_diag = 1;
+		    }
+		    if (!skip_diag && slope[nbr_ns[ct_dir]] > 0) {
+			/* slope to ns nbr > slope to center */
+			if (slope[ct_dir] <
+			    get_slope(ele_nbr[nbr_ns[ct_dir]],
+				       ele_nbr[ct_dir], ns_res))
+			    skip_diag = 1;
+		    }
+		}
+	    }
+
+	    if (!skip_diag) {
+		if (is_in_list == 0) {
+		    ele_up = ele_nbr[ct_dir];
+		    af.asp = drain[r_nbr - r + 1][c_nbr - c + 1];
+		    heap_add(r_nbr, c_nbr, ele_up);
+		    FLAG_SET(af.flag, INLISTFLAG);
+		    seg_put(&aspflag, (char *)&af, r_nbr, c_nbr);
+		}
+		else if (is_in_list && is_worked == 0) {
+		    if (FLAG_GET(af.flag, EDGEFLAG)) {
+			/* neighbour is edge in list, not yet worked */
+			if (af.asp < 0) {
+			    /* adjust flow direction for edge cell */
+			    af.asp = drain[r_nbr - r + 1][c_nbr - c + 1];
+			    seg_put(&aspflag, (char *)&af, r_nbr, c_nbr);
+			}
+		    }
+		    else if (FLAG_GET(af.flag, DEPRFLAG)) {
+			G_debug(3, "real depression");
+			/* neighbour is inside real depression, not yet worked */
+			if (af.asp == 0 && ele_val <= ele_nbr[ct_dir]) {
+			    af.asp = drain[r_nbr - r + 1][c_nbr - c + 1];
+			    FLAG_UNSET(af.flag, DEPRFLAG);
+			    seg_put(&aspflag, (char *)&af, r_nbr, c_nbr);
+			}
+		    }
+		}
+	    }
+	}    /* end neighbours */
+	/* add astar points to sorted list for flow accumulation and stream extraction */
+	first_cum--;
+	seg_put(&astar_pts, (char *)&heap_p.pnt, 0, first_cum);
+	seg_get(&aspflag, (char *)&af, r, c);
+	FLAG_SET(af.flag, WORKEDFLAG);
+	seg_put(&aspflag, (char *)&af, r, c);
+    }    /* end A* search */
+
+    G_percent(n_points, n_points, 1);	/* finish it */
+
+    return 1;
+}
+
+/*
+ * compare function for heap
+ * returns 1 if point1 < point2 else 0
+ */
+static int heap_cmp(HEAP_PNT *a, HEAP_PNT *b)
+{
+    if (a->ele < b->ele)
+	return 1;
+    else if (a->ele == b->ele) {
+	if (a->added < b->added)
+	    return 1;
+    }
+    return 0;
+}
+
+static int sift_up(GW_LARGE_INT start, HEAP_PNT child_p)
+{
+    GW_LARGE_INT parent, child;
+    HEAP_PNT heap_p;
+
+    child = start;
+
+    while (child > 1) {
+	parent = GET_PARENT(child);
+	seg_get(&search_heap, (char *)&heap_p, 0, parent);
+
+	/* push parent point down if child is smaller */
+	if (heap_cmp(&child_p, &heap_p)) {
+	    seg_put(&search_heap, (char *)&heap_p, 0, child);
+	    child = parent;
+	}
+	else
+	    /* no more sifting up, found slot for child */
+	    break;
+    }
+
+    /* add child to heap */
+    seg_put(&search_heap, (char *)&child_p, 0, child);
+
+    return 0;
+}
+
+/*
+ * add item to heap
+ * returns heap_size
+ */
+GW_LARGE_INT heap_add(int r, int c, CELL ele)
+{
+    HEAP_PNT heap_p;
+    
+    /* add point to next free position */
+
+    heap_size++;
+
+    heap_p.added = nxt_avail_pt;
+    heap_p.ele = ele;
+    heap_p.pnt.r = r;
+    heap_p.pnt.c = c;
+
+    nxt_avail_pt++;
+
+    /* sift up: move new point towards top of heap */
+    sift_up(heap_size, heap_p);
+
+    return heap_size;
+}
+
+/*
+ * drop item from heap
+ * returns heap size
+ */
+HEAP_PNT heap_drop(void)
+{
+    GW_LARGE_INT child, childr, parent;
+    GW_LARGE_INT i;
+    HEAP_PNT child_p, childr_p, last_p, root_p;
+
+    seg_get(&search_heap, (char *)&last_p, 0, heap_size);
+    seg_get(&search_heap, (char *)&root_p, 0, 1);
+
+    if (heap_size == 1) {
+	heap_size = 0;
+	return root_p;
+    }
+
+    parent = 1;
+    while ((child = GET_CHILD(parent)) < heap_size) {
+
+	seg_get(&search_heap, (char *)&child_p, 0, child);
+
+	if (child < heap_size) {
+	    childr = child + 1;
+	    i = child + 8;
+	    while (childr < heap_size && childr < i) {
+		seg_get(&search_heap, (char *)&childr_p, 0, childr);
+		if (heap_cmp(&childr_p, &child_p)) {
+		    child = childr;
+		    child_p = childr_p;
+		}
+		childr++;
+	    }
+	}
+
+	if (heap_cmp(&last_p, &child_p)) {
+	    break;
+	}
+
+	/* move hole down */
+	seg_put(&search_heap, (char *)&child_p, 0, parent);
+	parent = child;
+    }
+
+    /* fill hole */
+    if (parent < heap_size) {
+	seg_put(&search_heap, (char *)&last_p, 0, parent);
+    }
+
+    /* the actual drop */
+    heap_size--;
+
+    return root_p;
+}
+
+static double get_slope(CELL ele, CELL up_ele, double dist)
+{
+    if (ele >= up_ele)
+	return 0.0;
+    else
+	return (double)(up_ele - ele) / dist;
+}

raster/r.stream.extract/dseg.c

@@ -0,0 +1,154 @@
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/glocale.h>
+#include "seg.h"
+
+int dseg_open(DSEG *dseg, int srows, int scols, int nsegs_in_memory)
+{
+    char *filename;
+    int errflag;
+    int fd;
+
+    dseg->filename = NULL;
+    dseg->fd = -1;
+    dseg->name = NULL;
+    dseg->mapset = NULL;
+
+    filename = G_tempfile();
+    if (-1 == (fd = creat(filename, 0666))) {
+	G_warning(_("dseg_open(): unable to create segment file"));
+	return -2;
+    }
+    if (0 >
+	(errflag =
+	 segment_format(fd, Rast_window_rows(), Rast_window_cols(), srows, scols,
+			sizeof(DCELL)))) {
+	close(fd);
+	unlink(filename);
+	if (errflag == -1) {
+	    G_warning(_("dseg_open(): could not write segment file"));
+	    return -1;
+	}
+	else {
+	    G_warning(_("dseg_open(): illegal configuration parameter(s)"));
+	    return -3;
+	}
+    }
+    close(fd);
+    if (-1 == (fd = open(filename, 2))) {
+	unlink(filename);
+	G_warning(_("dseg_open(): unable to re-open segment file"));
+	return -4;
+    }
+    if (0 > (errflag = segment_init(&(dseg->seg), fd, nsegs_in_memory))) {
+	close(fd);
+	unlink(filename);
+	if (errflag == -1) {
+	    G_warning(_("dseg_open(): could not read segment file"));
+	    return -5;
+	}
+	else {
+	    G_warning(_("dseg_open(): out of memory"));
+	    return -6;
+	}
+    }
+    dseg->filename = filename;
+    dseg->fd = fd;
+    return 0;
+}
+
+int dseg_close(DSEG *dseg)
+{
+    segment_release(&(dseg->seg));
+    close(dseg->fd);
+    unlink(dseg->filename);
+    if (dseg->name) {
+	G_free(dseg->name);
+	dseg->name = NULL;
+    }
+    if (dseg->mapset) {
+	G_free(dseg->mapset);
+	dseg->mapset = NULL;
+    }
+    return 0;
+}
+
+int dseg_put(DSEG *dseg, DCELL *value, int row, int col)
+{
+    if (segment_put(&(dseg->seg), (DCELL *) value, row, col) < 0) {
+	G_warning(_("dseg_put(): could not write segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int dseg_put_row(DSEG *dseg, DCELL *value, int row)
+{
+    if (segment_put_row(&(dseg->seg), (DCELL *) value, row) < 0) {
+	G_warning(_("dseg_put(): could not write segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int dseg_get(DSEG *dseg, DCELL *value, int row, int col)
+{
+    if (segment_get(&(dseg->seg), (DCELL *) value, row, col) < 0) {
+	G_warning(_("dseg_get(): could not read segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int dseg_read_raster(DSEG *dseg, char *map_name, char *mapset)
+{
+    int row, nrows;
+    int map_fd;
+    DCELL *dbuffer;
+
+    dseg->name = NULL;
+    dseg->mapset = NULL;
+
+    map_fd = Rast_open_old(map_name, mapset);
+    nrows = Rast_window_rows();
+    dbuffer = Rast_allocate_d_buf();
+    for (row = 0; row < nrows; row++) {
+	Rast_get_d_row(map_fd, dbuffer, row);
+	if (segment_put_row(&(dseg->seg), (DCELL *) dbuffer, row) < 0) {
+	    G_free(dbuffer);
+	    Rast_close(map_fd);
+	    G_warning(_("dseg_read_raster(): unable to segment put row for <%s> in <%s>"),
+		    map_name, mapset);
+	    return (-1);
+	}
+    }
+
+    Rast_close(map_fd);
+    G_free(dbuffer);
+
+    dseg->name = G_store(map_name);
+    dseg->mapset = G_store(mapset);
+
+    return 0;
+}
+
+int dseg_write_cellfile(DSEG *dseg, char *map_name)
+{
+    int map_fd;
+    int row, nrows;
+    DCELL *dbuffer;
+
+    map_fd = Rast_open_new(map_name, DCELL_TYPE);
+    nrows = Rast_window_rows();
+    dbuffer = Rast_allocate_d_buf();
+    segment_flush(&(dseg->seg));
+    for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 1);
+	segment_get_row(&(dseg->seg), (DCELL *) dbuffer, row);
+	Rast_put_row(map_fd, dbuffer, DCELL_TYPE);
+    }
+    G_percent(row, nrows, 1);    /* finish it */
+    G_free(dbuffer);
+    Rast_close(map_fd);
+    return 0;
+}

raster/r.stream.extract/flag.h

@@ -0,0 +1,36 @@
+#ifndef __FLAG_H__
+#define __FLAG_H__
+
+/* a set of routines that allow the programmer to "flag" cells in a
+ * raster map. A flag is of type unsigned char, i.e. 8 bits can be set. 
+ *
+ * int flag_set(flag, bitno)
+ *     sets the flag at position bitno to one.
+ *
+ * int flag_unset(flag, bitno)
+ *     sets the flag at position bitno to zero.
+ *
+ * int flag_get(flag, bitno)
+ *     checks if the flag is set at postion bitno.
+ *
+ * Examples:
+ * set flag at position 0: FLAG_SET(flag, 0)
+ * unset (clear) flag at position 7: FLAG_UNSET(flag, 7)
+ * check flag at position 5: is_set_at_5 = FLAG_GET(flag, 5)
+ */
+
+/* flag positions */
+#define NULLFLAG         0      /* elevation is NULL */
+#define EDGEFLAG         1      /* edge cell */
+#define INLISTFLAG       2      /* in open A* list */
+#define WORKEDFLAG       3      /* in closed A* list/ accumulation done */
+#define STREAMFLAG       4      /* stream */
+#define DEPRFLAG         5      /* real depression */
+#define WORKED2FLAG      6      /* extraction done */
+/* last bit is unused */
+
+#define FLAG_SET(flag,bitno) ((flag) |= (1 << (bitno)))
+#define FLAG_UNSET(flag,bitno) ((flag) &= ~(1 << (bitno)))
+#define FLAG_GET(flag,bitno) ((flag) & (1 << (bitno)))
+
+#endif /* __FLAG_H__ */

raster/r.stream.extract/init_search.c

@@ -0,0 +1,122 @@
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+int init_search(int depr_fd)
+{
+    int r, c, r_nbr, c_nbr, ct_dir;
+    CELL *depr_buf, ele_value;
+    int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+    int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+    char asp_value, is_null;
+    WAT_ALT wa;
+    ASP_FLAG af, af_nbr;
+    GW_LARGE_INT n_depr_cells = 0;
+
+    nxt_avail_pt = heap_size = 0;
+
+    /* load edge cells and real depressions to A* heap */
+    if (depr_fd >= 0)
+	depr_buf = Rast_allocate_buf(CELL_TYPE);
+    else
+	depr_buf = NULL;
+
+    G_message(_("Initializing A* Search..."));
+    for (r = 0; r < nrows; r++) {
+	G_percent(r, nrows, 2);
+
+	if (depr_fd >= 0) {
+	    Rast_get_row(depr_fd, depr_buf, r, CELL_TYPE);
+	}
+
+	for (c = 0; c < ncols; c++) {
+
+	    seg_get(&aspflag, (char *)&af, r, c);
+	    is_null = FLAG_GET(af.flag, NULLFLAG);
+
+	    if (is_null)
+		continue;
+
+	    asp_value = 0;
+	    if (r == 0 || r == nrows - 1 || c == 0 || c == ncols - 1) {
+
+		if (r == 0 && c == 0)
+		    asp_value = -7;
+		else if (r == 0 && c == ncols - 1)
+		    asp_value = -5;
+		else if (r == nrows - 1 && c == 0)
+		    asp_value = -1;
+		else if (r == nrows - 1 && c == ncols - 1)
+		    asp_value = -3;
+		else if (r == 0)
+		    asp_value = -2;
+		else if (c == 0)
+		    asp_value = -4;
+		else if (r == nrows - 1)
+		    asp_value = -6;
+		else if (c == ncols - 1)
+		    asp_value = -8;
+
+		seg_get(&watalt, (char *)&wa, r, c);
+		ele_value = wa.ele;
+		heap_add(r, c, ele_value);
+		FLAG_SET(af.flag, INLISTFLAG);
+		FLAG_SET(af.flag, EDGEFLAG);
+		af.asp = asp_value;
+		seg_put(&aspflag, (char *)&af, r, c);
+		continue;
+	    }
+
+	    /* any neighbour NULL ? */
+	    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+		/* get r, c (r_nbr, c_nbr) for neighbours */
+		r_nbr = r + nextdr[ct_dir];
+		c_nbr = c + nextdc[ct_dir];
+
+		seg_get(&aspflag, (char *)&af_nbr, r_nbr, c_nbr);
+		is_null = FLAG_GET(af_nbr.flag, NULLFLAG);
+
+		if (is_null) {
+		    asp_value = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+		    seg_get(&watalt, (char *)&wa, r, c);
+		    ele_value = wa.ele;
+		    heap_add(r, c, ele_value);
+		    FLAG_SET(af.flag, INLISTFLAG);
+		    FLAG_SET(af.flag, EDGEFLAG);
+		    af.asp = asp_value;
+		    seg_put(&aspflag, (char *)&af, r, c);
+
+		    break;
+		}
+	    }
+	    if (asp_value) /* some neighbour was NULL, point added to list */
+		continue;
+	    
+	    /* real depression ? */
+	    if (depr_fd >= 0) {
+		if (!Rast_is_c_null_value(&depr_buf[c]) && depr_buf[c] != 0) {
+		    seg_get(&watalt, (char *)&wa, r, c);
+		    ele_value = wa.ele;
+		    heap_add(r, c, ele_value);
+		    FLAG_SET(af.flag, INLISTFLAG);
+		    FLAG_SET(af.flag, DEPRFLAG);
+		    af.asp = asp_value;
+		    seg_put(&aspflag, (char *)&af, r, c);
+		    n_depr_cells++;
+		}
+	    }
+	}
+    }
+    G_percent(nrows, nrows, 2);	/* finish it */
+
+    if (depr_fd >= 0) {
+	Rast_close(depr_fd);
+	G_free(depr_buf);
+    }
+
+    G_debug(1, "%lld edge cells", heap_size - n_depr_cells);
+    if (n_depr_cells)
+	G_debug(1, "%lld cells in depressions", n_depr_cells);
+
+    return 1;
+}

raster/r.stream.extract/load.c

@@ -0,0 +1,165 @@
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+/* need elevation map, do A* search on elevation like for r.watershed */
+
+int ele_round(double x)
+{
+    if (x >= 0.0)
+	return x + .5;
+    else
+	return x - .5;
+}
+
+/*
+ * loads elevation and optional flow accumulation map to memory and
+ * gets start points for A* Search
+ * start points are edges
+ */
+int load_maps(int ele_fd, int acc_fd)
+{
+    int r, c;
+    void *ele_buf, *ptr, *acc_buf = NULL, *acc_ptr = NULL;
+    CELL ele_value, *stream_id;
+    DCELL dvalue, acc_value;
+    size_t ele_size, acc_size = 0;
+    int ele_map_type, acc_map_type = 0;
+    WAT_ALT *wabuf;
+
+    ASP_FLAG *afbuf;
+
+    if (acc_fd < 0)
+	G_message(_("Loading elevation map..."));
+    else
+	G_message(_("Loading input maps..."));
+
+    n_search_points = n_points = 0;
+
+    ele_map_type = Rast_get_map_type(ele_fd);
+    ele_size = Rast_cell_size(ele_map_type);
+    ele_buf = Rast_allocate_buf(ele_map_type);
+
+    if (ele_buf == NULL) {
+	G_warning(_("Could not allocate memory"));
+	return -1;
+    }
+
+    if (acc_fd >= 0) {
+	acc_map_type = Rast_get_map_type(acc_fd);
+	acc_size = Rast_cell_size(acc_map_type);
+	acc_buf = Rast_allocate_buf(acc_map_type);
+	if (acc_buf == NULL) {
+	    G_warning(_("Could not allocate memory"));
+	    return -1;
+	}
+    }
+
+    ele_scale = 1;
+    if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
+	ele_scale = 1000;	/* should be enough to do the trick */
+
+    wabuf = G_malloc(ncols * sizeof(WAT_ALT));
+    afbuf = G_malloc(ncols * sizeof(ASP_FLAG));
+    stream_id = G_malloc(ncols * sizeof(CELL));
+
+    G_debug(1, "start loading %d rows, %d cols", nrows, ncols);
+    for (r = 0; r < nrows; r++) {
+
+	G_percent(r, nrows, 2);
+
+	Rast_get_row(ele_fd, ele_buf, r, ele_map_type);
+	ptr = ele_buf;
+
+	if (acc_fd >= 0) {
+	    Rast_get_row(acc_fd, acc_buf, r, acc_map_type);
+	    acc_ptr = acc_buf;
+	}
+
+	for (c = 0; c < ncols; c++) {
+
+	    afbuf[c].flag = 0;
+	    afbuf[c].asp = 0;
+	    stream_id[c] = 0;
+
+	    /* check for masked and NULL cells */
+	    if (Rast_is_null_value(ptr, ele_map_type)) {
+		FLAG_SET(afbuf[c].flag, NULLFLAG);
+		FLAG_SET(afbuf[c].flag, INLISTFLAG);
+		FLAG_SET(afbuf[c].flag, WORKEDFLAG);
+		FLAG_SET(afbuf[c].flag, WORKED2FLAG);
+		Rast_set_c_null_value(&ele_value, 1);
+		/* flow accumulation */
+		if (acc_fd >= 0) {
+		    if (!Rast_is_null_value(acc_ptr, acc_map_type))
+			G_fatal_error(_("Elevation map is NULL but accumulation map is not NULL!"));
+		}
+		Rast_set_d_null_value(&acc_value, 1);
+	    }
+	    else {
+		switch (ele_map_type) {
+		case CELL_TYPE:
+		    ele_value = *((CELL *) ptr);
+		    break;
+		case FCELL_TYPE:
+		    dvalue = *((FCELL *) ptr);
+		    dvalue *= ele_scale;
+		    ele_value = ele_round(dvalue);
+		    break;
+		case DCELL_TYPE:
+		    dvalue = *((DCELL *) ptr);
+		    dvalue *= ele_scale;
+		    ele_value = ele_round(dvalue);
+		    break;
+		}
+		if (acc_fd < 0)
+		    acc_value = 1;
+		else {
+		    if (Rast_is_null_value(acc_ptr, acc_map_type)) {
+			/* can this be ok after weighing ? */
+			G_fatal_error(_("Accumulation map is NULL but elevation map is not NULL!"));
+		    }
+
+		    switch (acc_map_type) {
+		    case CELL_TYPE:
+			acc_value = *((CELL *) acc_ptr);
+			break;
+		    case FCELL_TYPE:
+			acc_value = *((FCELL *) acc_ptr);
+			break;
+		    case DCELL_TYPE:
+			acc_value = *((DCELL *) acc_ptr);
+			break;
+		    }
+		}
+
+		n_points++;
+	    }
+
+	    wabuf[c].wat = acc_value;
+	    wabuf[c].ele = ele_value;
+	    ptr = G_incr_void_ptr(ptr, ele_size);
+	    if (acc_fd >= 0)
+		acc_ptr = G_incr_void_ptr(acc_ptr, acc_size);
+	}
+	seg_put_row(&watalt, (char *) wabuf, r);
+	seg_put_row(&aspflag, (char *) afbuf, r);
+	cseg_put_row(&stream, stream_id, r);
+    }
+    G_percent(nrows, nrows, 1);	/* finish it */
+
+    Rast_close(ele_fd);
+    G_free(ele_buf);
+    G_free(wabuf);
+    G_free(afbuf);
+    G_free(stream_id);
+
+    if (acc_fd >= 0) {
+	Rast_close(acc_fd);
+	G_free(acc_buf);
+    }
+    
+    G_debug(1, "%lld non-NULL cells", n_points);
+
+    return n_points;
+}

raster/r.stream.extract/local_proto.h

@@ -0,0 +1,95 @@
+
+#ifndef __LOCAL_PROTO_H__
+#define __LOCAL_PROTO_H__
+
+#include <grass/raster.h>
+#include "flag.h"
+#include "seg.h"
+
+#define GW_LARGE_INT off_t
+
+#define INDEX(r, c) ((r) * ncols + (c))
+#define MAXDEPTH 1000     /* maximum supported tree depth of stream network */
+
+#define POINT       struct a_point
+POINT {
+    int r, c;
+};
+
+#define HEAP_PNT    struct heap_point
+HEAP_PNT {
+   GW_LARGE_INT added;
+   CELL ele;
+   POINT pnt;
+};
+
+#define WAT_ALT    struct wat_altitude
+WAT_ALT {
+   CELL ele;
+   DCELL wat;
+};
+
+#define ASP_FLAG    struct aspect_flag
+ASP_FLAG {
+   char asp;
+   char flag;
+};
+
+struct snode
+{
+    int r, c;
+    int id;
+    int n_trib;           /* number of tributaries */
+    int n_trib_total;     /* number of all upstream stream segments */
+    int n_alloc;          /* n allocated tributaries */
+    int *trib;
+};
+
+/* extern variables */
+
+extern int nrows, ncols;
+extern GW_LARGE_INT n_search_points, n_points, nxt_avail_pt;
+extern GW_LARGE_INT heap_size;
+extern GW_LARGE_INT n_stream_nodes, n_alloc_nodes;
+extern POINT *outlets;
+extern struct snode *stream_node;
+extern GW_LARGE_INT n_outlets, n_alloc_outlets;
+extern char drain[3][3];
+extern char sides;
+extern int c_fac;
+extern int ele_scale;
+extern int have_depressions;
+
+extern SSEG search_heap;
+extern SSEG astar_pts;
+extern SSEG watalt, aspflag;
+/* extern BSEG bitflags, asp; */
+extern CSEG stream;
+
+/* load.c */
+int load_maps(int, int);
+
+/* init_search.c */
+int init_search(int);
+
+/* do_astar.c */
+int do_astar(void);
+GW_LARGE_INT heap_add(int, int, CELL);
+
+/* streams.c */
+int do_accum(double);
+int extract_streams(double, double, int);
+
+/* thin.c */
+int thin_streams(void);
+
+/* basins.c */
+int basin_borders(void);
+
+/* del_streams.c */
+int del_streams(int);
+
+/* close.c */
+int close_maps(char *, char *, char *);
+
+#endif /* __LOCAL_PROTO_H__ */

raster/r.stream.extract/main.c

@@ -0,0 +1,444 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.stream.extract
+ * AUTHOR(S):    Markus Metz <markus.metz.giswork gmail.com>
+ * PURPOSE:      Hydrological analysis
+ *               Extracts stream networks from accumulation raster with
+ *               given threshold
+ * COPYRIGHT:    (C) 1999-2009 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
+ *               for details.
+ *
+ *****************************************************************************/
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+#include <math.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+/* global variables */
+int nrows, ncols;
+GW_LARGE_INT n_search_points, n_points, nxt_avail_pt;
+GW_LARGE_INT heap_size;
+GW_LARGE_INT n_stream_nodes, n_alloc_nodes;
+POINT *outlets;
+struct snode *stream_node;
+GW_LARGE_INT n_outlets, n_alloc_outlets;
+char drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
+char sides;
+int c_fac;
+int ele_scale;
+int have_depressions;
+
+SSEG search_heap;
+SSEG astar_pts;
+SSEG watalt, aspflag;
+/* BSEG bitflags, asp; */
+CSEG stream;
+
+CELL *astar_order;
+
+int main(int argc, char *argv[])
+{
+    struct
+    {
+	struct Option *ele, *acc, *depression;
+	struct Option *threshold, *d8cut;
+	struct Option *mont_exp;
+	struct Option *min_stream_length;
+	struct Option *memory;
+    } input;
+    struct
+    {
+	struct Option *stream_rast;
+	struct Option *stream_vect;
+	struct Option *dir_rast;
+    } output;
+    struct GModule *module;
+    int ele_fd, acc_fd, depr_fd;
+    double threshold, d8cut, mont_exp;
+    int min_stream_length = 0, memory;
+    int seg_cols, seg_rows;
+    double seg2kb;
+    int num_open_segs, num_open_array_segs, num_seg_total;
+    double memory_divisor, heap_mem, disk_space;
+    const char *mapset;
+
+    G_gisinit(argv[0]);
+
+    /* Set description */
+    module = G_define_module();
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("hydrology"));
+    module->description = _("Stream network extraction");
+
+    input.ele = G_define_standard_option(G_OPT_R_INPUT);
+    input.ele->key = "elevation";
+    input.ele->label = _("Elevation map");
+    input.ele->description = _("Elevation on which entire analysis is based");
+
+    input.acc = G_define_standard_option(G_OPT_R_INPUT);
+    input.acc->key = "accumulation";
+    input.acc->label = _("Accumulation map");
+    input.acc->required = NO;
+    input.acc->description =
+	_("Stream extraction will use provided accumulation instead of calculating it anew");
+
+    input.depression = G_define_standard_option(G_OPT_R_INPUT);
+    input.depression->key = "depression";
+    input.depression->label = _("Map with real depressions");
+    input.depression->required = NO;
+    input.depression->description =
+	_("Streams will not be routed out of real depressions");
+
+    input.threshold = G_define_option();
+    input.threshold->key = "threshold";
+    input.threshold->label = _("Minimum flow accumulation for streams");
+    input.threshold->description = _("Must be > 0");
+    input.threshold->required = YES;
+    input.threshold->type = TYPE_DOUBLE;
+
+    input.d8cut = G_define_option();
+    input.d8cut->key = "d8cut";
+    input.d8cut->label = _("Use SFD above this threshold");
+    input.d8cut->description =
+	_("If accumulation is larger than d8cut, SFD is used instead of MFD."
+	  " Applies only if no accumulation map is given.");
+    input.d8cut->required = NO;
+    input.d8cut->answer = "infinity";
+    input.d8cut->type = TYPE_DOUBLE;
+
+    input.mont_exp = G_define_option();
+    input.mont_exp->key = "mexp";
+    input.mont_exp->type = TYPE_DOUBLE;
+    input.mont_exp->required = NO;
+    input.mont_exp->answer = "0";
+    input.mont_exp->label =
+	_("Montgomery exponent for slope, disabled with 0");
+    input.mont_exp->description =
+	_("Montgomery: accumulation is multiplied with pow(slope,mexp) and then compared with threshold.");
+
+    input.min_stream_length = G_define_option();
+    input.min_stream_length->key = "stream_length";
+    input.min_stream_length->type = TYPE_INTEGER;
+    input.min_stream_length->required = NO;
+    input.min_stream_length->answer = "0";
+    input.min_stream_length->label =
+	_("Delete stream segments shorter than stream_length cells.");
+    input.min_stream_length->description =
+	_("Applies only to first-order stream segments (springs/stream heads).");
+
+    input.memory = G_define_option();
+    input.memory->key = "memory";
+    input.memory->type = TYPE_INTEGER;
+    input.memory->required = NO;
+    input.memory->answer = "300";
+    input.memory->description = _("Maximum memory to be used in MB");
+
+    output.stream_rast = G_define_standard_option(G_OPT_R_OUTPUT);
+    output.stream_rast->key = "stream_rast";
+    output.stream_rast->description =
+	_("Output raster map with unique stream ids");
+    output.stream_rast->required = NO;
+    output.stream_rast->guisection = _("Output options");
+
+    output.stream_vect = G_define_standard_option(G_OPT_V_OUTPUT);
+    output.stream_vect->key = "stream_vect";
+    output.stream_vect->description =
+	_("Output vector with unique stream ids");
+    output.stream_vect->required = NO;
+    output.stream_vect->guisection = _("Output options");
+
+    output.dir_rast = G_define_standard_option(G_OPT_R_OUTPUT);
+    output.dir_rast->key = "direction";
+    output.dir_rast->description =
+	_("Output raster map with flow direction");
+    output.dir_rast->required = NO;
+    output.dir_rast->guisection = _("Output options");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    /***********************/
+    /*    check options    */
+    /***********************/
+
+    /* input maps exist ? */
+    if (!G_find_raster(input.ele->answer, ""))
+	G_fatal_error(_("Raster map <%s> not found"), input.ele->answer);
+
+    if (input.acc->answer) {
+	if (!G_find_raster(input.acc->answer, ""))
+	    G_fatal_error(_("Raster map <%s> not found"), input.acc->answer);
+    }
+
+    if (input.depression->answer) {
+	if (!G_find_raster(input.depression->answer, ""))
+	    G_fatal_error(_("Raster map <%s> not found"), input.depression->answer);
+	have_depressions = 1;
+    }
+    else
+	have_depressions = 0;
+
+    /* threshold makes sense */
+    threshold = atof(input.threshold->answer);
+    if (threshold <= 0)
+	G_fatal_error(_("Threshold must be > 0 but is %f"), threshold);
+
+    /* d8cut */
+    if (strcmp(input.d8cut->answer, "infinity") == 0) {
+	d8cut = DBL_MAX;
+    }
+    else {
+	d8cut = atof(input.d8cut->answer);
+	if (d8cut < 0)
+	    G_fatal_error(_("d8cut must be positive or zero but is %f"),
+			  d8cut);
+    }
+
+    /* Montgomery stream initiation */
+    if (input.mont_exp->answer) {
+	mont_exp = atof(input.mont_exp->answer);
+	if (mont_exp < 0)
+	    G_fatal_error(_("Montgomery exponent must be positive or zero but is %f"),
+			  mont_exp);
+	if (mont_exp > 3)
+	    G_warning(_("Montgomery exponent is %f, recommended range is 0.0 - 3.0"),
+		      mont_exp);
+    }
+    else
+	mont_exp = 0;
+
+    /* Minimum stream segment length */
+    if (input.min_stream_length->answer) {
+	min_stream_length = atoi(input.min_stream_length->answer);
+	if (min_stream_length < 0)
+	    G_fatal_error(_("Minimum stream length must be positive or zero but is %d"),
+			  min_stream_length);
+    }
+    else
+	min_stream_length = 0;
+	
+    if (input.memory->answer) {
+	memory = atoi(input.memory->answer);
+	if (memory <= 0)
+	    G_fatal_error(_("Memory must be positive but is %d"),
+			  memory);
+    }
+    else
+	memory = 300;
+
+    /* Check for some output map */
+    if ((output.stream_rast->answer == NULL)
+	&& (output.stream_vect->answer == NULL)
+	&& (output.dir_rast->answer == NULL)) {
+	G_fatal_error(_("Sorry, you must choose at least one output map."));
+    }
+
+    /*********************/
+    /*    preparation    */
+    /*********************/
+
+    /* open input maps */
+    mapset = G_find_raster2(input.ele->answer, "");
+    ele_fd = Rast_open_old(input.ele->answer, mapset);
+    if (ele_fd < 0)
+	G_fatal_error(_("Could not open input map %s"), input.ele->answer);
+
+    if (input.acc->answer) {
+	mapset = G_find_raster2(input.acc->answer, "");
+	acc_fd = Rast_open_old(input.acc->answer, mapset);
+	if (acc_fd < 0)
+	    G_fatal_error(_("Could not open input map %s"),
+			  input.acc->answer);
+    }
+    else
+	acc_fd = -1;
+
+    if (input.depression->answer) {
+	mapset = G_find_raster2(input.depression->answer, "");
+	depr_fd = Rast_open_old(input.depression->answer, mapset);
+	if (depr_fd < 0)
+	    G_fatal_error(_("Could not open input map %s"),
+			  input.depression->answer);
+    }
+    else
+	depr_fd = -1;
+
+    /* set global variables */
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+    sides = 8;	/* not a user option */
+    c_fac = 5;	/* not a user option, MFD covergence factor 5 gives best results */
+
+    /* segment structures */
+    seg_rows = seg_cols = 64;
+    seg2kb = seg_rows * seg_cols / 1024.;
+
+    /* balance segment files */
+    /* elevation + accumulation: * 2 */
+    memory_divisor = sizeof(WAT_ALT) * 2;
+    disk_space = sizeof(WAT_ALT);
+    /* stream ids: / 2 */
+    memory_divisor += sizeof(int) / 2.;
+    disk_space += sizeof(int);
+
+    /* aspect and flags: * 2 */
+    memory_divisor += sizeof(ASP_FLAG) * 4;
+    disk_space += sizeof(ASP_FLAG);
+
+    /* astar_points: / 16 */
+    /* ideally only a few but large segments */
+    memory_divisor += sizeof(POINT) / 16.;
+    disk_space += sizeof(POINT);
+    /* heap points: / 4 */
+    memory_divisor += sizeof(HEAP_PNT) / 4.;
+    disk_space += sizeof(HEAP_PNT);
+    
+    /* KB -> MB */
+    memory_divisor *= seg2kb / 1024.;
+    disk_space *= seg2kb / 1024.;
+
+    num_open_segs = memory / memory_divisor;
+    heap_mem = num_open_segs * seg2kb * sizeof(HEAP_PNT) /
+               (4. * 1024.);
+    num_seg_total = (ncols / seg_cols + 1) * (nrows / seg_rows + 1);
+    if (num_open_segs > num_seg_total) {
+	heap_mem += (num_open_segs - num_seg_total) * memory_divisor;
+	heap_mem -= (num_open_segs - num_seg_total) * seg2kb *
+		    sizeof(HEAP_PNT) / (4. * 1024.);
+	num_open_segs = num_seg_total;
+    }
+    if (num_open_segs < 16) {
+	num_open_segs = 16;
+	heap_mem = num_open_segs * seg2kb * sizeof(HEAP_PNT) /
+	           (4. * 1024.);
+    }
+    G_verbose_message(_("%.2f%% of data are kept in memory"),
+                      100. * num_open_segs / num_seg_total);
+    disk_space *= num_seg_total;
+    if (disk_space < 1024.0)
+	G_verbose_message(_("Will need up to %.2f MB of disk space"), disk_space);
+    else
+	G_verbose_message(_("Will need up to %.2f GB (%.0f MB) of disk space"),
+	           disk_space / 1024.0, disk_space);
+
+    /* open segment files */
+    G_verbose_message(_("Creating temporary files..."));
+    seg_open(&watalt, nrows, ncols, seg_rows, seg_cols, num_open_segs * 2,
+        sizeof(WAT_ALT), 1);
+    if (num_open_segs * 2 > num_seg_total)
+	heap_mem += (num_open_segs * 2 - num_seg_total) * seg2kb *
+	            sizeof(WAT_ALT) / 1024.;
+    cseg_open(&stream, seg_rows, seg_cols, num_open_segs / 2.);
+
+    seg_open(&aspflag, nrows, ncols, seg_rows, seg_cols, num_open_segs * 2,
+        sizeof(ASP_FLAG), 1);
+/*
+    bseg_open(&asp, seg_rows, seg_cols, num_open_segs);
+    bseg_open(&bitflags, seg_rows, seg_cols, num_open_segs * 4);
+*/
+
+    if (num_open_segs * 4 > num_seg_total)
+	heap_mem += (num_open_segs * 4 - num_seg_total) * seg2kb / 1024.;
+
+    /* load maps */
+    if (load_maps(ele_fd, acc_fd) < 0)
+	G_fatal_error(_("Could not load input map(s)"));
+    else if (!n_points)
+	G_fatal_error(_("No non-NULL cells in input map(s)"));
+
+    G_debug(1, "open segments for A* points");
+    /* columns per segment */
+    seg_cols = seg_rows * seg_rows;
+    num_seg_total = n_points / seg_cols;
+    if (n_points % seg_cols > 0)
+	num_seg_total++;
+    /* no need to have more segments open than exist */
+    num_open_array_segs = num_open_segs / 16.;
+    if (num_open_array_segs > num_seg_total)
+	num_open_array_segs = num_seg_total;
+    if (num_open_array_segs < 1)
+	num_open_array_segs = 1;
+    
+    G_debug(1, "segment size for A* points: %d", seg_cols);
+    seg_open(&astar_pts, 1, n_points, 1, seg_cols, num_open_array_segs,
+	     sizeof(POINT), 1);
+
+    /* one-based d-ary search_heap with astar_pts */
+    G_debug(1, "open segments for A* search heap");
+	
+    /* allowed memory for search heap in MB */
+    G_debug(1, "heap memory %.2f MB", heap_mem);
+    /* columns per segment */
+    /* larger is faster */
+    seg_cols = seg_rows * seg_rows * seg_rows;
+    num_seg_total = n_points / seg_cols;
+    if (n_points % seg_cols > 0)
+	num_seg_total++;
+    /* no need to have more segments open than exist */
+    num_open_array_segs = (1 << 20) * heap_mem / (seg_cols * sizeof(HEAP_PNT));
+    if (num_open_array_segs > num_seg_total)
+	num_open_array_segs = num_seg_total;
+    if (num_open_array_segs < 2)
+	num_open_array_segs = 2;
+
+    G_debug(1, "A* search heap open segments %d, total %d",
+            num_open_array_segs, num_seg_total);
+    G_debug(1, "segment size for heap points: %d", seg_cols);
+    /* the search heap will not hold more than 5% of all points at any given time ? */
+    /* chances are good that the heap will fit into one large segment */
+    seg_open(&search_heap, 1, n_points + 1, 1, seg_cols,
+	     num_open_array_segs, sizeof(HEAP_PNT), 1);
+
+    /********************/
+    /*    processing    */
+    /********************/
+
+    /* initialize A* search */
+    if (init_search(depr_fd) < 0)
+	G_fatal_error(_("Could not initialize search"));
+
+    /* sort elevation and get initial stream direction */
+    if (do_astar() < 0)
+	G_fatal_error(_("Could not sort elevation map"));
+    seg_close(&search_heap);
+
+    if (acc_fd < 0) {
+	/* accumulate surface flow */
+	if (do_accum(d8cut) < 0)
+	    G_fatal_error(_("Could not calculate flow accumulation"));
+    }
+
+    /* extract streams */
+    if (extract_streams(threshold, mont_exp, acc_fd < 0) < 0)
+	G_fatal_error(_("Could not extract streams"));
+
+    seg_close(&astar_pts);
+    seg_close(&watalt);
+
+    /* thin streams */
+    if (thin_streams() < 0)
+	G_fatal_error(_("Could not thin streams"));
+
+    /* delete short streams */
+    if (min_stream_length) {
+	if (del_streams(min_stream_length) < 0)
+	    G_fatal_error(_("Could not delete short stream segments"));
+    }
+
+    /* write output maps */
+    if (close_maps(output.stream_rast->answer, output.stream_vect->answer,
+		   output.dir_rast->answer) < 0)
+	G_fatal_error(_("Could not write output maps"));
+
+    cseg_close(&stream);
+    seg_close(&aspflag);
+
+    exit(EXIT_SUCCESS);
+}

raster/r.stream.extract/r.stream.extract.html

@@ -0,0 +1,240 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.stream.extract</em> extracts streams in both raster and vector
+format from a required input <em>elevation</em> map and optional input
+<em>accumulation</em> map.
+
+<h2>OPTIONS</h2>
+
+<dl>
+<dt><em>elevation</em> 
+<dd>Input map, required: Elevation on which the entire analysis is based.
+NULL (nodata) cells are ignored, zero and negative values are valid
+elevation data. Gaps in the elevation map that are located within the
+area of interest must be filled beforehand, e.g. with
+<em>r.fillnulls</em>, to avoid distortions.
+<p>
+<dt><em>accumulation</em>
+<dd>Input map, optional: Accumulation values of the provided
+<em>accumulation</em> map are used and not calculated from the input
+<em>elevation</em> map. If <em>accumulation</em> is given,
+<em>elevation</em> must be exactly the same map used to calculate
+<em>accumulation</em>. If <em>accumulation</em> was calculated with
+<a href="r.terraflow.html">r.terraflow</a>, the filled elevation output
+of r.terraflow must be used. Further on, the current region should be 
+aligned to the <em>accumulation map</em>. Flow direction is
+first calculated from <em>elevation</em> and then adjusted to
+<em>accumulation</em>. It is not necessary to provide <em>accumulation</em>
+as the number of cells, it can also be the optionally adjusted or
+weighed total contributing area in square meters or any other unit. 
+When an original flow accumulation map is adjusted or weighed, the 
+adjustment or weighing should not convert valid accumulation values to 
+NULL (nodata) values.
+<p>
+<dt><em>depression</em> 
+<dd>Input map, optional: All non-NULL and non-zero cells will be regarded
+as real depressions. Streams will not be routed out of depressions. If an
+area is marked as depression but the elevation model has no depression
+at this location, streams will not stop there. If a flow accumulation map
+and a map with real depressions are provided, the flow accumulation map
+must match the depression map such that flow is not distributed out of 
+the indicated depressions. It is recommended to use internally computed
+flow accumulation if a depression map is provided.
+<p>
+<dt><em>threshold</em>
+<dd>Required: <em>threshold</em> for stream initiation by overland flow:
+the minumum (optionally modifed) flow accumulation value that will initiate
+a new stream. If Montgomery's method for channel initiation is used, the
+cell value of the accumulation input map is multiplied by
+(tan(local slope))<sup>mexp</sup> and then compared to <em>threshold</em>.
+<p>
+<dt><em>d8cut</em>
+<dd>Minimum amount of overland flow (accumulation) when SFD (D8) will be
+used instead of MFD (FD8) to calculate flow accumulation. Only applies
+if no accumulation map is provided. Setting to 0 disables MFD completely.
+<p>
+<dt><em>mexp</em>
+<dd>Use the method of Montgomery and Foufoula-Georgiou (1993) to
+initiate a stream with exponent <em>mexp</em>. The cell value of the
+accumulation input map is multiplied by (tan(local slope))<sup>mexp</sup>
+and then compared to <em>threshold</em>. If threshold is reached or
+exceeded, a new stream is initiated. The default value 0 disables
+Montgomery. Montgomery and Foufoula-Georgiou (1993) generally recommend
+to use 2.0 as exponent. <em>mexp</em> values closer to 0 will produce
+streams more similar to streams extracted with Montgomery disabled.
+Larger <em>mexp</em> values decrease the number of streams in flat areas
+and increase the number of streams in steep areas. If <em>weight</em> is
+given, the weight is applied first.
+<p>
+<dt><em>stream_length</em>
+<dd>Minimum stream length in number of cells for first-order (head/spring)
+stream segments. All first-order stream segments shorter than
+<em>stream_length</em> will be deleted.
+
+<p>
+<dt><em>stream_rast</em>
+<dd>Output raster map with extracted streams. Cell values encode a 
+unique ID for each stream segment.
+<p>
+<dt><em>stream_vect</em>
+<dd>Output vector map with extracted stream segments and points. Points
+are written at the start location of each stream segment and at the
+outlet of a stream network. In layer 1, categories are unique IDs,
+identical to the cell value of the raster output. The attribute table
+for layer 1 holds information about the type of stream segment: start
+segment, or intermediate segment with tributaries. Columns are cat int,
+stream_type varchar(), type_code int. The encoding for type_code is 0 =
+start, 1 = intermediate. In layer 2, categories are identical to
+type_code in layer 1 with additional category 2 = outlet for outlet
+points. Points with category 1 = intermediate in layer 2 are at the
+location of confluences.
+<p>
+<dt><em>direction</em>
+<dd>Output raster map with flow direction for all non-NULL cells in
+input elevation. Flow direction is of D8 type with a range of 1 to 8.
+Multiplying values with 45 gives degrees CCW from East. Flow direction
+was adjusted during thinning, taking shortcuts and skipping cells that
+were eliminated by the thinning procedure.
+</dl>
+
+<h2>NOTES</h2>
+
+<h4>Stream extraction</h4>
+If no accumulation input map is provided, flow accumulation is
+determined with a hydrological anaylsis similar to
+<a href="r.watershed.html">r.watershed</a>. The algorithm is
+MFD (FD8) after Holmgren 1994, as for
+<a href="r.watershed.html">r.watershed</a>. The <em>threshold</em>
+option determines the number of streams and detail of stream networks.
+Whenever flow accumulation reaches <em>threshold</em>, a new stream is
+started and traced downstream to its outlet point. As for
+<a href="r.watershed.html">r.watershed</a>, flow accumulation is
+calculated as the number of cells draining through a cell.
+
+<h4>Weighed flow accumulation</h4>
+Flow accumulation can be calculated first, e.g. with
+<a href="r.watershed.html">r.watershed</a>, and then modified before
+using it as input for <em>r.stream.extract</em>. In its general form, a
+weighed accumulation map is generated by first creating a weighing map
+and then multiplying the accumulation map with the weighing map using
+<a href="r.mapcalc.html">r.mapcalc</a>. It is highly recommended to
+evaluate the weighed flow accumulation map first, before using it as
+input for <em>r.stream.extract</em>.
+<p>
+This allows e.g. to decrease the number of streams in dry areas and
+increase the number of streams in wet areas by setting <em>weight</em>
+to smaller than 1 in dry areas and larger than 1 in wet areas.
+<p>
+Another possibility is to restrict channel initiation to valleys
+determined from terrain morphology. Valleys can be determined with
+<a href="r.param.scale.html">r.param.scale</a> param=crosc
+(cross-sectional or tangential curvature). Curvature values &lt; 0
+indicate concave features, i.e. valleys. The size of the processing
+window determines whether narrow or broad valleys will be identified
+(See example below).
+
+<h4>Defining a region of interest</h4>
+The stream extraction procedure can be restricted to a certain region of 
+interest, e.g. a subbasin, by setting the computational region with 
+<em>g.region</em> and/or creating a MASK. Such region of interest should 
+be a complete catchment area, complete in the sense that the complete 
+area upstream of an outlet point is included and buffered with at least 
+one cell.
+
+<h4>Stream output</h4>
+The <em>stream_rast</em> output raster and vector contains stream
+segments with unique IDs. Note that these IDs are different from the IDs
+assigned by <a href="r.watershed.html">r.watershed</a>. The vector
+output also contains points at the location of the start of a stream
+segment, at confluences and at stream network outlet locations.
+<p>
+
+<h2>EXAMPLE</h2>
+This example is based on the elevation map <em>elev_ned_30m</em> in the
+North Carolina sample dataset and uses valleys determined with
+<a href="r.param.scale.html">r.param.scale</a> to weigh an accumulation
+map produced with <a href="r.watershed.html">r.watershed</a>.
+
+<pre>
+# set region
+g.region -p rast=elev_ned_30m@PERMANENT
+
+# calculate flow accumulation
+r.watershed ele=elev_ned_30m@PERMANENT acc=elevation.10m.acc
+
+# curvature to get narrow valleys
+r.param.scale input=elev_ned_30m@PERMANENT output=tangential_curv_5 size=5 param=crosc
+
+# curvature to get a bit broader valleys
+r.param.scale input=elev_ned_30m@PERMANENT output=tangential_curv_7 size=7 param=crosc
+
+# curvature to get broad valleys
+r.param.scale input=elev_ned_30m@PERMANENT output=tangential_curv_11 size=11 param=crosc
+
+# create weight map
+r.mapcalc "weight = if(tangential_curv_5 < 0, -100 * tangential_curv_5, \
+                    if(tangential_curv_7 < 0, -100 * tangential_curv_7, \
+		    if(tangential_curv_11 < 0, -100 * tangential_curv_11, 0.000001)))"
+
+# weigh accumulation map
+r.mapcalc "elev_ned_30m.acc.weighed = elev_ned_30m.acc * weight"
+
+# copy color table from original accumulation map
+r.colors map=elev_ned_30m.acc.weighed raster=elev_ned_30m.acc
+</pre>
+
+Display both the original and the weighed accumulation map.
+<br>
+Compare them and proceed if the weighed accumulation map makes sense.
+
+<pre>
+# extract streams
+r.stream.extract elevation=elev_ned_30m@PERMANENT \
+                 accumulation=elev_ned_30m.acc.weighed \
+		 threshold=1000 \
+		 stream_rast=elev_ned_30m.streams
+
+# extract streams using the original accumulation map
+r.stream.extract elevation=elev_ned_30m@PERMANENT \
+                 accumulation=elev_ned_30m.acc \
+		 threshold=1000 \
+		 stream_rast=elev_ned_30m.streams.noweight
+</pre>
+
+Now display both stream maps and decide which one is more realistic.
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.terraflow.html">r.terraflow</a>,
+<a href="r.param.scale.html">r.param.scale</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="r.thin.html">r.thin</a>,
+<a href="r.to.vect.html">r.to.vect</a>
+</em>
+
+<h2>REFERENCES</h2>
+Ehlschlaeger, C. (1989). <i>Using the A<sup>T</sup> Search Algorithm
+to Develop Hydrologic Models from Digital Elevation Data</i>,
+<b>Proceedings of International Geographic Information Systems (IGIS)
+Symposium '89</b>, pp 275-281 (Baltimore, MD, 18-19 March 1989).<br>
+URL: <a href="http://faculty.wiu.edu/CR-Ehlschlaeger2/older/IGIS/paper.html">
+http://faculty.wiu.edu/CR-Ehlschlaeger2/older/IGIS/paper.html</a>
+
+<p>
+Holmgren, P. (1994). <i>Multiple flow direction algorithms for runoff 
+modelling in grid based elevation models: An empirical evaluation.</i>
+<b>Hydrological Processes</b> Vol 8(4), pp 327-334.<br>
+DOI: <a href="http://dx.doi.org/10.1002/hyp.3360080405">10.1002/hyp.3360080405</a>
+
+<p>
+Montgomery, D.R., Foufoula-Georgiou, E. (1993). <i>Channel network source
+representation using digital elevation models.</i>
+<b>Water Resources Research</b> Vol 29(12), pp 3925-3934. 
+
+<h2>AUTHOR</h2>
+Markus Metz
+
+<p>
+<i>Last changed: $Date$</i>

raster/r.stream.extract/seg.c

@@ -0,0 +1,112 @@
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/glocale.h>
+#include "seg.h"
+
+int
+seg_open(SSEG *sseg, int nrows, int ncols, int row_in_seg, int col_in_seg,
+	 int nsegs_in_memory, int size_struct, int fill)
+{
+    char *filename;
+    int errflag;
+    int fd;
+
+    sseg->filename = NULL;
+    sseg->fd = -1;
+
+    filename = G_tempfile();
+    if (-1 == (fd = creat(filename, 0666))) {
+	G_warning(_("seg_open(): unable to create segment file"));
+	return -2;
+    }
+    if (fill)
+	errflag = segment_format(fd, nrows, ncols, row_in_seg,
+	                         col_in_seg, size_struct);
+    else
+	errflag = segment_format_nofill(fd, nrows, ncols, row_in_seg,
+	                         col_in_seg, size_struct);
+
+    if (0 > errflag) {
+	close(fd);
+	unlink(filename);
+	if (errflag == -1) {
+	    G_warning(_("seg_open(): could not write segment file"));
+	    return -1;
+	}
+	else {
+	    G_warning(_("seg_open(): illegal configuration parameter(s)"));
+	    return -3;
+	}
+    }
+    close(fd);
+    if (-1 == (fd = open(filename, 2))) {
+	unlink(filename);
+	G_warning(_("seg_open(): unable to re-open segment file"));
+	return -4;
+    }
+    if (0 > (errflag = segment_init(&(sseg->seg), fd, nsegs_in_memory))) {
+	close(fd);
+	unlink(filename);
+	if (errflag == -1) {
+	    G_warning(_("seg_open(): could not read segment file"));
+	    return -5;
+	}
+	else {
+	    G_warning(_("seg_open(): out of memory"));
+	    return -6;
+	}
+    }
+    sseg->filename = filename;
+    sseg->fd = fd;
+    return 0;
+}
+
+int seg_close(SSEG *sseg)
+{
+    segment_release(&(sseg->seg));
+    close(sseg->fd);
+    unlink(sseg->filename);
+    return 0;
+}
+
+int seg_put(SSEG *sseg, char *value, int row, int col)
+{
+    if (segment_put(&(sseg->seg), value, row, col) < 0) {
+	G_warning(_("seg_put(): could not write segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int seg_put_row(SSEG *sseg, char *value, int row)
+{
+    if (segment_put_row(&(sseg->seg), value, row) < 0) {
+	G_warning(_("seg_put_row(): could not write segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int seg_get(SSEG *sseg, char *value, int row, int col)
+{
+    if (segment_get(&(sseg->seg), value, row, col) < 0) {
+	G_warning(_("seg_get(): could not read segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int seg_get_row(SSEG *sseg, char *value, int row)
+{
+    if (segment_get_row(&(sseg->seg), value, row) < 0) {
+	G_warning(_("seg_get_row(): could not read segment file"));
+	return -1;
+    }
+    return 0;
+}
+
+int seg_flush(SSEG *sseg)
+{
+    segment_flush(&(sseg->seg));
+    return 0;
+}

raster/r.stream.extract/seg.h

@@ -0,0 +1,78 @@
+#ifndef __SEG_H__
+#define __SEG_H__
+
+#include <grass/raster.h>
+#include <grass/segment.h>
+
+#define CSEG struct _c_s_e_g_
+CSEG {
+    SEGMENT seg;		/* segment structure */
+    int fd;			/* fd for reading/writing segment file */
+    char *filename;		/* name of segment file */
+    char *name;			/* raster map read into segment file */
+    char *mapset;
+};
+
+#define DSEG struct _d_s_e_g_
+DSEG {
+    SEGMENT seg;		/* segment structure */
+    int fd;			/* fd for reading/writing segment file */
+    char *filename;		/* name of segment file */
+    char *name;			/* raster map read into segment file */
+    char *mapset;
+};
+
+#define BSEG struct _b_s_e_g_
+BSEG {
+    SEGMENT seg;		/* segment structure */
+    int fd;			/* fd for reading/writing segment file */
+    char *filename;		/* name of segment file */
+    char *name;			/* raster map read into segment file */
+    char *mapset;
+};
+
+#define SSEG struct _s_s_e_g_
+SSEG {
+    SEGMENT seg;		/* segment structure */
+    int fd;			/* fd for reading/writing segment file */
+    char *filename;		/* name of segment file */
+};
+
+/* bseg.c */
+int bseg_close(BSEG *);
+int bseg_get(BSEG *, char *, int, int);
+int bseg_open(BSEG *, int, int, int);
+int bseg_put(BSEG *, char *, int, int);
+int bseg_put_row(BSEG *, char *, int);
+int bseg_read_raster(BSEG *, char *, char *);
+int bseg_write_raster(BSEG *, char *);
+
+/* cseg.c */
+int cseg_close(CSEG *);
+int cseg_get(CSEG *, CELL *, int, int);
+int cseg_open(CSEG *, int, int, int);
+int cseg_put(CSEG *, CELL *, int, int);
+int cseg_put_row(CSEG *, CELL *, int);
+int cseg_read_raster(CSEG *, char *, char *);
+int cseg_write_raster(CSEG *, char *);
+
+/* dseg.c */
+int dseg_close(DSEG *);
+int dseg_get(DSEG *, double *, int, int);
+int dseg_open(DSEG *, int, int, int);
+int dseg_put(DSEG *, double *, int, int);
+int dseg_put_row(DSEG *, double *, int);
+int dseg_read_raster(DSEG *, char *, char *);
+int dseg_write_raster(DSEG *, char *);
+
+/* seg.c */
+int seg_close(SSEG *);
+int seg_get(SSEG *, char *, int, int);
+int seg_open(SSEG *, int, int, int, int, int, int, int);
+int seg_put(SSEG *, char *, int, int);
+int seg_put_row(SSEG *, char *, int);
+int seg_get(SSEG *, char *, int, int);
+int seg_get_row(SSEG *, char *, int);
+int seg_flush(SSEG *);
+
+#endif /* __SEG_H__ */

raster/r.stream.extract/streams.c

@@ -0,0 +1,713 @@
+#include <stdlib.h>
+#include <math.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+double mfd_pow(double base)
+{
+    int i;
+    double result = base;
+
+    for (i = 2; i <= c_fac; i++) {
+	result *= base;
+    }
+    return result;
+}
+
+static int continue_stream(CELL stream_id, int r_max, int c_max,
+		    int *stream_no)
+{
+    CELL curr_stream, stream_nbr, old_stream;
+    int r_nbr, c_nbr;
+    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    int stream_node_step = 1000;
+    ASP_FLAG af;
+
+    G_debug(3, "continue stream");
+    
+    cseg_get(&stream, &curr_stream, r_max, c_max);
+
+    if (curr_stream <= 0) {
+	/* no confluence, just continue */
+	G_debug(3, "no confluence, just continue stream");
+	curr_stream = stream_id;
+	cseg_put(&stream, &curr_stream, r_max, c_max);
+	seg_get(&aspflag, (char *)&af, r_max, c_max);
+	FLAG_SET(af.flag, STREAMFLAG);
+	seg_put(&aspflag, (char *)&af, r_max, c_max);
+	return 0;
+    }
+
+    G_debug(3, "confluence");
+	    
+    /* new confluence */
+    if (stream_node[curr_stream].r != r_max ||
+	stream_node[curr_stream].c != c_max) {
+	size_t new_size;
+	
+	G_debug(3, "new confluence");
+	/* set new stream id */
+	(*stream_no)++;
+	/* add stream node */
+	if (*stream_no >= n_alloc_nodes - 1) {
+	    n_alloc_nodes += stream_node_step;
+	    stream_node =
+		(struct snode *)G_realloc(stream_node,
+					  n_alloc_nodes *
+					  sizeof(struct snode));
+	}
+	stream_node[*stream_no].r = r_max;
+	stream_node[*stream_no].c = c_max;
+	stream_node[*stream_no].id = *stream_no;
+	stream_node[*stream_no].n_trib = 0;
+	stream_node[*stream_no].n_trib_total = 0;
+	stream_node[*stream_no].n_alloc = 0;
+	stream_node[*stream_no].trib = NULL;
+	n_stream_nodes++;
+
+	/* debug */
+	if (n_stream_nodes != *stream_no)
+	    G_warning(_("BUG: stream_no %d and n_stream_nodes %lld out of sync"),
+		      *stream_no, n_stream_nodes);
+
+	stream_node[*stream_no].n_alloc += 2;
+	new_size = stream_node[*stream_no].n_alloc * sizeof(int);
+	stream_node[*stream_no].trib =
+	    (int *)G_realloc(stream_node[*stream_no].trib, new_size);
+
+	/* add the two tributaries */
+	G_debug(3, "add tributaries");
+	stream_node[*stream_no].trib[stream_node[*stream_no].n_trib++] =
+	    curr_stream;
+	stream_node[*stream_no].trib[stream_node[*stream_no].n_trib++] =
+	    stream_id;
+
+	/* update stream IDs downstream */
+	G_debug(3, "update stream IDs downstream");
+	r_nbr = r_max;
+	c_nbr = c_max;
+	old_stream = curr_stream;
+	curr_stream = *stream_no;
+	cseg_put(&stream, &curr_stream, r_nbr, c_nbr);
+	seg_get(&aspflag, (char *)&af, r_nbr, c_nbr);
+
+	while (af.asp > 0) {
+	    r_nbr = r_nbr + asp_r[(int)af.asp];
+	    c_nbr = c_nbr + asp_c[(int)af.asp];
+	    cseg_get(&stream, &stream_nbr, r_nbr, c_nbr);
+	    if (stream_nbr != old_stream)
+		af.asp = -1;
+	    else {
+		cseg_put(&stream, &curr_stream, r_nbr, c_nbr);
+		seg_get(&aspflag, (char *)&af, r_nbr, c_nbr);
+	    }
+	}
+    }
+    else {
+	/* stream node already existing here */
+	G_debug(3, "existing confluence");
+	/* add new tributary to stream node */
+	if (stream_node[curr_stream].n_trib >=
+	    stream_node[curr_stream].n_alloc) {
+	    size_t new_size;
+
+	    stream_node[curr_stream].n_alloc += 2;
+	    new_size = stream_node[curr_stream].n_alloc * sizeof(int);
+	    stream_node[curr_stream].trib =
+		(int *)G_realloc(stream_node[curr_stream].trib, new_size);
+	}
+
+	stream_node[curr_stream].trib[stream_node[curr_stream].n_trib++] =
+	    stream_id;
+    }
+
+    G_debug(3, "%d tribs", stream_node[curr_stream].n_trib);
+    if (stream_node[curr_stream].n_trib == 1)
+	G_warning(_("BUG: stream node %d has only 1 tributary: %d"), curr_stream,
+		  stream_node[curr_stream].trib[0]);
+
+    return 1;
+}
+
+/*
+ * accumulate surface flow
+ */
+int do_accum(double d8cut)
+{
+    int r, c, dr, dc;
+    CELL ele_val, *ele_nbr;
+    DCELL value, *wat_nbr;
+    struct Cell_head window;
+    int mfd_cells, astar_not_set;
+    double *dist_to_nbr, *weight, sum_weight, max_weight;
+    double dx, dy;
+    int r_nbr, c_nbr, ct_dir, np_side;
+    int is_worked;
+    double prop;
+    int edge;
+    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+    int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+    GW_LARGE_INT workedon, killer;
+    char *flag_nbr;
+    POINT astarpoint;
+    WAT_ALT wa;
+    ASP_FLAG af, af_nbr;
+
+    G_message(_("Calculating flow accumulation..."));
+
+    /* distances to neighbours */
+    dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
+    weight = (double *)G_malloc(sides * sizeof(double));
+    flag_nbr = (char *)G_malloc(sides * sizeof(char));
+    wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
+    ele_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
+
+    G_get_set_window(&window);
+
+    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+	/* get r, c (r_nbr, c_nbr) for neighbours */
+	r_nbr = nextdr[ct_dir];
+	c_nbr = nextdc[ct_dir];
+	/* account for rare cases when ns_res != ew_res */
+	dy = abs(r_nbr) * window.ns_res;
+	dx = abs(c_nbr) * window.ew_res;
+	if (ct_dir < 4)
+	    dist_to_nbr[ct_dir] = dx + dy;
+	else
+	    dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+    }
+
+    /* distribute and accumulate */
+    for (killer = 0; killer < n_points; killer++) {
+
+	G_percent(killer, n_points, 1);
+	
+	seg_get(&astar_pts, (char *)&astarpoint, 0, killer);
+	r = astarpoint.r;
+	c = astarpoint.c;
+
+	seg_get(&aspflag, (char *)&af, r, c);
+
+	/* do not distribute flow along edges or out of real depressions */
+	if (af.asp <= 0) {
+	    FLAG_UNSET(af.flag, WORKEDFLAG);
+	    seg_put(&aspflag, (char *)&af, r, c);
+	    continue;
+	}
+
+	if (af.asp) {
+	    dr = r + asp_r[abs((int)af.asp)];
+	    dc = c + asp_c[abs((int)af.asp)];
+	}
+
+	seg_get(&watalt, (char *)&wa, r, c);
+	value = wa.wat;
+
+	/* WORKEDFLAG has been set during A* Search
+	 * reversed meaning here: 0 = done, 1 = not yet done */
+	FLAG_UNSET(af.flag, WORKEDFLAG);
+	seg_put(&aspflag, (char *)&af, r, c);
+
+	/***************************************/
+	/*  get weights for flow distribution  */
+	/***************************************/
+
+	max_weight = 0;
+	sum_weight = 0;
+	np_side = -1;
+	mfd_cells = 0;
+	astar_not_set = 1;
+	ele_val = wa.ele;
+	edge = 0;
+	/* this loop is needed to get the sum of weights */
+	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+	    /* get r_nbr, c_nbr for neighbours */
+	    r_nbr = r + nextdr[ct_dir];
+	    c_nbr = c + nextdc[ct_dir];
+	    weight[ct_dir] = -1;
+	    wat_nbr[ct_dir] = 0;
+	    ele_nbr[ct_dir] = 0;
+
+	    /* check that neighbour is within region */
+	    if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 && c_nbr < ncols) {
+
+		seg_get(&aspflag, (char *)&af_nbr, r_nbr, c_nbr);
+		flag_nbr[ct_dir] = af_nbr.flag;
+		if ((edge = FLAG_GET(flag_nbr[ct_dir], NULLFLAG)))
+		    break;
+		seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
+		wat_nbr[ct_dir] = wa.wat;
+		ele_nbr[ct_dir] = wa.ele;
+
+		/* WORKEDFLAG has been set during A* Search
+		 * reversed meaning here: 0 = done, 1 = not yet done */
+		is_worked = FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG) == 0;
+		if (is_worked == 0) {
+		    if (ele_nbr[ct_dir] <= ele_val) {
+			if (ele_nbr[ct_dir] < ele_val) {
+			    weight[ct_dir] =
+				mfd_pow((ele_val -
+					 ele_nbr[ct_dir]) / dist_to_nbr[ct_dir]);
+			}
+			if (ele_nbr[ct_dir] == ele_val) {
+			    weight[ct_dir] =
+				mfd_pow(0.5 / dist_to_nbr[ct_dir]);
+			}
+			sum_weight += weight[ct_dir];
+			mfd_cells++;
+
+			if (weight[ct_dir] > max_weight) {
+			    max_weight = weight[ct_dir];
+			}
+
+			if (dr == r_nbr && dc == c_nbr) {
+			    astar_not_set = 0;
+			}
+		    }
+		}
+		if (dr == r_nbr && dc == c_nbr)
+		    np_side = ct_dir;
+	    }
+	    else
+		edge = 1;
+	    if (edge)
+		break;
+	}
+
+	/* do not distribute flow along edges, this causes artifacts */
+	if (edge) {
+	    G_debug(3, "edge");
+	    continue;
+	}
+
+	/* honour A * path 
+	 * mfd_cells == 0: fine, SFD along A * path
+	 * mfd_cells == 1 && astar_not_set == 0: fine, SFD along A * path
+	 * mfd_cells > 0 && astar_not_set == 1: A * path not included, add to mfd_cells
+	 */
+
+	/************************************/
+	/*  distribute and accumulate flow  */
+	/************************************/
+
+	/* MFD, A * path not included, add to mfd_cells */
+	if (mfd_cells > 0 && astar_not_set == 1) {
+	    mfd_cells++;
+	    sum_weight += max_weight;
+	    weight[np_side] = max_weight;
+	}
+
+	/* use SFD (D8) if d8cut threshold exceeded */
+	if (fabs(value) > d8cut)
+	    mfd_cells = 0;
+
+	if (mfd_cells > 1) {
+	    prop = 0.0;
+	    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+		/* get r, c (r_nbr, c_nbr) for neighbours */
+		r_nbr = r + nextdr[ct_dir];
+		c_nbr = c + nextdc[ct_dir];
+
+		/* check that neighbour is within region */
+		if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
+		    c_nbr < ncols && weight[ct_dir] > -0.5) {
+		    is_worked = FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG) == 0;
+		    if (is_worked == 0) {
+
+			weight[ct_dir] = weight[ct_dir] / sum_weight;
+			/* check everything sums up to 1.0 */
+			prop += weight[ct_dir];
+
+			wa.wat = wat_nbr[ct_dir] + value * weight[ct_dir];
+			wa.ele = ele_nbr[ct_dir];
+			seg_put(&watalt, (char *)&wa, r_nbr, c_nbr);
+		    }
+		    else if (ct_dir == np_side) {
+			/* check for consistency with A * path */
+			workedon++;
+		    }
+		}
+	    }
+	    if (fabs(prop - 1.0) > 5E-6f) {
+		G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
+			  prop);
+	    }
+	}
+	/* get out of depression in SFD mode */
+	else {
+	    wa.wat = wat_nbr[np_side] + value;
+	    wa.ele = ele_nbr[np_side];
+	    seg_put(&watalt, (char *)&wa, dr, dc);
+	}
+    }
+    G_percent(1, 1, 2);
+
+    G_free(dist_to_nbr);
+    G_free(weight);
+    G_free(wat_nbr);
+    G_free(ele_nbr);
+    G_free(flag_nbr);
+
+    return 1;
+}
+
+/*
+ * extracts streams for threshold, accumulation is provided
+ */
+int extract_streams(double threshold, double mont_exp, int internal_acc)
+{
+    int r, c, dr, dc;
+    CELL is_swale, ele_val, *ele_nbr;
+    DCELL value, valued, *wat_nbr;
+    struct Cell_head window;
+    int mfd_cells, stream_cells, swale_cells;
+    double *dist_to_nbr;
+    double dx, dy;
+    int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side, max_side;
+    int is_worked;
+    double max_acc;
+    int edge, flat;
+    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+    int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+    /* sides */
+    /*
+     *  | 7 | 1 | 4 |
+     *  | 2 |   | 3 |
+     *  | 5 | 0 | 6 |
+     */
+    GW_LARGE_INT workedon, killer;
+    int stream_no = 0, stream_node_step = 1000;
+    double slope, diag;
+    char *flag_nbr;
+    POINT astarpoint;
+    WAT_ALT wa;
+    ASP_FLAG af, af_nbr;
+
+    G_message(_("Extracting streams..."));
+
+    /* init stream nodes */
+    n_alloc_nodes = stream_node_step;
+    stream_node =
+	(struct snode *)G_malloc(n_alloc_nodes * sizeof(struct snode));
+    n_stream_nodes = 0;
+
+    /* init outlet nodes */
+    n_alloc_outlets = stream_node_step;
+    outlets =
+	(POINT *)G_malloc(n_alloc_outlets * sizeof(POINT));
+    n_outlets = 0;
+
+    /* distances to neighbours */
+    dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
+    flag_nbr = (char *)G_malloc(sides * sizeof(char));
+    wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
+    ele_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
+
+    G_get_set_window(&window);
+
+    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+	/* get r, c (r_nbr, c_nbr) for neighbours */
+	r_nbr = nextdr[ct_dir];
+	c_nbr = nextdc[ct_dir];
+	/* account for rare cases when ns_res != ew_res */
+	dy = abs(r_nbr) * window.ns_res;
+	dx = abs(c_nbr) * window.ew_res;
+	if (ct_dir < 4)
+	    dist_to_nbr[ct_dir] = dx + dy;
+	else
+	    dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+    }
+
+    diag = sqrt(2);
+
+    workedon = 0;
+
+    /* extract streams */
+    for (killer =  0; killer < n_points; killer++) {
+	G_percent(killer, n_points, 1);
+	
+	seg_get(&astar_pts, (char *)&astarpoint, 0, killer);
+	r = astarpoint.r;
+	c = astarpoint.c;
+
+	seg_get(&aspflag, (char *)&af, r, c);
+	/* internal acc: SET, external acc: UNSET */
+	if (internal_acc)
+	    FLAG_SET(af.flag, WORKEDFLAG);
+	else
+	    FLAG_UNSET(af.flag, WORKEDFLAG);
+	seg_put(&aspflag, (char *)&af, r, c);
+
+	/* do not distribute flow along edges */
+	if (af.asp <= 0) {
+	    G_debug(3, "edge");
+	    is_swale = FLAG_GET(af.flag, STREAMFLAG);
+	    if (is_swale) {
+		G_debug(2, "edge outlet");
+		/* add outlet point */
+		if (n_outlets >= n_alloc_outlets) {
+		    n_alloc_outlets += stream_node_step;
+		    outlets =
+			(POINT *)G_realloc(outlets,
+						  n_alloc_outlets *
+						  sizeof(POINT));
+		}
+		outlets[n_outlets].r = r;
+		outlets[n_outlets].c = c;
+		n_outlets++;
+	    }
+
+	    if (af.asp == 0) {
+		/* can only happen with real depressions */
+		if (!have_depressions)
+		    G_fatal_error(_("Bug in stream extraction"));
+		G_debug(2, "bottom of real depression");
+	    } 
+	    continue;
+	}
+
+	if (af.asp) {
+	    dr = r + asp_r[abs((int)af.asp)];
+	    dc = c + asp_c[abs((int)af.asp)];
+	}
+	else {
+	    /* can only happen with real depressions,
+	     * but should not get to here */
+	    dr = r;
+	    dc = c;
+	}
+
+	r_nbr = r_max = dr;
+	c_nbr = c_max = dc;
+
+	seg_get(&watalt, (char *)&wa, r, c);
+	value = wa.wat;
+
+	/**********************************/
+	/*  find main drainage direction  */
+	/**********************************/
+
+	max_acc = -1;
+	max_side = np_side = -1;
+	mfd_cells = 0;
+	stream_cells = 0;
+	swale_cells = 0;
+	ele_val = wa.ele;
+	edge = 0;
+	flat = 1;
+	/* find main drainage direction */
+	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+	    /* get r_nbr, c_nbr for neighbours */
+	    r_nbr = r + nextdr[ct_dir];
+	    c_nbr = c + nextdc[ct_dir];
+	    wat_nbr[ct_dir] = 0;
+	    ele_nbr[ct_dir] = 0;
+	    flag_nbr[ct_dir] = 0;
+
+	    /* check that neighbour is within region */
+	    if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 && c_nbr < ncols) {
+
+		if (dr == r_nbr && dc == c_nbr)
+		    np_side = ct_dir;
+
+		seg_get(&aspflag, (char *)&af_nbr, r_nbr, c_nbr);
+		flag_nbr[ct_dir] = af_nbr.flag;
+		if ((edge = FLAG_GET(flag_nbr[ct_dir], NULLFLAG)))
+		    break;
+		seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
+		wat_nbr[ct_dir] = wa.wat;
+		ele_nbr[ct_dir] = wa.ele;
+
+		/* check for swale cells */
+		is_swale = FLAG_GET(flag_nbr[ct_dir], STREAMFLAG);
+		if (is_swale)
+		    swale_cells++;
+
+		/* check for stream cells */
+		valued = fabs(wat_nbr[ct_dir]);
+		/* check all upstream neighbours */
+		if (valued >= threshold && ct_dir != np_side &&
+		    ele_nbr[ct_dir] > ele_val)
+		    stream_cells++;
+
+		is_worked = FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG);
+		if (!internal_acc)
+		    is_worked = is_worked == 0;
+
+		if (is_worked == 0) {
+		    if (ele_nbr[ct_dir] != ele_val)
+			flat = 0;
+		    if (ele_nbr[ct_dir] <= ele_val) {
+
+			mfd_cells++;
+
+			/* set main drainage direction */
+			if (valued >= max_acc) {
+			    max_acc = valued;
+			    r_max = r_nbr;
+			    c_max = c_nbr;
+			    max_side = ct_dir;
+			}
+		    }
+		}
+		else if (ct_dir == np_side && !edge) {
+		    /* check for consistency with A * path */
+		    workedon++;
+		}
+	    }
+	    else
+		edge = 1;
+	    if (edge)
+		break;
+	}
+
+	is_swale = FLAG_GET(af.flag, STREAMFLAG);
+
+	/* do not continue streams along edges, these are artifacts */
+	if (edge) {
+	    G_debug(3, "edge");
+	    if (is_swale) {
+		G_debug(2, "edge outlet");
+		/* add outlet point */
+		if (n_outlets >= n_alloc_outlets) {
+		    n_alloc_outlets += stream_node_step;
+		    outlets =
+			(POINT *)G_realloc(outlets,
+						  n_alloc_outlets *
+						  sizeof(POINT));
+		}
+		outlets[n_outlets].r = r;
+		outlets[n_outlets].c = c;
+		n_outlets++;
+		if (af.asp > 0) {
+		    af.asp = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+		    seg_put(&aspflag, (char *)&af, r, c);
+		}
+	    }
+	    continue;
+	}
+
+	if (np_side < 0)
+	    G_fatal_error("np_side < 0");
+	    
+	/* set main drainage direction to A* path if possible */
+	if (mfd_cells > 0 && max_side != np_side) {
+	    if (fabs(wat_nbr[np_side] >= max_acc)) {
+		max_acc = fabs(wat_nbr[np_side]);
+		r_max = dr;
+		c_max = dc;
+		max_side = np_side;
+	    }
+	}
+	if (mfd_cells == 0) {
+	    flat = 0;
+	    r_max = dr;
+	    c_max = dc;
+	    max_side = np_side;
+	}
+
+	/* update aspect */
+	/* r_max == r && c_max == c should not happen */
+	if ((r_max != dr || c_max != dc) && (r_max != r || c_max != c)) {
+	    af.asp = drain[r - r_max + 1][c - c_max + 1];
+	    seg_put(&aspflag, (char *)&af, r, c);
+	}
+
+	/**********************/
+	/*  start new stream  */
+	/**********************/
+
+	/* Montgomery's stream initiation acc * (tan(slope))^mont_exp */
+	/* uses whatever unit is accumulation */
+	if (mont_exp > 0) {
+	    if (r_max == r && c_max == c)
+		G_warning
+		    (_("Can't use Montgomery's method, no stream direction found"));
+	    else {
+		slope = (double)(ele_val - ele_nbr[max_side]) / ele_scale;
+
+		if (max_side > 3)
+		    slope /= diag;
+
+		value *= pow(fabs(slope), mont_exp);
+	    }
+	}
+
+	if (!is_swale && fabs(value) >= threshold && stream_cells < 1 &&
+	    swale_cells < 1 && !flat) {
+	    G_debug(2, "start new stream");
+	    is_swale = ++stream_no;
+	    cseg_put(&stream, &is_swale, r, c);
+	    FLAG_SET(af.flag, STREAMFLAG);
+	    seg_put(&aspflag, (char *)&af, r, c);
+	    /* add stream node */
+	    if (stream_no >= n_alloc_nodes - 1) {
+		n_alloc_nodes += stream_node_step;
+		stream_node =
+		    (struct snode *)G_realloc(stream_node,
+					      n_alloc_nodes *
+					      sizeof(struct snode));
+	    }
+	    stream_node[stream_no].r = r;
+	    stream_node[stream_no].c = c;
+	    stream_node[stream_no].id = stream_no;
+	    stream_node[stream_no].n_trib = 0;
+	    stream_node[stream_no].n_trib_total = 0;
+	    stream_node[stream_no].n_alloc = 0;
+	    stream_node[stream_no].trib = NULL;
+	    n_stream_nodes++;
+
+	    /* debug */
+	    if (n_stream_nodes != stream_no)
+		G_warning(_("BUG: stream_no %d and n_stream_nodes %lld out of sync"),
+			  stream_no, n_stream_nodes);
+	}
+
+	/*********************/
+	/*  continue stream  */
+	/*********************/
+
+	if (is_swale > 0) {
+	    cseg_get(&stream, &is_swale, r, c);
+	    if (r_max == r && c_max == c) {
+		/* can't continue stream, add outlet point
+		 * r_max == r && c_max == c should not happen */
+		G_debug(1, "can't continue stream at r %d c %d", r, c);
+
+		if (n_outlets >= n_alloc_outlets) {
+		    n_alloc_outlets += stream_node_step;
+		    outlets =
+			(POINT *)G_malloc(n_alloc_outlets *
+						 sizeof(POINT));
+		}
+		outlets[n_outlets].r = r;
+		outlets[n_outlets].c = c;
+		n_outlets++;
+	    }
+	    else {
+		continue_stream(is_swale, r_max, c_max, &stream_no);
+	    }
+	}
+    }
+    G_percent(1, 1, 2);
+    if (workedon)
+	G_warning(_("MFD: A * path already processed when setting drainage direction: %lld of %lld cells"),
+		  workedon, n_points);
+
+    G_free(dist_to_nbr);
+    G_free(wat_nbr);
+    G_free(ele_nbr);
+    G_free(flag_nbr);
+
+    G_debug(1, "%lld outlets", n_outlets);
+    G_debug(1, "%lld nodes", n_stream_nodes);
+    G_debug(1, "%d streams", stream_no);
+
+    return 1;
+}

raster/r.stream.extract/thin.c

@@ -0,0 +1,172 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+int thin_seg(int stream_id)
+{
+    int thinned = 0;
+    int r, c, r_nbr, c_nbr, last_r, last_c;
+    CELL curr_stream, no_stream = 0;
+    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+    ASP_FLAG af;
+
+    r = stream_node[stream_id].r;
+    c = stream_node[stream_id].c;
+
+    cseg_get(&stream, &curr_stream, r, c);
+
+    seg_get(&aspflag, (char *)&af, r, c);
+    if (af.asp > 0) {
+	/* get downstream point */
+	last_r = r + asp_r[(int)af.asp];
+	last_c = c + asp_c[(int)af.asp];
+	cseg_get(&stream, &curr_stream, last_r, last_c);
+
+	if (curr_stream != stream_id)
+	    return thinned;
+
+	/* get next downstream point */
+	seg_get(&aspflag, (char *)&af, last_r, last_c);
+	while (af.asp > 0) {
+	    r_nbr = last_r + asp_r[(int)af.asp];
+	    c_nbr = last_c + asp_c[(int)af.asp];
+
+	    if (r_nbr == last_r && c_nbr == last_c)
+		return thinned;
+	    if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
+		return thinned;
+	    cseg_get(&stream, &curr_stream, r_nbr, c_nbr);
+	    if (curr_stream != stream_id)
+		return thinned;
+	    if (abs(r_nbr - r) < 2 && abs(c_nbr - c) < 2) {
+		/* eliminate last point */
+		cseg_put(&stream, &no_stream, last_r, last_c);
+		FLAG_UNSET(af.flag, STREAMFLAG);
+		seg_put(&aspflag, (char *)&af, last_r, last_c);
+		/* update start point */
+		seg_get(&aspflag, (char *)&af, r, c);
+		af.asp = drain[r - r_nbr + 1][c - c_nbr + 1];
+		seg_put(&aspflag, (char *)&af, r, c);
+
+		thinned = 1;
+	    }
+	    else {
+		/* nothing to eliminate, continue from last point */
+		r = last_r;
+		c = last_c;
+	    }
+	    last_r = r_nbr;
+	    last_c = c_nbr;
+	    seg_get(&aspflag, (char *)&af, last_r, last_c);
+	}
+    }
+
+    return thinned;
+}
+
+int thin_streams(void)
+{
+    int i, j, r, c, done;
+    CELL stream_id;
+    int next_node;
+    struct sstack
+    {
+	int stream_id;
+	int next_trib;
+    } *nodestack;
+    int top = 0, stack_step = 1000;
+    int n_trib_total;
+    int n_thinned = 0;
+
+    G_message(_("Thinning stream segments..."));
+
+    nodestack = (struct sstack *)G_malloc(stack_step * sizeof(struct sstack));
+
+    for (i = 0; i < n_outlets; i++) {
+	G_percent(i, n_outlets, 2);
+	r = outlets[i].r;
+	c = outlets[i].c;
+	cseg_get(&stream, &stream_id, r, c);
+
+	if (stream_id == 0)
+	    continue;
+
+	/* add root node to stack */
+	G_debug(2, "add root node");
+	top = 0;
+	nodestack[top].stream_id = stream_id;
+	nodestack[top].next_trib = 0;
+
+	/* depth first post order traversal */
+	G_debug(2, "traverse");
+	while (top >= 0) {
+
+	    done = 1;
+	    stream_id = nodestack[top].stream_id;
+	    G_debug(3, "stream_id %d, top %d", stream_id, top);
+	    if (nodestack[top].next_trib < stream_node[stream_id].n_trib) {
+		/* add to stack */
+		G_debug(3, "get next node");
+		next_node =
+		    stream_node[stream_id].trib[nodestack[top].next_trib];
+		G_debug(3, "add to stack: next %d, trib %d, n trib %d",
+			next_node, nodestack[top].next_trib,
+			stream_node[stream_id].n_trib);
+		nodestack[top].next_trib++;
+		top++;
+		if (top >= stack_step) {
+		    /* need more space */
+		    stack_step += 1000;
+		    nodestack =
+			(struct sstack *)G_realloc(nodestack,
+						   stack_step *
+						   sizeof(struct sstack));
+		}
+
+		nodestack[top].next_trib = 0;
+		nodestack[top].stream_id = next_node;
+		done = 0;
+		G_debug(3, "go further down");
+	    }
+	    if (done) {
+		/* thin stream segment */
+		G_debug(3, "thin stream segment %d", stream_id);
+
+		if (thin_seg(stream_id) == 0)
+		    G_debug(3, "segment %d not thinned", stream_id);
+		else {
+		    G_debug(3, "segment %d thinned", stream_id);
+		    n_thinned++;
+		}
+
+		top--;
+		/* count tributaries */
+		if (top >= 0) {
+		    n_trib_total = 0;
+		    stream_id = nodestack[top].stream_id;
+		    for (j = 0; j < stream_node[stream_id].n_trib; j++) {
+			/* intermediate */
+			if (stream_node[stream_node[stream_id].trib[j]].
+			    n_trib > 0)
+			    n_trib_total +=
+				stream_node[stream_node[stream_id].trib[j]].
+				n_trib_total;
+			/* start */
+			else
+			    n_trib_total++;
+		    }
+		    stream_node[stream_id].n_trib_total = n_trib_total;
+		}
+	    }
+	}
+    }
+    G_percent(n_outlets, n_outlets, 1);	/* finish it */
+
+    G_free(nodestack);
+    
+    G_verbose_message(_("%d of %lld stream segments were thinned"), n_thinned, n_stream_nodes);
+
+    return 1;
+}