A* Search performance improvement

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

raster/r.watershed/ram/Gwater.h

@@ -36,7 +36,7 @@
 #define POINT       struct points
 POINT {
     SHORT r, c; /* , downr, downc */
-    int nxt;
+    /* int nxt; */
 };
 
 extern struct Cell_head window;
@@ -51,7 +51,7 @@ 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;
 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;
@@ -115,6 +115,7 @@ int overland_cells(int, int, CELL, CELL, CELL *);
 
 /* ramseg.c */
 int size_array(int *, int, int);
+int seg_index_rc(int, int, int *, int *);
 
 /* sg_factor.c */
 int sg_factor(void);

raster/r.watershed/ram/do_astar.c

@@ -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 @@ int do_astar(void)
 
     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];
-
-	/* downhill path for flow accumulation is set here */
-	/* this path determines the order for flow accumulation calculation */
-	point->nxt = first_cum;
-	first_cum = doer;
+	/* add astar points to sorted list for flow accumulation */
+	astar_pts[first_cum] = index_doer;
+	first_cum--;
 
-	r = point->r;
-	c = point->c;
+	seg_index_rc(alt_seg, index_doer, &r, &c);
 
-	G_debug(3, "R:%2d C:%2d, ", r, c);
+	G_debug(3, "A* Search: row %d, column %d, ", 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 @@ int do_astar(void)
 		    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 @@ int do_astar(void)
 			    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 @@ int do_astar(void)
     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,14 +112,9 @@ int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
     if (heap_size > do_points)
 	G_fatal_error(_("heapsize too large"));
 
-    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; */
+    heap_index[heap_size] = nxt_avail_pt++;
 
-    nxt_avail_pt++;
+    astar_pts[heap_size] = SEG_INDEX(alt_seg, r, c);
 
     /* sift up: move new point towards top of heap */
 
@@ -155,29 +145,18 @@ int drop_pt(void)
     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 @@ int drop_pt(void)
 	/* move hole down */
 
 	heap_index[parent] = heap_index[child];
+	astar_pts[parent] = astar_pts[child];
 	parent = child;
 
     }
@@ -199,14 +179,11 @@ int drop_pt(void)
     /* 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 @@ int drop_pt(void)
 /* 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) {
-
-	    /* push parent point down */
-	    heap_index[child] = parentp;
-	    child = parent;
-
-	}
-	/* same ele, but parent is younger */
-	else if (elep == ele && parentp > childp) {
 
+	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;
-
 	}
 	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;
@@ -286,6 +253,7 @@ int replace(			/* ele was in there */
 /* 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 @@ int replace(			/* ele was in there */
 
     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;

raster/r.watershed/ram/do_cum.c

@@ -11,6 +11,7 @@ int do_cum(void)
     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."));
 
@@ -19,13 +20,11 @@ int do_cum(void)
 	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)];
@@ -33,10 +32,11 @@ int do_cum(void)
 	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;
@@ -49,14 +49,14 @@ int do_cum(void)
 		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) {
@@ -68,7 +68,6 @@ int do_cum(void)
 	    }
 	}
     }
-    G_percent(count, do_points, 1);	/* finish it */
     G_free(astar_pts);
 
     return 0;
@@ -113,6 +112,7 @@ int do_cum_mfd(void)
     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);
@@ -143,13 +143,11 @@ int do_cum_mfd(void)
     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)];
@@ -157,8 +155,9 @@ int do_cum_mfd(void)
 	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;
@@ -171,7 +170,7 @@ int do_cum_mfd(void)
 	    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 */
@@ -184,17 +183,19 @@ int do_cum_mfd(void)
 		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 = Rast_is_c_null_value(&ele_nbr);
 			edge = is_null;
 			if (!is_null && ele_nbr <= ele) {
@@ -234,7 +235,7 @@ int do_cum_mfd(void)
 		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;
 	    }
@@ -268,11 +269,13 @@ int do_cum_mfd(void)
 			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];
