[GRASS-SVN] r43756 - grass-addons/grass7/raster/r.stream/r.stream.segment

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Oct 1 07:21:56 EDT 2010


Author: jarekj71
Date: 2010-10-01 11:21:56 +0000 (Fri, 01 Oct 2010)
New Revision: 43756

Added:
   grass-addons/grass7/raster/r.stream/r.stream.segment/Makefile
   grass-addons/grass7/raster/r.stream/r.stream.segment/dirs.png
   grass-addons/grass7/raster/r.stream/r.stream.segment/io.c
   grass-addons/grass7/raster/r.stream/r.stream.segment/io.h
   grass-addons/grass7/raster/r.stream/r.stream.segment/local_proto.h
   grass-addons/grass7/raster/r.stream/r.stream.segment/local_vars.h
   grass-addons/grass7/raster/r.stream/r.stream.segment/main.c
   grass-addons/grass7/raster/r.stream/r.stream.segment/r.stream.segment.html
   grass-addons/grass7/raster/r.stream/r.stream.segment/sectors.png
   grass-addons/grass7/raster/r.stream/r.stream.segment/stream_segment.c
   grass-addons/grass7/raster/r.stream/r.stream.segment/stream_topology.c
   grass-addons/grass7/raster/r.stream/r.stream.segment/stream_vector.c
Log:
files


