← Back to branch summary

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

« back to all changes in this revision

Viewing changes to raster/r.viewshed/visibility.cpp

  • 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.viewshed/visibility.cpp
 
1
 
 
2
/****************************************************************************
 
3
 *
 
4
 * MODULE:       r.viewshed
 
5
 *
 
6
 * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
 
7
 *               Yi Zhuang - yzhuang@bowdoin.edu
 
8
 
 
9
 *               Ported to GRASS by William Richard -
 
10
 *               wkrichar@bowdoin.edu or willster3021@gmail.com
 
11
 *               Markus Metz: surface interpolation
 
12
 *
 
13
 * Date:         april 2011 
 
14
 * 
 
15
 * PURPOSE: To calculate the viewshed (the visible cells in the
 
16
 * raster) for the given viewpoint (observer) location.  The
 
17
 * visibility model is the following: Two points in the raster are
 
18
 * considered visible to each other if the cells where they belong are
 
19
 * visible to each other.  Two cells are visible to each other if the
 
20
 * line-of-sight that connects their centers does not intersect the
 
21
 * terrain. The terrain is NOT viewed as a tesselation of flat cells, 
 
22
 * i.e. if the line-of-sight does not pass through the cell center, 
 
23
 * elevation is determined using bilinear interpolation.
 
24
 * The viewshed algorithm is efficient both in
 
25
 * terms of CPU operations and I/O operations. It has worst-case
 
26
 * complexity O(n lg n) in the RAM model and O(sort(n)) in the
 
27
 * I/O-model.  For the algorithm and all the other details see the
 
28
 * paper: "Computing Visibility on * Terrains in External Memory" by
 
29
 * Herman Haverkort, Laura Toma and Yi Zhuang.
 
30
 *
 
31
 * COPYRIGHT: (C) 2008 by the GRASS Development Team
 
32
 *
 
33
 * This program is free software under the GNU General Public License
 
34
 * (>=v2). Read the file COPYING that comes with GRASS for details.
 
35
 *
 
36
 *****************************************************************************/
 
37
#include <stdio.h>
 
38
 
 
39
extern "C"
 
40
{
 
41
#include <grass/gis.h>
 
42
#include <grass/glocale.h>
 
43
}
 
44
 
 
45
#include "grid.h"
 
46
#include "visibility.h"
 
47
#include "grass.h"
 
48
 
 
49
 
 
50
 
 
51
/* ------------------------------------------------------------ */
 
52
/* viewpoint functions */
 
53
void print_viewpoint(Viewpoint vp)
 
54
{
 
55
    G_debug(3, "vp=(%d, %d, %.1f) ", vp.row, vp.col, vp.elev);
 
56
    return;
 
57
}
 
58
 
 
59
/* ------------------------------------------------------------ */
 
60
void set_viewpoint_coord(Viewpoint * vp, dimensionType row, dimensionType col)
 
61
{
 
62
    assert(vp);
 
63
    vp->row = row;
 
64
    vp->col = col;
 
65
    return;
 
66
}
 
67
 
 
68
/* ------------------------------------------------------------ */
 
69
void set_viewpoint_elev(Viewpoint * vp, float elev)
 
70
{
 
71
    assert(vp);
 
72
    vp->elev = elev;
 
73
    return;
 
74
}
 
75
 
 
76
/* ------------------------------------------------------------ */
 
77
/*copy from b to a */
 
78
void copy_viewpoint(Viewpoint * a, Viewpoint b)
 
79
{
 
80
    assert(a);
 
81
    a->row = b.row;
 
82
    a->col = b.col;
 
83
    a->elev = b.elev;
 
84
    return;
 
85
}
 
86
 
 
87
 
 
88
/* ------------------------------------------------------------ */
 
89
/* MemoryVisibilityGrid functions */
 
90
 
 
91
/* create and return a grid of the sizes specified in the header */
 
92
MemoryVisibilityGrid *create_inmem_visibilitygrid(GridHeader hd, Viewpoint vp)
 
