Ticket #344: r.w2.diff

./Makefile

@@ -1 +1 @@
 MODULE_TOPDIR = ../..
 
 SUBDIRS = \
+        front \
         ram \
-        seg \
-        front
+        seg
 # inter unused:
 #       shed
 

front/description.html

@@ -419 +419 @@
 </em>
 
 
-<h2>AUTHOR</h2>
+<h2>AUTHORS</h2>
 
+Original version:
 Charles Ehlschlaeger, U.S. Army Construction Engineering Research Laboratory
+<br>
+Speedups: Markus Metz <markus.metz.giswork at gmail.com>
+
 <p>
-<i>Last changed: $Date: 2008-05-17 07:11:45 +1200 (Sat, 17 May 2008) $</i>
+<i>Last changed: $Date: 2008-09-28 01:36:15 +0200 (Sun, 28 Sep 2008) $</i>

front/main.c

@@ -39 +39 @@
     struct Option *opt13;
     struct Option *opt14;
     struct Option *opt15;
+    struct Option *opt16;
     struct Flag *flag_flow;
     struct Flag *flag_seg;
     struct GModule *module;
@@ -176 +177 @@
     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 =
@@ -335 +343 @@
         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);
 

ram/Gwater.h

@@ -54 +54 @@
 
 GLOBAL struct Cell_head window;
 
+GLOBAL int *heap_index, heap_size;
 GLOBAL int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 GLOBAL SHORT nrows, ncols;
 GLOBAL double half_res, diag, max_length, dep_slope;
@@ -104 +105 @@
 /* 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);
 

ram/do_astar.c

@@ -1 +1 @@
 #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 +12 @@
     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 */
                         }
                     }
@@ -60 +95 @@
             }
         }
     }
-    G_percent(count, do_points, 3);     /* finish it */
+    G_percent(count, do_points, 1);     /* 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;
+
         }
-        p = astar_pts[p].nxt;
+        /* same ele, but parent is younger */
+        else if (elep == ele && parentp > childp) {
+
+            /* push parent point down */
+            heap_index[child] = parentp;
+            child = parent;
+
+        }
+        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 +265 @@
         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;
 }

ram/find_pour.c

@@ -8 +8 @@
 
     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 +34 @@
             }
         }
     }
+    G_percent(nrows, nrows, 1); /* finish it */
 
     return 0;
 }

ram/flag_create.c

@@ -15 +15 @@
         (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++) {

ram/init_vars.c

@@ -170 +170 @@
                 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) {
         pit_mapset = do_exist(pit_name);
@@ -247 +247 @@
                                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++) {

ram/main.c

@@ -5 +5 @@
  * 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_nospam yahoo.com>
+ *               Markus Metz <markus.metz.giswork gmail.com>
  * PURPOSE:      Watershed determination
  * COPYRIGHT:    (C) 1999-2006 by the GRASS Development Team
  *

seg/Gwater.h

@@ -18 +18 @@
 #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
@@ -49 +49 @@
 #define NEXTDCVAR
 #endif
 
+#define HEAP    struct heap_item
+HEAP {
+   int point;
+   CELL ele;
+};
+
 GLOBAL struct Cell_head window;
 
+GLOBAL SSEG heap_index;
+GLOBAL int heap_size;
 GLOBAL int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 GLOBAL SHORT nrows, ncols;
 GLOBAL double half_res, diag, max_length, dep_slope;
@@ -97 +105 @@
 /* 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);
 

seg/do_astar.c

@@ -1 +1 @@
 #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 +81 @@
                             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 +91 @@
     }
     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 +284 @@
                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 +300 @@
             seg_put(&astar_pts, (char *)&point, 0, now);
             return 0;
         }
-        now = point.nxt;
+        heap_run++;;
     }
     return 0;
 }

seg/find_pour.c

@@ -8 +8 @@
 
     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 +40 @@
             }
         }
     }
+    G_percent(nrows, nrows, 1); /* finish it */
 
     return 0;
 }

seg/init_vars.c

@@ -8 +8 @@
 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;
     extern FILE *fopen();
     char MASK_flag, *do_exist();
@@ -22 +28 @@
     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)
@@ -55 +61 @@
             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;
@@ -127 +138 @@
                 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 / 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);
+    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);
+
+    /* 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, ele_mapset);
     cseg_read_cell(&r_h, ele_name, ele_mapset);
-    cseg_open(&wat, SROW, SCOL, num_cseg_bytes);
+    cseg_open(&wat, seg_rows, seg_cols, num_open_segs);
 
     if (run_flag) {
         run_mapset = do_exist(run_name);
@@ -165 +190 @@
                     exit(EXIT_FAILURE);
         }
     }
-    cseg_open(&asp, SROW, SCOL, num_cseg_bytes);
+    cseg_open(&asp, seg_rows, seg_cols, num_open_segs);
     if (pit_flag) {
         pit_mapset = do_exist(pit_name);
         cseg_read_cell(&asp, pit_name, pit_mapset);
@@ -177 +202 @@
                     exit(EXIT_FAILURE);
         }
     }
-    bseg_open(&swale, SROW, SCOL, num_cseg_bytes);
+    bseg_open(&swale, seg_rows, seg_cols, num_open_segs);
     if (ob_flag) {
         ob_mapset = do_exist(ob_name);
         bseg_read_cell(&swale, ob_name, ob_mapset);
@@ -190 +215 @@
     }
     if (ril_flag) {
         ril_mapset = do_exist(ril_name);
-        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, ril_mapset);
     }
-    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())) {
@@ -218 +243 @@
             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);
@@ -364 +399 @@
                     }
                     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);
@@ -402 +437 @@
                 }
                 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;
 }

seg/main.c

@@ -7 +7 @@
  *               Roberto Flor <flor itc.it>,
  *               Brad Douglas <rez touchofmadness.com>, 
  *               Hamish Bowman <hamish_nospam 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
  *