speed up A* Search, fix half basin bug

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@39980 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.watershed/ram/Gwater.h

@@ -39,6 +39,11 @@ POINT {
     /* int nxt; */
 };
 
+#define OC_STACK struct overland_cells_stack
+OC_STACK {
+    int row, col;
+};
+
 extern struct Cell_head window;
 
 extern int mfd, c_fac, abs_acc, ele_scale;
@@ -47,6 +52,9 @@ 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;
 extern int bas_thres, tot_parts;
+extern CELL n_basins;
+extern OC_STACK *ocs;
+extern int ocs_alloced;
 extern FLAG *worked, *in_list, *s_b, *swale;
 extern RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
 extern RAMSEG r_h_seg, dep_seg;

raster/r.watershed/ram/close_maps2.c

@@ -24,13 +24,7 @@ int close_array_seg(void)
 	    theseg = haf;
 	    thesegseg = haf_seg;
 	}
-	max = -9;
-	for (r = 0; r < nrows; r++) {
-	    for (c = 0; c < ncols; c++) {
-		if (max < theseg[SEG_INDEX(thesegseg, r, c)])
-		    max = theseg[SEG_INDEX(thesegseg, r, c)];
-	    }
-	}
+	max = n_basins;
 	G_debug(1, "%d basins created", max);
 	Rast_init_colors(&colors);
 	Rast_make_random_colors(&colors, 1, max);

raster/r.watershed/ram/do_astar.c

@@ -84,6 +84,8 @@ int do_astar(void)
 		    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 */

raster/r.watershed/ram/find_pour.c

@@ -6,6 +6,9 @@ int find_pourpts(void)
     double easting, northing, stream_length;
     CELL old_elev, basin_num;
 
+    ocs_alloced = 2 * bas_thres;
+    ocs = (OC_STACK *)G_malloc(ocs_alloced * sizeof(OC_STACK));
+
     basin_num = 0;
     stream_length = old_elev = 0;
     for (row = 0; row < nrows; row++) {
@@ -36,6 +39,8 @@ int find_pourpts(void)
 	}
     }
     G_percent(nrows, nrows, 1);	/* finish it */
+    n_basins = basin_num;
+    G_free(ocs);
 
     return 0;
 }

raster/r.watershed/ram/main.c

@@ -30,6 +30,9 @@ int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 SHORT nrows, ncols;
 double half_res, diag, max_length, dep_slope;
 int bas_thres, tot_parts;
+CELL n_basins;
+OC_STACK *ocs;
+int ocs_alloced;
 FLAG *worked, *in_list, *s_b, *swale;
 RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
 RAMSEG r_h_seg, dep_seg;

raster/r.watershed/ram/no_stream.c

@@ -10,6 +10,7 @@ int no_stream(int row, int col, CELL basin_num,
     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++) {
@@ -33,29 +34,29 @@ int no_stream(int row, int col, CELL basin_num,
 	    downdir = asp[SEG_INDEX(asp_seg, row, col)];
 	    if (downdir < 0)
 		downdir = -downdir;
-	    if (sides == 8) {
-		if (uprow != row && upcol != col)
-		    stream_length += diag;
-		else if (uprow != row)
-		    stream_length += window.ns_res;
-		else
-		    stream_length += window.ew_res;
-	    }
-	    else {		/* sides == 4 */
-
-		asp_value = asp[SEG_INDEX(asp_seg, uprow, upcol)];
-		if (downdir == 2 || downdir == 6) {
-		    if (asp_value == 2 || asp_value == 6)
+	    if (arm_flag) {
+		if (sides == 8) {
+		    if (uprow != row && upcol != col)
+			stream_length += diag;
+		    else if (uprow != row)
 			stream_length += window.ns_res;
 		    else
-			stream_length += diag;
-		}
-		else {		/* downdir == 4,8 */
-
-		    if (asp_value == 4 || asp_value == 8)
 			stream_length += window.ew_res;
-		    else
-			stream_length += diag;
+		}
+		else {		/* sides == 4 */
+		    asp_value = asp[SEG_INDEX(asp_seg, uprow, upcol)];
+		    if (downdir == 2 || downdir == 6) {
+			if (asp_value == 2 || asp_value == 6)
+			    stream_length += window.ns_res;
+			else
+			    stream_length += diag;
+		    }
+		    else {		/* downdir == 4,8 */
+			if (asp_value == 4 || asp_value == 8)
+			    stream_length += window.ew_res;
+			else
+			    stream_length += diag;
+		    }
 		}
 	    }
 	    riteflag = leftflag = 0;
@@ -65,23 +66,24 @@ int no_stream(int row, int col, CELL basin_num,
 			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) {
-				overland_cells(r, c, basin_num, basin_num,
-					       &new_ele);
-				riteflag++;
-			    }
-			    else {
+			    switch (haf_basin_side
+			            (updir, (SHORT) downdir, thisdir)) {
+			    case LEFT:
 				overland_cells(r, c, basin_num, basin_num - 1,
 					       &new_ele);
 				leftflag++;
+				break;
+			    case RITE:
+				overland_cells(r, c, basin_num, basin_num,
+					       &new_ele);
+				riteflag++;
+				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 +98,7 @@ int no_stream(int row, int col, CELL basin_num,
 		    slope = MIN_SLOPE;
 		fprintf(fp, " %f %f\n", slope, stream_length);
 	    }
+	    haf[SEG_INDEX(haf_seg, row, col)] = basin_num;
 	    return 0;
 	}
     }

raster/r.watershed/ram/over_cells.c

@@ -2,7 +2,7 @@
 #define BIGNEG	-9999999
 
 int
-overland_cells(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
+overland_cells_recursive(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
 {
     int r, rr, c, cc;
     CELL new_ele, new_max_ele, value;
@@ -33,3 +33,59 @@ overland_cells(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
 
     return 0;
 }
+
+/* non-recursive version */
+int
+overland_cells(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
+{
+    int r, c, rr, cc, next_r, next_c;
+    int top = 0, idx;
+
+    /* put root on stack */
+    ocs[top].row = row;
+    ocs[top].col = col;
+    idx = SEG_INDEX(bas_seg, row, col);
+    bas[idx] = basin_num;
+    haf[idx] = haf_num;
+
+    top++;
+
+    while (top) {
+	top--;
+	next_r = ocs[top].row;
+	next_c = ocs[top].col;
+
+	for (r = next_r - 1, rr = 0; r <= next_r + 1; r++, rr++) {
+	    for (c = next_c - 1, cc = 0; c <= next_c + 1; c++, cc++) {
+		if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+		    idx = SEG_INDEX(bas_seg, r, c);
+		    if (asp[idx] == drain[rr][cc]) {
+			if (top >= ocs_alloced) {
+			    ocs_alloced += bas_thres;
+			    ocs = (OC_STACK *)G_realloc(ocs, ocs_alloced * sizeof(OC_STACK));
+			}
+			ocs[top].row = r;
+			ocs[top].col = c;
+			bas[idx] = basin_num;
+			haf[idx] = haf_num;
+			
+			top++;
+		    }
+		}
+	    }
+	}
+    }
+
+    /*
+    if (arm_flag) {
+	if (new_max_ele == BIGNEG) {
+	    cseg_get(&alt, hih_ele, row, col);
+	}
+	else {
+	    *hih_ele = new_max_ele;
+	}
+    }
+    */
+
+    return 0;
+}