1
2 /**********************************************************************
3 *
4 * flat area beautification after Garbrecht and Martz (1997)
5 *
6 * Garbrecht, J. & Martz, L. W. (1997)
7 * The assignment of drainage direction over flat surfaces in raster
8 * digital elevation models. J. Hydrol 193, 204-213.
9 *
10 * the method is modified for speed, only one pass is necessary to get
11 * the gradient away from higher terrain
12 *
13 *********************************************************************/
14
15
16 #include <limits.h>
17 #include <assert.h>
18 #include <grass/gis.h>
19 #include <grass/glocale.h>
20 #include <grass/rbtree.h>
21 #include "Gwater.h"
22 #include "do_astar.h"
23
24 struct pq_node
25 {
26 int idx;
27 struct pq_node *next;
28 };
29
30 struct pq
31 {
32 struct pq_node *first, *last;
33 int size;
34 };
35
pq_create(void)36 struct pq *pq_create(void)
37 {
38 struct pq *q = G_malloc(sizeof(struct pq));
39
40 q->first = G_malloc(sizeof(struct pq_node));
41 q->first->next = NULL;
42 q->first->idx = -1;
43 q->last = q->first;
44 q->size = 0;
45
46 return q;
47 }
48
49 /* dummy end must always be allocated and empty */
pq_add(int idx,struct pq * q)50 int pq_add(int idx, struct pq *q)
51 {
52 assert(q->last);
53 assert(q->last->idx == -1);
54
55 q->last->idx = idx;
56 if (q->last->next != NULL) {
57 G_fatal_error(_("Beautify flat areas: priority queue error"));
58 }
59
60 struct pq_node *n = (struct pq_node *)G_malloc(sizeof(struct pq_node));
61
62 n->next = NULL;
63 n->idx = -1;
64 q->last->next = n;
65 q->last = q->last->next;
66
67 assert(q->last != q->last->next);
68 assert(q->first != q->last);
69 q->size++;
70
71 return 0;
72 }
73
pq_drop(struct pq * q)74 int pq_drop(struct pq *q)
75 {
76 int idx = q->first->idx;
77 struct pq_node *n = q->first;
78
79 q->size--;
80
81 q->first = q->first->next;
82 assert(q->first);
83 assert(q->first != q->first->next);
84 assert(n != q->first);
85 G_free(n);
86
87 return idx;
88 }
89
pq_destroy(struct pq * q)90 int pq_destroy(struct pq *q)
91 {
92 struct pq_node *delme;
93
94 while (q->first) {
95 delme = q->first;
96 q->first = q->first->next;
97 G_free(delme);
98 }
99
100 G_free(q);
101
102 return 0;
103 }
104
105 struct orders
106 {
107 int index, uphill, downhill;
108 char flag;
109 };
110
cmp_orders(const void * a,const void * b)111 int cmp_orders(const void *a, const void *b)
112 {
113 struct orders *oa = (struct orders *)a;
114 struct orders *ob = (struct orders *)b;
115
116 return (oa->index < ob->index ? -1 : (oa->index > ob->index));
117 }
118
119 /*
120 * return 0 if nothing was modidied
121 * return 1 if elevation was modified
122 */
do_flatarea(int index,CELL ele,CELL * alt_org,CELL * alt_new)123 int do_flatarea(int index, CELL ele, CELL * alt_org, CELL * alt_new)
124 {
125 int upr, upc, r, c, ct_dir;
126 CELL is_in_list, is_worked, this_in_list;
127 int index_doer, index_up;
128 int n_flat_cells = 0, counter;
129 CELL ele_nbr, min_ele_diff;
130 int uphill_order, downhill_order, max_uphill_order, max_downhill_order;
131 int last_order;
132
133 struct pq *up_pq = pq_create();
134 struct pq *down_pq = pq_create();
135
136 struct orders inc_order, *order_found, *nbr_order_found;
137 struct RB_TREE *order_tree =
138 rbtree_create(cmp_orders, sizeof(struct orders));
139
140 pq_add(index, down_pq);
141 pq_add(index, up_pq);
142 inc_order.downhill = -1;
143 inc_order.uphill = 0;
144 inc_order.index = index;
145 inc_order.flag = 0;
146 rbtree_insert(order_tree, &inc_order);
147
148 n_flat_cells = 1;
149
150 min_ele_diff = INT_MAX;
151 max_uphill_order = max_downhill_order = 0;
152
153 /* get uphill start points */
154 G_debug(2, "get uphill start points");
155 counter = 0;
156 while (down_pq->size) {
157 if ((index_doer = pq_drop(down_pq)) == -1)
158 G_fatal_error("get start points: no more points in down queue");
159
160 seg_index_rc(alt_seg, index_doer, &r, &c);
161
162 FLAG_SET(flat_done, r, c);
163
164 /* check all neighbours, breadth first search */
165 for (ct_dir = 0; ct_dir < sides; ct_dir++) {
166 /* get r, c (upr, upc) for this neighbour */
167 upr = r + nextdr[ct_dir];
168 upc = c + nextdc[ct_dir];
169 /* check if r, c are within region */
170 if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
171 index_up = SEG_INDEX(alt_seg, upr, upc);
172 is_in_list = FLAG_GET(in_list, upr, upc);
173 is_worked = FLAG_GET(worked, upr, upc);
174 ele_nbr = alt_org[index_up];
175
176 if (ele_nbr == ele && !is_worked) {
177
178 inc_order.downhill = -1;
179 inc_order.uphill = -1;
180 inc_order.index = index_up;
181 inc_order.flag = 0;
182 /* not yet added to queue */
183 if ((order_found =
184 rbtree_find(order_tree, &inc_order)) == NULL) {
185 n_flat_cells++;
186
187 /* add to down queue if not yet in there */
188 pq_add(index_up, down_pq);
189
190 /* add to up queue if not yet in there */
191 if (is_in_list) {
192 pq_add(index_up, up_pq);
193 /* set uphill order to 0 */
194 inc_order.uphill = 0;
195 counter++;
196 }
197 rbtree_insert(order_tree, &inc_order);
198 }
199 }
200 }
201 }
202 }
203 /* flat area too small, not worth the effort */
204 if (n_flat_cells < 5) {
205 /* clean up */
206 pq_destroy(up_pq);
207 pq_destroy(down_pq);
208 rbtree_destroy(order_tree);
209
210 return 0;
211 }
212
213 G_debug(2, "%d flat cells, %d cells in tree, %d start cells",
214 n_flat_cells, (int)order_tree->count, counter);
215
216 pq_destroy(down_pq);
217 down_pq = pq_create();
218
219 /* got uphill start points, do uphill correction */
220 G_debug(2, "got uphill start points, do uphill correction");
221 counter = 0;
222 uphill_order = 1;
223 while (up_pq->size) {
224 int is_in_down_queue = 0;
225
226 if ((index_doer = pq_drop(up_pq)) == -1)
227 G_fatal_error("uphill order: no more points in up queue");
228
229 seg_index_rc(alt_seg, index_doer, &r, &c);
230 this_in_list = FLAG_GET(in_list, r, c);
231
232 /* get uphill order for this point */
233 inc_order.index = index_doer;
234 if ((order_found = rbtree_find(order_tree, &inc_order)) == NULL)
235 G_fatal_error(_("flat cell escaped for uphill correction"));
236
237 last_order = uphill_order - 1;
238 uphill_order = order_found->uphill;
239
240 if (last_order > uphill_order)
241 G_warning(_("queue error: last uphill order %d > current uphill order %d"),
242 last_order, uphill_order);
243
244 /* debug */
245 if (uphill_order == -1)
246 G_fatal_error(_("uphill order not set"));
247
248 if (max_uphill_order < uphill_order)
249 max_uphill_order = uphill_order;
250
251 uphill_order++;
252 counter++;
253
254 /* check all neighbours, breadth first search */
255 for (ct_dir = 0; ct_dir < sides; ct_dir++) {
256 /* get r, c (upr, upc) for this neighbour */
257 upr = r + nextdr[ct_dir];
258 upc = c + nextdc[ct_dir];
259 /* check if r, c are within region */
260 if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
261 index_up = SEG_INDEX(alt_seg, upr, upc);
262 is_in_list = FLAG_GET(in_list, upr, upc);
263 is_worked = FLAG_GET(worked, upr, upc);
264 ele_nbr = alt_org[index_up];
265
266 /* all cells that are in_list should have been added
267 * previously as uphill start points */
268 if (ele_nbr == ele && !is_worked) {
269
270 inc_order.index = index_up;
271 if ((nbr_order_found =
272 rbtree_find(order_tree, &inc_order)) == NULL) {
273 G_fatal_error(_("flat cell escaped in uphill correction"));
274 }
275
276 /* not yet added to queue */
277 if (nbr_order_found->uphill == -1) {
278 if (is_in_list)
279 G_warning("cell should be in queue");
280 /* add to up queue */
281 pq_add(index_up, up_pq);
282 /* set nbr uphill order = current uphill order + 1 */
283 nbr_order_found->uphill = uphill_order;
284 }
285 }
286 /* add focus cell to down queue */
287 if (!this_in_list && !is_in_down_queue &&
288 ele_nbr != ele && !is_in_list && !is_worked) {
289 pq_add(index_doer, down_pq);
290 /* set downhill order to 0 */
291 order_found->downhill = 0;
292 is_in_down_queue = 1;
293 }
294 if (ele_nbr > ele && min_ele_diff > ele_nbr - ele)
295 min_ele_diff = ele_nbr - ele;
296 }
297 }
298 }
299 /* debug: all flags should be set to 0 */
300
301 pq_destroy(up_pq);
302 up_pq = pq_create();
303
304 /* got downhill start points, do downhill correction */
305 G_debug(2, "got downhill start points, do downhill correction");
306 downhill_order = 1;
307 while (down_pq->size) {
308 if ((index_doer = pq_drop(down_pq)) == -1)
309 G_fatal_error(_("downhill order: no more points in down queue"));
310
311 seg_index_rc(alt_seg, index_doer, &r, &c);
312 this_in_list = FLAG_GET(in_list, r, c);
313
314 /* get downhill order for this point */
315 inc_order.index = index_doer;
316 if ((order_found = rbtree_find(order_tree, &inc_order)) == NULL)
317 G_fatal_error(_("flat cell escaped for downhill correction"));
318
319 last_order = downhill_order - 1;
320 downhill_order = order_found->downhill;
321
322 if (last_order > downhill_order)
323 G_warning(_("queue error: last downhill order %d > current downhill order %d"),
324 last_order, downhill_order);
325
326 /* debug */
327 if (downhill_order == -1)
328 G_fatal_error(_("downhill order: downhill order not set"));
329
330 if (max_downhill_order < downhill_order)
331 max_downhill_order = downhill_order;
332
333 downhill_order++;
334
335 /* check all neighbours, breadth first search */
336 for (ct_dir = 0; ct_dir < sides; ct_dir++) {
337 /* get r, c (upr, upc) for this neighbour */
338 upr = r + nextdr[ct_dir];
339 upc = c + nextdc[ct_dir];
340 /* check if r, c are within region */
341 if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
342 index_up = SEG_INDEX(alt_seg, upr, upc);
343 is_in_list = FLAG_GET(in_list, upr, upc);
344 is_worked = FLAG_GET(worked, upr, upc);
345 ele_nbr = alt_org[index_up];
346
347 if (ele_nbr == ele && !is_worked) {
348
349 inc_order.index = index_up;
350 if ((nbr_order_found =
351 rbtree_find(order_tree, &inc_order)) == NULL)
352 G_fatal_error(_("flat cell escaped in downhill correction"));
353
354 /* not yet added to queue */
355 if (nbr_order_found->downhill == -1) {
356
357 /* add to down queue */
358 pq_add(index_up, down_pq);
359 /* set nbr downhill order = current downhill order + 1 */
360 nbr_order_found->downhill = downhill_order;
361
362 /* add to up queue */
363 if (is_in_list) {
364 pq_add(index_up, up_pq);
365 /* set flag */
366 nbr_order_found->flag = 1;
367 }
368 }
369 }
370 }
371 }
372 }
373
374 /* got uphill and downhill order, adjust ele */
375
376 /* increment: ele += uphill_order + max_downhill_order - downhill_order */
377 /* decrement: ele += uphill_order - max_uphill_order - downhill_order */
378
379 G_debug(2, "adjust ele");
380 while (up_pq->size) {
381 if ((index_doer = pq_drop(up_pq)) == -1)
382 G_fatal_error("no more points in up queue");
383
384 seg_index_rc(alt_seg, index_doer, &r, &c);
385 this_in_list = FLAG_GET(in_list, r, c);
386
387 /* get uphill and downhill order for this point */
388 inc_order.index = index_doer;
389 if ((order_found = rbtree_find(order_tree, &inc_order)) == NULL)
390 G_fatal_error(_("flat cell escaped for adjustment"));
391
392 uphill_order = order_found->uphill;
393 downhill_order = order_found->downhill;
394
395 /* debug */
396 if (uphill_order == -1)
397 G_fatal_error(_("adjustment: uphill order not set"));
398 if (!this_in_list && downhill_order == -1)
399 G_fatal_error(_("adjustment: downhill order not set"));
400
401 /* increment */
402 if (this_in_list) {
403 downhill_order = max_downhill_order;
404 uphill_order = 0;
405 }
406 alt_new[index_doer] +=
407 (uphill_order +
408 (double)(max_downhill_order - downhill_order) / 2.0 +
409 0.5) / 2.0 + 0.5;
410
411 /* check all neighbours, breadth first search */
412 for (ct_dir = 0; ct_dir < sides; ct_dir++) {
413 /* get r, c (upr, upc) for this neighbour */
414 upr = r + nextdr[ct_dir];
415 upc = c + nextdc[ct_dir];
416 /* check if r, c are within region */
417 if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
418 index_up = SEG_INDEX(alt_seg, upr, upc);
419 is_in_list = FLAG_GET(in_list, upr, upc);
420 is_worked = FLAG_GET(worked, upr, upc);
421 ele_nbr = alt_org[index_up];
422
423 if (ele_nbr == ele && !is_worked) {
424
425 inc_order.index = index_up;
426 if ((nbr_order_found =
427 rbtree_find(order_tree, &inc_order)) == NULL)
428 G_fatal_error(_("flat cell escaped in adjustment"));
429
430 /* not yet added to queue */
431 if (nbr_order_found->flag == 0) {
432 if (is_in_list)
433 G_warning
434 ("adjustment: in_list cell should be in queue");
435 /* add to up queue */
436 pq_add(index_up, up_pq);
437 nbr_order_found->flag = 1;
438 }
439 }
440 }
441 }
442 }
443
444 /* clean up */
445 pq_destroy(up_pq);
446 pq_destroy(down_pq);
447 rbtree_destroy(order_tree);
448
449 return 1;
450 }
451