1 #include "Gwater.h"
2 #include "do_astar.h"
3 #include <grass/gis.h>
4 #include <grass/glocale.h>
5
6 static double get_slope2(CELL, CELL, double);
7
do_astar(void)8 int do_astar(void)
9 {
10 int count;
11 int upr, upc, r, c, ct_dir;
12 CELL alt_val, alt_nbr[8];
13 CELL is_in_list, is_worked, flat_is_done, nbr_flat_is_done;
14 int index_doer, index_up;
15
16 /* sides
17 * |7|1|4|
18 * |2| |3|
19 * |5|0|6|
20 */
21 int nbr_ew[8] = { 0, 1, 2, 3, 1, 0, 0, 1 };
22 int nbr_ns[8] = { 0, 1, 2, 3, 3, 2, 3, 2 };
23 double dx, dy, dist_to_nbr[8], ew_res, ns_res;
24 double slope[8];
25 int skip_diag;
26 CELL *alt_bak;
27
28 G_message(_("SECTION 2: A* Search."));
29
30 for (ct_dir = 0; ct_dir < sides; ct_dir++) {
31 /* get r, c (upr, upc) for neighbours */
32 upr = nextdr[ct_dir];
33 upc = nextdc[ct_dir];
34 /* account for rare cases when ns_res != ew_res */
35 dy = ABS(upr) * window.ns_res;
36 dx = ABS(upc) * window.ew_res;
37 if (ct_dir < 4)
38 dist_to_nbr[ct_dir] = dx + dy;
39 else
40 dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
41 }
42 ew_res = window.ew_res;
43 ns_res = window.ns_res;
44
45 count = 0;
46 first_astar = heap_index[1];
47 first_cum = do_points;
48
49 if (flat_flag) {
50 alt_bak =
51 (CELL *) G_malloc(sizeof(CELL) *
52 size_array(&alt_seg, nrows, ncols));
53 flat_done = flag_create(nrows, ncols);
54 flat_is_done = 0;
55
56 for (r = 0; r < nrows; r++) {
57 for (c = 0; c < ncols; c++) {
58 index_doer = SEG_INDEX(alt_seg, r, c);
59 alt_bak[index_doer] = alt[index_doer];
60
61 flag_unset(flat_done, r, c);
62 }
63 }
64 }
65 else {
66 alt_bak = NULL;
67 flat_done = NULL;
68 flat_is_done = 1;
69 }
70
71 /* A* Search: search uphill, get downhill paths */
72 while (heap_size > 0) {
73 G_percent(count++, do_points, 1);
74
75 /* get point with lowest elevation, in case of equal elevation
76 * of following points, oldest point = point added earliest */
77 index_doer = astar_pts[1];
78
79 drop_pt();
80
81 /* add astar points to sorted list for flow accumulation */
82 astar_pts[first_cum] = index_doer;
83 first_cum--;
84
85 seg_index_rc(alt_seg, index_doer, &r, &c);
86
87 G_debug(3, "A* Search: row %d, column %d, ", r, c);
88
89 alt_val = alt[index_doer];
90
91 if (flat_flag) {
92 flat_is_done = FLAG_GET(flat_done, r, c);
93 }
94
95 /* check all neighbours, breadth first search */
96 for (ct_dir = 0; ct_dir < sides; ct_dir++) {
97 /* get r, c (upr, upc) for this neighbour */
98 upr = r + nextdr[ct_dir];
99 upc = c + nextdc[ct_dir];
100 slope[ct_dir] = -1;
101 alt_nbr[ct_dir] = 0;
102 /* check if r, c are within region */
103 if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
104 index_up = SEG_INDEX(alt_seg, upr, upc);
105 is_in_list = FLAG_GET(in_list, upr, upc);
106 is_worked = FLAG_GET(worked, upr, upc);
107 skip_diag = 0;
108
109 alt_nbr[ct_dir] = alt[index_up];
110 if (flat_flag && !is_in_list && !is_worked) {
111 alt_val = alt_bak[index_doer];
112 alt_nbr[ct_dir] = alt_bak[index_up];
113 if (!flat_is_done && alt_nbr[ct_dir] == alt_val) {
114 do_flatarea(index_doer, alt_val, alt_bak, alt);
115 alt_nbr[ct_dir] = alt[index_up];
116 flat_is_done = 1;
117 nbr_flat_is_done = 1;
118 }
119 nbr_flat_is_done = FLAG_GET(flat_done, upr, upc);
120 if (!nbr_flat_is_done) {
121 /* use original ele values */
122 alt_val = alt_bak[index_doer];
123 alt_nbr[ct_dir] = alt_bak[index_up];
124 }
125 else {
126 /* use modified ele values */
127 alt_val = alt[index_doer];
128 alt_nbr[ct_dir] = alt[index_up];
129 }
130 }
131
132 /* avoid diagonal flow direction bias */
133 if (!is_worked) {
134 slope[ct_dir] =
135 get_slope2(alt_val, alt_nbr[ct_dir],
136 dist_to_nbr[ct_dir]);
137 }
138 if (!is_in_list || (!is_worked && asp[index_up] < 0)) {
139 if (ct_dir > 3 && slope[ct_dir] > 0) {
140 if (slope[nbr_ew[ct_dir]] >= 0) {
141 /* slope to ew nbr > slope to center */
142 if (slope[ct_dir] <
143 get_slope2(alt_nbr[nbr_ew[ct_dir]],
144 alt_nbr[ct_dir], ew_res))
145 skip_diag = 1;
146 }
147 if (!skip_diag && slope[nbr_ns[ct_dir]] >= 0) {
148 /* slope to ns nbr > slope to center */
149 if (slope[ct_dir] <
150 get_slope2(alt_nbr[nbr_ns[ct_dir]],
151 alt_nbr[ct_dir], ns_res))
152 skip_diag = 1;
153 }
154 }
155 }
156
157 if (!skip_diag) {
158 /* add neighbour as new point if not in the list */
159 if (!is_in_list) {
160 add_pt(upr, upc, alt_nbr[ct_dir]);
161 /* set flow direction */
162 asp[index_up] = drain[upr - r + 1][upc - c + 1];
163 }
164 else if (!is_worked) {
165 /* neighbour is edge in list, not yet worked */
166 if (asp[index_up] < 0 && slope[ct_dir] > 0) {
167 /* adjust flow direction for edge cell */
168 asp[index_up] = drain[upr - r + 1][upc - c + 1];
169
170 if (wat[index_doer] > 0)
171 wat[index_doer] = -1.0 * wat[index_doer];
172 }
173 /* neighbour is inside real depression, not yet worked */
174 else if (asp[index_up] == 0) {
175 asp[index_up] = drain[upr - r + 1][upc - c + 1];
176 }
177 }
178 }
179 } /* end if in region */
180 } /* end sides */
181 FLAG_SET(worked, r, c);
182 }
183 G_percent(count, do_points, 1); /* finish it */
184 if (mfd == 0)
185 flag_destroy(worked);
186
187 flag_destroy(in_list);
188 G_free(heap_index);
189
190 if (flat_flag) {
191 for (r = 0; r < nrows; r++) {
192 for (c = 0; c < ncols; c++) {
193 index_doer = SEG_INDEX(alt_seg, r, c);
194 alt[index_doer] = alt_bak[index_doer];
195 }
196 }
197 G_free(alt_bak);
198 flag_destroy(flat_done);
199 }
200
201 return 0;
202 }
203
204 /* compare two heap points */
205 /* return 1 if a < b else 0 */
cmp_pnt(CELL elea,CELL eleb,int addeda,int addedb)206 static int cmp_pnt(CELL elea, CELL eleb, int addeda, int addedb)
207 {
208 if (elea == eleb) {
209 return (addeda < addedb);
210 }
211 return (elea < eleb);
212 }
213
214 /* standard sift-up routine for d-ary min heap */
sift_up(int start,CELL ele)215 static int sift_up(int start, CELL ele)
216 {
217 register int parent, child, child_idx, child_added;
218 CELL elep;
219
220 child = start;
221 child_added = heap_index[child];
222 child_idx = astar_pts[child];
223
224 while (child > 1) {
225 GET_PARENT(parent, child);
226
227 elep = alt[astar_pts[parent]];
228 /* child smaller */
229 if (cmp_pnt(ele, elep, child_added, heap_index[parent])) {
230 /* push parent point down */
231 heap_index[child] = heap_index[parent];
232 astar_pts[child] = astar_pts[parent];
233 child = parent;
234 }
235 else
236 /* no more sifting up, found new slot for child */
237 break;
238 }
239
240 /* put point in new slot */
241 if (child < start) {
242 heap_index[child] = child_added;
243 astar_pts[child] = child_idx;
244 }
245
246 return 0;
247 }
248
249 /* add point routine for min heap */
add_pt(int r,int c,CELL ele)250 int add_pt(int r, int c, CELL ele)
251 {
252 FLAG_SET(in_list, r, c);
253
254 /* add point to next free position */
255 heap_size++;
256
257 if (heap_size > do_points)
258 G_fatal_error(_("heapsize too large"));
259
260 heap_index[heap_size] = nxt_avail_pt++;
261 astar_pts[heap_size] = SEG_INDEX(alt_seg, r, c);
262
263 /* sift up: move new point towards top of heap */
264 sift_up(heap_size, ele);
265
266 return 0;
267 }
268
269 /* drop point routine for min heap */
drop_pt(void)270 int drop_pt(void)
271 {
272 register int child, childr, parent;
273 CELL ele, eler;
274 register int i;
275
276 if (heap_size == 1) {
277 heap_index[1] = -1;
278 heap_size = 0;
279 return 0;
280 }
281
282 /* start with root */
283 parent = 1;
284
285 /* sift down: move hole back towards bottom of heap */
286
287 while (GET_CHILD(child, parent) <= heap_size) {
288 /* select child with lower ele, if both are equal, older child
289 * older child is older startpoint for flow path, important */
290 ele = alt[astar_pts[child]];
291 i = child + 3;
292 for (childr = child + 1; childr <= heap_size && childr < i; childr++) {
293 eler = alt[astar_pts[childr]];
294 if (cmp_pnt(eler, ele, heap_index[childr], heap_index[child])) {
295 child = childr;
296 ele = eler;
297 }
298 }
299
300 /* move hole down */
301 heap_index[parent] = heap_index[child];
302 astar_pts[parent] = astar_pts[child];
303 parent = child;
304 }
305
306 /* hole is in lowest layer, move to heap end */
307 if (parent < heap_size) {
308 heap_index[parent] = heap_index[heap_size];
309 astar_pts[parent] = astar_pts[heap_size];
310
311 ele = alt[astar_pts[parent]];
312 /* sift up last swapped point, only necessary if hole moved to heap end */
313 sift_up(parent, ele);
314 }
315
316 /* the actual drop */
317 heap_size--;
318
319 return 0;
320 }
321
get_slope(int r,int c,int downr,int downc,CELL ele,CELL downe)322 double get_slope(int r, int c, int downr, int downc, CELL ele, CELL downe)
323 {
324 double slope;
325
326 if (r == downr)
327 slope = (ele - downe) / window.ew_res;
328 else if (c == downc)
329 slope = (ele - downe) / window.ns_res;
330 else
331 slope = (ele - downe) / diag;
332
333 if (slope < MIN_SLOPE)
334 return (MIN_SLOPE);
335
336 return (slope);
337 }
338
get_slope2(CELL ele,CELL up_ele,double dist)339 static double get_slope2(CELL ele, CELL up_ele, double dist)
340 {
341 if (ele >= up_ele)
342 return 0.0;
343 else
344 return (double)(up_ele - ele) / dist;
345 }
346