93
{
 
94
 
 
95
    MemoryVisibilityGrid *visgrid;
 
96
 
 
97
    visgrid = (MemoryVisibilityGrid *) G_malloc(sizeof(MemoryVisibilityGrid));
 
98
 
 
99
    assert(visgrid);
 
100
 
 
101
 
 
102
    /* create the grid  */
 
103
    visgrid->grid = create_empty_grid();
 
104
    assert(visgrid->grid);
 
105
 
 
106
    /* create the header */
 
107
    visgrid->grid->hd = (GridHeader *) G_malloc(sizeof(GridHeader));
 
108
 
 
109
    assert(visgrid->grid->hd);
 
110
 
 
111
    /* set the header */
 
112
    copy_header(visgrid->grid->hd, hd);
 
113
 
 
114
    /* allocate the  Grid data */
 
115
    alloc_grid_data(visgrid->grid);
 
116
 
 
117
    /*allocate viewpoint */
 
118
    visgrid->vp = (Viewpoint *) G_malloc(sizeof(Viewpoint));
 
119
 
 
120
    assert(visgrid->vp);
 
121
    copy_viewpoint(visgrid->vp, vp);
 
122
 
 
123
    return visgrid;
 
124
}
 
125
 
 
126
 
 
127
 
 
128
 
 
129
/* ------------------------------------------------------------ */
 
130
void free_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid)
 
131
{
 
132
 
 
133
    assert(visgrid);
 
134
 
 
135
    if (visgrid->grid) {
 
136
        destroy_grid(visgrid->grid);
 
137
    }
 
138
    if (visgrid->vp) {
 
139
        G_free(visgrid->vp);
 
140
    }
 
141
    G_free(visgrid);
 
142
 
 
143
    return;
 
144
}
 
145
 
 
146
 
 
147
 
 
148
/* ------------------------------------------------------------ */
 
149
/*set all values of visgrid's Grid to the given value */
 
150
void set_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid, float val)
 
151
{
 
152
 
 
153
    assert(visgrid && visgrid->grid && visgrid->grid->hd &&
 
154
           visgrid->grid->grid_data);
 
155
 
 
156
    dimensionType i, j;
 
157
 
 
158
    for (i = 0; i < visgrid->grid->hd->nrows; i++) {
 
159
        assert(visgrid->grid->grid_data[i]);
 
160
        for (j = 0; j < visgrid->grid->hd->ncols; j++) {
 
161
            visgrid->grid->grid_data[i][j] = val;
 
162
        }
 
163
    }
 
164
    return;
 
165
}
 
166
 
 
167
 
 
168
 
 
169
/* ------------------------------------------------------------ */
 
170
/*set the (i,j) value of visgrid's Grid to the given value */
 
171
void add_result_to_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid,
 
172
                                        dimensionType i, dimensionType j,
 
173
                                        float val)
 
174
{
 
175
 
 
176
    assert(visgrid && visgrid->grid && visgrid->grid->hd &&
 
177
           visgrid->grid->grid_data);
 
178
    assert(i < visgrid->grid->hd->nrows);
 
179
    assert(j < visgrid->grid->hd->ncols);
 
180
    assert(visgrid->grid->grid_data[i]);
 
181
 
 
182
    visgrid->grid->grid_data[i][j] = val;
 
183
 
 
184
    return;
 
185
}
 
186
 
 
187
 
 
188
 
 
189
/* ------------------------------------------------------------ */
 
190
/*  The following functions are used to convert the visibility results
 
191
   recorded during the viewshed computation into the output grid into
 
192
   tehe output required by the user.  
 
193
 
 
194
   x is assumed to be the visibility value computed for a cell during the
 
195
   viewshed computation. 
 
196
 
 
197
   The value computed during the viewshed is the following:
 
198
 
 
199
   x is NODATA if the cell is NODATA; 
 
200
 
 
201
 
 
202
   x is INVISIBLE if the cell is invisible;
 
203
 
 
204
   x is the vertical angle of the cell wrt the viewpoint if the cell is
 
205
   visible---the angle is a value in (0,180).
 
206
 */
 
207
int is_visible(float x)
 
