r.stream.segment moved from grass-addons to trunk

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

raster/r.stream.segment/Makefile

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

raster/r.stream.segment/io.c

@@ -0,0 +1,564 @@
+#include "io.h"
+/* all in ram functions section */
+
+int ram_create_map(MAP * map, RASTER_MAP_TYPE data_type)
+{
+
+    /* 
+     * allocates 0 filled nrows*ncols map of type void;
+     * map parameters are stored in structure;
+     * map: map to be created;
+     * map type to be created must be CELL, FCELL, DCELL;
+     * */
+
+    int r;
+
+    if (data_type < 0 || data_type > 2)
+	G_fatal_error(_("Unable to create raster map of unrecognised type"));
+
+    map->data_type = data_type;
+    map->map_name = NULL;
+    map->nrows = Rast_window_rows();
+    map->ncols = Rast_window_cols();
+    map->data_size = Rast_cell_size(data_type);
+
+    /* preparing internal map */
+    switch (map->data_type) {
+    case CELL_TYPE:
+	map->map = G_calloc(map->nrows, sizeof(CELL *));
+	break;
+
+    case FCELL_TYPE:
+	map->map = G_calloc(map->nrows, sizeof(FCELL *));
+	break;
+
+    case DCELL_TYPE:
+	map->map = G_calloc(map->nrows, sizeof(DCELL *));
+	break;
+    }
+
+    for (r = 0; r < map->nrows; ++r)
+	(map->map)[r] = G_calloc(map->ncols, map->data_size);
+
+    return 0;
+}
+
+int ram_read_map(MAP * map, char *input_map_name, int check_res,
+		 RASTER_MAP_TYPE check_data_type)
+{
+    /*
+     * Function read external map and put it in MAP structure (created with create_map)
+     * map: map to be read can be of any data type, read map is converted to target map if necessary.
+     * input_map_name: name of the map to be read;
+     * map pointer to map structure (created with create_map);
+     * check_res: [1]: check res correspondence between region and map [0 no check];
+     * check_data_type [CELL, FCELL, DCELL] check if reading map is of particular type, [-1] no check;
+     */
+
+    int r, c;
+    char *mapset;
+    struct Cell_head cellhd, this_window;
+    char *maptypes[] = { "CELL", "FCELL", "DCELL" };
+    int input_map_fd;
+    RASTER_MAP_TYPE input_data_type;
+    size_t input_data_size;
+    void *input_buffer = NULL;
+    void *input_pointer;
+
+    /* checking if map exist */
+    mapset = (char *)G_find_raster2(input_map_name, "");
+    if (mapset == NULL)
+	G_fatal_error(_("Raster map <%s> not found"), input_map_name);
+
+    /* checking if region and input are the same */
+    G_get_window(&this_window);
+    Rast_get_cellhd(input_map_name, mapset, &cellhd);
+    if (check_res)
+	if (this_window.ew_res != cellhd.ew_res ||
+	    this_window.ns_res != cellhd.ns_res)
+          G_fatal_error(_("Region resolution and raster map <%s> resolution differs. "
+                          "Run 'g.region rast=%s' to set proper region resolution."),
+                        input_map_name, input_map_name);
+
+    /* checking if input map is of required type */
+    if (check_data_type != map->data_type)
+	G_debug(1,
+		"ram_open:required map type and internal map type differs: conversion forced!");
+    input_data_type = Rast_map_type(input_map_name, mapset);
+    if (check_data_type != -1)
+	if (input_data_type != check_data_type)
+	    G_fatal_error(_("Raster map <%s> is not of type '%s'"),
+			  input_map_name, maptypes[check_data_type]);
+
+    input_map_fd = Rast_open_old(input_map_name, mapset);
+    input_data_size = Rast_cell_size(input_data_type);
+
+    {				/* reading range */
+	struct Range map_range;
+	struct FPRange map_fp_range;
+	int min, max;
+
+	if (input_data_type == CELL_TYPE) {
+	    Rast_init_range(&map_range);
+	    Rast_read_range(input_map_name, mapset, &map_range);
+	    Rast_get_range_min_max(&map_range, &min, &max);
+	    map->min = (double)min;
+	    map->max = (double)max;
+	}
+	else {
+	    Rast_init_fp_range(&map_fp_range);
+	    Rast_read_fp_range(input_map_name, mapset, &map_fp_range);
+	    Rast_get_fp_range_min_max(&map_fp_range, &(map->min),
+				      &(map->max));
+	}
+    }
+    /* end opening and checking */
+
+    input_buffer = Rast_allocate_buf(input_data_type);
+
+    /* start reading */
+    G_message(_("Reading raster map <%s>..."), input_map_name);
+
+    for (r = 0; r < map->nrows; ++r) {
+	G_percent(r, map->nrows, 2);
+
+	Rast_get_row(input_map_fd, input_buffer, r, input_data_type);
+	input_pointer = input_buffer;
+
+	for (c = 0; c < map->ncols; ++c)
+	    if (!Rast_is_null_value
+		(input_pointer + c * input_data_size, input_data_type))
+		switch (map->data_type) {
+		case CELL_TYPE:
+		    ((CELL **) map->map)[r][c] =
+			Rast_get_c_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		case FCELL_TYPE:
+		    ((FCELL **) map->map)[r][c] =
+			Rast_get_f_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		case DCELL_TYPE:
+		    ((DCELL **) map->map)[r][c] =
+			Rast_get_d_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		default:
+		    G_fatal_error(_("Wrong internal data type"));
+		    break;
+		}
+    }				/*end for r */
+
+    G_free(input_buffer);
+    G_percent(r, map->nrows, 2);
+    Rast_close(input_map_fd);
+    return 0;
+}				/* end create floating point map */
+
+int ram_reset_map(MAP * map, int value)
+{
+    /*
+     * set all cells in the map to value
+     */
+    int r;
+
+    for (r = 0; r < map->nrows; ++r)
+	memset((map->map)[r], value, map->ncols * map->data_size);
+    return 0;
+}
+
+int ram_write_map(MAP * map, char *output_map_name,
+		  RASTER_MAP_TYPE output_data_type, int convert_to_null,
+		  double value)
+{
+    /* 
+     * write map to disk with output_map_name and output_data_type [CELL, FCELL, DCELL];
+     * if output_data_type = -1 than internal map type is used for output;
+     * if output map != -1 and types differ data_type, conversion is forced
+     * convert to null: check if convert to null a particular value in dataset;
+     */
+
+    int r, c;
+    int output_fd = 0;
+    struct History history;
+    void *row;
+
+    /* check for output format */
+    if (output_data_type == -1)
+	output_data_type = map->data_type;
+
+    if (output_data_type != map->data_type)
+	G_debug(1,
+		"ram_write:required map type and internal map type differs: conversion forced!");
+
+    G_message(_("Writing raster map <%s>..."), output_map_name);
+    output_fd = Rast_open_new(output_map_name, output_data_type);
+
+    /* writing */
+    for (r = 0; r < map->nrows; ++r) {
+	G_percent(r, map->nrows, 2);
+
+	if (convert_to_null) {
+	    row = map->map[r];
+	    switch (map->data_type) {
+	    case CELL_TYPE:
+		for (c = 0; c < map->ncols; ++c)
+		    if (((CELL *) row)[c] == (CELL) value)
+			Rast_set_c_null_value(row + c * (map->data_size), 1);
+		break;
+	    case FCELL_TYPE:
+		for (c = 0; c < map->ncols; ++c)
+		    if (((FCELL *) row)[c] == (FCELL) value)
+			Rast_set_f_null_value(row + c * (map->data_size), 1);
+		break;
+	    case DCELL_TYPE:
+		for (c = 0; c < map->ncols; ++c)
+		    if (((DCELL *) row)[c] == (DCELL) value)
+			Rast_set_d_null_value(row + c * (map->data_size), 1);
+		break;
+	    default:
+		G_debug(1, "ram_null:Cannot convert to null at: %d %d", r, c);
+	    }
+	}
+
+	Rast_put_row(output_fd, (map->map)[r], output_data_type);
+    }
+    G_percent(r, map->nrows, 2);
+    Rast_close(output_fd);
+    Rast_short_history(output_map_name, "raster", &history);
+    Rast_command_history(&history);
+    Rast_write_history(output_map_name, &history);
+    /* G_message(_("<%s> Done"), output_map_name); */
+    return 0;
+}
+
+int ram_release_map(MAP *map)
+{
+    /* 
+     * free memory allocated for map, set pointer to null;
+     */
+    int r;
+
+    for (r = 0; r < map->nrows; ++r)
+	G_free((map->map)[r]);
+    G_free(map->map);
+    map = NULL;
+    return 0;
+}
+
+
+/* memory swap functions section */
+
+
+int seg_create_map(SEG * seg, int srows, int scols, int number_of_segs,
+		   RASTER_MAP_TYPE data_type)
+{
+    /* create segment  and returns pointer to it;
+     * seg must be declared first;
+     * parameters are stored in structure;
+     * seg: segment to be created;
+     * srows, scols segment size
+     * number of segs max number of segs stored in memory
+     * data_type to be created must be CELL, FCELL, DCELL;
+     */
+
+    char *filename;
+    int fd;
+    int local_number_of_segs;
+
+    seg->fd = -1;
+    seg->filename = NULL;
+    seg->map_name = NULL;
+    seg->mapset = NULL;
+    seg->data_type = data_type;
+    seg->nrows = Rast_window_rows();
+    seg->ncols = Rast_window_cols();
+
+    local_number_of_segs =
+	(seg->nrows / srows + 1) * (seg->ncols / scols + 1);
+    number_of_segs =
+	(number_of_segs >
+	 local_number_of_segs) ? local_number_of_segs : number_of_segs;
+
+    G_debug(3, "seg_creat:number of segments %d", number_of_segs);
+
+    switch (seg->data_type) {
+    case CELL_TYPE:
+	seg->data_size = sizeof(CELL);
+	break;
+    case FCELL_TYPE:
+	seg->data_size = sizeof(FCELL);
+	break;
+    case DCELL_TYPE:
+	seg->data_size = sizeof(DCELL);
+	break;
+    default:
+	G_fatal_error(_("Unrecognisable data type"));
+    }
+
+    filename = G_tempfile();
+    fd = creat(filename, 0666);
+
+    if (0 >
+	segment_format(fd, seg->nrows, seg->ncols, srows, scols,
+		       seg->data_size)) {
+	close(fd);
+	unlink(filename);
+	G_fatal_error(_("Unable to format segment"));
+    }
+
+    close(fd);
+    if (0 > (fd = open(filename, 2))) {
+	unlink(filename);
+	G_fatal_error(_("Unable to re-open file '%s'"), filename);
+    }
+
+    if (0 > (fd = segment_init(&(seg->seg), fd, number_of_segs))) {
+	unlink(filename);
+	G_fatal_error(_("Unable to init segment file or out of memory"));
+    }
+
+    seg->filename = G_store(filename);
+    seg->fd = fd;
+    return 0;
+}
+
+int seg_read_map(SEG * seg, char *input_map_name, int check_res,
+		 RASTER_MAP_TYPE check_data_type)
+{
+
+    /*
+     * Function read external map and put it in SEG structure (created with seg_create_map)
+     * map to be read can be of any data type, read map is converted if necessary.
+     * input_map_name: name of the map to be read;
+     * seg: pointer to map structure (created with create_map);
+     * check_res: [1]: check res correspondence between region and map [0 no check];
+     * check_data_type [CELL, FCELL, DCELL] check if reading map is of particular type, [-1] no check;
+     */
+
+    int input_fd;
+    int r, c;
+    char *mapset;
+    struct Cell_head cellhd, this_window;
+    char *maptypes[] = { "CELL", "FCELL", "DCELL" };
+    RASTER_MAP_TYPE input_data_type;
+    size_t input_data_size;
+    void *input_buffer = NULL;
+    void *target_buffer = NULL;
+    void *input_pointer = NULL;
+
+    /* checking if map exist */
+    mapset = (char *)G_find_raster2(input_map_name, "");
+    if (mapset == NULL)
+	G_fatal_error(_("Raster map <%s> not found"),
+		      input_map_name);
+    seg->mapset = mapset;
+
+    /* checking if region and input are the same */
+    G_get_window(&this_window);
+    Rast_get_cellhd(input_map_name, mapset, &cellhd);
+
+    /* check resolution equal any integer check;  equal 0 no check */
+    if (check_res)
+	if (this_window.ew_res != cellhd.ew_res ||
+	    this_window.ns_res != cellhd.ns_res)
+          G_fatal_error(_("Region resolution and raster map <%s> resolution differs. "
+                          "Run 'g.region rast=%s' to set proper region resolution."),
+			  input_map_name, input_map_name);
+
+    if (check_data_type != seg->data_type)
+	G_debug(1,
+		"ram_open:required map type and internal map type differs: conversion forced!");
+    input_data_type = Rast_map_type(input_map_name, mapset);
+    if (check_data_type != -1)
+	if (input_data_type != check_data_type)
+	    G_fatal_error(_("Raster map <%s> is not of type '%s'"),
+			  input_map_name, maptypes[check_data_type]);
+
+    input_fd = Rast_open_old(input_map_name, mapset);
+    input_data_size = Rast_cell_size(input_data_type);
+
+    {				/* reading range */
+	struct Range map_range;
+	struct FPRange map_fp_range;
+	int min, max;
+
+	if (input_data_type == CELL_TYPE) {
+	    Rast_init_range(&map_range);
+	    Rast_read_range(input_map_name, mapset, &map_range);
+	    Rast_get_range_min_max(&map_range, &min, &max);
+	    seg->min = (double)min;
+	    seg->max = (double)max;
+	}
+	else {
+	    Rast_init_fp_range(&map_fp_range);
+	    Rast_read_fp_range(input_map_name, mapset, &map_fp_range);
+	    Rast_get_fp_range_min_max(&map_fp_range, &(seg->min),
+				      &(seg->max));
+	}
+    }
+
+    /* end opening and checking */
+
+    G_message(_("Reading raster map <%s>..."), input_map_name);
+    input_buffer = Rast_allocate_buf(input_data_type);
+
+    target_buffer = Rast_allocate_buf(seg->data_type);
+
+    for (r = 0; r < seg->nrows; ++r) {
+	G_percent(r, seg->nrows, 2);
+	Rast_get_row(input_fd, input_buffer, r, input_data_type);
+	input_pointer = input_buffer;
+	memset(target_buffer, 0, seg->ncols * seg->data_size);
+
+	for (c = 0; c < seg->ncols; ++c)
+	    if (!Rast_is_null_value
+		(input_pointer + c * input_data_size, input_data_type)) {
+		switch (seg->data_type) {
+		case CELL_TYPE:
+		    ((CELL *) target_buffer)[c] =
+			Rast_get_c_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		case FCELL_TYPE:
+		    ((FCELL *) target_buffer)[c] =
+			Rast_get_f_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		case DCELL_TYPE:
+		    ((DCELL *) target_buffer)[c] =
+			Rast_get_d_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		default:
+		    G_fatal_error(_("Wrong internal data type"));
+		    break;
+		}
+	    }
+
+	if (0 > segment_put_row(&(seg->seg), target_buffer, r)) {
+	    G_free(input_buffer);
+	    G_free(target_buffer);
+	    Rast_close(input_fd);
+	    G_fatal_error(_("seg_read: Cannot segment put row %d for map %s"),
+			  r, input_map_name);
+	}
+    }				/* end for row */
+
+    G_percent(r, seg->nrows, 2);
+    Rast_close(input_fd);
+    G_free(input_buffer);
+    G_free(target_buffer);
+
+    seg->map_name = G_store(input_map_name);
+    seg->mapset = G_store(mapset);
+
+    return 0;
+}
+
+int seg_reset_map(SEG * seg, int value)
+{
+    /*
+     * set all cells in the map to value
+     */
+    int r, c;
+
+    for (r = 0; r < seg->nrows; ++r)
+	for (c = 0; c < seg->ncols; ++c)
+	    segment_put(&(seg->seg), &value, r, c);
+
+    return 0;
+}
+
+int seg_write_map(SEG * seg, char *output_map_name,
+		  RASTER_MAP_TYPE output_data_type, int convert_to_null,
+		  double value)
+{
+    /* 
+     * write seg to disk with output_map_name and output_data_type [CELL, FCELL, DCELL];
+     * if output_data_type = -1 than internal map type is used for output;
+     * if output map != -1 and types differ data_type, conversion is forced
+     * convert to null: check if convert to null a particular value in dataset;
+     */
+    int output_fd;
+    int r, c;
+    void *output_buffer;
+    void *row;
+    struct History history;
+
+    /* check for output format */
+    if (output_data_type == -1)
+	output_data_type = seg->data_type;
+
+    if (output_data_type != seg->data_type)
+	G_debug(1,
+		"ram_write:required map type and internal map type differs: conversion forced!");
+
+    G_message(_("Writing raster map <%s>..."), output_map_name);
+    output_fd = Rast_open_new(output_map_name, output_data_type);
+    output_buffer = Rast_allocate_buf(output_data_type);
+    segment_flush(&(seg->seg));
+
+    /* writing */
+    for (r = 0; r < seg->nrows; ++r) {
+
+	G_percent(r, seg->nrows, 2);
+	if (0 > segment_get_row(&(seg->seg), output_buffer, r))
+	    G_warning(_("Unable to segment read row %d for raster map <%s>"),
+		      r, output_map_name);
+
+	if (convert_to_null) {
+
+	    row = output_buffer;
+	    switch (seg->data_type) {
+	    case CELL_TYPE:
+		for (c = 0; c < seg->ncols; ++c)
+		    if (((CELL *) output_buffer)[c] == (CELL) value)
+			Rast_set_c_null_value(row + c * (seg->data_size), 1);
+		break;
+	    case FCELL_TYPE:
+		for (c = 0; c < seg->ncols; ++c)
+		    if (((FCELL *) output_buffer)[c] == (FCELL) value)
+			Rast_set_f_null_value(row + c * (seg->data_size), 1);
+		break;
+	    case DCELL_TYPE:
+		for (c = 0; c < seg->ncols; ++c)
+		    if (((DCELL *) output_buffer)[c] == (DCELL) value)
+			Rast_set_d_null_value(row + c * (seg->data_size), 1);
+		break;
+	    default:
+		G_warning(_("Unable to convert to NULL at: %d %d"), r,
+			  c);
+	    }
+	}
+	Rast_put_row(output_fd, output_buffer, output_data_type);
+    }
+
+    G_percent(r, seg->nrows, 2);
+    G_free(output_buffer);
+    Rast_close(output_fd);
+    Rast_short_history(output_map_name, "raster", &history);
+    Rast_command_history(&history);
+    Rast_write_history(output_map_name, &history);
+    /* G_message(_("%s Done"), output_map_name); */
+
+    return 0;
+}
+
+int seg_release_map(SEG * seg)
+{
+    /* 
+     * release segment close files, set pointers to null;
+     */
+    segment_release(&(seg->seg));
+    close(seg->fd);
+    unlink(seg->filename);
+
+    if (seg->map_name)
+	G_free(seg->map_name);
+    if (seg->mapset)
+	G_free(seg->mapset);
+
+    return 0;
+}

