[GRASS-SVN] r39727 -
grass/branches/develbranch_6/raster/r.watershed/seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Nov 14 06:30:00 EST 2009
Author: mmetz
Date: 2009-11-14 06:30:00 -0500 (Sat, 14 Nov 2009)
New Revision: 39727
Modified:
grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c
Log:
fix for ticket #807 backported from trunk
Modified: grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c 2009-11-14 11:24:33 UTC (rev 39726)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c 2009-11-14 11:30:00 UTC (rev 39727)
@@ -5,20 +5,45 @@
#include "Gwater.h"
#include "do_astar.h"
+double get_slope2(CELL, CELL, double);
+
int do_astar(void)
{
POINT point;
int doer, count;
SHORT upr, upc, r, c, ct_dir;
- CELL work_val, alt_val, alt_up, asp_up;
+ CELL work_val, alt_val, alt_nbr[8], alt_up, asp_up;
DCELL wat_val;
CELL in_val, drain_val;
HEAP heap_pos;
+ /* sides
+ * |7|1|4|
+ * |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 };
+ double dx, dy, dist_to_nbr[8], ew_res, ns_res;
- /* double slope; */
+ double slope[8];
G_message(_("SECTION 2: A * Search."));
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (r_nbr, c_nbr) for neighbours */
+ upr = nextdr[ct_dir];
+ upc = nextdc[ct_dir];
+ /* account for rare cases when ns_res != ew_res */
+ dy = ABS(upr) * window.ns_res;
+ dx = ABS(upc) * window.ew_res;
+ if (ct_dir < 4)
+ dist_to_nbr[ct_dir] = dx + dy;
+ else
+ dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
+ }
+ ew_res = window.ew_res;
+ ns_res = window.ns_res;
+
count = 0;
seg_get(&heap_index, (char *)&heap_pos, 0, 1);
first_astar = heap_pos.point;
@@ -54,11 +79,37 @@
/* get r, c (upr, upc) for this neighbour */
upr = r + nextdr[ct_dir];
upc = c + nextdc[ct_dir];
+ slope[ct_dir] = alt_nbr[ct_dir] = 0;
/* check that upr, upc are within region */
if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
- /* check if neighbour is in the list */
- /* if not, add as new point */
+ /* slope */
bseg_get(&in_list, &in_val, upr, upc);
+ bseg_get(&worked, &work_val, upr, upc);
+ if (!work_val) {
+ cseg_get(&alt, &alt_up, upr, upc);
+ alt_nbr[ct_dir] = alt_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 to ew nbr > slope to center */
+ if (slope[ct_dir] <=
+ get_slope2(alt_nbr[nbr_ew[ct_dir]],
+ alt_nbr[ct_dir], ew_res))
+ in_val = 1;
+ }
+ if (!in_val && slope[nbr_ns[ct_dir]]) {
+ /* 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))
+ in_val = 1;
+ }
+ }
+ }
+ /* put neighbour in search list if not yet in */
if (in_val == 0) {
cseg_get(&alt, &alt_up, upr, upc);
add_pt(upr, upc, alt_up, alt_val);
@@ -69,7 +120,6 @@
else {
/* check if neighbour has not been worked on,
* update values for asp and wat */
- bseg_get(&worked, &work_val, upr, upc);
if (!work_val) {
cseg_get(&asp, &asp_up, upr, upc);
if (asp_up < 0) {
@@ -80,15 +130,11 @@
wat_val = -wat_val;
dseg_put(&wat, &wat_val, r, c);
}
- /* cseg_get(&alt, &alt_up, upr, upc); */
- /* replace(upr, upc, r, c); */ /* alt_up used to be */
- /* slope = get_slope (upr, upc, r, c, alt_up, alt_val);
- dseg_put (&slp, &slope, upr, upc); */
}
}
}
}
- }
+ } /* end sides */
}
if (mfd == 0)
@@ -107,10 +153,6 @@
POINT point;
HEAP heap_pos;
- /* double slp_value; */
-
- /* slp_value = get_slope(r, c, downr, downc, ele, downe);
- dseg_put (&slp, &slp_value, r, c); */
bseg_put(&in_list, &one, r, c);
/* add point to next free position */
@@ -126,8 +168,6 @@
point.r = r;
point.c = c;
-/* point.downr = downr;
- point.downc = downc; */
seg_put(&astar_pts, (char *)&point, 0, nxt_avail_pt);
@@ -146,8 +186,8 @@
int drop_pt(void)
{
int child, childr, parent;
- int childp; /* , childrp */
- CELL ele; /* , eler */
+ int childp;
+ CELL ele;
int i;
HEAP heap_pos;
@@ -177,8 +217,6 @@
i = child + 3;
while (childr <= heap_size && childr < i) {
seg_get(&heap_index, (char *)&heap_pos, 0, childr);
- /* childrp = heap_pos.point;
- eler = heap_pos.ele; */
if (heap_pos.ele < ele) {
child = childr;
childp = heap_pos.point;
@@ -190,7 +228,6 @@
childp = heap_pos.point;
}
childr++;
- /* i++; */
}
}
@@ -211,7 +248,6 @@
/* sift up last swapped point, only necessary if hole moved to heap end */
sift_up(parent, ele);
-
}
/* the actual drop */
@@ -224,8 +260,7 @@
/* standard sift-up routine for d-ary min heap */
int sift_up(int start, CELL ele)
{
- int parent, child, childp; /* parentp, */
- /* CELL elep; */
+ int parent, child, childp;
HEAP heap_pos;
child = start;
@@ -235,8 +270,6 @@
while (child > 1) {
parent = GET_PARENT(child);
seg_get(&heap_index, (char *)&heap_pos, 0, parent);
- /* parentp = heap_pos.point;
- elep = heap_pos.ele; */
/* parent ele higher */
if (heap_pos.ele > ele) {
@@ -244,7 +277,6 @@
/* push parent point down */
seg_put(&heap_index, (char *)&heap_pos, 0, child);
child = parent;
-
}
/* same ele, but parent is younger */
else if (heap_pos.ele == ele && heap_pos.point > childp) {
@@ -252,7 +284,6 @@
/* push parent point down */
seg_put(&heap_index, (char *)&heap_pos, 0, child);
child = parent;
-
}
else
/* no more sifting up, found new slot for child */
@@ -284,6 +315,14 @@
return (slope);
}
+double get_slope2(CELL ele, CELL up_ele, double dist)
+{
+ if (ele == up_ele)
+ return 0.5 / dist;
+ else
+ return (double)(up_ele - ele) / dist;
+}
+
/* no longer needed */
int replace(SHORT upr, SHORT upc, SHORT r, SHORT c)
{ /* ele was in there */
@@ -299,7 +338,7 @@
G_warning("pnt_index incorrect!");
return 1;
}
-/* point.downr = r;
+ /* point.downr = r;
point.downc = c; */
seg_put(&astar_pts, (char *)&point, 0, now);
More information about the grass-commit
mailing list