Merge r.watershed2 from Markus Metz into r.watershed.

(merge from devbr6 https://trac.osgeo.org/grass/changeset/34645)

fast version of r.watershed:
Implemented is a new sorting algorithm for both the ram and the 
segmented version, resulting in significant speed improvements,
with identical results.

The new sorting algorithm is implemented for both the ram and the
segmented version.
The new segmented version is still much slower than the new ram 
version, but the new segmented version is already much faster than 
the original ram version. A new option is available for the segmented 
mode specifying the maximum amount of memory to be used in MB.


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

raster/r.watershed/front/main.c

@@ -6,7 +6,7 @@
  *               Brad Douglas <rez touchofmadness.com>,
  *		 Hamish Bowman <hamish_b yahoo.com>
  * PURPOSE:      Watershed determination
- * COPYRIGHT:    (C) 1999-2006 by the GRASS Development Team
+ * COPYRIGHT:    (C) 1999-2008 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
@@ -40,6 +40,7 @@ int main(int argc, char *argv[])
     struct Option *opt13;
     struct Option *opt14;
     struct Option *opt15;
+    struct Option *opt16;
     struct Flag *flag_flow;
     struct Flag *flag_seg;
     struct GModule *module;
@@ -177,6 +178,13 @@ int main(int argc, char *argv[])
     opt15->gisprompt = "new,cell,raster";
     opt15->guisection = _("Output_options");
 
+    opt16 = G_define_option() ;
+    opt16->key         = "memory";
+    opt16->type        = TYPE_INTEGER;
+    opt16->required    = NO;
+    opt16->answer      = "300"; /* 300MB default value, please keep in sync with r.terraflow */
+    opt16->description = _("Maximum memory to be used with -m flag (in MB)");
+
     flag_flow = G_define_flag();
     flag_flow->key = '4';
     flag_flow->description =
@@ -336,6 +344,13 @@ int main(int argc, char *argv[])
 	strcat(command, "\"");
     }
 
+    if (flag_seg->answer && opt16->answer) {
+	strcat(command, " mb=");
+	strcat(command, "\"");
+	strcat(command, opt16->answer);
+	strcat(command, "\"");
+    }
+
     G_debug(1, "Mode: %s", flag_seg->answer ? "Segmented" : "All in RAM");
     G_debug(1, "Running: %s", command);
 

raster/r.watershed/front/r.watershed.html

@@ -419,8 +419,13 @@ Food Production In Arid Zones</b> (Tucson, AZ, 12-16 Oct. 1987).
 </em>
 
 
-<h2>AUTHOR</h2>
+<h2>AUTHORS</h2>
 
+Original version:
 Charles Ehlschlaeger, U.S. Army Construction Engineering Research Laboratory
+<br>
+Faster sorting algorithm:
+Markus Metz &lt;markus.metz.giswork at gmail.com&gt;
+
 <p>
 <i>Last changed: $Date$</i>

raster/r.watershed/ram/Gwater.h