raster/r.stream.segment/io.h

@@ -0,0 +1,58 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <string.h>
+#include <math.h>
+#include <grass/glocale.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/segment.h>
+
+#define NOT_IN_REGION(x) (r + nextr[(x)] < 0 || r + nextr[(x)] > (nrows - 1) || \
+                          c + nextc[(x)] < 0 || c + nextc[(x)] > (ncols - 1))
+#define NR(x) (r + nextr[(x)])
+#define NC(x) (c + nextc[(x)])
+#define INDEX(r,c) ((r) * ncols + (c))
+#define DIAG(x) (((x) + 4) > 8 ? ((x) - 4) : ((x) + 4))
+
+#define SROWS 256
+#define SCOLS 256
+
+typedef struct {
+	void **map; /* matrix of data */
+	double min, max; /* data range : may require casting */
+	int nrows, ncols;
+	char *map_name; /* map name, unused */
+	RASTER_MAP_TYPE data_type; /* type of data */
+	size_t data_size; /* type of data */
+} MAP;
+
+typedef struct {
+	SEGMENT seg;		/* segmented data store */
+	int fd;					/* segment temporary file name descriptor */
+	char *filename; /* segment temporary file name */
+	char *map_name; /* map name converted to segment */
+	char *mapset;
+	int nrows, ncols; /* store nrows and rcols */
+	RASTER_MAP_TYPE data_type; /* data type of the map */
+	size_t data_size; /* size of cell returned by sizeof */
+	double min, max; /* data range */
+} SEG;
+
+
+/* all in ram functions */
+int ram_create_map(MAP *, RASTER_MAP_TYPE);
+int ram_read_map(MAP *, char *, int, RASTER_MAP_TYPE);
+int ram_reset_map(MAP *, int);
+int ram_write_map(MAP *, char *, RASTER_MAP_TYPE, int, double);
+int ram_release_map(MAP *);
+int ram_destory_map(MAP *);
+
+/* memory swap functions */
+int seg_create_map(SEG *, int, int, int, RASTER_MAP_TYPE);
+int seg_read_map(SEG *, char *, int, RASTER_MAP_TYPE);
+int seg_reset_map (SEG *, int);
+int seg_write_map(SEG *, char *, RASTER_MAP_TYPE, int, double);
+int seg_release_map(SEG *);

raster/r.stream.segment/local_proto.h

