[GRASS-SVN] r50104 - grass-addons/grass7/raster/r.stream.extract
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jan 9 07:05:52 EST 2012
Author: mmetz
Date: 2012-01-09 04:05:52 -0800 (Mon, 09 Jan 2012)
New Revision: 50104
Modified:
grass-addons/grass7/raster/r.stream.extract/do_astar.c
grass-addons/grass7/raster/r.stream.extract/init_search.c
grass-addons/grass7/raster/r.stream.extract/main.c
grass-addons/grass7/raster/r.stream.extract/streams.c
Log:
r.stream.extract: fix memory and disk space calculations, code cleanup for A* search
Modified: grass-addons/grass7/raster/r.stream.extract/do_astar.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/do_astar.c 2012-01-09 11:38:48 UTC (rev 50103)
+++ grass-addons/grass7/raster/r.stream.extract/do_astar.c 2012-01-09 12:05:52 UTC (rev 50104)
@@ -40,7 +40,7 @@
G_message(_("A* Search..."));
- G_get_set_window(&window);
+ Rast_get_window(&window);
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (r_nbr, c_nbr) for neighbours */
@@ -81,68 +81,68 @@
c_nbr = c + nextdc[ct_dir];
slope[ct_dir] = ele_nbr[ct_dir] = 0;
skip_me = 0;
+
/* check that neighbour is within region */
- if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 && c_nbr < ncols) {
+ if (r_nbr < 0 || r_nbr >= nrows || c_nbr < 0 || c_nbr >= ncols)
+ continue;
- bseg_get(&bitflags, &flag_value, r_nbr, c_nbr);
- is_in_list = FLAG_GET(flag_value, INLISTFLAG);
- is_worked = FLAG_GET(flag_value, WORKEDFLAG);
- if (!is_worked) {
- seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
- ele_nbr[ct_dir] = wa.ele;
- slope[ct_dir] =
- get_slope(ele_val, ele_nbr[ct_dir],
- dist_to_nbr[ct_dir]);
- }
- /* avoid diagonal flow direction bias */
- if (!is_in_list) {
- if (ct_dir > 3 && slope[ct_dir] > 0) {
- if (slope[nbr_ew[ct_dir]] > 0) {
- /* slope to ew nbr > slope to center */
- if (slope[ct_dir] <
- get_slope(ele_nbr[nbr_ew[ct_dir]],
- ele_nbr[ct_dir], ew_res))
- skip_me = 1;
- }
- if (!skip_me && slope[nbr_ns[ct_dir]] > 0) {
- /* slope to ns nbr > slope to center */
- if (slope[ct_dir] <
- get_slope(ele_nbr[nbr_ns[ct_dir]],
- ele_nbr[ct_dir], ns_res))
- skip_me = 1;
- }
+ bseg_get(&bitflags, &flag_value, r_nbr, c_nbr);
+ is_in_list = FLAG_GET(flag_value, INLISTFLAG);
+ is_worked = FLAG_GET(flag_value, WORKEDFLAG);
+ if (!is_worked) {
+ seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
+ ele_nbr[ct_dir] = wa.ele;
+ slope[ct_dir] = get_slope(ele_val, ele_nbr[ct_dir],
+ dist_to_nbr[ct_dir]);
+ }
+ /* avoid diagonal flow direction bias */
+ if (!is_in_list) {
+ if (ct_dir > 3 && slope[ct_dir] > 0) {
+ if (slope[nbr_ew[ct_dir]] > 0) {
+ /* slope to ew nbr > slope to center */
+ if (slope[ct_dir] <
+ get_slope(ele_nbr[nbr_ew[ct_dir]],
+ ele_nbr[ct_dir], ew_res))
+ skip_me = 1;
}
+ if (!skip_me && slope[nbr_ns[ct_dir]] > 0) {
+ /* slope to ns nbr > slope to center */
+ if (slope[ct_dir] <
+ get_slope(ele_nbr[nbr_ns[ct_dir]],
+ ele_nbr[ct_dir], ns_res))
+ skip_me = 1;
+ }
}
+ }
- if (is_in_list == 0 && skip_me == 0) {
- ele_up = ele_nbr[ct_dir];
- asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
- bseg_put(&asp, &asp_val, r_nbr, c_nbr);
- heap_add(r_nbr, c_nbr, ele_up);
- FLAG_SET(flag_value, INLISTFLAG);
- bseg_put(&bitflags, &flag_value, r_nbr, c_nbr);
+ if (is_in_list == 0 && skip_me == 0) {
+ ele_up = ele_nbr[ct_dir];
+ asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
+ bseg_put(&asp, &asp_val, r_nbr, c_nbr);
+ heap_add(r_nbr, c_nbr, ele_up);
+ FLAG_SET(flag_value, INLISTFLAG);
+ bseg_put(&bitflags, &flag_value, r_nbr, c_nbr);
+ }
+ else if (is_in_list && is_worked == 0) {
+ if (FLAG_GET(flag_value, EDGEFLAG)) {
+ /* neighbour is edge in list, not yet worked */
+ bseg_get(&asp, &asp_val, r_nbr, c_nbr);
+ if (asp_val < 0) {
+ /* adjust flow direction for edge cell */
+ asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
+ bseg_put(&asp, &asp_val, r_nbr, c_nbr);
+ }
}
- else if (is_in_list && is_worked == 0) {
- if (FLAG_GET(flag_value, EDGEFLAG)) {
- /* neighbour is edge in list, not yet worked */
- bseg_get(&asp, &asp_val, r_nbr, c_nbr);
- if (asp_val < 0) {
- /* adjust flow direction for edge cell */
- asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
- bseg_put(&asp, &asp_val, r_nbr, c_nbr);
- }
+ else if (FLAG_GET(flag_value, DEPRFLAG)) {
+ G_debug(3, "real depression");
+ /* neighbour is inside real depression, not yet worked */
+ bseg_get(&asp, &asp_val, r_nbr, c_nbr);
+ if (asp_val == 0 && ele_val <= ele_nbr[ct_dir]) {
+ asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
+ bseg_put(&asp, &asp_val, r_nbr, c_nbr);
+ FLAG_UNSET(flag_value, DEPRFLAG);
+ bseg_put(&bitflags, &flag_value, r_nbr, c_nbr);
}
- else if (FLAG_GET(flag_value, DEPRFLAG)) {
- G_debug(3, "real depression");
- /* neighbour is inside real depression, not yet worked */
- bseg_get(&asp, &asp_val, r_nbr, c_nbr);
- if (asp_val == 0 && ele_val <= ele_nbr[ct_dir]) {
- asp_val = drain[r_nbr - r + 1][c_nbr - c + 1];
- bseg_put(&asp, &asp_val, r_nbr, c_nbr);
- FLAG_UNSET(flag_value, DEPRFLAG);
- bseg_put(&bitflags, &flag_value, r_nbr, c_nbr);
- }
- }
}
}
} /* end neighbours */
@@ -152,7 +152,7 @@
bseg_get(&bitflags, &flag_value, r, c);
FLAG_SET(flag_value, WORKEDFLAG);
bseg_put(&bitflags, &flag_value, r, c);
- } /* end A* search */
+ } /* end A* search */
G_percent(n_points, n_points, 1); /* finish it */
@@ -190,8 +190,8 @@
seg_put(&search_heap, (char *)&heap_p, 0, child);
child = parent;
}
- /* no more sifting up, found slot for child */
else
+ /* no more sifting up, found slot for child */
break;
}
@@ -271,7 +271,10 @@
parent = child;
}
- seg_put(&search_heap, (char *)&last_p, 0, parent);
+ /* fill hole */
+ if (parent < heap_size) {
+ seg_put(&search_heap, (char *)&last_p, 0, parent);
+ }
/* the actual drop */
heap_size--;
Modified: grass-addons/grass7/raster/r.stream.extract/init_search.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/init_search.c 2012-01-09 11:38:48 UTC (rev 50103)
+++ grass-addons/grass7/raster/r.stream.extract/init_search.c 2012-01-09 12:05:52 UTC (rev 50104)
@@ -37,9 +37,9 @@
if (is_null)
continue;
+ asp_value = 0;
if (r == 0 || r == nrows - 1 || c == 0 || c == ncols - 1) {
- asp_value = 0;
if (r == 0 && c == 0)
asp_value = -7;
else if (r == 0 && c == ncols - 1)
@@ -68,7 +68,6 @@
}
/* any neighbour NULL ? */
- asp_value = 0;
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (r_nbr, c_nbr) for neighbours */
r_nbr = r + nextdr[ct_dir];
Modified: grass-addons/grass7/raster/r.stream.extract/main.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/main.c 2012-01-09 11:38:48 UTC (rev 50103)
+++ grass-addons/grass7/raster/r.stream.extract/main.c 2012-01-09 12:05:52 UTC (rev 50104)
@@ -65,8 +65,9 @@
double threshold, d8cut, mont_exp;
int min_stream_length = 0, memory;
int seg_cols, seg_rows;
+ double seg2kb;
int num_open_segs, num_open_array_segs, num_seg_total;
- double memory_divisor, heap_mem;
+ double memory_divisor, heap_mem, disk_space;
const char *mapset;
G_gisinit(argv[0]);
@@ -278,6 +279,7 @@
/* segment structures */
seg_rows = seg_cols = 64;
+ seg2kb = seg_rows * seg_cols / 1024.;
/* elevation + accumulation: 12 byte -> 48 KB / segment
* aspect: 1 byte -> 4 KB / segment
* stream: 4 byte -> 16 KB / segment
@@ -293,48 +295,65 @@
/* balance segment files */
/* elevation + accumulation: * 2 */
- memory_divisor = 48 * 2;
+ memory_divisor = sizeof(WAT_ALT) * 2;
+ disk_space = sizeof(WAT_ALT);
/* aspect: as is */
- memory_divisor += 4;
+ memory_divisor += sizeof(char);
+ disk_space += sizeof(char);
/* stream ids: / 2 */
- memory_divisor += 16 / 2;
+ memory_divisor += sizeof(int) / 2.;
+ disk_space += sizeof(int);
/* flags: * 4 */
- memory_divisor += 4 * 4;
+ memory_divisor += sizeof(char) * 4;
+ disk_space += sizeof(char);
/* astar_points: / 16 */
/* ideally only a few but large segments */
- memory_divisor += 32 / 16;
- /* heap points: / 5 */
- memory_divisor += 64. / 5.;
+ memory_divisor += sizeof(POINT) / 16.;
+ disk_space += sizeof(POINT);
+ /* heap points: / 4 */
+ memory_divisor += sizeof(HEAP_PNT) / 4.;
+ disk_space += sizeof(HEAP_PNT);
/* KB -> MB */
- memory_divisor /= 1024.;
+ memory_divisor *= seg2kb / 1024.;
+ disk_space *= seg2kb / 1024.;
num_open_segs = memory / memory_divisor;
- heap_mem = num_open_segs * 64. / (5. * 1024.);
+ heap_mem = num_open_segs * seg2kb * sizeof(HEAP_PNT) /
+ (4. * 1024.);
num_seg_total = (ncols / seg_cols + 1) * (nrows / seg_rows + 1);
if (num_open_segs > num_seg_total) {
heap_mem += (num_open_segs - num_seg_total) * memory_divisor;
- heap_mem -= (num_open_segs - num_seg_total) * 64. / (5. * 1024.);
+ heap_mem -= (num_open_segs - num_seg_total) * seg2kb *
+ sizeof(HEAP_PNT) / (4. * 1024.);
num_open_segs = num_seg_total;
}
if (num_open_segs < 16) {
num_open_segs = 16;
- heap_mem = num_open_segs * 64. / 5.;
+ heap_mem = num_open_segs * seg2kb * sizeof(HEAP_PNT) /
+ (4. * 1024.);
}
- G_verbose_message(_("%d of %d segments are kept in memory"),
- num_open_segs, num_seg_total);
+ G_verbose_message(_("%.2f of data are kept in memory"),
+ 100. * num_open_segs / num_seg_total);
+ disk_space *= num_seg_total;
+ if (disk_space < 1024.0)
+ G_verbose_message(_("Will need up to %.2f MB of disk space"), disk_space);
+ else
+ G_verbose_message(_("Will need up to %.2f GB (%.0f MB) of disk space"),
+ disk_space / 1024.0, disk_space);
/* open segment files */
G_verbose_message(_("Create temporary files..."));
seg_open(&watalt, nrows, ncols, seg_rows, seg_cols, num_open_segs * 2,
sizeof(WAT_ALT), 1);
if (num_open_segs * 2 > num_seg_total)
- heap_mem += (num_open_segs * 2 - num_seg_total) * 48. * 2. / 1024.;
+ heap_mem += (num_open_segs * 2 - num_seg_total) * seg2kb *
+ sizeof(WAT_ALT) / 1024.;
cseg_open(&stream, seg_rows, seg_cols, num_open_segs / 2.);
bseg_open(&asp, seg_rows, seg_cols, num_open_segs);
bseg_open(&bitflags, seg_rows, seg_cols, num_open_segs * 4);
if (num_open_segs * 4 > num_seg_total)
- heap_mem += (num_open_segs * 4 - num_seg_total) * 4. * 4. / 1024.;
+ heap_mem += (num_open_segs * 4 - num_seg_total) * seg2kb / 1024.;
/* load maps */
if (load_maps(ele_fd, acc_fd) < 0)
@@ -371,7 +390,7 @@
if (n_points % seg_cols > 0)
num_seg_total++;
/* no need to have more segments open than exist */
- num_open_array_segs = 1024. * 1024. * heap_mem / (seg_cols * 16.);
+ num_open_array_segs = (1 << 20) * heap_mem / (seg_cols * sizeof(HEAP_PNT));
if (num_open_array_segs > num_seg_total)
num_open_array_segs = num_seg_total;
if (num_open_array_segs < 2)
Modified: grass-addons/grass7/raster/r.stream.extract/streams.c
===================================================================
--- grass-addons/grass7/raster/r.stream.extract/streams.c 2012-01-09 11:38:48 UTC (rev 50103)
+++ grass-addons/grass7/raster/r.stream.extract/streams.c 2012-01-09 12:05:52 UTC (rev 50104)
@@ -6,14 +6,10 @@
double mfd_pow(double base)
{
- int x;
- double result;
+ int i;
+ double result = base;
- result = base;
- if (c_fac == 1)
- return result;
-
- for (x = 2; x <= c_fac; x++) {
+ for (i = 2; i <= c_fac; i++) {
result *= base;
}
return result;
More information about the grass-commit
mailing list