208
{
 
209
    /* if GRASS is on, we cannot guarantee that NODATA is negative; so
 
210
       we need to check */
 
211
    int isnull = Rast_is_null_value(&x, G_SURFACE_TYPE);
 
212
 
 
213
    if (isnull)
 
214
        return 0;
 
215
    else
 
216
        return (x >= 0);
 
217
}
 
218
int is_invisible_not_nodata(float x)
 
219
{
 
220
 
 
221
    return ((int)x == (int)INVISIBLE);
 
222
}
 
223
 
 
224
int is_invisible_nodata(float x)
 
225
{
 
226
 
 
227
    return (!is_visible(x)) && (!is_invisible_not_nodata(x));
 
228
}
 
229
 
 
230
/* ------------------------------------------------------------ */
 
231
/* This function is called when the program runs in
 
232
   viewOptions.outputMode == OUTPUT_BOOL. */
 
233
float booleanVisibilityOutput(float x)
 
234
{
 
235
    /* NODATA and INVISIBLE are both negative values */
 
236
    if (is_visible(x))
 
237
        return BOOL_VISIBLE;
 
238
    else
 
239
        return BOOL_INVISIBLE;
 
240
}
 
241
 
 
242
/* ------------------------------------------------------------ */
 
243
/* This function is called when the program runs in
 
244
   viewOptions.outputMode == OUTPUT_ANGLE. In this case x represents
 
245
   the right value.  */
 
246
float angleVisibilityOutput(float x)
 
247
{
 
248
 
 
249
    return x;
 
250
}
 
251
 
 
252
 
 
253
/* ------------------------------------------------------------ */
 
254
/* visgrid is the structure that records the visibility information
 
255
   after the sweep is done.  Use it to write the visibility output
 
256
   grid and then distroy it.
 
257
 */
 
258
void save_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid,
 
259
                               ViewOptions viewOptions, Viewpoint vp)
 
260
{
 
261
 
 
262
    if (viewOptions.outputMode == OUTPUT_BOOL)
 
263
        save_grid_to_GRASS(visgrid->grid, viewOptions.outputfname, CELL_TYPE,
 
264
                           OUTPUT_BOOL);
 
265
    else if (viewOptions.outputMode == OUTPUT_ANGLE)
 
266
        save_grid_to_GRASS(visgrid->grid, viewOptions.outputfname, FCELL_TYPE,
 
267
                           OUTPUT_ANGLE);
 
268
    else
 
269
        /* elevation  output */
 
270
        save_vis_elev_to_GRASS(visgrid->grid, viewOptions.inputfname,
 
271
                               viewOptions.outputfname,
 
272
                               vp.elev + viewOptions.obsElev);
 
273
 
 
274
    free_inmem_visibilitygrid(visgrid);
 
275
 
 
276
    return;
 
277
}
 
278
 
 
279
 
 
280
 
 
281
/* ------------------------------------------------------------ */
 
282
/* IOVisibilityGrid functions */
 
283
/* ------------------------------------------------------------ */
 
284
 
 
285
/* ------------------------------------------------------------ */
 
286
/*create grid from given header and viewpoint */
 
287
IOVisibilityGrid *init_io_visibilitygrid(GridHeader hd, Viewpoint vp)
 
288
{
 
289
    IOVisibilityGrid *visgrid;
 
290
 
 
291
    visgrid = (IOVisibilityGrid *) G_malloc(sizeof(IOVisibilityGrid));
 
292
 
 
293
    assert(visgrid);
 
294
 
 
295
    /*header */
 
296
    visgrid->hd = (GridHeader *) G_malloc(sizeof(GridHeader));
 
297
 
 
298
    assert(visgrid->hd);
 
299
    copy_header(visgrid->hd, hd);
 
300
 
 
301
    /*viewpoint */
 
302
    visgrid->vp = (Viewpoint *) G_malloc(sizeof(Viewpoint));
 
303
 
 
304
    assert(visgrid->vp);
 
305
    copy_viewpoint(visgrid->vp, vp);
 
306
 
 
307
    /*stream */
 
308
    visgrid->visStr = new AMI_STREAM < VisCell > ();
 
309
    assert(visgrid->visStr);
 
310
 
 
311
    return visgrid;
 
312
}
 