@@ -0,0 +1,30 @@
+#include "io.h"
+#include "local_vars.h"
+
+int free_attributes(int number_of_streams);
+int convert_border_dir(int r, int c, int dir);
+
+int create_sectors(STREAM *cur_stream, int seg_length, int seg_skip,
+		   double seg_treshold);
+int calc_tangents(STREAM *cur_stream, int seg_length, int seg_skip,
+		  int number_streams);
+
+int create_sector_vector(char *out_vector, int number_of_streams, int radians);
+int create_segment_vector(char *out_vector, int number_of_streams,
+			  int radians);
+
+
+int ram_build_streamlines(CELL **streams, CELL **dirs, FCELL **elevation,
+			  int number_of_streams);
+int ram_fill_streams(CELL **unique_streams, int number_of_streams);
+int ram_find_contributing_cell(int r, int c, CELL **dirs, FCELL **elevation);
+int ram_identify_next_stream(CELL **streams, int number_of_streams);
+int ram_number_of_streams(CELL **streams, CELL **dirs, int *ordered);
+
+int seg_build_streamlines(SEGMENT *streams, SEGMENT *dirs,
+			  SEGMENT *elevation, int number_of_streams);
+int seg_fill_streams(SEGMENT *unique_streams, int number_of_streams);
+int seg_find_contributing_cell(int r, int c, SEGMENT *dirs,
+			       SEGMENT *elevation);
+int seg_identify_next_stream(SEGMENT *streams, int number_of_streams);
+int seg_number_of_streams(SEGMENT *streams, SEGMENT *dirs, int *ordered);

raster/r.stream.segment/local_vars.h

@@ -0,0 +1,71 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/glocale.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+
+#define SQRT2 1.414214
+	
+typedef struct {
+	int stream;
+	int next_stream;
+	int number_of_cells;
+	int order;
+	unsigned long int * points;
+	float * elevation;
+	double * distance;
+	unsigned long int init;
+	unsigned long int outlet; /* outlet is cell from next stream */
+	int last_cell_dir; /* to add outlet to vector */
+	float direction;
+	float length;
+	float stright;
+	float drop;
+	float tangent;
+	float continuation;
+	int number_of_sectors;
+	int* sector_breakpoints; /* index of breakpoints in *points vector */
+	int* sector_cats;
+	float* sector_directions;
+	float* sector_strights;
+	double* sector_lengths;
+	float* sector_drops; /* gradient calculated at the end */
+} STREAM;	
+
+typedef struct {
+	float long_dir_diff;
+	float short_dir_diff;
+	int long_break;
+	int decision;
+} DIRCELLS;
+
+#ifdef MAIN
+#	define GLOBAL
+#else
+#	define GLOBAL extern
+#endif
+
+#ifdef MAIN
+#	define GLOBAL
+#else
+#	define GLOBAL extern
+#endif
+
+#ifndef PI
+ #define PI (4*atan(1))
+#endif
+
+#define DEG2RAD(d) ((d)*PI/180)
+#define RAD2DEG(r) ((r)*180/PI)
+
+GLOBAL int nextr[9];
+GLOBAL int nextc[9];
+
+GLOBAL int nrows, ncols; 
+GLOBAL STREAM* stream_attributes;
+
+GLOBAL struct Cell_head window;

raster/r.stream.segment/main.c

@@ -0,0 +1,291 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.stream.segment
+ * AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl
+ *							 
+ * PURPOSE:	 Calculate geometrical attributes for segments of current order, 
+ * 			 divide segments on near straight line portions and 
+ * 			 segment orientation and angles between streams and its
+ *           tributaries. For stream direction it use algorithm to divide
+ *           particular streams of the same order into near-straight line
+ *           portions.
+ * 				
+ *							
+ *
+ * COPYRIGHT:		 (C) 2002,2010-2014 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.
+ *
+ *****************************************************************************/
+#define MAIN
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+
+int main(int argc, char *argv[])
+{
+
+    struct GModule *module;	/* GRASS module for parsing arguments */
+    struct Option *in_dir_opt,	/* options */
+     *in_stm_opt,
+	*in_elev_opt,
+	*out_segment_opt,
+	*out_sector_opt,
+	*opt_length, *opt_skip, *opt_threshold, *opt_swapsize;
+
+    struct Flag *flag_radians, *flag_segmentation;	/* segmentation library */
+
+    int i;
+    int seg_length, seg_skip;
+    int radians, segmentation;	/* flags */
+    double seg_treshold;
+    int number_of_streams, ordered;
+
+    /* initialize GIS environment */
+    G_gisinit(argv[0]);		/* reads grass env, stores program name to G_program_name() */
+
+    /* initialize module */
+    module = G_define_module();
+    module->description =
+	_("Divides network into near straight-line segments and calculate its order.");
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("hydrology"));
+    G_add_keyword(_("stream network"));
+    G_add_keyword(_("stream divide"));
+
+    in_stm_opt = G_define_standard_option(G_OPT_R_INPUT);
+    in_stm_opt->key = "stream_raster";
+    in_stm_opt->description = _("Name of input streams mask raster map");
+
+    in_dir_opt = G_define_standard_option(G_OPT_R_INPUT);
+    in_dir_opt->key = "direction";
+    in_dir_opt->description = _("Name for input raster map with flow direction");
+
+    in_elev_opt = G_define_standard_option(G_OPT_R_ELEV);
+
+    out_segment_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+    out_segment_opt->key = "segments";
+    out_segment_opt->description =
+	_("Name for output vector map to write segment attributes");
+
+    out_sector_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+    out_sector_opt->key = "sectors";
+    out_sector_opt->description =
+	_("Name for output vector map to write segment attributes");
+
+    opt_length = G_define_option();
+    opt_length->key = "length";
+    opt_length->label = _("Search length to calculate direction");
+    opt_length->description = _("Must be > 0");
+    opt_length->answer = "15";
+    opt_length->type = TYPE_INTEGER;
+
+    opt_skip = G_define_option();
+    opt_skip->key = "skip";
+    opt_skip->label = _("Skip segments shorter than");
+    opt_skip->description = _("Must be >= 0");
+    opt_skip->answer = "5";
+    opt_skip->type = TYPE_INTEGER;
+
+    opt_threshold = G_define_option();
+    opt_threshold->key = "threshold";
+    opt_threshold->label = _("Max angle (degrees) between stream segments");
+    opt_threshold->description = _("Must be > 0");
+    opt_threshold->answer = "160";
+    opt_threshold->type = TYPE_DOUBLE;
+
+    opt_swapsize = G_define_option();
+    opt_swapsize->key = "memory";
+    opt_swapsize->type = TYPE_INTEGER;
+    opt_swapsize->answer = "300";
+    opt_swapsize->description = _("Max memory used in memory swap mode (MB)");
+    opt_swapsize->guisection = _("Memory setings");
+
+    flag_radians = G_define_flag();
+    flag_radians->key = 'r';
+    flag_radians->description =
+	_("Output angles in radians (default: degrees)");
+
+    flag_segmentation = G_define_flag();
+    flag_segmentation->key = 'm';
+    flag_segmentation->description = _("Use memory swap (operation is slow)");
+    flag_segmentation->guisection = _("Memory setings");
+
+    if (G_parser(argc, argv))	/* parser */
+	exit(EXIT_FAILURE);
+
+    seg_length = atoi(opt_length->answer);
+    seg_treshold = atof(opt_threshold->answer);
+    seg_skip = atoi(opt_skip->answer);
+    radians = (flag_radians->answer != 0);
+    segmentation = (flag_segmentation->answer != 0);
+
+    if (seg_length <= 0)
+        G_fatal_error(_("Search's length must be > 0"));
+    if (seg_treshold < 0 || seg_treshold > 180)
+	G_fatal_error(_("Threshold must be between 0 and 180"));
+    if (seg_skip < 0)
+	G_fatal_error(_("Segment's length must be >= 0"));
+
+    seg_treshold = DEG2RAD(seg_treshold);
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+    Rast_get_window(&window);
+    G_begin_distance_calculations();
+
+
+
+    if (!segmentation) {
+	MAP map_dirs, map_streams, map_elevation, map_unique_streams;
+	CELL **streams, **dirs, **unique_streams = NULL;
+	FCELL **elevation;
+
+	G_message(_("All in RAM calculation..."));
+	ram_create_map(&map_streams, CELL_TYPE);
+	ram_read_map(&map_streams, in_stm_opt->answer, 1, CELL_TYPE);
+	ram_create_map(&map_dirs, CELL_TYPE);
+	ram_read_map(&map_dirs, in_dir_opt->answer, 1, CELL_TYPE);
+	ram_create_map(&map_elevation, FCELL_TYPE);
+	ram_read_map(&map_elevation, in_elev_opt->answer, 0, -1);
+
+	streams = (CELL **) map_streams.map;
+	dirs = (CELL **) map_dirs.map;
+	elevation = (FCELL **) map_elevation.map;
+
+	number_of_streams =
+	    ram_number_of_streams(streams, dirs, &ordered) + 1;
+	ram_build_streamlines(streams, dirs, elevation, number_of_streams);
+
+	if (ordered) {
+	    ram_create_map(&map_unique_streams, CELL_TYPE);
+	    unique_streams = (CELL **) map_unique_streams.map;
+	    ram_fill_streams(unique_streams, number_of_streams);
+	    ram_identify_next_stream(unique_streams, number_of_streams);
+	    ram_release_map(&map_unique_streams);
+	}
+	else
+	    ram_identify_next_stream(streams, number_of_streams);
+
+	ram_release_map(&map_streams);
+	ram_release_map(&map_dirs);
+	ram_release_map(&map_elevation);
+    }
+
+
+    if (segmentation) {
+	SEG map_dirs, map_streams, map_elevation, map_unique_streams;
+	SEGMENT *streams, *dirs, *unique_streams = NULL;
+	SEGMENT *elevation;
+	int number_of_segs;
+
+        G_message(_("Memory swap calculation (may take some time)..."));
+	number_of_segs = (int)atof(opt_swapsize->answer);
+	number_of_segs = number_of_segs < 32 ? (int)(32 / 0.18) : number_of_segs / 0.18;
+
+	seg_create_map(&map_streams, SROWS, SCOLS, number_of_segs, CELL_TYPE);
+	seg_read_map(&map_streams, in_stm_opt->answer, 1, CELL_TYPE);
+	seg_create_map(&map_dirs, SROWS, SCOLS, number_of_segs, CELL_TYPE);
+	seg_read_map(&map_dirs, in_dir_opt->answer, 1, CELL_TYPE);
+	seg_create_map(&map_elevation, SROWS, SCOLS, number_of_segs,
+		       FCELL_TYPE);
+	seg_read_map(&map_elevation, in_elev_opt->answer, 0, -1);
+
+	streams = &map_streams.seg;
+	dirs = &map_dirs.seg;
+	elevation = &map_elevation.seg;
+
+	number_of_streams =
+	    seg_number_of_streams(streams, dirs, &ordered) + 1;
+	seg_build_streamlines(streams, dirs, elevation, number_of_streams);
+
+	if (ordered) {
+	    seg_create_map(&map_unique_streams, SROWS, SCOLS, number_of_segs,
+			   CELL_TYPE);
+	    unique_streams = &map_unique_streams.seg;
+	    seg_fill_streams(unique_streams, number_of_streams);
+	    seg_identify_next_stream(unique_streams, number_of_streams);
+	    seg_release_map(&map_unique_streams);
+	}
+	else
+	    seg_identify_next_stream(streams, number_of_streams);
+
+	seg_release_map(&map_streams);
+	seg_release_map(&map_dirs);
+	seg_release_map(&map_elevation);
+    }
+
+
+    for (i = 1; i < number_of_streams; ++i)
+	G_message("%d %d %d", stream_attributes[i].stream,
+		  stream_attributes[i].next_stream,
+		  stream_attributes[i].last_cell_dir);
+
+
+    /*
+       for(i=1;i<number_of_streams;++i)
+       printf("STREAM %d NEXT_STREAM %d POINT %d \n",
+       stream_attributes[i].stream,
+       stream_attributes[i].next_stream,
+       stream_attributes[i].outlet);
+
+     */
+    G_message(_("Creating sectors and calculating attributes..."));
+
+    for (i = 1; i < number_of_streams; ++i) {
+	create_sectors(&stream_attributes[i], seg_length, seg_skip,
+		       seg_treshold);
+	calc_tangents(&stream_attributes[i], seg_length, seg_skip,
+		      number_of_streams);
+    }
+
+
+    /*
+
+
+       for(j=1;j<number_of_streams;++j)
+       G_message("STREAM %d   ncells %d len %f str %f sin %f",
+       stream_attributes[j].stream,
+       stream_attributes[j].number_of_cells,
+       stream_attributes[j].length,
+       stream_attributes[j].stright,
+       stream_attributes[j].length/stream_attributes[j].stright);
+       G_message("%d",j);
+
+
+       for(j=1;j<number_of_streams;++j)
+       printf("STREAM %d   max %f min %f drop %f\n",
+       stream_attributes[j].stream,
+       stream_attributes[j].elevation[1],
+       stream_attributes[j].elevation[stream_attributes[j].number_of_cells-1],
+       stream_attributes[j].drop);
+
+       for(j=1;j<number_of_streams;++j)
+       for (i = 0; i < stream_attributes[j].number_of_sectors; ++i)
+       printf("%d  cat %d |BRAEK %d | dir %.2f | len %.2f | drop %.3f  \n" ,j,
+       stream_attributes[j].sector_cats[i],
+       stream_attributes[j].sector_breakpoints[i],
+       stream_attributes[j].sector_directions[i],
+       stream_attributes[j].sector_lengths[i],
+       stream_attributes[j].sector_drops[i]);
+
+       for(j=1;j<number_of_streams;++j) {
+       printf("        %d %d \n" ,j,stream_attributes[j].number_of_cells);
+       for (i = 0; i <= stream_attributes[j].number_of_cells; ++i)
+       printf("%d %f \n" ,i,stream_attributes[j].elevation[i]);
+       }
+     */
+
+    create_segment_vector(out_segment_opt->answer, number_of_streams,
+			  radians);
+    create_sector_vector(out_sector_opt->answer, number_of_streams, radians);
+
+    free_attributes(number_of_streams);
+    G_message("Done");
+    exit(EXIT_SUCCESS);
+}

