2
/****************************************************************************
6
* AUTHOR(S): Laura Toma, Bowdoin College - ltoma@bowdoin.edu
7
* Yi Zhuang - yzhuang@bowdoin.edu
9
* Ported to GRASS by William Richard -
10
* wkrichar@bowdoin.edu or willster3021@gmail.com
11
* Markus Metz: surface interpolation
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.
31
* COPYRIGHT: (C) 2008 by the GRASS Development Team
33
* This program is free software under the GNU General Public License
34
* (>=v2). Read the file COPYING that comes with GRASS for details.
36
*****************************************************************************/
47
#include <grass/config.h>
48
#include <grass/gis.h>
49
#include <grass/glocale.h>
57
/* ------------------------------------------------------------ */
59
void copy_header(GridHeader * a, GridHeader b)
64
a->xllcorner = b.xllcorner;
65
a->yllcorner = b.yllcorner;
68
a->nodata_value = b.nodata_value;
73
/* ------------------------------------------------------------ */
74
/*returns 1 if value is Nodata, 0 if it is not */
75
int is_nodata(GridHeader * hd, float value)
79
return Rast_is_null_value(&value, G_SURFACE_TYPE);
83
/* ------------------------------------------------------------ */
84
/*returns 1 if value is Nodata, 0 if it is not */
85
int is_nodata(Grid * grid, float value)
88
return is_nodata(grid->hd, value);
93
/* ------------------------------------------------------------ */
94
/* create an empty grid and return it. The header and the data are set
96
Grid *create_empty_grid()
99
Grid *ptr_grid = (Grid *) G_malloc(sizeof(Grid));
103
/*initialize structure */
105
ptr_grid->grid_data = NULL;
108
printf("**DEBUG: createEmptyGrid \n");
118
/* ------------------------------------------------------------ */
119
/* allocate memroy for grid_data; grid must have a header that gives
121
void alloc_grid_data(Grid * pgrid)
126
pgrid->grid_data = (float **)G_malloc(pgrid->hd->nrows * sizeof(float *));
128
assert(pgrid->grid_data);
132
for (i = 0; i < pgrid->hd->nrows; i++) {
133
pgrid->grid_data[i] =
134
(float *)G_malloc(pgrid->hd->ncols * sizeof(float));
136
assert(pgrid->grid_data[i]);
140
printf("**DEBUG: allocGridData\n");
148
/* ------------------------------------------------------------ */
149
/*destroy the structure and reclaim all memory allocated */
150
void destroy_grid(Grid * grid)
154
/*free grid data if its allocated */
155
if (grid->grid_data) {
158
for (i = 0; i < grid->hd->nrows; i++) {
159
if (!grid->grid_data[i])
160
G_free((float *)grid->grid_data[i]);
163
G_free((float **)grid->grid_data);
170
printf("**DEBUG: grid destroyed.\n");