@@ -285,7 +288,7 @@ int do_cum_mfd(void)
 				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) {
@@ -307,7 +310,7 @@ int do_cum_mfd(void)
 	    }
 
 	    if (mfd_cells < 2) {
-		valued = wat[SEG_INDEX(wat_seg, dr, dc)];
+		valued = wat[down_index];
 		if (value > 0) {
 		    if (valued > 0)
 			valued += value;
@@ -320,22 +323,22 @@ int do_cum_mfd(void)
 		    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 */
@@ -356,7 +359,6 @@ int do_cum_mfd(void)
 	    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);

raster/r.watershed/ram/init_vars.c

@@ -13,9 +13,10 @@ int init_vars(int argc, char *argv[])
     CELL *buf, alt_value, wat_value, asp_value, block_value;
     DCELL dvalue;
     void *elebuf, *ptr;
-    int fd, index, ele_map_type;
+    int fd, ele_map_type;
     size_t ele_size;
     char MASK_flag;
+    int seg_idx;
 
     G_gisinit(argv[0]);
     ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
@@ -159,7 +160,7 @@ int init_vars(int argc, char *argv[])
 	Rast_get_row(fd, elebuf, r, ele_map_type);
 	ptr = elebuf;
 	for (c = 0; c < ncols; c++) {
-	    index = SEG_INDEX(alt_seg, r, c);
+	    seg_idx = SEG_INDEX(alt_seg, r, c);
 
 	    /* all flags need to be manually set to zero */
 	    flag_unset(swale, r, c);
@@ -188,9 +189,9 @@ int init_vars(int argc, char *argv[])
 		    alt_value = ele_round(dvalue);
 		}
 	    }
-	    alt[index] = alt_value;
+	    alt[seg_idx] = alt_value;
 	    if (er_flag) {
-		r_h[index] = alt_value;
+		r_h[seg_idx] = alt_value;
 	    }
 	    ptr = G_incr_void_ptr(ptr, ele_size);
 	}
