[GRASS-SVN] r34646 - in grass/trunk/raster/r.watershed: front ram
seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Dec 1 02:10:43 EST 2008
Author: hamish
Date: 2008-12-01 02:10:39 -0500 (Mon, 01 Dec 2008)
New Revision: 34646
Added:
grass/trunk/raster/r.watershed/ram/do_astar.h
grass/trunk/raster/r.watershed/seg/do_astar.h
Modified:
grass/trunk/raster/r.watershed/front/main.c
grass/trunk/raster/r.watershed/front/r.watershed.html
grass/trunk/raster/r.watershed/ram/Gwater.h
grass/trunk/raster/r.watershed/ram/do_astar.c
grass/trunk/raster/r.watershed/ram/find_pour.c
grass/trunk/raster/r.watershed/ram/flag_create.c
grass/trunk/raster/r.watershed/ram/init_vars.c
grass/trunk/raster/r.watershed/ram/main.c
grass/trunk/raster/r.watershed/seg/Gwater.h
grass/trunk/raster/r.watershed/seg/do_astar.c
grass/trunk/raster/r.watershed/seg/find_pour.c
grass/trunk/raster/r.watershed/seg/init_vars.c
grass/trunk/raster/r.watershed/seg/main.c
Log:
Merge r.watershed2 from Markus Metz into r.watershed.
(merge from devbr6 r34645)
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.
Modified: grass/trunk/raster/r.watershed/front/main.c
===================================================================
--- grass/trunk/raster/r.watershed/front/main.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/front/main.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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 @@
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 @@
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 @@
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);
Modified: grass/trunk/raster/r.watershed/front/r.watershed.html
===================================================================
--- grass/trunk/raster/r.watershed/front/r.watershed.html 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/front/r.watershed.html 2008-12-01 07:10:39 UTC (rev 34646)
@@ -419,8 +419,13 @@
</em>
-<h2>AUTHOR</h2>
+<h2>AUTHORS</h2>
+Original version:
Charles Ehlschlaeger, U.S. Army Construction Engineering Research Laboratory
+<br>
+Faster sorting algorithm:
+Markus Metz <markus.metz.giswork at gmail.com>
+
<p>
<i>Last changed: $Date$</i>
Modified: grass/trunk/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/Gwater.h 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/ram/Gwater.h 2008-12-01 07:10:39 UTC (rev 34646)
@@ -40,6 +40,7 @@
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 @@
/* 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);
Modified: grass/trunk/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/ram/do_astar.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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 @@
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 = astar_pts[doer].r = point->r;
- c = astar_pts[doer].c = point->c;
+ r = point->r;
+ 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 @@
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;
+
}
- 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;
- now = first_astar;
- while (now != -1) {
+ /* find the current neighbour point and
+ * set flow direction to focus point */
+
+ heap_run = 0;
+
+ 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;
}
Copied: grass/trunk/raster/r.watershed/ram/do_astar.h (from rev 34645, grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.h)
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.h (rev 0)
+++ grass/trunk/raster/r.watershed/ram/do_astar.h 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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__ */
Modified: grass/trunk/raster/r.watershed/ram/find_pour.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/find_pour.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/ram/find_pour.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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;
}
Modified: grass/trunk/raster/r.watershed/ram/flag_create.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/flag_create.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/ram/flag_create.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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++) {
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -148,7 +148,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) {
fd = G_open_cell_old(pit_name, "");
@@ -224,8 +224,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++) {
Modified: grass/trunk/raster/r.watershed/ram/main.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/main.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/ram/main.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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;
Modified: grass/trunk/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/Gwater.h 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/seg/Gwater.h 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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 @@
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 @@
/* 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);
Modified: grass/trunk/raster/r.watershed/seg/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_astar.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/seg/do_astar.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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 @@
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;
- now = first_astar;
- while (now != -1) {
+ heap_run = 0;
+
+ 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;
}
Copied: grass/trunk/raster/r.watershed/seg/do_astar.h (from rev 34645, grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.h)
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_astar.h (rev 0)
+++ grass/trunk/raster/r.watershed/seg/do_astar.h 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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__ */
Modified: grass/trunk/raster/r.watershed/seg/find_pour.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/find_pour.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/seg/find_pour.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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;
}
Modified: grass/trunk/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/init_vars.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/seg/init_vars.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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 @@
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 @@
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 @@
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);
- cseg_open(&alt, SROW, SCOL, num_cseg_bytes);
- cseg_open(&r_h, SROW, SCOL, 4);
+ /* 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, 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 @@
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 @@
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 @@
}
}
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 @@
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 @@
}
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 @@
}
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;
}
Modified: grass/trunk/raster/r.watershed/seg/main.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/main.c 2008-12-01 05:59:51 UTC (rev 34645)
+++ grass/trunk/raster/r.watershed/seg/main.c 2008-12-01 07:10:39 UTC (rev 34646)
@@ -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;
More information about the grass-commit
mailing list