[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