@@ -216,8 +217,8 @@ int init_vars(int argc, char *argv[])
 	    Rast_get_c_row(fd, buf, r);
 	    for (c = 0; c < ncols; c++) {
 		if (MASK_flag) {
-		    index = FLAG_GET(worked, r, c);
-		    if (!index)
+		    block_value = FLAG_GET(worked, r, c);
+		    if (!block_value)
 			wat[SEG_INDEX(wat_seg, r, c)] = buf[c];
 		    else
 			wat[SEG_INDEX(wat_seg, r, c)] = 0.0;
@@ -233,8 +234,8 @@ int init_vars(int argc, char *argv[])
 	for (r = 0; r < nrows; r++) {
 	    for (c = 0; c < ncols; c++) {
 		if (MASK_flag) {
-		    index = FLAG_GET(worked, r, c);
-		    if (!index)
+		    block_value = FLAG_GET(worked, r, c);
+		    if (!block_value)
 			wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
 		}
 		else
@@ -305,7 +306,7 @@ int init_vars(int argc, char *argv[])
 			       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 */
@@ -321,22 +322,23 @@ int init_vars(int argc, char *argv[])
 	for (r = 0; r < nrows; r++) {
 	    G_percent(r, nrows, 3);
 	    for (c = 0; c < ncols; c++) {
+		seg_idx = SEG_INDEX(wat_seg, r, c);
 		if (FLAG_GET(worked, r, c)) {
-		    wat[SEG_INDEX(wat_seg, r, c)] = 0;
+		    wat[seg_idx] = 0;
 		}
 		else {
 		    if (er_flag)
-			s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
-		    asp_value = asp[SEG_INDEX(asp_seg, r, c)];
+			s_l[seg_idx] = half_res;
+		    asp_value = asp[seg_idx];
 		    if (r == 0 || c == 0 || r == nrows - 1 ||
 			c == ncols - 1 || asp_value != 0) {
-			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+			wat_value = wat[seg_idx];
 			if (wat_value > 0)
-			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			    wat[seg_idx] = -wat_value;
 			/* set depression */
 			if (asp_value) {
 			    asp_value = 0;
-			    wat[SEG_INDEX(wat_seg, r, c)] = ABS(wat_value);
+			    wat[seg_idx] = ABS(wat_value);
 			}
 			else if (r == 0)
 			    asp_value = -2;
@@ -346,73 +348,73 @@ int init_vars(int argc, char *argv[])
 			    asp_value = -6;
 			else if (c == ncols - 1)
 			    asp_value = -8;
-			asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
-			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+			asp[seg_idx] = asp_value;
+			alt_value = alt[seg_idx];
 			add_pt(r, c, alt_value, alt_value);
 		    }
 		    else if (FLAG_GET(worked, r - 1, c)) {
-			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+			alt_value = alt[seg_idx];
 			add_pt(r, c, alt_value, alt_value);
-			asp[SEG_INDEX(asp_seg, r, c)] = -2;
-			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+			asp[seg_idx] = -2;
+			wat_value = wat[seg_idx];
 			if (wat_value > 0)
-			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			    wat[seg_idx] = -wat_value;
 		    }
 		    else if (FLAG_GET(worked, r + 1, c)) {
-			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+			alt_value = alt[seg_idx];
 			add_pt(r, c, alt_value, alt_value);
-			asp[SEG_INDEX(asp_seg, r, c)] = -6;
-			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+			asp[seg_idx] = -6;
+			wat_value = wat[seg_idx];
 			if (wat_value > 0)
-			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			    wat[seg_idx] = -wat_value;
 		    }
 		    else if (FLAG_GET(worked, r, c - 1)) {
-			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+			alt_value = alt[seg_idx];
 			add_pt(r, c, alt_value, alt_value);
-			asp[SEG_INDEX(asp_seg, r, c)] = -4;
-			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+			asp[seg_idx] = -4;
+			wat_value = wat[seg_idx];
 			if (wat_value > 0)
-			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			    wat[seg_idx] = -wat_value;
 		    }
 		    else if (FLAG_GET(worked, r, c + 1)) {
-			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+			alt_value = alt[seg_idx];
 			add_pt(r, c, alt_value, alt_value);
-			asp[SEG_INDEX(asp_seg, r, c)] = -8;
-			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+			asp[seg_idx] = -8;
+			wat_value = wat[seg_idx];
 			if (wat_value > 0)
-			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			    wat[seg_idx] = -wat_value;
 		    }
 		    else if (sides == 8 && FLAG_GET(worked, r - 1, c - 1)) {
-			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+			alt_value = alt[seg_idx];
 			add_pt(r, c, alt_value, alt_value);
-			asp[SEG_INDEX(asp_seg, r, c)] = -3;
-			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+			asp[seg_idx] = -3;
+			wat_value = wat[seg_idx];
 			if (wat_value > 0)
-			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			    wat[seg_idx] = -wat_value;
 		    }
 		    else if (sides == 8 && FLAG_GET(worked, r - 1, c + 1)) {
-			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+			alt_value = alt[seg_idx];
 			add_pt(r, c, alt_value, alt_value);
-			asp[SEG_INDEX(asp_seg, r, c)] = -1;
-			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+			asp[seg_idx] = -1;
+			wat_value = wat[seg_idx];
 			if (wat_value > 0)
-			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			    wat[seg_idx] = -wat_value;
 		    }
 		    else if (sides == 8 && FLAG_GET(worked, r + 1, c - 1)) {
-			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+			alt_value = alt[seg_idx];
 			add_pt(r, c, alt_value, alt_value);
-			asp[SEG_INDEX(asp_seg, r, c)] = -5;
-			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+			asp[seg_idx] = -5;
+			wat_value = wat[seg_idx];
 			if (wat_value > 0)
-			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			    wat[seg_idx] = -wat_value;
 		    }
 		    else if (sides == 8 && FLAG_GET(worked, r + 1, c + 1)) {
-			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+			alt_value = alt[seg_idx];
 			add_pt(r, c, alt_value, alt_value);
-			asp[SEG_INDEX(asp_seg, r, c)] = -7;
-			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+			asp[seg_idx] = -7;
+			wat_value = wat[seg_idx];
 			if (wat_value > 0)
-			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			    wat[seg_idx] = -wat_value;
 		    }
 		}
 	    }
@@ -422,19 +424,20 @@ int init_vars(int argc, char *argv[])
 	for (r = 0; r < nrows; r++) {
 	    G_percent(r, nrows, 3);
 	    for (c = 0; c < ncols; c++) {
+		seg_idx = SEG_INDEX(wat_seg, r, c);
 		if (er_flag)
-		    s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
-		asp_value = asp[SEG_INDEX(asp_seg, r, c)];
+		    s_l[seg_idx] = half_res;
+		asp_value = asp[seg_idx];
 		if (r == 0 || c == 0 || r == nrows - 1 ||
 		    c == ncols - 1 || asp_value != 0) {
-		    wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+		    wat_value = wat[seg_idx];
 		    if (wat_value > 0) {
-			wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+			wat[seg_idx] = -wat_value;
 		    }
 		    /* set depression */
 		    if (asp_value) {
 			asp_value = 0;
-			wat[SEG_INDEX(wat_seg, r, c)] = ABS(wat_value);
+			wat[seg_idx] = ABS(wat_value);
 		    }
 		    else if (r == 0)
 			asp_value = -2;
@@ -444,8 +447,8 @@ int init_vars(int argc, char *argv[])
 			asp_value = -6;
 		    else if (c == ncols - 1)
 			asp_value = -8;
-		    asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
-		    alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+		    asp[seg_idx] = asp_value;
+		    alt_value = alt[seg_idx];
 		    add_pt(r, c, alt_value, alt_value);
 		}
 	    }

raster/r.watershed/ram/main.c

@@ -34,7 +34,7 @@ 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;
 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;

raster/r.watershed/ram/ramseg.c

@@ -13,3 +13,16 @@ int size_array(int *ram_seg, int nrows, int ncols)
     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;
+}