313
 
 
314
 
 
315
 
 
316
/* ------------------------------------------------------------ */
 
317
/*free the grid */
 
318
void free_io_visibilitygrid(IOVisibilityGrid * grid)
 
319
{
 
320
    assert(grid);
 
321
 
 
322
    if (grid->hd)
 
323
        G_free(grid->hd);
 
324
    if (grid->vp)
 
325
        G_free(grid->vp);
 
326
    if (grid->visStr)
 
327
        delete grid->visStr;
 
328
 
 
329
    G_free(grid);
 
330
 
 
331
    return;
 
332
}
 
333
 
 
334
 
 
335
 
 
336
/* ------------------------------------------------------------ */
 
337
/*write cell to stream */
 
338
void add_result_to_io_visibilitygrid(IOVisibilityGrid * visgrid,
 
339
                                     VisCell * cell)
 
340
{
 
341
 
 
342
    assert(visgrid && cell);
 
343
 
 
344
    AMI_err ae;
 
345
 
 
346
    assert(visgrid->visStr);
 
347
    ae = visgrid->visStr->write_item(*cell);
 
348
    assert(ae == AMI_ERROR_NO_ERROR);
 
349
    return;
 
350
}
 
351
 
 
352
 
 
353
/* ------------------------------------------------------------ */
 
354
/*compare function, (i,j) grid order */
 
355
int IJCompare::compare(const VisCell & a, const VisCell & b)
 
356
{
 
357
    if (a.row > b.row)
 
358
        return 1;
 
359
    if (a.row < b.row)
 
360
        return -1;
 
361
 
 
362
    /*a.row==b.row */
 
363
    if (a.col > b.col)
 
364
        return 1;
 
365
    if (a.col < b.col)
 
366
        return -1;
 
367
    /*all equal */
 
368
    return 0;
 
369
}
 
370
 
 
371
 
 
372
/* ------------------------------------------------------------ */
 
373
/*sort stream in grid order */
 
374
void sort_io_visibilitygrid(IOVisibilityGrid * visGrid)
 
375
{
 
376
 
 
377
    assert(visGrid);
 
378
    assert(visGrid->visStr);
 
379
    if (visGrid->visStr->stream_len() == 0)
 
380
        return;
 
381
 
 
382
    AMI_STREAM < VisCell > *sortedStr;
 
383
    AMI_err ae;
 
384
    IJCompare cmpObj;
 
385
 
 
386
    ae = AMI_sort(visGrid->visStr, &sortedStr, &cmpObj, 1);
 
387
    assert(ae == AMI_ERROR_NO_ERROR);
 
388
    assert(sortedStr);
 
389
    sortedStr->seek(0);
 
390
 
 
391
    visGrid->visStr = sortedStr;
 
392
    return;
 
393
}
 
394
 
 
395
 
 
396
/* ------------------------------------------------------------ */
 
397
void
 
398
save_io_visibilitygrid(IOVisibilityGrid * visgrid,
 
399
                       ViewOptions viewOptions, Viewpoint vp)
 
400
{
 
401
 
 
402
    if (viewOptions.outputMode == OUTPUT_BOOL)
 
403
        save_io_visibilitygrid_to_GRASS(visgrid, viewOptions.outputfname,
 
404
                                        CELL_TYPE, booleanVisibilityOutput);
 
405
 
 
406
    else if (viewOptions.outputMode == OUTPUT_ANGLE)
 
407
        save_io_visibilitygrid_to_GRASS(visgrid, viewOptions.outputfname,
 
408
                                        FCELL_TYPE, angleVisibilityOutput);
 
409
    else
 
410
        /* elevation  output */
 
411
        save_io_vis_and_elev_to_GRASS(visgrid, viewOptions.inputfname,
 
412
                                      viewOptions.outputfname,
 
413
                                      vp.elev + viewOptions.obsElev);
 
414
 
 
415
    /*free visibiliyty grid */
 
416
    free_io_visibilitygrid(visgrid);
 
417
 
 
418
    return;
 
419
}