Added: grass-addons/grass7/raster/r.stream/r.stream.segment/Makefile
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/Makefile	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/Makefile	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../../..
+
+PGM = r.stream.segment
+
+LIBES = $(GISLIB) $(RASTERLIB) $(SEGMENTLIB) $(VECTLIB) $(DBMILIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP) $(SEGMENTDEP) $(VECTDEP) $(DBMIDEP)
+
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/dirs.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.stream/r.stream.segment/dirs.png
___________________________________________________________________
Added: svn:mime-type
   + application/octet-stream

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/io.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/io.c	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/io.c	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,516 @@
+#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, c;
+		
+		if(data_type < 0 || data_type > 2)
+	G_fatal_error(_("ram_creat: Cannot create 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) {
+/*
+ * Funciton 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 neccesary.
+ * input_map_name: name of the map to be read;
+ * map pointer to map stucture (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 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(_("<%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 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(_("ram_open: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) {
+	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 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(_("seg_create: unrecognisabe 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(_("seg_create: cannot format segment"));
+	}
+	
+	close(fd);
+	if(0 > (fd = open(filename,2))) {
+		unlink(filename);
+		G_fatal_error(_("seg_create: cannot re-open file"));
+	}
+	
+	if(0>(fd = segment_init(&(seg->seg),fd,number_of_segs))) {
+		unlink(filename);
+		G_fatal_error(_("seg_create: cannot 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) {
+	
+/*
+ * Funciton 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 neccesary.
+ * input_map_name: name of the map to be read;
+ * seg: pointer to map stucture (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(_("seg_read: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 anyinteger 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 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(_("<%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 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);
+ }
+
+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 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(_("seg_write: Cannot segment read row %d for 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(_("ram_null:Cannot 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;
+}

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/io.h
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/io.h	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/io.h	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,56 @@
+#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/segment.h>
+#include <grass/vector.h>
+#include <grass/dbmi.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(x,y) (x)*ncols+(y)
+#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 requre 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_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*);
+

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/local_proto.h	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/local_proto.h	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,3 @@
+#include "io.h"
+#include "local_vars.h"
+

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/local_vars.h
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/local_vars.h	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/local_vars.h	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,75 @@
+#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;
+
+

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/main.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/main.c	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/main.c	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,294 @@
+
+/****************************************************************************
+ *
+ * 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 stright line portons and 
+ * 							 and segment orientation and angles between streams and its
+ *               tributuaries. For stream direction it use algorithim to divide
+ *               particular streams of the same order  into near-stright line
+ *               portions.
+ * 				
+ *							
+ *
+ * COPYRIGHT:		 (C) 2002,2010 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,j;
+	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 =
+	_("Divide network into near strigh-line segments and calculate its order");
+  G_add_keyword("Stream divide");
+  G_add_keyword("Stream direction");
+  G_add_keyword("Stream gradient");
+
+  in_stm_opt = G_define_standard_option(G_OPT_R_INPUT);	
+  in_stm_opt->key = "streams";
+  in_stm_opt->description = "Name of streams mask input map";
+
+  in_dir_opt = G_define_standard_option(G_OPT_R_INPUT);	
+  in_dir_opt->key = "dirs";
+  in_dir_opt->description = "Name of flow direction input map";
+
+  in_elev_opt = G_define_standard_option(G_OPT_R_INPUT);	
+  in_elev_opt->key = "elevation";
+  in_elev_opt->description = "Name of elevation map";
+    
+	out_segment_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+	out_segment_opt->key="segments";
+	out_segment_opt->description =_("OUTPUT vector file to write segment atributes");
+	
+	out_sector_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+	out_sector_opt->key="sectors";
+	out_sector_opt->description =_("OUTPUT vector file to write segment atributes");
+	
+  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_length->guisection=_("Optional");	
+		
+	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_skip->guisection=_("Optional");	
+		
+	opt_threshold = G_define_option();
+  opt_threshold->key = "treshold";
+  opt_threshold->label = _("Max angle (degrees) beetwen stream segments");
+  opt_threshold->description = _("Must be > 0");
+	opt_threshold->answer = "160";
+	opt_threshold->type = TYPE_DOUBLE;
+	opt_threshold->guisection=_("Optional");	
+	
+  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=_("Optional");	
+    
+  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)");  
+
+      
+    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("Treshold 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_output, 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_output, 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 < 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);
+}
+
+
+
+
+
+/*	
+		int i;
+		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);
+}

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/r.stream.segment.html
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/r.stream.segment.html	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/r.stream.segment.html	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,109 @@
+<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</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 slighty 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>dirs</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 resoultion must be the same. 
+Also <em>stream</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 streams and dirs.</DD>
+<DT><b>length</b></DT>
+<DD>Integer values indicating the search length (in cells) to determine stright 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 neigbour. 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 requre longer length.</DD>
+
+<DT><b>treshold</b></DT>
+<DD>real value indicates the internal angle between upstream and downsteam direction to treat actual cell as lying on the stright line. Grater value (up to 180 degrees) produces more segments. Lesser values produced less segments. Values below 90 in most cases will not produce any addational segments to these resulting from ordering
+</DL>
+<DL>
+<h2>OUTPUTS</h2>
+<P>The module produces two vector maps: one representing original segments (where segment is a streamline where its order remains unchanged) and second divided into near stright line sectors resulting form segmentation proccess. Most of 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>stright</B>: double precision, length of stright line between segment nodes
+<li><B>sinusoid</B>: double precision, sinusoid (length/stright)
+<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, deifference 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 estabilsh 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>stright</B>: double precision, length of stright line between sector nodes
+<li><B>sinusoid</B>: double precision, sinusoid (length/stright)
+<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, deifference 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>
+<h2>DESCRIPTION</h2>
+<P>
+The main idea comes from works of Horton (1932) and Howard (1971, 1990). The module is designed to inverstigate 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>
+Module can work only if direction map, stream map and region map has same settings. It is also required that stream map and direction map come from the same source. For lots of reason this limitation probably cannot be omitted.   this means if stream map comes from r.stream.extract also direction map from r.stream.extract must be used. If stream network was generated with MFD method also MFD direction map must be used.
+
+
+<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.

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/sectors.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.stream/r.stream.segment/sectors.png
___________________________________________________________________
Added: svn:mime-type
   + application/octet-stream

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/stream_segment.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/stream_segment.c	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/stream_segment.c	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,267 @@
+#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, dir_diff_short;
+		
+	float local_minimum=PI;
+	int number_of_sectors=0, cells_in_segment=1;
+	int category_num=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 attributies */
+	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;
+	int outlet=cur_stream->outlet;
+	int last_cell=cur_stream->number_of_cells-1;
+	int reached_end=1;
+	
+		/*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)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;
+	
+	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;
+}
+

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/stream_topology.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/stream_topology.c	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/stream_topology.c	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,572 @@
+#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 tributuaries */
+
+  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 tributuaries */
+
+  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 min=100000000, max=-99999999;
+	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 min=100000000, max=-99999999;
+	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 next_r, next_c;
+  int prev_r, prev_c;
+  int stream_num=1, cell_num=0;
+  int contrib_cell;
+  STREAM* SA;
+	int tmp=0;
+  stream_num = 1;
+	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 next_r, next_c;
+  int prev_r, prev_c;
+  int stream_num=1, cell_num=0;
+  int contrib_cell;
+  STREAM* SA;
+	int tmp=0;
+  stream_num = 1;
+	int border_dir;
+	int streams_cell, dirs_cell;
+	int streams_next_cell, dirs_next_cell;
+	int streams_prev_cell, dirs_prev_cell;
+	float elevation_cell, elevation_next_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++;
+	}
+}
+	
+/* buildning 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,j;
+	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,j;
+	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 broder 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;
+}		

Added: grass-addons/grass7/raster/r.stream/r.stream.segment/stream_vector.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.segment/stream_vector.c	                        (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.segment/stream_vector.c	2010-10-01 11:21:56 UTC (rev 43756)
@@ -0,0 +1,367 @@
+#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=SA[i].sector_breakpoints[j];
+		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);
+
+	//printf("%s \n",buf);
+
+	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("Cannot create table %s", db_get_string(&db_sql));
+    }
+
+    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);
+
+	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*/
+	
+	//printf("%s  \n",buf); }}
+
+
+	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(_("Cannot 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,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];
+  
+  /* 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);
+
+	//printf("%s \n",buf);
+
+	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("Cannot create table %s", db_get_string(&db_sql));
+    }
+
+    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);
+  
+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(direction);
+	out_azimuth=RAD2DEG(azimuth);		
+	next_direction=RAD2DEG(direction);
+	next_azimuth=RAD2DEG(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*/
+	
+	//printf("%s  \n",buf);
+
+
+	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(_("Cannot 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;
+}



More information about the grass-commit mailing list