[GRASS-SVN] r39995 - grass/branches/develbranch_6/raster/r.watershed/ram

svn_grass at osgeo.org svn_grass at osgeo.org
Mon Dec 14 05:03:47 EST 2009


Author: mmetz
Date: 2009-12-14 05:03:46 -0500 (Mon, 14 Dec 2009)
New Revision: 39995

Modified:
   grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c
   grass/branches/develbranch_6/raster/r.watershed/ram/no_stream.c
Log:
backporting fixes for flow direction and basin delineation from trunk

Modified: grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c	2009-12-14 09:52:38 UTC (rev 39994)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c	2009-12-14 10:03:46 UTC (rev 39995)
@@ -17,10 +17,11 @@
      * |2| |3|
      * |5|0|6|
      */
-    int nbr_ew[8] = { 0, 1, 2, 3, 1, 0, 0, 1};
-    int nbr_ns[8] = { 0, 1, 2, 3, 3, 2, 3, 2};
+    int nbr_ew[8] = { 0, 1, 2, 3, 1, 0, 0, 1 };
+    int nbr_ns[8] = { 0, 1, 2, 3, 3, 2, 3, 2 };
     double dx, dy, dist_to_nbr[8], ew_res, ns_res;
     double slope[8];
+    int skip_diag;
 
     G_message(_("SECTION 2: A * Search."));
 
@@ -39,6 +40,8 @@
     ew_res = window.ew_res;
     ns_res = window.ns_res;
 
+    G_debug(0, "ew_res %e, ns_res %e", ew_res, ns_res);
+
     count = 0;
     first_astar = heap_index[1];
     first_cum = do_points;
@@ -78,39 +81,46 @@
 		/* if not, add as new point */
 		is_in_list = FLAG_GET(in_list, upr, upc);
 		is_worked = FLAG_GET(worked, upr, upc);
+		skip_diag = 0;
+		/* avoid diagonal flow direction bias */
 		if (!is_worked) {
 		    alt_nbr[ct_dir] = alt[index_up];
-		    slope[ct_dir] = get_slope2(alt_val, alt_nbr[ct_dir], dist_to_nbr[ct_dir]);
-		    /* avoid diagonal flow direction bias */
-		    if (ct_dir > 3) {
-			if (slope[nbr_ew[ct_dir]]) {
+		    slope[ct_dir] =
+			get_slope2(alt_val, alt_nbr[ct_dir],
+				   dist_to_nbr[ct_dir]);
+		}
+		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_slope2(alt_nbr[nbr_ew[ct_dir]], alt_nbr[ct_dir], ew_res))
-				is_in_list = 1;
+			    if (slope[ct_dir] <
+				get_slope2(alt_nbr[nbr_ew[ct_dir]],
+					   alt_nbr[ct_dir], ew_res))
+				skip_diag = 1;
 			}
-			if (!is_in_list && slope[nbr_ns[ct_dir]]) {
+			if (!skip_diag && slope[nbr_ns[ct_dir]] > 0) {
 			    /* slope to ns nbr > slope to center */
-			    if (slope[ct_dir] < get_slope2(alt_nbr[nbr_ns[ct_dir]], alt_nbr[ct_dir], ns_res))
-				is_in_list = 1;
+			    if (slope[ct_dir] <
+				get_slope2(alt_nbr[nbr_ns[ct_dir]],
+					   alt_nbr[ct_dir], ns_res))
+				skip_diag = 1;
 			}
 		    }
 		}
-		
+
 		/* add neighbour as new point if not in the list */
-		if (is_in_list == 0) {
+		if (is_in_list == 0 && skip_diag == 0) {
 		    add_pt(upr, upc, alt_nbr[ct_dir], alt_val);
 		    /* set flow direction */
 		    asp[index_up] = drain[upr - r + 1][upc - c + 1];
 		}
-		else {
-		    if (is_worked == 0) {
-			/* neighbour is edge in list, not yet worked */
-			if (asp[index_up] < 0) {
-			    asp[index_up] = drain[upr - r + 1][upc - c + 1];
+		else if (is_in_list && is_worked == 0) {
+		    /* neighbour is edge in list, not yet worked */
+		    if (asp[index_up] < 0) {
+			asp[index_up] = drain[upr - r + 1][upc - c + 1];
 
-			    if (wat[index_doer] > 0)
-				wat[index_doer] = -wat[index_doer];
-			}
+			if (wat[index_doer] > 0)
+			    wat[index_doer] = -wat[index_doer];
 		    }
 		}
 	    }  /* end if in region */
@@ -142,7 +152,6 @@
 /* new add point routine for min heap */
 int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
 {
-
     FLAG_SET(in_list, r, c);
 
     /* add point to next free position */
@@ -291,7 +300,7 @@
     if (ele == up_ele)
 	return 0.5 / dist;
     else
-	return (double) (up_ele - ele) / dist;
+	return (double)(up_ele - ele) / dist;
 }
 
 /* new replace */
@@ -312,8 +321,8 @@
 	/* 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; */
+	    /* astar_pts[now].downr = r;
+	       astar_pts[now].downc = c; */
 	    return 0;
 	}
 	heap_run++;

Modified: grass/branches/develbranch_6/raster/r.watershed/ram/no_stream.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/no_stream.c	2009-12-14 09:52:38 UTC (rev 39994)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/no_stream.c	2009-12-14 10:03:46 UTC (rev 39995)
@@ -10,6 +10,7 @@
     SHORT updir, riteflag, leftflag, thisdir;
 
     while (1) {
+	bas[SEG_INDEX(bas_seg, row, col)] = basin_num;
 	max_drain = -1;
 	for (r = row - 1, rr = 0; r <= row + 1; r++, rr++) {
 	    for (c = col - 1, cc = 0; c <= col + 1; c++, cc++) {
@@ -65,23 +66,25 @@
 			aspect = asp[SEG_INDEX(asp_seg, r, c)];
 			if (aspect == drain[rr][cc]) {
 			    thisdir = updrain[rr][cc];
-			    if (haf_basin_side(updir,
-					       (SHORT) downdir,
-					       thisdir) == RITE) {
+			    switch (haf_basin_side(updir,
+			                          (SHORT) downdir,
+						  thisdir)) {
+			    case RITE:
 				overland_cells(r, c, basin_num, basin_num,
 					       &new_ele);
 				riteflag++;
-			    }
-			    else {
+				break;
+			    case LEFT:
 				overland_cells(r, c, basin_num, basin_num - 1,
 					       &new_ele);
 				leftflag++;
+				break;
 			    }
 			}
 		    }
 		}
 	    }
-	    if (leftflag >= riteflag)
+	    if (leftflag > riteflag)
 		haf[SEG_INDEX(haf_seg, row, col)] = basin_num - 1;
 	    else
 		haf[SEG_INDEX(haf_seg, row, col)] = basin_num;
@@ -96,6 +99,7 @@
 		    slope = MIN_SLOPE;
 		fprintf(fp, " %f %f\n", slope, stream_length);
 	    }
+	    haf[SEG_INDEX(haf_seg, row, col)] = basin_num;
 	    return 0;
 	}
     }



More information about the grass-commit mailing list