[GRASS-SVN] r46199 - grass-addons/raster/r.viewshed
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri May 6 03:30:56 EDT 2011
Author: mmetz
Date: 2011-05-06 00:30:56 -0700 (Fri, 06 May 2011)
New Revision: 46199
Removed:
grass-addons/raster/r.viewshed/print_message.cc
grass-addons/raster/r.viewshed/print_message.h
Modified:
grass-addons/raster/r.viewshed/Makefile
grass-addons/raster/r.viewshed/distribute.cc
grass-addons/raster/r.viewshed/distribute.h
grass-addons/raster/r.viewshed/eventlist.cc
grass-addons/raster/r.viewshed/eventlist.h
grass-addons/raster/r.viewshed/grass.cc
grass-addons/raster/r.viewshed/grass.h
grass-addons/raster/r.viewshed/grid.cc
grass-addons/raster/r.viewshed/grid.h
grass-addons/raster/r.viewshed/main.cc
grass-addons/raster/r.viewshed/r.viewshed.html
grass-addons/raster/r.viewshed/rbbst.cc
grass-addons/raster/r.viewshed/rbbst.h
grass-addons/raster/r.viewshed/statusstructure.cc
grass-addons/raster/r.viewshed/statusstructure.h
grass-addons/raster/r.viewshed/viewshed.cc
grass-addons/raster/r.viewshed/viewshed.h
grass-addons/raster/r.viewshed/visibility.cc
grass-addons/raster/r.viewshed/visibility.h
Log:
fix for ticket #390
Modified: grass-addons/raster/r.viewshed/Makefile
===================================================================
--- grass-addons/raster/r.viewshed/Makefile 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/Makefile 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,10 +8,9 @@
include $(MODULE_TOPDIR)/include/Make/Module.make
SOURCES = main.cc distribute.cc eventlist.cc grid.cc grass.cc viewshed.cc \
- rbbst.cc statusstructure.cc visibility.cc \
- print_message.cc
+ rbbst.cc statusstructure.cc visibility.cc
HEADERS = distribute.h eventlist.h grid.h grass.h viewshed.h \
- rbbst.h statusstructure.h visibility.h print_message.h
+ rbbst.h statusstructure.h visibility.h
OBJARCH=OBJ.$(ARCH)
Modified: grass-addons/raster/r.viewshed/distribute.cc
===================================================================
--- grass-addons/raster/r.viewshed/distribute.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/distribute.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -42,7 +41,6 @@
#include <stdio.h>
#include <math.h>
-#ifdef __GRASS__
extern "C"
{
#include <grass/config.h>
@@ -53,14 +51,8 @@
/* include IOSTREAM header */
#include <grass/iostream/ami.h>
-#else
-/* if GRASS is not defined */
-#include <ami.h>
-#endif
-
-
-/*note: MAX_STREAM_OPEN defined in IOStrea/include/ami_stream.h, which
+/*note: MAX_STREAM_OPEN defined in IOStream/include/ami_stream.h, which
is included by ami.h */
#include "distribute.h"
@@ -68,10 +60,7 @@
#include "visibility.h"
#include "eventlist.h"
#include "statusstructure.h"
-#include "print_message.h"
-#ifdef __GRASS__
#include "grass.h"
-#endif
#define DISTRIBDEBUG if(0)
@@ -82,7 +71,7 @@
#define ANGLE_FACTOR 1
-#define EPSILON .00000001
+#define EPSILON GRASS_EPSILON
#define PRINT_DISTRIBUTE if(0)
@@ -97,15 +86,14 @@
{
assert(inputfname && hd && vp);
- print_message("Start distributed sweeping.\n");
- fflush(stdout);
+ G_message(_("Start distributed sweeping."));
/* ------------------------------ */
/*initialize the visibility grid */
IOVisibilityGrid *visgrid;
visgrid = init_io_visibilitygrid(*hd, *vp);
- printf("distribute_and_sweep: visgrid=%s\n", visgrid->visStr->name());
+ G_debug(1, "distribute_and_sweep: visgrid=%s", visgrid->visStr->name());
/* ------------------------------ */
@@ -114,17 +102,14 @@
rt_start(initEventTime);
AMI_STREAM < AEvent > *eventList;
-#ifdef __GRASS__
- eventList = grass_init_event_list(inputfname, vp, hd,
+
+ eventList = init_event_list(inputfname, vp, hd,
viewOptions, NULL, visgrid);
-#else
- eventList =
- init_event_list(inputfname, vp, hd, viewOptions, NULL, visgrid);
-#endif
+
assert(eventList);
eventList->seek(0);
rt_stop(initEventTime);
- printf("distribute_and_sweep: eventlist=%s\n", eventList->sprint());
+ G_debug(1, "distribute_and_sweep: eventlist=%s", eventList->sprint());
/* ------------------------------ */
@@ -132,22 +117,16 @@
Rtimer sortEventTime;
rt_start(sortEventTime);
- PRINT_DISTRIBUTE {
- print_message("sorting events by distance from viewpoint..");
- fflush(stdout);
- }
+ G_debug(1, "Sorting events by distance from viewpoint..");
sort_event_list_by_distance(&eventList, *vp);
- PRINT_DISTRIBUTE {
- print_message("..sorting done.\n");
- fflush(stdout);
- }
+ G_debug(1, "..sorting done.");
/* debugging */
/*sortCheck(eventList, *vp); */
eventList->seek(0); /*this does not seem to be ensured by sort */
rt_stop(sortEventTime);
- printf("distribute_and_sweep: eventlist=%s\n", eventList->sprint());
+ G_debug(1, "distribute_and_sweep: eventlist=%s", eventList->sprint());
@@ -161,16 +140,15 @@
/*distribute recursively the events and write results to visgrid.
invariant: distribute_sector deletes its eventlist */
nvis = distribute_sector(eventList, NULL, 0, ANGLE_FACTOR * 2 * M_PI,
- visgrid, vp, viewOptions);
+ visgrid, vp, hd, viewOptions);
rt_stop(sweepTime);
/* ------------------------------ */
/*cleanup */
- print_message("Distribution sweeping done.\n");
- fflush(stdout);
+ G_message(_("Distribution sweeping done."));
- printf("Total cells %ld, visible cells %ld (%.1f percent).\n",
+ G_verbose_message("Total cells %ld, visible cells %ld (%.1f percent).",
(long)visgrid->hd->nrows * visgrid->hd->ncols,
nvis,
(float)((float)nvis * 100 /
@@ -202,19 +180,19 @@
AMI_STREAM < AEvent > *enterBndEvents,
double start_angle, double end_angle,
IOVisibilityGrid * visgrid, Viewpoint * vp,
- ViewOptions viewOptions)
+ GridHeader *hd, ViewOptions viewOptions)
{
assert(eventList && visgrid && vp);
/*enterBndEvents may be NULL first time */
- PRINT_DISTRIBUTE printf("***DISTRIBUTE sector [%.4f, %.4f]***\n",
+ G_debug(2, "*** DISTRIBUTE sector [%.4f, %.4f] ***",
start_angle, end_angle);
- printf("initial_gradient: %lf\n", SMALLEST_GRADIENT);
- printf("eventlist: %s\n", eventList->sprint());
+ G_debug(2, "initial_gradient: %lf", SMALLEST_GRADIENT);
+ G_debug(2, "eventlist: %s", eventList->sprint());
if (enterBndEvents)
- printf("BndEvents: %s\n", enterBndEvents->sprint());
+ G_debug(2, "BndEvents: %s", enterBndEvents->sprint());
PRINT_DISTRIBUTE LOG_avail_memo();
unsigned long nvis = 0;
@@ -226,8 +204,8 @@
MM_manager.memory_available()) {
if (enterBndEvents) {
nvis += solve_in_memory(eventList, enterBndEvents,
- start_angle, end_angle, visgrid, vp,
- viewOptions);
+ start_angle, end_angle, visgrid, hd,
+ vp, viewOptions);
return nvis;
}
else {
@@ -240,7 +218,7 @@
}
/*else, must recurse */
- PRINT_DISTRIBUTE print_message("in EXTERNAL memory\n");
+ PRINT_DISTRIBUTE G_debug(2, "In EXTERNAL memory");
/*compute number of sectors */
int nsect = compute_n_sectors();
@@ -299,20 +277,19 @@
double ssize = (end_angle - start_angle) / nsect;
boundaryEvents++;
- PRINTWARNING {
- print_message("WARNING!event ");
- print_event(*e);
- print_message("CLOSE TO BOUNDARY\n");
- printf("angle=%f close to sector boundaries=[%f, %f]\n",
- e->angle, s * ssize, (s + 1) * ssize);
- }
+
+ G_debug(2, "WARNING! event ");
+ print_event(*e, 3);
+ G_debug(2, "close to boundary");
+ G_debug(2, "angle=%f close to sector boundaries=[%f, %f]",
+ e->angle, s * ssize, (s + 1) * ssize);
}
- DISTRIBDEBUG printf("event %7lu: ", (unsigned long)i);
- DISTRIBDEBUG print_event(*e);
- DISTRIBDEBUG printf("d=%8.1f, ",
+ G_debug(2, "event %7lu: ", (unsigned long)i);
+ print_event(*e, 2);
+ G_debug(2, "d=%8.1f, ",
get_square_distance_from_viewpoint(*e, *vp));
- DISTRIBDEBUG printf("s=%3d ", s);
+ G_debug(2, "s=%3d ", s);
assert(is_inside(s, nsect));
total[s]++;
@@ -329,8 +306,8 @@
exit_angle = calculate_exit_angle(e->row, e->col, vp);
exit_s =
get_event_sector(exit_angle, start_angle, end_angle, nsect);
- DISTRIBDEBUG
- printf(" ENTER (a=%.2f,s=%3d)---> EXIT (a=%.2f,s=%3d) ",
+
+ G_debug(2, " ENTER (a=%.2f,s=%3d)---> EXIT (a=%.2f,s=%3d) ",
e->angle, s, exit_angle, exit_s);
/*note: exit_s can be -1 (outside) */
if (exit_s == s) {
@@ -359,8 +336,8 @@
enter_angle = calculate_enter_angle(e->row, e->col, vp);
enter_s =
get_event_sector(enter_angle, start_angle, end_angle, nsect);
- DISTRIBDEBUG
- printf(" EXIT (a=%.2f,s=%3d)--->ENTER (a=%.2f,s=%3d) ",
+
+ G_debug(2, " EXIT (a=%.2f,s=%3d)--->ENTER (a=%.2f,s=%3d) ",
e->angle, s, enter_angle, enter_s);
/*don't need to check spanned sectors because it is done on its
@@ -373,12 +350,12 @@
/*semi-short event
//the corresponding ENTER event must insert itself in sectorBnd[s] */
e->eventType = ENTERING_EVENT;
- BND_DEBUG {
- print_message("BND event ");
- print_event(*e);
- printf("in bndSector %d\n", s);
- fflush(stdout);
- }
+
+ G_debug(2, "BND event ");
+ print_event(*e, 2);
+ G_debug(2, "in bndSector %d", s);
+
+
insert_event_in_sector(e, s, §orBnd[s], high[s],
vp, bndInsert, bndDrop);
@@ -389,12 +366,11 @@
longEvents++;
/*the corresponding ENTER event must insert itself in sectorBnd[s] */
e->eventType = ENTERING_EVENT;
- BND_DEBUG {
- print_message("BND event ");
- print_event(*e);
- printf("in bndSector %d\n", s);
- fflush(stdout);
- }
+
+ G_debug(2, "BND event ");
+ print_event(*e, 2);
+ G_debug(2, "in bndSector %d", s);
+
insert_event_in_sector(e, s, §orBnd[s], high[s],
vp, bndInsert, bndDrop);
}
@@ -402,7 +378,7 @@
} /*switch event-type */
- DISTRIBDEBUG print_message("\n");
+ G_debug(2, "\n");
} /*for event i */
@@ -414,7 +390,7 @@
bndDrop);
/*sanity checks */
- PRINT_DISTRIBUTE printf("boundary events in distribution: %ld\n",
+ G_debug(2, "boundary events in distribution: %ld",
boundaryEvents);
print_sector_stats(nbEvents, nsect, high, total, insert, drop, sector,
sectorBnd, bndInsert,
@@ -434,21 +410,17 @@
delete enterBndEvents;
/*save stream names of new sectors */
-#ifdef __GRASS__
char **sectorName = (char **)G_malloc(nsect * sizeof(char *));
char **sectorBndName = (char **)G_malloc(nsect * sizeof(char *));
-#else
- char **sectorName = (char **)malloc(nsect * sizeof(char *));
- char **sectorBndName = (char **)malloc(nsect * sizeof(char *));
-#endif
+
assert(sectorName && sectorBndName);
for (int i = 0; i < nsect; i++) {
sector[i].name(§orName[i]);
- PRINT_DISTRIBUTE printf("saving stream %d: %s\t", i, sectorName[i]);
+ G_debug(2, "saving stream %d: %s\t", i, sectorName[i]);
sector[i].persist(PERSIST_PERSISTENT);
sectorBnd[i].name(§orBndName[i]);
- PRINT_DISTRIBUTE printf("saving BndStr %d: %s\n", i,
+ G_debug(2, "saving BndStr %d: %s", i,
sectorBndName[i]);
sectorBnd[i].persist(PERSIST_PERSISTENT);
}
@@ -465,20 +437,20 @@
for (int i = 0; i < nsect; i++) {
/*recover stream */
- PRINT_DISTRIBUTE printf("\nopening sector stream %s ", sectorName[i]);
+ G_debug(3, "opening sector stream %s ", sectorName[i]);
AMI_STREAM < AEvent > *str =
new AMI_STREAM < AEvent > (sectorName[i]);
assert(str);
- PRINT_DISTRIBUTE printf(" len=%lu\n",
+ G_debug(3, " len=%lu",
(unsigned long)str->stream_len());
/*recover boundary stream */
- PRINT_DISTRIBUTE printf("opening boundary sector stream %s ",
+ G_debug(3, "opening boundary sector stream %s ",
sectorBndName[i]);
AMI_STREAM < AEvent > *bndStr =
new AMI_STREAM < AEvent > (sectorBndName[i]);
assert(str);
- PRINT_DISTRIBUTE printf(" len=%lu\n",
+ G_debug(3, " len=%lu",
(unsigned long)bndStr->stream_len());
@@ -488,19 +460,14 @@
start_angle + (i +
1) * ((end_angle -
start_angle) / nsect),
- visgrid, vp, viewOptions);
+ visgrid, vp, hd, viewOptions);
}
/*cleanup */
-#ifdef __GRASS__
G_free(sectorName);
G_free(sectorBndName);
-#else
- free(sectorName);
- free(sectorBndName);
-#endif
- PRINT_DISTRIBUTE printf("Distribute sector [ %.4f, %.4f] done.\n",
+ G_debug(2, "Distribute sector [ %.4f, %.4f] done.",
start_angle, end_angle);
return nvis;
@@ -524,7 +491,7 @@
double *high, long *insert, long *drop)
{
- PRINT_DISTRIBUTE printf("Distribute boundary of sector [ %.4f, %.4f] ",
+ G_debug(3, "Distribute boundary of sector [ %.4f, %.4f] ",
start_angle, end_angle);
assert(bndEvents && sectorBnd && vp && high && insert && drop);
AEvent *e;
@@ -560,8 +527,7 @@
} /*for i */
- PRINT_DISTRIBUTE
- printf("Distribute boundary of sector [ %.4f, %.4f] done.\n",
+ G_debug(3, "Distribute boundary of sector [ %.4f, %.4f] done.",
start_angle, end_angle);
return;
@@ -571,7 +537,7 @@
//***********************************************************************
-/* Solves a segment it inemory. it is called by distribute() when
+/* Solves a segment in memory. it is called by distribute() when
sector fits in memory. eventList is the list of events in
increasing order of distance from the viewpoint; enterBndEvents is
the list of ENTER events that are outside the sector, whose
@@ -581,19 +547,19 @@
solved by switching to radial sweep. */
unsigned long solve_in_memory(AMI_STREAM < AEvent > *eventList,
AMI_STREAM < AEvent > *enterBndEvents,
- double start_angle,
- double end_angle, IOVisibilityGrid * visgrid,
+ double start_angle, double end_angle,
+ IOVisibilityGrid * visgrid, GridHeader *hd,
Viewpoint * vp, ViewOptions viewOptions)
{
assert(eventList && visgrid && vp);
- PRINT_DISTRIBUTE print_message("solve INTERNAL memory\n");
+ G_debug(2, "solve INTERNAL memory");
unsigned long nvis = 0; /*number of visible cells */
- printf("solve_in_memory: eventlist: %s\n", eventList->sprint());
+ G_debug(2, "solve_in_memory: eventlist: %s", eventList->sprint());
if (enterBndEvents)
- printf("BndEvents: %s\n", enterBndEvents->sprint());
+ G_debug(2, "BndEvents: %s", enterBndEvents->sprint());
if (eventList->stream_len() == 0) {
delete eventList;
@@ -646,31 +612,49 @@
assert(ae == AMI_ERROR_END_OF_STREAM);
*/
if (enterBndEvents) {
+ double ax, ay;
+
enterBndEvents->seek(0);
inevents = enterBndEvents->stream_len();
for (off_t i = 0; i < inevents; i++) {
ae = enterBndEvents->read_item(&e);
assert(ae == AMI_ERROR_NO_ERROR);
- DEBUGINIT {
- print_message("INMEM init: initializing boundary ");
- print_event(*e);
- print_message("\n");
- }
+
+ G_debug(3, "INMEM init: initializing boundary ");
+ print_event(*e, 3);
+ G_debug(3, "\n");
+
/*this must span the first boundary of this sector; insert it
//in status structure */
sn.col = e->col;
sn.row = e->row;
- sn.elev = e->elev;
- calculate_dist_n_gradient(&sn, vp);
+
+ e->eventType = ENTERING_EVENT;
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_event_gradient(&sn, 0, ay, ax, e->elev[0], vp, *hd);
+
+ e->eventType = CENTER_EVENT;
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_dist_n_gradient(&sn, e->elev[1], vp, *hd);
+
+ e->eventType = EXITING_EVENT;
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_event_gradient(&sn, 2, ay, ax, e->elev[2], vp, *hd);
+
+ if (sn.angle[0] > sn.angle[1])
+ sn.angle[0] -= 2 * M_PI;
+
+ //calculate_dist_n_gradient(&sn, vp);
insert_into_status_struct(sn, status_struct);
}
}
- PRINT_DISTRIBUTE {
- printf("initialized active structure with %d events\n", inevents);
- fflush(stdout);
- }
+ G_debug(2, "initialized active structure with %d events", inevents);
+
/*sweep the event list */
VisCell viscell;
off_t nbEvents = eventList->stream_len();
@@ -678,42 +662,73 @@
/*printf("nbEvents = %ld\n", (long) nbEvents); */
eventList->seek(0);
for (off_t i = 0; i < nbEvents; i++) {
-
/*get out one event at a time and process it according to its type */
ae = eventList->read_item(&e);
assert(ae == AMI_ERROR_NO_ERROR);
- SOLVEINMEMDEBUG {
- print_message("INMEM sweep: next event: ");
- print_event(*e);
- }
+ G_debug(3, "INMEM sweep: next event: ");
+ print_event(*e, 3);
+
sn.col = e->col;
sn.row = e->row;
- sn.elev = e->elev;
+ //sn.elev = e->elev;
/*calculate Distance to VP and Gradient */
- calculate_dist_n_gradient(&sn, vp);
+ calculate_dist_n_gradient(&sn, e->elev[1] + vp->target_offset, vp, *hd);
switch (e->eventType) {
case ENTERING_EVENT:
+ double ax, ay;
+
/*insert node into structure */
- SOLVEINMEMDEBUG {
- print_message("..ENTER-EVENT: insert\n");
- }
+ G_debug(3, "..ENTER-EVENT: insert");
+
/*don't insert if its close to the boundary---the segment was
//already inserted in initialization above
//if (!is_almost_on_boundary(e->angle, start_angle))
//insertIntoStatusStruct(sn,status_struct); */
+ sn.angle[0] = calculate_enter_angle(sn.row, sn.col, vp);
+ sn.angle[1] = calculate_angle(sn.col, sn.row, vp->col, vp->row);
+ sn.angle[2] = calculate_exit_angle(sn.row, sn.col, vp);
+
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ //sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row);
+ sn.angle[0] = e->angle;
+ calculate_event_gradient(&sn, 0, ay, ax, e->elev[0], vp, *hd);
+
+ e->eventType = CENTER_EVENT;
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_dist_n_gradient(&sn, e->elev[1], vp, *hd);
+
+ e->eventType = EXITING_EVENT;
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_event_gradient(&sn, 2, ay, ax, e->elev[2], vp, *hd);
+
+ e->eventType = ENTERING_EVENT;
+
+ if (e->angle < M_PI) {
+ if (sn.angle[0] > sn.angle[1])
+ sn.angle[0] -= 2 * M_PI;
+ }
+ else {
+ if (sn.angle[0] > sn.angle[1]) {
+ sn.angle[1] += 2 * M_PI;
+ sn.angle[2] += 2 * M_PI;
+ }
+ }
+
insert_into_status_struct(sn, status_struct);
break;
case EXITING_EVENT:
/*delete node out of status structure */
SOLVEINMEMDEBUG {
- print_message("..EXIT-EVENT: delete\n");
+ G_debug(3, "..EXIT-EVENT: delete");
/*find its corresponding enter event and its sector */
double enter_angle =
calculate_enter_angle(e->row, e->col, vp);
- printf(" EXIT (a=%f)--->ENTER (a=%f) ", e->angle,
+ G_debug(3, " EXIT (a=%f)--->ENTER (a=%f) ", e->angle,
enter_angle);
}
delete_from_status_struct(status_struct, sn.dist2vp);
@@ -721,7 +736,7 @@
case CENTER_EVENT:
SOLVEINMEMDEBUG {
- print_message("..QUERY-EVENT: query\n");
+ G_debug(3, "..QUERY-EVENT: query");
}
/*calculate visibility
@@ -733,15 +748,16 @@
double max;
max =
- find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
+ find_max_gradient_in_status_struct(status_struct, sn.dist2vp,
+ e->angle, sn.gradient[1]);
viscell.row = sn.row;
viscell.col = sn.col;
- if (max <= sn.gradient_offset) {
+ if (max <= sn.gradient[1]) {
/*the point is visible */
viscell.angle =
- get_vertical_angle(*vp, sn, viewOptions.doCurv);
+ get_vertical_angle(*vp, sn, e->elev[1] + vp->target_offset, viewOptions.doCurv);
assert(viscell.angle > 0);
/* viscell.vis = VISIBLE; */
add_result_to_io_visibilitygrid(visgrid, &viscell);
@@ -760,10 +776,9 @@
} /* for each event */
- PRINT_DISTRIBUTE print_message("in memory sweeping done.\n");
+ G_debug(2, "in memory sweeping done.");
- PRINT_DISTRIBUTE
- printf("Total cells %lu, visible cells %lu (%.1f percent).\n",
+ G_debug(2, "Total cells %lu, visible cells %lu (%.1f percent).",
(unsigned long)eventList->stream_len(), nvis,
(float)((float)nvis * 100 / (float)(eventList->stream_len())));
@@ -868,7 +883,7 @@
Viewpoint * vp, AEvent * e, double *high)
{
- DISTRIBDEBUG printf("LONG CELL: spans [%3d, %3d] ", start_s, end_s);
+ G_debug(4, "LONG CELL: spans [%3d, %3d] ", start_s, end_s);
double ctrgrad = calculate_center_gradient(e, vp);
/*ENTER event is outside */
@@ -954,49 +969,27 @@
}
assert(totalSector == nevents);
-#ifdef __GRASS__
PRINT_DISTRIBUTE {
- G_message("-----nsectors=%d\n", nsect);
+ G_debug(3, "-----nsectors=%d", nsect);
for (int i = 0; i < nsect; i++) {
- G_message("\ts=%3d ", i);
- G_message("[%.4f, %.4f] ",
+ G_debug(3, "\ts=%3d ", i);
+ G_debug(3, "[%.4f, %.4f] ",
get_sector_start(i, start_angle, end_angle, nsect),
get_sector_end(i, start_angle, end_angle, nsect));
- G_message("high = %9.1f, ", high[i]);
- G_message("total = %10ld, ", total[i]);
- G_message("inserted = %10ld, ", insert[i]);
- G_message("dropped = %10ld, ", drop[i]);
- G_message("BOUNDARY = %5ld", bndInsert[i]);
- G_message("\n");
+ G_debug(3, "high = %9.1f, ", high[i]);
+ G_debug(3, "total = %10ld, ", total[i]);
+ G_debug(3, "inserted = %10ld, ", insert[i]);
+ G_debug(3, "dropped = %10ld, ", drop[i]);
+ G_debug(3, "BOUNDARY = %5ld", bndInsert[i]);
+ G_debug(3, "\n");
}
}
- G_message("Distribute [%.4f, %.4f]: nsect=%d, ", start_angle, end_angle,
- nsect);
- G_message
- ("total events %lu, inserted %lu, dropped %lu, long events=%ld\n",
+ G_debug(3, "Distribute [%.4f, %.4f]: nsect=%d, ",
+ start_angle, end_angle, nsect);
+ G_debug(3,
+ "total events %lu, inserted %lu, dropped %lu, long events=%ld",
totalSector, totalInsert, totalDrop, longEvents);
-#else
- PRINT_DISTRIBUTE {
- printf("-----nsectors=%d\n", nsect);
- for (int i = 0; i < nsect; i++) {
- printf("\ts=%3d ", i);
- printf("[%.4f, %.4f] ",
- get_sector_start(i, start_angle, end_angle, nsect),
- get_sector_end(i, start_angle, end_angle, nsect));
- printf("high = %9.1f, ", high[i]);
- printf("total = %10ld, ", total[i]);
- printf("inserted = %10ld, ", insert[i]);
- printf("dropped = %10ld, ", drop[i]);
- printf("BOUNDARY = %5ld", bndInsert[i]);
- printf("\n");
- }
- }
- printf("Distribute [%.4f, %.4f]: nsect=%d, ", start_angle, end_angle,
- nsect);
- printf("total events %lu, inserted %lu, dropped %lu, long events=%ld\n",
- totalSector, totalInsert, totalDrop, longEvents);
- fflush(stdout);
-#endif
+
return;
}
@@ -1033,8 +1026,8 @@
if (2 * nsect > MAX_STREAMS_OPEN - 10)
nsect = (MAX_STREAMS_OPEN - 10) / 2;
}
- printf("nsectors set to %d\n", nsect);
- fflush(stdout);
+ G_debug(1, "nsectors set to %d", nsect);
+
return nsect;
}
@@ -1101,8 +1094,8 @@
/*note: if on boundary, PRECISION ISSUES?? should insert both sectors? */
DISTRIBDEBUG {
- print_event(*e);
- printf(" insert in sector %3d\n", s);
+ print_event(*e, 2);
+ G_debug(2, " insert in sector %3d", s);
}
AMI_err ae;
@@ -1127,8 +1120,8 @@
if (!is_center_gradient_occluded(e, high_s, vp)) {
insert[s]++;
DISTRIBDEBUG {
- print_event(*e);
- printf(" insert in sector %3d\n", s);
+ print_event(*e, 2);
+ G_debug(2, " insert in sector %3d", s);
}
AMI_err ae;
@@ -1171,7 +1164,7 @@
assert(e && vp);
double eg = calculate_center_gradient(e, vp);
- printf(" dropping grad=%.2f, high=%.2f\n", eg, high);
+ G_debug(3, " dropping grad=%.2f, high=%.2f", eg, high);
return;
}
Modified: grass-addons/raster/r.viewshed/distribute.h
===================================================================
--- grass-addons/raster/r.viewshed/distribute.h 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/distribute.h 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -74,7 +73,7 @@
AMI_STREAM < AEvent > *enterBndEvents,
double start_angle, double end_angle,
IOVisibilityGrid * visgrid, Viewpoint * vp,
- ViewOptions viewOptions);
+ GridHeader *hd, ViewOptions viewOptions);
/* bndEvents is a stream of events that cross into the sector's
(first) boundary; they must be distributed to the boundary streams
@@ -101,8 +100,8 @@
unsigned long solve_in_memory(AMI_STREAM < AEvent > *eventList,
AMI_STREAM < AEvent > *enterBndEvents,
double start_angle, double end_angle,
- IOVisibilityGrid * visgrid, Viewpoint * vp,
- ViewOptions viewOptions);
+ IOVisibilityGrid * visgrid, GridHeader *hd,
+ Viewpoint * vp, ViewOptions viewOptions);
/*returns 1 if enter angle is within epsilon from boundary angle */
Modified: grass-addons/raster/r.viewshed/eventlist.cc
===================================================================
--- grass-addons/raster/r.viewshed/eventlist.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/eventlist.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -72,9 +71,9 @@
sqdist = (e->row - vp->row) * (e->row - vp->row) +
(e->col - vp->col) * (e->col - vp->col);
- gradient = (e->elev - vp->elev) * (e->elev - vp->elev) / sqdist;
+ gradient = (e->elev[1] - vp->elev) * (e->elev[1] - vp->elev) / sqdist;
/*maintain sign */
- if (e->elev < vp->elev)
+ if (e->elev[1] < vp->elev)
gradient = -gradient;
return gradient;
}
@@ -121,7 +120,8 @@
double x, y;
e.eventType = EXITING_EVENT;
- e.angle = e.elev = 0;
+ e.angle = 0;
+ e.elev[0] = e.elev[1] = e.elev[2] = 0;
e.row = row;
e.col = col;
calculate_event_position(e, vp->row, vp->col, &y, &x);
@@ -138,7 +138,8 @@
double x, y;
e.eventType = ENTERING_EVENT;
- e.angle = e.elev = 0;
+ e.angle = 0;
+ e.elev[0] = e.elev[1] = e.elev[2] = 0;
e.row = row;
e.col = col;
calculate_event_position(e, vp->row, vp->col, &y, &x);
@@ -151,7 +152,7 @@
double viewpointX, double viewpointY)
{
double angle = atan(fabs(eventY - viewpointY) / fabs(eventX - viewpointX));
-
+
/*M_PI is defined in math.h to represent 3.14159... */
if (viewpointY == eventY && eventX > viewpointX) {
return 0; /*between 1st and 4th quadrant */
@@ -168,7 +169,7 @@
}
else if (eventX < viewpointX && eventY < viewpointY) {
/*second quadrant */
- return ((M_PI) - angle);
+ return (M_PI - angle);
}
else if (viewpointY == eventY && eventX < viewpointX) {
@@ -183,11 +184,11 @@
}
else if (viewpointX == eventX && viewpointY < eventY) {
/*between 3rd and 4th quadrant */
- return ((M_PI) + (M_PI) / 2);
+ return (M_PI * 3.0 / 2.0);
}
else if (eventX > viewpointX && eventY > viewpointY) {
/*4th quadrant */
- return ((M_PI) * 2.0 - angle);
+ return (M_PI * 2.0 - angle);
}
assert(eventX == viewpointX && eventY == viewpointY);
return 0;
@@ -338,19 +339,161 @@
}
assert(fabs(*x - e.col) < 1 && fabs(*y - e.row) < 1);
- /*if ((fabs(*x -e.col) >=1) || (fabs(*y -e.row) >=1)) {
- //printf("x-e.col=%f, y-e.row=%f ", fabs(*x -e.col), fabs(*y -e.row));
- //print_event(e);
- //printf("vp=(%d, %d), x=%.3f, y=%.3f", viewpointRow, viewpointCol, *x, *y);
- //assert(0);
- //exit(1);
- // } */
+ /*
+ if ((fabs(*x -e.col) >=1) || (fabs(*y -e.row) >=1)) {
+ G_warning("x-e.col=%f, y-e.row=%f ", fabs(*x -e.col), fabs(*y -e.row));
+ print_event(e, 0);
+ G_warning("vp=(%d, %d), x=%.3f, y=%.3f", viewpointRow, viewpointCol, *x, *y);
+ exit(1);
+ }
+ */
return;
}
+void
+calculate_event_row_col(AEvent e, dimensionType viewpointRow,
+ dimensionType viewpointCol, int *y, int *x)
+{
+ assert(x && y);
+ *x = 0;
+ *y = 0;
+ if (e.eventType == CENTER_EVENT) {
+ G_fatal_error("calculate_event_row_col() must not be called for CENTER events");
+ }
+
+ if (e.row < viewpointRow && e.col < viewpointCol) {
+ /*first quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row - 1;
+ *x = e.col + 1;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row + 1;
+ *x = e.col - 1;
+ }
+
+ }
+ else if (e.col == viewpointCol && e.row < viewpointRow) {
+ /*between the first and second quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row + 1;
+ *x = e.col + 1;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row + 1;
+ *x = e.col - 1;
+ }
+
+ }
+ else if (e.col > viewpointCol && e.row < viewpointRow) {
+ /*second quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row + 1;
+ *x = e.col + 1;
+ }
+ else { /*otherwise it is EXITING_EVENT */
+ *y = e.row - 1;
+ *x = e.col - 1;
+ }
+
+ }
+ else if (e.row == viewpointRow && e.col > viewpointCol) {
+ /*between the second and the fourth quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row + 1;
+ *x = e.col - 1;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row - 1;
+ *x = e.col - 1;
+ }
+
+ }
+ else if (e.col > viewpointCol && e.row > viewpointRow) {
+ /*fourth quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row + 1;
+ *x = e.col - 1;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row - 1;
+ *x = e.col + 1;
+ }
+
+ }
+ else if (e.col == viewpointCol && e.row > viewpointRow) {
+ /*between the third and fourth quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row - 1;
+ *x = e.col - 1;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row - 1;
+ *x = e.col + 1;
+ }
+
+ }
+ else if (e.col < viewpointCol && e.row > viewpointRow) {
+ /*third quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row - 1;
+ *x = e.col - 1;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row + 1;
+ *x = e.col + 1;
+ }
+
+ }
+ else if (e.row == viewpointRow && e.col < viewpointCol) {
+ /*between first and third quadrant */
+ if (e.eventType == ENTERING_EVENT) { /*if it is ENTERING_EVENT */
+ *y = e.row - 1;
+ *x = e.col + 1;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row + 1;
+ *x = e.col + 1;
+ }
+ }
+ else {
+ /*must be the viewpoint cell itself */
+ G_debug(1, "calculate_event_row_col() called for viewpoint cell itself");
+ assert(e.row == viewpointRow && e.col == viewpointCol);
+ *x = e.col;
+ *y = e.row;
+ }
+
+ /* assert(fabs(*x - e.col) <= 1 && fabs(*y - e.row) <= 1); */
+
+ if ((abs(*x - e.col) > 1) || (abs(*y - e.row) > 1)) {
+ G_warning("calculate_event_row_col() :");
+ G_warning("x-e.col=%d, y-e.row=%d", abs(*x - e.col), abs(*y - e.row));
+ print_event(e, 0);
+ G_warning("vp=(%d, %d), x=%d, y=%d", viewpointRow, viewpointCol, *x, *y);
+ exit(1);
+ }
+
+ return;
+}
+
/* ------------------------------------------------------------ */
-void print_event(AEvent a)
+void print_event(AEvent a, int debug_level)
{
char c = '0';
@@ -360,8 +503,13 @@
c = 'X';
if (a.eventType == CENTER_EVENT)
c = 'Q';
- printf("ev=[(%3d, %3d), e=%8.1f a=%4.2f t=%c] ",
- a.row, a.col, a.elev, a.angle, c);
+
+ if (debug_level < 1)
+ G_warning("ev=[(%3d, %3d), e=%8.1f a=%4.2f t=%c] ",
+ a.row, a.col, a.elev[1], a.angle, c);
+ else
+ G_debug(debug_level, "ev=[(%3d, %3d), e=%8.1f a=%4.2f t=%c] ",
+ a.row, a.col, a.elev[1], a.angle, c);
return;
}
@@ -398,8 +546,8 @@
(eventy - vp.row) * (eventy - vp.row);
/*don't take sqrt, it is expensive; suffices for comparison */
- print_event(a);
- printf(" pos= (%.3f. %.3f) sqdist=%.3f\n", eventx, eventy, dist);
+ print_event(a, 2);
+ G_debug(2, " pos= (%.3f. %.3f) sqdist=%.3f", eventx, eventy, dist);
return dist;
}
@@ -417,408 +565,18 @@
if ((int)maxDist == INFINITY_DISTANCE)
return 0;
-#ifdef _GRASS_
if (maxDist < G_distance(G_col_to_easting(vp.col + 0.5, &hd.window),
G_row_to_northing(vp.row + 0.5, &hd.window),
G_col_to_easting(col + 0.5, &hd.window),
G_row_to_northing(row + 0.5, &hd.window))) {
return 1;
}
- else {
- return 0;
- }
-#else
-
- double dif_x, dif_y, sqdist;
-
- dif_x = (vp.col - col);
- dif_y = (vp.row - row);
-
- /* expensive to take squareroots so use squares */
- sqdist = (dif_x * dif_x + dif_y * dif_y) * hd.cellsize * hd.cellsize;
-
- if (sqdist > maxDist * maxDist) {
- return 1;
- }
- else {
- return 0;
- }
-#endif
+ return 0;
}
-void testeventlist(AEvent * elist, size_t n)
-{
- printf("testing event list..%lu ", n);
- fflush(stdout);
- AEvent e = { 0, 0, 0 };
- for (size_t i = 0; i < n; i++)
- elist[i] = e;
- printf("ok ");
- fflush(stdout);
-}
-
-/*///////////////////////////////////////////////////////////
- ------------------------------------------------------------
-
- input: an array capable to hold the max number of events, an arcascii
- file, a grid header and a viewpoint; action: figure out all events
- in the input file, and write them to the event list. data is
- allocated and initialized with all the cells on the same row as the
- viewpoint. it returns the number of events.
- */
-size_t
-init_event_list_in_memory(AEvent * eventList, char *inputfname,
- Viewpoint * vp,
- GridHeader * hd, ViewOptions viewOptions,
- double **data, MemoryVisibilityGrid * visgrid)
-{
- printf("computing events..");
- fflush(stdout);
- assert(eventList && inputfname && hd && vp && visgrid);
-
-
- /* data is used to store all the cells on the same row as the
- viewpoint. */
- *data = (double *)malloc(hd->ncols * sizeof(double));
- assert(*data);
-
- /* open input raster */
- FILE *grid_fp;
-
- grid_fp = fopen(inputfname, "r");
- assert(grid_fp);
- /* we do this just to position the pointer after header for reading
- the data */
- read_header_from_arcascii_file(grid_fp);
-
-
- /* scan throught the arcascii file */
- size_t nevents = 0;
- dimensionType i, j;
- double ax, ay;
- AEvent e;
-
- e.angle = -1;
- for (i = 0; i < hd->nrows; i++) {
- for (j = 0; j < hd->ncols; j++) {
-
- e.row = i;
- e.col = j;
- fscanf(grid_fp, "%f", &(e.elev));
- //printf("(i=%3d, j=%3d): e=%.1f\n", i,j,e.elev); fflush(stdout);
-
- /*write the row of data going through the viewpoint */
- if (i == vp->row) {
- (*data)[j] = e.elev;
- if (j == vp->col) {
- set_viewpoint_elev(vp, e.elev + viewOptions.obsElev);
- /*what to do when viewpoint is NODATA ? */
- if (is_nodata(hd, e.elev)) {
- printf("WARNING: viewpoint is NODATA. ");
- printf("Will assume its elevation is %.f\n", e.elev);
- }
- }
- }
-
- /*don't insert the viewpoint events in the list */
- if (i == vp->row && j == vp->col)
- continue;
-
- /*don't insert the nodata cell events */
- if (is_nodata(hd, e.elev)) {
- /* record this cell as being NODATA; this is necessary so
- that we can distingush invisible events, from nodata
- events in the output */
- add_result_to_inmem_visibilitygrid(visgrid, i, j,
- hd->nodata_value);
- continue;
- }
-
- /* if point is outside maxDist, do NOT include it as an
- event */
- if (is_point_outside_max_dist
- (*vp, *hd, i, j, viewOptions.maxDist))
- continue;
-
- e.eventType = ENTERING_EVENT;
- calculate_event_position(e, vp->row, vp->col, &ay, &ax);
- e.angle = calculate_angle(ax, ay, vp->col, vp->row);
- eventList[nevents] = e;
- nevents++;
-
- e.eventType = CENTER_EVENT;
- calculate_event_position(e, vp->row, vp->col, &ay, &ax);
- e.angle = calculate_angle(ax, ay, vp->col, vp->row);
- eventList[nevents] = e;
- nevents++;
-
- e.eventType = EXITING_EVENT;
- calculate_event_position(e, vp->row, vp->col, &ay, &ax);
- e.angle = calculate_angle(ax, ay, vp->col, vp->row);
- eventList[nevents] = e;
- nevents++;
-
- } /* for j */
- } /* for i */
- fclose(grid_fp);
-
- printf("..done\n");
- fflush(stdout);
- printf("Event array size: %lu x %dB (%lu MB)\n",
- (unsigned long)nevents, (int)sizeof(AEvent),
- (unsigned long)(((long long)(nevents * sizeof(AEvent))) >> 20));
- fflush(stdout);
-
- return nevents;
-}
-
-
-
-
-
-
-
-/*///////////////////////////////////////////////////////////
- ------------------------------------------------------------
- input: an arcascii file, a grid header and a viewpoint; action:
- figure out all events in the input file, and write them to the
- stream. It assumes the file pointer is positioned rigth after the
- grid header so that this function can read all grid elements.
-
- if data is not NULL, it creates an array that stores all events on
- the same row as the viewpoint.
- */
-AMI_STREAM < AEvent > *init_event_list(char *inputfname, Viewpoint * vp,
- GridHeader * hd,
- ViewOptions viewOptions, double **data,
- IOVisibilityGrid * visgrid)
-{
- printf("computing events..");
- fflush(stdout);
- assert(inputfname && hd && vp && visgrid);
-
- /*create the event stream that will hold the events */
- AMI_STREAM < AEvent > *eventList = new AMI_STREAM < AEvent > ();
- assert(eventList);
-
- if (data != NULL) {
- /*data is used to store all the cells on the same row as the
- //viewpoint. */
- *data = (double *)malloc(hd->ncols * sizeof(double));
- assert(*data);
- }
-
- FILE *grid_fp = fopen(inputfname, "r");
-
- assert(grid_fp);
- /*we do this just to position the pointer after header for reading
- //the data
- //GridHeader* foo = */
- read_header_from_arcascii_file(grid_fp);
-
- /*scan throught the arcascii file */
- dimensionType i, j;
- double ax, ay;
- AEvent e;
-
- e.angle = -1;
- for (i = 0; i < hd->nrows; i++) {
- for (j = 0; j < hd->ncols; j++) {
-
- e.row = i;
- e.col = j;
- fscanf(grid_fp, "%f", &(e.elev));
-
- if (data != NULL) {
- /*write the row of data going through the viewpoint */
- if (i == vp->row) {
- (*data)[j] = e.elev;
- }
- }
-
- if (i == vp->row && j == vp->col) {
- set_viewpoint_elev(vp, e.elev + viewOptions.obsElev);
- /*what to do when viewpoint is NODATA */
- if (is_nodata(hd, e.elev)) {
- printf("WARNING: viewpoint is NODATA. ");
- printf("Will assume its elevation is %.f\n", e.elev);
- };
- }
-
- /*don't insert the viewpoint events in the list */
- if (i == vp->row && j == vp->col)
- continue;
-
-
- /*don't insert the nodata cell events */
- if (is_nodata(hd, e.elev)) {
- /* record this cell as being NODATA. ; this is necessary so
- that we can distingush invisible events, from nodata
- events in the output */
- VisCell visCell = { i, j, hd->nodata_value };
- add_result_to_io_visibilitygrid(visgrid, &visCell);
- continue;
- }
-
- /* if point is outside maxDist, do NOT include it as an
- event */
- if (is_point_outside_max_dist
- (*vp, *hd, i, j, viewOptions.maxDist))
- continue;
-
-
- e.eventType = ENTERING_EVENT;
- calculate_event_position(e, vp->row, vp->col, &ay, &ax);
- e.angle = calculate_angle(ax, ay, vp->col, vp->row);
- eventList->write_item(e);
-
- e.eventType = CENTER_EVENT;
- calculate_event_position(e, vp->row, vp->col, &ay, &ax);
- e.angle = calculate_angle(ax, ay, vp->col, vp->row);
- eventList->write_item(e);
-
- e.eventType = EXITING_EVENT;
- calculate_event_position(e, vp->row, vp->col, &ay, &ax);
- e.angle = calculate_angle(ax, ay, vp->col, vp->row);
- eventList->write_item(e);
- } /* for j */
- } /* for i */
-
- fclose(grid_fp);
-
- printf("..done\n");
- printf("nbEvents = %lu\n", (unsigned long)eventList->stream_len());
- printf("Event stream length: %lu x %dB (%lu MB)\n",
- (unsigned long)eventList->stream_len(), (int)sizeof(AEvent),
- (unsigned
- long)(((long long)(eventList->stream_len() *
- sizeof(AEvent))) >> 20));
- fflush(stdout);
-
- return eventList;
-}
-
-
-
-
-
-
-
-// /*///////////////////////////////////////////////////////////
-// ------------------------------------------------------------
-// input: an arcascii file, a grid header and a viewpoint; action:
-// figure out all events in the input file, and write them to the
-// stream. data is allocated and initialized with all the cells on the
-// same row as the viewpoint. It assumes the file pointer is positioned
-// rigth after the grid header so that this function can read all grid
-// elements.
-// */
-// AMI_STREAM < AEvent > *
-// init_event_list(FILE * grid_fp, Viewpoint * vp,
-// GridHeader * hd, ViewOptions viewOptions,
-// IOVisibilityGrid *visgrid) {
-
-// printf("computing events..");
-// fflush(stdout);
-// assert(grid_fp && hd && vp && visgrid);
-
-// /*create the event stream that will hold the events */
-// AMI_STREAM < AEvent > *eventList = new AMI_STREAM < AEvent > ();
-// assert(eventList);
-
-
-// /*we do this just to position the pointer after header for reading
-// //the data GridHeader* foo = */
-// read_header_from_arcascii_file(grid_fp);
-
-// /*scan throught the arcascii file */
-// dimensionType i, j;
-// double ax, ay;
-// AEvent e;
-
-// e.angle = -1;
-// for (i = 0; i < hd->nrows; i++) {
-// for (j = 0; j < hd->ncols; j++) {
-
-// e.row = i;
-// e.col = j;
-// fscanf(grid_fp, "%f", &(e.elev));
-
-// if (i == vp->row && j == vp->col) {
-// set_viewpoint_elev(vp, e.elev+viewOptions.obsElev);
-// /*what to do when viewpoint is NODATA */
-// if (is_nodata(hd, e.elev)) {
-// printf("WARNING: viewpoint is NODATA. ");
-// printf("Will assume its elevation is %.f\n", e.elev);
-// };
-// }
-
-// /*don't insert the viewpoint events in the list */
-// if (i == vp->row && j == vp->col)
-// continue;
-
-
-// /*don't insert the nodata cell events */
-// if (is_nodata(hd, e.elev)) {
-// /*printf("(%d, %d) dropping\n", i, j); */
-// /* record this cell as being NODATA; this is necessary so
-// that we can distingush invisible events, from nodata
-// events in the output */
-// VisCell viscell = {i,j, hd->nodata_value};
-// add_result_to_io_visibilitygrid(visgrid, &viscell);
-// continue;
-// }
-
-// /* if point is outside maxDist, do NOT include it as an
-// event */
-// if(is_point_outside_max_dist(*vp, *hd, i, j, viewOptions.maxDist))
-// continue;
-
-
-// e.eventType = ENTERING_EVENT;
-// calculate_event_position(e, vp->row, vp->col, &ay, &ax);
-// e.angle = calculate_angle(ax, ay, vp->col, vp->row);
-// eventList->write_item(e);
-// /*d = get_square_distance_from_viewpoint(e, *vp);
-// //printf("(%d, %d) insert ENTER (%.2f,%.2f) d=%.2f\n", i, j, ay, ax, d); */
-
-// e.eventType = CENTER_EVENT;
-// calculate_event_position(e, vp->row, vp->col, &ay, &ax);
-// e.angle = calculate_angle(ax, ay, vp->col, vp->row);
-// eventList->write_item(e);
-// /*d = get_square_distance_from_viewpoint(e, *vp);
-// //printf("(%d, %d) insert CENTER (%.2f,%.2f) d=%.2f\n", i,j, ay, ax, d); */
-
-// e.eventType = EXITING_EVENT;
-// calculate_event_position(e, vp->row, vp->col, &ay, &ax);
-// e.angle = calculate_angle(ax, ay, vp->col, vp->row);
-// eventList->write_item(e);
-// /*d = get_square_distance_from_viewpoint(e, *vp);
-// //printf("(%d, %d) insert EXIT (%.2f,%.2f) d=%.2f\n", i, j, ay, ax, d); */
-// }
-// }
-// printf("..done\n");
-// printf("nbEvents = %lu\n", (unsigned long)eventList->stream_len());
-// printf("Event stream length: %lu x %dB (%lu MB)\n",
-// (unsigned long)eventList->stream_len(), (int)sizeof(AEvent),
-// (unsigned
-// long)(((long long)(eventList->stream_len() *
-// sizeof(AEvent))) >> 20));
-// fflush(stdout);
-// return eventList;
-// }
-
-
-
-
-
-
-
-
-
/* ------------------------------------------------------------
//note: this is expensive because distance is not storedin the event
//and must be computed on the fly */
@@ -881,7 +639,7 @@
}
/* ------------------------------------------------------------ */
-/* a copy of the function above is needed by qsort, when teh
+/* a copy of the function above is needed by qsort, when the
computation runs in memory */
int radial_compare_events(const void *x, const void *y)
@@ -958,34 +716,3 @@
return;
}
-
-
-
-void sort_check(AMI_STREAM < AEvent > *eventList, Viewpoint vp)
-{
- printf("checking sort..");
- fflush(stdout);
- assert(eventList);
-
- size_t i, nbe = eventList->stream_len();
-
- eventList->seek(0);
- AEvent *crt, *next;
- double crtd, nextd;
- AMI_err ae;
-
- ae = eventList->read_item(&crt);
- assert(ae == AMI_ERROR_NO_ERROR);
- crtd = get_square_distance_from_viewpoint(*crt, vp);
-
- for (i = 0; i < nbe - 1; i++) {
- ae = eventList->read_item(&next);
- assert(ae == AMI_ERROR_NO_ERROR);
- nextd = get_square_distance_from_viewpoint(*next, vp);
- if (crtd > nextd)
- assert(0);
- crtd = nextd;
- }
- printf("..sort test passed\n");
- return;
-}
Modified: grass-addons/raster/r.viewshed/eventlist.h
===================================================================
--- grass-addons/raster/r.viewshed/eventlist.h 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/eventlist.h 2011-05-06 07:30:56 UTC (rev 46199)
@@ -5,11 +5,12 @@
*
* AUTHOR(S): Laura Toma, Bowdoin College - ltoma at bowdoin.edu
* Yi Zhuang - yzhuang at bowdoin.edu
-
+ *
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,19 +18,17 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
* paper: "Computing Visibility on * Terrains in External Memory" by
* Herman Haverkort, Laura Toma and Yi Zhuang.
*
- * COPYRIGHT: (C) 2008 by the GRASS Development Team
+ * COPYRIGHT: (C) 2008 - 2011 by the GRASS Development Team
*
* This program is free software under the GNU General Public License
* (>=v2). Read the file COPYING that comes with GRASS for details.
@@ -43,11 +42,7 @@
#include "grid.h"
#include "visibility.h"
-#ifdef __GRASS__
#include <grass/iostream/ami.h>
-#else
-#include <ami.h>
-#endif
#define ENTERING_EVENT 1
@@ -57,7 +52,11 @@
typedef struct event_
{
dimensionType row, col; //location of the center of cell
- float elev; //elevation here
+ // 3 elevation values:
+ // elev[0]: entering
+ // elev[1]: center
+ // elev[2]: exiting
+ surface_type elev[3]; //elevation here
double angle;
char eventType;
@@ -77,37 +76,6 @@
-/* ------------------------------------------------------------ */
-/* input: an array capable to hold the max number of events, an
- arcascii file, a grid header and a viewpoint; action: figure out
- all events in the input file, and write them to the event
- list. data is allocated and initialized with all the cells on the
- same row as the viewpoint. it returns the number of events.
- */
-size_t
-init_event_list_in_memory(AEvent * eventList, char *inputfname,
- Viewpoint * vp,
- GridHeader * hd, ViewOptions viewOptions,
- double **data, MemoryVisibilityGrid * visgrid);
-
-
-
-
-
-/* ------------------------------------------------------------ */
-/* input: an arcascii file, an eventlist stream, and a viewpoint.
- action: figure out all events in the input file, and write them to
- the stream. hd is initialized from file. if data is not NULL, data
- is allocated and initialized with all the cells on the same row as
- the viewpoint.
- */
-AMI_STREAM < AEvent > *init_event_list(char *inputfname, Viewpoint * vp,
- GridHeader * hd,
- ViewOptions viewOptions, double **data,
- IOVisibilityGrid * visgrid);
-
-
-
/*sort the event list by the angle around the viewpoint) */
void sort_event_list(AMI_STREAM < AEvent > **eventList);
@@ -123,7 +91,7 @@
{
public:int compare(const AEvent &, const AEvent &);
};
-void print_event(AEvent a);
+void print_event(AEvent a, int debug_level);
/*computes the distance from the event to the viewpoint. Note: all 3
@@ -135,7 +103,6 @@
/*sort the event list in distance order */
void sort_event_list_by_distance(AMI_STREAM < AEvent > **eventList,
Viewpoint vp);
-void sort_check(AMI_STREAM < AEvent > *eventList, Viewpoint vp);
/*return the angle from this event wrt viewpoint; the type of the
@@ -164,6 +131,13 @@
void calculate_event_position(AEvent e, dimensionType viewpointRow,
dimensionType viewpointCol, double *y,
double *x);
+ /* calculate the neighbouring row, col of the given event, and store them in x
+ //and y. */
+void
+calculate_event_row_col(AEvent e, dimensionType viewpointRow,
+ dimensionType viewpointCol, int *y, int *x);
+
+
double calculate_angle(double eventX, double eventY, double viewpointX,
double viewpointY);
Modified: grass-addons/raster/r.viewshed/grass.cc
===================================================================
--- grass-addons/raster/r.viewshed/grass.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/grass.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -5,11 +5,12 @@
*
* AUTHOR(S): Laura Toma, Bowdoin College - ltoma at bowdoin.edu
* Yi Zhuang - yzhuang at bowdoin.edu
-
+ *
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -40,8 +39,6 @@
#include <stdlib.h>
#include <stdio.h>
-#ifdef __GRASS__
-
extern "C"
{
#include <grass/config.h>
@@ -61,8 +58,8 @@
curvature of the earth; otherwise return the passed height
unchanged.
*/
-float adjust_for_curvature(Viewpoint vp, dimensionType row,
- dimensionType col, float h,
+surface_type adjust_for_curvature(Viewpoint vp, double row,
+ double col, surface_type h,
ViewOptions viewOptions, GridHeader *hd)
{
@@ -72,7 +69,6 @@
assert(viewOptions.ellps_a != 0);
/* distance must be in meters because ellps_a is in meters */
-
double dist = G_distance(G_col_to_easting(vp.col + 0.5, &(hd->window)),
G_row_to_northing(vp.row + 0.5, &(hd->window)),
G_col_to_easting(col + 0.5, &(hd->window)),
@@ -86,9 +82,8 @@
/* ************************************************************ */
-/*return a GridHeader that has all the relevant data filled in from
- GRASS */
-GridHeader *read_header_from_GRASS(char *rastName, Cell_head * region)
+/*return a GridHeader that has all the relevant data filled in */
+GridHeader *read_header(char *rastName, Cell_head * region)
{
assert(rastName);
@@ -123,17 +118,46 @@
"this case this may result in innacuracies."));
// exit(EXIT_FAILURE);
}
- hd->cellsize = (float)region->ew_res;
+ hd->ew_res = region->ew_res;
+ hd->ns_res = region->ns_res;
//store the null value of the map
- G_set_f_null_value(&(hd->nodata_value), 1);
+ G_set_null_value(&(hd->nodata_value), 1, G_SURFACE_TYPE);
G_message("Nodata value set to %f", hd->nodata_value);
+
+
+
return hd;
}
+/* calculate ENTER and EXIT event elevation (bilinear interpolation) */
+surface_type calculate_event_elevation(AEvent e, int nrows, int ncols,
+ dimensionType vprow, dimensionType vpcol,
+ G_SURFACE_T **inrast, RASTER_MAP_TYPE data_type)
+{
+ int row1, col1;
+ surface_type event_elev;
+ G_SURFACE_T elev1, elev2, elev3, elev4;
+
+ calculate_event_row_col(e, vprow, vpcol, &row1, &col1);
+ if (row1 >= 0 && row1 < nrows && col1 >= 0 && col1 < ncols) {
+ elev1 = inrast[row1 - e.row + 1][col1];
+ elev2 = inrast[row1 - e.row + 1][e.col];
+ elev3 = inrast[1][col1];
+ elev4 = inrast[1][e.col];
+ if (G_is_null_value(&elev1, data_type) || G_is_null_value(&elev2, data_type) ||
+ G_is_null_value(&elev3, data_type) || G_is_null_value(&elev4, data_type))
+ event_elev = inrast[1][e.col];
+ else {
+ event_elev = (elev1 + elev2 + elev3 + elev4) / 4.;
+ }
+ }
+ else
+ event_elev = inrast[1][e.col];
+ return event_elev;
+}
-
/* ************************************************************ */
/* input: an array capable to hold the max number of events, a raster
name, a viewpoint and the viewOptions; action: figure out all events
@@ -143,9 +167,9 @@
AEvent* with all the events for the map. Used when solving in
memory, so the AEvent* should fit in memory. */
size_t
-grass_init_event_list_in_memory(AEvent * eventList, char *rastName,
+init_event_list_in_memory(AEvent * eventList, char *rastName,
Viewpoint * vp, GridHeader * hd,
- ViewOptions viewOptions, double **data,
+ ViewOptions viewOptions, surface_type ***data,
MemoryVisibilityGrid * visgrid)
{
@@ -155,8 +179,12 @@
/*alloc data ; data is used to store all the cells on the same row
as the viewpoint. */
- *data = (double *)malloc(G_window_cols() * sizeof(double));
+ *data = (surface_type **)G_malloc(3 * sizeof(surface_type *));
assert(*data);
+ (*data)[0] = (surface_type *)G_malloc(3 * G_window_cols() * sizeof(surface_type));
+ assert((*data)[0]);
+ (*data)[1] = (*data)[0] + G_window_cols();
+ (*data)[2] = (*data)[1] + G_window_cols();
/*get the mapset name */
char *mapset;
@@ -174,13 +202,26 @@
/*get the data_type */
RASTER_MAP_TYPE data_type;
- data_type = G_raster_map_type(rastName, mapset);
+ /* data_type = G_raster_map_type(rastName, mapset); */
+ data_type = G_SURFACE_TYPE;
- /*buffer to hold a row */
- void *inrast;
+ /*buffer to hold 3 rows */
+ G_SURFACE_T **inrast;
+ int nrows = G_window_rows();
+ int ncols = G_window_rows();
- inrast = G_allocate_raster_buf(data_type);
+ inrast = (G_SURFACE_T **)G_malloc(3 * sizeof(G_SURFACE_T *));
assert(inrast);
+ inrast[0] = (G_SURFACE_T *)G_allocate_raster_buf(data_type);
+ assert(inrast[0]);
+ inrast[1] = (G_SURFACE_T *)G_allocate_raster_buf(data_type);
+ assert(inrast[1]);
+ inrast[2] = (G_SURFACE_T *)G_allocate_raster_buf(data_type);
+ assert(inrast[2]);
+
+ G_set_null_value(inrast[0], ncols, data_type);
+ G_set_null_value(inrast[1], ncols, data_type);
+ G_set_null_value(inrast[2], ncols, data_type);
/*DCELL to check for loss of prescion- haven't gotten that to work
yet though */
@@ -192,15 +233,29 @@
/*scan through the raster data */
dimensionType i, j, k;
+ int row1, col1;
double ax, ay;
+ G_SURFACE_T *elev1, *elev2, *elev3, *elev4;
AEvent e;
- int nrows = G_window_rows();
+
+ /* read first row */
+ G_get_raster_row(infd, inrast[2], 0, data_type);
e.angle = -1;
for (i = 0; i < nrows; i++) {
/*read in the raster row */
- int rasterRowResult = G_get_raster_row(infd, inrast, i, data_type);
+ int rasterRowResult = 1;
+
+ G_SURFACE_T *tmprast = inrast[0];
+ inrast[0] = inrast[1];
+ inrast[1] = inrast[2];
+ inrast[2] = tmprast;
+ if (i < nrows - 1)
+ rasterRowResult = G_get_raster_row(infd, inrast[2], i + 1, data_type);
+ else
+ G_set_null_value(inrast[2], ncols, data_type);
+
if (rasterRowResult <= 0)
G_fatal_error(_("Coord not read from row %d of <%s>"), i,
rastName);
@@ -212,33 +267,23 @@
e.row = i;
e.col = j;
- /*read the elevation value into the event, depending on data_type */
- switch (data_type) {
- case CELL_TYPE:
- isnull = G_is_c_null_value(&(((CELL *) inrast)[j]));
- e.elev = (float)(((CELL *) inrast)[j]);
- break;
- case FCELL_TYPE:
- isnull = G_is_f_null_value(&(((FCELL *) inrast)[j]));
- e.elev = (float)(((FCELL *) inrast)[j]);
- break;
- case DCELL_TYPE:
- isnull = G_is_d_null_value(&(((DCELL *) inrast)[j]));
- e.elev = (float)(((DCELL *) inrast)[j]);
- break;
- }
+ /*read the elevation value into the event */
+ isnull = G_is_null_value(&(inrast[1][j]), data_type);
+ e.elev[1] = inrast[1][j];
/* adjust for curvature */
- e.elev = adjust_for_curvature(*vp, i, j, e.elev, viewOptions, hd);
+ e.elev[1] = adjust_for_curvature(*vp, i, j, e.elev[1], viewOptions, hd);
/*write it into the row of data going through the viewpoint */
if (i == vp->row) {
- (*data)[j] = e.elev;
+ (*data)[0][j] = e.elev[1];
+ (*data)[1][j] = e.elev[1];
+ (*data)[2][j] = e.elev[1];
}
/* set the viewpoint, and don't insert it into eventlist */
if (i == vp->row && j == vp->col) {
- set_viewpoint_elev(vp, e.elev + viewOptions.obsElev);
+ set_viewpoint_elev(vp, e.elev[1] + viewOptions.obsElev);
if (viewOptions.tgtElev > 0)
vp->target_offset = viewOptions.tgtElev;
else
@@ -249,6 +294,9 @@
G_message(_("Will assume its elevation is = %f"),
vp->elev);
}
+
+ add_result_to_inmem_visibilitygrid(visgrid, i, j,
+ 180);
continue;
}
@@ -271,6 +319,34 @@
/* if it got here it is not the viewpoint, not NODATA, and
within max distance from viewpoint; generate its 3 events
and insert them */
+
+ /* get ENTER elevation */
+ e.eventType = ENTERING_EVENT;
+ e.elev[0] = calculate_event_elevation(e, nrows, ncols,
+ vp->row, vp->col, inrast, data_type);
+ /* adjust for curvature */
+ if (viewOptions.doCurv) {
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ e.elev[0] = adjust_for_curvature(*vp, ay, ax, e.elev[0], viewOptions, hd);
+ }
+
+ /* get EXIT elevation */
+ e.eventType = EXITING_EVENT;
+ e.elev[2] = calculate_event_elevation(e, nrows, ncols,
+ vp->row, vp->col, inrast, data_type);
+ /* adjust for curvature */
+ if (viewOptions.doCurv) {
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ e.elev[2] = adjust_for_curvature(*vp, ay, ax, e.elev[2], viewOptions, hd);
+ }
+
+ /*write adjusted elevation into the row of data going through the viewpoint */
+ if (i == vp->row) {
+ (*data)[0][j] = e.elev[0];
+ (*data)[1][j] = e.elev[1];
+ (*data)[2][j] = e.elev[2];
+ }
+
/*put event into event list */
e.eventType = ENTERING_EVENT;
calculate_event_position(e, vp->row, vp->col, &ay, &ax);
@@ -294,8 +370,13 @@
}
G_percent(nrows, nrows, 2);
- G_message(_("...done creating event list"));
G_close_cell(infd);
+
+ G_free(inrast[0]);
+ G_free(inrast[1]);
+ G_free(inrast[2]);
+ G_free(inrast);
+
return nevents;
}
@@ -312,20 +393,25 @@
if data is not NULL, it creates an array that stores all events on
the same row as the viewpoint.
*/
-AMI_STREAM < AEvent > *grass_init_event_list(char *rastName, Viewpoint * vp,
+AMI_STREAM < AEvent > *init_event_list(char *rastName, Viewpoint * vp,
GridHeader * hd,
ViewOptions viewOptions,
- double **data,
+ surface_type ***data,
IOVisibilityGrid * visgrid)
{
+
G_message(_("computing events ..."));
assert(rastName && vp && hd && visgrid);
if (data != NULL) {
/*data is used to store all the cells on the same row as the
//viewpoint. */
- *data = (double *)G_malloc(G_window_cols() * sizeof(double));
+ *data = (surface_type **)G_malloc(3 * sizeof(surface_type *));
assert(*data);
+ (*data)[0] = (surface_type *)G_malloc(3 * G_window_cols() * sizeof(surface_type));
+ assert((*data)[0]);
+ (*data)[1] = (*data)[0] + G_window_cols();
+ (*data)[2] = (*data)[1] + G_window_cols();
}
/*create the event stream that will hold the events */
@@ -347,20 +433,37 @@
RASTER_MAP_TYPE data_type;
- data_type = G_raster_map_type(rastName, mapset);
- void *inrast;
+ /* data_type = G_raster_map_type(rastName, mapset); */
+ data_type = G_SURFACE_TYPE;
+ G_SURFACE_T **inrast;
+ int nrows = G_window_rows();
+ int ncols = G_window_rows();
- inrast = G_allocate_raster_buf(data_type);
+ inrast = (G_SURFACE_T **)G_malloc(3 * sizeof(G_SURFACE_T *));
assert(inrast);
+ inrast[0] = (G_SURFACE_T *)G_allocate_raster_buf(data_type);
+ assert(inrast[0]);
+ inrast[1] = (G_SURFACE_T *)G_allocate_raster_buf(data_type);
+ assert(inrast[1]);
+ inrast[2] = (G_SURFACE_T *)G_allocate_raster_buf(data_type);
+ assert(inrast[2]);
+
+ G_set_null_value(inrast[0], ncols, data_type);
+ G_set_null_value(inrast[1], ncols, data_type);
+ G_set_null_value(inrast[2], ncols, data_type);
/*scan through the raster data */
DCELL d;
int isnull = 0;
dimensionType i, j, k;
+ int row1, col1;
double ax, ay;
+ G_SURFACE_T *elev1, *elev2, *elev3, *elev4;
AEvent e;
- int nrows = G_window_rows();
+ /* read first row */
+ G_get_raster_row(infd, inrast[2], 0, data_type);
+
e.angle = -1;
/*start scanning through the grid */
@@ -369,62 +472,68 @@
G_percent(i, nrows, 2);
/*read in the raster row */
- if (G_get_raster_row(infd, inrast, i, data_type) <= 0)
- G_fatal_error(_("Coord not read from row %d of %s"), i, rastName);
+ int rasterRowResult = 1;
+
+ G_SURFACE_T *tmprast = inrast[0];
+ inrast[0] = inrast[1];
+ inrast[1] = inrast[2];
+ inrast[2] = tmprast;
+ if (i < nrows - 1)
+ rasterRowResult = G_get_raster_row(infd, inrast[2], i + 1, data_type);
+ else
+ G_set_null_value(inrast[2], ncols, data_type);
+
+ if (rasterRowResult <= 0)
+ G_fatal_error(_("Coord not read from row %d of <%s>"), i,
+ rastName);
+
/*fill event list with events from this row */
- for (j = 0; j < G_window_cols(); j++) {
+ for (j = 0; j < ncols; j++) {
e.row = i;
e.col = j;
- /*read the elevation value into the event, depending on data_type */
- switch (data_type) {
- case CELL_TYPE:
- isnull = G_is_c_null_value(&(((CELL *) inrast)[j]));
- e.elev = (float)(((CELL *) inrast)[j]);
- break;
- case FCELL_TYPE:
- isnull = G_is_f_null_value(&(((FCELL *) inrast)[j]));
- e.elev = (float)(((FCELL *) inrast)[j]);
- break;
- case DCELL_TYPE:
- isnull = G_is_d_null_value(&(((DCELL *) inrast)[j]));
- e.elev = (float)(((DCELL *) inrast)[j]);
- break;
- }
+ /*read the elevation value into the event */
+ isnull = G_is_null_value(&(inrast[1][j]), data_type);
+ e.elev[1] = inrast[1][j];
/* adjust for curvature */
- e.elev = adjust_for_curvature(*vp, i, j, e.elev, viewOptions, hd);
+ e.elev[1] = adjust_for_curvature(*vp, i, j, e.elev[1], viewOptions, hd);
if (data != NULL) {
- /**write the row of data going through the viewpoint */
+ /*write the row of data going through the viewpoint */
if (i == vp->row) {
- (*data)[j] = e.elev;
+ (*data)[0][j] = e.elev[1];
+ (*data)[1][j] = e.elev[1];
+ (*data)[2][j] = e.elev[1];
}
}
/* set the viewpoint */
if (i == vp->row && j == vp->col) {
- set_viewpoint_elev(vp, e.elev + viewOptions.obsElev);
+ set_viewpoint_elev(vp, e.elev[1] + viewOptions.obsElev);
+ /*what to do when viewpoint is NODATA */
+ if (is_nodata(hd, e.elev[1])) {
+ G_warning("Viewpoint is NODATA.");
+ G_message("Will assume its elevation is %.f", e.elev[1]);
+ };
if (viewOptions.tgtElev > 0)
vp->target_offset = viewOptions.tgtElev;
else
vp->target_offset = 0.;
- /*what to do when viewpoint is NODATA */
- if (is_nodata(hd, e.elev)) {
- G_warning("Viewpoint is NODATA.");
- G_message("Will assume its elevation is %.f", e.elev);
- };
- }
- /*don't insert viewpoint into eventlist */
- if (i == vp->row && j == vp->col)
+ /* add viewpoint to visibility grid */
+ VisCell visCell = { i, j, 180 };
+ add_result_to_io_visibilitygrid(visgrid, &visCell);
+
+ /*don't insert viewpoint into eventlist */
continue;
+ }
/*don't insert the nodata cell events */
- if (is_nodata(hd, e.elev)) {
+ if (is_nodata(hd, e.elev[1])) {
/* record this cell as being NODATA. ; this is necessary so
that we can distingush invisible events, from nodata
events in the output */
@@ -439,6 +548,36 @@
(*vp, *hd, i, j, viewOptions.maxDist))
continue;
+ /* get ENTER elevation */
+ e.eventType = ENTERING_EVENT;
+ e.elev[0] = calculate_event_elevation(e, nrows, ncols,
+ vp->row, vp->col, inrast, data_type);
+ /* adjust for curvature */
+ if (viewOptions.doCurv) {
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ e.elev[0] = adjust_for_curvature(*vp, ay, ax, e.elev[0], viewOptions, hd);
+ }
+
+ /* get EXIT elevation */
+ e.eventType = EXITING_EVENT;
+ e.elev[2] = calculate_event_elevation(e, nrows, ncols,
+ vp->row, vp->col, inrast, data_type);
+ /* adjust for curvature */
+ if (viewOptions.doCurv) {
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ e.elev[2] = adjust_for_curvature(*vp, ay, ax, e.elev[2], viewOptions, hd);
+ }
+
+ if (data != NULL) {
+
+ /*write the row of data going through the viewpoint */
+ if (i == vp->row) {
+ (*data)[0][j] = e.elev[0];
+ (*data)[1][j] = e.elev[1];
+ (*data)[2][j] = e.elev[2];
+ }
+ }
+
/*put event into event list */
e.eventType = ENTERING_EVENT;
calculate_event_position(e, vp->row, vp->col, &ay, &ax);
@@ -459,14 +598,19 @@
} /* for i */
G_percent(nrows, nrows, 2);
- G_message(_("...done creating event list\n"));
G_close_cell(infd);
- G_message("nbEvents = %lu", (unsigned long)eventList->stream_len());
- G_message("Event stream length: %lu x %dB (%lu MB)",
+
+ G_free(inrast[0]);
+ G_free(inrast[1]);
+ G_free(inrast[2]);
+ G_free(inrast);
+
+ G_debug(1, "nbEvents = %lu", (unsigned long)eventList->stream_len());
+ G_debug(1, "Event stream length: %lu x %dB (%lu MB)",
(unsigned long)eventList->stream_len(), (int)sizeof(AEvent),
(unsigned long)(((long long)(eventList->stream_len() *
sizeof(AEvent))) >> 20));
- fflush(stdout);
+
return eventList;
}
@@ -868,7 +1012,3 @@
G_close_cell(visfd);
return;
}
-
-
-
-#endif
Modified: grass-addons/raster/r.viewshed/grass.h
===================================================================
--- grass-addons/raster/r.viewshed/grass.h 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/grass.h 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -37,15 +36,12 @@
*****************************************************************************/
-#ifdef __GRASS__
-
#ifndef _GRASS_H
#define _GRASS_H
#include <math.h>
extern "C"
{
-#include <grass/config.h>
#include <grass/gis.h>
#include <grass/glocale.h>
}
@@ -60,8 +56,8 @@
curvature of the earth; otherwise return the passed height
unchanged.
*/
-float adjust_for_curvature(Viewpoint vp, dimensionType row,
- dimensionType col, float h,
+float adjust_for_curvature(Viewpoint vp, double row,
+ double col, float h,
ViewOptions viewOptions);
@@ -72,8 +68,14 @@
/*return a GridHeader with all the relevant data filled in from GRASS */
-GridHeader *read_header_from_GRASS(char *rastName, Cell_head * region);
+GridHeader *read_header(char *rastName, Cell_head * region);
+/* calculate ENTER and EXIT event elevation */
+surface_type calculate_event_elevation(AEvent e, int nrows, int ncols,
+ dimensionType vprow, dimensionType vpcol,
+ G_SURFACE_T **inrast, RASTER_MAP_TYPE data_type);
+
+
/* ************************************************************ */
/* input: an array capable to hold the max number of events, a raster
name, a viewpoint and the viewOptions; action: figure out all events
@@ -83,9 +85,9 @@
AEvent* with all the events for the map. Used when solving in
memory, so the AEvent* should fit in memory. */
size_t
-grass_init_event_list_in_memory(AEvent * eventList, char *rastName,
+init_event_list_in_memory(AEvent * eventList, char *rastName,
Viewpoint * vp, GridHeader * hd,
- ViewOptions viewOptions, double **data,
+ ViewOptions viewOptions, surface_type ***data,
MemoryVisibilityGrid * visgrid);
@@ -99,10 +101,10 @@
if data is not NULL, it creates an array that stores all events on
the same row as the viewpoint.
*/
-AMI_STREAM < AEvent > *grass_init_event_list(char *rastName, Viewpoint * vp,
+AMI_STREAM < AEvent > *init_event_list(char *rastName, Viewpoint * vp,
GridHeader * hd,
ViewOptions viewOptions,
- double **data,
+ surface_type ***data,
IOVisibilityGrid * visgrid);
@@ -149,5 +151,3 @@
char *visfname, float vp_elev);
#endif/*_GRASS_H*/
-
-#endif /*__GRASS__*/
Modified: grass-addons/raster/r.viewshed/grid.cc
===================================================================
--- grass-addons/raster/r.viewshed/grid.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/grid.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -43,72 +42,19 @@
#include <stdlib.h>
#include <assert.h>
-#ifdef __GRASS__
extern "C"
{
#include <grass/config.h>
#include <grass/gis.h>
#include <grass/glocale.h>
}
-#endif
#include "grid.h"
-/* ------------------------------------------------------------ */
-/*read header from file and return it; */
-GridHeader *read_header_from_arcascii_file(char *fname)
-{
- assert(fname);
- FILE *fp;
- fp = fopen(fname, "r");
- assert(fp);
- GridHeader *hd = read_header_from_arcascii_file(fp);
-
- fclose(fp);
- return hd;
-}
-
-
-
/* ------------------------------------------------------------ */
-/*read header from file and return it; */
-GridHeader *read_header_from_arcascii_file(FILE * fp)
-{
- assert(fp);
- rewind(fp);
-
- GridHeader *hd = (GridHeader *) malloc(sizeof(GridHeader));
-
- assert(hd);
-
- int nrows, ncols;
-
- fscanf(fp, "%*s%d\n", &ncols);
- fscanf(fp, "%*s%d\n", &nrows);
- /*check that you dont lose precision */
- if (nrows <= maxDimension && ncols <= maxDimension) {
- hd->nrows = (dimensionType) nrows;
- hd->ncols = (dimensionType) ncols;
- }
- else {
- fprintf(stderr, "grid dimension too big for current precision\n");
- printf("change type and re-compile\n");
- exit(1);
- }
- fscanf(fp, "%*s%f\n", &(hd->xllcorner));
- fscanf(fp, "%*s%f\n", &(hd->yllcorner));
- fscanf(fp, "%*s%f\n", &(hd->cellsize));
- fscanf(fp, "%*s%f\n", &(hd->nodata_value));
-
- return hd;
-}
-
-
-
-/* ------------------------------------------------------------ */
/*copy from b to a */
void copy_header(GridHeader * a, GridHeader b)
{
@@ -117,56 +63,25 @@
a->ncols = b.ncols;
a->xllcorner = b.xllcorner;
a->yllcorner = b.yllcorner;
- a->cellsize = b.cellsize;
+ a->ns_res = b.ns_res;
+ a->ew_res = b.ew_res;
a->nodata_value = b.nodata_value;
return;
}
/* ------------------------------------------------------------ */
-/*print header */
-void print_grid_header(GridHeader * hd)
+/*returns 1 if value is Nodata, 0 if it is not */
+int is_nodata(GridHeader * hd, float value)
{
-
assert(hd);
- print_grid_header(stdout, hd);
- return;
-}
+ return G_is_null_value(&value, G_SURFACE_TYPE);
-/* ------------------------------------------------------------ */
-void print_grid_header(FILE * fp, GridHeader * hd)
-{
- assert(fp && hd);
- fprintf(fp, "ncols\t%d\n", hd->ncols);
- fprintf(fp, "nrows\t%d\n", hd->nrows);
- fprintf(fp, "xllcorner\t%f\n", hd->xllcorner);
- fprintf(fp, "yllcorner\t%f\n", hd->yllcorner);
- fprintf(fp, "cellsize\t%f\n", hd->cellsize);
- fprintf(fp, "NODATA_value\t%f\n", hd->nodata_value);
- return;
}
-
-
-
/* ------------------------------------------------------------ */
/*returns 1 if value is Nodata, 0 if it is not */
-int is_nodata(GridHeader * hd, float value)
-{
- assert(hd);
-#ifdef __GRASS__
- return G_is_f_null_value(&value);
-#else
- if (fabs(value - hd->nodata_value) < 0.000001) {
- return 1;
- }
- return 0;
-#endif
-}
-
-/* ------------------------------------------------------------ */
-/*returns 1 if value is Nodata, 0 if it is not */
int is_nodata(Grid * grid, float value)
{
assert(grid);
@@ -181,11 +96,7 @@
Grid *create_empty_grid()
{
-#ifdef __GRASS__
Grid *ptr_grid = (Grid *) G_malloc(sizeof(Grid));
-#else
- Grid *ptr_grid = (Grid *) malloc(sizeof(Grid));
-#endif
assert(ptr_grid);
@@ -212,23 +123,15 @@
assert(pgrid);
assert(pgrid->hd);
-#ifdef __GRASS__
pgrid->grid_data = (float **)G_malloc(pgrid->hd->nrows * sizeof(float *));
-#else
- pgrid->grid_data = (float **)malloc(pgrid->hd->nrows * sizeof(float *));
-#endif
+
assert(pgrid->grid_data);
dimensionType i;
for (i = 0; i < pgrid->hd->nrows; i++) {
-#ifdef __GRASS__
pgrid->grid_data[i] =
(float *)G_malloc(pgrid->hd->ncols * sizeof(float));
-#else
- pgrid->grid_data[i] =
- (float *)malloc(pgrid->hd->ncols * sizeof(float));
-#endif
assert(pgrid->grid_data[i]);
}
@@ -241,63 +144,7 @@
return;
}
-/* ------------------------------------------------------------ */
-/*reads header and data from file */
-Grid *read_grid_from_arcascii_file(char *filename)
-{
- assert(filename);
- FILE *fp = fopen(filename, "r");
-
- assert(fp);
-
- Grid *grid = create_empty_grid();
-
- grid->hd = read_header_from_arcascii_file(fp);
- alloc_grid_data(grid);
-
- /*READ DATA */
- dimensionType i, j;
- int first_flag = 1;
- float value = 0;
-
- for (i = 0; i < grid->hd->nrows; i++) {
- for (j = 0; j < grid->hd->ncols; j++) {
- fscanf(fp, "%f", &value);
- grid->grid_data[i][j] = value;
-
- if (first_flag) {
- if (is_nodata(grid, value))
- continue;
- grid->minvalue = grid->maxvalue = value;
- first_flag = 0;
- }
- else {
- if (is_nodata(grid, value))
- continue;
- if (value > grid->maxvalue)
- grid->maxvalue = value;
- if (value < grid->minvalue)
- grid->minvalue = value;
- }
- }
- fscanf(fp, "\n");
- }
-
- fclose(fp);
-
-#ifdef DEBUG_ON
- printf("**DEBUG: readGridFromArcasciiFile():\n");
- fflush(stdout);
-#endif
-
- return grid;
-}
-
-
-
-
-
/* ------------------------------------------------------------ */
/*destroy the structure and reclaim all memory allocated */
void destroy_grid(Grid * grid)
@@ -305,7 +152,6 @@
assert(grid);
/*free grid data if its allocated */
-#ifdef __GRASS__
if (grid->grid_data) {
dimensionType i;
@@ -320,23 +166,6 @@
G_free(grid->hd);
G_free(grid);
-
-#else
- if (grid->grid_data) {
- dimensionType i;
-
- for (i = 0; i < grid->hd->nrows; i++) {
- if (!grid->grid_data[i])
- free((float *)grid->grid_data[i]);
- }
-
- free((float **)grid->grid_data);
- }
-
- free(grid->hd);
- free(grid);
-#endif
-
#ifdef _DEBUG_ON
printf("**DEBUG: grid destroyed.\n");
fflush(stdout);
@@ -344,60 +173,3 @@
return;
}
-
-
-
-
-
-
-
-/* ------------------------------------------------------------ */
-/*save the grid into an arcascii file. Loops through all elements x
- in row-column order and writes fun(x) to file */
-void
-save_grid_to_arcascii_file(Grid * grid, char *filename, float (*fun) (float))
-{
-
- assert(filename && grid);
- printf("saving grid to %s\n", filename);
- fflush(stdout);
-
- FILE *outfile = fopen(filename, "r");
-
- if (outfile) { /*outfile already exists */
- printf("The output file already exists. It will be overwritten\n");
- fclose(outfile);
- int ret = remove(filename); /*delete the existing file */
-
- if (ret != 0) {
- printf("unable to overwrite the existing output file.\n!");
- exit(1);
- }
- }
- FILE *fp = fopen(filename, "a");
-
- assert(fp);
-
- /*print header */
- print_grid_header(fp, grid->hd);
-
- /*print data */
- dimensionType i, j;
-
- for (i = 0; i < grid->hd->nrows; i++) {
- for (j = 0; j < grid->hd->ncols; j++) {
-
- /* call fun() on this element and write it to file */
- fprintf(fp, "%.1f ", fun(grid->grid_data[i][j]));
-
- }
- fprintf(fp, "\n");
- }
-
-#ifdef _DEBUG_ON
- printf("**DEBUG: saveGridToArcasciiFile: saved to %s\n", filename);
- fflush(stdout);
-#endif
-
- return;
-}
Modified: grass-addons/raster/r.viewshed/grid.h
===================================================================
--- grass-addons/raster/r.viewshed/grid.h 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/grid.h 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -47,13 +46,15 @@
#include <stdio.h>
#include <limits.h>
-#ifdef __GRASS__
extern "C"
{
#include <grass/gis.h>
}
-#endif
+#define G_SURFACE_TYPE FCELL_TYPE
+typedef float surface_type;
+typedef FCELL G_SURFACE_T;
+
/* this should accomodate grid sizes up to 2^16-1=65,535
If this is not enough, change type and recompile */
typedef unsigned short int dimensionType;
@@ -64,15 +65,13 @@
{
dimensionType ncols; /*number of columns in the grid */
dimensionType nrows; /*number of rows in the grid */
- float xllcorner; /*xllcorner refers to the western edge of grid */
- float yllcorner; /*yllcorner refers to the southern edge of grid */
- float cellsize; /*the resolution of the grid */
- float nodata_value; /*the value that represents missing data */
-
-#ifdef __GRASS__
- /* GRASS */
+ double xllcorner; /*xllcorner refers to the western edge of grid */
+ double yllcorner; /*yllcorner refers to the southern edge of grid */
+ double ew_res; /*the ew resolution of the grid */
+ double ns_res; /*the ns resolution of the grid */
+ surface_type nodata_value; /*the value that represents missing data */
+
struct Cell_head window;
-#endif
} GridHeader;
@@ -90,22 +89,10 @@
-
-/* create and return the header of the grid stored in this file; */
-GridHeader *read_header_from_arcascii_file(char *fname);
-
-/* create and return the header of the grid stored in this file; */
-GridHeader *read_header_from_arcascii_file(FILE * fp);
-
/*copy from b to a */
void copy_header(GridHeader * a, GridHeader b);
-/*print header */
-void print_grid_header(GridHeader * hd);
-void print_grid_header(FILE * fp, GridHeader * hd);
-
-
/*returns 1 if value is Nodata, 0 if it is not */
int is_nodata(GridHeader * hd, float value);
int is_nodata(Grid * grid, float value);
@@ -116,17 +103,8 @@
/*allocate memory for grid data, grid must have a header */
void alloc_grid_data(Grid * grid);
-/*scan an arcascii file and fill the information in the given structure */
-Grid *read_grid_from_arcascii_file(char *filename);
-
/*destroy the structure and reclaim all memory allocated */
void destroy_grid(Grid * grid);
-/*save grid into an arcascii file; Loops through all elements x in
- row-column order and writes fun(x) to file */
-void save_grid_to_arcascii_file(Grid * grid, char *filename,
- float (*fun) (float));
-
-
#endif
Modified: grass-addons/raster/r.viewshed/main.cc
===================================================================
--- grass-addons/raster/r.viewshed/main.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/main.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -43,7 +42,6 @@
#include <ctype.h>
#include <unistd.h>
-#ifdef __GRASS__
extern "C"
{
#include <grass/config.h>
@@ -53,17 +51,12 @@
#include "grass.h"
#include <grass/iostream/ami.h>
-#else
-#include <ami.h>
-#endif
-
#include "viewshed.h"
#include "visibility.h"
#include "grid.h"
#include "rbbst.h"
#include "statusstructure.h"
#include "distribute.h"
-#include "print_message.h"
@@ -102,18 +95,10 @@
Rtimer totalTime);
void print_timings_external_memory(Rtimer totalTime, Rtimer viewshedTime,
Rtimer outputTime, Rtimer sortOutputTime);
-void print_status(Viewpoint vp, ViewOptions viewOptions,
- long long memSizeBytes);
-void print_usage();
-#ifdef __GRASS__
void parse_args(int argc, char *argv[], int *vpRow, int *vpCol,
ViewOptions * viewOptions, long long *memSizeBytes,
Cell_head * window);
-#else
-void parse_args(int argc, char *argv[], int *vpRow, int *vpCol,
- ViewOptions * viewOptions, long long *memSizeBytes);
-#endif
@@ -125,7 +110,6 @@
int main(int argc, char *argv[])
{
-#ifdef __GRASS__
/* GRASS initialization stuff */
struct GModule *module;
@@ -142,10 +126,8 @@
if (G_get_set_window(®ion) == -1)
G_fatal_error("Error getting current region");
-#endif
-
/* ************************************************************ */
/* parameters set up */
long long memSizeBytes = DEFAULT_MEMORY;
@@ -172,14 +154,9 @@
viewOptions.outputMode = OUTPUT_ANGLE;
viewOptions.doCurv = 0;
-#ifdef __GRASS__
parse_args(argc, argv, &vpRow, &vpCol, &viewOptions, &memSizeBytes,
®ion);
-#else
- parse_args(argc, argv, &vpRow, &vpCol, &viewOptions, &memSizeBytes);
-#endif
-
/* set viewpoint with the coordinates specified by user. The
height of the viewpoint is not known at this point---it will be
set during the execution of the algorithm */
@@ -187,16 +164,13 @@
set_viewpoint_coord(&vp, vpRow, vpCol);
- print_status(vp, viewOptions, memSizeBytes);
-
/* ************************************************************ */
/* set up the header of the raster with all raster info and make
sure the requested viewpoint is on the map */
GridHeader *hd;
-#ifdef __GRASS__
- hd = read_header_from_GRASS(viewOptions.inputfname, ®ion);
+ hd = read_header(viewOptions.inputfname, ®ion);
assert(hd);
G_get_set_window(&(hd->window));
@@ -208,28 +182,9 @@
G_warning(_("viewpont: (row=%d, col=%d)"), vp.row, vp.col);
G_fatal_error(_("grid: (rows=%d, cols=%d)"), hd->nrows, hd->ncols);
}
-#else
- /*open file input file and read grid header from grid ascii file */
- hd = read_header_from_arcascii_file(viewOptions.inputfname);
- assert(hd);
- printf("input grid: (rows=%d, cols=%d)\n", hd->nrows, hd->ncols);
- fflush(stdout);
- /*sanity check */
- if (!(vp.row < hd->nrows && vp.col < hd->ncols)) {
- printf("Viewpoint outside grid\n");
- printf("viewpont: (row=%d, col=%d)\n", vp.row, vp.col);
- printf("grid: (rows=%d, cols=%d)\n", hd->nrows, hd->ncols);
- exit(1);
- /* LT: there is no need to exit if viewpoint is outside grid,
- the algorithm will work correctly in theory. But this
- requires some changes. To do. */
- }
-#endif /*__GRASS__*/
-
/* set curvature params */
-#ifdef __GRASS__
viewOptions.cellsize = region.ew_res;
double e2;
@@ -242,42 +197,39 @@
}
G_begin_distance_calculations();
-
-#else
- /* in standalone mode we do not know how to adjust for curvature */
- assert(viewOptions.doCurv == 0);
- viewOptions.ellps_a = 0;
-#endif
+
/* ************************************************************ */
/* decide whether the computation of the viewshed will take place
in-memory or in external memory */
int IN_MEMORY;
long long inmemSizeBytes = get_viewshed_memory_usage(hd);
- printf("In-memory memory usage is %lld B (%d MB), \
-max mem allowed=%lld B(%dMB)\n", inmemSizeBytes, (int)(inmemSizeBytes >> 20), memSizeBytes, (int)(memSizeBytes >> 20));
+ G_verbose_message(_("In-memory memory usage is %lld B (%d MB), \
+ max mem allowed=%lld B(%dMB)"), inmemSizeBytes,
+ (int)(inmemSizeBytes >> 20), memSizeBytes,
+ (int)(memSizeBytes >> 20));
if (inmemSizeBytes < memSizeBytes) {
IN_MEMORY = 1;
- print_message("*************\nIN_MEMORY MODE\n*************\n");
+ G_verbose_message("***** IN_MEMORY MODE *****");
}
else {
- print_message("*************\nEXTERNAL_MEMORY MODE\n**********\n");
+ G_verbose_message("***** EXTERNAL_MEMORY MODE *****");
IN_MEMORY = 0;
}
- fflush(stdout);
+
/* the mode can be forced to in memory or external if the user
wants to test or debug a specific mode */
#ifdef FORCE_EXTERNAL
IN_MEMORY = 0;
- print_message("FORCED EXTERNAL\n");
+ G_debug(1, "FORCED EXTERNAL");
#endif
#ifdef FORCE_INTERNAL
IN_MEMORY = 1;
- print_message("FORCED INTERNAL\n");
+ G_debug(1, "FORCED INTERNAL");
#endif
@@ -324,27 +276,27 @@
if (getenv(STREAM_TMPDIR) != NULL) {
/*if already set */
- fprintf(stderr, "%s=%s\n", STREAM_TMPDIR, getenv(STREAM_TMPDIR));
- printf("Intermediate stream location: %s\n",
+ G_debug(1, "%s=%s", STREAM_TMPDIR, getenv(STREAM_TMPDIR));
+ G_debug(1, "Intermediate stream location: %s",
getenv(STREAM_TMPDIR));
}
else {
/*set it */
sprintf(buf, "%s=%s", STREAM_TMPDIR, viewOptions.streamdir);
- fprintf(stderr, "setting %s ", buf);
+ G_debug(1, "setting %s ", buf);
putenv(buf);
if (getenv(STREAM_TMPDIR) == NULL) {
- fprintf(stderr, ", not set\n");
+ G_fatal_error(_("%s not set"), "STREAM_TMPDIR");
exit(1);
}
else {
- fprintf(stderr, ", ok.\n");
+ G_debug(1, "are ok.");
}
- printf("Intermediate stream location: %s\n", viewOptions.streamdir);
+ G_debug(1, "Intermediate stream location: %s", viewOptions.streamdir);
}
- fprintf(stderr, "Intermediate files will not be deleted "
- "in case of abnormal termination.\n");
- fprintf(stderr, "To save space delete these files manually!\n");
+ G_important_message(_("Intermediate files will not be deleted \
+ in case of abnormal termination."));
+ G_important_message(_("To save space delete these files manually!"));
/* initialize IOSTREAM memory manager */
@@ -377,8 +329,8 @@
/* external memory, base case */
/* ************************************************************ */
if (BASE_CASE) {
- print_message
- ("---Active structure small, starting base case---\n");
+ G_debug
+ (1, "---Active structure small, starting base case---");
Rtimer totalTime, viewshedTime, outputTime, sortOutputTime;
@@ -419,9 +371,9 @@
/************************************************************ */
else { /* if not BASE_CASE */
#ifndef FORCE_DISTRIBUTION
- print_message("---Active structure does not fit in memory,");
+ G_debug(1, "---Active structure does not fit in memory,");
#else
- print_message("FORCED DISTRIBUTION\n");
+ G_debug(1, "FORCED DISTRIBUTION");
#endif
Rtimer totalTime, sweepTime, outputTime, sortOutputTime;
@@ -459,15 +411,11 @@
/*end external memory, distribution sweep */
+ /**************************************/
+ /* FINISH UP, ALL CASES */
+ /**************************************/
-
-
-
- /* ************************************************************ */
- /* FINISH UP, ALL CASES */
- /* ************************************************************ */
/*close input file and free grid header */
-#ifdef __GRASS__
G_free(hd);
/*following GRASS's coding standards for history and exiting */
struct History history;
@@ -476,20 +424,13 @@
G_command_history(&history);
G_write_history(viewOptions.outputfname, &history);
exit(EXIT_SUCCESS);
-#else
- free(hd);
-#endif
}
-
-
-
/* ------------------------------------------------------------ */
-/*There are two versions of the function parse_args, one for when it is in GRASS, one for independent compilition. */
-#ifdef __GRASS__
+/* parse arguments */
void
parse_args(int argc, char *argv[], int *vpRow, int *vpCol,
ViewOptions * viewOptions, long long *memSizeBytes,
@@ -544,7 +485,7 @@
elevationFlag = G_define_flag();
elevationFlag->key = 'e';
elevationFlag->description =
- _("Output format is {NODATA, -1 (invisible), elev-viewpoint_elev (visible)}");
+ _("Output format is invisible = NULL, else current elev - viewpoint_elev");
/* viewpoint coordinates */
struct Option *viewLocOpt;
@@ -602,12 +543,12 @@
struct Option *memAmountOpt;
memAmountOpt = G_define_option();
- memAmountOpt->key = "memory_usage";
+ memAmountOpt->key = "memory";
memAmountOpt->type = TYPE_INTEGER;
memAmountOpt->required = NO;
memAmountOpt->key_desc = "value";
memAmountOpt->description =
- _("The amount of main memory in MB to be used");
+ _("Amount of memory to be used in MB");
memAmountOpt->answer = "500";
/* temporary STREAM path */
@@ -641,7 +582,7 @@
viewOptions->maxDist = atof(maxDistOpt->answer);
if (viewOptions->maxDist < 0 && viewOptions->maxDist != INFINITY_DISTANCE) {
- G_fatal_error(_("Negative max distance value is not valid"));
+ G_fatal_error(_("A negative max distance value is not allowed"));
}
@@ -657,8 +598,9 @@
int memSizeMB = atoi(memAmountOpt->answer);
if (memSizeMB < 0) {
- printf("Memory cannot be negative.\n");
- exit(1);
+ G_warning(_("Amount of memory cannot be negative."));
+ G_warning(_(" Converting %d to %d MB"), memSizeMB, -memSizeMB);
+ memSizeMB = -memSizeMB;
}
*memSizeBytes = (long long)memSizeMB;
*memSizeBytes = (*memSizeBytes) << 20;
@@ -672,13 +614,13 @@
if (row_col->answer) {
*vpRow = atoi(viewLocOpt->answers[1]);
*vpCol = atoi(viewLocOpt->answers[0]);
- printf("viewpoint in row-col mode: (%d,%d)\n", *vpRow, *vpCol);
+ G_debug(1, "viewpoint in row-col mode: (%d,%d)", *vpRow, *vpCol);
}
else {
*vpRow = (int)G_northing_to_row(atof(viewLocOpt->answers[1]), window);
*vpCol = (int)G_easting_to_col(atof(viewLocOpt->answers[0]), window);
- printf("viewpoint converted from lat-lon mode: (%d,%d)\n", *vpRow,
- *vpCol);
+ G_debug(1, "viewpoint converted from current projection: (%.3f, %.3f) to row, col (%d, %d)",
+ atof(viewLocOpt->answers[1]), atof(viewLocOpt->answers[0]), *vpRow, *vpCol);
}
@@ -690,105 +632,6 @@
/* ------------------------------------------------------------ */
-#else /* if not in GRASS mode */
-void
-parse_args(int argc, char *argv[], int *vpRow, int *vpCol,
- ViewOptions * viewOptions, long long *memSizeBytes)
-{
-
- assert(vpRow && vpCol && viewOptions && memSizeBytes);
-
- int gotrow = 0, gotcol = 0, gotinput = 0, gotoutput = 0;
- int c;
-
- /*deal with flags for the output using getopt */
- while ((c = getopt(argc, argv, "i:o:r:c:v:e:d:m:")) != -1) {
- switch (c) {
- case 'i':
- /* inputfile name */
- //*inputfname = optarg;
- strcpy(viewOptions->inputfname, optarg);
- gotinput = 1;
- break;
- case 'o':
- //*outputfname = optarg;
- strcpy(viewOptions->outputfname, optarg);
- gotoutput = 1;
- break;
- case 'r':
- *vpRow = atoi(optarg);
- gotrow = 1;
- break;
- case 'c':
- *vpCol = atoi(optarg);
- gotcol = 1;
- break;
- case 'v':
- /* output mode */
- if (strcmp(optarg, "angle") == 0)
- viewOptions->outputMode = OUTPUT_ANGLE;
- else if (strcmp(optarg, "bool") == 0)
- viewOptions->outputMode = OUTPUT_BOOL;
- else if (strcmp(optarg, "elev") == 0)
- viewOptions->outputMode = OUTPUT_ELEV;
- else {
- printf("unknown option %s: use -v: [angle|bool|elev]\n",
- optarg);
- exit(1);
- }
- break;
- case 'e':
- viewOptions->obsElev = atof(optarg);
- break;
- case 'd':
- viewOptions->maxDist = atof(optarg);
- if (viewOptions->maxDist < 0) {
- printf("max distance needs to be positive\n");
- exit(EXIT_FAILURE);
- }
- break;
- case 'm':
- int memSizeMB;
-
- memSizeMB = atoi(optarg);
- if (memSizeMB < 0) {
- printf("Memory cannot be negative.\n");
- exit(1);
- }
- *memSizeBytes = (long long)memSizeMB;
- *memSizeBytes = (*memSizeBytes) << 20;
- break;
- case '?':
- if (optopt == 'i' || optopt == 'o' || optopt == 'r' ||
- optopt == 'c' || optopt == 'e' || optopt == 'd' ||
- optopt == 'm')
- fprintf(stderr, "Option -%c requires an argument.\n", optopt);
- else if (isprint(optopt))
- fprintf(stderr, "Unknown option '-%c'.\n", optopt);
- else
- fprintf(stderr, "Unknown option character '\\x%x.\n", optopt);
- print_usage();
- exit(1);
- //default:
- //exit(1);
- }
- } /* while getopt */
-
- /*check to make sure the required options are set. */
- if (!(gotinput && gotoutput && gotrow && gotcol)) {
- printf("Not all required options set.\n");
- print_usage();
- exit(1);
- }
- printf("viewpoint: (%d, %d)\n", *vpRow, *vpCol);
- return;
-}
-#endif
-
-
-
-
-/* ------------------------------------------------------------ */
/*print the timings for the internal memory method of computing the
viewshed */
void
@@ -797,21 +640,19 @@
char timeused[100];
- printf("TOTAL TIMING: \n");
+ G_verbose_message("TOTAL TIMING:");
+
rt_sprint_safe(timeused, sweepTime);
- printf("\t%30s", "sweep:");
- printf(timeused);
- printf("\n");
+ G_verbose_message("sweep: %s", timeused);
+ G_verbose_message("\n");
rt_sprint_safe(timeused, outputTime);
- printf("\t%30s", "output:");
- printf(timeused);
- printf("\n");
+ G_verbose_message("output: %s", timeused);
+ G_verbose_message("\n");
rt_sprint_safe(timeused, totalTime);
- printf("\t%30s", "total:");
- printf(timeused);
- printf("\n");
+ G_verbose_message("total: %s", timeused);
+ G_verbose_message("\n");
}
@@ -825,123 +666,19 @@
/*print timings */
char timeused[100];
- printf("\n\nTOTAL TIMING: \n");
+ G_verbose_message("\n\nTOTAL TIMING:");
+
rt_sprint_safe(timeused, viewshedTime);
- printf("\t%30s", "total sweep:");
- printf(timeused);
- printf("\n");
+ G_verbose_message("total sweep: %s", timeused);
+
rt_sprint_safe(timeused, sortOutputTime);
- printf("\t%30s", "sort output:");
- printf(timeused);
- printf("\n");
+ G_verbose_message("sort output: %s", timeused);
+
rt_sprint_safe(timeused, outputTime);
- printf("\t%30s", "Write result grid:");
- printf(timeused);
- printf("\n");
+ G_verbose_message("Write result grid: %s", timeused);
+
rt_sprint_safe(timeused, totalTime);
- printf("\t%30s", "Total Time:");
- printf(timeused);
- printf("\n\n");
+ G_verbose_message("Total Time: %s", timeused);
+ G_verbose_message("\n\n");
return;
}
-
-/* ------------------------------------------------------------ */
-void print_status(Viewpoint vp, ViewOptions viewOptions,
- long long memSizeBytes)
-{
-
-
-#ifdef __GRASS__
- G_message(_("Options set as:\n"));
- G_message(_("---input: %s \n---output: %s \n---viewpoint: (%d, %d)"),
- viewOptions.inputfname, viewOptions.outputfname,
- vp.row, vp.col);
- if (viewOptions.outputMode == OUTPUT_ANGLE) {
- G_message(_("---outputting viewshed in angle mode:"));
- G_message(_("---The output is {NODATA, %d(invisible),angle(visible)}.\n"),
- INVISIBLE);
- }
- if (viewOptions.outputMode == OUTPUT_BOOL) {
- G_message(_("---outputting viewshed in boolean mode: "));
- G_message(_("---The output is {%d (invisible), %d (visible)}.\n"),
- BOOL_INVISIBLE, BOOL_VISIBLE);
- }
- if (viewOptions.outputMode == OUTPUT_ELEV) {
- G_message(_("---outputting viewshed in elevation mode: "));
- G_message(_("---The output is {NODATA, %d (invisible), elev (visible)}.\n"),
- INVISIBLE);
- }
- G_message(_("---observer elevation above terrain: %f\n"),
- viewOptions.obsElev);
-
- if (viewOptions.maxDist == INFINITY_DISTANCE)
- G_message(_("---max distance: infinity\n"));
- else
- G_message(_("---max distance: %f\n"), viewOptions.maxDist);
-
- G_message(_("---consider earth curvature: %d\n"), viewOptions.doCurv);
-
- G_message(_("---max memory = %d MB\n"), (int)(memSizeBytes >> 20));
- G_message(_("---------------------------------\n"));
- G_message(_("---temporary files streamdir: %s\n"), viewOptions.streamdir);
-
-#else
- printf("---------------------------------\nOptions set as:\n");
- printf("input: %s \noutput: %s\n",
- viewOptions.inputfname, viewOptions.outputfname);
- printf("viewpoint: (%d, %d)\n", vp.row, vp.col);
- if (viewOptions.outputMode == OUTPUT_ANGLE) {
- printf("outputting viewshed in angle mode: ");
- printf("The output is {NODATA, %d (invisible), angle (visible)}.\n",
- INVISIBLE);
- }
- if (viewOptions.outputMode == OUTPUT_BOOL) {
- printf("outputting viewshed in boolean mode: ");
- printf("The output is {%d (invisible), %d (visible)}.\n",
- BOOL_INVISIBLE, BOOL_VISIBLE);
- }
- if (viewOptions.outputMode == OUTPUT_ELEV) {
- printf("outputting viewshed in elevation mode: ");
- printf("The output is {NODATA, %d (invisible), elev (visible)}.\n",
- INVISIBLE);
- }
-
- printf("observer elevation above terrain: %f\n", viewOptions.obsElev);
-
- if (viewOptions.maxDist == INFINITY_DISTANCE)
- printf("max distance: infinity\n");
- else
- printf("max distance: %f\n", viewOptions.maxDist);
-
- printf("consider earth curvature: %d\n", viewOptions.doCurv);
-
- printf("max memory = %d MB\n", (int)(memSizeBytes >> 20));
- printf("---------------------------------\n");
- printf("---temporary files streamdir: %s\n", viewOptions.streamdir);
-
-#endif
-
-}
-
-/* ------------------------------------------------------------ */
-/*print the usage information. Only used in the stand-alone version */
-void print_usage()
-{
- printf
- ("\nusage: ioviewshed -i <input name> -o <output name> -r <row number> -c <column number> [-v <angle | bool | elev>] [-e <observer elevation>] [-d <max distance>] [-m <memory usage MB>]\n\n");
-
- printf("OPTIONS\n");
- printf("-i \t input map name.\n");
- printf("-o \t output map name.\n");
- printf("-r \t row number.\n");
- printf("-c \t column number.\n");
- printf("-v \t iutput mode. Default is angle.\n");
- printf
- (" \t\t angle: output is {NODATA, -1 (invisible), angle (visible)}\n\t\t\t angle is a value in [0,180] and represents the vertical angle wrt viewpoint.\n");
- printf(" \t\t bool: output is {0 (invisible), 1 (visible)}.\n");
- printf
- (" \t\t elev: output is {NODATA, -1 (invisible), elev (visible)}. This is not implemented in the standalone version.\n");
- printf("-e \t observer elevation. Default is 0.\n");
- printf("-d \t maximum distance. Default is infinity.\n");
- printf("-m \t memory usage in MB. Default is 500.\n");
-}
Deleted: grass-addons/raster/r.viewshed/print_message.cc
===================================================================
--- grass-addons/raster/r.viewshed/print_message.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/print_message.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -1,67 +0,0 @@
-
-/****************************************************************************
- *
- * MODULE: r.viewshed
- *
- * AUTHOR(S): Laura Toma, Bowdoin College - ltoma at bowdoin.edu
- * Yi Zhuang - yzhuang at bowdoin.edu
- *
- * Ported to GRASS by William Richard -
- * wkrichar at bowdoin.edu or willster3021 at gmail.com
- *
- * Date: July 2008
- *
- * PURPOSE: To calculate the viewshed (the visible cells in the
- * raster) for the given viewpoint (observer) location. The
- * visibility model is the following: Two points in the raster are
- * considered visible to each other if the cells where they belong are
- * visible to each other. Two cells are visible to each other if the
- * line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
- * terms of CPU operations and I/O operations. It has worst-case
- * complexity O(n lg n) in the RAM model and O(sort(n)) in the
- * I/O-model. For the algorithm and all the other details see the
- * paper: "Computing Visibility on * Terrains in External Memory" by
- * Herman Haverkort, Laura Toma and Yi Zhuang.
- *
- * COPYRIGHT: (C) 2008 by the GRASS Development Team
- *
- * This program is free software under the GNU General Public License
- * (>=v2). Read the file COPYING that comes with GRASS for details.
- *
- *****************************************************************************/
-
-
-#include <assert.h>
-#include <stdio.h>
-
-#ifdef __GRASS__
-extern "C"
-{
-#include <grass/config.h>
-#include <grass/gis.h>
-#include <grass/glocale.h>
-}
-#endif
-
-#include "print_message.h"
-
-
-/* Used to make printing messages following the GRASS standards cleaner
- * and easier to read. Needs to be in seperate file so that both GRASS
- * and non-GRASS modes can access it.
- */
-void print_message(char *message)
-{
- assert(message);
-#ifdef __GRASS__
- G_message(_(message));
-#else
- printf(message);
-#endif
-}
Deleted: grass-addons/raster/r.viewshed/print_message.h
===================================================================
--- grass-addons/raster/r.viewshed/print_message.h 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/print_message.h 2011-05-06 07:30:56 UTC (rev 46199)
@@ -1,46 +0,0 @@
-
-/****************************************************************************
- *
- * MODULE: r.viewshed
- *
- * AUTHOR(S): Laura Toma, Bowdoin College - ltoma at bowdoin.edu
- * Yi Zhuang - yzhuang at bowdoin.edu
-
- * Ported to GRASS by William Richard -
- * wkrichar at bowdoin.edu or willster3021 at gmail.com
- *
- * Date: july 2008
- *
- * PURPOSE: To calculate the viewshed (the visible cells in the
- * raster) for the given viewpoint (observer) location. The
- * visibility model is the following: Two points in the raster are
- * considered visible to each other if the cells where they belong are
- * visible to each other. Two cells are visible to each other if the
- * line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
- * terms of CPU operations and I/O operations. It has worst-case
- * complexity O(n lg n) in the RAM model and O(sort(n)) in the
- * I/O-model. For the algorithm and all the other details see the
- * paper: "Computing Visibility on * Terrains in External Memory" by
- * Herman Haverkort, Laura Toma and Yi Zhuang.
- *
- * COPYRIGHT: (C) 2008 by the GRASS Development Team
- *
- * This program is free software under the GNU General Public License
- * (>=v2). Read the file COPYING that comes with GRASS for details.
- *
- *****************************************************************************/
-
-
-#ifndef __PRINT_MESSAGE_H
-#define __PRINT_MESSAGE_H
-
-/*used to make printing messages following the GRASS standards cleaner and easier to read. Needs to be in seperate file so that both GRASS and non-GRASS modes can access it. */
-void print_message(char *message);
-
-#endif
Modified: grass-addons/raster/r.viewshed/r.viewshed.html
===================================================================
--- grass-addons/raster/r.viewshed/r.viewshed.html 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/r.viewshed.html 2011-05-06 07:30:56 UTC (rev 46199)
@@ -32,8 +32,9 @@
<li> -1, if the point is not visible
-<li> a value in [0, 180] representing the vertical angle wrt
-viewpoint, in degrees, if the point is visible.
+<li> a value in [0, 180] representing the vertical angle with regard to
+the viewpoint, in degrees, if the point is visible. 0 degrees: directly
+above the observer, 180 degrees: directly below the observer
</ul>
@@ -75,6 +76,11 @@
<em>obs_elev=...</em>. The value entered is in the same units
as the elevation.
+<p>Target elevation: By default the target is assumed to have
+height=0 above the terrain. The user can change this using option
+<em>tgt_elev=...</em> to determine if objects of a given height would be
+visible. The value entered is in the same units as the elevation.
+
<p>Maximum visibility distance: By default there is no restriction on
the maximum distance to which the observer can see. The user can set
a maximum distance of visibility using option <em>max_dist=...</em>.
@@ -109,20 +115,16 @@
<h3>The algorithm:</h3>
<p>
<em>r.viewshed</em> uses the following model for determining
-visibility: Each point in the raster is assumed to represent the
-elevation of a flat cell centered at that point (the "circum-cell" of
-the grid point). The height of a cell is assumed to be constant,
-and thus the terrain is viewed as a tesselation of flat cells. Two
-points are visible to each other if their line-of-sight does not
-intersect the terrain. Since for an arbitrary point x in the terrain,
-the closest grid point to it is the one whose "circum-cell" contains
-x, this means that this model does a nearest-neighbor interpolation of
-heights.
+visibility: The height of a cell is assumed to be variable, and the
+actual height of a point falling into a cell, but not identical the cell
+center, is interpolated. Thus the terrain is viewed as a smooth surface.
+Two points are visible to each other if their line-of-sight does not
+intersect the terrain. The height for an arbitrary point x in the terrain
+is interpolated from the 4 surrounding neighbours. This means that this
+model does a bilinear interpolation of heights.
-This model is suitable for high resolution rasters; it may not be
-accurate for low resolution rasters, where it may be better to
-interpolate the height at a point based on several neighbors, rather
-than just nearest neighbor.
+This model is suitable for both low and high resolution rasters as well
+as terrain with flat and steep slopes.
<p>The core of the algorithm is determining, for each cell, the
line-of-sight and its intersections with the cells in the terrain. For
@@ -150,9 +152,7 @@
cell's gradient is higher, it is marked as visible, whereas if it is
lower, it is marked as invisible.
-<p>The advantage of considering the cells to be flat is that the
-algorithm can keep track and query efficiently the active events.
-For a (square) raster of <em>n</em> point in total, the standard
+<p>For a (square) raster of <em>n</em> point in total, the standard
viewshed algorithm uses <em>O(n sqrt(n))= O(n<sup>3/2</sup>)</em>
time, while the sweep-line algorithm uses <em>O(n lg n)</em> time.
This algorithm is efficient in terms of CPU operations and can be also
@@ -210,5 +210,7 @@
<p>William Richard (Bowdoin College): <tt>willster3021 at gmail.com </tt>
+<p>Markus Metz
+
<p><i>Last changed: $Date$</i>
Modified: grass-addons/raster/r.viewshed/rbbst.cc
===================================================================
--- grass-addons/raster/r.viewshed/rbbst.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/rbbst.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -74,13 +73,11 @@
#include <math.h>
#include "rbbst.h"
-#ifdef __GRASS__
extern "C"
{
-#include <grass/config.h>
#include <grass/gis.h>
+#include <grass/glocale.h>
}
-#endif
@@ -93,13 +90,8 @@
RBTree *create_tree(TreeValue tv)
{
init_nil_node();
-#ifdef __GRASS__
RBTree *rbt = (RBTree *) G_malloc(sizeof(RBTree));
TreeNode *root = (TreeNode *) G_malloc(sizeof(TreeNode));
-#else
- RBTree *rbt = (RBTree *) malloc(sizeof(RBTree));
- TreeNode *root = (TreeNode *) malloc(sizeof(TreeNode));
-#endif
rbt->root = root;
rbt->root->value = tv;
@@ -131,11 +123,7 @@
destroy_sub_tree(node->left);
destroy_sub_tree(node->right);
-#ifdef __GRASS__
G_free(node);
-#else
- free(node);
-#endif
return;
}
@@ -157,23 +145,25 @@
}
/*------------The following is designed for viewshed's algorithm-------*/
-double find_max_gradient_within_key(RBTree * rbt, double key)
+double find_max_gradient_within_key(RBTree * rbt, double key, double angle, double gradient)
{
- return find_max_value_within_key(rbt->root, key);
+ return find_max_value_within_key(rbt->root, key, angle, gradient);
}
/*<--------------------------------->
//Private below this line */
void init_nil_node()
{
-#ifdef __GRASS__
NIL = (TreeNode *) G_malloc(sizeof(TreeNode));
-#else
- NIL = (TreeNode *) malloc(sizeof(TreeNode));
-#endif
NIL->color = RB_BLACK;
- NIL->value.gradient = SMALLEST_GRADIENT;
+ NIL->value.angle[0] = 0;
+ NIL->value.angle[1] = 0;
+ NIL->value.angle[2] = 0;
+ NIL->value.gradient[0] = SMALLEST_GRADIENT;
+ NIL->value.gradient[1] = SMALLEST_GRADIENT;
+ NIL->value.gradient[2] = SMALLEST_GRADIENT;
NIL->value.maxGradient = SMALLEST_GRADIENT;
+ NIL->value.key = 0;
NIL->parent = NULL;
NIL->left = NULL;
@@ -188,9 +178,9 @@
//2: v1 > v2 */
char compare_values(TreeValue * v1, TreeValue * v2)
{
- if (v1->gradient > v2->gradient)
+ if (v1->gradient[1] > v2->gradient[1])
return 1;
- if (v1->gradient < v2->gradient)
+ if (v1->gradient[1] < v2->gradient[1])
return -1;
return 0;
@@ -215,11 +205,7 @@
{
TreeNode *ret;
-#ifdef __GRASS__
ret = (TreeNode *) G_malloc(sizeof(TreeNode));
-#else
- ret = (TreeNode *) malloc(sizeof(TreeNode));
-#endif
ret->color = RB_RED;
@@ -278,7 +264,7 @@
TreeNode *inserted = nextNode;
/*update augmented maxGradient */
- nextNode->value.maxGradient = nextNode->value.gradient;
+ nextNode->value.maxGradient = nextNode->value.gradient[1];
while (nextNode->parent != NIL) {
if (nextNode->parent->value.maxGradient < nextNode->value.maxGradient)
nextNode->parent->value.maxGradient = nextNode->value.maxGradient;
@@ -383,6 +369,8 @@
while (y != NIL && x == y->right) {
x = y;
+ if (y->parent == NIL)
+ return y;
y = y->parent;
}
return y;
@@ -401,10 +389,8 @@
z = search_for_node(*root, key);
if (z == NIL) {
- printf("ATTEMPT to delete key=%f failed\n", key);
- fprintf(stderr, "Node not found. Deletion fails.\n");
- exit(1);
- return; /*node to delete is not found */
+ /*node to delete is not found */
+ G_fatal_error(_("Attempt to delete node with key=%f failed"), key);
}
/*1-3 */
@@ -412,6 +398,10 @@
y = z;
else
y = tree_successor(z);
+
+ if (y == NIL) {
+ G_fatal_error(_("Successor node not found. Deletion fails."));
+ }
/*4-6 */
if (y->left != NIL)
@@ -442,7 +432,7 @@
double left, right;
while (curNode->parent != NIL) {
- if (curNode->parent->value.maxGradient == y->value.gradient) {
+ if (curNode->parent->value.maxGradient == y->value.gradient[1]) {
left = find_max_value(curNode->parent->left);
right = find_max_value(curNode->parent->right);
@@ -451,10 +441,10 @@
else
curNode->parent->value.maxGradient = right;
- if (curNode->parent->value.gradient >
+ if (curNode->parent->value.gradient[1] >
curNode->parent->value.maxGradient)
curNode->parent->value.maxGradient =
- curNode->parent->value.gradient;
+ curNode->parent->value.gradient[1];
}
else {
break;
@@ -468,17 +458,22 @@
toFix->left->value.maxGradient >
toFix->right->value.maxGradient ? toFix->left->value.
maxGradient : toFix->right->value.maxGradient;
- if (tmpMax > toFix->value.gradient)
+ if (tmpMax > toFix->value.gradient[1])
toFix->value.maxGradient = tmpMax;
else
- toFix->value.maxGradient = toFix->value.gradient;
+ toFix->value.maxGradient = toFix->value.gradient[1];
/*13-15 */
- if (y != z) {
- double zGradient = z->value.gradient;
+ if (y != NIL && y != z) {
+ double zGradient = z->value.gradient[1];
z->value.key = y->value.key;
- z->value.gradient = y->value.gradient;
+ z->value.gradient[0] = y->value.gradient[0];
+ z->value.gradient[1] = y->value.gradient[1];
+ z->value.gradient[2] = y->value.gradient[2];
+ z->value.angle[0] = y->value.angle[0];
+ z->value.angle[1] = y->value.angle[1];
+ z->value.angle[2] = y->value.angle[2];
toFix = z;
@@ -487,14 +482,14 @@
toFix->left->value.maxGradient >
toFix->right->value.maxGradient ? toFix->left->value.
maxGradient : toFix->right->value.maxGradient;
- if (tmpMax > toFix->value.gradient)
+ if (tmpMax > toFix->value.gradient[1])
toFix->value.maxGradient = tmpMax;
else
- toFix->value.maxGradient = toFix->value.gradient;
+ toFix->value.maxGradient = toFix->value.gradient[1];
while (z->parent != NIL) {
if (z->parent->value.maxGradient == zGradient) {
- if (z->parent->value.gradient != zGradient &&
+ if (z->parent->value.gradient[1] != zGradient &&
(!(z->parent->left->value.maxGradient == zGradient &&
z->parent->right->value.maxGradient == zGradient))) {
@@ -506,10 +501,10 @@
else
z->parent->value.maxGradient = right;
- if (z->parent->value.gradient >
+ if (z->parent->value.gradient[1] >
z->parent->value.maxGradient)
z->parent->value.maxGradient =
- z->parent->value.gradient;
+ z->parent->value.gradient[1];
}
@@ -528,11 +523,7 @@
rb_delete_fixup(root, x);
/*18 */
-#ifdef __GRASS__
G_free(y);
-#else
- free(y);
-#endif
return;
}
@@ -633,7 +624,7 @@
/* find max within the max key */
-double find_max_value_within_key(TreeNode * root, double maxKey)
+double find_max_value_within_key(TreeNode * root, double maxKey, double angle, double gradient)
{
TreeNode *keyNode = search_for_node(root, maxKey);
@@ -644,16 +635,79 @@
exit(1);
}
- double max = find_max_value(keyNode->left);
+ TreeNode *currNode = keyNode;
+ double max = SMALLEST_GRADIENT;
double tmpMax;
+ double curr_gradient;
+ /* traverse all nodes with smaller distance */
+ while (currNode != NIL) {
+ int checkme = (currNode->value.angle[0] <= angle &&
+ currNode->value.angle[2] >= angle);
+
+ if (!checkme && currNode->value.key > 0) {
+ G_warning(_("\nangles outside angle %.4f"), angle);
+ G_warning(_("ENTER angle %.4f"), currNode->value.angle[0]);
+ G_warning(_("CENTER angle %.4f"), currNode->value.angle[1]);
+ G_warning(_("EXIT angle %.4f"), currNode->value.angle[2]);
+ G_warning(_("\nENTER gradient %.4f"), currNode->value.gradient[0]);
+ G_warning(_("CENTER gradient %.4f"), currNode->value.gradient[1]);
+ G_warning(_("EXIT gradient %.4f"), currNode->value.gradient[2]);
+ }
+
+ if (currNode->value.key > maxKey) {
+ G_fatal_error(_("current dist too large %.4f > %.4f"), currNode->value.key, maxKey);
+ }
+
+
+ if (checkme && currNode != keyNode) {
+ if (angle < currNode->value.angle[1]) {
+ curr_gradient = currNode->value.gradient[1] +
+ (currNode->value.gradient[0] - currNode->value.gradient[1]) *
+ (currNode->value.angle[1] - angle) /
+ (currNode->value.angle[1] - currNode->value.angle[0]);
+ }
+ else if (angle > currNode->value.angle[1]) {
+ curr_gradient = currNode->value.gradient[1] +
+ (currNode->value.gradient[2] - currNode->value.gradient[1]) *
+ (angle - currNode->value.angle[1]) /
+ (currNode->value.angle[2] - currNode->value.angle[1]);
+ }
+ else /* angle == currNode->value.angle[1] */
+ curr_gradient = currNode->value.gradient[1];
+
+ if (curr_gradient > max)
+ max = curr_gradient;
+
+ if (max > gradient)
+ return max;
+ }
+ /* get next smaller key */
+ if (currNode->left != NIL) {
+ currNode = currNode->left;
+ while (currNode->right != NIL)
+ currNode = currNode->right;
+ }
+ else {
+ /* at smallest item in this branch, go back up */
+ TreeNode *lastNode;
+
+ do {
+ lastNode = currNode;
+ currNode = currNode->parent;
+ } while (currNode != NIL && lastNode == currNode->left);
+ }
+ }
+ return max;
+
+ /* old code assuming flat cells */
while (keyNode->parent != NIL) {
if (keyNode == keyNode->parent->right) { /*its the right node of its parent; */
tmpMax = find_max_value(keyNode->parent->left);
if (tmpMax > max)
max = tmpMax;
- if (keyNode->parent->value.gradient > max)
- max = keyNode->parent->value.gradient;
+ if (keyNode->parent->value.gradient[1] > max)
+ max = keyNode->parent->value.gradient[1];
}
keyNode = keyNode->parent;
}
@@ -675,20 +729,20 @@
tmpMax = x->left->value.maxGradient > y->left->value.maxGradient ?
x->left->value.maxGradient : y->left->value.maxGradient;
- if (tmpMax > x->value.gradient)
+ if (tmpMax > x->value.gradient[1])
x->value.maxGradient = tmpMax;
else
- x->value.maxGradient = x->value.gradient;
+ x->value.maxGradient = x->value.gradient[1];
/*fix y */
tmpMax = x->value.maxGradient > y->right->value.maxGradient ?
x->value.maxGradient : y->right->value.maxGradient;
- if (tmpMax > y->value.gradient)
+ if (tmpMax > y->value.gradient[1])
y->value.maxGradient = tmpMax;
else
- y->value.maxGradient = y->value.gradient;
+ y->value.maxGradient = y->value.gradient[1];
/*left rotation
//see pseudocode on page 278 in CLRS */
@@ -727,19 +781,19 @@
tmpMax = x->right->value.maxGradient > y->right->value.maxGradient ?
x->right->value.maxGradient : y->right->value.maxGradient;
- if (tmpMax > y->value.gradient)
+ if (tmpMax > y->value.gradient[1])
y->value.maxGradient = tmpMax;
else
- y->value.maxGradient = y->value.gradient;
+ y->value.maxGradient = y->value.gradient[1];
/*fix x */
tmpMax = x->left->value.maxGradient > y->value.maxGradient ?
x->left->value.maxGradient : y->value.maxGradient;
- if (tmpMax > x->value.gradient)
+ if (tmpMax > x->value.gradient[1])
x->value.maxGradient = tmpMax;
else
- x->value.maxGradient = x->value.gradient;
+ x->value.maxGradient = x->value.gradient[1];
/*ratation */
y->left = x->right;
Modified: grass-addons/raster/r.viewshed/rbbst.h
===================================================================
--- grass-addons/raster/r.viewshed/rbbst.h 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/rbbst.h 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -59,7 +58,8 @@
double key;
/*anything else below this line is optional */
- double gradient;
+ double gradient[3];
+ double angle[3];
double maxGradient;
} TreeValue;
@@ -93,7 +93,7 @@
/*------------The following is designed for kreveld's algorithm-------*/
-double find_max_gradient_within_key(RBTree * rbt, double key);
+double find_max_gradient_within_key(RBTree * rbt, double key, double angle, double gradient);
/*LT: not sure if this is correct */
int is_empty(RBTree * t);
@@ -152,6 +152,6 @@
double find_max_value(TreeNode * root);
/*find max within the max key */
-double find_max_value_within_key(TreeNode * root, double maxKey);
+double find_max_value_within_key(TreeNode * root, double maxKey, double angle, double gradient);
#endif
Modified: grass-addons/raster/r.viewshed/statusstructure.cc
===================================================================
--- grass-addons/raster/r.viewshed/statusstructure.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/statusstructure.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -42,15 +41,12 @@
#include <stdlib.h>
#include <assert.h>
-#ifdef __GRASS__
extern "C"
{
-#include <grass/config.h>
#include <grass/gis.h>
#include <grass/glocale.h>
}
#include "grass.h"
-#endif
#include "statusstructure.h"
@@ -62,23 +58,31 @@
/*find the vertical angle in degrees between the viewpoint and the
point represented by the StatusNode. Assumes all values (except
gradient) in sn have been filled. The value returned is in [0,
- 180]. if doCurv is set we need to consider the curvature of the
+ 180]. A value of 0 is directly below the specified viewing position,
+ 90 is due horizontal, and 180 is directly above the observer.
+ If doCurv is set we need to consider the curvature of the
earth */
-float get_vertical_angle(Viewpoint vp, StatusNode sn, int doCurv)
+float get_vertical_angle(Viewpoint vp, StatusNode sn, surface_type elev, int doCurv)
{
/*determine the difference in elevation, based on the curvature */
- float diffElev;
+ double diffElev;
+ diffElev = vp.elev - elev;
- diffElev = vp.elev - sn.elev;
-
/*calculate and return the angle in degrees */
assert(fabs(sn.dist2vp) > 0.001);
+ /* 0 above, 180 below */
if (diffElev >= 0.0)
return (atan(sqrt(sn.dist2vp) / diffElev) * (180 / M_PI));
else
- return ((atan(fabs(diffElev) / sqrt(sn.dist2vp)) * (180 / M_PI)) + 90);
+ return (atan(fabs(diffElev) / sqrt(sn.dist2vp)) * (180 / M_PI) + 90);
+
+ /* 180 above, 0 below */
+ if (diffElev >= 0.0)
+ return (atan(diffElev / sqrt(sn.dist2vp)) * (180 / M_PI) + 90);
+ else
+ return (atan(sqrt(sn.dist2vp) / fabs(diffElev)) * (180 / M_PI));
}
@@ -90,21 +94,20 @@
long long sizeBytes;
- printf("Estimated size active structure:");
- printf(" (key=%d, ptr=%d, total node=%d B)",
+ G_verbose_message(_("Estimated size active structure:"));
+ G_verbose_message(_(" (key=%d, ptr=%d, total node=%d B)"),
(int)sizeof(TreeValue),
(int)sizeof(TreeNode *), (int)sizeof(TreeNode));
sizeBytes = sizeof(TreeNode) * max(hd->ncols, hd->nrows);
- printf(" Total= %lld B\n", sizeBytes);
+ G_verbose_message(_(" Total= %lld B"), sizeBytes);
return sizeBytes;
}
-
-
/* ------------------------------------------------------------ */
/*given a StatusNode, fill in its dist2vp and gradient */
-void calculate_dist_n_gradient(StatusNode * sn, Viewpoint * vp)
+void calculate_dist_n_gradient(StatusNode * sn, double elev,
+ Viewpoint * vp, GridHeader hd)
{
assert(sn && vp);
/*sqrt is expensive
@@ -112,44 +115,98 @@
// pow(sn->col - vp->col,2.0)));
//sn->gradient = (sn->elev - vp->elev)/(sn->dist2vp); */
- double elevDiff = (double)sn->elev - vp->elev;
+ double diffElev = elev - vp->elev;
+ double dx = ((double)sn->col - vp->col) * hd.ew_res;
+ double dy = ((double)sn->row - vp->row) * hd.ns_res;
+
+ sn->dist2vp = (dx * dx) + (dy * dy);
- sn->dist2vp = (sn->row - vp->row) * (sn->row - vp->row) +
- (sn->col - vp->col) * (sn->col - vp->col);
- sn->gradient = elevDiff * elevDiff / (sn->dist2vp);
- sn->gradient = atan(sqrt(sn->gradient));
- /*maintain sign */
- if (sn->elev < vp->elev)
- sn->gradient = -sn->gradient;
+ if (diffElev == 0) {
+ sn->gradient[1] = 0;
+ return;
+ }
- elevDiff += vp->target_offset;
- sn->gradient_offset = elevDiff * elevDiff / (sn->dist2vp);
- sn->gradient_offset = atan(sqrt(sn->gradient_offset));
+ /* PI / 2 above, - PI / 2 below, like r.los */
+ sn->gradient[1] = atan(diffElev / sqrt(sn->dist2vp));
+
+ return;
+
+ /* PI above, 0 below. slower than r.los - like */
+ if (diffElev >= 0.0)
+ sn->gradient[1] = (atan(diffElev / sqrt(sn->dist2vp)) + M_PI / 2);
+ else
+ sn->gradient[1] = (atan(sqrt(sn->dist2vp) / fabs(diffElev)));
+
+ return;
+
+ /* a little bit faster but not accurate enough */
+ sn->gradient[1] = (diffElev * diffElev) / (sn->dist2vp);
/*maintain sign */
- if (sn->elev + vp->target_offset < vp->elev)
- sn->gradient_offset = -sn->gradient_offset;
+ if (elev < vp->elev)
+ sn->gradient[1] = -sn->gradient[1];
+
return;
}
+/* ------------------------------------------------------------ */
+/* calculate gradient for ENTERING or EXITING event */
+void calculate_event_gradient(StatusNode * sn, int e_idx,
+ double row, double col, double elev,
+ Viewpoint * vp, GridHeader hd)
+{
+ assert(sn && vp);
+ /*sqrt is expensive
+ //sn->dist2vp = sqrt((float) ( pow(sn->row - vp->row,2.0) +
+ // pow(sn->col - vp->col,2.0)));
+ //sn->gradient = (sn->elev - vp->elev)/(sn->dist2vp); */
+
+ double diffElev = elev - vp->elev;
+ double dx = (col - vp->col) * hd.ew_res;
+ double dy = (row - vp->row) * hd.ns_res;
+ double dist2vp = (dx * dx) + (dy * dy);
+ /* PI / 2 above, - PI / 2 below */
+ sn->gradient[e_idx] = atan(diffElev / sqrt(dist2vp));
+
+ return;
+
+ /* PI above, 0 below. slower than r.los - like */
+ if (diffElev >= 0.0)
+ sn->gradient[e_idx] = (atan(diffElev / sqrt(dist2vp)) + M_PI / 2);
+ else
+ sn->gradient[e_idx] = (atan(sqrt(dist2vp) / fabs(diffElev)));
+
+ return;
+
+ /* faster but not accurate enough */
+ sn->gradient[e_idx] = (diffElev * diffElev) / (dist2vp);
+ /*maintain sign */
+ if (elev < vp->elev)
+ sn->gradient[e_idx] = -sn->gradient[e_idx];
+
+ return;
+}
+
+
/* ------------------------------------------------------------ */
/*create an empty status list */
StatusList *create_status_struct()
{
StatusList *sl;
-#ifdef __GRASS__
sl = (StatusList *) G_malloc(sizeof(StatusList));
-#else
- sl = (StatusList *) malloc(sizeof(StatusList));
-#endif
assert(sl);
TreeValue tv;
- tv.gradient = SMALLEST_GRADIENT;
+ tv.gradient[0] = SMALLEST_GRADIENT;
+ tv.gradient[1] = SMALLEST_GRADIENT;
+ tv.gradient[2] = SMALLEST_GRADIENT;
+ tv.angle[0] = 0;
+ tv.angle[1] = 0;
+ tv.angle[2] = 0;
tv.key = 0;
tv.maxGradient = SMALLEST_GRADIENT;
@@ -165,11 +222,8 @@
{
assert(sl);
delete_tree(sl->rbt);
-#ifdef __GRASS__
G_free(sl);
-#else
- free(sl);
-#endif
+
return;
}
@@ -194,9 +248,15 @@
TreeValue tv;
tv.key = sn.dist2vp;
- tv.gradient = sn.gradient;
+ tv.gradient[0] = sn.gradient[0];
+ tv.gradient[1] = sn.gradient[1];
+ tv.gradient[2] = sn.gradient[2];
+ tv.angle[0] = sn.angle[0];
+ tv.angle[1] = sn.angle[1];
+ tv.angle[2] = sn.angle[2];
tv.maxGradient = SMALLEST_GRADIENT;
insert_into(sl->rbt, tv);
+
return;
}
@@ -204,7 +264,7 @@
/* ------------------------------------------------------------ */
/*find the node with max Gradient within the distance (from viewpoint)
//given */
-double find_max_gradient_in_status_struct(StatusList * sl, double dist)
+double find_max_gradient_in_status_struct(StatusList * sl, double dist, double angle, double gradient)
{
assert(sl);
/*note: if there is nothing in the status struccture, it means this
@@ -214,10 +274,10 @@
/*it is also possible that the status structure is not empty, but
there are no events with key < dist ---in this case it returns
SMALLEST_GRADIENT; */
- return find_max_gradient_within_key(sl->rbt, dist);
+ return find_max_gradient_within_key(sl->rbt, dist, angle, gradient);
}
-/*returns true is it is empty */
+/*returns true if it is empty */
int is_empty(StatusList * sl)
{
assert(sl);
Modified: grass-addons/raster/r.viewshed/statusstructure.h
===================================================================
--- grass-addons/raster/r.viewshed/statusstructure.h 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/statusstructure.h 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -44,9 +43,7 @@
This header file defines the status structure and related functions.
*/
-#ifdef __GRASS__
-#include <grass/config.h>
-#endif
+#include <grass/gis.h>
#include "grid.h"
#include "rbbst.h"
@@ -55,10 +52,12 @@
typedef struct statusnode_
{
dimensionType row, col; /*position of the cell */
- float elev; /*elevation of cell */
+
+ /* float elev; */ /*elevation of cell */
double dist2vp; /*distance to the viewpoint */
- double gradient; /*gradient of the Line of Sight */
- double gradient_offset; /*gradient of the Line of Sight with local elevation offset */
+ double gradient[3]; /*ENTER, CENTER, EXIT gradients of the Line of Sight */
+ double angle[3]; /*ENTER, CENTER, EXIT angles of the Line of Sight */
+ /* double gradient_offset; */ /*gradient of the Line of Sight with local elevation offset */
} StatusNode;
@@ -75,9 +74,14 @@
long long get_active_str_size_bytes(GridHeader * hd);
-//given a StatusNode, fill in its dist2vp and gradient
-void calculate_dist_n_gradient(StatusNode * sn, Viewpoint * vp);
+/*given a StatusNode, fill in its dist2vp and gradient */
+void calculate_dist_n_gradient(StatusNode * sn, double elev,
+ Viewpoint * vp, GridHeader hd);
+/* calculate gradient for ENTERING or EXITING event */
+void calculate_event_gradient(StatusNode * sn, int e_idx,
+ double row, double col, double elev,
+ Viewpoint * vp, GridHeader hd);
/*create an empty status list. */
StatusList *create_status_struct();
@@ -96,12 +100,12 @@
/*find the node with max Gradient. The node must be
//within the distance (from viewpoint) given */
-double find_max_gradient_in_status_struct(StatusList * sl, double dist);
+double find_max_gradient_in_status_struct(StatusList * sl, double dist, double angle, double gradient);
/*find the vertical angle in degrees between the viewpoint and the
point represented by the StatusNode. Assumes all values (except
gradient) in sn have been filled. */
-float get_vertical_angle(Viewpoint vp, StatusNode sn, int doCurv);
+float get_vertical_angle(Viewpoint vp, StatusNode sn, surface_type elev, int doCurv);
#endif
Modified: grass-addons/raster/r.viewshed/viewshed.cc
===================================================================
--- grass-addons/raster/r.viewshed/viewshed.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/viewshed.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -5,11 +5,12 @@
*
* AUTHOR(S): Laura Toma, Bowdoin College - ltoma at bowdoin.edu
* Yi Zhuang - yzhuang at bowdoin.edu
-
+ *
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: oct 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -41,13 +40,17 @@
#include <stdlib.h>
#include <stdio.h>
+extern "C"
+{
+#include "grass/gis.h"
+#include "grass/glocale.h"
+}
+
#include "viewshed.h"
#include "visibility.h"
#include "eventlist.h"
#include "statusstructure.h"
-#ifdef __GRASS
#include "grass.h"
-#endif
#define VIEWSHEDDEBUG if(0)
@@ -65,27 +68,25 @@
/* the output visibility grid */
long long totalcells = (long long)hd->nrows * (long long)hd->ncols;
- printf("rows=%d, cols=%d, total = %lld\n", hd->nrows, hd->ncols,
+ G_verbose_message(_("rows=%d, cols=%d, total = %lld"), hd->nrows, hd->ncols,
totalcells);
long long gridMemUsage = totalcells * sizeof(float);
- VIEWSHEDDEBUG {
- printf("grid usage=%lld\n", gridMemUsage);
- }
+ G_debug(1, "grid usage=%lld", gridMemUsage);
/* the event array */
long long eventListMemUsage = totalcells * 3 * sizeof(AEvent);
- VIEWSHEDDEBUG {
- printf("memory_usage: eventList=%lld\n", eventListMemUsage);
- }
+ G_debug(1, "memory_usage: eventList=%lld", eventListMemUsage);
/* the double array <data> that stores all the cells in the same row
as the viewpoint */
long long dataMemUsage = (long long)(hd->ncols * sizeof(double));
- printf("viewshed memory usage: size AEvent=%dB, nevents=%lld, \
- total=%lld B (%d MB)\n", (int)sizeof(AEvent), totalcells * 3, gridMemUsage + eventListMemUsage + dataMemUsage, (int)((gridMemUsage + eventListMemUsage + dataMemUsage) >> 20));
+ G_debug(1, "viewshed memory usage: size AEvent=%dB, nevents=%lld, \
+ total=%lld B (%d MB)", (int)sizeof(AEvent), totalcells * 3,
+ gridMemUsage + eventListMemUsage + dataMemUsage,
+ (int)((gridMemUsage + eventListMemUsage + dataMemUsage) >> 20));
return (gridMemUsage + eventListMemUsage + dataMemUsage);
@@ -100,22 +101,16 @@
char timeused[1000];
- printf("sweep timings:\n");
+ G_verbose_message(_("sweep timings:"));
rt_sprint_safe(timeused, initEventTime);
- printf("\t%30s", "init events: ");
- printf(timeused);
- printf("\n");
+ G_verbose_message("init events: %s", timeused);
rt_sprint_safe(timeused, sortEventTime);
- printf("\t%30s", "sort events: ");
- printf(timeused);
- printf("\n");
+ G_verbose_message("sort events: %s", timeused);
rt_sprint_safe(timeused, sweepTime);
- printf("\t%30s", "process events: ");
- printf(timeused);
- printf("\n");
- fflush(stdout);
+ G_verbose_message("process events: %s", timeused);
+
return;
}
@@ -123,8 +118,8 @@
/* ------------------------------------------------------------ */
static void print_statusnode(StatusNode sn)
{
- printf("processing (row=%d, col=%d, elev=%f, dist=%f, grad=%f)",
- sn.row, sn.col, sn.elev, sn.dist2vp, sn.gradient);
+ G_debug(3, "processing (row=%d, col=%d, dist=%f, grad=%f)",
+ sn.row, sn.col, sn.dist2vp, sn.gradient[1]);
return;
}
@@ -143,40 +138,37 @@
long long totalsize = hd->ncols * hd->nrows * 3;
totalsize *= sizeof(AEvent);
- printf("total size of eventlist is %lld B (%d MB); ",
+ G_debug(1, "total size of eventlist is %lld B (%d MB); ",
totalsize, (int)(totalsize >> 20));
/* what's the size of size_t on this machine? */
int sizet_size;
sizet_size = (int)sizeof(size_t);
- printf("size_t is %lu B\n", sizeof(size_t));
+ G_debug(1, "size_t is %d B", sizet_size);
if (sizet_size >= 8) {
- printf("64-bit platform, great.\n");
+ G_debug(1, "64-bit platform, great.");
}
else {
/* this is the max value of size_t */
long long maxsizet = ((long long)1 << (sizeof(size_t) * 8)) - 1;
- printf("max size_t is %lld\n", maxsizet);
+ G_debug(1, "max size_t is %lld", maxsizet);
/* checking whether allocating totalsize causes an overflow */
if (totalsize > maxsizet) {
- printf
- ("running the program in-memory mode requires memory beyond the capability of the platform. Use external mode, or 64-bit platform.\n");
- exit(1);
+ G_fatal_error(_("Running the program in-memory mode requires " \
+ "memory beyond the capability of the platform. " \
+ "Use external mode, or a 64-bit platform."));
}
}
- printf("allocating..");
-#ifdef __GRASS__
+ G_debug(1, "allocating eventList...");
eventList = (AEvent *) G_malloc(totalsize);
-#else
- eventList = (AEvent *) malloc(totalsize * sizeof(char));
-#endif
+
assert(eventList);
- printf("..ok\n");
+ G_debug(1, "...ok");
return eventList;
}
@@ -206,8 +198,7 @@
{
assert(inputfname && hd && vp);
- printf("Start sweeping.\n");
- fflush(stdout);
+ G_verbose_message(_("Start sweeping."));
/* ------------------------------ */
/* create the visibility grid */
@@ -217,20 +208,17 @@
/* set everything initially invisible */
set_inmem_visibilitygrid(visgrid, INVISIBLE);
assert(visgrid);
- VIEWSHEDDEBUG {
- printf("visibility grid size: %d x %d x %d B (%d MB)\n",
+ G_debug(1, "visibility grid size: %d x %d x %d B (%d MB)",
hd->nrows, hd->ncols, (int)sizeof(float),
(int)(((long long)(hd->nrows * hd->ncols *
sizeof(float))) >> 20));
- }
-
/* ------------------------------ */
/* construct the event list corresponding to the given input file
and viewpoint; this creates an array of all the cells on the
same row as the viewpoint */
- double *data;
+ surface_type **data;
size_t nevents;
Rtimer initEventTime;
@@ -239,23 +227,19 @@
AEvent *eventList = allocate_eventlist(hd);
-#ifdef __GRASS__
- nevents = grass_init_event_list_in_memory(eventList, inputfname, vp, hd,
+ nevents = init_event_list_in_memory(eventList, inputfname, vp, hd,
viewOptions, &data, visgrid);
-#else
- nevents = init_event_list_in_memory(eventList, inputfname, vp, hd,
- viewOptions, &data, visgrid);
-#endif
+
assert(data);
rt_stop(initEventTime);
- printf("actual nb events is %lu\n", nevents);
+ G_debug(1, "actual nb events is %lu", (long unsigned int)nevents);
/* ------------------------------ */
/*sort the events radially by angle */
Rtimer sortEventTime;
rt_start(sortEventTime);
- printf("sorting events..");
+ G_verbose_message(_("sorting events..."));
fflush(stdout);
/*this is recursive and seg faults for large arrays
@@ -268,13 +252,11 @@
RadialCompare cmpObj;
quicksort(eventList, nevents, cmpObj);
- printf("done\n");
+ G_verbose_message(_("done"));
fflush(stdout);
rt_stop(sortEventTime);
-
-
/* ------------------------------ */
/*create the status structure */
StatusList *status_struct = create_status_struct();
@@ -285,28 +267,52 @@
rt_start(sweepTime);
for (dimensionType i = vp->col + 1; i < hd->ncols; i++) {
+ AEvent e;
+ double ax, ay;
+
sn.col = i;
sn.row = vp->row;
- sn.elev = data[i];
- if (!is_nodata(visgrid->grid->hd, sn.elev) &&
+ e.col = i;
+ e.row = vp->row;
+ e.elev[0] = data[0][i];
+ e.elev[1] = data[1][i];
+ e.elev[2] = data[2][i];
+
+ if (!is_nodata(visgrid->grid->hd, data[1][i]) &&
!is_point_outside_max_dist(*vp, *hd, sn.row, sn.col,
viewOptions.maxDist)) {
/*calculate Distance to VP and Gradient, store them into sn */
- calculate_dist_n_gradient(&sn, vp);
- VIEWSHEDDEBUG {
- printf("inserting: ");
- print_statusnode(sn);
- printf("\n");
- }
+ /* need either 3 elevation values or
+ * 3 gradients calculated from 3 elevation values */
+ /* need also 3 angles */
+ e.eventType = ENTERING_EVENT;
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_event_gradient(&sn, 0, ay, ax, e.elev[0], vp, *hd);
+
+ e.eventType = CENTER_EVENT;
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_dist_n_gradient(&sn, e.elev[1], vp, *hd);
+
+ e.eventType = EXITING_EVENT;
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_event_gradient(&sn, 2, ay, ax, e.elev[2], vp, *hd);
+
+ assert(sn.angle[1] == 0);
+
+ if (sn.angle[0] > sn.angle[1])
+ sn.angle[0] -= 2 * M_PI;
+
+ G_debug(2, "inserting: ");
+ print_statusnode(sn);
/*insert sn into the status structure */
insert_into_status_struct(sn, status_struct);
}
}
-#ifdef __GRASS__
+ G_free(data[0]);
G_free(data);
-#else
- free(data);
-#endif
@@ -315,66 +321,93 @@
long nvis = 0; /*number of visible cells */
AEvent *e;
-#ifdef __GRASS__
- G_message("Determine visibility...");
+ G_message(_("Determine visibility..."));
G_percent(0, 100, 2);
-#endif
+
for (size_t i = 0; i < nevents; i++) {
-#ifdef __GRASS__
- int perc = (int)((double)i / nevents * 1000000000.);
- if (perc > 0)
- G_percent(perc, 1000000000, 2);
-#endif
+ int perc = (int)((double)i / nevents * 1000000.);
+ if (perc > 0 && perc < 1000000)
+ G_percent(perc, 1000000, 2);
+
/*get out one event at a time and process it according to its type */
e = &(eventList[i]);
sn.col = e->col;
sn.row = e->row;
- sn.elev = e->elev;
+ //sn.elev = e->elev;
/*calculate Distance to VP and Gradient */
- calculate_dist_n_gradient(&sn, vp);
- INMEMORY_DEBUG {
- printf("event: ");
- print_event(*e);
- printf("sn.dist=%f, sn.gradient=%f\n", sn.dist2vp, sn.gradient);
- }
+ calculate_dist_n_gradient(&sn, e->elev[1] + vp->target_offset, vp, *hd);
+ G_debug(3, "event: ");
+ print_event(*e, 3);
+ G_debug(3, "sn.dist=%f, sn.gradient=%f", sn.dist2vp, sn.gradient[1]);
switch (e->eventType) {
case ENTERING_EVENT:
+ double ax, ay;
/*insert node into structure */
- VIEWSHEDDEBUG {
- printf("..ENTER-EVENT: insert\n");
+ G_debug(3, "..ENTER-EVENT: insert");
+
+ /* need either 3 elevation values or
+ * 3 gradients calculated from 3 elevation values */
+ /* need also 3 angles */
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ //sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row);
+ sn.angle[0] = e->angle;
+ calculate_event_gradient(&sn, 0, ay, ax, e->elev[0], vp, *hd);
+
+ e->eventType = CENTER_EVENT;
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_dist_n_gradient(&sn, e->elev[1], vp, *hd);
+
+ e->eventType = EXITING_EVENT;
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_event_gradient(&sn, 2, ay, ax, e->elev[2], vp, *hd);
+
+ e->eventType = ENTERING_EVENT;
+
+ if (e->angle < M_PI) {
+ if (sn.angle[0] > sn.angle[1])
+ sn.angle[0] -= 2 * M_PI;
}
+ else {
+ if (sn.angle[0] > sn.angle[1]) {
+ sn.angle[1] += 2 * M_PI;
+ sn.angle[2] += 2 * M_PI;
+ }
+ }
+
insert_into_status_struct(sn, status_struct);
break;
case EXITING_EVENT:
/*delete node out of status structure */
- VIEWSHEDDEBUG {
- printf("..EXIT-EVENT: delete\n");
- }
+ G_debug(3, "..EXIT-EVENT: delete");
+ /* need only distance */
delete_from_status_struct(status_struct, sn.dist2vp);
break;
case CENTER_EVENT:
- VIEWSHEDDEBUG {
- printf("..QUERY-EVENT: query\n");
- }
+ G_debug(3, "..QUERY-EVENT: query");
/*calculate visibility */
double max;
+ /* consider current angle and gradient */
max =
- find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
+ find_max_gradient_in_status_struct(status_struct, sn.dist2vp,
+ e->angle, sn.gradient[1]);
/*the point is visible: store its vertical angle */
- if (max <= sn.gradient_offset) {
+ if (max <= sn.gradient[1]) {
+ float vert_angle = get_vertical_angle(*vp, sn, e->elev[1] + vp->target_offset,
+ viewOptions.doCurv);
+
add_result_to_inmem_visibilitygrid(visgrid, sn.row, sn.col,
- get_vertical_angle(*vp, sn,
- viewOptions.
- doCurv));
- assert(get_vertical_angle(*vp, sn, viewOptions.doCurv) >= 0);
+ vert_angle);
+ assert(vert_angle >= 0);
/* when you write the visibility grid you assume that
visible values are positive */
nvis++;
@@ -388,9 +421,10 @@
}
}
rt_stop(sweepTime);
+ G_percent(1, 1, 1);
- printf("Sweeping done.\n");
- printf("Total cells %ld, visible cells %ld (%.1f percent).\n",
+ G_verbose_message(_("Sweeping done."));
+ G_verbose_message(_("Total cells %ld, visible cells %ld (%.1f percent)."),
(long)visgrid->grid->hd->nrows * visgrid->grid->hd->ncols,
nvis,
(float)((float)nvis * 100 /
@@ -400,14 +434,9 @@
print_viewshed_timings(initEventTime, sortEventTime, sweepTime);
/*cleanup */
-#ifdef __GRASS__
G_free(eventList);
-#else
- free(eventList);
-#endif
return visgrid;
-
}
@@ -430,8 +459,7 @@
{
assert(inputfname && hd && vp);
- printf("Start sweeping.\n");
- fflush(stdout);
+ G_message(_("Start sweeping."));
/* ------------------------------ */
@@ -448,16 +476,13 @@
Rtimer initEventTime, sortEventTime, sweepTime;
AMI_STREAM < AEvent > *eventList;
- double *data;
+ surface_type **data;
rt_start(initEventTime);
-#ifdef __GRASS__
- eventList = grass_init_event_list(inputfname, vp, hd, viewOptions,
+
+ eventList = init_event_list(inputfname, vp, hd, viewOptions,
&data, visgrid);
-#else
- eventList =
- init_event_list(inputfname, vp, hd, viewOptions, &data, visgrid);
-#endif
+
assert(eventList && data);
eventList->seek(0);
rt_stop(initEventTime);
@@ -466,9 +491,8 @@
/* ------------------------------ */
/*sort the events radially by angle */
-#ifdef __GRASS__
- G_message("Sort the events...");
-#endif
+ G_verbose_message(_("Sorting events..."));
+
rt_start(sortEventTime);
sort_event_list(&eventList);
eventList->seek(0); /*this does not seem to be ensured by sort?? */
@@ -483,38 +507,59 @@
structure */
StatusNode sn;
-#ifdef __GRASS__
- G_message("Calculate distances...");
-#endif
+ G_message(_("Calculating distances..."));
+
rt_start(sweepTime);
for (dimensionType i = vp->col + 1; i < hd->ncols; i++) {
+ AEvent e;
+ double ax, ay;
-#ifdef __GRASS__
G_percent(i, hd->ncols, 2);
-#endif
+
sn.col = i;
sn.row = vp->row;
- sn.elev = data[i];
- if (!is_nodata(visgrid->hd, sn.elev) &&
+ e.col = i;
+ e.row = vp->row;
+ e.elev[0] = data[0][i];
+ e.elev[1] = data[1][i];
+ e.elev[2] = data[2][i];
+ if (!is_nodata(visgrid->hd, data[1][i]) &&
!is_point_outside_max_dist(*vp, *hd, sn.row, sn.col,
viewOptions.maxDist)) {
/*calculate Distance to VP and Gradient, store them into sn */
- calculate_dist_n_gradient(&sn, vp);
- VIEWSHEDDEBUG {
- printf("inserting: ");
- print_statusnode(sn);
- printf("\n");
- }
+ /* need either 3 elevation values or
+ * 3 gradients calculated from 3 elevation values */
+ /* need also 3 angles */
+ e.eventType = ENTERING_EVENT;
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_event_gradient(&sn, 0, ay, ax, e.elev[0], vp, *hd);
+
+ e.eventType = CENTER_EVENT;
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_dist_n_gradient(&sn, e.elev[1], vp, *hd);
+
+ e.eventType = EXITING_EVENT;
+ calculate_event_position(e, vp->row, vp->col, &ay, &ax);
+ sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_event_gradient(&sn, 2, ay, ax, e.elev[2], vp, *hd);
+
+ assert(sn.angle[1] == 0);
+
+ if (sn.angle[0] > sn.angle[1])
+ sn.angle[0] -= 2 * M_PI;
+
+ G_debug(3, "inserting: ");
+ print_statusnode(sn);
+
/*insert sn into the status structure */
insert_into_status_struct(sn, status_struct);
}
}
-#ifdef __GRASS__
G_percent(hd->ncols, hd->ncols, 2);
+ G_free(data[0]);
G_free(data);
-#else
- free(data);
-#endif
/* ------------------------------ */
@@ -527,64 +572,93 @@
/*printf("nbEvents = %ld\n", (long) nbEvents); */
-#ifdef __GRASS__
- G_message("Determine visibility...");
+ G_message(_("Determine visibility..."));
G_percent(0, 100, 2);
-#endif
+
for (off_t i = 0; i < nbEvents; i++) {
-
-#ifdef __GRASS__
- int perc = (int)((double)i / nbEvents * 1000000000.);
+
+ int perc = (int)((double)i / nbEvents * 1000000.);
if (perc > 0)
- G_percent(perc, 1000000000, 2);
-#endif
+ G_percent(perc, 1000000, 2);
+
/*get out one event at a time and process it according to its type */
ae = eventList->read_item(&e);
assert(ae == AMI_ERROR_NO_ERROR);
sn.col = e->col;
sn.row = e->row;
- sn.elev = e->elev;
+ //sn.elev = e->elev;
/*calculate Distance to VP and Gradient */
- calculate_dist_n_gradient(&sn, vp);
- VIEWSHEDDEBUG {
- printf("next event: ");
- print_statusnode(sn);
- }
+ calculate_dist_n_gradient(&sn, e->elev[1] + vp->target_offset, vp, *hd);
+ G_debug(3, "next event: ");
+ print_statusnode(sn);
+
switch (e->eventType) {
case ENTERING_EVENT:
+ double ax, ay;
+
/*insert node into structure */
- VIEWSHEDDEBUG {
- printf("..ENTER-EVENT: insert\n");
+ /* need either 3 elevation values or
+ * 3 gradients calculated from 3 elevation values */
+ /* need also 3 angles */
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ //sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row);
+ sn.angle[0] = e->angle;
+ calculate_event_gradient(&sn, 0, ay, ax, e->elev[0], vp, *hd);
+
+ e->eventType = CENTER_EVENT;
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_dist_n_gradient(&sn, e->elev[1], vp, *hd);
+
+ e->eventType = EXITING_EVENT;
+ calculate_event_position(*e, vp->row, vp->col, &ay, &ax);
+ sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row);
+ calculate_event_gradient(&sn, 2, ay, ax, e->elev[2], vp, *hd);
+
+ e->eventType = ENTERING_EVENT;
+
+ if (e->angle < M_PI) {
+ if (sn.angle[0] > sn.angle[1])
+ sn.angle[0] -= 2 * M_PI;
}
+ else {
+ if (sn.angle[0] > sn.angle[1]) {
+ sn.angle[1] += 2 * M_PI;
+ sn.angle[2] += 2 * M_PI;
+ }
+ }
+
+ G_debug(3, "..ENTER-EVENT: insert");
+
insert_into_status_struct(sn, status_struct);
break;
case EXITING_EVENT:
/*delete node out of status structure */
- VIEWSHEDDEBUG {
- printf("..EXIT-EVENT: delete\n");
- }
+
+ G_debug(3, "..EXIT-EVENT: delete");
+
delete_from_status_struct(status_struct, sn.dist2vp);
break;
case CENTER_EVENT:
- VIEWSHEDDEBUG {
- printf("..QUERY-EVENT: query\n");
- }
+ G_debug(3, "..QUERY-EVENT: query");
+
/*calculate visibility */
viscell.row = sn.row;
viscell.col = sn.col;
double max;
max =
- find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
+ find_max_gradient_in_status_struct(status_struct, sn.dist2vp,
+ e->angle, sn.gradient[1]);
/*the point is visible */
- if (max <= sn.gradient_offset) {
+ if (max <= sn.gradient[1]) {
viscell.angle =
- get_vertical_angle(*vp, sn, viewOptions.doCurv);
+ get_vertical_angle(*vp, sn, e->elev[1] + vp->target_offset, viewOptions.doCurv);
assert(viscell.angle >= 0);
/* viscell.vis = VISIBLE; */
add_result_to_io_visibilitygrid(visgrid, &viscell);
@@ -601,9 +675,10 @@
}
} /* for each event */
rt_stop(sweepTime);
+ G_percent(1, 1, 1);
- printf("Sweeping done.\n");
- printf("Total cells %ld, visible cells %ld (%.1f percent).\n",
+ G_message(_("Sweeping done."));
+ G_verbose_message(_("Total cells %ld, visible cells %ld (%.1f percent)."),
(long)visgrid->hd->nrows * visgrid->hd->ncols,
nvis,
(float)((float)nvis * 100 /
Modified: grass-addons/raster/r.viewshed/viewshed.h
===================================================================
--- grass-addons/raster/r.viewshed/viewshed.h 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/viewshed.h 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -44,12 +43,8 @@
#include "grid.h"
#include "eventlist.h"
-#ifdef __GRASS__
#include "grass.h"
#include <grass/iostream/ami.h>
-#else
-#include <ami.h>
-#endif
/* ------------------------------------------------------------ */
Modified: grass-addons/raster/r.viewshed/visibility.cc
===================================================================
--- grass-addons/raster/r.viewshed/visibility.cc 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/visibility.cc 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -37,14 +36,11 @@
*****************************************************************************/
#include <stdio.h>
-#ifdef __GRASS__
extern "C"
{
-#include <grass/config.h>
#include <grass/gis.h>
#include <grass/glocale.h>
}
-#endif
#include "grid.h"
#include "visibility.h"
@@ -56,7 +52,7 @@
/* viewpoint functions */
void print_viewpoint(Viewpoint vp)
{
- printf("vp=(%d, %d, %.1f) ", vp.row, vp.col, vp.elev);
+ G_debug(3, "vp=(%d, %d, %.1f) ", vp.row, vp.col, vp.elev);
return;
}
@@ -98,11 +94,8 @@
MemoryVisibilityGrid *visgrid;
-#ifdef __GRASS__
visgrid = (MemoryVisibilityGrid *) G_malloc(sizeof(MemoryVisibilityGrid));
-#else
- visgrid = (MemoryVisibilityGrid *) malloc(sizeof(MemoryVisibilityGrid));
-#endif
+
assert(visgrid);
@@ -111,11 +104,8 @@
assert(visgrid->grid);
/* create the header */
-#ifdef __GRASS__
visgrid->grid->hd = (GridHeader *) G_malloc(sizeof(GridHeader));
-#else
- visgrid->grid->hd = (GridHeader *) malloc(sizeof(GridHeader));
-#endif
+
assert(visgrid->grid->hd);
/* set the header */
@@ -125,11 +115,8 @@
alloc_grid_data(visgrid->grid);
/*allocate viewpoint */
-#ifdef __GRASS__
visgrid->vp = (Viewpoint *) G_malloc(sizeof(Viewpoint));
-#else
- visgrid->vp = (Viewpoint *) malloc(sizeof(Viewpoint));
-#endif
+
assert(visgrid->vp);
copy_viewpoint(visgrid->vp, vp);
@@ -148,18 +135,11 @@
if (visgrid->grid) {
destroy_grid(visgrid->grid);
}
-#ifdef __GRASS__
if (visgrid->vp) {
G_free(visgrid->vp);
}
G_free(visgrid);
-#else
- if (visgrid->vp) {
- free(visgrid->vp);
- }
- free(visgrid);
-#endif
return;
}
@@ -228,16 +208,12 @@
{
/* if GRASS is on, we cannot guarantee that NODATA is negative; so
we need to check */
-#ifdef __GRASS__
- int isnull = G_is_f_null_value(&x);
+ int isnull = G_is_null_value(&x, G_SURFACE_TYPE);
if (isnull)
return 0;
else
return (x >= 0);
-#else
- return (x >= 0);
-#endif
}
int is_invisible_not_nodata(float x)
{
@@ -274,10 +250,6 @@
}
-
-
-
-
/* ------------------------------------------------------------ */
/* visgrid is the structure that records the visibility information
after the sweep is done. Use it to write the visibility output
@@ -287,7 +259,6 @@
ViewOptions viewOptions, Viewpoint vp)
{
-#ifdef __GRASS__
if (viewOptions.outputMode == OUTPUT_BOOL)
save_grid_to_GRASS(visgrid->grid, viewOptions.outputfname, CELL_TYPE,
booleanVisibilityOutput);
@@ -299,22 +270,6 @@
save_vis_elev_to_GRASS(visgrid->grid, viewOptions.inputfname,
viewOptions.outputfname,
vp.elev + viewOptions.obsElev);
-#else
- if (viewOptions.outputMode == OUTPUT_BOOL)
- save_grid_to_arcascii_file(visgrid->grid, viewOptions.outputfname,
- booleanVisibilityOutput);
- else if (viewOptions.outputMode == OUTPUT_ANGLE)
- save_grid_to_arcascii_file(visgrid->grid, viewOptions.outputfname,
- angleVisibilityOutput);
- else {
- /* elevation output */
- printf("Elevation output not implemented in the standalone version.");
- printf("Output in angle mode\n");
- save_grid_to_arcascii_file(visgrid->grid, viewOptions.outputfname,
- angleVisibilityOutput);
- /* if you want elevation output, use grass */
- }
-#endif
free_inmem_visibilitygrid(visgrid);
@@ -323,42 +278,29 @@
-
-
/* ------------------------------------------------------------ */
/* IOVisibilityGrid functions */
/* ------------------------------------------------------------ */
-
-
/* ------------------------------------------------------------ */
/*create grid from given header and viewpoint */
IOVisibilityGrid *init_io_visibilitygrid(GridHeader hd, Viewpoint vp)
{
IOVisibilityGrid *visgrid;
-#ifdef __GRASS__
visgrid = (IOVisibilityGrid *) G_malloc(sizeof(IOVisibilityGrid));
-#else
- visgrid = (IOVisibilityGrid *) malloc(sizeof(IOVisibilityGrid));
-#endif
+
assert(visgrid);
/*header */
-#ifdef __GRASS__
visgrid->hd = (GridHeader *) G_malloc(sizeof(GridHeader));
-#else
- visgrid->hd = (GridHeader *) malloc(sizeof(GridHeader));
-#endif
+
assert(visgrid->hd);
copy_header(visgrid->hd, hd);
/*viewpoint */
-#ifdef __GRASS__
visgrid->vp = (Viewpoint *) G_malloc(sizeof(Viewpoint));
-#else
- visgrid->vp = (Viewpoint *) malloc(sizeof(Viewpoint));
-#endif
+
assert(visgrid->vp);
copy_viewpoint(visgrid->vp, vp);
@@ -376,7 +318,7 @@
void free_io_visibilitygrid(IOVisibilityGrid * grid)
{
assert(grid);
-#ifdef __GRASS__
+
if (grid->hd)
G_free(grid->hd);
if (grid->vp)
@@ -385,24 +327,12 @@
delete grid->visStr;
G_free(grid);
-#else
- if (grid->hd)
- free(grid->hd);
- if (grid->vp)
- free(grid->vp);
- if (grid->visStr)
- delete grid->visStr;
- free(grid);
-#endif
return;
}
-
-
-
/* ------------------------------------------------------------ */
/*write cell to stream */
void add_result_to_io_visibilitygrid(IOVisibilityGrid * visgrid,
@@ -420,7 +350,6 @@
}
-
/* ------------------------------------------------------------ */
/*compare function, (i,j) grid order */
int IJCompare::compare(const VisCell & a, const VisCell & b)
@@ -464,14 +393,12 @@
}
-
/* ------------------------------------------------------------ */
void
save_io_visibilitygrid(IOVisibilityGrid * visgrid,
ViewOptions viewOptions, Viewpoint vp)
{
-#ifdef __GRASS__
if (viewOptions.outputMode == OUTPUT_BOOL)
save_io_visibilitygrid_to_GRASS(visgrid, viewOptions.outputfname,
CELL_TYPE, booleanVisibilityOutput);
@@ -484,97 +411,9 @@
save_io_vis_and_elev_to_GRASS(visgrid, viewOptions.inputfname,
viewOptions.outputfname,
vp.elev + viewOptions.obsElev);
-#else
- if (viewOptions.outputMode == OUTPUT_BOOL)
- save_io_visibilitygrid_to_arcascii(visgrid, viewOptions.outputfname,
- booleanVisibilityOutput);
- else if (viewOptions.outputMode == OUTPUT_ANGLE)
- save_io_visibilitygrid_to_arcascii(visgrid, viewOptions.outputfname,
- angleVisibilityOutput);
- else {
- /* elevation output */
- printf("Elevation output not implemented in the standalone version.");
- printf("Output in angle mode\n");
- save_io_visibilitygrid_to_arcascii(visgrid, viewOptions.outputfname,
- angleVisibilityOutput);
- /* if you want elevation output, use grass */
- }
-#endif
/*free visibiliyty grid */
free_io_visibilitygrid(visgrid);
return;
}
-
-
-
-
-/* ------------------------------------------------------------ */
-/*write visibility grid to arcascii file. assume all cells that are
- not in stream are NOT visible. assume stream is sorted in (i,j)
- order. for each value x it writes to grass fun(x) */
-void
-save_io_visibilitygrid_to_arcascii(IOVisibilityGrid * visgrid,
- char *outputfname, float (*fun) (float))
-{
-
- assert(visgrid && outputfname);
-
- /*open file */
- FILE *fp = fopen(outputfname, "w");
-
- assert(fp);
-
- /*write header */
- print_grid_header(fp, visgrid->hd);
-
- /*sort the stream in (i,j) order */
- /*sortVisGrid(visgrid); */
-
- /* set up visibility stream */
- AMI_STREAM < VisCell > *vstr = visgrid->visStr;
- off_t streamLen, counter = 0;
-
- streamLen = vstr->stream_len();
- vstr->seek(0);
-
-
- /*read the first element */
- AMI_err ae;
- VisCell *curResult = NULL;
-
- if (streamLen > 0) {
- ae = vstr->read_item(&curResult);
- assert(ae == AMI_ERROR_NO_ERROR);
- counter++;
- }
- for (dimensionType i = 0; i < visgrid->hd->nrows; i++) {
- for (dimensionType j = 0; j < visgrid->hd->ncols; j++) {
-
- if (curResult->row == i && curResult->col == j) {
- /* this cell is present in the visibility stream; it must
- be either visible, or NODATA */
- fprintf(fp, "%.1f ", fun(curResult->angle));
-
- /*read next element of stream */
- if (counter < streamLen) {
- ae = vstr->read_item(&curResult);
- assert(ae == AMI_ERROR_NO_ERROR);
- counter++;
- }
- }
- else {
- /* this cell is not in stream, then it is invisible */
- fprintf(fp, "%.1f ", fun(INVISIBLE));
- }
- } /* for j */
- fprintf(fp, "\n");
- } /* for i */
-
- /* assert that we read teh entire file, otherwise something is
- wrong */
- assert(counter == streamLen);
- fclose(fp);
- return;
-}
Modified: grass-addons/raster/r.viewshed/visibility.h
===================================================================
--- grass-addons/raster/r.viewshed/visibility.h 2011-05-05 12:24:33 UTC (rev 46198)
+++ grass-addons/raster/r.viewshed/visibility.h 2011-05-06 07:30:56 UTC (rev 46199)
@@ -8,8 +8,9 @@
* Ported to GRASS by William Richard -
* wkrichar at bowdoin.edu or willster3021 at gmail.com
+ * Markus Metz: surface interpolation
*
- * Date: july 2008
+ * Date: april 2011
*
* PURPOSE: To calculate the viewshed (the visible cells in the
* raster) for the given viewpoint (observer) location. The
@@ -17,12 +18,10 @@
* considered visible to each other if the cells where they belong are
* visible to each other. Two cells are visible to each other if the
* line-of-sight that connects their centers does not intersect the
- * terrain. The height of a cell is assumed to be constant, and the
- * terrain is viewed as a tesselation of flat cells. This model is
- * suitable for high resolution rasters; it may not be accurate for
- * low resolution rasters, where it may be better to interpolate the
- * height at a point based on the neighbors, rather than assuming
- * cells are "flat". The viewshed algorithm is efficient both in
+ * terrain. The terrain is NOT viewed as a tesselation of flat cells,
+ * i.e. if the line-of-sight does not pass through the cell center,
+ * elevation is determined using bilinear interpolation.
+ * The viewshed algorithm is efficient both in
* terms of CPU operations and I/O operations. It has worst-case
* complexity O(n lg n) in the RAM model and O(sort(n)) in the
* I/O-model. For the algorithm and all the other details see the
@@ -39,13 +38,9 @@
#ifndef visibility_h
#define visibility_h
-#ifdef __GRASS__
#include <grass/config.h>
#include <grass/iostream/ami.h>
-#else
-#include <ami.h>
-#endif
#include "grid.h"
@@ -249,15 +244,7 @@
save_io_visibilitygrid(IOVisibilityGrid * visgrid,
ViewOptions viewoptions, Viewpoint vp);
-/* write visibility grid to arcascii file. assume all cells that are
- not in stream are NOT visible. assume stream is sorted. calls fun
- on every individual cell */
-void
-save_io_visibilitygrid_to_arcascii(IOVisibilityGrid * visgrid,
- char *outfname, float (*fun) (float));
-
-
/*sort stream in grid (i,j) order */
void sort_io_visibilitygrid(IOVisibilityGrid * visGrid);
More information about the grass-commit
mailing list