← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.watershed/ram/do_astar.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

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