[GRASS-SVN] r41270 - grass-addons/raster/r.fuzzy
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Mar 3 01:46:46 EST 2010
Author: jarekj71
Date: 2010-03-03 01:46:46 -0500 (Wed, 03 Mar 2010)
New Revision: 41270
Added:
grass-addons/raster/r.fuzzy/Makefile
grass-addons/raster/r.fuzzy/description.html
grass-addons/raster/r.fuzzy/global.h
grass-addons/raster/r.fuzzy/init.c
grass-addons/raster/r.fuzzy/io.c
grass-addons/raster/r.fuzzy/main.c
grass-addons/raster/r.fuzzy/order.c
grass-addons/raster/r.fuzzy/seg_line.c
grass-addons/raster/r.fuzzy/tangent.c
Log:
files for r.fuzzy
Added: grass-addons/raster/r.fuzzy/Makefile
===================================================================
--- grass-addons/raster/r.fuzzy/Makefile (rev 0)
+++ grass-addons/raster/r.fuzzy/Makefile 2010-03-03 06:46:46 UTC (rev 41270)
@@ -0,0 +1,14 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = r.stream.angle
+
+LIBES = $(VECTLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/raster/r.fuzzy/Makefile
___________________________________________________________________
Added: svn:executable
+
Added: grass-addons/raster/r.fuzzy/description.html
===================================================================
--- grass-addons/raster/r.fuzzy/description.html (rev 0)
+++ grass-addons/raster/r.fuzzy/description.html 2010-03-03 06:46:46 UTC (rev 41270)
@@ -0,0 +1,92 @@
+<h2>OPTIONS</h2>
+<DL>
+<DT><b>-r</b></DT>
+<DD>Outputs in radians (Default is degrees). Input is always in degrees</DD>
+<p>
+<DT><b>-e</b></DT>
+<DD>Extended mode. For now supports only direction of streams of higher order and elevation difference between starting and finishig point of the segment.
+</DD>
+<p>
+</DD>
+<p>
+
+<DT><b>stream</b></DT>
+<DD>Stream network: name of input stream map on which segmentation will be performed produced by r.watershed or r.stream.extract. 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>
+<p>
+<DT><b>dir</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>
+<P>
+<DT><b>elev</b></DT>
+<DD>Elevation map. Must be the same as used to calculate streams and dirs. Map can by of any type, but for calculation speed FCELL is the best choice.</DD>
+
+<P>
+<DT><b>order</b></DT>
+<DD>Ordering method. If different than none ordering is choosen, the segmentation proccess ignores junctions and create segments along streams belonging to the same order.</DD>
+
+<P>
+<DT><b>length</b></DT>
+<DD>Integer values indicating the search length (in cells) to determine stright line. The longest length parameter the more tolerant the module treats local stream undulation and inequalities. Default value of 15 is suitable for 30 meters DEMS. More details DEMS may requre longer length.</DD>
+
+<DT><b>skip</b></DT>
+<DD>Integer values indicating the length (in cells) to skip local short segment and join them with the longer neigbor. The shortest length parameter the more short segments will be produced by the module due to undulation and inequalities. Default value of 15 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 that resulting from ordering</DD>
+
+
+<h2>OUTPUTS</h2>
+<P>The module produces one vector map divided into segments resulting form segmentation proccess. Every segmement has ascribed following parameters:</p>
+<ul>
+<li>c_stream: current stream identifier;
+<li>c_order: current stream order according ordering system, 0 if no ordering have been chooded;
+<li>direction: flow downstream direction (degrees or radians) from 0 to 360 (north).
+<li>length: current segment length
+Extended mode (will be improved in the future)
+<li>ncells: number of cells
+<li>n_stream: identifier of the next (topological) stream (not segment)
+<li>n_tangent: direction of tangent line in the point where stream (not segment) reaches its higer order strem. If ordering is NONE it is direction of the upstream part of next toplogical stream in the network.
+</ul>
+</DL>
+
+<h2>DESCRIPTION</h2>
+<P>
+The main idea comes from works of Horton (1932) and Howard (1971, 1990). The module is designed to 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.
+<P>
+Module is still in experimental phase. It means that something may change in the future in the module core.
+
+
+<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.stats.html">r.stream.stats</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="r.reclass.html">r.reclass</a>,
+<a href="r.patch.html">r.patch</a>
+</em>
+
+<h2>AUTHOR</h2>
+Jarek Jasiewicz
+<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.
+
+<HR>
+<P><a href="index.html">Main index</a> - <a href="raster.html">raster index</a> - <a href="full_index.html">Full index</a></P>
+<P>© 2003-2009 <a href="http://grass.osgeo.org">GRASS Development Team</a></p>
+</body>
+</html>
Property changes on: grass-addons/raster/r.fuzzy/description.html
___________________________________________________________________
Added: svn:executable
+
Added: grass-addons/raster/r.fuzzy/global.h
===================================================================
--- grass-addons/raster/r.fuzzy/global.h (rev 0)
+++ grass-addons/raster/r.fuzzy/global.h 2010-03-03 06:46:46 UTC (rev 41270)
@@ -0,0 +1,172 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/dbmi.h>
+#include <grass/Vect.h>
+
+ /* define */
+#define NODES struct node
+#define STREAMS struct stream
+#define HACK 4
+#define HORTON 2
+#define STRAHLER 1
+#define NONE 0
+
+#ifndef PI
+ #define PI (4*atan(1))
+#endif
+
+#define deg2rad(d) ((d)*PI/180)
+#define rad2deg(r) ((r)*180/PI)
+
+
+
+#ifdef MAIN
+# define GLOBAL
+#else
+# define GLOBAL extern
+#endif
+
+/* dirs
+ *** next ***
+
+ 3 | 2 | 1
+ --- --- ---
+ 4 | 0 | 8
+ --- --- ---
+ 5 | 6 | 7
+
+ *** prev ***
+
+ 7 | 6 | 5
+ --- --- ---
+ 8 | 0 | 4
+ --- --- ---
+ 1 | 2 | 3
+ */
+
+#define STREAM struct _stream //zostaje
+STREAM {
+ int stream; /* index */
+ int strahler, shreeve, horton, hack; /* orders */
+ int init; /* stream segment identifier (if ordering is in use) */
+ int next_stream, out_stream;
+ int trib_num;
+ int cell_num;
+ int trib[5];
+ double accum;
+ int init_r, init_c;
+ int out_r, out_c;
+ float tangent_dir; /* tangent direction at init */
+ float stream_dir; /* final stream direction */
+ };
+
+
+#define DIRCELLS struct _dircells
+DIRCELLS {
+ int r;
+ int c;
+ float dir_diff;
+ float small_dir_diff;
+ int candidate;
+ int tunning;
+ int decision;
+ int category;
+};
+
+#define SEGMENTS struct _segments //modyfikacja
+SEGMENTS {
+ int out_stream; /* for streams == stream */
+ int init_stream; /* for streams == stream */
+ int cell_num;
+ int seg_num;
+ float dir_angle; /* final direction */
+ float dir_init, dir_middle, dir_last, dir_full;
+ float *angles;
+ float *lengths;
+ float *drops;
+ int *cellnums;
+ int *cats;
+};
+
+ /* functions.c */
+
+
+/* io.c */
+int open_raster(char *mapname);
+int create_base_maps(void);
+int stream_number(void);
+int close_vector(void);
+
+/* init.c */
+int trib_nums(int r, int c);
+int init_streams(void);
+int find_nodes(void);
+int do_cum_length (void);
+int cur_orders(int cur_stream, int ordering);
+float calc_dir(int rp,int cp,int rn,int cn);
+float calc_drop(int rp,int cp,int rn,int cn);
+float calc_length(int rp,int cp,int rn,int cn);
+
+/* order.c */
+int strahler(void);
+int horton(void);
+int shreeve(void);
+int hack(void);
+
+/* tangent.c */
+int calc_tangent(int ordering);
+int add_missing_dirs(int ordering);
+
+/* seg_line */
+int do_segments(SEGMENTS *segments, int ordering);
+int do_decision(SEGMENTS *segments);
+
+
+
+ /* variables */
+
+
+GLOBAL struct Cell_head window;
+GLOBAL char *in_dirs, *in_streams, *in_elev; /* input dirrection and accumulation raster names */
+
+char *out_vector;
+char *out_table_rels;
+
+struct Map_info Out;
+
+GLOBAL int seg_length;
+GLOBAL int seg_outlet;
+GLOBAL float seg_treshold;
+GLOBAL int radians, extended;
+
+GLOBAL CELL **dirs, **streams; /* matrix with input data */
+GLOBAL FCELL **elevation;
+
+GLOBAL int stack_max;
+GLOBAL int nrows, ncols;
+GLOBAL int ordering;
+
+GLOBAL STREAM *s_streams; /* stream structure all parameters we have here */
+GLOBAL SEGMENTS *seg_common;
+
+/*
+GLOBAL SEGMENTS *seg_strahler;
+GLOBAL SEGMENTS *seg_horton;
+GLOBAL SEGMENTS *seg_hack;
+*/
+
+int *springs, *outlets;
+int springs_num, outlets_num;
+int stream_num;
+
+struct line_pnts *Segments; //
+struct line_cats *Cats; //
+
+
+
+
+
+
Property changes on: grass-addons/raster/r.fuzzy/global.h
___________________________________________________________________
Added: svn:executable
+
Added: grass-addons/raster/r.fuzzy/init.c
===================================================================
--- grass-addons/raster/r.fuzzy/init.c (rev 0)
+++ grass-addons/raster/r.fuzzy/init.c 2010-03-03 06:46:46 UTC (rev 41270)
@@ -0,0 +1,339 @@
+#include <grass/glocale.h>
+#include "global.h"
+
+int trib_nums(int r, int c)
+{ /* calcualte number of tributuaries */
+
+
+ int trib = 0;
+ int i, j;
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ for (i = 1; i < 9; ++i) {
+ if (r + nextr[i] < 0 || r + nextr[i] > (nrows - 1) || c + nextc[i] < 0
+ || c + nextc[i] > (ncols - 1))
+ continue;
+ j = (i + 4) > 8 ? i - 4 : i + 4;
+ if (streams[r + nextr[i]][c + nextc[i]] > 0 &&
+ dirs[r + nextr[i]][c + nextc[i]] == j)
+ trib++;
+ }
+ if (trib > 5)
+ G_fatal_error("Error finding nodes. Stream and direction maps probably do not match...");
+ if (trib > 3)
+ G_warning("Stream network may be too dense...");
+
+ return trib;
+} /* end trib_num */
+
+int init_streams(void)
+{
+ int i;
+ s_streams = (STREAM *) G_malloc((stream_num + 1) * sizeof(STREAM));
+ seg_common=(SEGMENTS *)G_malloc((stream_num + 1) *sizeof(SEGMENTS));
+
+ /*
+ seg_strahler=(SEGMENTS *)G_malloc((stream_num + 1) *sizeof(SEGMENTS));
+ seg_horton=(SEGMENTS *)G_malloc((stream_num + 1) *sizeof(SEGMENTS));
+ seg_hack=(SEGMENTS *)G_malloc((stream_num + 1) *sizeof(SEGMENTS));
+ */
+
+ for (i = 0; i <= stream_num; ++i) {
+ s_streams[i].stream = -1;
+ s_streams[i].next_stream = -1;
+ s_streams[i].out_stream = -1;
+ s_streams[i].strahler = -1;
+ s_streams[i].shreeve = -1;
+ s_streams[i].hack = -1;
+ s_streams[i].horton = -1;
+ s_streams[i].init = -1;
+ s_streams[i].trib_num = 0;
+ s_streams[i].cell_num = 0;
+ s_streams[i].accum = 0.;
+ s_streams[i].trib[0] = 0;
+ s_streams[i].trib[1] = 0;
+ s_streams[i].trib[2] = 0;
+ s_streams[i].trib[3] = 0;
+ s_streams[i].trib[4] = 0;
+ s_streams[i].init_r = -1;
+ s_streams[i].init_c = -1;
+ s_streams[i].out_r = -1;
+ s_streams[i].out_c = -1;
+ s_streams[i].stream_dir=-999;
+ s_streams[i].tangent_dir=-999;
+
+ seg_common[i].init_stream=-1;
+ seg_common[i].out_stream=-1;
+ seg_common[i].cell_num=-1;
+ seg_common[i].seg_num=-1;
+ seg_common[i].dir_angle=-1;
+ seg_common[i].dir_init=-1;
+ seg_common[i].dir_middle=-1;
+ seg_common[i].dir_last=-1;
+ seg_common[i].dir_full=-1;
+ seg_common[i].angles=NULL;
+ seg_common[i].lengths=NULL;
+ seg_common[i].drops=NULL;
+ seg_common[i].cellnums=NULL;
+ seg_common[i].cats=NULL;
+
+ }
+ return 0;
+}
+
+int find_nodes(void)
+{
+
+ int d, i, j; /* d: direction, i: iteration */
+ int r, c;
+ int trib_num, trib = 0;
+ int next_stream = -1, cur_stream;
+
+ springs_num = 0, outlets_num = 0;
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ G_message(_("Finding nodes..."));
+
+ outlets = (int *)G_malloc((stream_num) * sizeof(int));
+ springs = (int *)G_malloc((stream_num) * sizeof(int));
+
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+ if (streams[r][c] > 0) {
+ trib_num = trib_nums(r, c);
+ trib = 0;
+ d = abs(dirs[r][c]); /* abs */
+ if (r + nextr[d] < 0 || r + nextr[d] > (nrows - 1) ||
+ c + nextc[d] < 0 || c + nextc[d] > (ncols - 1)) {
+ next_stream = -1;
+ }
+ else {
+ next_stream = streams[r + nextr[d]][c + nextc[d]];
+ if (next_stream < 1)
+ next_stream = -1;
+ }
+ if (dirs[r][c]==0)
+ next_stream=-1;
+ cur_stream = streams[r][c];
+
+ if (cur_stream != next_stream) { /* building hierarchy */
+
+ if (outlets_num > (stream_num - 1))
+ G_fatal_error("Error finding nodes. Stream and direction maps probably do not match...");
+
+ s_streams[cur_stream].stream = cur_stream;
+ s_streams[cur_stream].next_stream = next_stream;
+ if (next_stream < 0)
+ outlets[outlets_num++] = cur_stream; /* collecting outlets */
+ }
+
+ if (trib_num == 0) { /* adding springs */
+
+ if (springs_num > (stream_num - 1))
+ G_fatal_error("Error finding nodes. Stream and direction maps probably do not match...");
+
+ s_streams[cur_stream].trib_num = 0;
+ springs[springs_num++] = cur_stream; /* collecting springs */
+ s_streams[cur_stream].init_r = r;
+ s_streams[cur_stream].init_c = c;
+ }
+
+ if (trib_num > 1) { /* adding tributuaries */
+ s_streams[cur_stream].trib_num = trib_num;
+ s_streams[cur_stream].init_r = r;
+ s_streams[cur_stream].init_c = c;
+
+ for (i = 1; i < 9; ++i) {
+ if (trib > 4)
+ G_fatal_error("Error finding nodes. Stream and direction maps probably do not match...");
+
+ if (r + nextr[i] < 0 || r + nextr[i] > (nrows - 1) ||
+ c + nextc[i] < 0 || c + nextc[i] > (ncols - 1))
+ continue;
+
+ j = (i + 4) > 8 ? i - 4 : i + 4;
+ if (streams[r + nextr[i]][c + nextc[i]] > 0 &&
+ dirs[r + nextr[i]][c + nextc[i]] == j) {
+ s_streams[cur_stream].trib[trib++] =
+ streams[r + nextr[i]][c + nextc[i]];
+ }
+ } /* end for i... */
+
+ } /* if trib_num>1 */
+ } /* end if streams */
+ } /* c */
+ } /* r */
+ return 0;
+}
+
+int do_cum_length (void) {
+
+ 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 i, s, d; /* s - streams index; d - direction */
+ int done = 1;
+ int r, c;
+ int next_r, next_c;
+ int cur_stream;
+ float cur_northing, cur_easting;
+ float next_northing, next_easting;
+ double cur_length=0.;
+ double cur_accum=0.;
+ G_message("Finding longests streams...");
+
+ for (s=0; s<springs_num; ++s) { /* main loop on springs */
+ r = s_streams[springs[s]].init_r;
+ c = s_streams[springs[s]].init_c;
+ cur_stream=s_streams[springs[s]].stream;
+ cur_length=0;
+ done=1;
+
+ while(done) {
+
+ cur_northing = window.north - (r + .5) * window.ns_res;
+ cur_easting = window.west + (c + .5) * window.ew_res;
+ d=dirs[r][c];
+ next_r=r+nextr[d];
+ next_c=c+nextc[d];
+
+ if (d<1|| /* end of route: sink */
+ r + nextr[d] < 0 || r + nextr[d] > (nrows - 1) || /*border*/
+ c + nextc[d] < 0 || c + nextc[d] > (ncols - 1) ||
+ streams[next_r][next_c]<1) { /* mask */
+
+ cur_length=(window.ns_res+window.ew_res)/2;
+ s_streams[cur_stream].accum += cur_length;
+ cur_length=G_distance(next_easting, next_northing, cur_easting, cur_northing);
+ s_streams[cur_stream].accum += cur_length;
+ s_streams[cur_stream].cell_num++;
+ s_streams[cur_stream].out_r=r;
+ s_streams[cur_stream].out_c=c;
+ break;
+ }
+
+
+ next_northing = window.north - (next_r + .5) * window.ns_res;
+ next_easting = window.west + (next_c + .5) * window.ew_res;
+ cur_length = G_distance(next_easting, next_northing, cur_easting, cur_northing);
+ s_streams[cur_stream].accum += cur_length;
+ s_streams[cur_stream].cell_num++;
+
+ if (streams[next_r][next_c] != cur_stream) {
+ s_streams[cur_stream].out_r=r;
+ s_streams[cur_stream].out_c=c;
+ cur_stream=streams[next_r][next_c];
+ cur_accum=0;
+
+ for (i = 0; i < s_streams[cur_stream].trib_num; ++i) {
+ if (s_streams[s_streams[cur_stream].trib[i]].accum == 0) {
+ done=0;
+ cur_accum=0;
+ break; /* do not pass accum*/
+ }
+ if(s_streams[s_streams[cur_stream].trib[i]].accum>cur_accum)
+ cur_accum=s_streams[s_streams[cur_stream].trib[i]].accum;
+ } /* end for i */
+ s_streams[cur_stream].accum=cur_accum;
+ }/* end if */
+
+ r=next_r;
+ c=next_c;
+
+ } /* end while */
+ } /* end for s*/
+return 0;
+}
+
+int cur_orders(int cur_stream, int ordering) {
+ switch (ordering) {
+ case NONE:
+ return s_streams[cur_stream].stream;
+ break;
+ case HACK:
+ return s_streams[cur_stream].hack;
+ break;
+ case STRAHLER:
+ return s_streams[cur_stream].strahler;
+ break;
+ case HORTON:
+ return s_streams[cur_stream].horton;
+ break;
+ }
+}
+
+/*
+ * rp: previous row (most upstream row)
+ * cp: previous col (most upstream col)
+ * nr: next row (most downstream row)
+ * nc: next col (most downstream col)
+ */
+
+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_drop(int rp,int cp,int rn,int cn) {
+ if ((elevation[rp][cp]-elevation[rn][cn])<0.)
+ return 0.;
+ else
+ return elevation[rp][cp]-elevation[rn][cn];
+}
+
+float calc_length(int rp,int cp,int rn,int cn) {
+
+ 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 i, d; /* s - streams index; d - direction */
+ int done = 1;
+ int r, c;
+ int next_r, next_c;
+ float cur_northing, cur_easting;
+ float next_northing, next_easting;
+ float length;
+
+ r=rp;
+ c=cp;
+ length=0.;
+ done=1;
+
+ while(done) {
+
+ cur_northing = window.north - (r + .5) * window.ns_res;
+ cur_easting = window.west + (c + .5) * window.ew_res;
+ d=dirs[r][c];
+ next_r=r+nextr[d];
+ next_c=c+nextc[d];
+
+ if (d<1|| /* end of route: sink */
+ r + nextr[d] < 0 || r + nextr[d] > (nrows - 1) || /*border*/
+ c + nextc[d] < 0 || c + nextc[d] > (ncols - 1))
+ {
+ length += (window.ns_res+window.ew_res)/2;
+ return length;
+ }
+
+
+ if (r==rn && c==cn) {
+ length += (dirs[r][c]%2) ? 1.41 * (window.ns_res+window.ew_res)/2 : (window.ns_res+window.ew_res)/2;
+ return length;
+ }
+ next_northing = window.north - (next_r + .5) * window.ns_res;
+ next_easting = window.west + (next_c + .5) * window.ew_res;
+ length += G_distance(next_easting, next_northing, cur_easting, cur_northing);
+
+ r=next_r;
+ c=next_c;
+
+
+ } /* end while */
+
+}
Property changes on: grass-addons/raster/r.fuzzy/init.c
___________________________________________________________________
Added: svn:executable
+
Added: grass-addons/raster/r.fuzzy/io.c
===================================================================
--- grass-addons/raster/r.fuzzy/io.c (rev 0)
+++ grass-addons/raster/r.fuzzy/io.c 2010-03-03 06:46:46 UTC (rev 41270)
@@ -0,0 +1,314 @@
+#include <grass/glocale.h>
+#include "global.h"
+
+
+int open_raster(char *mapname)
+{
+ int fd = 0;
+ char *mapset;
+ struct Cell_head cellhd; /* it stores region information, and header information of rasters */
+
+ mapset = G_find_cell2(mapname, ""); /* checking if map exist */
+
+ if (mapset == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), mapname);
+
+ if (G_get_cellhd(mapname, mapset, &cellhd) < 0)
+ G_fatal_error(_("Unable to read file header of <%s>"), mapname);
+
+ if (mapname != in_elev) {
+ /* checking if the input maps type is CELL and region is valid */
+ if (G_raster_map_type(mapname, mapset) != CELL_TYPE)
+ G_fatal_error(_("<%s> is not of type CELL"), mapname);
+ if (window.ew_res != cellhd.ew_res || window.ns_res != cellhd.ns_res)
+ G_fatal_error(_("Region resolution and map %s resolution differs. Run g.region rast=%s to set proper resolution"),
+ mapname, mapname);
+ } /* end if map not elev */
+
+ if ((fd = G_open_cell_old(mapname, mapset)) < 0) /* file descriptor */
+ G_fatal_error(_("Unable to open raster map <%s>"), mapname);
+
+ if (G_get_cellhd(mapname, mapset, &cellhd) < 0)
+ G_fatal_error(_("Unable to read file header of <%s>"), mapname);
+
+ return fd;
+}
+
+int create_base_maps(void)
+{
+ int r, c;
+ int elev_type;
+ int in_dir_fd, in_stm_fd, in_dem_fd;
+ CELL *r_dirs, *r_streams;
+ FCELL *r_dem_f = NULL;
+ DCELL *r_dem_d = NULL;
+ CELL *r_dem_c = NULL;
+
+ in_dir_fd = open_raster(in_dirs);
+ in_stm_fd = open_raster(in_streams);
+
+ r_dirs = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ r_streams = (CELL *) G_malloc(sizeof(CELL) * ncols);
+
+ dirs = (CELL **) G_malloc(sizeof(CELL *) * nrows);
+ streams = (CELL **) G_malloc(sizeof(CELL *) * nrows);
+
+
+ if (extended) {
+ in_dem_fd = open_raster(in_elev);
+
+ elev_type = G_get_raster_map_type(in_dem_fd);
+ switch (elev_type) {
+ case CELL_TYPE:
+ r_dem_c = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ break;
+ case FCELL_TYPE:
+ r_dem_f = (FCELL *) G_malloc(sizeof(FCELL) * ncols);
+ break;
+ case DCELL_TYPE:
+ r_dem_d = (DCELL *) G_malloc(sizeof(DCELL) * ncols);
+ break;
+ }
+ elevation = (FCELL **) G_malloc(sizeof(FCELL *) * nrows);
+ }
+
+
+
+ G_message(_("Reading maps..."));
+
+ for (r = 0; r < nrows; ++r) {
+ G_percent(r, nrows, 2);
+
+ /* dirs & streams */
+ dirs[r] = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ streams[r] = (CELL *) G_malloc(sizeof(CELL) * ncols);
+
+ if (G_get_c_raster_row(in_dir_fd, r_dirs, r) < 0 ||
+ G_get_c_raster_row(in_stm_fd, r_streams, r) < 0) {
+ G_close_cell(in_dir_fd);
+ G_close_cell(in_stm_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"), r);
+ }
+
+if (extended) {
+ elevation[r] = (FCELL *) G_malloc(sizeof(FCELL) * ncols);
+
+ switch (elev_type) {
+ case CELL_TYPE: /* CELL */
+ if (G_get_c_raster_row(in_dem_fd, r_dem_c, r) < 0) {
+ G_close_cell(in_dem_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"),
+ r);
+ }
+ for (c = 0; c < ncols; ++c) {
+ if (G_is_c_null_value(&r_dem_c[c])) {
+ G_set_f_null_value(&elevation[r][c], 1);
+ }
+ else {
+ elevation[r][c] = (FCELL) r_dem_c[c];
+ }
+ } /* end for */
+ break; /* CELL */
+
+ case DCELL_TYPE: /* DCELL */
+ if (G_get_d_raster_row(in_dem_fd, r_dem_d, r) < 0) {
+ G_close_cell(in_dem_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"),
+ r);
+ }
+ for (c = 0; c < ncols; ++c) {
+ if (G_is_d_null_value(&r_dem_d[c])) {
+ G_set_f_null_value(&elevation[r][c], 1);
+ }
+ else {
+ elevation[r][c] = (FCELL) r_dem_d[c];
+
+ }
+ } /* end for */
+ break; /* DCELL */
+
+ case FCELL_TYPE: /* FCELL */
+ if (G_get_f_raster_row(in_dem_fd, elevation[r], r) < 0) {
+ G_close_cell(in_dem_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"),
+ r);
+ }
+ break; /* FCELL */
+
+ } /* end switch */
+ } /* end if elev */
+
+
+ for (c = 0; c < ncols; ++c) {
+ if (G_is_c_null_value(&r_dirs[c])) {
+ dirs[r][c] = 0;
+ }
+ else {
+ dirs[r][c] = r_dirs[c];
+ }
+
+ if (G_is_c_null_value(&r_streams[c])) {
+ streams[r][c] = 0;
+ }
+ else {
+ streams[r][c] = r_streams[c];
+ }
+ } /* end for c */
+
+ /* END dirs & streams & accums */
+
+ } /*end for r */
+
+ switch (elev_type) {
+ case CELL_TYPE:
+ G_free(r_dem_c);
+ break;
+ case DCELL_TYPE:
+ G_free(r_dem_d);
+ break;
+ }
+
+ G_free(r_streams);
+ G_free(r_dirs);
+ G_percent(r, nrows, 2);
+ G_close_cell(in_dir_fd);
+ G_close_cell(in_stm_fd);
+
+ return 0;
+} /* end create base maps */
+
+
+int stream_number(void)
+{
+ char *cur_mapset;
+ CELL c_min, c_max;
+ struct Range stream_range;
+ G_init_range(&stream_range);
+ cur_mapset = G_find_cell2(in_streams, "");
+ G_read_range(in_streams,cur_mapset,&stream_range);
+ G_get_range_min_max(&stream_range, &c_min, &c_max);
+ stream_num=c_max;
+ return 0;
+}
+
+
+int close_vector () {
+
+ int i,j;
+ int cur_order;
+ int index_cat=0;
+ float angle, tangent;
+ struct field_info *Fi;
+ dbDriver *driver;
+ dbHandle handle;
+ dbString table_name, db_sql, val_string;
+ char *cat_col_name="cat";
+ char buf[1000];
+
+ 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);
+ }
+
+ if (extended) /* calculate extended topology */
+ sprintf(buf,"create table %s (%s integer, \
+ c_stream integer, \
+ segment integer, \
+ c_order integer, \
+ direction double precision, \
+ length double precision, \
+ drop double precision, \
+ ncells integer, \
+ n_stream integer, \
+ n_tangent double precision\
+ )",
+ Fi->table,cat_col_name);
+ else /* only simple results */
+ sprintf(buf,"create table %s (%s integer, \
+ c_stream integer, \
+ c_order integer, \
+ direction double precision, \
+ length 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("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<stream_num+1;++i) {
+
+ for (j=0;j<seg_common[i].seg_num;++j) {
+ index_cat=seg_common[i].cats[j];
+
+ angle=(radians) ? seg_common[i].angles[j] : rad2deg(seg_common[i].angles[j]);
+ tangent=(radians) ?
+ s_streams[s_streams[i].next_stream].tangent_dir :
+ rad2deg(s_streams[s_streams[i].next_stream].tangent_dir);
+
+ cur_order=(ordering) ? cur_orders(s_streams[i].stream,ordering) : 0;
+
+ if(extended)
+ sprintf(buf,"insert into %s values (%d, %d, %d, %d, %f, %f, %f, %d, %d, %f)",
+ Fi->table,
+ index_cat,
+ s_streams[i].stream,
+ j,
+ cur_order,
+ angle,
+ seg_common[i].lengths[j],
+ seg_common[i].drops[j],
+ seg_common[i].cellnums[j],
+ s_streams[i].next_stream,
+ tangent
+ );
+ else
+ sprintf(buf,"insert into %s values (%d, %d, %d, %f, %f)",
+ Fi->table,
+ index_cat,
+ s_streams[i].stream,
+ cur_order,
+ angle,
+ seg_common[i].lengths[j]
+ );
+
+ 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);
+}
Property changes on: grass-addons/raster/r.fuzzy/io.c
___________________________________________________________________
Added: svn:executable
+
Added: grass-addons/raster/r.fuzzy/main.c
===================================================================
--- grass-addons/raster/r.fuzzy/main.c (rev 0)
+++ grass-addons/raster/r.fuzzy/main.c 2010-03-03 06:46:46 UTC (rev 41270)
@@ -0,0 +1,199 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.stream.angle
+ * AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl
+ *
+ * PURPOSE: Calculate streams orientation and angles between streams and its
+ * tributuaries. For stream direction it use algorithim to divide
+ * particular streams of same order into near-stright line
+ * segments.
+ *
+ *
+ *
+ * COPYRIGHT: (C) 2002,2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#define MAIN
+#include <grass/glocale.h>
+#include "global.h"
+
+/*
+ * main function *
+ */
+int main(int argc, char *argv[])
+{
+
+ struct GModule *module; /* GRASS module for parsing arguments */
+ struct Option *in_dir_opt,
+ *in_stm_opt, /* streams */
+ *in_elev_opt, /* streams */
+ *in_length,
+ *in_outlet,
+ *in_threshold,
+ *in_ordering,
+ *out_segments; /* options */
+
+ struct Flag *out_rad,
+ *out_ext;
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ module->keywords =
+ _("stream topology, route azimuth, route direction");
+ module->description =
+ _("Route azimuth, direction and relation to streams of higher order");
+
+ /*input option direction is reqired, acummulation is optional */
+
+ in_stm_opt = G_define_option();
+ in_stm_opt->key = "stream";
+ in_stm_opt->type = TYPE_STRING;
+ in_stm_opt->required = YES;
+ in_stm_opt->gisprompt = "old,cell,raster";
+ in_stm_opt->description =
+ "Name of stream mask input map (r.watershed or r.stream.extract)";
+
+ in_dir_opt = G_define_option(); /* input directon file */
+ in_dir_opt->key = "dirs";
+ in_dir_opt->type = TYPE_STRING;
+ in_dir_opt->required = YES;
+ in_dir_opt->gisprompt = "old,cell,raster";
+ in_dir_opt->description =
+ "Name of direction input map (r.watershed or r.stream.extract)";
+
+ in_elev_opt = G_define_option(); /* input directon file */
+ in_elev_opt->key = "elev";
+ in_elev_opt->type = TYPE_STRING;
+ in_elev_opt->required = NO;
+ in_elev_opt->gisprompt = "old,cell,raster";
+ in_elev_opt->answer = NULL;
+ in_elev_opt->description =
+ "Name of elevation input map";
+
+ in_ordering = G_define_option();
+ in_ordering->key = "order";
+ in_ordering->description ="Stream ordering method";
+ in_ordering->type = TYPE_STRING;
+ in_ordering->required = YES;
+ in_ordering->options ="none,hack,horton,strahler";
+ in_ordering->answer = "none";
+
+ in_length = G_define_option();
+ in_length->key = "length";
+ in_length->label = _("Search length to calculat direction");
+ in_length->description = _("Must be > 0");
+ in_length->required = YES;
+ in_length->answer = "15";
+ in_length->type = TYPE_INTEGER;
+
+ in_outlet = G_define_option();
+ in_outlet->key = "skip";
+ in_outlet->label = _("Skip segments shorter than");
+ in_outlet->description = _("Must be >= 0");
+ in_outlet->required = YES;
+ in_outlet->answer = "5";
+ in_outlet->type = TYPE_INTEGER;
+
+ in_threshold = G_define_option();
+ in_threshold->key = "treshold";
+ in_threshold->label = _("Max angle (degrees) beetwen stream segments to ");
+ in_threshold->description = _("Must be > 0");
+ in_threshold->required = YES;
+ in_threshold->answer = "150";
+ in_threshold->type = TYPE_DOUBLE;
+
+ out_segments = G_define_option();
+ out_segments->key = "seg_vector";
+ out_segments->label = _("Vector to store new network with segments");
+ out_segments->required = YES;
+ out_segments->answer = NULL;
+ out_segments->gisprompt = "new,vector,vector";
+
+ out_rad = G_define_flag();
+ out_rad->key = 'r';
+ out_rad->description = _("Output angles in radians (default: degrees)");
+
+ out_ext = G_define_flag();
+ out_ext->key = 'e';
+ out_ext->description = _("Extended topology (default: calculate direction only)");
+
+ if (G_parser(argc, argv)) /* parser */
+ exit(EXIT_FAILURE);
+
+ G_get_window(&window);
+
+ /* stores input options to variables */
+ in_dirs = in_dir_opt->answer;
+ in_streams = in_stm_opt->answer;
+ in_elev = in_elev_opt->answer;
+
+ out_vector=out_segments->answer;
+
+ seg_length = atoi(in_length->answer);
+ seg_treshold = atof(in_threshold->answer);
+ seg_outlet = atoi(in_outlet->answer);
+
+ radians = (out_rad->answer != 0);
+ extended = (out_ext->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_outlet < 0)
+ G_fatal_error("Segment's length must be >= 0");
+ if (extended && !in_elev)
+ G_fatal_error("For extended topology elevation map is required");
+
+ seg_treshold=deg2rad(seg_treshold);
+
+
+
+ /* stores output options to variables */
+
+ if (!strcmp(in_ordering->answer, "none"))
+ ordering = NONE;
+ else if (!strcmp(in_ordering->answer, "hack"))
+ ordering = HACK;
+ else if (!strcmp(in_ordering->answer, "horton"))
+ ordering = HORTON;
+ else if (!strcmp(in_ordering->answer, "strahler"))
+ ordering = STRAHLER;
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+
+ create_base_maps();
+ stream_number();
+
+ stack_max = stream_num; /* stack's size depends on number of streams */
+ G_begin_distance_calculations();
+ init_streams();
+ find_nodes();
+ do_cum_length();
+ strahler();
+ if (ordering==HORTON)
+ horton();
+ if (ordering==HACK)
+ hack();
+
+ if(extended) {
+ G_message("Calculate tangents at stream joins...");
+ calc_tangent(ordering);
+ }
+
+ G_message("Calculate segments...");
+ do_segments(seg_common,ordering);
+ if (extended)
+ add_missing_dirs(ordering);
+ close_vector();
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/raster/r.fuzzy/main.c
___________________________________________________________________
Added: svn:executable
+
Added: grass-addons/raster/r.fuzzy/order.c
===================================================================
--- grass-addons/raster/r.fuzzy/order.c (rev 0)
+++ grass-addons/raster/r.fuzzy/order.c 2010-03-03 06:46:46 UTC (rev 41270)
@@ -0,0 +1,258 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+int strahler(void)
+{
+
+ int i, j, done = 1;
+ int cur_stream, next_stream;
+ int max_strahler = 0, max_strahler_num;
+
+ G_message("Calculating Strahler's stream order ...");
+
+ for (j = 0; j < springs_num; ++j) { /* main loop on springs */
+
+ cur_stream = s_streams[springs[j]].stream;
+ do { /* we must go at least once, if stream is of first order and is outlet */
+ max_strahler_num = 1;
+ max_strahler = 0;
+ next_stream = s_streams[cur_stream].next_stream;
+
+ if (s_streams[cur_stream].trib_num == 0) { /* assign 1 for spring stream */
+ s_streams[cur_stream].strahler = 1;
+ cur_stream = next_stream;
+ done = 1;
+ }
+ else {
+ done = 1;
+
+ for (i = 0; i < s_streams[cur_stream].trib_num; ++i) { /* loop for determining strahler */
+ if (s_streams[s_streams[cur_stream].trib[i]].strahler < 0) {
+ done = 0;
+ break; /* strahler is not determined, break for loop */
+ }
+ else if (s_streams[s_streams[cur_stream].trib[i]].
+ strahler > max_strahler) {
+ max_strahler =
+ s_streams[s_streams[cur_stream].trib[i]].strahler;
+ max_strahler_num = 1;
+ }
+ else if (s_streams[s_streams[cur_stream].trib[i]].
+ strahler == max_strahler) {
+ ++max_strahler_num;
+ }
+ } /* end determining strahler */
+
+ if (done == 1) {
+ s_streams[cur_stream].strahler =
+ (max_strahler_num >
+ 1) ? ++max_strahler : max_strahler;
+ cur_stream = next_stream; /* if next_stream<0 we in outlet stream */
+ }
+
+ }
+ } while (done && next_stream > 0);
+ } /* end for of main loop */
+ return 0;
+} /* end strahler */
+
+int shreeve(void)
+{
+
+ int i, j, done = 1;
+ int cur_stream, next_stream;
+ int max_shreeve = 0;
+
+ G_message("Calculating Shreeve's stream magnitude ...");
+
+ for (j = 0; j < springs_num; ++j) { /* main loop on springs */
+
+ cur_stream = s_streams[springs[j]].stream;
+ do { /* we must go at least once, if stream is of first order and is outlet */
+
+ max_shreeve = 0;
+ next_stream = s_streams[cur_stream].next_stream;
+
+ if (s_streams[cur_stream].trib_num == 0) { /* assign 1 for spring stream */
+
+ s_streams[cur_stream].shreeve = 1;
+ cur_stream = next_stream;
+ done = 1;
+
+ }
+ else {
+ done = 1;
+
+ for (i = 0; i < s_streams[cur_stream].trib_num; ++i) { /* loop for determining strahler */
+ if (s_streams[s_streams[cur_stream].trib[i]].shreeve < 0) {
+ done = 0;
+ break; /* shreeve is not determined, break for loop */
+ }
+ else {
+ max_shreeve +=
+ s_streams[s_streams[cur_stream].trib[i]].shreeve;
+ }
+ } /* end determining strahler */
+
+ if (done == 1) {
+ s_streams[cur_stream].shreeve = max_shreeve;
+ cur_stream = next_stream; /* if next_stream<0 we in outlet stream */
+ }
+ }
+
+ } while (done && next_stream > 0);
+ } /* end main loop */
+ return 0;
+} /* end shreeve */
+
+int horton(void)
+{
+
+ int *stack;
+ int top, i, j;
+ int cur_stream, cur_horton;
+ int max_strahler;
+ double max_accum;
+ int up_stream = 0;
+
+ G_message("Calculating Hortons's stream order ...");
+ stack = (int *)G_malloc(stack_max * sizeof(int));
+
+ for (j = 0; j < outlets_num; ++j) {
+ cur_stream = s_streams[outlets[j]].stream; /* outlet: init */
+ cur_horton = s_streams[cur_stream].strahler;
+ stack[0] = 0;
+ stack[1] = cur_stream;
+ top = 1;
+
+ do { /* on every stream */
+ max_strahler = 0;
+ max_accum = 0;
+
+ if (s_streams[cur_stream].trib_num == 0) { /* spring: go back on stack */
+
+ s_streams[cur_stream].horton = cur_horton;
+ cur_stream = stack[--top];
+
+ }
+ else if (s_streams[cur_stream].trib_num > 1) { /* node */
+
+ up_stream = 0; /* calculating up_stream */
+ for (i = 0; i < s_streams[cur_stream].trib_num; ++i) {
+ if (s_streams[s_streams[cur_stream].trib[i]].horton < 0) {
+
+ if (s_streams[s_streams[cur_stream].trib[i]].
+ strahler > max_strahler) {
+ max_strahler =
+ s_streams[s_streams[cur_stream].trib[i]].
+ strahler;
+ max_accum =
+ s_streams[s_streams[cur_stream].trib[i]].
+ accum;
+ up_stream = s_streams[cur_stream].trib[i];
+
+ }
+ else if (s_streams[s_streams[cur_stream].trib[i]].
+ strahler == max_strahler) {
+
+ if (s_streams[s_streams[cur_stream].trib[i]].
+ accum > max_accum) {
+ max_accum =
+ s_streams[s_streams[cur_stream].trib[i]].
+ accum;
+ up_stream = s_streams[cur_stream].trib[i];
+ }
+ }
+ }
+ } /* end determining up_stream */
+
+ if (up_stream) { /* at least one branch is not assigned */
+ if (s_streams[cur_stream].horton < 0) {
+ s_streams[cur_stream].horton = cur_horton;
+ }
+ else {
+ cur_horton = s_streams[up_stream].strahler;
+ }
+ cur_stream = up_stream;
+ stack[++top] = cur_stream;
+
+ }
+ else { /* all asigned, go downstream */
+ cur_stream = stack[--top];
+
+ } /* end up_stream */
+ } /* end spring/node */
+ } while (cur_stream);
+ } /* end for outlets */
+ G_free(stack);
+ return 0;
+}
+
+int hack(void)
+{
+
+ int *stack;
+ int top, i, j;
+ int cur_stream, cur_hack;
+ double max_accum;
+ int up_stream = 0;
+
+ G_message("Calculating Hack's main streams ...");
+ stack = (int *)G_malloc(stack_max * sizeof(int));
+
+ for (j = 0; j < outlets_num; ++j) {
+ cur_stream = s_streams[outlets[j]].stream; /* outlet: init */
+ cur_hack = 1;
+ stack[0] = 0;
+ stack[1] = cur_stream;
+ top = 1;
+
+ do {
+ max_accum = 0;
+
+ if (s_streams[cur_stream].trib_num == 0) { /* spring: go back on stack */
+
+ s_streams[cur_stream].hack = cur_hack;
+ cur_stream = stack[--top];
+
+ }
+ else if (s_streams[cur_stream].trib_num > 1) { /* node */
+ up_stream = 0; /* calculating up_stream */
+
+ for (i = 0; i < s_streams[cur_stream].trib_num; ++i) { /* determining upstream */
+ if (s_streams[s_streams[cur_stream].trib[i]].hack < 0) {
+ if (s_streams[s_streams[cur_stream].trib[i]].accum >
+ max_accum) {
+ max_accum =
+ s_streams[s_streams[cur_stream].trib[i]].
+ accum;
+ up_stream = s_streams[cur_stream].trib[i];
+ }
+ }
+ } /* end determining up_stream */
+
+ if (up_stream) { /* at least one branch is not assigned */
+ if (s_streams[cur_stream].hack < 0) {
+ s_streams[cur_stream].hack = cur_hack;
+ }
+ else {
+ cur_hack = s_streams[cur_stream].hack;
+ ++cur_hack;
+ }
+ cur_stream = up_stream;
+ stack[++top] = cur_stream;
+
+ }
+ else { /* all asigned, go downstream */
+
+ cur_stream = stack[--top];
+
+ } /* end up_stream */
+ } /* end spring/node */
+ } while (cur_stream);
+ } /* end for outlets */
+ G_free(stack);
+ return 0;
+}
Property changes on: grass-addons/raster/r.fuzzy/order.c
___________________________________________________________________
Added: svn:executable
+
Added: grass-addons/raster/r.fuzzy/seg_line.c
===================================================================
--- grass-addons/raster/r.fuzzy/seg_line.c (rev 0)
+++ grass-addons/raster/r.fuzzy/seg_line.c 2010-03-03 06:46:46 UTC (rev 41270)
@@ -0,0 +1,368 @@
+#include <grass/glocale.h>
+#include "global.h"
+
+/* the idea of algorithm
+ *
+ */
+
+int do_segments(SEGMENTS *segments, int ordering) {
+
+ 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 i,j,k,s;
+ int done, in_loop; /*flags */
+ int r,c;
+ int init, out, cur_stream, next_stream, cell_num;
+ int strahler_springs;
+ int cur_order, next_order;
+ int init_stream,out_stream;
+ DIRCELLS* cells;
+ int tmp=0;
+
+ /* variables for line segmentation */
+ int small_seg_length;
+ int cell_down, cell_up, small_cell_down, small_cell_up;
+ float dir_down, dir_up, small_dir_down, small_dir_up;
+ int out_segs, cat_num;
+ int local_minimum_point;
+ float local_minimum=PI;
+ float local_minimum_tunning=PI;
+ int seg_nums=0;
+ int prev_index=0;
+ int cells_in_segment=1;
+ int seg_cats=0;
+
+ double north_offset, west_offset, ns_res, ew_res;
+
+ ns_res = window.ns_res;
+ ew_res = window.ew_res;
+ north_offset = window.north - 0.5 * ns_res;
+ west_offset = window.west + 0.5 * ew_res;
+
+ /* the rest of vector code is in function close_vect in io.c file */
+
+ Segments = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ if (Vect_open_new(&Out, out_vector, 0)<0)
+ G_fatal_error(_("Unable to create new vector map <%s>"), out_vector);
+
+ /* global variables:
+ * seg_length
+ * seg_outlet
+ * seg_treshold
+ */
+
+ /* for STRAHLER ordering only*/
+ if (ordering==STRAHLER) {
+ for (i=1;i<stream_num+1;++i) {
+ if (s_streams[i].strahler>1) {
+ done=1;
+ for (j=0;j<s_streams[i].trib_num;++j) {
+ if (s_streams[s_streams[i].trib[j]].strahler==s_streams[i].strahler) {
+ done=0;
+ break;
+ }
+ }
+ if(done)
+ springs[springs_num++]=s_streams[i].stream;
+ }
+ }
+ }
+
+ /* for NONE ordering only*/
+ if (ordering==NONE) {
+ for (i=1;i<stream_num+1;++i) {
+ if (s_streams[i].strahler>1)
+ springs[springs_num++]=s_streams[i].stream;
+ }
+ }
+
+ /* end for strahler and none ordering only*/
+
+ for(s=0;s<springs_num;++s) {
+
+ init=springs[s];
+ cur_stream=init;
+ next_stream=s_streams[cur_stream].next_stream;
+ cell_num=s_streams[cur_stream].cell_num;
+
+ /* find outlet stream */
+ while (cur_orders(cur_stream,ordering)==cur_orders(next_stream,ordering)) {
+ cur_stream = next_stream;
+ next_stream = s_streams[cur_stream].next_stream;
+ cell_num += s_streams[cur_stream].cell_num;
+ }
+
+ /* reset */
+
+ out=cur_stream;
+ //s_streams[init].out_stream=out;
+ cur_stream=init;
+ next_stream=s_streams[cur_stream].next_stream;
+
+
+ while (cur_orders(cur_stream,ordering)==cur_orders(next_stream,ordering)) {
+ segments[cur_stream].init_stream=init;
+ segments[cur_stream].out_stream=out;
+ cur_stream = next_stream;
+ next_stream = s_streams[cur_stream].next_stream;
+ }
+
+ segments[cur_stream].init_stream=init;
+ segments[cur_stream].out_stream=cur_stream;
+ segments[cur_stream].cell_num=cell_num;
+
+
+ r=s_streams[cur_stream].out_r;
+ c=s_streams[cur_stream].out_c;
+
+ cells=(DIRCELLS *) G_malloc((cell_num+1) * sizeof(DIRCELLS));
+
+
+ for (j=0; j<cell_num; ++j) {
+
+ cells[j].r=r;
+ cells[j].c=c;
+ cells[j].dir_diff=-1;
+ cells[j].small_dir_diff=-1;
+ cells[j].candidate=0;
+ cells[j].tunning=0;
+ cells[j].decision=0;
+
+
+ for (i = 1; i < 9; ++i) {
+ if (r + nextr[i] < 0 || r + nextr[i] > (nrows - 1)
+ || c + nextc[i] < 0 || c + nextc[i] > (ncols - 1))
+ continue;
+
+ k = (i + 4) > 8 ? i - 4 : i + 4;
+
+ if (streams[r + nextr[i]][c + nextc[i]] > 0 &&
+ dirs[r + nextr[i]][c + nextc[i]] == k) {
+
+ cur_order=cur_orders(streams[r][c],ordering);
+ next_order=cur_orders(streams[r + nextr[i]][c + nextc[i]],ordering);
+
+ if (next_order==cur_order) {
+
+ r+=nextr[i];
+ c+=nextc[i];
+ break;
+ }
+ }
+ }
+ } /* end for j, cells initialized! */
+
+ if(cell_num<seg_outlet+1) {
+
+ segments[cur_stream].dir_full=
+ calc_dir(cells[cell_num-1].r, cells[cell_num-1].c,
+ cells[0].r,cells[0].c);
+
+ } else {
+
+ segments[cur_stream].dir_init=
+ calc_dir(cells[seg_outlet].r, cells[seg_outlet].c,
+ cells[0].r,cells[0].c);
+
+ segments[cur_stream].dir_full=
+ calc_dir(cells[cell_num-1].r, cells[cell_num-1].c,
+ cells[0].r,cells[0].c);
+
+ segments[cur_stream].dir_middle=
+ calc_dir(cells[(cell_num/2)+(seg_outlet/2)].r,
+ cells[(cell_num/2)+(seg_outlet/2)].c,
+ cells[(cell_num/2)-(seg_outlet/2)].r,
+ cells[(cell_num/2)-(seg_outlet/2)].c);
+
+ segments[cur_stream].dir_last=
+ calc_dir(cells[cell_num-1].r, cells[cell_num-1].c,
+ cells[(cell_num-1)-seg_outlet].r,cells[(cell_num-1)-seg_outlet].c);
+
+ }
+
+ small_seg_length = seg_length / 3;
+
+ for (i=(seg_outlet+1); i<(cell_num-seg_outlet);++i) {
+
+ cell_down = (i<seg_length) ? i : seg_length;
+ cell_up = (i>((cell_num-1)-seg_length)) ? ((cell_num-1)-i) : seg_length;
+
+ dir_down =
+ calc_dir(cells[i].r, cells[i].c,
+ cells[i-cell_down].r,cells[i-cell_down].c);
+
+ dir_up =
+ calc_dir(cells[i+cell_up].r,cells[i+cell_up].c,
+ cells[i].r,cells[i].c);
+
+ dir_up = (dir_up >= PI) ?
+ dir_up - PI : dir_up + PI; /* dir_up: direction upstream */
+
+ cells[i].dir_diff =
+ fabs(dir_down-dir_up) > PI ?
+ (PI * 2)-fabs(dir_down-dir_up) : fabs(dir_down-dir_up);
+
+ small_cell_down =
+ (i<small_seg_length) ? i : small_seg_length;
+
+ small_cell_up = (i>((cell_num-1)-small_seg_length)) ?
+ ((cell_num-1)-i) : small_seg_length;
+
+ small_dir_down =
+ calc_dir(cells[i].r, cells[i].c,
+ cells[i-small_cell_down].r,cells[i-small_cell_down].c);
+
+ small_dir_up =
+ calc_dir(cells[i+small_cell_up].r,cells[i+small_cell_up].c,
+ cells[i].r,cells[i].c);
+
+ small_dir_up = (small_dir_up >= PI) ?
+ small_dir_up - PI :
+ small_dir_up + PI; /* small_dir_up: direction upstream */
+
+ cells[i].small_dir_diff=
+ fabs(small_dir_down-small_dir_up) > PI ?
+ (PI * 2)-fabs(small_dir_down-small_dir_up) :
+ fabs(small_dir_down-small_dir_up);
+
+ cells[i].candidate = (cells[i].dir_diff<seg_treshold) ? 1 : 0;
+ cells[i].tunning = (cells[i].small_dir_diff<(seg_treshold*0.8)) ? 1 : 0;
+
+ } /* end for i */
+
+ /* decision system: prototype */
+
+
+ /* fine tunning */
+ local_minimum=PI;
+ cat_num=0;
+ in_loop=0;
+
+ for (i=0;i<cell_num;++i) {
+ if(cells[i].candidate) {
+
+ out_segs=0;
+ if (local_minimum>cells[i].small_dir_diff) {
+ local_minimum=cells[i].small_dir_diff;
+ local_minimum_point=i;
+ in_loop=1;
+ } /* end local minimum */
+
+ } else if (!cells[i].candidate && in_loop) {
+
+ out_segs++;
+ if (out_segs==(seg_length/5)) {
+ cells[local_minimum_point].decision=1;
+ local_minimum=PI;
+ in_loop=0;
+ }
+ }
+ } /* end for: fine tunning */
+
+ /* cleaning */
+ for (i=0,out_segs=0;i<cell_num;++i,out_segs++) {
+
+ if(cells[i].decision) {
+
+ if(out_segs<seg_outlet && i>seg_outlet) { /* was: seg_length/5 */
+ cells[i].decision=0;
+ i=local_minimum_point;
+ } else {
+ local_minimum_point=i;
+ }
+ out_segs=0;
+ }
+ }
+
+ seg_nums=0;
+ prev_index=0;
+ cells_in_segment=1;
+
+ for (i=0;i<cell_num;++i) {
+ if (cells[i].decision==1 || i==(cell_num-1))
+ seg_nums++;
+ }
+
+ segments[cur_stream].seg_num=seg_nums;
+ segments[cur_stream].angles=(float*)G_malloc(seg_nums*sizeof(float));
+ segments[cur_stream].lengths=(float*)G_malloc(seg_nums*sizeof(float));
+ segments[cur_stream].drops=(float*)G_malloc(seg_nums*sizeof(float));
+ segments[cur_stream].cellnums=(int*)G_malloc(seg_nums*sizeof(int));
+ segments[cur_stream].cats=(int*)G_malloc(seg_nums*sizeof(int));
+
+ seg_nums=0;
+
+ /* we go from segment outlet upstream */
+ for (i=0;i<cell_num;++i) {
+ if (cells[i].decision==1 || i==(cell_num-1)) {
+
+ segments[cur_stream].angles[seg_nums]=
+ calc_dir(cells[i].r,cells[i].c,
+ cells[prev_index].r,cells[prev_index].c);
+
+ segments[cur_stream].lengths[seg_nums]=
+ calc_length(cells[i].r,cells[i].c,
+ cells[prev_index].r,cells[prev_index].c);
+
+ if(extended) {
+ segments[cur_stream].drops[seg_nums]=
+ calc_drop(cells[i].r,cells[i].c,
+ cells[prev_index].r,cells[prev_index].c);
+
+ segments[cur_stream].cellnums[seg_nums]=
+ cells_in_segment;
+ }
+ segments[cur_stream].cats[seg_nums++]=
+ ++seg_cats;
+
+ cells_in_segment=0;
+
+ if (i<(cell_num-1))
+ prev_index=i+1;
+ }
+ cells_in_segment++;
+
+ }
+
+ Vect_reset_line(Segments);
+ Vect_reset_cats(Cats);
+
+ for(i=cell_num-1;i>-1;--i) {
+ Vect_append_point(Segments, west_offset + cells[i].c * ew_res,
+ north_offset - cells[i].r * ns_res, 0);
+
+ if(cells[i].decision==1 || i==0) {/* start new line */
+
+ if(dirs[cells[i].r][cells[i].c]>0) /*add point to create network */
+
+ Vect_append_point(Segments, west_offset +
+ (cells[i].c+nextc[dirs[cells[i].r][cells[i].c]]) * ew_res,
+ north_offset -
+ (cells[i].r+nextr[dirs[cells[i].r][cells[i].c]]) * ns_res, 0);
+
+ Vect_cat_set(Cats, 1, segments[cur_stream].cats[--seg_nums]);
+
+ Vect_write_line(&Out, GV_LINE, Segments, Cats);
+ Vect_reset_line(Segments);
+ Vect_reset_cats(Cats);
+ }
+ }
+
+ G_free(cells);
+
+ }
+
+return 0;
+}
+
+
+
+
+
+
+
+
+
Property changes on: grass-addons/raster/r.fuzzy/seg_line.c
___________________________________________________________________
Added: svn:executable
+
Added: grass-addons/raster/r.fuzzy/tangent.c
===================================================================
--- grass-addons/raster/r.fuzzy/tangent.c (rev 0)
+++ grass-addons/raster/r.fuzzy/tangent.c 2010-03-03 06:46:46 UTC (rev 41270)
@@ -0,0 +1,192 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+
+
+/* avg_direction calculate direction of stream in the confluence point as
+ * avaraged direction of current segment and previous segment of the same order
+ * the number of cell used for this can be determinde by user by parameter step_num */
+
+int calc_tangent(int ordering) {
+ 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 cur_stream;
+ int prev_stream=0;
+ int cur_order, prev_order;
+ int r, c, d;
+ int rn, cn; /* row, col cur stream */
+ int next_rn, next_cn;
+ int rp, cp; /* row, col prev stream */
+ int next_rp, next_cp;
+ int init_stream;
+ int i, j, k;
+ int step_num;
+ int *tribstack;
+ int top=1;
+ int done;
+
+ int out_dir;
+
+ tribstack=(int *)G_malloc(stack_max * sizeof(int));
+
+ for (j=0;j<outlets_num;++j) {
+ cur_stream = s_streams[outlets[j]].stream;
+ tribstack[0]=0;
+ top=1;
+
+ do {
+ /* prev streams and tribs */
+ prev_stream=0;
+ cur_order=cur_orders(cur_stream,ordering);
+
+ for (i=0;i<s_streams[cur_stream].trib_num;++i) {
+ prev_order=cur_orders(s_streams[cur_stream].trib[i],ordering);
+ if (prev_order==cur_order) {
+ prev_stream=s_streams[cur_stream].trib[i]; /* is continuation */
+ } else {
+ tribstack[top++]=s_streams[cur_stream].trib[i];
+
+ } /* end else */
+
+ } /* for i */
+
+ if (!prev_stream) {
+ cur_stream=tribstack[--top];
+ if(cur_stream) /* what to do with the stack */
+ continue;
+ else
+ break;
+ }
+
+ rn = r = s_streams[cur_stream].init_r;
+ cn = c = s_streams[cur_stream].init_c;
+
+ for (i = 1; i < 9; ++i) {
+ if (r + nextr[i] < 0 || r + nextr[i] > (nrows - 1)
+ || c + nextc[i] < 0 || c + nextc[i] > (ncols - 1))
+ continue;
+ if (streams[r + nextr[i]][c + nextc[i]] == prev_stream) {
+ rp=r + nextr[i];
+ cp=c + nextc[i];
+ }
+ } /* for i */
+
+ /* calculate direction for cur_stream */
+ step_num=seg_length;
+
+ while(step_num) {
+
+ step_num--;
+
+ //downstream
+ d=dirs[rn][cn];
+ next_rn=rn+nextr[d];
+ next_cn=cn+nextc[d];
+
+ if (d<1|| //no more cells
+ next_rn < 0 || next_rn > (nrows - 1) ||
+ next_cn < 0 || next_cn > (ncols - 1))
+
+ break; //leave internal while
+
+ rn=next_rn;
+ cn=next_cn;
+
+ //upstream
+ prev_order=0;
+ done=0;
+
+ for (i = 1; i < 9; ++i) {
+ if (rp + nextr[i] < 0 || rp + nextr[i] > (nrows - 1)
+ || cp + nextc[i] < 0 || cp + nextc[i] > (ncols - 1))
+ continue;
+
+ k = (i + 4) > 8 ? i - 4 : i + 4;
+
+ if (streams[rp + nextr[i]][cp + nextc[i]] > 0 &&
+ dirs[rp + nextr[i]][cp + nextc[i]] == k) {
+ prev_order=cur_orders(streams[rp + nextr[i]][cp + nextc[i]],ordering);
+
+ if (prev_order==cur_order) {
+ rp+=nextr[i];
+ cp+=nextc[i];
+ done=1;
+ break; /* leave for */
+ }
+ }
+ } /* end for i */
+
+ if (!done)
+ break; /* leave internal while */
+ } /* end internal while */
+
+ s_streams[cur_stream].tangent_dir = calc_dir(rp,cp,rn,cn);
+
+ cur_stream=prev_stream;
+ } while (top);
+
+ } /* for j */
+G_free(tribstack);
+return 0;
+}
+
+/* for strahler and none ordering
+ * add missing dirs for streams of higher order
+ */
+
+int add_missing_dirs(int ordering) {
+ int i, j;
+ int out;
+ int seg_num;
+ float last_angle, pre_last_angle;
+ float last_angle_diff;
+ int last_length, pre_last_length;
+
+ for (j = 1; j <= stream_num; ++j) {
+
+ if (s_streams[j].strahler==1 || seg_common[j].init_stream!=j)
+ continue;
+
+ out=seg_common[j].out_stream;
+ seg_num=seg_common[out].seg_num;
+
+ if(seg_num==1) {
+ s_streams[j].tangent_dir=seg_common[out].angles[0];
+ continue;
+ }
+
+ last_angle=seg_common[out].angles[seg_num-1];
+ last_length=seg_common[out].lengths[seg_num-1];
+
+ if(last_length>seg_outlet) {
+ s_streams[j].tangent_dir=seg_common[out].angles[seg_num-1];
+ continue;
+ }
+
+ pre_last_angle=seg_common[out].angles[seg_num-2];
+
+ last_angle_diff=(abs(pre_last_angle-last_angle)>PI) ?
+ 2*PI-abs(pre_last_angle-last_angle) : abs(pre_last_angle-last_angle);
+
+ if(last_angle_diff>PI-seg_outlet) {
+ s_streams[j].tangent_dir=seg_common[out].angles[seg_num-2];
+ } else {
+ s_streams[j].tangent_dir=
+ calc_dir(s_streams[j].init_r,
+ s_streams[j].init_c,
+ s_streams[s_streams[j].next_stream].out_r,
+ s_streams[s_streams[j].next_stream].out_c);
+
+
+ }
+ /* calc stream direction using coordinates
+ * of two last segments of current stream segment */
+
+ } /* for j */
+return 0;
+}
+
Property changes on: grass-addons/raster/r.fuzzy/tangent.c
___________________________________________________________________
Added: svn:executable
+
More information about the grass-commit
mailing list