+
-Last changed: $Date: 2008-05-17 07:11:45 +1200 (Sat, 17 May 2008) $
+Last changed: $Date: 2008-09-28 01:36:15 +0200 (Sun, 28 Sep 2008) $
diff -u -u r.watershed/front/main.c r.watershed2/front/main.c
--- r.watershed/front/main.c 2008-11-28 19:24:52.000000000 +1300
+++ r.watershed2/front/main.c 2008-11-28 19:19:42.000000000 +1300
@@ -39,6 +39,7 @@
struct Option *opt13;
struct Option *opt14;
struct Option *opt15;
+ struct Option *opt16;
struct Flag *flag_flow;
struct Flag *flag_seg;
struct GModule *module;
@@ -176,6 +177,13 @@
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,6 +343,13 @@
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);
diff -u -u r.watershed/ram/Gwater.h r.watershed2/ram/Gwater.h
--- r.watershed/ram/Gwater.h 2008-11-28 19:50:45.000000000 +1300
+++ r.watershed2/ram/Gwater.h 2008-11-28 19:56:10.000000000 +1300
@@ -54,6 +54,7 @@
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,6 +105,7 @@
/* 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);
diff -u -u r.watershed/ram/do_astar.c r.watershed2/ram/do_astar.c
--- r.watershed/ram/do_astar.c 2008-08-05 18:39:03.000000000 +1200
+++ r.watershed2/ram/do_astar.c 2008-11-28 20:08:40.000000000 +1300
@@ -1,7 +1,9 @@
#include "Gwater.h"
+#include "do_astar.h"
#include
#include
+int sift_up(int, CELL);
int do_astar(void)
{
@@ -10,49 +12,82 @@
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,68 +95,163 @@
}
}
}
- 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,25 +265,34 @@
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;
}
Only in r.watershed2/ram: do_astar.h
diff -u -u r.watershed/ram/find_pour.c r.watershed2/ram/find_pour.c
--- r.watershed/ram/find_pour.c 2008-08-05 18:39:03.000000000 +1200
+++ r.watershed2/ram/find_pour.c 2008-11-29 14:37:47.000000000 +1300
@@ -8,6 +8,7 @@
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 @@
}
}
}
+ G_percent(nrows, nrows, 1); /* finish it */
return 0;
}
diff -u -u r.watershed/ram/flag_create.c r.watershed2/ram/flag_create.c
--- r.watershed/ram/flag_create.c 2008-08-05 18:39:03.000000000 +1200
+++ r.watershed2/ram/flag_create.c 2008-11-28 20:08:40.000000000 +1300
@@ -15,7 +15,7 @@
(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++) {
diff -u -u r.watershed/ram/init_vars.c r.watershed2/ram/init_vars.c
--- r.watershed/ram/init_vars.c 2008-08-05 18:39:03.000000000 +1200
+++ r.watershed2/ram/init_vars.c 2008-11-29 14:41:54.000000000 +1300
@@ -170,7 +170,7 @@
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,8 +247,15 @@
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++) {
diff -u -u r.watershed/ram/main.c r.watershed2/ram/main.c
--- r.watershed/ram/main.c 2008-08-05 18:39:03.000000000 +1200
+++ r.watershed2/ram/main.c 2008-11-28 20:08:40.000000000 +1300
@@ -5,6 +5,7 @@
* AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
* Markus Neteler , Roberto Flor ,
* Brad Douglas , Hamish Bowman
+ * Markus Metz
* PURPOSE: Watershed determination
* COPYRIGHT: (C) 1999-2006 by the GRASS Development Team
*
diff -u -u r.watershed/seg/Gwater.h r.watershed2/seg/Gwater.h
--- r.watershed/seg/Gwater.h 2008-11-28 20:03:37.000000000 +1300
+++ r.watershed2/seg/Gwater.h 2008-11-29 14:49:44.000000000 +1300
@@ -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
@@ -49,8 +49,16 @@
#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,6 +105,7 @@
/* 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);
diff -u -u r.watershed/seg/do_astar.c r.watershed2/seg/do_astar.c
--- r.watershed/seg/do_astar.c 2008-08-05 18:39:04.000000000 +1200
+++ r.watershed2/seg/do_astar.c 2008-11-28 20:08:53.000000000 +1300
@@ -1,47 +1,74 @@
#include
#include
-#include "Gwater.h"
#include
#include
+#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 @@
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 @@
}
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 @@
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 @@
seg_put(&astar_pts, (char *)&point, 0, now);
return 0;
}
- now = point.nxt;
+ heap_run++;;
}
return 0;
}
Only in r.watershed2/seg: do_astar.h
diff -u -u r.watershed/seg/find_pour.c r.watershed2/seg/find_pour.c
--- r.watershed/seg/find_pour.c 2008-08-05 18:39:04.000000000 +1200
+++ r.watershed2/seg/find_pour.c 2008-11-29 14:38:02.000000000 +1300
@@ -8,6 +8,7 @@
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 @@
}
}
}
+ G_percent(nrows, nrows, 1); /* finish it */
return 0;
}
diff -u -u r.watershed/seg/init_vars.c r.watershed2/seg/init_vars.c
--- r.watershed/seg/init_vars.c 2008-08-08 16:25:24.000000000 +1200
+++ r.watershed2/seg/init_vars.c 2008-11-29 14:57:03.000000000 +1300
@@ -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;
extern FILE *fopen();
char MASK_flag, *do_exist();
@@ -22,7 +28,7 @@
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,6 +61,11 @@
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,32 +138,46 @@
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,7 +190,7 @@
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,7 +202,7 @@
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,11 +215,11 @@
}
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,17 +243,27 @@
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,15 +399,15 @@
}
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,11 +437,12 @@
}
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;
}
diff -u -u r.watershed/seg/main.c r.watershed2/seg/main.c
--- r.watershed/seg/main.c 2008-11-29 14:59:29.000000000 +1300
+++ r.watershed2/seg/main.c 2008-11-28 20:08:53.000000000 +1300
@@ -7,6 +7,7 @@
* Roberto Flor ,
* Brad Douglas ,
* Hamish Bowman
+ * Markus Metz
* PURPOSE: Watershed determination using the GRASS segmentation lib
* COPYRIGHT: (C) 1999-2006 by the GRASS Development Team
*