[GRASS-SVN] r39322 -
grass/branches/develbranch_6/raster/r.watershed/ram
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Sep 29 03:31:56 EDT 2009
Author: mmetz
Date: 2009-09-29 03:31:55 -0400 (Tue, 29 Sep 2009)
New Revision: 39322
Modified:
grass/branches/develbranch_6/raster/r.watershed/ram/Gwater.h
grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c
grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c
grass/branches/develbranch_6/raster/r.watershed/ram/init_vars.c
grass/branches/develbranch_6/raster/r.watershed/ram/main.c
grass/branches/develbranch_6/raster/r.watershed/ram/ramseg.c
Log:
A* Search performance improvement
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/Gwater.h 2009-09-29 07:28:49 UTC (rev 39321)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/Gwater.h 2009-09-29 07:31:55 UTC (rev 39322)
@@ -50,7 +50,7 @@
extern RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
extern RAMSEG r_h_seg, dep_seg;
extern RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
-extern POINT *astar_pts;
+extern int *astar_pts;
extern CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
extern DCELL *wat;
extern int ril_fd;
@@ -114,6 +114,7 @@
/* ramseg.c */
int size_array(int *, int, int);
+int seg_index_rc(int, int, int *, int *);
/* sg_factor.c */
int sg_factor(void);
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c 2009-09-29 07:28:49 UTC (rev 39321)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c 2009-09-29 07:31:55 UTC (rev 39322)
@@ -5,8 +5,7 @@
int do_astar(void)
{
- POINT *point;
- int doer, count;
+ int count;
SHORT upr, upc, r, c, ct_dir;
CELL alt_val, alt_up, asp_up, wat_val;
CELL in_val, drain_val;
@@ -17,44 +16,31 @@
count = 0;
first_astar = heap_index[1];
+ first_cum = do_points;
- /* A * Search: search uphill, get downhill path */
- while (first_astar != -1) {
+ /* A* Search: search uphill, get downhill paths */
+ while (heap_size > 0) {
G_percent(count++, do_points, 1);
/* 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]] */
+ index_doer = astar_pts[1];
- doer = heap_index[1];
-
- point = &(astar_pts[doer]);
-
- /* 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: replace first_astar with heap_index[1] in line 22 */
- first_astar = heap_index[1];
+ /* add astar points to sorted list for flow accumulation */
+ astar_pts[first_cum] = index_doer;
+ first_cum--;
- /* downhill path for flow accumulation is set here */
- /* this path determines the order for flow accumulation calculation */
- point->nxt = first_cum;
- first_cum = doer;
+ seg_index_rc(alt_seg, index_doer, &r, &c);
- r = point->r;
- c = point->c;
+ G_debug(3, "A* Search: row %d, column %d, ", r, c);
- G_debug(3, "R:%2d C:%2d, ", r, c);
-
- index_doer = SEG_INDEX(alt_seg, r, c);
alt_val = alt[index_doer];
FLAG_SET(worked, r, c);
- /* check all neighbours */
+ /* check neighbours */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (upr, upc) for this neighbour */
upr = r + nextdr[ct_dir];
@@ -73,9 +59,8 @@
asp[index_up] = drain_val;
}
else {
- /* check if neighbour has been worked on,
- * if not, update values for asp and wat */
in_val = FLAG_GET(worked, upr, upc);
+ /* neighbour is edge in list, not yet worked */
if (in_val == 0) {
asp_up = asp[index_up];
if (asp_up < 0) {
@@ -84,17 +69,14 @@
wat_val = wat[index_doer];
if (wat_val > 0)
wat[index_doer] = -wat_val;
-
- /* replace(upr, upc, r, c); */ /* alt_up used to be */
}
}
}
}
}
}
- /* this was a lot of indentation, all in the sake of speed... */
- /* improve code aesthetics? */
- G_percent(count, do_points, 3); /* finish it */
+ /* this was a lot of indentation, improve code aesthetics? */
+ G_percent(count, do_points, 1); /* finish it */
if (mfd == 0)
flag_destroy(worked);
@@ -104,6 +86,19 @@
return 0;
}
+/* compare two heap points */
+/* return 1 if a < b else 0 */
+int cmp_pnt(CELL elea, CELL eleb, int addeda, int addedb)
+{
+ if (elea < eleb)
+ return 1;
+ else if (elea == eleb) {
+ if (addeda < addedb)
+ return 1;
+ }
+ return 0;
+}
+
/* new add point routine for min heap */
int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
{
@@ -117,15 +112,10 @@
if (heap_size > do_points)
G_fatal_error(_("heapsize too large"));
- heap_index[heap_size] = nxt_avail_pt;
+ 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; */
+ astar_pts[heap_size] = SEG_INDEX(alt_seg, r, c);
- nxt_avail_pt++;
-
/* sift up: move new point towards top of heap */
sift_up(heap_size, ele);
@@ -155,29 +145,18 @@
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 */
- ele =
- alt[SEG_INDEX
- (alt_seg, astar_pts[heap_index[child]].r,
- astar_pts[heap_index[child]].c)];
+ ele = alt[astar_pts[child]];
if (child < heap_size) {
childr = child + 1;
- i = child + 3; /* change the number, GET_CHILD() and GET_PARENT() to play with different d-ary heaps */
+ i = child + 3;
while (childr <= heap_size && childr < i) {
- eler =
- alt[SEG_INDEX
- (alt_seg, astar_pts[heap_index[childr]].r,
- astar_pts[heap_index[childr]].c)];
- if (eler < ele) {
+ eler = alt[astar_pts[childr]];
+ /* get smallest child */
+ if (cmp_pnt(eler, ele, heap_index[childr], heap_index[child])) {
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++; */
}
/* break if childr > last entry? that saves sifting up again
* OTOH, this is another comparison
@@ -192,6 +171,7 @@
/* move hole down */
heap_index[parent] = heap_index[child];
+ astar_pts[parent] = astar_pts[child];
parent = child;
}
@@ -199,14 +179,11 @@
/* hole is in lowest layer, move to heap end */
if (parent < heap_size) {
heap_index[parent] = heap_index[heap_size];
+ astar_pts[parent] = astar_pts[heap_size];
- ele =
- alt[SEG_INDEX
- (alt_seg, astar_pts[heap_index[parent]].r,
- astar_pts[heap_index[parent]].c)];
+ ele = alt[astar_pts[parent]];
/* sift up last swapped point, only necessary if hole moved to heap end */
sift_up(parent, ele);
-
}
/* the actual drop */
@@ -218,43 +195,33 @@
/* standard sift-up routine for d-ary min heap */
int sift_up(int start, CELL ele)
{
- register int parent, parentp, child, childp;
+ register int parent, child, child_idx, child_added;
CELL elep;
child = start;
- childp = heap_index[child];
+ child_added = heap_index[child];
+ child_idx = astar_pts[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) {
-
+ elep = alt[astar_pts[parent]];
+ /* child smaller */
+ if (cmp_pnt(ele, elep, child_added, heap_index[parent])) {
/* push parent point down */
- heap_index[child] = parentp;
+ heap_index[child] = heap_index[parent];
+ astar_pts[child] = astar_pts[parent];
child = parent;
-
}
- /* 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 */
+ /* put point in new slot */
if (child < start) {
- heap_index[child] = childp;
+ heap_index[child] = child_added;
+ astar_pts[child] = child_idx;
}
return 0;
@@ -280,12 +247,13 @@
}
-/* no longer needed */
+/* new replace */
int replace( /* ele was in there */
SHORT upr, SHORT upc, SHORT r, SHORT c)
/* CELL ele; */
{
int now, heap_run;
+ int r2, c2;
/* find the current neighbour point and
* set flow direction to focus point */
@@ -294,7 +262,9 @@
while (heap_run <= heap_size) {
now = heap_index[heap_run];
- if (astar_pts[now].r == upr && astar_pts[now].c == upc) {
+ /* if (astar_pts[now].r == upr && astar_pts[now].c == upc) { */
+ seg_index_rc(alt_seg, astar_pts[now], &r2, &c2);
+ if (r2 == upr && c2 == upc) {
/*astar_pts[now].downr = r;
astar_pts[now].downc = c; */
return 0;
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c 2009-09-29 07:28:49 UTC (rev 39321)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c 2009-09-29 07:31:55 UTC (rev 39322)
@@ -10,6 +10,7 @@
int killer, threshold, count;
SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int this_index, down_index;
G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
@@ -18,13 +19,11 @@
threshold = 60;
else
threshold = bas_thres;
- while (first_cum != -1) {
- G_percent(count++, do_points, 2);
- killer = first_cum;
- first_cum = astar_pts[killer].nxt;
- r = astar_pts[killer].r;
- c = astar_pts[killer].c;
- aspect = asp[SEG_INDEX(asp_seg, r, c)];
+ for (killer = 1; killer <= do_points; killer++) {
+ G_percent(killer, do_points, 1);
+ this_index = astar_pts[killer];
+ aspect = asp[this_index];
+ seg_index_rc(alt_seg, this_index, &r, &c);
if (aspect) {
dr = r + asp_r[ABS(aspect)];
dc = c + asp_c[ABS(aspect)];
@@ -32,10 +31,11 @@
else
dr = dc = -1;
if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
- value = wat[SEG_INDEX(wat_seg, r, c)];
+ down_index = SEG_INDEX(wat_seg, dr, dc);
+ value = wat[this_index];
if ((int)(ABS(value) + 0.5) >= threshold)
FLAG_SET(swale, r, c);
- valued = wat[SEG_INDEX(wat_seg, dr, dc)];
+ valued = wat[down_index];
if (value > 0) {
if (valued > 0)
valued += value;
@@ -48,14 +48,14 @@
else
valued = value - valued;
}
- wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
+ wat[down_index] = valued;
valued = ABS(valued) + 0.5;
is_swale = FLAG_GET(swale, r, c);
/* update asp for depression */
if (is_swale && pit_flag) {
- if (aspect > 0 && asp[SEG_INDEX(asp_seg, dr, dc)] == 0) {
+ if (aspect > 0 && asp[down_index] == 0) {
aspect = -aspect;
- asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ asp[this_index] = aspect;
}
}
if (is_swale || ((int)valued) >= threshold) {
@@ -67,7 +67,6 @@
}
}
}
- G_percent(count, do_points, 1); /* finish it */
G_free(astar_pts);
return 0;
@@ -112,6 +111,7 @@
int workedon, edge;
SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int this_index, down_index, nbr_index;
G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
G_debug(1, "MFD convergence factor set to %d.", c_fac);
@@ -142,13 +142,11 @@
else
threshold = bas_thres;
- while (first_cum != -1) {
- G_percent(count++, do_points, 2);
- killer = first_cum;
- first_cum = astar_pts[killer].nxt;
- r = astar_pts[killer].r;
- c = astar_pts[killer].c;
- aspect = asp[SEG_INDEX(asp_seg, r, c)];
+ for (killer = 1; killer <= do_points; killer++) {
+ G_percent(killer, do_points, 1);
+ this_index = astar_pts[killer];
+ seg_index_rc(alt_seg, this_index, &r, &c);
+ aspect = asp[this_index];
if (aspect) {
dr = r + asp_r[ABS(aspect)];
dc = c + asp_c[ABS(aspect)];
@@ -156,8 +154,9 @@
else
dr = dc = -1;
if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
- value = wat[SEG_INDEX(wat_seg, r, c)];
- valued = wat[SEG_INDEX(wat_seg, dr, dc)];
+ value = wat[this_index];
+ down_index = SEG_INDEX(wat_seg, dr, dc);
+ valued = wat[down_index];
r_max = dr;
c_max = dc;
@@ -170,7 +169,7 @@
stream_cells = 0;
swale_cells = 0;
astar_not_set = 1;
- ele = alt[SEG_INDEX(alt_seg, r, c)];
+ ele = alt[this_index];
is_null = 0;
edge = 0;
/* this loop is needed to get the sum of weights */
@@ -183,17 +182,19 @@
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols) {
+ nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
+
/* check for swale or stream cells */
is_swale = FLAG_GET(swale, r_nbr, c_nbr);
if (is_swale)
swale_cells++;
- valued = wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)];
+ valued = wat[nbr_index];
if ((ABS(valued) + 0.5) >= threshold)
stream_cells++;
is_worked = FLAG_GET(worked, r_nbr, c_nbr);
if (is_worked == 0) {
- ele_nbr = alt[SEG_INDEX(alt_seg, r_nbr, c_nbr)];
+ ele_nbr = alt[nbr_index];
is_null = G_is_c_null_value(&ele_nbr);
edge = is_null;
if (!is_null && ele_nbr <= ele) {
@@ -233,7 +234,7 @@
is_swale = FLAG_GET(swale, r, c);
if (is_swale && aspect > 0) {
aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
- asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ asp[this_index] = aspect;
}
continue;
}
@@ -267,11 +268,13 @@
is_worked = FLAG_GET(worked, r_nbr, c_nbr);
if (is_worked == 0) {
+ nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
+
weight[ct_dir] = weight[ct_dir] / sum_weight;
/* check everything sums up to 1.0 */
prop += weight[ct_dir];
- valued = wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)];
+ valued = wat[nbr_index];
if (value > 0) {
if (valued > 0)
valued += value * weight[ct_dir];
@@ -284,7 +287,7 @@
else
valued = value * weight[ct_dir] - valued;
}
- wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)] = valued;
+ wat[nbr_index] = valued;
/* get main drainage direction */
if (ABS(valued) >= max_acc) {
@@ -306,7 +309,7 @@
}
if (mfd_cells < 2) {
- valued = wat[SEG_INDEX(wat_seg, dr, dc)];
+ valued = wat[down_index];
if (value > 0) {
if (valued > 0)
valued += value;
@@ -319,22 +322,22 @@
else
valued = value - valued;
}
- wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
+ wat[down_index] = valued;
}
/* update asp */
if (dr != r_max || dc != c_max) {
aspect = drain[r - r_max + 1][c - c_max + 1];
- if (asp[SEG_INDEX(asp_seg, r, c)] < 0)
+ if (asp[this_index] < 0)
aspect = -aspect;
- asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ asp[this_index] = aspect;
}
is_swale = FLAG_GET(swale, r, c);
/* update asp for depression */
if (is_swale && pit_flag) {
if (aspect > 0 && asp[SEG_INDEX(asp_seg, r_max, c_max)] == 0) {
aspect = -aspect;
- asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+ asp[this_index] = aspect;
}
}
/* start new stream */
@@ -355,7 +358,6 @@
FLAG_SET(worked, r, c);
}
}
- G_percent(count, do_points, 1); /* finish it */
if (workedon)
G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
workedon, do_points);
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/init_vars.c 2009-09-29 07:28:49 UTC (rev 39321)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/init_vars.c 2009-09-29 07:31:55 UTC (rev 39322)
@@ -304,7 +304,7 @@
sizeof(double));
}
- astar_pts = (POINT *) G_malloc(do_points * sizeof(POINT));
+ astar_pts = (int *) G_malloc((do_points + 1) * sizeof(int));
/* heap_index will track astar_pts in ternary min-heap */
/* heap_index is one-based */
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/main.c 2009-09-29 07:28:49 UTC (rev 39321)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/main.c 2009-09-29 07:31:55 UTC (rev 39322)
@@ -34,7 +34,7 @@
RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
RAMSEG r_h_seg, dep_seg;
RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
-POINT *astar_pts;
+int *astar_pts;
CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
DCELL *wat;
int ril_fd;
Modified: grass/branches/develbranch_6/raster/r.watershed/ram/ramseg.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/ramseg.c 2009-09-29 07:28:49 UTC (rev 39321)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/ramseg.c 2009-09-29 07:31:55 UTC (rev 39322)
@@ -13,3 +13,16 @@
size -= (*ram_seg << RAMSEGBITS) - ncols;
return (size);
}
+
+/* get r, c from seg_index */
+int seg_index_rc(int ramseg, int seg_index, int *r, int *c)
+{
+ int seg_no, seg_remainder;
+
+ seg_no = seg_index >> DOUBLEBITS;
+ seg_remainder = seg_index - (seg_no << DOUBLEBITS);
+ *r = ((seg_no / ramseg) << RAMSEGBITS) + (seg_remainder >> RAMSEGBITS);
+ *c = ((seg_no - ((*r) >> RAMSEGBITS) * ramseg) << RAMSEGBITS) +
+ seg_remainder - (((*r) & SEGLENLESS) << RAMSEGBITS);
+ return seg_no;
+}
More information about the grass-commit
mailing list