raster/r.stream.segment/r.stream.segment.html

@@ -0,0 +1,214 @@
+<h2>DESCRIPTION</h2>
+
+<h2>OPTIONS</h2>
+<dl>
+<dt><b>-r</b></dt>
+<dd>Directions and azimut output in radians. Default is degrees.</dd>
+
+<dt><b>-m</b></dt>
+<dd>Only for very large data sets. Use segment library to optimize memory
+consumption during analysis</dd>
+
+<dt><b>stream_rast</b></dt>
+<dd>Stream network: name of input stream map. Streams shall be ordered according
+one of the r.stream.order ordering system as well as unordered (with original
+stream identifiers)  Because streams network produced by r.watershed and
+r.stream.extract may slightly differ in detail, it is required to use both stream
+and direction map produced by the same module. Stream background shall have NULL
+value or zero value. Background values of NULL are by default produced by
+r.watershed and r.stream.extract. If not 0 or NULL use <a
+href="r.mapcalc.html">r.mapcalc</a> to set background values to null.  
+</dd>
+
+<dt><b>direction</b></dt>
+<dd>Flow direction: name of input direction map produced by r.watershed or
+r.stream.extract. If r.stream.extract output map is used, it only has non-NULL
+values in places where streams occur. NULL (nodata) cells are ignored, zero and
+negative values are valid direction data if they vary from -8 to 8 (CCW from
+East in steps of 45 degrees). Direction map shall be of type CELL values. Region
+resolution and map resolution must be the same. 
+Also <em>stream_rast</em> network map must have the same resolution. It is checked by
+default. If resolutions differ the module informs about it and stops. Region
+boundary and maps boundary may be differ but it may lead to unexpected
+results.</dd>
+
+<dt><b>elevation</b></dt>
+<dd>Elevation: name of input elevation map. Map can be of type CELL, FCELL or
+DCELL. It is not restricted to resolution of region settings as stream_rast and
+direction.</dd>
+
+<dt><b>length</b></dt>
+<dd>Integer values indicating the search length (in cells) to determine straight
+line. The longest length parameter the module treats more tolerant local stream
+undulation and inequalities. Default value of 15 is suitable for 30 meters
+DEMs. More detail DEMs may requre longer length.</dd>
+
+<dt><b>skip</b></dt>
+<dd>Integer values indicating the length (in cells) local short segment to skip
+and join them to the longer neighbour. The shortest length parameter the more
+short segments will be produced by the module due to undulation and
+inequalities. Default value of 5 is suitable for 30 meters DEMS. More details
+DEMS may require longer length.</dd>
+
+<dt><b>threshold</b></dt>
+<dd>real value indicates the internal angle between upstream and downstream
+direction to treat actual cell as lying on the straight line. Greater values (up to
+180 degrees) produce more segments. Lower values produce less segments.
+Values below 90 in most cases will not produce any additional segments to those
+resulting from ordering.</dl>
+
+<dl>
+<h2>OUTPUTS</h2>
+<p>The module produces two vector maps: one representing original segments
+(where a segment is a streamline where its order remains unchanged) and the second
+divided into near straight line sectors resulting form segmentation process.
+Most of the segment and sectors attributes are the same as in r.stream.order vector
+output.</p>
+
+<dl>
+<dt><b>segments</b></dt>
+<dd>
+Vector map where every segment has its own category and following attributes:
+<ul>
+<li><b>segment</b>: integer, segment identifier
+<li><b>next_segment</b>: integer, topological next segment identifier
+<li><b>s_order</b>: integer, segment order
+<li><b>next_order</b>: integer, topological next segment order
+<li><b>direction</b>: double precision, full segment direction (0-360)
+<li><b>azimuth</b>: double precision, full segment azimuth (0-180) 
+<li><b>length</b>: double precision, segment length
+<li><b>straight</b>: double precision, length of straight line between segment
+nodes
+<li><b>sinusoid</b>: double precision, sinusoid (length/straight)
+<li><b>elev_min</b>: double precision, minimum elevation (elevation at segment
+start)
+<li><b>elev_max</b>: double precision, maximum elevation (elevation at segment
+end)
+<li><b>s_drop</b>: double precision, difference between start and end of the
+segment
+<li><b>gradient</b>: double precision, drop/length
+<li><b>out_direction</b>: double precision, direction (0-360) of segment end
+sector
+<li><b>out_azimuth</b>: double precision,  azimuth (0-180) of segment end sector
+<li><b>out_length</b>: double precision, length of segment end sector
+<li><b>out_drop</b>: double precision, drop of segment end sector
+<li><b>out_gradient</b>: double precision, gradient of segment end sector
+<li><b>tangent_dir</b>: double precision, direction of tangent in segment outlet
+to the next stream 
+<li><b>tangent_azimuth</b>: double precision, azimuth of tangent in segment
+outlet to the next stream 
+<li><b>next_direction</b>: double precision, direction of next stream in join
+with current segment 
+<li><b>next_azimuth</b>: double precision, azimuth of next stream in join with
+current segment 
+</ul>
+<img src="dirs.png">
+</dd>
+<dt><b>sectors</b></dt>
+<dd>Vector map where every sector has its own category and following attributes:
+<ul>
+<li><b>sector</b>: integer, sector category
+<li><b>segment</b>: integer, segment category (to establish relationship)
+<li><b>s_order</b>: integer, segment order
+<li><b>direction</b>: double precision, sector direction
+<li><b>azimuth</b>: double precision, sector azimuth
+<li><b>length</b>: double precision, sector length
+<li><b>straight</b>: double precision, length of straight line between sector
+nodes
+<li><b>sinusoid</b>: double precision, sinusoid (length/straight)
+<li><b>elev_min</b>: double precision, minimum elevation (elevation at sector
+start)
+<li><b>elev_max</b>: double precision, minimum elevation (elevation at sector
+end)
+<li><b>s_drop</b>: double precision, difference between start and end of the
+sector
+<li><b>gradient</b>: double precision, drop/length
+</ul>
+<img src="sectors.png">
+Relation between segments and sector may be set up by segment key.
+</dd>
+</dl>
+<p>
+The main idea comes from works of Horton (1932) and Howard (1971, 1990). The
+module is designed to investigate network lineaments and calculate angle
+relations between tributaries and its major streams. The main problem in
+calculating directional parameters is that streams usually are not straight
+lines. Therefore as the first step of the procedure, partitioning of streams
+into near-straight-line segments is required.
+<p>
+The segmentation process uses a method similar to the one used by Van & Ventura
+(1997) to detect corners and partition curves into straight lines and gentle
+arcs. Because it is almost impossible to determine exactly straight sections
+without creating numerous very short segments, the division process requires
+some approximation. The approximation is made based on three parameters: (1) the
+downstream/upstream search length, (2) the short segment skipping threshold, and
+(3) the maximum angle between downstream/upstream segments to be considered as a
+straight line. In order to designate straight sections of the streams, the
+algorithm is searching for those points where curves significantly change their
+direction.
+The definition of stream segments depends on the ordering method selected by the
+user, Strahler's, Horton's or Hack's main stream, or the network may remain
+unordered. All junctions of streams to streams of higher order are always split
+points, but for ordered networks, streams of higher order may be divided into
+sections which ignore junctions with streams of lower order. In unordered
+networks all junctions are always split points.
+In extended mode the module also calculates the direction of a stream to its
+higher order stream. If the higher order stream begins at the junction with the
+current stream (Strahler's ordering only) or if the network is unordered, the
+direction is calculated as the direction of the line between junction point and
+downstream point (Howard 1971) within the user-defined global search distance.
+If a higher order stream continues at the junction, its direction is calculated
+as the direction of the tangent line to the stream of higher order at the
+junction point. To avoid local fluctuation, the tangent line is approximated as
+a secant line joining downstream/upstream points at a distance globally defined
+by the search length parameter (1). Such a definition of the angle between
+streams is not fully compatible with Horton's original criterion.
+
+<h2>NOTES</h2>
+<p>
+The module can work only if direction map, stream_rast map and region have the same
+settings. It is also required that stream_rast map and direction map come from the
+same source. For lots of reason this limitation probably cannot be omitted.  
+This means that if stream_rast map comes from r.stream.extract, also direction map from
+r.stream.extract must be used. If stream_rast network was generated with MFD method,
+also MFD direction map must be used.
+
+<h2>EXAMPLE</h2>
+
+<div class="code"><pre>
+g.region -p -a rast=elevation
+r.watershed elevation=elevation threshold=10000 drainage=direction stream=streams
+r.stream.order stream_vect=streams direction=direction strahler=riverorder_strahler
+r.stream.segment stream_rast=riverorder_strahler direction=direction \
+  elevation=elevation segments=river_segment sectors=river_sector
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.stream.extract.html">r.stream.extract</a>,
+<a href="r.stream.order.html">r.stream.order</a>,
+<a href="r.stream.basins.html">r.stream.basins</a>,
+<a href="r.stream.slope.html">r.stream.slope</a>,
+<a href="r.stream.stats.html">r.stream.stats</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+</em>
+
+<h2>REFERENCES</h2>
+<p>Horton, R. E., (1932). Drainage basin characteristics: Am. Geophys. Union
+Trans., (3), 350-361.
+<p>Howard, A.D. (1971). Optimal angles of stream junction: Geometric, Stability
+to capture and Minimum Power Criteria, Water Resour. Res. 7(4), 863-873.
+<p>Howard, A.D. (1990). Theoretical model of optimal drainage networks Water
+Resour. Res., 26(9),  2107-2117.
+<p>Van, W., Ventura, J.A. (1997). Segmentation of Planar Curves into
+Straight-Line Segments and Elliptical Arcs, Graphical Models and Image
+Processing 59(6), 484-494.
+
+<h2>AUTHOR</h2>
+Jarek Jasiewicz, Adam Mickiewicz University, Geoecology and Geoinformation
+Institute.
+
+<p><i>Last changed: $Date$</i>
+

raster/r.stream.segment/stream_segment.c

@@ -0,0 +1,284 @@
+#include "local_proto.h"
+static int sector_cat = 0;
+
+float calc_dir(int rp, int cp, int rn, int cn)
+{
+    return
+	(cp - cn) == 0 ?
+	(rp - rn) > 0 ? 0 : PI :
+	(cp - cn) < 0 ?
+	PI / 2 + atan((rp - rn) / (float)(cp - cn)) :
+	3 * PI / 2 + atan((rp - rn) / (float)(cp - cn));
+}
+
+float calc_length(double *distance, int start, int stop)
+{
+    float cum_length = 0;
+    int i;
+
+    for (i = start; i < stop; ++i)
+	cum_length += distance[i];
+    return cum_length;
+}
+
+double calc_drop(float *elevation, int start, int stop)
+{
+    float result;
+
+    result = elevation[start] - elevation[stop];
+    return result < 0 ? 0 : result;
+}
+
+double calc_stright(int rp, int cp, int rn, int cn)
+{
+    double northing, easting, next_northing, next_easting;
+
+    northing = window.north - (rp + .5) * window.ns_res;
+    easting = window.west + (cp + .5) * window.ew_res;
+    next_northing = window.north - (rn + .5) * window.ns_res;
+    next_easting = window.west + (cn + .5) * window.ew_res;
+    return G_distance(easting, northing, next_easting, next_northing);
+}
+
+int create_sectors(STREAM *cur_stream, int seg_length, int seg_skip,
+		   double seg_treshold)
+{
+    DIRCELLS *streamline;
+    unsigned long int *P;	/* alias for points */
+
+    int i, prev_i = 0;
+    int number_of_cells;
+    int cell_down, cell_up;
+    int r, c, r_up, c_up, r_down, c_down;
+    int seg_length_short = seg_length / 3;
+    float dir_down, dir_up, dir_diff;
+
+    float local_minimum = PI;
+    int number_of_sectors = 0;
+    int local_minimum_point = 0;
+    int sector_index = 0;
+    int in_loop = 0;
+    int num_of_points = 0, num_of_breakpoints = 0;
+
+    number_of_cells = cur_stream->number_of_cells - 1;
+    P = cur_stream->points;
+
+    streamline =
+	(DIRCELLS *) G_malloc((number_of_cells + 1) * sizeof(DIRCELLS));
+
+    /* init cells */
+    for (i = 0; i < number_of_cells + 1; ++i) {
+	streamline[i].long_dir_diff = 0;
+	streamline[i].short_dir_diff = 0;
+	streamline[i].long_break = 0;
+	streamline[i].decision = 0;
+    }
+
+    /* upstream: to init, downstream: to outlet */
+    for (i = seg_skip; i < number_of_cells - seg_skip; ++i) {
+	cell_up = i < seg_length ? i : seg_length;
+	cell_down = i > number_of_cells - 1 - seg_length ?
+	    number_of_cells - 1 - i : seg_length;
+
+	r = (int)P[i] / ncols;
+	c = (int)P[i] % ncols;
+	r_up = (int)P[i - cell_up] / ncols;
+	c_up = (int)P[i - cell_up] % ncols;
+	r_down = (int)P[i + cell_down] / ncols;
+	c_down = (int)P[i + cell_down] % ncols;
+
+	dir_down = calc_dir(r, c, r_down, c_down);
+	dir_up = calc_dir(r, c, r_up, c_up);
+	dir_diff = fabs(dir_up - dir_down);
+	streamline[i].long_dir_diff =
+	    dir_diff > PI ? PI * 2 - dir_diff : dir_diff;
+	streamline[i].long_break =
+	    (streamline[i].long_dir_diff < seg_treshold) ? 1 : 0;
+
+	cell_up = i < seg_length_short ? i : seg_length_short;
+	cell_down = i > number_of_cells - 1 - seg_length_short ?
+	    number_of_cells - 1 - i : seg_length_short;
+
+	r = (int)P[i] / ncols;
+	c = (int)P[i] % ncols;
+	r_up = (int)P[i - cell_up] / ncols;
+	c_up = (int)P[i - cell_up] % ncols;
+	r_down = (int)P[i + cell_down] / ncols;
+	c_down = (int)P[i + cell_down] % ncols;
+
+	dir_down = calc_dir(r, c, r_down, c_down);
+	dir_up = calc_dir(r, c, r_up, c_up);
+	dir_diff = fabs(dir_up - dir_down);
+	streamline[i].short_dir_diff =
+	    dir_diff > PI ? PI * 2 - dir_diff : dir_diff;
+    }
+
+    /* look for breakpoints */
+    for (i = 0; i < number_of_cells; ++i) {
+
+	if (streamline[i].long_break) {
+	    num_of_breakpoints = 0;
+	    if (local_minimum > streamline[i].short_dir_diff) {
+		local_minimum = streamline[i].short_dir_diff;
+		local_minimum_point = i;
+		in_loop = 1;
+	    }			/* end local minimum */
+
+	}
+	else if (!streamline[i].long_break && in_loop) {
+	    num_of_breakpoints++;
+	    if (num_of_breakpoints == (seg_length / 5)) {
+		streamline[local_minimum_point].decision = 1;
+		local_minimum = PI;
+		in_loop = 0;
+	    }
+	}
+    }
+
+    /* cleaning breakpoints */
+    for (i = 0, num_of_points = 0; i < number_of_cells; ++i, ++num_of_points) {
+
+	if (streamline[i].decision) {
+	    //printf("       BEFORE  %d %d\n",i,num_of_points);
+	    if (i < seg_skip || (i > seg_skip && num_of_points < seg_skip)) {
+		streamline[i].decision = 0;
+		i = local_minimum_point;
+	    }
+	    else {
+		local_minimum_point = i;
+	    }
+	    num_of_points = 0;
+	}
+    }
+
+    /* number of segment in streamline */
+    for (i = 0; i < number_of_cells + 1; ++i)
+	if (streamline[i].decision == 1 || i == (number_of_cells - 1))
+	    number_of_sectors++;
+
+
+    cur_stream->number_of_sectors = number_of_sectors;
+    cur_stream->sector_breakpoints =
+	(int *)G_malloc(number_of_sectors * sizeof(int));
+    cur_stream->sector_cats =
+	(int *)G_malloc(number_of_sectors * sizeof(int));
+    cur_stream->sector_directions =
+	(float *)G_malloc(number_of_sectors * sizeof(float));
+    cur_stream->sector_strights =
+	(float *)G_malloc(number_of_sectors * sizeof(float));
+    cur_stream->sector_lengths =
+	(double *)G_malloc(number_of_sectors * sizeof(double));
+    cur_stream->sector_drops =
+	(float *)G_malloc(number_of_sectors * sizeof(float));
+
+    /* add attributes */
+    for (i = 0, prev_i = 0; i < number_of_cells + 1; ++i) {
+	if (streamline[i].decision == 1 || i == (number_of_cells - 1)) {
+
+	    r = (int)P[i] / ncols;
+	    c = (int)P[i] % ncols;
+	    r_up = (int)P[prev_i] / ncols;
+	    c_up = (int)P[prev_i] % ncols;
+
+	    cur_stream->sector_breakpoints[sector_index] = i;
+
+	    cur_stream->sector_directions[sector_index] =
+		calc_dir(r_up, c_up, r, c);
+
+	    cur_stream->sector_lengths[sector_index] =
+		calc_length(cur_stream->distance, prev_i, i);
+
+	    cur_stream->sector_strights[sector_index] =
+		calc_stright(r_up, c_up, r, c);
+
+	    cur_stream->sector_drops[sector_index] =
+		calc_drop(cur_stream->elevation, prev_i, i);
+
+	    cur_stream->sector_cats[sector_index] = ++sector_cat;
+	    sector_index++;
+	    if (i < (number_of_cells - 1))
+		prev_i = i;
+	}
+    }
+
+
+    /*
+       for (i = 0; i < number_of_cells; ++i)
+       printf("%d | %f  %f |break %d | Dec %d  \n" ,i,
+       streamline[i].long_dir_diff,
+       streamline[i].short_dir_diff,
+       streamline[i].long_break,
+       streamline[i].decision); 
+
+     */
+    G_free(streamline);
+    return 0;
+}
+
+int calc_tangents(STREAM *cur_stream, int seg_length, int seg_skip,
+		  int number_streams)
+{
+
+    int i;
+    int cell_up, cell_down;
+    int r, c, r_up, c_up, r_down, c_down;
+    STREAM *SA = stream_attributes;
+    unsigned long int *P = cur_stream->points;
+    int next_stream = cur_stream->next_stream;
+    unsigned int outlet = cur_stream->outlet;
+    int last_cell = cur_stream->number_of_cells - 1;
+    int reached_end = 1;
+
+    G_debug(3, "calc_tangents(): number_streams=%d", number_streams);
+    /*before calc tangents add rest of streamline attributes */
+    r_up = (int)P[1] / ncols;
+    c_up = (int)P[1] % ncols;
+    r_down = (int)P[last_cell] / ncols;
+    c_down = (int)P[last_cell] % ncols;
+
+    cur_stream->direction = calc_dir(r_up, c_up, r_down, c_down);
+    cur_stream->length = calc_length(cur_stream->distance, 1, last_cell);
+    cur_stream->stright = calc_stright(r_up, c_up, r_down, c_down);
+    cur_stream->drop = calc_drop(cur_stream->elevation, 1, last_cell);
+
+    if (next_stream < 1) {
+	cur_stream->tangent = -1;
+	cur_stream->continuation = -1;
+	return 0;
+    }
+
+    /* find location of outlet in next stream */
+    for (i = 1; i < SA[next_stream].number_of_cells; ++i) {
+	if (SA[next_stream].points[i] == outlet) {
+	    reached_end = 0;
+	    break;
+	}
+    }
+
+    /* outlet not lies on the next stream */
+    if (reached_end) {
+	G_warning(_("Network topology error: cannot identify stream join for stream %d"),
+		  cur_stream->stream);
+	cur_stream->tangent = -1;
+	cur_stream->continuation = -1;
+	return 0;
+    }
+
+    cell_up = i <= seg_length ? i - 1 : seg_length;
+    cell_down = i >= (SA[next_stream].number_of_cells - seg_length) ?
+	SA[next_stream].number_of_cells - seg_length - 1 : seg_length;
+
+    r = (int)SA[next_stream].points[i] / ncols;
+    c = (int)SA[next_stream].points[i] % ncols;
+    r_up = (int)SA[next_stream].points[i - cell_up] / ncols;
+    c_up = (int)SA[next_stream].points[i - cell_up] % ncols;
+    r_down = (int)SA[next_stream].points[i + cell_down] / ncols;
+    c_down = (int)SA[next_stream].points[i + cell_down] % ncols;
+
+    cur_stream->continuation = calc_dir(r, c, r_down, c_down);
+    cur_stream->tangent = i == 1 ? -1 :
+	i < seg_skip ? cur_stream->continuation : calc_dir(r_up, c_up, r_down,
+							   c_down);
+
+    return 0;
+}

raster/r.stream.segment/stream_topology.c

@@ -0,0 +1,570 @@
+#include "local_proto.h"
+double get_distance(int r, int c, int d)
+{
+    double northing, easting, next_northing, next_easting;
+    int next_r, next_c;
+
+    next_r = NR(d);
+    next_c = NC(d);
+    northing = window.north - (r + .5) * window.ns_res;
+    easting = window.west + (c + .5) * window.ew_res;
+    next_northing = window.north - (next_r + .5) * window.ns_res;
+    next_easting = window.west + (next_c + .5) * window.ew_res;
+    return G_distance(easting, northing, next_easting, next_northing);
+}
+
+int ram_trib_nums(int r, int c, CELL ** streams, CELL ** dirs)
+{				/* calculate number of tributaries */
+
+    int trib_num = 0;
+    int i, j;
+    int next_r, next_c;
+
+    for (i = 1; i < 9; ++i) {
+	if (NOT_IN_REGION(i))
+	    continue;
+
+	j = DIAG(i);
+	next_r = NR(i);
+	next_c = NC(i);
+
+	if (streams[next_r][next_c] > 0 && dirs[next_r][next_c] == j)
+	    trib_num++;
+    }
+
+    if (trib_num > 1)
+	for (i = 1; i < 9; ++i) {
+	    if (NOT_IN_REGION(i))
+		continue;
+
+	    j = DIAG(i);
+	    next_r = NR(i);
+	    next_c = NC(i);
+
+	    if (streams[next_r][next_c] == streams[r][c] &&
+		dirs[next_r][next_c] == j)
+		trib_num--;
+	}
+
+    if (trib_num > 5)
+	G_fatal_error(_("Error finding inits. Stream and direction maps probably do not match."));
+    if (trib_num > 3)
+	G_warning(_("Stream network may be too dense"));
+
+    return trib_num;
+}				/* end trib_num */
+
+
+int seg_trib_nums(int r, int c, SEGMENT *streams, SEGMENT *dirs)
+{				/* calculate number of tributaries */
+
+    int trib_num = 0;
+    int i, j;
+    int next_r, next_c;
+    int streams_cell, streams_next_cell, dirs_next_cell;
+
+    segment_get(streams, &streams_cell, r, c);
+    for (i = 1; i < 9; ++i) {
+	if (NOT_IN_REGION(i))
+	    continue;
+
+	j = DIAG(i);
+	next_r = NR(i);
+	next_c = NC(i);
+
+	segment_get(streams, &streams_next_cell, next_r, next_c);
+	segment_get(dirs, &dirs_next_cell, next_r, next_c);
+
+	if (streams_next_cell > 0 && dirs_next_cell == j)
+	    trib_num++;
+    }
+
+    if (trib_num > 1)
+	for (i = 1; i < 9; ++i) {
+	    if (NOT_IN_REGION(i))
+		continue;
+
+	    j = DIAG(i);
+	    next_r = NR(i);
+	    next_c = NC(i);
+
+	    segment_get(streams, &streams_next_cell, next_r, next_c);
+	    segment_get(dirs, &dirs_next_cell, next_r, next_c);
+
+	    if (streams_next_cell == streams_cell && dirs_next_cell == j)
+		trib_num--;
+	}
+
+    if (trib_num > 5)
+	G_fatal_error(_("Error finding inits. Stream and direction maps probably do not match."));
+    if (trib_num > 3)
+	G_warning(_("Stream network may be too dense"));
+
+    return trib_num;
+}				/* end trib_num */
+
+
+int ram_number_of_streams(CELL **streams, CELL **dirs, int *ordered)
+{
+    int r, c;
+    int stream_num = 0;
+    int one = 0, two = 0;
+
+    for (r = 0; r < nrows; ++r)
+	for (c = 0; c < ncols; ++c)
+	    if (streams[r][c] > 0)
+		if (ram_trib_nums(r, c, streams, dirs) != 1) {
+		    stream_num++;
+		    if (streams[r][c] == 1)
+			one++;
+		    if (streams[r][c] == 2)
+			two++;
+		}
+    *ordered = (one > 1 || two > 1) ? 1 : 0;
+    /* if there is more than 1 stream with identifier 1 or 2  network is ordered */
+
+    return stream_num;
+}
+
+int seg_number_of_streams(SEGMENT *streams, SEGMENT *dirs, int *ordered)
+{
+    int r, c;
+    int stream_num = 0;
+    int one = 0, two = 0;
+    int streams_cell;
+
+    for (r = 0; r < nrows; ++r)
+	for (c = 0; c < ncols; ++c) {
+	    segment_get(streams, &streams_cell, r, c);
+	    if (streams_cell > 0)
+		if (seg_trib_nums(r, c, streams, dirs) != 1) {
+		    stream_num++;
+		    if (streams_cell == 1)
+			one++;
+		    if (streams_cell == 2)
+			two++;
+		}
+	}
+    *ordered = (one > 1 || two > 1) ? 1 : 0;
+    /* if there is more than 1 stream with identifier 1 or 2  network is ordered */
+
+    return stream_num;
+}
+
+int ram_build_streamlines(CELL **streams, CELL **dirs, FCELL **elevation,
+			  int number_of_streams)
+{
+    int r, c, i;
+    int d, next_d;
+    int prev_r, prev_c;
+    int stream_num = 1, cell_num = 0;
+    int contrib_cell;
+    STREAM *SA;
+    int border_dir;
+
+    stream_attributes =
+	(STREAM *) G_malloc(number_of_streams * sizeof(STREAM));
+    G_message(_("Finding inits..."));
+    SA = stream_attributes;
+
+    for (r = 0; r < nrows; ++r)
+	for (c = 0; c < ncols; ++c)
+	    if (streams[r][c])
+		if (ram_trib_nums(r, c, streams, dirs) != 1) {	/* adding inits */
+		    if (stream_num > number_of_streams)
+			G_fatal_error(_("Error finding inits. Stream and direction maps probably do not match."));
+
+		    SA[stream_num].stream = stream_num;
+		    SA[stream_num].init = INDEX(r, c);
+		    stream_num++;
+		}
+
+    for (i = 1; i < stream_num; ++i) {
+
+
+	r = (int)SA[i].init / ncols;
+	c = (int)SA[i].init % ncols;
+	SA[i].order = streams[r][c];
+	SA[i].number_of_cells = 0;
+	do {
+
+	    SA[i].number_of_cells++;
+	    d = abs(dirs[r][c]);
+	    if (NOT_IN_REGION(d) || d == 0)
+		break;
+	    r = NR(d);
+	    c = NC(d);
+	} while (streams[r][c] == SA[i].order);
+
+	SA[i].number_of_cells += 2;	/* add two extra points for init+ and outlet+ */
+    }
+
+    for (i = 1; i < number_of_streams; ++i) {
+
+	SA[i].points = (unsigned long int *)
+	    G_malloc((SA[i].number_of_cells) * sizeof(unsigned long int));
+	SA[i].elevation = (float *)
+	    G_malloc((SA[i].number_of_cells) * sizeof(float));
+	SA[i].distance = (double *)
+	    G_malloc((SA[i].number_of_cells) * sizeof(double));
+
+	r = (int)SA[i].init / ncols;
+	c = (int)SA[i].init % ncols;
+	contrib_cell = ram_find_contributing_cell(r, c, dirs, elevation);
+	prev_r = NR(contrib_cell);
+	prev_c = NC(contrib_cell);
+
+	/* add one point contributing to init to calculate parameters */
+	/* what to do if there is no contributing points? */
+	SA[i].points[0] = (contrib_cell == 0) ? -1 : INDEX(prev_r, prev_c);
+	SA[i].elevation[0] = (contrib_cell == 0) ? -99999 :
+	    elevation[prev_r][prev_c];
+	d = (contrib_cell == 0) ? dirs[r][c] : dirs[prev_r][prev_c];
+	SA[i].distance[0] = (contrib_cell == 0) ? get_distance(r, c, d) :
+	    get_distance(prev_r, prev_c, d);
+
+	SA[i].points[1] = INDEX(r, c);
+	SA[i].elevation[1] = elevation[r][c];
+	d = abs(dirs[r][c]);
+	SA[i].distance[1] = get_distance(r, c, d);
+
+	cell_num = 2;
+	do {
+	    d = abs(dirs[r][c]);
+
+	    if (NOT_IN_REGION(d) || d == 0) {
+		SA[i].points[cell_num] = -1;
+		SA[i].distance[cell_num] = SA[i].distance[cell_num - 1];
+		SA[i].elevation[cell_num] =
+		    2 * SA[i].elevation[cell_num - 1] -
+		    SA[i].elevation[cell_num - 2];
+		border_dir = convert_border_dir(r, c, dirs[r][c]);
+		SA[i].last_cell_dir = border_dir;
+		break;
+	    }
+	    r = NR(d);
+	    c = NC(d);
+	    SA[i].last_cell_dir = dirs[r][c];
+	    SA[i].points[cell_num] = INDEX(r, c);
+	    SA[i].elevation[cell_num] = elevation[r][c];
+	    next_d = (abs(dirs[r][c]) == 0) ? d : abs(dirs[r][c]);
+	    SA[i].distance[cell_num] = get_distance(r, c, next_d);
+	    cell_num++;
+	    if (cell_num > SA[i].number_of_cells)
+		G_fatal_error(_("To much points in stream line"));
+	} while (streams[r][c] == SA[i].order);
+
+	if (SA[i].elevation[0] == -99999)
+	    SA[i].elevation[0] = 2 * SA[i].elevation[1] - SA[i].elevation[2];
+    }
+    return 0;
+}
+
+
+int seg_build_streamlines(SEGMENT *streams, SEGMENT *dirs,
+			  SEGMENT *elevation, int number_of_streams)
+{
+    int r, c, i;
+    int d, next_d;
+    int prev_r, prev_c;
+    int stream_num = 1, cell_num = 0;
+    int contrib_cell;
+    STREAM *SA;
+    int border_dir;
+    int streams_cell, dirs_cell;
+    int dirs_prev_cell;
+    float elevation_cell, elevation_prev_cell;
+
+    stream_attributes =
+	(STREAM *) G_malloc(number_of_streams * sizeof(STREAM));
+    G_message(_("Finding inits..."));
+    SA = stream_attributes;
+
+    /* finding inits */
+    for (r = 0; r < nrows; ++r)
+	for (c = 0; c < ncols; ++c) {
+	    segment_get(streams, &streams_cell, r, c);
+
+	    if (streams_cell)
+		if (seg_trib_nums(r, c, streams, dirs) != 1) {	/* adding inits */
+		    if (stream_num > number_of_streams)
+			G_fatal_error(_("Error finding inits. Stream and direction maps probably do not match."));
+
+		    SA[stream_num].stream = stream_num;
+		    SA[stream_num].init = INDEX(r, c);
+		    stream_num++;
+		}
+	}
+
+    /* building streamline */
+    for (i = 1; i < stream_num; ++i) {
+
+	r = (int)SA[i].init / ncols;
+	c = (int)SA[i].init % ncols;
+	segment_get(streams, &streams_cell, r, c);
+	SA[i].order = streams_cell;
+	SA[i].number_of_cells = 0;
+
+	do {
+
+	    SA[i].number_of_cells++;
+	    segment_get(dirs, &dirs_cell, r, c);
+
+	    d = abs(dirs_cell);
+	    if (NOT_IN_REGION(d) || d == 0)
+		break;
+	    r = NR(d);
+	    c = NC(d);
+	    segment_get(streams, &streams_cell, r, c);
+	} while (streams_cell == SA[i].order);
+
+	SA[i].number_of_cells += 2;	/* add two extra points for point before init and after outlet */
+    }
+
+    for (i = 1; i < number_of_streams; ++i) {
+
+	SA[i].points = (unsigned long int *)
+	    G_malloc((SA[i].number_of_cells) * sizeof(unsigned long int));
+	SA[i].elevation = (float *)
+	    G_malloc((SA[i].number_of_cells) * sizeof(float));
+	SA[i].distance = (double *)
+	    G_malloc((SA[i].number_of_cells) * sizeof(double));
+
+	r = (int)SA[i].init / ncols;
+	c = (int)SA[i].init % ncols;
+	contrib_cell = seg_find_contributing_cell(r, c, dirs, elevation);
+	prev_r = NR(contrib_cell);
+	prev_c = NC(contrib_cell);
+
+	/* add one point contributing to init to calculate parameters */
+	/* what to do if there is no contributing points? */
+
+	segment_get(dirs, &dirs_cell, r, c);
+	segment_get(dirs, &dirs_prev_cell, prev_r, prev_c);
+	segment_get(elevation, &elevation_prev_cell, prev_r, prev_c);
+	segment_get(elevation, &elevation_cell, r, c);
+
+	SA[i].points[0] = (contrib_cell == 0) ? -1 : INDEX(prev_r, prev_c);
+	SA[i].elevation[0] = (contrib_cell == 0) ? -99999 :
+	    elevation_prev_cell;
+	d = (contrib_cell == 0) ? dirs_cell : dirs_prev_cell;
+	SA[i].distance[0] = (contrib_cell == 0) ? get_distance(r, c, d) :
+	    get_distance(prev_r, prev_c, d);
+
+	SA[i].points[1] = INDEX(r, c);
+	SA[i].elevation[1] = elevation_cell;
+	d = abs(dirs_cell);
+	SA[i].distance[1] = get_distance(r, c, d);
+
+	cell_num = 2;
+	do {
+	    segment_get(dirs, &dirs_cell, r, c);
+	    d = abs(dirs_cell);
+
+	    if (NOT_IN_REGION(d) || d == 0) {
+		SA[i].points[cell_num] = -1;
+		SA[i].distance[cell_num] = SA[i].distance[cell_num - 1];
+		SA[i].elevation[cell_num] =
+		    2 * SA[i].elevation[cell_num - 1] -
+		    SA[i].elevation[cell_num - 2];
+		border_dir = convert_border_dir(r, c, dirs_cell);
+		SA[i].last_cell_dir = border_dir;
+		break;
+	    }
+	    r = NR(d);
+	    c = NC(d);
+	    segment_get(dirs, &dirs_cell, r, c);
+	    SA[i].last_cell_dir = dirs_cell;
+	    SA[i].points[cell_num] = INDEX(r, c);
+	    segment_get(elevation, &SA[i].elevation[cell_num], r, c);
+	    next_d = (abs(dirs_cell) == 0) ? d : abs(dirs_cell);
+	    SA[i].distance[cell_num] = get_distance(r, c, next_d);
+	    cell_num++;
+	    if (cell_num > SA[i].number_of_cells)
+		G_fatal_error(_("To much points in stream line"));
+	    segment_get(streams, &streams_cell, r, c);
+	} while (streams_cell == SA[i].order);
+
+	if (SA[i].elevation[0] == -99999)
+	    SA[i].elevation[0] = 2 * SA[i].elevation[1] - SA[i].elevation[2];
+    }
+    return 0;
+}
+
+int ram_find_contributing_cell(int r, int c, CELL **dirs, FCELL **elevation)
+{
+    int i, j = 0;
+    int next_r, next_c;
+    float elev_min = 9999;
+
+    for (i = 1; i < 9; ++i) {
+	if (NOT_IN_REGION(i))
+	    continue;
+	next_r = NR(i);
+	next_c = NC(i);
+	if (dirs[next_r][next_c] == DIAG(i) &&
+	    elevation[next_r][next_c] < elev_min) {
+	    elev_min = elevation[next_r][next_c];
+	    j = i;
+	}
+    }
+
+    return j;
+}
+
+int seg_find_contributing_cell(int r, int c, SEGMENT *dirs,
+			       SEGMENT *elevation)
+{
+    int i, j = 0;
+    int next_r, next_c;
+    float elev_min = 9999;
+    int dirs_next_cell;
+    float elevation_next_cell;
+
+    for (i = 1; i < 9; ++i) {
+	if (NOT_IN_REGION(i))
+	    continue;
+	next_r = NR(i);
+	next_c = NC(i);
+	segment_get(elevation, &elevation_next_cell, next_r, next_c);
+	segment_get(dirs, &dirs_next_cell, next_r, next_c);
+
+	if (dirs_next_cell == DIAG(i) && elevation_next_cell < elev_min) {
+	    elev_min = elevation_next_cell;
+	    j = i;
+	}
+    }
+    return j;
+}
+
+int ram_fill_streams(CELL **unique_streams, int number_of_streams)
+{
+    int r, c;
+    int i, j;
+    STREAM *SA;
+
+    SA = stream_attributes;
+
+    for (i = 1; i < number_of_streams; ++i) {
+	for (j = 1; j < SA[i].number_of_cells - 1; ++j) {
+	    r = (int)SA[i].points[j] / ncols;
+	    c = (int)SA[i].points[j] % ncols;
+	    unique_streams[r][c] = SA[i].stream;
+	}
+    }
+    return 0;
+}
+
+int seg_fill_streams(SEGMENT *unique_streams, int number_of_streams)
+{
+    int r, c;
+    int i, j;
+    STREAM *SA;
+
+    SA = stream_attributes;
+
+    for (i = 1; i < number_of_streams; ++i) {
+	for (j = 1; j < SA[i].number_of_cells - 1; ++j) {
+	    r = (int)SA[i].points[j] / ncols;
+	    c = (int)SA[i].points[j] % ncols;
+	    segment_put(unique_streams, &SA[i].stream, r, c);
+	}
+    }
+    return 0;
+}
+
+int ram_identify_next_stream(CELL **streams, int number_of_streams)
+{
+    int r, c;
+    int i;
+    STREAM *SA;
+
+    SA = stream_attributes;
+
+    for (i = 1; i < number_of_streams; ++i) {
+	if (SA[i].points[SA[i].number_of_cells - 1] == -1) {
+	    SA[i].next_stream = -1;
+	    SA[i].outlet = -1;
+	}
+	else {
+	    r = (int)SA[i].points[SA[i].number_of_cells - 1] / ncols;
+	    c = (int)SA[i].points[SA[i].number_of_cells - 1] % ncols;
+	    SA[i].next_stream = streams[r][c];
+	    SA[i].outlet = SA[i].points[SA[i].number_of_cells - 1];
+	}
+    }
+    return 0;
+}
+
+int seg_identify_next_stream(SEGMENT *streams, int number_of_streams)
+{
+    int r, c;
+    int i;
+    STREAM *SA;
+
+    SA = stream_attributes;
+
+    for (i = 1; i < number_of_streams; ++i) {
+	if (SA[i].points[SA[i].number_of_cells - 1] == -1) {
+	    SA[i].next_stream = -1;
+	    SA[i].outlet = -1;
+	}
+	else {
+	    r = (int)SA[i].points[SA[i].number_of_cells - 1] / ncols;
+	    c = (int)SA[i].points[SA[i].number_of_cells - 1] % ncols;
+	    segment_get(streams, &SA[i].next_stream, r, c);
+	    SA[i].outlet = SA[i].points[SA[i].number_of_cells - 1];
+	}
+    }
+    return 0;
+}
+
+
+int free_attributes(int number_of_streams)
+{
+    int i;
+    STREAM *SA;
+
+    SA = stream_attributes;
+
+    for (i = 1; i < number_of_streams; ++i) {
+	G_free(SA[i].points);
+	G_free(SA[i].elevation);
+	G_free(SA[i].distance);
+	G_free(SA[i].sector_breakpoints);
+	G_free(SA[i].sector_cats);
+	G_free(SA[i].sector_directions);
+	G_free(SA[i].sector_lengths);
+	G_free(SA[i].sector_drops);
+    }
+    G_free(stream_attributes);
+    return 0;
+}
+
+
+int convert_border_dir(int r, int c, int dir)
+{
+    /* this function must be added to other modules */
+    /* this is added to fix r.stream.extract issue with broader cell direction */
+    if (dir)
+	return dir;
+
+    if (r == 0 && c == 0)
+	return -3;
+    else if (r == 0 && c == ncols - 1)
+	return -1;
+    else if (r == nrows - 1 && c == ncols - 1)
+	return -7;
+    else if (r == nrows - 1 && c == 0)
+	return -5;
+    else if (r == 0)
+	return -2;
+    else if (r == nrows - 1)
+	return -6;
+    else if (c == 0)
+	return -4;
+    else if (c == ncols - 1)
+	return -8;
+    else
+	return 0;
+}

raster/r.stream.segment/stream_vector.c

@@ -0,0 +1,334 @@
+#include "local_proto.h"
+
+int create_sector_vector(char *out_vector, int number_of_streams, int radians)
+{
+    int i, j, k;
+    int r, c, d;
+    int start, stop;
+    float northing, easting;
+    struct Map_info Out;
+    struct line_pnts *Segments;
+    struct line_cats *Cats;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    struct field_info *Fi;
+    dbString table_name, db_sql, val_string;
+    dbDriver *driver;
+    dbHandle handle;
+    char *cat_col_name = "cat";
+    char buf[1000];
+
+    int sector_category, segment, sector, order;
+    double direction, azimuth, length, stright, sinusoid;
+    double elev_min, elev_max, drop, gradient;
+
+    Segments = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+    Vect_open_new(&Out, out_vector, 0);
+
+    Vect_reset_line(Segments);
+    Vect_reset_cats(Cats);
+
+    for (i = 1; i < number_of_streams; ++i) {
+	stop = 1;
+	for (j = 0; j < SA[i].number_of_sectors; ++j) {
+	    start = stop;
+	    stop = (j == SA[i].number_of_sectors - 1) ?
+		SA[i].sector_breakpoints[j] +
+		1 : SA[i].sector_breakpoints[j] + 1;
+	    Vect_cat_set(Cats, 1, SA[i].sector_cats[j]);
+	    for (k = start; k <= stop; ++k) {
+		if (SA[i].points[k] == -1) {
+		    d = abs(SA[i].last_cell_dir);
+		    r = NR(d);
+		    c = NC(d);
+		}
+		else {
+		    r = (int)SA[i].points[k] / ncols;
+		    c = (int)SA[i].points[k] % ncols;
+		}
+		easting = window.west + (c + .5) * window.ew_res;
+		northing = window.north - (r + .5) * window.ns_res;
+		Vect_append_point(Segments, easting, northing, 0);
+	    }
+	    Vect_write_line(&Out, GV_LINE, Segments, Cats);
+	    Vect_reset_line(Segments);
+	    Vect_reset_cats(Cats);
+	}
+    }
+
+    /* add attributes */
+    db_init_string(&db_sql);
+    db_init_string(&val_string);
+    db_init_string(&table_name);
+    db_init_handle(&handle);
+
+    Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
+    driver = db_start_driver_open_database(Fi->driver, Fi->database);
+    if (driver == NULL) {
+	G_fatal_error(_("Unable to start driver <%s>"), Fi->driver);
+    }
+
+    /* create table */
+    sprintf(buf, "create table %s (%s integer, "
+            "segment integer, "
+            "sector integer, "
+            "s_order integer, "
+            "direction double precision, "
+            "azimuth double precision, "
+            "length double precision, "
+            "stright double precision, "
+            "sinusoid double precision, "
+            "elev_min double precision, "
+            "elev_max double precision, "
+            "s_drop double precision, "
+            "gradient double precision)", Fi->table, cat_col_name);
+    
+    db_set_string(&db_sql, buf);
+
+    if (db_execute_immediate(driver, &db_sql) != DB_OK) {
+	db_close_database(driver);
+	db_shutdown_driver(driver);
+	G_fatal_error(_("Unable to create table: '%s'"), db_get_string(&db_sql));
+    }
+
+    if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
+        G_warning(_("Unable to create create index on table <%s>"), Fi->table);
+
+    if (db_grant_on_table(driver, Fi->table,
+			  DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK)
+	G_fatal_error(_("Unable to grant privileges on table <%s>"), Fi->table);
+
+    db_begin_transaction(driver);
+
+    for (i = 1; i < number_of_streams; ++i) {
+	stop = 1;
+	for (j = 0; j < SA[i].number_of_sectors; ++j) {
+	    start = stop;
+	    stop = SA[i].sector_breakpoints[j];
+
+	    /* calculate and add parameters */
+	    sector_category = SA[i].sector_cats[j];
+	    segment = SA[i].stream;
+	    sector = j + 1;
+	    order = SA[i].order;
+	    direction = SA[i].sector_directions[j];
+	    azimuth = direction <= PI ? direction : direction - PI;
+	    length = SA[i].sector_lengths[j];
+	    stright = SA[i].sector_strights[j];
+	    sinusoid = length / stright;
+	    elev_max = SA[i].elevation[start];
+	    elev_min = SA[i].elevation[stop];
+	    drop = elev_max - elev_min;
+	    gradient = drop / length;
+
+	    if (!radians) {
+		direction = RAD2DEG(direction);
+		azimuth = RAD2DEG(azimuth);
+	    }
+
+	    sprintf(buf, "insert into %s values( %d, %d, %d, %d, "
+                    "%f, %f, %f, %f, %f, "
+                    "%f, %f, %f, %f)", Fi->table, sector_category, segment, sector, order,	/*4 */
+		    direction, azimuth, length, stright, sinusoid,	/*9 */
+		    elev_max, elev_min, drop, gradient);	/*13 */
+            
+	    db_set_string(&db_sql, buf);
+
+	    if (db_execute_immediate(driver, &db_sql) != DB_OK) {
+		db_close_database(driver);
+		db_shutdown_driver(driver);
+		G_fatal_error(_("Unable to inset new row: '%s'"),
+			      db_get_string(&db_sql));
+	    }
+	}
+    }
+    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);
+
+    Vect_hist_command(&Out);
+    Vect_build(&Out);
+    Vect_close(&Out);
+    return 0;
+
+}
+
+int create_segment_vector(char *out_vector, int number_of_streams,
+			  int radians)
+{
+    int i, k;
+    int r, c, d;
+    float northing, easting;
+    struct Map_info Out;
+    struct line_pnts *Segments;
+    struct line_cats *Cats;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    struct field_info *Fi;
+    dbString table_name, db_sql, val_string;
+    dbDriver *driver;
+    dbHandle handle;
+    char *cat_col_name = "cat";
+    char buf[1000];
+
+    /* variables to store table attributes */
+    int last;
+    int segment, next_segment, order, next_order;
+    double direction, azimuth, length, stright, sinusoid;
+    double elev_min, elev_max, drop, gradient;
+    double out_direction, out_azimuth, out_length, out_drop, out_gradient;
+    double tangent_dir, tangent_azimuth, next_direction, next_azimuth;
+
+    Segments = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+    Vect_open_new(&Out, out_vector, 0);
+
+    Vect_reset_line(Segments);
+    Vect_reset_cats(Cats);
+
+    for (i = 1; i < number_of_streams; ++i) {
+	Vect_cat_set(Cats, 1, SA[i].stream);
+	for (k = 1; k < SA[i].number_of_cells; ++k) {
+	    if (SA[i].points[k] == -1) {
+		d = abs(SA[i].last_cell_dir);
+		r = NR(d);
+		c = NC(d);
+	    }
+	    else {
+		r = (int)SA[i].points[k] / ncols;
+		c = (int)SA[i].points[k] % ncols;
+	    }
+	    easting = window.west + (c + .5) * window.ew_res;
+	    northing = window.north - (r + .5) * window.ns_res;
+	    Vect_append_point(Segments, easting, northing, 0);
+	}
+	Vect_write_line(&Out, GV_LINE, Segments, Cats);
+	Vect_reset_line(Segments);
+	Vect_reset_cats(Cats);
+    }
+
+    /* add attributes */
+    db_init_string(&db_sql);
+    db_init_string(&val_string);
+    db_init_string(&table_name);
+    db_init_handle(&handle);
+
+    Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
+    driver = db_start_driver_open_database(Fi->driver, Fi->database);
+    if (driver == NULL) {
+	G_fatal_error(_("Unable to start driver <%s>"), Fi->driver);
+    }
+
+    /* create table */
+    sprintf(buf, "create table %s (%s integer, "
+            "segment integer, "
+            "next_segment integer, "
+            "s_order integer, "
+            "next_order integer, "
+            "direction double precision, "
+            "azimuth double precision, "
+            "length double precision, "
+            "stright double precision, "
+            "sinusoid double precision, "
+            "elev_min double precision, "
+            "elev_max double precision, "
+            "s_drop double precision, "
+            "gradient double precision, "
+            "out_direction double precision, "
+            "out_azimuth double precision, "
+            "out_length double precision, "
+            "out_drop double precision, "
+            "out_gradient double precision, "
+            "tangent_dir double precision, "
+            "tangent_azimuth double precision, "
+            "next_direction double precision, "
+            "next_azimuth double precision)", Fi->table, cat_col_name);
+    
+    db_set_string(&db_sql, buf);
+
+    if (db_execute_immediate(driver, &db_sql) != DB_OK) {
+	db_close_database(driver);
+	db_shutdown_driver(driver);
+	G_fatal_error(_("Unable to create table '%s'"), db_get_string(&db_sql));
+    }
+
+    if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
+	G_warning(_("Unable to create index on table <%s>"), Fi->table);
+
+    if (db_grant_on_table(driver, Fi->table,
+			  DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK)
+	G_fatal_error(_("Unable grant privileges on table <%s>"), Fi->table);
+
+    db_begin_transaction(driver);
+
+    for (i = 1; i < number_of_streams; ++i) {
+	/* calculate and add parameters */
+	segment = SA[i].stream;
+	next_segment = SA[i].next_stream;
+	order = SA[i].order;
+	next_order = next_segment == -1 ? -1 : SA[next_segment].order;
+	direction = SA[i].direction;
+	azimuth = direction <= PI ? direction : direction - PI;
+	length = SA[i].length;
+	stright = SA[i].stright;
+	sinusoid = length / stright;
+	elev_max = SA[i].elevation[1];
+	elev_min = SA[i].elevation[SA[i].number_of_cells - 1];
+	drop = SA[i].drop;
+	gradient = drop / length;
+	last = SA[i].number_of_sectors - 1;
+	out_direction = SA[i].sector_directions[last];
+	out_azimuth =
+	    out_direction <= PI ? out_direction : out_direction - PI;
+	out_length = SA[i].sector_lengths[last];
+	out_drop = SA[i].sector_drops[last];
+	out_gradient = out_drop / out_length;
+	tangent_dir = SA[i].tangent;
+	tangent_azimuth = tangent_dir <= PI ? tangent_dir : tangent_dir - PI;
+	next_direction = SA[i].continuation;
+	next_azimuth =
+	    next_direction <= PI ? next_direction : next_direction - PI;
+
+	if (!radians) {
+	    direction = RAD2DEG(direction);
+	    azimuth = RAD2DEG(azimuth);
+	    out_direction = RAD2DEG(out_direction);
+	    out_azimuth = RAD2DEG(out_azimuth);
+	    tangent_dir = RAD2DEG(tangent_dir);
+	    tangent_azimuth = RAD2DEG(tangent_azimuth);
+	    next_direction = RAD2DEG(next_direction);
+	    next_azimuth = RAD2DEG(next_azimuth);
+	}
+
+	sprintf(buf, "insert into %s values( %d, %d, %d, %d, %d, "
+                "%f, %f, %f, %f, %f, "
+                "%f, %f, %f, %f, "
+                "%f, %f, %f, %f, %f, "
+                "%f, %f, %f, %f)", Fi->table, i, segment, next_segment, order, next_order,	/*5 */
+		direction, azimuth, length, stright, sinusoid,	/*10 */
+		elev_max, elev_min, drop, gradient,	/*14 */
+		out_direction, out_azimuth, out_length, out_drop, out_gradient,	/*19 */
+		tangent_dir, tangent_azimuth, next_direction, next_azimuth);	/*23 */
+
+	db_set_string(&db_sql, buf);
+
+	if (db_execute_immediate(driver, &db_sql) != DB_OK) {
+	    db_close_database(driver);
+	    db_shutdown_driver(driver);
+	    G_fatal_error(_("Unable to insert new row: '%s'"),
+			  db_get_string(&db_sql));
+	}
+
+    }
+    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);
+
+    Vect_hist_command(&Out);
+    Vect_build(&Out);
+    Vect_close(&Out);
+    return 0;
+}