@@ -40,6 +40,7 @@ POINT {
 
 extern struct Cell_head window;
 
+extern int *heap_index, heap_size;
 extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 extern SHORT nrows, ncols;
 extern double half_res, diag, max_length, dep_slope;
@@ -85,6 +86,7 @@ CELL def_basin(int, int, CELL, double, CELL);
 /* do_astar.c */
 int do_astar(void);
 int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int drop_pt(void);
 double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
 int replace(SHORT, SHORT, SHORT, SHORT);
 

raster/r.watershed/ram/do_astar.c

@@ -1,7 +1,9 @@
 #include "Gwater.h"
+#include "do_astar.h"
 #include <grass/gis.h>
 #include <grass/glocale.h>
 
+int sift_up(int, CELL);
 
 int do_astar(void)
 {
@@ -10,49 +12,82 @@ int do_astar(void)
     SHORT upr, upc, r, c, ct_dir;
     CELL alt_val, alt_up, asp_up, wat_val;
     CELL in_val, drain_val;
+    int index_doer, index_up;
+
 
     G_message(_("SECTION 2: A * Search."));
 
     count = 0;
+    first_astar = heap_index[1];
+
+    /* A * Search: search uphill, get downhill path */
     while (first_astar != -1) {
 	G_percent(count++, do_points, 1);
-	doer = first_astar;
+
+	/* start with point with lowest elevation, in case of equal elevation
+	 * of following points, oldest point = point added earliest */
+	/* old routine: astar_pts[first_astar] (doer = first_astar) */
+	/* new routine: astar_pts[heap_index[1]] */
+
+	doer = heap_index[1];
+
 	point = &(astar_pts[doer]);
-	first_astar = point->nxt;
+
+	/* drop astar_pts[doer] from heap */
+	/* necessary to protect the current point (doer) from modification */
+	/* equivalent to first_astar = point->next in old code */
+	drop_pt();
+
+	/* can go, dragged on from old code */
+	first_astar = heap_index[1];
+
+	/* downhill path for flow accumulation is set here */
+	/* this path determines the order for flow accumulation calculation */
 	point->nxt = first_cum;
+	first_cum = doer;
+
+	r = point->r;
+	c = point->c;
 
-	r = astar_pts[doer].r = point->r;
-	c = astar_pts[doer].c = point->c;
 	G_debug(3, "R:%2d C:%2d, ", r, c);
 
-	astar_pts[doer].downr = point->downr;
-	astar_pts[doer].downc = point->downc;
-	astar_pts[doer].nxt = point->nxt;
-	first_cum = doer;
+	index_doer = SEG_INDEX(alt_seg, r, c);
+	alt_val = alt[index_doer];
+	wat_val = wat[index_doer];
+
 	FLAG_SET(worked, r, c);
-	alt_val = alt[SEG_INDEX(alt_seg, r, c)];
+
+	/* check all neighbours */
 	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+	    /* get r, c (upr, upc) for this neighbour */
 	    upr = r + nextdr[ct_dir];
 	    upc = c + nextdc[ct_dir];
+	    index_up = SEG_INDEX(alt_seg, upr, upc);
+	    /* check that r, c are within region */
 	    if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
+		/* check if neighbour is in the list */
+		/* if not, add as new point */
 		in_val = FLAG_GET(in_list, upr, upc);
 		if (in_val == 0) {
-		    alt_up = alt[SEG_INDEX(alt_seg, upr, upc)];
+		    alt_up = alt[index_up];
+		    /* flow direction is set here */
 		    add_pt(upr, upc, r, c, alt_up, alt_val);
 		    drain_val = drain[upr - r + 1][upc - c + 1];
-		    asp[SEG_INDEX(asp_seg, upr, upc)] = drain_val;
+		    asp[index_up] = drain_val;
+
 		}
 		else {
+		    /* check if neighbour has not been worked on,
+		     * update values for asp and wat */
 		    in_val = FLAG_GET(worked, upr, upc);
 		    if (in_val == 0) {
-			asp_up = asp[SEG_INDEX(asp_seg, upr, upc)];
+			asp_up = asp[index_up];
 			if (asp_up < -1) {
-			    asp[SEG_INDEX(asp_seg, upr, upc)] =
-				drain[upr - r + 1][upc - c + 1];
-			    wat_val = wat[SEG_INDEX(wat_seg, r, c)];
+			    asp[index_up] = drain[upr - r + 1][upc - c + 1];
+
 			    if (wat_val > 0)
-				wat[SEG_INDEX(wat_seg, r, c)] = -wat_val;
-			    alt_up = alt[SEG_INDEX(alt_seg, upr, upc)];
+				wat[index_doer] = -wat_val;
+
 			    replace(upr, upc, r, c);	/* alt_up used to be */
 			}
 		    }
@@ -63,65 +98,160 @@ int do_astar(void)
     G_percent(count, do_points, 3);	/* finish it */
     flag_destroy(worked);
     flag_destroy(in_list);
+    G_free(heap_index);
 
     return 0;
 }
 
+/* new add point routine for min heap */
 int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
 {
-    int p;
-    CELL check_ele;
-    POINT point;
 
     FLAG_SET(in_list, r, c);
-    if (first_astar == -1) {
-	astar_pts[nxt_avail_pt].r = r;
-	astar_pts[nxt_avail_pt].c = c;
-	astar_pts[nxt_avail_pt].downr = downr;
-	astar_pts[nxt_avail_pt].downc = downc;
-	astar_pts[nxt_avail_pt].nxt = -1;
-	first_astar = nxt_avail_pt;
-	nxt_avail_pt++;
+
+    /* add point to next free position */
+
+    heap_size++;
+
+    if (heap_size > do_points)
+	G_fatal_error(_("heapsize too large"));
+
+    heap_index[heap_size] = nxt_avail_pt;
+
+    astar_pts[nxt_avail_pt].r = r;
+    astar_pts[nxt_avail_pt].c = c;
+    astar_pts[nxt_avail_pt].downr = downr;
+    astar_pts[nxt_avail_pt].downc = downc;
+
+    nxt_avail_pt++;
+
+    /* sift up: move new point towards top of heap */
+
+    sift_up(heap_size, ele);
+
+    return 0;
+}
+
+/* new drop point routine for min heap */
+int drop_pt(void)
+{
+    int child, childr, parent;
+    CELL ele, eler;
+    int i;
+
+    if (heap_size == 1) {
+	heap_index[1] = -1;
+	heap_size = 0;
 	return 0;
     }
-    p = first_astar;;
-    while (1) {
-	point.r = astar_pts[p].r;
-	point.c = astar_pts[p].c;
-	check_ele = alt[SEG_INDEX(alt_seg, point.r, point.c)];
-	if (check_ele > ele) {
-	    point.downr = astar_pts[p].downr;
-	    point.downc = astar_pts[p].downc;
-	    point.nxt = astar_pts[p].nxt;
-	    astar_pts[p].r = r;
-	    astar_pts[p].c = c;
-	    astar_pts[p].downr = downr;
-	    astar_pts[p].downc = downc;
-	    astar_pts[p].nxt = nxt_avail_pt;
-	    astar_pts[nxt_avail_pt].r = point.r;
-	    astar_pts[nxt_avail_pt].c = point.c;
-	    astar_pts[nxt_avail_pt].downr = point.downr;
-	    astar_pts[nxt_avail_pt].downc = point.downc;
-	    astar_pts[nxt_avail_pt].nxt = point.nxt;
-	    nxt_avail_pt++;
-	    return 0;
+
+    /* start with root */
+    parent = 1;
+
+    /* sift down: move hole back towards bottom of heap */
+    /* sift-down routine customised for A * Search logic */
+
+    while ((child = GET_CHILD(parent)) <= heap_size) {
+	/* select child with lower ele, if equal, older child
+	 * older child is older startpoint for flow path, important
+	 * chances are good that the left child will be the older child,
+	 * but just to make sure we test */
+	ele =
+	    alt[SEG_INDEX
+		(alt_seg, astar_pts[heap_index[child]].r,
+		 astar_pts[heap_index[child]].c)];
+	if (child < heap_size) {
+	    childr = child + 1;
+	    i = 1;
+	    while (childr <= heap_size && i < 3) {
+		eler =
+		    alt[SEG_INDEX
+			(alt_seg, astar_pts[heap_index[childr]].r,
+			 astar_pts[heap_index[childr]].c)];
+		if (eler < ele) {
+		    child = childr;
+		    ele = eler;
+		    /* make sure we get the oldest child */
+		}
+		else if (ele == eler &&
+			 heap_index[child] > heap_index[childr]) {
+		    child = childr;
+		}
+		childr++;
+		i++;
+	    }
 	}
-	point.nxt = astar_pts[p].nxt;
-	if (point.nxt == -1) {
-	    astar_pts[p].nxt = nxt_avail_pt;
-	    astar_pts[nxt_avail_pt].r = r;
-	    astar_pts[nxt_avail_pt].c = c;
-	    astar_pts[nxt_avail_pt].downr = downr;
-	    astar_pts[nxt_avail_pt].downc = downc;
-	    astar_pts[nxt_avail_pt].nxt = -1;
-	    nxt_avail_pt++;
 
-	    return 0;
+	/* move hole down */
+
+	heap_index[parent] = heap_index[child];
+	parent = child;
+
+    }
+
+    /* hole is in lowest layer, move to heap end */
+    if (parent < heap_size) {
+	heap_index[parent] = heap_index[heap_size];
+
+	ele =
+	    alt[SEG_INDEX
+		(alt_seg, astar_pts[heap_index[parent]].r,
+		 astar_pts[heap_index[parent]].c)];
+	/* sift up last swapped point, only necessary if hole moved to heap end */
+	sift_up(parent, ele);
+
+    }
+
+    /* the actual drop */
+    heap_size--;
+
+    return 0;
+}
+
+/* standard sift-up routine for d-ary min heap */
+int sift_up(int start, CELL ele)
+{
+    int parent, parentp, child, childp;
+    CELL elep;
+
+    child = start;
+    childp = heap_index[child];
+
+    while (child > 1) {
+	parent = GET_PARENT(child);
+	parentp = heap_index[parent];
+
+	elep =
+	    alt[SEG_INDEX
+		(alt_seg, astar_pts[parentp].r, astar_pts[parentp].c)];
+	/* parent ele higher */
+	if (elep > ele) {
+
+	    /* push parent point down */
+	    heap_index[child] = parentp;
+	    child = parent;
+
+	}
+	/* same ele, but parent is younger */
+	else if (elep == ele && parentp > childp) {
+
+	    /* push parent point down */
+	    heap_index[child] = parentp;
+	    child = parent;
+
 	}
-	p = astar_pts[p].nxt;
+	else
+	    /* no more sifting up, found new slot for child */
+	    break;
+    }
+
+    /* set heap_index for child */
+    if (child < start) {
+	heap_index[child] = childp;
     }
 
     return 0;
+
 }
 
 double
@@ -135,25 +265,34 @@ get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
 	slope = (ele - downe) / window.ns_res;
     else
 	slope = (ele - downe) / diag;
+
     if (slope < MIN_SLOPE)
 	return (MIN_SLOPE);
+
     return (slope);
 }
 
+
+/* new replace */
 int replace(			/* ele was in there */
 	       SHORT upr, SHORT upc, SHORT r, SHORT c)
 /* CELL ele;  */
 {
-    int now;
+    int now, heap_run;
+
+    /* find the current neighbour point and 
+     * set flow direction to focus point */
+
+    heap_run = 0;
 
-    now = first_astar;
-    while (now != -1) {
+    while (heap_run <= heap_size) {
+	now = heap_index[heap_run];
 	if (astar_pts[now].r == upr && astar_pts[now].c == upc) {
 	    astar_pts[now].downr = r;
 	    astar_pts[now].downc = c;
 	    return 0;
 	}
-	now = astar_pts[now].nxt;
+	heap_run++;
     }
     return 0;
 }

raster/r.watershed/ram/do_astar.h

@@ -0,0 +1,7 @@
+#ifndef __DO_ASTAR_H__
+#define __DO_ASTAR__
+
+#define GET_PARENT(c) (int) (((c) - 2) / 3 + 1)
+#define GET_CHILD(p) (int) (((p) * 3) - 1)
+
+#endif /* __DO_ASTAR_H__ */

raster/r.watershed/ram/find_pour.c

@@ -8,6 +8,7 @@ int find_pourpts(void)
 
     basin_num = 0;
     for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 3);
 	northing = window.north - (row + .5) * window.ns_res;
 	for (col = 0; col < ncols; col++) {
 	    value = FLAG_GET(swale, row, col);
@@ -33,6 +34,7 @@ int find_pourpts(void)
 	    }
 	}
     }
+    G_percent(nrows, nrows, 1);	/* finish it */
 
     return 0;
 }

raster/r.watershed/ram/flag_create.c

@@ -15,7 +15,7 @@ FLAG *flag_create(int nrows, int ncols)
 	(unsigned char **)G_malloc(nrows * sizeof(unsigned char *));
 
     temp =
-	(unsigned char *)G_calloc(nrows * new_flag->leng,
+	(unsigned char *)G_malloc(nrows * new_flag->leng *
 				  sizeof(unsigned char));
 
     for (i = 0; i < nrows; i++) {

raster/r.watershed/ram/init_vars.c

@@ -148,7 +148,7 @@ int init_vars(int argc, char *argv[])
 		wat[SEG_INDEX(wat_seg, r, c)] = 1;
 	}
     }
-    asp = (CELL *) G_calloc(size_array(&asp_seg, nrows, ncols), sizeof(CELL));
+    asp = (CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
 
     if (pit_flag) {
 	fd = G_open_cell_old(pit_name, "");
@@ -224,8 +224,15 @@ int init_vars(int argc, char *argv[])
 			       sizeof(double));
     }
 
+    /* heap_index will track astar_pts in the binary min-heap */
+    /* heap_index is one-based */
+    heap_index = (int *)G_malloc((do_points + 1) * sizeof(int));
+
     G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
 
+    /* heap is empty */
+    heap_size = 0;
+
     first_astar = first_cum = -1;
     if (MASK_flag) {
 	for (r = 0; r < nrows; r++) {

raster/r.watershed/ram/main.c

@@ -5,9 +5,10 @@
  * AUTHOR(S):    Charles Ehlschlaeger, CERL (original contributor)
  *               Markus Neteler <neteler itc.it>, Roberto Flor <flor itc.it>, 
  *               Brad Douglas <rez touchofmadness.com>,
- *		 Hamish Bowman <hamish_b yahoo com>
+ *		 Hamish Bowman <hamish_b yahoo com>,
+ *               Markus Metz <markus.metz.giswork gmail.com>
  * PURPOSE:      Watershed determination
- * COPYRIGHT:    (C) 1999-2006 by the GRASS Development Team
+ * COPYRIGHT:    (C) 1999-2008 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
@@ -23,6 +24,7 @@
 
 struct Cell_head window;
 
+int *heap_index, heap_size;
 int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 SHORT nrows, ncols;
 double half_res, diag, max_length, dep_slope;

raster/r.watershed/seg/Gwater.h

@@ -18,10 +18,10 @@
 #define MIN_GRADIENT_DEGREES	1
 #define DEG_TO_RAD		((2 * M_PI) / 360.)
 #define METER_TO_FOOT		(1 / 0.3048)
-#define MAX_BYTES		2000000
-#define PAGE_BLOCK		512
-#define SROW			9
-#define SCOL   			13
+#define MAX_BYTES		10485760
+#define PAGE_BLOCK		1024
+#define SROW			200
+#define SCOL   			200
 #define RITE			1
 #define LEFT			2
 #define NEITHER			0
@@ -35,8 +35,16 @@ POINT {
     int nxt;
 };
 
+#define HEAP    struct heap_item
+HEAP {
+   int point;
+   CELL ele;
+};
+
 extern struct Cell_head window;
 
+extern SSEG heap_index;
+extern int heap_size;
 extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 extern SHORT nrows, ncols;
 extern double half_res, diag, max_length, dep_slope;
@@ -78,6 +86,7 @@ CELL def_basin(int, int, CELL, double, CELL);
 /* do_astar.c */
 int do_astar(void);
 int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int drop_pt(void);
 double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
 int replace(SHORT, SHORT, SHORT, SHORT);
 

raster/r.watershed/seg/do_astar.c

@@ -1,47 +1,74 @@
 #include <stdlib.h>
 #include <unistd.h>
-#include "Gwater.h"
 #include <grass/gis.h>
 #include <grass/glocale.h>
+#include "Gwater.h"
+#include "do_astar.h"
 
+int sift_up(int, CELL);
 
 int do_astar(void)
 {
     POINT point;
     int doer, count;
-    double get_slope();
     SHORT upr, upc, r, c, ct_dir;
     CELL work_val, alt_val, alt_up, asp_up, wat_val;
     CELL in_val, drain_val;
-    double slope;
+    HEAP heap_pos;
+
+    /* double slope; */
 
     G_message(_("SECTION 2: A * Search."));
 
     count = 0;
+    seg_get(&heap_index, (char *)&heap_pos, 0, 1);
+    first_astar = heap_pos.point;
+
+    /* A * Search: downhill path through all not masked cells */
     while (first_astar != -1) {
 	G_percent(count++, do_points, 1);
-	doer = first_astar;
+
+	seg_get(&heap_index, (char *)&heap_pos, 0, 1);
+	doer = heap_pos.point;
+
 	seg_get(&astar_pts, (char *)&point, 0, doer);
-	first_astar = point.nxt;
+
+	/* drop astar_pts[doer] from heap */
+	drop_pt();
+
+	seg_get(&heap_index, (char *)&heap_pos, 0, 1);
+	first_astar = heap_pos.point;
+
 	point.nxt = first_cum;
 	seg_put(&astar_pts, (char *)&point, 0, doer);
+
 	first_cum = doer;
 	r = point.r;
 	c = point.c;
+
 	bseg_put(&worked, &one, r, c);
 	cseg_get(&alt, &alt_val, r, c);
+
+	/* check all neighbours, breadth first search */
 	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+	    /* get r, c (upr, upc) for this neighbour */
 	    upr = r + nextdr[ct_dir];
 	    upc = c + nextdc[ct_dir];
+	    /* check that upr, upc are within region */
 	    if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
+		/* check if neighbour is in the list */
+		/* if not, add as new point */
 		bseg_get(&in_list, &in_val, upr, upc);
 		if (in_val == 0) {
 		    cseg_get(&alt, &alt_up, upr, upc);
 		    add_pt(upr, upc, r, c, alt_up, alt_val);
+		    /* flow direction is set here */
 		    drain_val = drain[upr - r + 1][upc - c + 1];
 		    cseg_put(&asp, &drain_val, upr, upc);
 		}
 		else {
+		    /* check if neighbour has not been worked on,
+		     * update values for asp, wat and slp */
 		    bseg_get(&worked, &work_val, upr, upc);
 		    if (!work_val) {
 			cseg_get(&asp, &asp_up, upr, upc);
@@ -54,9 +81,8 @@ int do_astar(void)
 			    cseg_put(&wat, &wat_val, r, c);
 			    cseg_get(&alt, &alt_up, upr, upc);
 			    replace(upr, upc, r, c);	/* alt_up used to be */
-			    slope =
-				get_slope(upr, upc, r, c, alt_up, alt_val);
-			    dseg_put(&slp, &slope, upr, upc);
+			    /* slope = get_slope (upr, upc, r, c, alt_up, alt_val);
+			       dseg_put (&slp, &slope, upr, upc); */
 			}
 		    }
 		}
@@ -65,58 +91,176 @@ int do_astar(void)
     }
     bseg_close(&worked);
     bseg_close(&in_list);
+    seg_close(&heap_index);
 
     G_percent(count, do_points, 1);	/* finish it */
     return 0;
 }
 
+/* new add point routine for min heap */
 int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
 {
-    int p;
-    CELL check_ele;
-    POINT point, new_point;
-    double slp_value, get_slope();
+    POINT point;
+    HEAP heap_pos;
 
-    if (nxt_avail_pt == do_points)
-	G_fatal_error(_("problem w/ astar algorithm"));
+    /* double slp_value; */
 
-    slp_value = get_slope(r, c, downr, downc, ele, downe);
-    dseg_put(&slp, &slp_value, r, c);
+    /* slp_value = get_slope(r, c, downr, downc, ele, downe);
+       dseg_put (&slp, &slp_value, r, c); */
     bseg_put(&in_list, &one, r, c);
-    new_point.r = r;
-    new_point.c = c;
-    new_point.downr = downr;
-    new_point.downc = downc;
-    if (first_astar == -1) {
-	new_point.nxt = -1;
-	seg_put(&astar_pts, (char *)&new_point, 0, nxt_avail_pt);
-	first_astar = nxt_avail_pt;
-	nxt_avail_pt++;
+
+    /* add point to next free position */
+
+    heap_size++;
+
+    if (heap_size > do_points)
+	G_fatal_error(_("heapsize too large"));
+
+    heap_pos.point = nxt_avail_pt;
+    heap_pos.ele = ele;
+    seg_put(&heap_index, (char *)&heap_pos, 0, heap_size);
+
+    point.r = r;
+    point.c = c;
+    point.downr = downr;
+    point.downc = downc;
+
+    seg_put(&astar_pts, (char *)&point, 0, nxt_avail_pt);
+
+    nxt_avail_pt++;
+
+    /* sift up: move new point towards top of heap */
+
+    sift_up(heap_size, ele);
+
+    return 0;
+}
+
+/* drop point routine for min heap */
+int drop_pt(void)
+{
+    int child, childr, parent;
+    int childp, childrp, parentp;
+    CELL ele, eler;
+    int i;
+    POINT point;
+    HEAP heap_pos;
+
+    if (heap_size == 1) {
+	parent = -1;
+	heap_pos.point = -1;
+	heap_pos.ele = 0;
+	seg_put(&heap_index, (char *)&heap_pos, 0, 1);
+	heap_size = 0;
 	return 0;
     }
-    p = first_astar;;
-    while (1) {
-	seg_get(&astar_pts, (char *)&point, 0, p);
-	cseg_get(&alt, &check_ele, point.r, point.c);
-	if (check_ele > ele) {
-	    new_point.nxt = nxt_avail_pt;
-	    seg_put(&astar_pts, (char *)&new_point, 0, p);
-	    seg_put(&astar_pts, (char *)&point, 0, nxt_avail_pt);
-	    nxt_avail_pt++;
-	    return 0;
+
+    /* start with heap root */
+    parent = 1;
+
+    /* sift down: move hole back towards bottom of heap */
+    /* sift-down routine customised for A * Search logic */
+
+    while ((child = GET_CHILD(parent)) <= heap_size) {
+	/* select child with lower ele, if both are equal, older child
+	 * older child is older startpoint for flow path, important */
+	seg_get(&heap_index, (char *)&heap_pos, 0, child);
+	childp = heap_pos.point;
+	ele = heap_pos.ele;
+	if (child < heap_size) {
+	    childr = child + 1;
+	    i = 1;
+	    while (childr <= heap_size && i < 3) {
+		seg_get(&heap_index, (char *)&heap_pos, 0, childr);
+		childrp = heap_pos.point;
+		eler = heap_pos.ele;
+		if (eler < ele) {
+		    child = childr;
+		    childp = childrp;
+		    ele = eler;
+		}
+		/* make sure we get the oldest child */
+		else if (ele == eler && childp > childrp) {
+		    child = childr;
+		    childp = childrp;
+		}
+		childr++;
+		i++;
+	    }
 	}
-	if (point.nxt == -1) {
-	    point.nxt = nxt_avail_pt;
-	    seg_put(&astar_pts, (char *)&point, 0, p);
-	    new_point.nxt = -1;
-	    seg_put(&astar_pts, (char *)&new_point, 0, nxt_avail_pt);
-	    nxt_avail_pt++;
 
-	    return 0;
+	/* move hole down */
+	heap_pos.point = childp;
+	heap_pos.ele = ele;
+	seg_put(&heap_index, (char *)&heap_pos, 0, parent);
+	parent = child;
+
+    }
+
+    /* hole is in lowest layer, move to heap end */
+    if (parent < heap_size) {
+	seg_get(&heap_index, (char *)&heap_pos, 0, heap_size);
+	seg_put(&heap_index, (char *)&heap_pos, 0, parent);
+
+	ele = heap_pos.ele;
+
+	/* sift up last swapped point, only necessary if hole moved to heap end */
+	sift_up(parent, ele);
+
+    }
+
+    /* the actual drop */
+    heap_size--;
+
+    return 0;
+
+}
+
+/* standard sift-up routine for d-ary min heap */
+int sift_up(int start, CELL ele)
+{
+    int parent, parentp, child, childp;
+    POINT point;
+    CELL elep;
+    HEAP heap_pos;
+
+    child = start;
+    seg_get(&heap_index, (char *)&heap_pos, 0, child);
+    childp = heap_pos.point;
+
+    while (child > 1) {
+	parent = GET_PARENT(child);
+	seg_get(&heap_index, (char *)&heap_pos, 0, parent);
+	parentp = heap_pos.point;
+	elep = heap_pos.ele;
+
+	/* parent ele higher */
+	if (elep > ele) {
+
+	    /* push parent point down */
+	    seg_put(&heap_index, (char *)&heap_pos, 0, child);
+	    child = parent;
+
 	}
-	p = point.nxt;
+	/* same ele, but parent is younger */
+	else if (elep == ele && parentp > childp) {
+
+	    /* push parent point down */
+	    seg_put(&heap_index, (char *)&heap_pos, 0, child);
+	    child = parent;
+
+	}
+	else
+	    /* no more sifting up, found new slot for child */
+	    break;
     }
 
+    /* no more sifting up, found new slot for child */
+    if (child < start) {
+	heap_pos.point = childp;
+	heap_pos.ele = ele;
+	seg_put(&heap_index, (char *)&heap_pos, 0, child);
+    }
     return 0;
 }
 
@@ -140,11 +284,15 @@ int replace(			/* ele was in there */
 	       SHORT upr, SHORT upc, SHORT r, SHORT c)
 /* CELL ele;  */
 {
-    int now;
+    int now, heap_run;
     POINT point;
+    HEAP heap_pos;
+
+    heap_run = 0;
 
-    now = first_astar;
-    while (now != -1) {
+    while (heap_run <= heap_size) {
+	seg_get(&heap_index, (char *)&heap_pos, 0, heap_run);
+	now = heap_pos.point;
 	seg_get(&astar_pts, (char *)&point, 0, now);
 	if (point.r == upr && point.c == upc) {
 	    point.downr = r;
@@ -152,7 +300,7 @@ int replace(			/* ele was in there */
 	    seg_put(&astar_pts, (char *)&point, 0, now);
 	    return 0;
 	}
-	now = point.nxt;
+	heap_run++;;
     }
     return 0;
 }

raster/r.watershed/seg/do_astar.h

@@ -0,0 +1,7 @@
+#ifndef __DO_ASTAR_H__
+#define __DO_ASTAR__
+
+#define GET_PARENT(c) ((int) (((c) - 2) / 3 + 1))
+#define GET_CHILD(p) ((int) ((p) * 3 - 1))
+
+#endif /* __DO_ASTAR_H__ */

raster/r.watershed/seg/find_pour.c

@@ -8,6 +8,7 @@ int find_pourpts(void)
 
     basin_num = 0;
     for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 3);
 	northing = window.north - (row + .5) * window.ns_res;
 	for (col = 0; col < ncols; col++) {
 	    /* cseg_get (&wat, &value, row, col);
@@ -39,6 +40,7 @@ int find_pourpts(void)
 	    }
 	}
     }
+    G_percent(nrows, nrows, 1);	/* finish it */
 
     return 0;
 }

raster/r.watershed/seg/init_vars.c

@@ -8,7 +8,13 @@
 int init_vars(int argc, char *argv[])
 {
     SHORT r, c;
-    int fd, num_cseg_total, num_cseg, num_cseg_bytes;
+    int fd, num_cseg_total, num_open_segs;
+    int seg_rows, seg_cols;
+    double segs_mb;
+    char *mb_opt;
+
+    /* int page_block, num_cseg; */
+    int max_bytes;
     CELL *buf, alt_value, wat_value, asp_value, worked_value;
     char MASK_flag;
 
@@ -21,7 +27,7 @@ int init_vars(int argc, char *argv[])
     max_length = dzero = 0.0;
     ril_value = -1.0;
     /* dep_slope = 0.0; */
-    num_cseg_bytes = MAX_BYTES - 4 * PAGE_BLOCK;
+    max_bytes = 0;
     sides = 8;
     for (r = 1; r < argc; r++) {
 	if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
@@ -54,6 +60,11 @@ int init_vars(int argc, char *argv[])
 	    ls_flag++;
 	else if (sscanf(argv[r], "ob=%[^\n]", ob_name) == 1)
 	    ob_flag++;
+	else if (sscanf(argv[r], "mb=%[^\n]", mb_opt) == 1) {
+	    if (sscanf(mb_opt, "%lf", &segs_mb) == 0) {
+		segs_mb = 300;
+	    }
+	}
 	else if (sscanf(argv[r], "r=%[^\n]", ril_name) == 1) {
 	    if (sscanf(ril_name, "%lf", &ril_value) == 0) {
 		ril_value = -1.0;
@@ -110,32 +121,46 @@ int init_vars(int argc, char *argv[])
 		window.ns_res * window.ns_res);
     if (sides == 4)
 	diag *= 0.5;
-    if (ls_flag)
-	num_cseg_bytes -= PAGE_BLOCK;
-    if (sg_flag)
-	num_cseg_bytes -= PAGE_BLOCK;
-    if (ril_flag)
-	num_cseg_bytes -= PAGE_BLOCK;
-    /* if (dep_flag == -1) num_cseg_bytes -= PAGE_BLOCK; */
-    if (sl_flag)
-	num_cseg_bytes -= sizeof(double) * SROW * SCOL * 4;
-    num_cseg = sizeof(CELL) * 3 + sizeof(double);
-    num_cseg_bytes /= num_cseg * 4 * SROW * SCOL;
+
+    /* segment parameters: one size fits all. Fine tune? */
+    /* Segment rows and cols: 200 */
+    /* 1 segment open for all rasters: 2.86 MB */
+    /* num_open_segs = segs_mb / 2.86 */
+
+    seg_rows = SROW;
+    seg_cols = SCOL;
+
+    if (segs_mb < 3.0) {
+	segs_mb = 300;
+	G_warning(_("Maximum memory to be used was smaller than 3 MB, set to default = 300 MB."));
+    }
+
+    num_open_segs = segs_mb / 2.86;
+
+    G_debug(1, "segs MB: %.0f", segs_mb);
+    G_debug(1, "seg cols: %d", seg_cols);
+    G_debug(1, "seg rows: %d", seg_rows);
+
     num_cseg_total = nrows / SROW + 1;
-    G_debug(1, "    segments in row:\t%d", num_cseg_total);
+    G_debug(1, "   row segments:\t%d", num_cseg_total);
 
     num_cseg_total = ncols / SCOL + 1;
-    G_debug(1, "segments in columns:\t%d", num_cseg_total);
+    G_debug(1, "column segments:\t%d", num_cseg_total);
+
+    num_cseg_total = (ncols / seg_cols + 1) * (nrows / seg_rows + 1);
+    G_debug(1, " total segments:\t%d", num_cseg_total);
+    G_debug(1, "  open segments:\t%d", num_open_segs);
 
-    num_cseg_total = (ncols / SCOL + 1) * (nrows / SROW + 1);
-    G_debug(1, "     total segments:\t%d", num_cseg_total);
-    G_debug(1, "      open segments:\t%d", num_cseg_bytes);
+    /* nonsense to have more segments open than exist */
+    if (num_open_segs > nrows)
+	num_open_segs = nrows;
+    G_debug(1, "  open segments after adjusting:\t%d", num_open_segs);
 
-    cseg_open(&alt, SROW, SCOL, num_cseg_bytes);
-    cseg_open(&r_h, SROW, SCOL, 4);
+    cseg_open(&alt, seg_rows, seg_cols, num_open_segs);
+    cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
     cseg_read_cell(&alt, ele_name, "");
     cseg_read_cell(&r_h, ele_name, "");
-    cseg_open(&wat, SROW, SCOL, num_cseg_bytes);
+    cseg_open(&wat, seg_rows, seg_cols, num_open_segs);
 
     if (run_flag) {
 	cseg_read_cell(&wat, run_name, "");
@@ -147,7 +172,7 @@ int init_vars(int argc, char *argv[])
 		    exit(EXIT_FAILURE);
 	}
     }
-    cseg_open(&asp, SROW, SCOL, num_cseg_bytes);
+    cseg_open(&asp, seg_rows, seg_cols, num_open_segs);
     if (pit_flag) {
 	cseg_read_cell(&asp, pit_name, "");
     }
@@ -158,7 +183,7 @@ int init_vars(int argc, char *argv[])
 		    exit(EXIT_FAILURE);
 	}
     }
-    bseg_open(&swale, SROW, SCOL, num_cseg_bytes);
+    bseg_open(&swale, seg_rows, seg_cols, num_open_segs);
     if (ob_flag) {
 	bseg_read_cell(&swale, ob_name, "");
     }
@@ -169,11 +194,11 @@ int init_vars(int argc, char *argv[])
 	}
     }
     if (ril_flag) {
-	dseg_open(&ril, 1, (int)PAGE_BLOCK / sizeof(double), 1);
+	dseg_open(&ril, 1, seg_rows * seg_cols, num_open_segs);
 	dseg_read_cell(&ril, ril_name, "");
     }
-    bseg_open(&in_list, SROW, SCOL, num_cseg_bytes);
-    bseg_open(&worked, SROW, SCOL, num_cseg_bytes);
+    bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
+    bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
     MASK_flag = 0;
     do_points = nrows * ncols;
     if (NULL != G_find_file("cell", "MASK", G_mapset())) {
@@ -197,17 +222,27 @@ int init_vars(int argc, char *argv[])
 	    G_free(buf);
 	}
     }
-    dseg_open(&slp, SROW, SCOL, num_cseg_bytes);
-    dseg_open(&s_l, SROW, SCOL, 4);
+    /* dseg_open(&slp, SROW, SCOL, num_open_segs); */
+    dseg_open(&s_l, seg_rows, seg_cols, num_open_segs);
     if (sg_flag)
-	dseg_open(&s_g, 1, (int)PAGE_BLOCK / sizeof(double), 1);
+	dseg_open(&s_g, 1, seg_rows * seg_cols, num_open_segs);
     if (ls_flag)
-	dseg_open(&l_s, 1, (int)PAGE_BLOCK / sizeof(double), 1);
-    seg_open(&astar_pts, 1, do_points, 1, PAGE_BLOCK / sizeof(POINT),
-	     4, sizeof(POINT));
-    first_astar = first_cum = -1;
+	dseg_open(&l_s, 1, seg_rows * seg_cols, num_open_segs);
+    seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols,
+	     num_open_segs, sizeof(POINT));
+
+    /* heap_index will track astar_pts in the binary min-heap */
+    /* heap_index is one-based */
+    seg_open(&heap_index, 1, do_points + 1, 1, seg_rows * seg_cols,
+	     num_open_segs, sizeof(HEAP));
+
     G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
 
+    /* heap is empty */
+    heap_size = 0;
+
+    first_astar = first_cum = -1;
+
     if (MASK_flag) {
 	for (r = 0; r < nrows; r++) {
 	    G_percent(r, nrows, 3);
@@ -343,15 +378,15 @@ int init_vars(int argc, char *argv[])
 		    }
 		    else {
 			bseg_put(&in_list, &zero, r, c);
-			dseg_put(&slp, &dzero, r, c);
+			/* dseg_put(&slp, &dzero, r, c); */
 		    }
 		}
 	    }
 	}
-	G_percent(r, nrows, 3);	/* finish it */
     }
     else {
 	for (r = 0; r < nrows; r++) {
+	    G_percent(r, nrows, 3);
 	    for (c = 0; c < ncols; c++) {
 		bseg_put(&worked, &zero, r, c);
 		dseg_put(&s_l, &half_res, r, c);
@@ -381,11 +416,12 @@ int init_vars(int argc, char *argv[])
 		}
 		else {
 		    bseg_put(&in_list, &zero, r, c);
-		    dseg_put(&slp, &dzero, r, c);
+		    /* dseg_put(&slp, &dzero, r, c); */
 		}
 	    }
 	}
     }
+    G_percent(r, nrows, 3);	/* finish it */
 
     return 0;
 }

raster/r.watershed/seg/main.c

@@ -6,9 +6,10 @@
  *               Markus Neteler <neteler itc.it>, 
  *               Roberto Flor <flor itc.it>,
  *               Brad Douglas <rez touchofmadness.com>, 
- *               Hamish Bowman <hamish_b yahoo.com>
+ *               Hamish Bowman <hamish_b yahoo.com>,
+ *               Markus Metz <markus.metz.giswork gmail.com>
  * PURPOSE:      Watershed determination using the GRASS segmentation lib
- * COPYRIGHT:    (C) 1999-2006 by the GRASS Development Team
+ * COPYRIGHT:    (C) 1999-2008 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
@@ -24,6 +25,8 @@
 
 struct Cell_head window;
 
+SSEG heap_index;
+int heap_size;
 int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 SHORT nrows, ncols;
 double half_res, diag, max_length, dep_slope;