[GRASS-SVN] r32752 - in grass/trunk/raster: . r.viewshed
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Aug 13 15:57:41 EDT 2008
Author: ltoma
Date: 2008-08-13 15:57:41 -0400 (Wed, 13 Aug 2008)
New Revision: 32752
Added:
grass/trunk/raster/r.viewshed/
grass/trunk/raster/r.viewshed/Makefile
grass/trunk/raster/r.viewshed/distribute.cc
grass/trunk/raster/r.viewshed/distribute.h
grass/trunk/raster/r.viewshed/eventlist.cc
grass/trunk/raster/r.viewshed/eventlist.h
grass/trunk/raster/r.viewshed/grass.cc
grass/trunk/raster/r.viewshed/grass.h
grass/trunk/raster/r.viewshed/grid.cc
grass/trunk/raster/r.viewshed/grid.h
grass/trunk/raster/r.viewshed/main.cc
grass/trunk/raster/r.viewshed/print_message.cc
grass/trunk/raster/r.viewshed/print_message.h
grass/trunk/raster/r.viewshed/r.viewshed.html
grass/trunk/raster/r.viewshed/rbbst.cc
grass/trunk/raster/r.viewshed/rbbst.h
grass/trunk/raster/r.viewshed/statusstructure.cc
grass/trunk/raster/r.viewshed/statusstructure.h
grass/trunk/raster/r.viewshed/sweep1.png
grass/trunk/raster/r.viewshed/sweep2.png
grass/trunk/raster/r.viewshed/viewshed.cc
grass/trunk/raster/r.viewshed/viewshed.h
grass/trunk/raster/r.viewshed/visibility.cc
grass/trunk/raster/r.viewshed/visibility.h
Log:
First working version. Uses grass/lib/iostrean, Makefile and include changes from Paul kelly
Added: grass/trunk/raster/r.viewshed/Makefile
===================================================================
--- grass/trunk/raster/r.viewshed/Makefile (rev 0)
+++ grass/trunk/raster/r.viewshed/Makefile 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,52 @@
+# Path to the module directory, i.e. the grass sourcecode main path.
+# Relative or absolute.
+MODULE_TOPDIR = ../../
+
+PGM = r.viewshed
+
+# Includes the grass make-System.
+# For further information look in the Multi.make file.
+include $(MODULE_TOPDIR)/include/Make/Multi.make
+
+
+SOURCES = main.cc distribute.cc eventlist.cc grid.cc grass.cc viewshed.cc \
+ rbbst.cc statusstructure.cc visibility.cc \
+ print_message.cc
+HEADERS = distribute.h eventlist.h grid.h grass.h viewshed.h \
+ rbbst.h statusstructure.h visibility.h print_message.h
+
+
+OBJARCH=OBJ.$(ARCH)
+OBJ := $(patsubst %.cc,$(OBJARCH)/%.o,$(SOURCES))
+
+LIBS = $(GISLIB) $(IOSTREAMLIB)
+DEPLIBS = $(DEPGISLIB) $(IOSTREAMDEP)
+
+# Using GNU c++ compiler.
+CXX = g++
+
+# Set compiler and load flags.
+# (See 'man g++' for help).
+CXXFLAGS += -O3 -DNO_STATS #-DNDEBUG
+CXXFLAGS += $(WARNING_FLAGS) -DUSER=\"$(USER)\" -D__GRASS__
+# CXXFLAGS += -DPEARL
+CXXFLAGS += -D_FILE_OFFSET_BITS=64 -D_LARGEFILE_SOURCE -fmessage-length=0
+CXXFLAGS += -ffast-math -funroll-loops
+
+
+WARNING_FLAGS = -Wall -Wformat -Wparentheses -Wpointer-arith -Wno-conversion \
+ -Wreturn-type \
+ -Wcomment -Wno-sign-compare -Wno-unused
+#-Wimplicit-int -Wimplicit-function-declaration
+
+# Make rules.
+# Default is the 'default' target for the grass make system.
+$(OBJARCH)/%.o:%.cc
+ $(CXX) -c $(CXXFLAGS) $< -o $@
+
+
+default: $(BIN)/$(PGM)$(EXE)
+ $(MAKE) htmlcmd
+
+$(BIN)/$(PGM)$(EXE): $(OBJ) $(DEPLIBS)
+ $(CXX) $(LDFLAGS) -o $@ $(OBJ) $(LIBS)
Property changes on: grass/trunk/raster/r.viewshed/Makefile
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/distribute.cc
===================================================================
--- grass/trunk/raster/r.viewshed/distribute.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/distribute.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,1159 @@
+/****************************************************************************
+ *
+ * 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 <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#ifdef __GRASS__
+extern "C"
+{
+#include <grass/config.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+}
+
+/* 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
+ is included by ami.h */
+
+#include "distribute.h"
+#include "viewshed.h"
+#include "visibility.h"
+#include "eventlist.h"
+#include "statusstructure.h"
+#include "print_message.h"
+#ifdef __GRASS__
+#include "grass.h"
+#endif
+
+
+#define DISTRIBDEBUG if(0)
+#define SOLVEINMEMDEBUG if(0)
+#define DEBUGINIT if(0)
+#define PRINTWARNING if(0)
+#define BND_DEBUG if(0)
+
+
+#define ANGLE_FACTOR 1
+#define EPSILON .00000001
+#define PRINT_DISTRIBUTE if(0)
+
+
+
+////////////////////////////////////////////////////////////////////////
+/* distribution sweep: write results to visgrid.
+ */
+IOVisibilityGrid *distribute_and_sweep(char* inputfname,
+ GridHeader * hd,
+ Viewpoint *vp,
+ ViewOptions viewOptions) {
+
+ assert(inputfname && hd && vp);
+ print_message("Start distributed sweeping.\n");
+ fflush(stdout);
+
+ /* ------------------------------ */
+ /*initialize the visibility grid */
+ IOVisibilityGrid *visgrid;
+ visgrid = init_io_visibilitygrid(*hd, *vp);
+ printf("distribute_and_sweep: visgrid=%s\n", visgrid->visStr->name());
+
+
+ /* ------------------------------ */
+ /*construct event list corresp to the input file and viewpoint */
+ Rtimer initEventTime;
+ rt_start(initEventTime);
+ AMI_STREAM < AEvent > *eventList;
+#ifdef __GRASS__
+ eventList = grass_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());
+
+
+ /* ------------------------------ */
+ /*sort the events concentrically */
+ Rtimer sortEventTime;
+ rt_start(sortEventTime);
+ PRINT_DISTRIBUTE {
+ print_message("sorting events by distance from viewpoint..");
+ fflush(stdout);
+ }
+
+ sort_event_list_by_distance(&eventList, *vp);
+ PRINT_DISTRIBUTE {
+ print_message("..sorting done.\n");
+ fflush(stdout);
+ }
+
+ /* 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());
+
+
+
+ /* ------------------------------ */
+ /* start distribution */
+ long nvis; /*number of cells visible. Returned by distribute_sector */
+
+ Rtimer sweepTime;
+ rt_start(sweepTime);
+ /*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);
+ rt_stop(sweepTime);
+
+
+ /* ------------------------------ */
+ /*cleanup */
+ print_message("Distribution sweeping done.\n");
+ fflush(stdout);
+
+ printf("Total cells %ld, visible cells %ld (%.1f percent).\n",
+ (long)visgrid->hd->nrows * visgrid->hd->ncols,
+ nvis,
+ (float)((float)nvis * 100 /
+ (float)(visgrid->hd->nrows * visgrid->hd->ncols)));
+
+ print_viewshed_timings(initEventTime, sortEventTime, sweepTime);
+
+ return visgrid;
+}
+
+
+
+
+
+
+//***********************************************************************
+/* distribute recursively the events and write results to
+ visgrid. eventList is a list of events sorted by distance that fall
+ within angle boundaries start_angle and end_angle; enterBndEvents
+ is a stream that contains all the ENTER events that are not in this
+ sector, but their corresponding Q or X events are in this sector.
+
+ When problem is small enough, solve it in memory and write results to
+ visgrid.
+
+ invariant: distribute_sector deletes eventList and enterBndEvents
+ */
+unsigned long distribute_sector(AMI_STREAM < AEvent > *eventList,
+ AMI_STREAM < AEvent > *enterBndEvents,
+ double start_angle, double end_angle,
+ IOVisibilityGrid * visgrid, Viewpoint * vp,
+ ViewOptions viewOptions) {
+
+
+ assert(eventList && visgrid && vp);
+ /*enterBndEvents may be NULL first time */
+
+ PRINT_DISTRIBUTE printf("***DISTRIBUTE sector [%.4f, %.4f]***\n",
+ start_angle, end_angle);
+ printf("initial_gradient: %lf\n", SMALLEST_GRADIENT);
+ printf("eventlist: %s\n", eventList->sprint());
+ if (enterBndEvents)
+ printf("BndEvents: %s\n", enterBndEvents->sprint());
+ PRINT_DISTRIBUTE LOG_avail_memo();
+
+ unsigned long nvis = 0;
+
+ /*******************************************************
+ BASE CASE
+ *******************************************************/
+ if (eventList->stream_len() * sizeof(AEvent) <
+ MM_manager.memory_available()) {
+ if (enterBndEvents) {
+ nvis += solve_in_memory(eventList, enterBndEvents,
+ start_angle,end_angle,visgrid,vp,viewOptions);
+ return nvis;
+ }
+ else {
+ /*we are here if the problem fits in memory, and
+ //enterBNdEvents==NULL---this only happens the very first time;
+ //technically at this point we should do one pass though the
+ //data and collect the events that cross the 1st and 4th
+ //quadrant; instead we force recursion, nect= 2 */
+ }
+ }
+
+ /*else, must recurse */
+ PRINT_DISTRIBUTE print_message("in EXTERNAL memory\n");
+
+ /*compute number of sectors */
+ int nsect = compute_n_sectors();
+ assert(nsect > 1);
+
+ /*an array of streams, one for each sector; sector[i] will keep all
+ the cells that are completely inside sector i */
+ AMI_STREAM < AEvent > *sector = new AMI_STREAM < AEvent >[nsect];
+
+ /*the array of gradient values, one for each sector; the gradient is
+ the gradient of the center of a cell that spans the sector
+ completely */
+ double *high = new double[nsect];
+
+ for (int i = 0; i < nsect; i++)
+ high[i] = SMALLEST_GRADIENT;
+
+ /*an array of streams, one for each stream boundary; sectorBnd[0]
+ will keep all the cells crossing into sector 0 from below; and so on. */
+ AMI_STREAM < AEvent > *sectorBnd = new AMI_STREAM < AEvent >[nsect];
+
+ /*keeps stats for each sector */
+ long *total = new long[nsect];
+ long *insert = new long[nsect];
+ long *drop = new long[nsect];
+ long *bndInsert = new long[nsect];
+ long *bndDrop = new long[nsect];
+
+ for (int i = 0; i < nsect; i++)
+ total[i] = insert[i] = drop[i] = bndInsert[i] = bndDrop[i] = 0;
+ long longEvents = 0;
+
+ /*******************************************************
+ CONCENTRIC SWEEP
+ *******************************************************/
+ AEvent *e;
+ AMI_err ae;
+ double exit_angle, enter_angle;
+ int exit_s, s, enter_s;
+ long boundaryEvents = 0;
+ off_t nbEvents = eventList->stream_len();
+
+ 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);
+ assert(is_inside(e, start_angle, end_angle));
+
+ /*compute its sector */
+ s = get_event_sector(e->angle, start_angle, end_angle, nsect);
+ /*detect boundary cases ==> precision issues */
+ if (is_almost_on_boundary(e->angle, s, start_angle, end_angle, nsect)) {
+ 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);
+ }
+ }
+
+ DISTRIBDEBUG printf("event %7lu: ", (unsigned long)i);
+ DISTRIBDEBUG print_event(*e);
+ DISTRIBDEBUG printf("d=%8.1f, ",
+ get_square_distance_from_viewpoint(*e, *vp));
+ DISTRIBDEBUG printf("s=%3d ", s);
+
+ assert(is_inside(s, nsect));
+ total[s]++;
+
+ /*insert event in sector if not occluded */
+ insert_event_in_sector(e, s, §or[s], high[s], vp, insert, drop);
+
+ switch (e->eventType) {
+ case CENTER_EVENT:
+ break;
+
+ case ENTERING_EVENT:
+ /*find its corresponding exit event and its sector */
+ 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) ",
+ e->angle, s, exit_angle, exit_s);
+ /*note: exit_s can be -1 (outside) */
+ if (exit_s == s) {
+ /*short event, fit in sector s */
+
+ }
+ else if (exit_s == (s + 1) % nsect || (exit_s + 1) % nsect == s) {
+ /*semi-short event; insert in sector s, and in sector boundary s+1
+ NOTE: to avoid precision issues, the events are inserted
+ when processing the EXIT_EVENT
+ insertEventInSector(e, (s+1)%nsect, §orBnd[(s+1)%nsect],
+ high[(s+1)%nsect], vp,bndInsert, bndDrop); */
+
+ }
+ else {
+ /*long event; insert in sector s, and in sector boundary exit_s */
+ process_long_cell(s, exit_s, nsect, vp, e, high);
+ longEvents++;
+ /*insertEventInSector(e, exit_s, §orBnd[exit_s], high[exit_s],
+ // vp, bndInsert, bndDrop); */
+ }
+ break;
+
+ case EXITING_EVENT:
+ /*find its corresponding enter event and its sector */
+ 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) ",
+ e->angle, s, enter_angle, enter_s);
+
+ /*don't need to check spanned sectors because it is done on its
+ //ENTER event; actually...you do, because its enter event may
+ //not be in this sector==> enter_s = -1 (outside) */
+ if (enter_s == s) {
+ /*short event, fit in sector */
+ }
+ else if (enter_s == (s + 1) % nsect || (enter_s + 1) % nsect == s) {
+ /*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);
+ }
+ insert_event_in_sector(e, s, §orBnd[s], high[s],
+ vp, bndInsert, bndDrop);
+
+ }
+ else {
+ /*long event */
+ process_long_cell(enter_s, s, nsect, vp, e, high);
+ 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);
+ }
+ insert_event_in_sector(e, s, §orBnd[s], high[s],
+ vp, bndInsert, bndDrop);
+ }
+ break;
+
+ } /*switch event-type */
+
+ DISTRIBDEBUG print_message("\n");
+ } /*for event i */
+
+
+ /*distribute the enterBnd events to the boundary streams of the
+ //sectors; note: the boundary streams are not sorted by distance. */
+ if (enterBndEvents)
+ distribute_bnd_events(enterBndEvents, sectorBnd, nsect, vp,
+ start_angle, end_angle, high, bndInsert,
+ bndDrop);
+
+ /*sanity checks */
+ PRINT_DISTRIBUTE printf("boundary events in distribution: %ld\n",
+ boundaryEvents);
+ print_sector_stats(nbEvents, nsect, high, total, insert, drop, sector,
+ sectorBnd, bndInsert,
+ longEvents, start_angle, end_angle);
+ /*cleanup after sector stats */
+ delete[]total;
+ delete[]insert;
+ delete[]drop;
+ delete[]high;
+ delete[]bndInsert;
+ delete[]bndDrop;
+
+ /*we finished processing this sector, delete the event list */
+ delete eventList;
+
+ if (enterBndEvents)
+ 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]);
+
+ sector[i].persist(PERSIST_PERSISTENT);
+ sectorBnd[i].name(§orBndName[i]);
+ PRINT_DISTRIBUTE printf("saving BndStr %d: %s\n", i,
+ sectorBndName[i]);
+ sectorBnd[i].persist(PERSIST_PERSISTENT);
+ }
+
+ /*delete [] sector;
+ //does this call delete on every single stream?
+ //for (int i=0; i< nsect; i++ ) delete sector[i]; */
+ delete[]sector;
+ delete[]sectorBnd;
+ /*LOG_avail_memo(); */
+
+
+ /*solve recursively each sector */
+ for (int i = 0; i < nsect; i++) {
+
+ /*recover stream */
+ PRINT_DISTRIBUTE printf("\nopening sector stream %s ", sectorName[i]);
+
+ AMI_STREAM < AEvent > *str =
+ new AMI_STREAM < AEvent > (sectorName[i]);
+ assert(str);
+ PRINT_DISTRIBUTE printf(" len=%lu\n",
+ (unsigned long)str->stream_len());
+ /*recover boundary stream */
+ PRINT_DISTRIBUTE printf("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",
+ (unsigned long)bndStr->stream_len());
+
+
+ nvis += distribute_sector(str, bndStr,
+ start_angle+i*((end_angle-start_angle)/nsect),
+ start_angle+(i+1)*((end_angle-start_angle)/nsect),
+ visgrid, vp, 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",
+ start_angle, end_angle);
+
+ return nvis;
+}
+
+
+
+
+/***********************************************************************
+ enterBndEvents is a stream of events that cross into the sector's
+ (first) boundary; they must be distributed to the boundary streams
+ of the sub-sectors of this sector. Note: the boundary streams of
+ the sub-sectors may not be empty; as a result, events get appended
+ at the end, and they will not be sorted by distance from the
+ vp.
+ */
+void
+distribute_bnd_events(AMI_STREAM < AEvent > *bndEvents,
+ AMI_STREAM < AEvent > *sectorBnd, int nsect,
+ Viewpoint * vp, double start_angle, double end_angle,
+ double *high, long *insert, long *drop)
+{
+
+ PRINT_DISTRIBUTE printf("Distribute boundary of sector [ %.4f, %.4f] ",
+ start_angle, end_angle);
+ assert(bndEvents && sectorBnd && vp && high && insert && drop);
+ AEvent *e;
+ AMI_err ae;
+ double exit_angle;
+ int exit_s;
+ off_t nbEvents = bndEvents->stream_len();
+
+ bndEvents->seek(0);
+ for (off_t i = 0; i < nbEvents; i++) {
+
+ /*get out one event at a time */
+ ae = bndEvents->read_item(&e);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ /*must be outside, but better not check to avoid precision issues
+ //assert(!is_inside(e, start_angle, end_angle));
+
+ //each event must be an ENTER event that falls in a diff sector
+ //than its EXIT */
+ assert(e->eventType == ENTERING_EVENT);
+
+ /*find its corresponding exit event and its sector */
+ exit_angle = calculate_exit_angle(e->row, e->col, vp);
+ exit_s = get_event_sector(exit_angle, start_angle, end_angle, nsect);
+
+ /*exit_s cannot be outside sector; though we have to be careful
+ //with precision */
+ assert(is_inside(exit_s, nsect));
+
+ /*insert this event in the boundary stream of this sub-sector */
+ insert_event_in_sector(e, exit_s, §orBnd[exit_s], high[exit_s],
+ vp, insert, drop);
+
+ } /*for i */
+
+ PRINT_DISTRIBUTE
+ printf("Distribute boundary of sector [ %.4f, %.4f] done.\n",
+ start_angle, end_angle);
+
+ return;
+}
+
+
+
+
+//***********************************************************************
+/* Solves a segment it inemory. 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
+ corresponding EXIT events are inside the sector. start_angle and
+ end_angle are the boundaries of the sector, and visgrid is the grid
+ to which the visible/invisible cells must be written. The sector is
+ 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,
+ Viewpoint * vp, ViewOptions viewOptions)
+{
+
+ assert(eventList && visgrid && vp);
+ PRINT_DISTRIBUTE print_message("solve INTERNAL memory\n");
+
+ unsigned long nvis = 0; /*number of visible cells */
+
+ printf("solve_in_memory: eventlist: %s\n", eventList->sprint());
+ if (enterBndEvents)
+ printf("BndEvents: %s\n", enterBndEvents->sprint());
+
+ if (eventList->stream_len() == 0) {
+ delete eventList;
+ if (enterBndEvents) delete enterBndEvents;
+ return nvis;
+ }
+
+ /*sort the events radially */
+ sort_event_list(&eventList);
+
+ /*create the status structure */
+ StatusList *status_struct = create_status_struct();
+
+ /*initialize status structure with all ENTER events whose EXIT
+ //events is inside the sector */
+ AEvent *e;
+ AMI_err ae;
+ StatusNode sn;
+ int inevents = 0;
+
+ /*
+ eventList->seek(0);
+ double enter_angle;
+ ae = eventList->read_item(&e);
+ while (ae == AMI_ERROR_NO_ERROR) {
+ if (e->eventType == EXITING_EVENT) {
+ enter_angle = calculateEnterAngle(e->row, e->col, vp);
+ //to get around precision problems, insert cell in active
+ //structure when close to the boundary --- these may cause a
+ //double insert, but subsequent ENTER events close to the
+ //boundary will be careful not to insert
+ if (!is_inside(enter_angle, start_angle, end_angle) ||
+ is_almost_on_boundary(enter_angle, start_angle)) {
+ DEBUGINIT {printf("inserting "); print_event(*e); printf("\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;
+ calculateDistNGradient(&sn, vp);
+ insertIntoStatusStruct(sn, status_struct);
+ inevents++;
+ }
+ }
+ ae = eventList->read_item(&e);
+ }
+ assert(ae == AMI_ERROR_END_OF_STREAM);
+ */
+ if (enterBndEvents) {
+ 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");
+ }
+ /*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);
+ insert_into_status_struct(sn, status_struct);
+ }
+ }
+ PRINT_DISTRIBUTE {
+ printf("initialized active structure with %d events\n", inevents);
+ fflush(stdout);
+ }
+
+
+ /*sweep the event list */
+ VisCell viscell;
+ off_t nbEvents = eventList->stream_len();
+
+ /*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);
+ }
+
+ sn.col = e->col;
+ sn.row = e->row;
+ sn.elev = e->elev;
+ /*calculate Distance to VP and Gradient */
+ calculate_dist_n_gradient(&sn, vp);
+
+ switch (e->eventType) {
+ case ENTERING_EVENT:
+ /*insert node into structure */
+ SOLVEINMEMDEBUG {
+ print_message("..ENTER-EVENT: insert\n");
+ }
+ /*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); */
+ insert_into_status_struct(sn, status_struct);
+ break;
+
+ case EXITING_EVENT:
+ /*delete node out of status structure */
+ SOLVEINMEMDEBUG {
+ print_message("..EXIT-EVENT: delete\n");
+ /*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,
+ enter_angle);
+ }
+ delete_from_status_struct(status_struct, sn.dist2vp);
+ break;
+
+ case CENTER_EVENT:
+ SOLVEINMEMDEBUG {
+ print_message("..QUERY-EVENT: query\n");
+ }
+ /*calculate visibility
+
+ //note: if there is nothing in the status structure, it means
+ //there is no prior event to occlude it, so this cell is
+ //VISIBLE. this is taken care of in the status structure --- if
+ //a query event comes when the structure is empty, the query
+ //returns minimum gradient */
+ double max;
+
+ max = find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
+
+ viscell.row = sn.row;
+ viscell.col = sn.col;
+
+ if (max <= sn.gradient) {
+ /*the point is visible */
+ viscell.angle = get_vertical_angle(*vp, sn, viewOptions.doCurv);
+ assert(viscell.angle >0);
+ /* viscell.vis = VISIBLE; */
+ add_result_to_io_visibilitygrid(visgrid, &viscell);
+ /*make sure nvis is correct*/
+ nvis++;
+ }
+ else {
+ /* else the point is invisible; we do not write it to the
+ //visibility stream, because we only record visible and nodata
+ //values to the stream */
+ /* viscell.vis = INVISIBLE; */
+ /* add_result_to_io_visibilitygrid(visgrid, &viscell); */
+ }
+ break;
+ }
+ } /* for each event */
+
+
+ PRINT_DISTRIBUTE print_message("in memory sweeping done.\n");
+
+ PRINT_DISTRIBUTE
+ printf("Total cells %lu, visible cells %lu (%.1f percent).\n",
+ (unsigned long)eventList->stream_len(), nvis,
+ (float)((float)nvis * 100 / (float)(eventList->stream_len())));
+
+ /*cleanup */
+ delete_status_structure(status_struct);
+
+
+
+ /*invariant: must delete its eventList */
+ delete eventList;
+
+ if (enterBndEvents)
+ delete enterBndEvents;
+
+ return nvis;
+}
+
+
+
+/***********************************************************************
+ //returns 1 if enter angle is within epsilon from boundary angle*/
+int is_almost_on_boundary(double angle, double boundary_angle)
+{
+ /*printf("is_almost_on_boundary: %f (%f) %f\n", angle, angle-2*M_PI, boundary_angle); */
+ return (fabs(angle - boundary_angle) < EPSILON) ||
+ (fabs(angle - 2 * M_PI - boundary_angle) < EPSILON);
+}
+
+
+/***********************************************************************
+ // returns 1 if angle is within epsilon the boundaries of sector s*/
+int
+is_almost_on_boundary(double angle, int s, double start_angle,
+ double end_angle, int nsect)
+{
+ /*the boundaries of sector s */
+ double ssize = (end_angle - start_angle) / nsect;
+
+ return is_almost_on_boundary(angle, s * ssize) ||
+ is_almost_on_boundary(angle, (s + 1) * ssize);
+}
+
+
+/***********************************************************************
+ returns true if the event is inside the given sector */
+int is_inside(AEvent * e, double start_angle, double end_angle)
+{
+ assert(e);
+ return (e->angle >= (start_angle - EPSILON) &&
+ e->angle <= (end_angle + EPSILON));
+}
+
+/***********************************************************************
+ returns true if this angle is inside the given sector */
+int is_inside(double angle, double start_angle, double end_angle)
+{
+ return (angle >= (start_angle - EPSILON) &&
+ angle <= (end_angle + EPSILON));
+}
+
+
+
+/***********************************************************************
+ return the start angle of the i-th sector. Assuming that
+[start..end] is split into nsectors */
+double
+get_sector_start(int i, double start_angle, double end_angle, int nsect)
+{
+ assert(is_inside(i, nsect));
+ return start_angle + i * ((end_angle - start_angle) / nsect);
+}
+
+
+
+/***********************************************************************
+ return the start angle of the i-th sector. Assuming that
+[start..end] is split into nsectors */
+double get_sector_end(int i, double start_angle, double end_angle, int nsect)
+{
+ assert(is_inside(i, nsect));
+ return start_angle + (i + 1) * ((end_angle - start_angle) / nsect);
+}
+
+
+
+/***********************************************************************
+ return 1 is s is inside sector; that is, if it is not -1 */
+int is_inside(int s, int nsect)
+{
+ return (s >= 0 && s < nsect);
+}
+
+
+/***********************************************************************
+ the event e spans sectors from start_s to end_s; Action: update
+ high[] for each spanned sector. start_s and both_s can be -1, which
+ means outside given sector---in that case long cell spans to the
+ boundary of the sector.
+ */
+void
+process_long_cell(int start_s, int end_s, int nsect,
+ Viewpoint * vp, AEvent * e, double *high)
+{
+
+ DISTRIBDEBUG printf("LONG CELL: spans [%3d, %3d] ", start_s, end_s);
+ double ctrgrad = calculate_center_gradient(e, vp);
+
+ /*ENTER event is outside */
+ if (start_s == -1) {
+ assert(e->eventType == EXITING_EVENT);
+ assert(is_inside(end_s, nsect));
+ /*span from 0 to end_s */
+ for (int j = 0; j < end_s; j++) {
+ if (high[j] < ctrgrad) {
+ /*printf("update high[%d] from %.2f to %.2f ", j, high[j], ctrgrad); */
+ high[j] = ctrgrad;
+ }
+ }
+ return;
+ }
+
+ /*EXIT event is outside */
+ if (end_s == -1) {
+ assert(e->eventType == ENTERING_EVENT);
+ assert(is_inside(start_s, nsect));
+ /*span from start_s to nsect */
+ for (int j = start_s + 1; j < nsect; j++) {
+ if (high[j] < ctrgrad) {
+ /*printf("update high[%d] from %.2f to %.2f ", j, high[j], ctrgrad); */
+ high[j] = ctrgrad;
+ }
+ }
+ return;
+ }
+
+ /*the usual scenario, both inside sector */
+ if (start_s < end_s) {
+ /*we must update high[] in start_s+1..end_s-1 */
+ for (int j = start_s + 1; j < end_s; j++) {
+ if (high[j] < ctrgrad) {
+ /*printf("update high[%d] from %.2f to %.2f ", j, high[j], ctrgrad); */
+ high[j] = ctrgrad;
+ }
+ }
+ return;
+ }
+ else {
+ /*start_s > end_s: we must insert in [start_s..nsect] and [0, end_s] */
+ for (int j = start_s + 1; j < nsect; j++) {
+ if (high[j] < ctrgrad) {
+ /*printf("update high[%d] from %.2f to %.2f ", j, high[j], ctrgrad); */
+ high[j] = ctrgrad;
+ }
+ }
+ for (int j = 0; j < end_s; j++) {
+ if (high[j] < ctrgrad) {
+ /*printf("update high[%d] from %.2f to %.2f ", j, high[j], ctrgrad); */
+ high[j] = ctrgrad;
+ }
+ }
+ }
+ return;
+}
+
+/***********************************************************************
+ prints how many events were inserted and dropped in each sector */
+void
+print_sector_stats(off_t nevents, int nsect, double *high,
+ long *total,
+ long *insert, long *drop, AMI_STREAM < AEvent > *sector,
+ AMI_STREAM < AEvent > *bndSector, long *bndInsert,
+ long longEvents, double start_angle, double end_angle)
+{
+
+
+ unsigned long totalSector = 0, totalDrop = 0, totalInsert = 0;
+
+ for (int i = 0; i < nsect; i++) {
+ assert(total[i] == insert[i] + drop[i]);
+ assert(insert[i] == sector[i].stream_len());
+
+ assert(bndInsert[i] == bndSector[i].stream_len());
+
+ totalSector += total[i];
+ totalDrop += drop[i];
+ totalInsert += insert[i];
+
+ }
+ assert(totalSector == nevents);
+
+#ifdef __GRASS__
+ PRINT_DISTRIBUTE {
+ G_message("-----nsectors=%d\n", nsect);
+ for (int i = 0; i < nsect; i++) {
+ G_message("\ts=%3d ", i);
+ G_message("[%.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_message("Distribute [%.4f, %.4f]: nsect=%d, ", start_angle, end_angle,
+ nsect);
+ G_message
+ ("total events %lu, inserted %lu, dropped %lu, long events=%ld\n",
+ 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;
+}
+
+
+
+/***********************************************************************
+ computes the number of sector for the distribution sweep; technically
+ M/2B because you need 2 streams per sector */
+int compute_n_sectors()
+{
+
+ long memSizeBytes = MM_manager.memory_available();
+ unsigned int blockSizeBytes = UntypedStream::get_block_length();
+
+ /*printf("computeNSect: block=%d, mem=%d\n", blockSizeBytes, (int)memSizeBytes); */
+ int nsect = (int)(memSizeBytes / (2 * blockSizeBytes));
+
+ /*be safe */
+ if (nsect > 4)
+ nsect = nsect / 2;
+
+ /*if it happens that we are at the end of memory, set nsect=2;
+ //technically, if there is not enough memory to hold two
+ //blocks, the function should enter solve_in_memory; so there
+ //is not enough memory to solve in memory nor to
+ //distribute...this shoudl happen only under tests with very
+ //very little memory. just set nsect=2 and hope that it
+ //works */
+
+ if (nsect == 0 || nsect == 1)
+ nsect = 2;
+ else {
+ /*we'll have 2 streams for each sector open; subtract 10 to be safe */
+ if (2 * nsect > MAX_STREAMS_OPEN - 10)
+ nsect = (MAX_STREAMS_OPEN - 10) / 2;
+ }
+ printf("nsectors set to %d\n", nsect);
+ fflush(stdout);
+ return nsect;
+}
+
+
+
+/***********************************************************************
+ compute the sector that contains this angle; there are nsect
+ sectors that span the angle interval [sstartAngle, sendAngle]. if
+ angle falls outside, return -1*/
+int
+get_event_sector(double angle, double sstartAngle, double sendAngle,
+ int nsect)
+{
+
+ int s = -1;
+
+ /*first protect against rounding errors
+ //in the last sector */
+ if (fabs(angle - sendAngle) < EPSILON)
+ return nsect - 1;
+
+ /*in the first sector */
+ if (fabs(angle - sstartAngle) < EPSILON)
+ return 0;
+
+ double ssize = fabs(sstartAngle - sendAngle) / nsect;
+
+ s = (int)((angle - sstartAngle) / ssize);
+ /*printf("getsector: fit %.2f in (%.2f. %.2f)", angle, sstartAngle, sendAngle);
+ //printf("ssize = %.2f, s=%d", ssize, s);
+ //assert (s >= 0 && s < nsect); */
+ if (s < 0 || s >= nsect) {
+ /*falls outside --- this can happen when finding sector of pair
+ //event; e.g. ENTER is inside sector, and its pair EXIT event
+ //falls outside sector. */
+ s = -1;
+ }
+ return s;
+}
+
+
+
+/***********************************************************************
+ insert event in this sector */
+void insert_event_in_sector(AMI_STREAM < AEvent > *str, AEvent * e)
+{
+
+ assert(str && e);
+ AMI_err ae;
+
+ ae = str->write_item(*e);
+ assert(ae == AMI_ERROR_NO_ERROR);
+
+ return;
+}
+
+
+/**********************************************************************
+ insert event e into sector */
+void
+insert_event_in_sector_no_drop(AEvent * e, int s, AMI_STREAM < AEvent > *str,
+ long *insert)
+{
+
+ /*note: if on boundary, PRECISION ISSUES?? should insert both sectors? */
+ DISTRIBDEBUG {
+ print_event(*e);
+ printf(" insert in sector %3d\n", s);
+ }
+ AMI_err ae;
+
+ ae = str->write_item(*e);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ insert[s]++;
+
+ return;
+}
+
+/**********************************************************************
+ insert event e into sector if it is not occluded by high_s */
+void
+insert_event_in_sector(AEvent * e, int s, AMI_STREAM < AEvent > *str,
+ double high_s, Viewpoint * vp, long *insert,
+ long *drop)
+{
+
+ /*note: if on boundary, PRECISION ISSUES?? should insert both sectors?
+
+ //if not occluded by high_s insert it in its sector */
+ if (!is_center_gradient_occluded(e, high_s, vp)) {
+ insert[s]++;
+ DISTRIBDEBUG {
+ print_event(*e);
+ printf(" insert in sector %3d\n", s);
+ }
+ AMI_err ae;
+
+ ae = str->write_item(*e);
+ assert(ae == AMI_ERROR_NO_ERROR);
+
+ }
+ else {
+ /*assert(calculateCenterGradient(e, vp) <= high[s]);
+ //technically, if its a QUERY we should write it as invisible in
+ //vis stream; but as an optimization, our vis stream only
+ //records visible stuff. */
+ DISTRIBDEBUG print_dropped(e, vp, high_s);
+
+ drop[s]++;
+ }
+
+ return;
+}
+
+
+/***********************************************************************
+ returns 1 if the center of event is occluded by the gradient, which
+ is assumed to be in line with the event */
+int is_center_gradient_occluded(AEvent * e, double gradient, Viewpoint * vp)
+{
+
+ assert(e && vp);
+ double eg = calculate_center_gradient(e, vp);
+
+ return (eg < gradient);
+}
+
+/***********************************************************************
+called when dropping an event e, high is the highest gradiant value
+//in its sector*/
+void print_dropped(AEvent * e, Viewpoint * vp, double high)
+{
+
+ assert(e && vp);
+ double eg = calculate_center_gradient(e, vp);
+
+ printf(" dropping grad=%.2f, high=%.2f\n", eg, high);
+
+ return;
+}
Property changes on: grass/trunk/raster/r.viewshed/distribute.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/distribute.h
===================================================================
--- grass/trunk/raster/r.viewshed/distribute.h (rev 0)
+++ grass/trunk/raster/r.viewshed/distribute.h 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,184 @@
+/****************************************************************************
+ *
+ * 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 __DISTRIBUTE_H
+#define __DISTRIBUTE_H
+
+
+#include "visibility.h"
+#include "grid.h"
+#include "eventlist.h"
+
+
+
+/* distribution sweep: write results to visgrid.
+ */
+IOVisibilityGrid *distribute_and_sweep(char* inputfname,
+ GridHeader * hd,
+ Viewpoint * vp,
+ ViewOptions viewOptions);
+
+
+
+
+/* distribute recursively the events and write results to
+ visgrid. eventList is a list of events sorted by distance that fall
+ within angle boundaries start_angle and end_angle; enterBndEvents
+ is a stream that contains all the ENTER events that are not in this
+ sector, but their corresponding Q or X events are in this sector.
+
+ when problem is small enough, solve it in memory and write results
+ to visgrid.
+
+ invariant: distribute_sector deletes eventList and enterBndEvents
+
+ returns the number of visible cells.
+ */
+unsigned long distribute_sector(AMI_STREAM < AEvent > *eventList,
+ AMI_STREAM < AEvent > *enterBndEvents,
+ double start_angle, double end_angle,
+ IOVisibilityGrid * visgrid, Viewpoint * vp,
+ ViewOptions viewOptions);
+
+/* bndEvents is a stream of events that cross into the sector's
+ (first) boundary; they must be distributed to the boundary streams
+ of the sub-sectors of this sector. Note: the boundary streams of
+ the sub-sectors may not be empty; as a result, events get appended
+ at the end, and they will not be sorted by distance from the
+ vp.
+ */
+void distribute_bnd_events(AMI_STREAM < AEvent > *bndEvents,
+ AMI_STREAM < AEvent > *SectorBnd, int nsect,
+ Viewpoint * vp, double start_angle,
+ double end_angle, double *high, long *insert,
+ long *drop);
+
+
+/* same as above, but does it inemory. it is called 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 corresponding EXIT events
+ are inside the sector. start_angle and end_angle are the
+ boundaries of the sector, and visgrid is the grid to which the
+ visible/invisible cells must be written. The sector is solved by
+ switching to radial sweep. Returns the number of visible cells. */
+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);
+
+
+/*returns 1 if enter angle is within epsilon from boundary angle */
+int is_almost_on_boundary(double angle, double boundary_angle);
+
+/* returns 1 if angle is within epsilon the boundaries of sector s */
+int is_almost_on_boundary(double angle, int s, double start_angle,
+ double end_angle, int nsect);
+
+/* computes the number of sector for the distribution sweep;
+ technically M/B */
+int compute_n_sectors();
+
+/* return 1 is s is inside sector; that is, if it is not -1 */
+int is_inside(int s, int nsect);
+
+/* compute the sector that contains this angle; there are nsect
+ sectors that span the angle interval [sstartAngle, sendAngle] */
+int get_event_sector(double angle, double sstartAngle, double sendAngle,
+ int nsect);
+
+
+/* insert event in this sector */
+void insert_event_in_sector(AMI_STREAM < AEvent > *str, AEvent * e);
+
+/* insert event e into sector if it is not occluded by high_s */
+void insert_event_in_sector(AEvent * e, int s, AMI_STREAM < AEvent > *str,
+ double high_s, Viewpoint * vp, long *insert,
+ long *drop);
+
+/**********************************************************************
+ insert event e into sector, no occlusion check */
+void insert_event_in_sector_no_drop(AEvent * e, int s,
+ AMI_STREAM < AEvent > *str, long *insert);
+
+
+
+/* returns 1 if the center of event is occluded by the gradient, which
+ is assumed to be in line with the event */
+int is_center_gradient_occluded(AEvent * e, double gradient, Viewpoint * vp);
+
+/* called when dropping an event e, high is the highest gradiant value
+ in its sector */
+void print_dropped(AEvent * e, Viewpoint * vp, double high);
+
+
+/* prints how many events were inserted and dropped in each sector */
+void print_sector_stats(off_t nevents, int nsect, double *high, long *total,
+ long *insert, long *drop,
+ AMI_STREAM < AEvent > *sector,
+ AMI_STREAM < AEvent > *bndSector, long *bndInsert,
+ long longEvents, double start_angle,
+ double end_angle);
+
+
+/* the event e spans sectors from start_s to end_s; Action: update
+ high[] for each spanned sector.
+ */
+void process_long_cell(int start_s, int end_s, int nsect,
+ Viewpoint * vp, AEvent * e, double *high);
+
+
+/* return the start angle of the i-th sector. Assuming that
+ [start..end] is split into nsectors */
+double get_sector_start(int i, double start_angle, double end_angle,
+ int nsect);
+
+
+/* return the start angle of the i-th sector. Assuming that
+ [start..end] is split into nsectors */
+double get_sector_end(int i, double start_angle, double end_angle, int nsect);
+
+/* returns true if the event is inside the given sector */
+int is_inside(AEvent * e, double start_angle, double end_angle);
+
+/* returns true if this angle is inside the given sector */
+int is_inside(double angle, double start_angle, double end_angle);
+
+
+#endif
Property changes on: grass/trunk/raster/r.viewshed/distribute.h
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/eventlist.cc
===================================================================
--- grass/trunk/raster/r.viewshed/eventlist.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/eventlist.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,961 @@
+/****************************************************************************
+ *
+ * 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 <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+
+#include "eventlist.h"
+
+
+
+/* forced to use this because DistanceCompare::compare has troubles if
+ i put it inside the class */
+Viewpoint globalVP;
+
+
+
+/* ------------------------------------------------------------
+ compute the gradient of the CENTER of this event wrt viewpoint. For
+ efficiency it does not compute the gradient, but the square of the
+ arctan of the gradient. Assuming all gradients are computed the same
+ way, this is correct. */
+double calculate_center_gradient(AEvent * e, Viewpoint * vp)
+{
+
+ assert(e && vp);
+ double gradient, sqdist;
+
+ /*square of the distance from the center of this event to vp */
+ 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;
+ /*maintain sign */
+ if (e->elev < vp->elev)
+ gradient = -gradient;
+ return gradient;
+}
+
+
+
+
+
+/* ------------------------------------------------------------
+ //calculate the angle at which the event is. Return value is the angle.
+
+ angle quadrants:
+ 2 1
+ 3 4
+ ----->x
+ |
+ |
+ |
+ V y
+
+ */
+
+/*/////////////////////////////////////////////////////////////////////
+ //return the angle from this event wrt viewpoint; the type of the
+ //event is taken into position to compute a different amngle for each
+ //event associated with a cell */
+double calculate_event_angle(AEvent * e, Viewpoint * vp)
+{
+
+ assert(e && vp);
+ double ex, ey;
+
+ calculate_event_position(*e, vp->row, vp->col, &ey, &ex);
+ return calculate_angle(ex, ey, vp->col, vp->row);
+}
+
+
+/*/////////////////////////////////////////////////////////////////////
+ //calculate the exit angle corresponding to this cell */
+double
+calculate_exit_angle(dimensionType row, dimensionType col, Viewpoint * vp)
+{
+ AEvent e;
+ double x, y;
+
+ e.eventType = EXITING_EVENT;
+ e.angle = e.elev = 0;
+ e.row = row;
+ e.col = col;
+ calculate_event_position(e, vp->row, vp->col, &y, &x);
+ return calculate_angle(x, y, vp->col, vp->row);
+}
+
+
+/*/////////////////////////////////////////////////////////////////////
+ //calculate the enter angle corresponding to this cell */
+double
+calculate_enter_angle(dimensionType row, dimensionType col, Viewpoint * vp)
+{
+ AEvent e;
+ double x, y;
+
+ e.eventType = ENTERING_EVENT;
+ e.angle = e.elev = 0;
+ e.row = row;
+ e.col = col;
+ calculate_event_position(e, vp->row, vp->col, &y, &x);
+ return calculate_angle(x, y, vp->col, vp->row);
+}
+
+/*///////////////////////////////////////////////////////////////////// */
+double
+calculate_angle(double eventX, double eventY,
+ double viewpointX, double viewpointY)
+{
+
+ /*M_PI is defined in math.h to represent 3.14159... */
+ if (viewpointY == eventY && eventX > viewpointX) {
+ return 0; /*between 1st and 4th quadrant */
+ }
+ else if (eventX > viewpointX && eventY < viewpointY) {
+ /*first quadrant */
+ return atan((viewpointY - eventY) / (eventX - viewpointX));
+
+ }
+ else if (viewpointX == eventX && viewpointY > eventY) {
+ /*between 1st and 2nd quadrant */
+ return (M_PI) / 2;
+
+ }
+ else if (eventX < viewpointX && eventY < viewpointY) {
+ /*second quadrant */
+ return ((M_PI) / 2 +
+ atan((viewpointX - eventX) / (viewpointY - eventY)));
+
+ }
+ else if (viewpointY == eventY && eventX < viewpointX) {
+ /*between 1st and 3rd quadrant */
+ return M_PI;
+
+ }
+ else if (eventY > viewpointY && eventX < viewpointX) {
+ /*3rd quadrant */
+ return (M_PI + atan((eventY - viewpointY) / (viewpointX - eventX)));
+
+ }
+ else if (viewpointX == eventX && viewpointY < eventY) {
+ /*between 3rd and 4th quadrant */
+ return (M_PI * 3.0 / 2.0);
+ }
+ else if (eventX > viewpointX && eventY > viewpointY) {
+ /*4th quadrant */
+ return (M_PI * 3.0 / 2.0 +
+ atan((eventX - viewpointX) / (eventY - viewpointY)));
+ }
+ assert(eventX == viewpointX && eventY == viewpointY);
+ return 0;
+}
+
+
+
+/* ------------------------------------------------------------ */
+/* calculate the exact position of the given event, and store them in x
+ and y.
+ quadrants: 1 2
+ 3 4
+ ----->x
+ |
+ |
+ |
+ V y
+ */
+void
+calculate_event_position(AEvent e, dimensionType viewpointRow,
+ dimensionType viewpointCol, double *y, double *x)
+{
+
+ assert(x && y);
+ *x = 0;
+ *y = 0;
+
+ if (e.eventType == CENTER_EVENT) {
+ /*FOR CENTER_EVENTS */
+ *y = e.row;
+ *x = e.col;
+ return;
+ }
+
+ if (e.row < viewpointRow && e.col < viewpointCol) {
+ /*first quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row - 0.5;
+ *x = e.col + 0.5;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row + 0.5;
+ *x = e.col - 0.5;
+ }
+
+ }
+ 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 + 0.5;
+ *x = e.col + 0.5;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row + 0.5;
+ *x = e.col - 0.5;
+ }
+
+ }
+ else if (e.col > viewpointCol && e.row < viewpointRow) {
+ /*second quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row + 0.5;
+ *x = e.col + 0.5;
+ }
+ else { /*otherwise it is EXITING_EVENT */
+ *y = e.row - 0.5;
+ *x = e.col - 0.5;
+ }
+
+ }
+ 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 + 0.5;
+ *x = e.col - 0.5;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row - 0.5;
+ *x = e.col - 0.5;
+ }
+
+ }
+ else if (e.col > viewpointCol && e.row > viewpointRow) {
+ /*fourth quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row + 0.5;
+ *x = e.col - 0.5;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row - 0.5;
+ *x = e.col + 0.5;
+ }
+
+ }
+ 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 - 0.5;
+ *x = e.col - 0.5;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row - 0.5;
+ *x = e.col + 0.5;
+ }
+
+ }
+ else if (e.col < viewpointCol && e.row > viewpointRow) {
+ /*third quadrant */
+ if (e.eventType == ENTERING_EVENT) {
+ /*if it is ENTERING_EVENT */
+ *y = e.row - 0.5;
+ *x = e.col - 0.5;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row + 0.5;
+ *x = e.col + 0.5;
+ }
+
+ }
+ 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 - 0.5;
+ *x = e.col + 0.5;
+ }
+ else {
+ /*otherwise it is EXITING_EVENT */
+ *y = e.row + 0.5;
+ *x = e.col + 0.5;
+ }
+ }
+ else {
+ /*must be the 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 ((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);
+ // } */
+ return;
+}
+
+
+/* ------------------------------------------------------------ */
+void print_event(AEvent a)
+{
+ char c = '0';
+
+ if (a.eventType == ENTERING_EVENT)
+ c = 'E';
+ if (a.eventType == EXITING_EVENT)
+ 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);
+ return;
+}
+
+
+/* ------------------------------------------------------------ */
+/*computes the distance from the event to the viewpoint. Note: all 3
+ //events associate to a cell are considered at the same distance, from
+ //the center of the cell to the viewpoint */
+double
+get_square_distance_from_viewpoint(const AEvent & a, const Viewpoint & vp)
+{
+
+ double eventy, eventx;
+
+ calculate_event_position(a, vp.row, vp.col, &eventy, &eventx);
+
+ double dist = (eventx - vp.col) * (eventx - vp.col) +
+ (eventy - vp.row) * (eventy - vp.row);
+ /*don't take sqrt, it is expensive; suffices for comparison */
+ return dist;
+}
+
+/* ------------------------------------------------------------ */
+/* a duplicate of get_square_distance_from_viewpoint() needed for debug */
+double
+get_square_distance_from_viewpoint_with_print(const AEvent & a,
+ const Viewpoint & vp)
+{
+
+ double eventy, eventx;
+
+ calculate_event_position(a, vp.row, vp.col, &eventy, &eventx);
+ double dist = (eventx - vp.col) * (eventx - vp.col) +
+ (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);
+
+ return dist;
+}
+
+
+/* ------------------------------------------------------------ */
+/*determines if the point at row,col is outside the maximum distance
+ limit. Return 1 if the point is outside limit, 0 if point is inside
+ limit. */
+int is_point_outside_max_dist(Viewpoint vp, GridHeader hd,
+ dimensionType row, dimensionType col,
+ float maxDist)
+{
+ /* it is not too smart to compare floats */
+ if ((int)maxDist == INFINITY_DISTANCE)
+ return 0;
+
+ 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;
+ }
+}
+
+
+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 */
+int DistanceCompare::compare(const AEvent & a, const AEvent & b)
+{
+
+ /*calculate distance from viewpoint
+ //don't take sqrt, it is expensive; suffices for comparison */
+ double da, db;
+
+ /*da = get_square_distance_from_viewpoint(a, globalVP);
+ //db = get_square_distance_from_viewpoint(b, globalVP); */
+
+ /*in the event these are not inlined */
+ double eventy, eventx;
+
+ calculate_event_position(a, globalVP.row, globalVP.col, &eventy, &eventx);
+ da = (eventx - globalVP.col) * (eventx - globalVP.col) +
+ (eventy - globalVP.row) * (eventy - globalVP.row);
+ calculate_event_position(b, globalVP.row, globalVP.col, &eventy, &eventx);
+ db = (eventx - globalVP.col) * (eventx - globalVP.col) +
+ (eventy - globalVP.row) * (eventy - globalVP.row);
+
+ if (da > db) {
+ return 1;
+ }
+ else if (da < db) {
+ return -1;
+ }
+ else {
+ return 0;
+ }
+ return 0;
+}
+
+
+/* ------------------------------------------------------------ */
+int RadialCompare::compare(const AEvent & a, const AEvent & b)
+{
+
+ if (a.row == b.row && a.col == b.col && a.eventType == b.eventType)
+ return 0;
+
+ assert(a.angle >= 0 && b.angle >= 0);
+
+ if (a.angle > b.angle) {
+ return 1;
+ }
+ else if (a.angle < b.angle) {
+ return -1;
+ }
+ else {
+ /*a.angle == b.angle */
+ if (a.eventType == EXITING_EVENT)
+ return -1;
+ else if (a.eventType == ENTERING_EVENT)
+ return 1;
+ return 0;
+ }
+}
+
+/* ------------------------------------------------------------ */
+/* a copy of the function above is needed by qsort, when teh
+ computation runs in memory */
+
+int radial_compare_events(const void *x, const void *y)
+{
+
+ AEvent *a, *b;
+
+ a = (AEvent *) x;
+ b = (AEvent *) y;
+ if (a->row == b->row && a->col == b->col && a->eventType == b->eventType)
+ return 0;
+
+ assert(a->angle >= 0 && b->angle >= 0);
+
+ if (a->angle > b->angle) {
+ return 1;
+ }
+ else if (a->angle < b->angle) {
+ return -1;
+ }
+ else {
+ /*a->angle == b->angle */
+ if (a->eventType == EXITING_EVENT)
+ return -1;
+ else if (a->eventType == ENTERING_EVENT)
+ return 1;
+ return 0;
+ }
+}
+
+
+
+/* ------------------------------------------------------------ */
+/*sort the event list in radial order */
+void sort_event_list(AMI_STREAM < AEvent > **eventList)
+{
+
+ /*printf("sorting events.."); fflush(stdout); */
+ assert(*eventList);
+
+ AMI_STREAM < AEvent > *sortedStr;
+ RadialCompare cmpObj;
+ AMI_err ae;
+
+ ae = AMI_sort(*eventList, &sortedStr, &cmpObj, 1);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ *eventList = sortedStr;
+ /*printf("..done.\n"); fflush(stdout); */
+ return;
+}
+
+
+/* ------------------------------------------------------------ */
+/*sort the event list in distance order */
+void
+sort_event_list_by_distance(AMI_STREAM < AEvent > **eventList, Viewpoint vp)
+{
+
+ /*printf("sorting events by distance from viewpoint.."); fflush(stdout); */
+ assert(*eventList);
+
+ AMI_STREAM < AEvent > *sortedStr;
+ DistanceCompare cmpObj;
+
+ globalVP.row = vp.row;
+ globalVP.col = vp.col;
+ /*printViewpoint(globalVP); */
+ AMI_err ae;
+
+ ae = AMI_sort(*eventList, &sortedStr, &cmpObj, 1);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ *eventList = sortedStr;
+ /*printf("..sorting done.\n"); fflush(stdout); */
+ 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;
+}
Property changes on: grass/trunk/raster/r.viewshed/eventlist.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/eventlist.h
===================================================================
--- grass/trunk/raster/r.viewshed/eventlist.h (rev 0)
+++ grass/trunk/raster/r.viewshed/eventlist.h 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,170 @@
+
+/****************************************************************************
+ *
+ * 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 _EVENTLIST_H
+#define _EVENTLIST_H
+
+#include "grid.h"
+#include "visibility.h"
+
+#ifdef __GRASS__
+#include <grass/iostream/ami.h>
+#else
+#include <ami.h>
+#endif
+
+
+#define ENTERING_EVENT 1
+#define EXITING_EVENT -1
+#define CENTER_EVENT 0
+
+typedef struct event_
+{
+ dimensionType row, col; //location of the center of cell
+ float elev; //elevation here
+ double angle;
+ char eventType;
+
+ //type of the event: ENTERING_EVENT, EXITING_EVENT, CENTER_EVENT
+} AEvent;
+
+
+
+/* ------------------------------------------------------------ */
+/*determines if the point at row,col is outside the maximum distance
+ limit wrt viewpoint. Return 1 if the point is outside
+ limit, 0 if point is inside limit. */
+int
+is_point_outside_max_dist(Viewpoint vp, GridHeader hd,
+ dimensionType row, dimensionType col,
+ float maxDist);
+
+
+
+/* ------------------------------------------------------------ */
+/* 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);
+
+class RadialCompare
+{
+ public:int compare(const AEvent &, const AEvent &);
+};
+int radial_compare_events(const void *a, const void *b);
+
+
+ /*sort the event list by the distance from the viewpoint */
+class DistanceCompare
+{
+ public:int compare(const AEvent &, const AEvent &);
+};
+void print_event(AEvent a);
+
+
+ /*computes the distance from the event to the viewpoint. Note: all 3
+ //events associate to a cell are considered at the same distance, from
+ //the center of teh cell to the viewpoint */
+double get_square_distance_from_viewpoint(const AEvent & a,
+ const Viewpoint & vp);
+
+ /*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
+ //event is taken into position to compute a different amngle for each
+ //event associated with a cell */
+double calculate_event_angle(AEvent * e, Viewpoint * vp);
+
+
+ /*compute the gradient of the CENTER of this event wrt viewpoint. For
+ //efficiency it does not compute the gradient, but the square of the
+ //tan of the gradient. Assuming all gradients are computed the same
+ //way, this is correct. */
+double calculate_center_gradient(AEvent * e, Viewpoint * vp);
+
+
+ /*calculate the exit angle corresponding to this cell */
+double calculate_exit_angle(dimensionType row, dimensionType col,
+ Viewpoint * vp);
+
+ /*calculate the enter angle corresponding to this cell */
+double calculate_enter_angle(dimensionType row, dimensionType col,
+ Viewpoint * vp);
+
+ /*calculate the exact position of the given event, and store them in x
+ //and y. */
+void calculate_event_position(AEvent e, dimensionType viewpointRow,
+ dimensionType viewpointCol, double *y,
+ double *x);
+double calculate_angle(double eventX, double eventY, double viewpointX,
+ double viewpointY);
+
+#endif
Property changes on: grass/trunk/raster/r.viewshed/eventlist.h
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/grass.cc
===================================================================
--- grass/trunk/raster/r.viewshed/grass.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/grass.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,822 @@
+/****************************************************************************
+ *
+ * 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 <stdlib.h>
+#include <stdio.h>
+
+#ifdef __GRASS__
+
+extern "C"
+{
+#include <grass/config.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+}
+
+#include "grass.h"
+#include "visibility.h"
+
+
+
+
+
+/* ------------------------------------------------------------ */
+/* if viewOptions.doCurv is on then adjust the passed height for
+ curvature of the earth; otherwise return the passed height
+ unchanged.
+*/
+float adjust_for_curvature(Viewpoint vp, dimensionType row,
+ dimensionType col, float h,
+ ViewOptions viewOptions) {
+
+ if (!viewOptions.doCurv) return h;
+
+ assert(viewOptions.ellps_a != 0);
+ double dif_x, dif_y, sqdist;
+ dif_x = (vp.col - col);
+ dif_y = (vp.row - row);
+ sqdist = (dif_x * dif_x + dif_y * dif_y) * viewOptions.cellsize * viewOptions.cellsize;
+
+ return h - (sqdist / (2 * viewOptions.ellps_a));
+}
+
+
+
+
+
+/* ************************************************************ */
+/*return a GridHeader that has all the relevant data filled in from
+ GRASS */
+GridHeader *read_header_from_GRASS(char *rastName, Cell_head *region) {
+
+ assert(rastName);
+
+ /*allocate the grid header to fill */
+ GridHeader *hd = (GridHeader *) G_malloc(sizeof(GridHeader));
+ assert(hd);
+
+ /*get the num rows and cols from GRASS */
+ int nrows, ncols;
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ /*check for loss of prescion */
+ if (nrows <= maxDimension && ncols <= maxDimension) {
+ hd->nrows = (dimensionType) nrows;
+ hd->ncols = (dimensionType) ncols;
+ }
+ else {
+ G_fatal_error(_("grid dimension too big for current precision"));
+ exit(EXIT_FAILURE);
+ }
+
+ /*fill in rest of header */
+ hd->xllcorner = G_col_to_easting(0, region);
+ hd->yllcorner = G_row_to_northing(0, region);
+ /*Cell_head stores 2 resolutions, while GridHeader only stores 1
+ //make sure the two Cell_head resolutions are equal */
+ if (fabs(region->ew_res - region->ns_res) > .001) {
+ G_warning(_("east-west resolution does not equal north-south resolutio. The viewshed computation assumes the cells are square, so in this case this may result in innacuracies."));
+ // exit(EXIT_FAILURE);
+ //
+ }
+ hd->cellsize = (float) region->ew_res;
+ //store the null value of the map
+ G_set_f_null_value(& (hd->nodata_value), 1);
+ printf("nodata value set to %f\n", hd->nodata_value);
+ return hd;
+}
+
+
+
+
+
+
+
+
+/* ************************************************************ */
+/* input: an array capable to hold the max number of events, a raster
+ name, a viewpoint and the viewOptions; 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. initialize and fill
+ 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,
+ Viewpoint * vp, GridHeader* hd,
+ ViewOptions viewOptions, double **data,
+ MemoryVisibilityGrid* visgrid ) {
+
+ G_verbose_message(_("computing events ..."));
+ assert(eventList && vp && visgrid);
+ //GRASS should be defined
+
+ /*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));
+ assert(*data);
+
+ /*get the mapset name */
+ char *mapset;
+ mapset = G_find_cell(rastName, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("raster map [%s] not found"), rastName);
+ exit(EXIT_FAILURE);
+ }
+
+ /*open map */
+ int infd;
+ if ((infd = G_open_cell_old(rastName, mapset)) < 0) {
+ G_fatal_error(_("Cannot open raster file [%s]"), rastName);
+ exit(EXIT_FAILURE);
+ }
+
+ /*get the data_type */
+ RASTER_MAP_TYPE data_type;
+ data_type = G_raster_map_type(rastName, mapset);
+
+ /*buffer to hold a row */
+ void *inrast;
+ inrast = G_allocate_raster_buf(data_type);
+ assert(inrast);
+
+ /*DCELL to check for loss of prescion- haven't gotten that to work
+ yet though */
+ DCELL d;
+ int isnull = 0;
+
+ /*keep track of the number of events added, to be returned later */
+ size_t nevents = 0;
+
+ /*scan through the raster data */
+ dimensionType i, j, k;
+ double ax, ay;
+ AEvent e;
+
+ e.angle = -1;
+ for (i = 0; i < G_window_rows(); i++) {
+ /*read in the raster row */
+ int rasterRowResult = G_get_raster_row(infd, inrast, i, data_type);
+ if (rasterRowResult <= 0) {
+ G_fatal_error(_("coord not read from row %d of %s"), i, rastName);
+ exit(EXIT_FAILURE);
+ }
+
+ /*fill event list with events from this row */
+ for (j = 0; j < G_window_cols(); 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;
+ }
+
+ /* adjust for curvature */
+ e.elev = adjust_for_curvature(*vp, i, j, e.elev, viewOptions);
+
+ /*write it into the row of data going through the viewpoint */
+ if (i == vp->row) {
+ (*data)[j] = e.elev;
+ }
+
+ /* 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);
+ if (isnull) {
+ /*what to do when viewpoint is NODATA ? */
+ G_message(_("WARNING: viewpoint is NODATA. "));
+ G_message(_("Will assume its elevation is = %f"), vp->elev);
+ }
+ continue;
+ }
+
+ /*don't insert in eventlist nodata cell events */
+ if (isnull) {
+ /* 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;
+
+ /* if it got here it is not the viewpoint, not NODATA, and
+ within max distance from viewpoint; generate its 3 events
+ and insert them */
+ /*put event into event list */
+ 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++;
+
+ }
+ }
+
+ G_message(_("...done creating event list"));
+ G_close_cell(infd);
+ 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 > *
+grass_init_event_list(char *rastName, Viewpoint* vp, GridHeader* hd,
+ ViewOptions viewOptions, double **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));
+ assert(*data);
+ }
+
+ /*create the event stream that will hold the events */
+ AMI_STREAM < AEvent > *eventList = new AMI_STREAM < AEvent > ();
+ assert(eventList);
+
+ /*determine which mapset we are in */
+ char *mapset;
+ mapset = G_find_cell(rastName, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("raster map [%s] not found"), rastName);
+ exit(EXIT_FAILURE);
+ }
+
+ /*open map */
+ int infd;
+ if ((infd = G_open_cell_old(rastName, mapset)) < 0) {
+ G_fatal_error(_("Cannot open raster file [%s]"), rastName);
+ exit(EXIT_FAILURE);
+ }
+ RASTER_MAP_TYPE data_type;
+ data_type = G_raster_map_type(rastName, mapset);
+ void *inrast;
+ inrast = G_allocate_raster_buf(data_type);
+ assert(inrast);
+
+ /*scan through the raster data */
+ DCELL d;
+ int isnull = 0;
+ dimensionType i, j, k;
+ double ax, ay;
+ AEvent e;
+ e.angle = -1;
+
+ /*start scanning through the grid */
+ for (i = 0; i < G_window_rows(); i++) {
+
+ /*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);
+ exit(EXIT_FAILURE);
+ }
+ /*fill event list with events from this row */
+ for (j = 0; j < G_window_cols(); 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;
+ }
+
+ /* adjust for curvature */
+ e.elev = adjust_for_curvature(*vp, i, j, e.elev, viewOptions);
+
+ if (data!= NULL) {
+ /**write the row of data going through the viewpoint */
+ if (i == vp->row) {
+ (*data)[j] = e.elev;
+ }
+ }
+
+ /* set the viewpoint */
+ 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 viewpoint into eventlist */
+ 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;
+
+ /*put event into event list */
+ 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 */
+
+ G_message(_("...done creating event list\n"));
+ G_close_cell(infd);
+ 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;
+}
+
+
+
+
+
+
+/* ************************************************************ */
+/* saves the grid into a GRASS raster. Loops through all elements x
+ in row-column order and writes fun(x) to file.*/
+void
+save_grid_to_GRASS(Grid* grid, char* filename, RASTER_MAP_TYPE type,
+ float(*fun)(float)) {
+
+ G_message(_("saving grid to %s"), filename);
+ assert(grid && filename);
+
+ /*open the new raster */
+ int outfd;
+ outfd = G_open_raster_new(filename, type);
+
+ /*get the buffer to store values to read and write to each row */
+ void *outrast;
+ outrast = G_allocate_raster_buf(type);
+ assert(outrast);
+
+ dimensionType i, j;
+ for (i = 0; i < G_window_rows(); i++) {
+ for (j = 0; j < G_window_cols(); j++) {
+
+ switch (type) {
+ case CELL_TYPE:
+ ((CELL *) outrast)[j] = (CELL) fun(grid->grid_data[i][j]);
+ break;
+ case FCELL_TYPE:
+ ((FCELL *) outrast)[j] = (FCELL) fun(grid->grid_data[i][j]);
+ break;
+ case DCELL_TYPE:
+ ((DCELL *) outrast)[j] = (DCELL) fun(grid->grid_data[i][j]);
+ break;
+ }
+ } /* for j */
+ G_put_raster_row(outfd, outrast, type);
+ } /* for i */
+
+ G_close_cell(outfd);
+ return;
+}
+
+
+
+
+
+/* ************************************************************ */
+/* using the visibility information recorded in visgrid, it creates an
+ output viewshed raster with name outfname; for every point p that
+ is visible in the grid, the corresponding value in the output
+ raster is elevation(p) - viewpoint_elevation(p); the elevation
+ values are read from elevfname raster */
+
+void
+save_vis_elev_to_GRASS(Grid* visgrid, char* elevfname, char* visfname,
+ float vp_elev) {
+
+ G_message(_("saving grid to %s"), visfname);
+ assert(visgrid && elevfname && visfname);
+
+ /*get the mapset name */
+ char *mapset;
+ mapset = G_find_cell(elevfname, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("raster map [%s] not found"), elevfname);
+ exit(EXIT_FAILURE);
+ }
+ /*open elevation map */
+ int elevfd;
+ if ((elevfd = G_open_cell_old(elevfname, mapset)) < 0) {
+ G_fatal_error(_("Cannot open raster file [%s]"), elevfname);
+ exit(EXIT_FAILURE);
+ }
+
+ /*get elevation data_type */
+ RASTER_MAP_TYPE elev_data_type;
+ elev_data_type = G_raster_map_type(elevfname, mapset);
+
+ /* create the visibility raster of same type */
+ int visfd;
+ visfd = G_open_raster_new(visfname, elev_data_type);
+
+ /*get the buffers to store each row */
+ void *elevrast;
+ elevrast = G_allocate_raster_buf(elev_data_type);
+ assert(elevrast);
+ void *visrast;
+ visrast = G_allocate_raster_buf(elev_data_type);
+ assert(visrast);
+
+ dimensionType i, j;
+ double elev=0, viewshed_value;
+ for (i = 0; i < G_window_rows(); i++) {
+ /* get the row from elevation */
+ if (G_get_raster_row(elevfd, elevrast, i, elev_data_type) <= 0) {
+ G_fatal_error(_("save_vis_elev_to_GRASS: could not read row %d"), i);
+ }
+ for (j = 0; j < G_window_cols(); j++) {
+
+ /* read the current elevation value */
+ int isNull = 0;
+ switch (elev_data_type) {
+ case CELL_TYPE:
+ isNull = G_is_c_null_value(&((CELL *) elevrast)[j]);
+ elev = (double)(((CELL *) elevrast)[j]);
+ break;
+ case FCELL_TYPE:
+ isNull = G_is_f_null_value(&((FCELL *) elevrast)[j]);
+ elev = (double)(((FCELL *) elevrast)[j]);
+ break;
+ case DCELL_TYPE:
+ isNull = G_is_d_null_value(&((DCELL *) elevrast)[j]);
+ elev = (double)(((DCELL *) elevrast)[j]);
+ break;
+ }
+
+ if (is_visible(visgrid->grid_data[i][j])) {
+ /* elevation cannot be null */
+ assert(!isNull);
+ /* write elev - viewpoint_elevation*/
+ viewshed_value = elev - vp_elev;
+ writeValue(visrast,j,viewshed_value, elev_data_type);
+ } else if (is_invisible_not_nodata(visgrid->grid_data[i][j])) {
+ /* elevation cannot be null */
+ assert(!isNull);
+ /* write INVISIBLE */
+ viewshed_value = INVISIBLE;
+ writeValue(visrast,j,viewshed_value, elev_data_type);
+ } else {
+ /* nodata */
+ assert(isNull);
+ /* write NODATA */
+ writeNodataValue(visrast, j, elev_data_type);
+ }
+
+
+ } /* for j*/
+ G_put_raster_row(visfd, visrast, elev_data_type);
+ } /* for i */
+
+ G_close_cell(elevfd);
+ G_close_cell(visfd);
+ return;
+}
+
+
+
+
+/* helper function to deal with GRASS writing to a row buffer */
+void writeValue(void* bufrast, int j, double x, RASTER_MAP_TYPE data_type) {
+
+ switch (data_type) {
+ case CELL_TYPE:
+ ((CELL *) bufrast)[j] = (CELL) x;
+ break;
+ case FCELL_TYPE:
+ ((FCELL *) bufrast)[j] = (FCELL) x;
+ break;
+ case DCELL_TYPE:
+ ((DCELL *) bufrast)[j] = (DCELL) x;
+ break;
+ default:
+ G_fatal_error(_("unknown type"));
+ }
+}
+
+void writeNodataValue(void* bufrast, int j, RASTER_MAP_TYPE data_type) {
+
+ switch (data_type) {
+ case CELL_TYPE:
+ G_set_c_null_value(& ((CELL *) bufrast)[j], 1);
+ break;
+ case FCELL_TYPE:
+ G_set_f_null_value(& ((FCELL *) bufrast)[j], 1);
+ break;
+ case DCELL_TYPE:
+ G_set_d_null_value(& ((DCELL *) bufrast)[j], 1);
+ break;
+ default:
+ G_fatal_error(_("unknown type"));
+ exit(EXIT_FAILURE);
+ }
+}
+
+
+/* ************************************************************ */
+/* write the visibility grid to GRASS. 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_GRASS(IOVisibilityGrid * visgrid,
+ char* fname, RASTER_MAP_TYPE type,
+ float (*fun)(float)) {
+
+ G_message(_("saving grid to %s"), fname);
+ assert(fname && visgrid);
+
+ /* open the output raster and set up its row buffer */
+ int visfd;
+ visfd = G_open_raster_new(fname, type);
+ void * visrast;
+ visrast = G_allocate_raster_buf(type);
+ assert(visrast);
+
+ /*set up reading data from 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++;
+ }
+
+ dimensionType i, j;
+ for (i = 0; i < G_window_rows(); i++) {
+ for (j = 0; j < G_window_cols(); j++) {
+
+ if (curResult->row == i && curResult->col == j) {
+ /*cell is recodred in the visibility stream: it must be
+ either visible, or NODATA */
+ writeValue(visrast,j, fun(curResult->angle), type);
+
+ /*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, so it is invisible */
+ writeValue(visrast, j, fun(INVISIBLE), type);
+ }
+ } /* for j */
+
+ G_put_raster_row(visfd, visrast, type);
+ } /* for i */
+
+ G_close_cell(visfd);
+}
+
+
+
+
+
+
+
+/* ************************************************************ */
+/* using the visibility information recorded in visgrid, it creates
+ an output viewshed raster with name outfname; for every point p
+ that is visible in the grid, the corresponding value in the output
+ raster is elevation(p) - viewpoint_elevation(p); the elevation
+ values are read from elevfname raster. assume stream is sorted in
+ (i,j) order. */
+void
+save_io_vis_and_elev_to_GRASS(IOVisibilityGrid * visgrid, char* elevfname,
+ char* visfname, float vp_elev) {
+
+ G_message(_("saving grid to %s"), visfname);
+ assert(visfname && visgrid);
+
+ /*get mapset name and data type */
+ char *mapset;
+ mapset = G_find_cell(elevfname, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("opening %s: cannot find raster"), elevfname);
+ exit(EXIT_FAILURE);
+ }
+
+ /*open elevation map */
+ int elevfd;
+ if ((elevfd = G_open_cell_old(elevfname, mapset)) < 0) {
+ G_fatal_error(_("Cannot open raster file [%s]"), elevfname);
+ exit(EXIT_FAILURE);
+ }
+
+ /*get elevation data_type */
+ RASTER_MAP_TYPE elev_data_type;
+ elev_data_type = G_raster_map_type(elevfname, mapset);
+
+ /* open visibility raster */
+ int visfd;
+ visfd = G_open_raster_new(visfname, elev_data_type);
+
+ /*get the buffers to store each row */
+ void *elevrast;
+ elevrast = G_allocate_raster_buf(elev_data_type);
+ assert(elevrast);
+ void *visrast;
+ visrast = G_allocate_raster_buf(elev_data_type);
+ assert(visrast);
+
+ /*set up stream reading stuff */
+ off_t streamLen, counter = 0;
+ AMI_err ae;
+ VisCell *curResult = NULL;
+
+ AMI_STREAM < VisCell > *vstr = visgrid->visStr;
+ streamLen = vstr->stream_len();
+ vstr->seek(0);
+
+ /*read the first element */
+ if (streamLen > 0) {
+ ae = vstr->read_item(&curResult);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ counter++;
+ }
+
+ dimensionType i, j;
+ double elev=0, viewshed_value;
+ for (i = 0; i < G_window_rows(); i++) {
+
+ if (G_get_raster_row(elevfd, elevrast, i, elev_data_type) <= 0) {
+ G_fatal_error(_("could not read row %d"), i);
+ exit(EXIT_FAILURE);
+ }
+
+ for (j = 0; j < G_window_cols(); j++) {
+
+ /* read the current elevation value */
+ int isNull = 0;
+ switch (elev_data_type) {
+ case CELL_TYPE:
+ isNull = G_is_c_null_value(&((CELL *) elevrast)[j]);
+ elev = (double)(((CELL *) elevrast)[j]);
+ break;
+ case FCELL_TYPE:
+ isNull = G_is_f_null_value(&((FCELL *) elevrast)[j]);
+ elev = (double)(((FCELL *) elevrast)[j]);
+ break;
+ case DCELL_TYPE:
+ isNull = G_is_d_null_value(&((DCELL *) elevrast)[j]);
+ elev = (double)(((DCELL *) elevrast)[j]);
+ break;
+ }
+
+
+ if (curResult->row == i && curResult->col == j) {
+ /*cell is recodred in the visibility stream: it must be
+ either visible, or NODATA */
+ if (is_visible(curResult->angle))
+ writeValue(visrast,j, elev -vp_elev, elev_data_type);
+ else
+ writeNodataValue(visrast,j, elev_data_type);
+ /*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, so it is invisible */
+ writeValue(visrast, j, INVISIBLE, elev_data_type);
+ }
+ } /* for j */
+
+ G_put_raster_row(visfd, visrast, elev_data_type);
+ } /* for i */
+
+ G_close_cell(elevfd);
+ G_close_cell(visfd);
+ return;
+}
+
+
+
+#endif
Property changes on: grass/trunk/raster/r.viewshed/grass.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/grass.h
===================================================================
--- grass/trunk/raster/r.viewshed/grass.h (rev 0)
+++ grass/trunk/raster/r.viewshed/grass.h 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,151 @@
+/****************************************************************************
+ *
+ * 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.
+ *
+ *****************************************************************************/
+
+
+#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>
+}
+
+#include "eventlist.h"
+#include "grid.h"
+#include "visibility.h"
+
+
+/* ------------------------------------------------------------ */
+/* if viewOptions.doCurv is on then adjust the passed height for
+ curvature of the earth; otherwise return the passed height
+ unchanged.
+*/
+float adjust_for_curvature(Viewpoint vp, dimensionType row,
+ dimensionType col, float h,
+ ViewOptions viewOptions);
+
+
+/* helper function to deal with GRASS writing to a row buffer */
+void writeValue(void* ptr, int j, double x, RASTER_MAP_TYPE data_type);
+void writeNodataValue(void* ptr, int j, RASTER_MAP_TYPE data_type);
+
+
+
+/*return a GridHeader with all the relevant data filled in from GRASS */
+GridHeader *read_header_from_GRASS(char *rastName, Cell_head * region);
+
+/* ************************************************************ */
+/* input: an array capable to hold the max number of events, a raster
+ name, a viewpoint and the viewOptions; 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. initialize and fill
+ 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,
+ Viewpoint * vp, GridHeader* hd,
+ ViewOptions viewOptions, double **data,
+ MemoryVisibilityGrid* visgrid );
+
+
+
+/* ************************************************************ */
+/* 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 > *
+grass_init_event_list(char *rastName, Viewpoint* vp, GridHeader* hd,
+ ViewOptions viewOptions, double **data,
+ IOVisibilityGrid* visgrid);
+
+
+/* ************************************************************ */
+/* saves the grid into a GRASS raster. Loops through all elements x
+ in row-column order and writes fun(x) to file.*/
+void
+save_grid_to_GRASS(Grid* grid, char* filename, RASTER_MAP_TYPE type,
+ float(*fun)(float));
+
+
+/* ************************************************************ */
+/* using the visibility information recorded in visgrid, it creates an
+ output viewshed raster with name outfname; for every point p that
+ is visible in the grid, the corresponding value in the output
+ raster is elevation(p) - viewpoint_elevation(p); the elevation
+ values are read from elevfname raster */
+
+void
+save_vis_elev_to_GRASS(Grid* visgrid, char* elevfname, char* visfname,
+ float vp_elev);
+
+
+/* ************************************************************ */
+/* write the visibility grid to GRASS. 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_GRASS(IOVisibilityGrid * visgrid,
+ char* outfname, RASTER_MAP_TYPE type,
+ float (*fun)(float));
+
+
+
+/* ************************************************************ */
+/* using the visibility information recorded in visgrid, it creates
+ an output viewshed raster with name outfname; for every point p
+ that is visible in the grid, the corresponding value in the output
+ raster is elevation(p) - viewpoint_elevation(p); the elevation
+ values are read from elevfname raster. assume stream is sorted in
+ (i,j) order. */
+void
+save_io_vis_and_elev_to_GRASS(IOVisibilityGrid * visgrid, char* elevfname,
+ char* visfname, float vp_elev);
+
+#endif/*_GRASS_H*/
+
+#endif /*__GRASS__*/
Property changes on: grass/trunk/raster/r.viewshed/grass.h
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/grid.cc
===================================================================
--- grass/trunk/raster/r.viewshed/grid.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/grid.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,397 @@
+/****************************************************************************
+ *
+ * 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 <stdio.h>
+#include <math.h>
+#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)
+{
+ assert(a);
+ a->nrows = b.nrows;
+ a->ncols = b.ncols;
+ a->xllcorner = b.xllcorner;
+ a->yllcorner = b.yllcorner;
+ a->cellsize = b.cellsize;
+ a->nodata_value = b.nodata_value;
+ return;
+}
+
+
+/* ------------------------------------------------------------ */
+/*print header */
+void print_grid_header(GridHeader * hd)
+{
+
+ assert(hd);
+ print_grid_header(stdout, hd);
+ return;
+}
+
+
+/* ------------------------------------------------------------ */
+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);
+ return is_nodata(grid->hd, value);
+}
+
+
+
+/* ------------------------------------------------------------ */
+/* create an empty grid and return it. The header and the data are set
+ to NULL. */
+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);
+
+ /*initialize structure */
+ ptr_grid->hd = NULL;
+ ptr_grid->grid_data = NULL;
+
+#ifdef _DEBUG_ON
+ printf("**DEBUG: createEmptyGrid \n");
+ fflush(stdout);
+#endif
+
+ return ptr_grid;
+}
+
+
+
+
+/* ------------------------------------------------------------ */
+/* allocate memroy for grid_data; grid must have a header that gives
+ the dimensions */
+void alloc_grid_data(Grid * pgrid)
+{
+ 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]);
+ }
+
+#ifdef _DEBUG_ON
+ printf("**DEBUG: allocGridData\n");
+ fflush(stdout);
+#endif
+
+ 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)
+{
+ assert(grid);
+
+ /*free grid data if its allocated */
+#ifdef __GRASS__
+ if (grid->grid_data) {
+ dimensionType i;
+
+ for (i = 0; i < grid->hd->nrows; i++) {
+ if (!grid->grid_data[i])
+ G_free((float *)grid->grid_data[i]);
+ }
+
+ G_free((float **)grid->grid_data);
+ }
+
+ 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);
+#endif
+
+ 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;
+}
Property changes on: grass/trunk/raster/r.viewshed/grid.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/grid.h
===================================================================
--- grass/trunk/raster/r.viewshed/grid.h (rev 0)
+++ grass/trunk/raster/r.viewshed/grid.h 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,119 @@
+/****************************************************************************
+ *
+ * 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.
+ *
+ *****************************************************************************/
+
+/*
+ A grid in ArcInfo Ascii Grid Format
+ */
+
+
+#ifndef __GRID_H
+#define __GRID_H
+
+#include <stdio.h>
+#include <limits.h>
+
+
+
+/* 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;
+static const dimensionType maxDimension = USHRT_MAX - 1;
+
+
+typedef struct grid_header {
+ 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 */
+} GridHeader;
+
+
+
+typedef struct grid_ {
+ GridHeader *hd;
+
+ /*two dimensional array holding all the values in the grid */
+ float **grid_data;
+
+ float minvalue; /*the minimum value in the grid */
+ float maxvalue; /*the maximum value in the grid */
+} Grid;
+
+
+
+
+/* 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);
+
+/* create and return an empty grid */
+Grid *create_empty_grid();
+
+/*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
Property changes on: grass/trunk/raster/r.viewshed/grid.h
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/main.cc
===================================================================
--- grass/trunk/raster/r.viewshed/main.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/main.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,863 @@
+
+/****************************************************************************
+ *
+ * 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 <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <ctype.h>
+#include <unistd.h>
+
+#ifdef __GRASS__
+extern "C"
+{
+#include <grass/config.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+}
+#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"
+
+
+
+
+/* if the user does not specify how much memory is available for the
+ program, this is the default value used (in bytes) */
+#define DEFAULT_MEMORY 500<<20
+
+
+/* observer elevation above the terrain */
+#define DEFAULT_OBS_ELEVATION 0
+
+
+/* All these flags are used for debugging */
+
+/* if this flag is set, it always runs in memory */
+//#define FORCE_INTERNAL
+
+/* if this is set, it runs in external memory, even if the problem is
+ small enough to fit in memory. In external memory it first tries
+ to run the base-case, then recursion. */
+//#define FORCE_EXTERNAL
+
+/* if this flag is set it runs in external memory, and starts
+ recursion without checking the base-case. */
+//#define FORCE_DISTRIBUTION
+
+
+
+
+
+/* ------------------------------------------------------------ */
+/* forward declarations */
+/* ------------------------------------------------------------ */
+void print_timings_internal(Rtimer sweepTime, Rtimer outputTime, 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
+
+
+
+
+
+
+
+/* ------------------------------------------------------------ */
+int main(int argc, char *argv[]) {
+
+#ifdef __GRASS__
+ /* GRASS initialization stuff */
+ struct GModule *module;
+
+ /*initialize GIS environment */
+ G_gisinit(argv[0]);
+
+ /*initialize module */
+ module = G_define_module();
+ module->keywords = _("raster, viewshed, line of sight");
+ module->description = _("IO-efficient viewshed algorithm");
+
+ struct Cell_head region;
+ if (G_get_set_window(®ion) == -1) {
+ G_fatal_error("error getting current region");
+ exit(EXIT_FAILURE);
+ }
+#endif
+
+
+
+ /* ************************************************************ */
+ /* parameters set up */
+ long long memSizeBytes = DEFAULT_MEMORY;
+ /* the maximum size of main memory that the program ca use. The
+ user can specify it, otherwise the default value of 500MB is
+ used. The program uses this value to decied in which mode to
+ run --- in internal memory, or external memory. */
+
+ int vpRow, vpCol;
+ /* the coordinates of the viewpoint in the raster; right now the
+ algorithm assumes that the viewpoint is inside the grid, though
+ this is not necessary; some changes will be needed to make it
+ work with a viewpoint outside the terrain */
+
+ ViewOptions viewOptions;
+ //viewOptions.inputfname = (char*)malloc(500);
+ //viewOptions.outputfname = (char*)malloc(500);
+ //assert(inputfname && outputfname);
+ viewOptions.obsElev = DEFAULT_OBS_ELEVATION;
+ viewOptions.maxDist = INFINITY_DISTANCE;
+ 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 */
+ Viewpoint vp;
+ 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);
+ assert(hd);
+
+ /* 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.*/
+ if (!(vp.row < hd->nrows && vp.col < hd->ncols)) {
+ G_fatal_error(_("Viewpoint outside grid"));
+ G_fatal_error(_("viewpont: (row=%d, col=%d)"), vp.row, vp.col);
+ G_fatal_error(_("grid: (rows=%d, cols=%d)"), hd->nrows, hd->ncols);
+ exit(EXIT_FAILURE);
+ }
+#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;
+ G_get_ellipsoid_parameters(&viewOptions.ellps_a, &e2);
+ if (viewOptions.ellps_a == 0) {
+ /*according to r.los, this can be
+ problematic, so we'll have a backup, hardcoded radius :-( */
+ G_warning(_
+ ("Problems obtaining current ellipsoid parameters, usting sphere (6370997.0)"));
+ viewOptions.ellps_a = 6370997.00;
+ }
+#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));
+ if (inmemSizeBytes < memSizeBytes) {
+ IN_MEMORY = 1;
+ print_message("*************\nIN_MEMORY MODE\n*************\n");
+ }
+ else {
+ print_message("*************\nEXTERNAL_MEMORY MODE\n**********\n");
+ 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");
+#endif
+
+#ifdef FORCE_INTERNAL
+ IN_MEMORY = 1;
+ print_message("FORCED INTERNAL\n");
+#endif
+
+
+ /* ************************************************************ */
+ /* compute viewshed in memory */
+ /* ************************************************************ */
+ if (IN_MEMORY) {
+ /*//////////////////////////////////////////////////// */
+ /*/viewshed in internal memory */
+ /*//////////////////////////////////////////////////// */
+ Rtimer totalTime, outputTime, sweepTime;
+ MemoryVisibilityGrid *visgrid;
+
+ rt_start(totalTime);
+
+ /*compute the viewshed and store it in visgrid */
+ rt_start(sweepTime);
+ visgrid = viewshed_in_memory(viewOptions.inputfname,hd,&vp, viewOptions);
+ rt_stop(sweepTime);
+
+ /* write the output */
+ rt_start(outputTime);
+ save_inmem_visibilitygrid(visgrid, viewOptions, vp);
+ rt_stop(outputTime);
+
+ rt_stop(totalTime);
+
+ print_timings_internal(sweepTime, outputTime, totalTime);
+ }
+
+
+
+
+ /* ************************************************************ */
+ /* compute viewshed in external memory */
+ /* ************************************************************ */
+ else {
+
+ /* ************************************************************ */
+ /* set up external memory mode */
+ /* setup STREAM_DIR if not already set */
+ char buf[1000];
+ if (getenv(STREAM_TMPDIR) != NULL) {
+ /*if already set */
+ fprintf(stderr, "%s=%s\n", STREAM_TMPDIR, getenv(STREAM_TMPDIR));
+ printf("Intermediate stream location: %s\n", getenv(STREAM_TMPDIR));
+ }
+ else {
+ /*set it */
+ sprintf(buf, "%s=%s", STREAM_TMPDIR, "/var/tmp/");
+ fprintf(stderr, "setting %s ", buf);
+ putenv(buf);
+ if (getenv(STREAM_TMPDIR) == NULL) {
+ fprintf(stderr, ", not set\n");
+ exit(1);
+ }
+ else {
+ fprintf(stderr, ", ok.\n");
+ }
+ printf("Intermediate stream location: %s\n", "/var/tmp/" );
+ }
+ fprintf(stderr, "Intermediate files will not be deleted "
+ "in case of abnormal termination.\n");
+ fprintf(stderr, "To save space delete these files manually!\n");
+
+
+ /* initialize IOSTREAM memory manager */
+ MM_manager.set_memory_limit(memSizeBytes);
+ MM_manager.warn_memory_limit();
+ MM_manager.print_limit_mode();
+
+
+
+ /* ************************************************************ */
+ /* BASE CASE OR DISTRIBUTION */
+ /* determine whether base-case of external algorithm is enough,
+ or recursion is necessary */
+ int BASE_CASE = 0;
+
+ if (get_active_str_size_bytes(hd) < memSizeBytes)
+ BASE_CASE = 1;
+
+ /*if the user set the FORCE_DISTRIBUTION flag, then the
+ algorithm runs in the fuly recursive mode (even if this is
+ not necessary). This is used solely for debugging purpses */
+#ifdef FORCE_DISTRIBUTION
+ BASE_CASE = 0;
+#endif
+
+
+
+
+ /* ************************************************************ */
+ /* external memory, base case */
+ /* ************************************************************ */
+ if (BASE_CASE) {
+ print_message
+ ("---Active structure small, starting base case---\n");
+
+ Rtimer totalTime, viewshedTime, outputTime, sortOutputTime;
+
+ rt_start(totalTime);
+
+ /*run viewshed's algorithm */
+ IOVisibilityGrid *visgrid;
+
+ rt_start(viewshedTime);
+ visgrid = viewshed_external(viewOptions.inputfname,hd,&vp,viewOptions);
+ rt_stop(viewshedTime);
+
+ /*sort output */
+ rt_start(sortOutputTime);
+ sort_io_visibilitygrid(visgrid);
+ rt_stop(sortOutputTime);
+
+ /*save output stream to file. */
+ rt_start(outputTime);
+ save_io_visibilitygrid(visgrid, viewOptions, vp);
+ rt_stop(outputTime);
+
+ rt_stop(totalTime);
+
+ print_timings_external_memory(totalTime, viewshedTime,
+ outputTime, sortOutputTime);
+ }
+
+
+
+
+
+ /************************************************************/
+ /* external memory, recursive distribution sweeping recursion */
+ /************************************************************ */
+ else { /* if not BASE_CASE */
+#ifndef FORCE_DISTRIBUTION
+ print_message("---Active structure does not fit in memory,");
+#else
+ print_message("FORCED DISTRIBUTION\n");
+#endif
+
+ Rtimer totalTime, sweepTime, outputTime, sortOutputTime;
+
+ rt_start(totalTime);
+
+ /*get the viewshed solution by distribution */
+ IOVisibilityGrid *visgrid;
+
+ rt_start(sweepTime);
+ visgrid = distribute_and_sweep(viewOptions.inputfname,hd,&vp,viewOptions);
+
+ rt_stop(sweepTime);
+
+ /*sort the visibility grid so that it is in order when it is
+ outputted */
+ rt_start(sortOutputTime);
+ sort_io_visibilitygrid(visgrid);
+ rt_stop(sortOutputTime);
+
+ rt_start(outputTime);
+ save_io_visibilitygrid(visgrid, viewOptions, vp);
+ rt_stop(outputTime);
+
+
+ rt_stop(totalTime);
+
+ print_timings_external_memory(totalTime, sweepTime,
+ outputTime, sortOutputTime);
+
+ }
+ }
+ /*end external memory, distribution sweep */
+
+
+
+
+
+
+ /* ************************************************************ */
+ /* 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;
+ G_short_history(viewOptions.outputfname, "raster", &history);
+ 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__
+void
+parse_args(int argc, char *argv[], int *vpRow, int *vpCol,
+ ViewOptions* viewOptions, long long *memSizeBytes,
+ Cell_head *window) {
+
+ assert(vpRow && vpCol && memSizeBytes && window);
+
+ /* the input */
+ struct Option *inputOpt;
+ inputOpt = G_define_standard_option(G_OPT_R_ELEV);
+ inputOpt->key = "input";
+
+ /* the output */
+ struct Option *outputOpt;
+ outputOpt = G_define_standard_option(G_OPT_R_OUTPUT);
+ outputOpt->description = "Name of output viewshed raser map\n\t\t\tdefault format: {NODATA, -1 (invisible), vertical angle wrt viewpoint (visible)}";
+
+
+ /* row-column flag */
+ struct Flag *row_col;
+ row_col = G_define_flag();
+ row_col->key = 'r';
+ row_col->description =
+ _("Use row-column location rather than latitude-longitude location");
+
+ /* curvature flag */
+ struct Flag *curvature;
+ curvature = G_define_flag();
+ curvature->key = 'c';
+ curvature->description =
+ _("Consider the curvature of the earth (current ellipsoid)");
+
+
+ /* boolean output flag */
+ struct Flag *booleanOutput;
+ booleanOutput = G_define_flag();
+ booleanOutput->key = 'b';
+ booleanOutput->description =
+ _("Output format is {0 (invisible) 1 (visible)}");
+
+ /* output mode = elevation flag */
+ struct Flag *elevationFlag;
+ elevationFlag = G_define_flag();
+ elevationFlag->key = 'e';
+ elevationFlag->description =
+ ("Output format is {NODATA, -1 (invisible), elev-viewpoint_elev (visible)}");
+
+ /* viewpoint coordinates */
+ struct Option *viewLocOpt;
+ viewLocOpt = G_define_option();
+ viewLocOpt->key = "viewpoint_location";
+ viewLocOpt->type = TYPE_STRING;
+ viewLocOpt->required = YES;
+ viewLocOpt->key_desc = "lat,long";
+ viewLocOpt->description =
+ ("Coordinates of viewing position in latitude-longitude (if -r flag is present, then coordinates are row-column)");
+
+ /* observer elevation */
+ struct Option *obsElevOpt;
+ obsElevOpt = G_define_option();
+ obsElevOpt->key = "observer_elevation";
+ obsElevOpt->type = TYPE_DOUBLE;
+ obsElevOpt->required = NO;
+ obsElevOpt->key_desc = "value";
+ obsElevOpt->description = _("Viewing elevation above the ground");
+ obsElevOpt->answer = "0.0";
+
+ /* max distance */
+ struct Option *maxDistOpt;
+ maxDistOpt = G_define_option();
+ maxDistOpt->key = "max_dist";
+ maxDistOpt->type = TYPE_DOUBLE;
+ maxDistOpt->required = NO;
+ maxDistOpt->key_desc = "value";
+ maxDistOpt->description =
+ ("Maximum visibility radius. By default infinity (-1).");
+ char infdist[10];
+ sprintf(infdist, "%d", INFINITY_DISTANCE);
+ maxDistOpt->answer = infdist;
+
+
+ /* memory size */
+ struct Option *memAmountOpt;
+ memAmountOpt = G_define_option();
+ memAmountOpt->key = "memory_usage";
+ memAmountOpt->type = TYPE_INTEGER;
+ memAmountOpt->required = NO;
+ memAmountOpt->key_desc = "value";
+ memAmountOpt->description = _("The amount of main memory in MB to be used");
+ memAmountOpt->answer = "500";
+
+
+ /*fill the options and flags with G_parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+
+ /* store the parameters into a structure to be used along the way */
+ strcpy(viewOptions->inputfname, inputOpt->answer);
+ strcpy(viewOptions->outputfname, outputOpt->answer);
+
+ viewOptions->obsElev = atof(obsElevOpt->answer);
+
+ viewOptions->maxDist = atof(maxDistOpt->answer);
+ if (viewOptions->maxDist < 0 &&
+ viewOptions->maxDist!= INFINITY_DISTANCE) {
+ G_fatal_error(_("negative max distance value is not valid"));
+ exit(EXIT_FAILURE);
+ }
+
+
+
+ viewOptions->doCurv = curvature->answer;
+ if (booleanOutput->answer)
+ viewOptions->outputMode = OUTPUT_BOOL;
+ else if (elevationFlag->answer)
+ viewOptions->outputMode = OUTPUT_ELEV;
+ else viewOptions->outputMode = OUTPUT_ANGLE;
+
+ int memSizeMB = atoi(memAmountOpt->answer);
+ if (memSizeMB <0) {
+ printf("Memory cannot be negative.\n");
+ exit(1);
+ }
+ *memSizeBytes = (long long)memSizeMB;
+ *memSizeBytes = (*memSizeBytes) << 20;
+
+ /*The algorithm runs with the viewpoint row and col, so depending
+ on if the row_col flag is present we either need to store the
+ row and col, or convert the lat-lon coordinates to row and
+ column format */
+ 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);
+ }
+ 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);
+
+ }
+
+
+ return;
+}
+
+
+
+
+/* ------------------------------------------------------------ */
+#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
+print_timings_internal(Rtimer sweepTime, Rtimer outputTime, Rtimer totalTime) {
+
+ char timeused[100];
+ printf("TOTAL TIMING: \n");
+ rt_sprint_safe(timeused, sweepTime);
+ printf("\t%30s", "sweep:");
+ printf(timeused); printf("\n");
+
+ rt_sprint_safe(timeused, outputTime);
+ printf("\t%30s", "output:");
+ printf(timeused); printf("\n");
+
+ rt_sprint_safe(timeused, totalTime);
+ printf("\t%30s", "total:");
+ printf(timeused);
+ printf("\n");
+}
+
+
+/* ------------------------------------------------------------ */
+/*print the timings for the external memory method of solving the viewshed */
+void
+print_timings_external_memory(Rtimer totalTime, Rtimer viewshedTime,
+ Rtimer outputTime, Rtimer sortOutputTime)
+{
+
+ /*print timings */
+ char timeused[100];
+
+ printf("\n\nTOTAL TIMING: \n");
+ rt_sprint_safe(timeused, viewshedTime);
+ printf("\t%30s", "total sweep:");
+ printf(timeused);
+ printf("\n");
+ rt_sprint_safe(timeused, sortOutputTime);
+ printf("\t%30s", "sort output:");
+ printf(timeused);
+ printf("\n");
+ rt_sprint_safe(timeused, outputTime);
+ printf("\t%30s", "Write result grid:");
+ printf(timeused);
+ printf("\n");
+ rt_sprint_safe(timeused, totalTime);
+ printf("\t%30s", "Total Time:");
+ printf(timeused);
+ printf("\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"));
+
+#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");
+
+#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");
+}
Property changes on: grass/trunk/raster/r.viewshed/main.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/print_message.cc
===================================================================
--- grass/trunk/raster/r.viewshed/print_message.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/print_message.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,64 @@
+
+/****************************************************************************
+ *
+ * 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
+}
Property changes on: grass/trunk/raster/r.viewshed/print_message.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/print_message.h
===================================================================
--- grass/trunk/raster/r.viewshed/print_message.h (rev 0)
+++ grass/trunk/raster/r.viewshed/print_message.h 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,46 @@
+
+/****************************************************************************
+ *
+ * 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
Property changes on: grass/trunk/raster/r.viewshed/print_message.h
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/r.viewshed.html
===================================================================
--- grass/trunk/raster/r.viewshed/r.viewshed.html (rev 0)
+++ grass/trunk/raster/r.viewshed/r.viewshed.html 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,225 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS: r.viewshed</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+<h2>NAME</h2>
+<em><b>r.viewshed</b></em>
+
+<H2>DESCRIPTION</H2>
+
+<p><em>r.viewshed</em> is a module that computes the viewshed of a
+point on a raster terrain. That is, given an elevation raster, and the
+location of an observer, it generates a raster output map showing
+which cells are visible from the given location.
+
+The algorithm underlying <em>r.viewshed</em> minimizes both the CPU
+operations and the transfer of data between main memory and disk; as a
+result <em>r.viewshed</em> runs fast on very large rasters.
+
+<h3>Options and flags:</h3>
+
+To run <em>r.viewshed</em>, the user must specify an input map name,
+an output map name, and the location of the viewpoint.
+
+<p>Viewpoint: For the time being the viewpoint is assumed to be
+located inside the terrain. The viewpoint location is, by default, in
+coordinate format, but if the <b>-r</b> flag is set, it is in
+row-column format instead.
+
+<p>Output: The output map may have one of three possible formats,
+based on which flags are set.
+
+<ol>
+<li>By default, if no flag is set, the output is in angle-mode, and
+each point in the output map is marked as:
+
+<ul>
+<li> Nodata/null, if the respective point in the elevation map is
+Nodata/null
+
+<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.
+
+</ul>
+
+
+<li>If the -b flag is set, the output is in boolean-mode, and each point
+in the output map is marked as:
+
+<ul>
+
+<li> 0 if the point is Nodata/null or not visible
+
+<li> 1 if the point is visible.
+
+</ul>
+
+
+<li>If the -e flag is set, the output is in elevation-mode, and each point
+in the output map is marked as:
+
+<ul>
+<li> Nodata (null), if the respective point in the elevation map is
+Nodata (null)
+
+<li> -1, if the point is not visible
+
+<li> the difference in elevation between the point and the viewpoint,
+if the point is visible.
+
+</ul>
+
+</ol>
+
+<p>Curvature of the earth: By default the elevations are not adjusted for
+the curvature of the earth. The user can turn this on with flag
+<b>-c</b>.
+
+<p>Observer elevation: By default the observer is assumed to have
+height=0 above the terrain. The user can change this using option
+<em>observer_elevation=...</em>. 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>.
+The value entered is in the same units as the cell size of the raster.
+
+
+<p>Main memory usage: By default <em>r.viewshed</em> assumes it has
+500MB of main memory, and sets up its internal data structures so that
+it does not require more than this amount of RAM. The user can set
+the amount of memory used by the program by setting the
+<em>memory_usage</em> to the number of MB of memory they would like to
+be used.
+
+<h3>Memory mode</h3>
+<p> The algorithm can run in two modes: in internal memory, which
+means that it keeps all necessary data structures in memory during the
+computation. And in external memory, which means that the data
+structures are external, i.e. on disk. <em>r.viewshed</em> decides
+which mode to run in using the amount of main memory specified by the
+user. The internal mode is (much) faster than the external mode.
+
+<p>Ideally, the user should specify on the command line the amount of
+physical memory that is free for the program to use. Underestimating
+the memory may result in <em>r.viewshed</em> running in external mode
+instead of internal, which is slower. Overestimating the amount of
+free memory may result in <em>r.viewshed</em> running in internal mode
+and using virtual memory, which is slower than the external mode.
+
+
+
+
+<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.
+
+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.
+
+<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
+a (square) grid of <em>n</em> cells, there can be <em>O(n
+<sup>1/2</sup>)</em> cells that intersect the LOS. If we test every
+single such cell for every point in the grid, this adds up to
+<em>O(n<sup>3/2</sup>)</em> tests. We can do all these tests faster if
+we re-use information from one point to the next (two grid points that
+are close to each other will be intersected by a lot of the same
+points) and organize the computation differently.
+
+<p>More precisely, the algorithm uses a technique called <em>line
+sweeping</em>: It considers a half-line centered at the viewpoint, and
+rotates it radially around the viewpoint, 360 degrees. During the
+sweep it keeps track of all the cells that intersect the sweep line at
+that time; These are called the <em>active</em> cells. A cell has 3
+associated events: when it is first met by the sweep line and inserted
+into the active structure; when it is last met by the sweep line and
+deleted from the active structure; and when the sweep line passes over
+its centerpoint, at which time its visibility is determined. To
+determine the visibility of a cell all cells that intersect the
+line-of-sight must be active, so they are in the active structure.
+The algorithm looks at all the active cells that are between the point
+and the viewpoint, and finds the maximum gradient among these. If the
+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
+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
+made efficient in terms of I/O-operations. For all details see the
+REFERENCES below.
+
+
+<p>
+<table width=50% align=center>
+ <tr>
+ <th><img src="sweep1.jpg" width=200 alt="[SDF]" border=0></th>
+ <th><img src="sweep2.jpg" width=200 alt="[SDF]" border=0></th>
+ </tr>
+ <tr>
+ <th>The sweep-line.</th>
+ <th>The active cells.</th>
+ </tr>
+</table>
+
+
+<h3>An example:</h3>
+
+
+Using the Spearfish dataset: calculating the viewpoint from the top
+of a moutain:
+
+<div class="code"><pre> r.viewshed input=elevation output=viewshed viewpoint_location=598869,4916642 mem=800 </pre></div>
+
+<h3>REFERENCES</h3>
+
+<ul>
+
+ <li>Computing Visibility on Terrains in External Memory. Herman
+ Haverkort, Laura Toma and Yi Zhuang. To appear in <em>ACM Journal
+ on Experimental Algorithmics</em>.
+
+ <li><a
+ href="http://www.siam.org/proceedings/alenex/2007/alx07_002haverkorth.pdf">
+ Computing Visibility on Terrains in External Memory</a>. Herman
+ Haverkort, Laura Toma and Yi Zhuang. In the <em>Proceedings of
+ the 9th Workshop on Algorithm Engineering and Experiments /
+ Workshop on Analytic Algorithms and Combinatorics (ALENEX/ANALCO
+ 2007)</em>.</li>
+
+</ul>
+
+<h3>AUTHORS</h3>
+
+<p>Laura Toma (Bowdoin College): <tt>ltoma at bowdoin.edu</tt>
+
+<p> Yi Zhuang (Carnegie-Mellon University): <tt>yzhuang at andrew.cmu.edu</tt>
+
+<p>William Richard (Bowdoin College): <tt>willster3021 at gmail.com </tt>
+<HR>
+<P><a href="index.html">Main index</a> - <a href="raster.html">raster index</a> - <a href="full_index.html">Full index</a></P>
+<P>© 2003-2008 <a href="http://grass.osgeo.org">GRASS Development Team</a></p>
+</body>
+</html>
Added: grass/trunk/raster/r.viewshed/rbbst.cc
===================================================================
--- grass/trunk/raster/r.viewshed/rbbst.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/rbbst.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,764 @@
+
+/****************************************************************************
+ *
+ * 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.
+ *
+ *****************************************************************************/
+
+/*
+
+ A R/B BST. Always call initNILnode() before using the tree.
+ Version 0.0.0
+
+ Version 0.0.1
+ Rewrote BST Deletion to improve efficiency
+
+ Version 0.0.2
+ Bug fixed in deletion.
+ CLRS pseudocode forgot to make sure that x is not NIL before
+ calling rbDeleteFixup(root,x).
+
+ Version 0.0.3
+ Some Cleanup. Separated the public portion and the
+ private porthion of the interface in the header
+
+
+ =================================
+ This is based on BST 1.0.4
+ BST change log
+ <---------------->
+ find max is implemented in this version.
+ Version 1.0.2
+
+ Version 1.0.4
+ Major bug fix in deletion (when the node has two children,
+ one of them has a wrong parent pointer after the rotation in the deletion.)
+ <----------------->
+ */
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include "rbbst.h"
+
+#ifdef __GRASS__
+extern "C"
+{
+#include <grass/config.h>
+#include <grass/gis.h>
+}
+#endif
+
+
+
+TreeNode *NIL;
+
+#define EPSILON 0.0000001
+
+
+/*public:--------------------------------- */
+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;
+ rbt->root->left = NIL;
+ rbt->root->right = NIL;
+ rbt->root->parent = NIL;
+ rbt->root->color = RB_BLACK;
+
+ return rbt;
+}
+
+/*LT: not sure if this is correct */
+int is_empty(RBTree * t)
+{
+ assert(t);
+ return (t->root == NIL);
+}
+
+void delete_tree(RBTree * t)
+{
+ destroy_sub_tree(t->root);
+ return;
+}
+
+void destroy_sub_tree(TreeNode * node)
+{
+ if (node == NIL)
+ return;
+
+ destroy_sub_tree(node->left);
+ destroy_sub_tree(node->right);
+#ifdef __GRASS__
+ G_free(node);
+#else
+ free(node);
+#endif
+ return;
+}
+
+void insert_into(RBTree * rbt, TreeValue value)
+{
+ insert_into_tree(&(rbt->root), value);
+ return;
+}
+
+void delete_from(RBTree * rbt, double key)
+{
+ delete_from_tree(&(rbt->root), key);
+ return;
+}
+
+TreeNode *search_for_node_with_key(RBTree * rbt, double key)
+{
+ return search_for_node(rbt->root, key);
+}
+
+/*------------The following is designed for viewshed's algorithm-------*/
+double find_max_gradient_within_key(RBTree * rbt, double key)
+{
+ return find_max_value_within_key(rbt->root, key);
+}
+
+/*<--------------------------------->
+ //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.maxGradient = SMALLEST_GRADIENT;
+
+ NIL->parent = NULL;
+ NIL->left = NULL;
+ NIL->right = NULL;
+ return;
+}
+
+/*you can write change this compare function, depending on your TreeValue struct
+ //compare function used by findMaxValue
+ //-1: v1 < v2
+ //0: v1 = v2
+ //2: v1 > v2 */
+char compare_values(TreeValue * v1, TreeValue * v2)
+{
+ if (v1->gradient > v2->gradient)
+ return 1;
+ if (v1->gradient < v2->gradient)
+ return -1;
+
+ return 0;
+}
+
+
+/*a function used to compare two doubles */
+char compare_double(double a, double b)
+{
+ if (fabs(a - b) < EPSILON)
+ return 0;
+ if (a - b < 0)
+ return -1;
+
+ return 1;
+}
+
+
+
+/*create a tree node */
+TreeNode *create_tree_node(TreeValue value)
+{
+ TreeNode *ret;
+
+#ifdef __GRASS__
+ ret = (TreeNode *) G_malloc(sizeof(TreeNode));
+#else
+ ret = (TreeNode *) malloc(sizeof(TreeNode));
+#endif
+
+ ret->color = RB_RED;
+
+ ret->left = NIL;
+ ret->right = NIL;
+ ret->parent = NIL;
+
+ ret->value = value;
+ ret->value.maxGradient = SMALLEST_GRADIENT;
+ return ret;
+}
+
+/*create node with its value set to the value given
+ //and insert the node into the tree
+ //rbInsertFixup may change the root pointer, so TreeNode** is passed in */
+void insert_into_tree(TreeNode ** root, TreeValue value)
+{
+ TreeNode *curNode;
+ TreeNode *nextNode;
+
+ curNode = *root;
+
+ if (compare_double(value.key, curNode->value.key) == -1) {
+ nextNode = curNode->left;
+ }
+ else {
+ nextNode = curNode->right;
+ }
+
+
+ while (nextNode != NIL) {
+ curNode = nextNode;
+
+ if (compare_double(value.key, curNode->value.key) == -1) {
+ nextNode = curNode->left;
+ }
+ else {
+ nextNode = curNode->right;
+ }
+ }
+
+ /*create a new node
+ //and place it at the right place
+ //created node is RED by default */
+ nextNode = create_tree_node(value);
+
+ nextNode->parent = curNode;
+
+ if (compare_double(value.key, curNode->value.key) == -1) {
+ curNode->left = nextNode;
+ }
+ else {
+ curNode->right = nextNode;
+ }
+
+ TreeNode *inserted = nextNode;
+
+ /*update augmented maxGradient */
+ nextNode->value.maxGradient = nextNode->value.gradient;
+ while (nextNode->parent != NIL) {
+ if (nextNode->parent->value.maxGradient < nextNode->value.maxGradient)
+ nextNode->parent->value.maxGradient = nextNode->value.maxGradient;
+
+ if (nextNode->parent->value.maxGradient > nextNode->value.maxGradient)
+ break;
+ nextNode = nextNode->parent;
+ }
+
+ /*fix rb tree after insertion */
+ rb_insert_fixup(root, inserted);
+
+ return;
+}
+
+void rb_insert_fixup(TreeNode ** root, TreeNode * z)
+{
+ /*see pseudocode on page 281 in CLRS */
+ TreeNode *y;
+
+ while (z->parent->color == RB_RED) {
+ if (z->parent == z->parent->parent->left) {
+ y = z->parent->parent->right;
+ if (y->color == RB_RED) { /*case 1 */
+ z->parent->color = RB_BLACK;
+ y->color = RB_BLACK;
+ z->parent->parent->color = RB_RED;
+ z = z->parent->parent;
+ }
+ else {
+ if (z == z->parent->right) { /*case 2 */
+ z = z->parent;
+ left_rotate(root, z); /*convert case 2 to case 3 */
+ }
+ z->parent->color = RB_BLACK; /*case 3 */
+ z->parent->parent->color = RB_RED;
+ right_rotate(root, z->parent->parent);
+ }
+
+ }
+ else { /*(z->parent == z->parent->parent->right) */
+ y = z->parent->parent->left;
+ if (y->color == RB_RED) { /*case 1 */
+ z->parent->color = RB_BLACK;
+ y->color = RB_BLACK;
+ z->parent->parent->color = RB_RED;
+ z = z->parent->parent;
+ }
+ else {
+ if (z == z->parent->left) { /*case 2 */
+ z = z->parent;
+ right_rotate(root, z); /*convert case 2 to case 3 */
+ }
+ z->parent->color = RB_BLACK; /*case 3 */
+ z->parent->parent->color = RB_RED;
+ left_rotate(root, z->parent->parent);
+ }
+ }
+ }
+ (*root)->color = RB_BLACK;
+
+ return;
+}
+
+
+
+
+/*search for a node with the given key */
+TreeNode *search_for_node(TreeNode * root, double key)
+{
+ TreeNode *curNode = root;
+
+ while (curNode != NIL && compare_double(key, curNode->value.key) != 0) {
+
+ if (compare_double(key, curNode->value.key) == -1) {
+ curNode = curNode->left;
+ }
+ else {
+ curNode = curNode->right;
+ }
+
+ }
+
+ return curNode;
+}
+
+/*function used by treeSuccessor */
+TreeNode *tree_minimum(TreeNode * x)
+{
+ while (x->left != NIL)
+ x = x->left;
+
+ return x;
+}
+
+/*function used by deletion */
+TreeNode *tree_successor(TreeNode * x)
+{
+ if (x->right != NIL)
+ return tree_minimum(x->right);
+ TreeNode *y = x->parent;
+
+ while (y != NIL && x == y->right) {
+ x = y;
+ y = y->parent;
+ }
+ return y;
+}
+
+
+/*delete the node out of the tree */
+void delete_from_tree(TreeNode ** root, double key)
+{
+ double tmpMax;
+ TreeNode *z;
+ TreeNode *x;
+ TreeNode *y;
+ TreeNode *toFix;
+
+ 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 */
+ }
+
+ /*1-3 */
+ if (z->left == NIL || z->right == NIL)
+ y = z;
+ else
+ y = tree_successor(z);
+
+ /*4-6 */
+ if (y->left != NIL)
+ x = y->left;
+ else
+ x = y->right;
+
+ /*7 */
+ x->parent = y->parent;
+
+ /*8-12 */
+ if (y->parent == NIL) {
+ *root = x;
+
+ toFix = *root; /*augmentation to be fixed */
+ }
+ else {
+ if (y == y->parent->left)
+ y->parent->left = x;
+ else
+ y->parent->right = x;
+
+ toFix = y->parent; /*augmentation to be fixed */
+ }
+
+ /*fix augmentation for removing y */
+ TreeNode *curNode = y;
+ double left, right;
+
+ while (curNode->parent != NIL) {
+ if (curNode->parent->value.maxGradient == y->value.gradient) {
+ left = find_max_value(curNode->parent->left);
+ right = find_max_value(curNode->parent->right);
+
+ if (left > right)
+ curNode->parent->value.maxGradient = left;
+ else
+ curNode->parent->value.maxGradient = right;
+
+ if (curNode->parent->value.gradient >
+ curNode->parent->value.maxGradient)
+ curNode->parent->value.maxGradient =
+ curNode->parent->value.gradient;
+ }
+ else {
+ break;
+ }
+ curNode = curNode->parent;
+ }
+
+
+ /*fix augmentation for x */
+ tmpMax =
+ toFix->left->value.maxGradient >
+ toFix->right->value.maxGradient ? toFix->left->value.
+ maxGradient : toFix->right->value.maxGradient;
+ if (tmpMax > toFix->value.gradient)
+ toFix->value.maxGradient = tmpMax;
+ else
+ toFix->value.maxGradient = toFix->value.gradient;
+
+ /*13-15 */
+ if (y != z) {
+ double zGradient = z->value.gradient;
+
+ z->value.key = y->value.key;
+ z->value.gradient = y->value.gradient;
+
+
+ toFix = z;
+ /*fix augmentation */
+ tmpMax =
+ toFix->left->value.maxGradient >
+ toFix->right->value.maxGradient ? toFix->left->value.
+ maxGradient : toFix->right->value.maxGradient;
+ if (tmpMax > toFix->value.gradient)
+ toFix->value.maxGradient = tmpMax;
+ else
+ toFix->value.maxGradient = toFix->value.gradient;
+
+ while (z->parent != NIL) {
+ if (z->parent->value.maxGradient == zGradient) {
+ if (z->parent->value.gradient != zGradient &&
+ (!(z->parent->left->value.maxGradient == zGradient &&
+ z->parent->right->value.maxGradient == zGradient))) {
+
+ left = find_max_value(z->parent->left);
+ right = find_max_value(z->parent->right);
+
+ if (left > right)
+ z->parent->value.maxGradient = left;
+ else
+ z->parent->value.maxGradient = right;
+
+ if (z->parent->value.gradient >
+ z->parent->value.maxGradient)
+ z->parent->value.maxGradient =
+ z->parent->value.gradient;
+
+ }
+
+ }
+ else {
+ if (z->value.maxGradient > z->parent->value.maxGradient)
+ z->parent->value.maxGradient = z->value.maxGradient;
+ }
+ z = z->parent;
+ }
+
+ }
+
+ /*16-17 */
+ if (y->color == RB_BLACK && x != NIL)
+ rb_delete_fixup(root, x);
+
+ /*18 */
+#ifdef __GRASS__
+ G_free(y);
+#else
+ free(y);
+#endif
+
+ return;
+}
+
+/*fix the rb tree after deletion */
+void rb_delete_fixup(TreeNode ** root, TreeNode * x)
+{
+ TreeNode *w;
+
+ while (x != *root && x->color == RB_BLACK) {
+ if (x == x->parent->left) {
+ w = x->parent->right;
+ if (w->color == RB_RED) {
+ w->color = RB_BLACK;
+ x->parent->color = RB_RED;
+ left_rotate(root, x->parent);
+ w = x->parent->right;
+ }
+
+ if (w == NIL) {
+ x = x->parent;
+ continue;
+ }
+
+ if (w->left->color == RB_BLACK && w->right->color == RB_BLACK) {
+ w->color = RB_RED;
+ x = x->parent;
+ }
+ else {
+ if (w->right->color == RB_BLACK) {
+ w->left->color = RB_BLACK;
+ w->color = RB_RED;
+ right_rotate(root, w);
+ w = x->parent->right;
+ }
+
+ w->color = x->parent->color;
+ x->parent->color = RB_BLACK;
+ w->right->color = RB_BLACK;
+ left_rotate(root, x->parent);
+ x = *root;
+ }
+
+ }
+ else { /*(x==x->parent->right) */
+ w = x->parent->left;
+ if (w->color == RB_RED) {
+ w->color = RB_BLACK;
+ x->parent->color = RB_RED;
+ right_rotate(root, x->parent);
+ w = x->parent->left;
+ }
+
+ if (w == NIL) {
+ x = x->parent;
+ continue;
+ }
+
+ if (w->right->color == RB_BLACK && w->left->color == RB_BLACK) {
+ w->color = RB_RED;
+ x = x->parent;
+ }
+ else {
+ if (w->left->color == RB_BLACK) {
+ w->right->color = RB_BLACK;
+ w->color = RB_RED;
+ left_rotate(root, w);
+ w = x->parent->left;
+ }
+
+ w->color = x->parent->color;
+ x->parent->color = RB_BLACK;
+ w->left->color = RB_BLACK;
+ right_rotate(root, x->parent);
+ x = *root;
+ }
+
+ }
+ }
+ x->color = RB_BLACK;
+
+ return;
+}
+
+/*find the max value in the given tree
+ //you need to provide a compare function to compare the nodes */
+double find_max_value(TreeNode * root)
+{
+ if (!root)
+ return SMALLEST_GRADIENT;
+ assert(root);
+ /*assert(root->value.maxGradient != SMALLEST_GRADIENT);
+ //LT: this shoudl be fixed
+ //if (root->value.maxGradient != SMALLEST_GRADIENT) */
+ return root->value.maxGradient;
+}
+
+
+
+/*find max within the max key */
+double find_max_value_within_key(TreeNode * root, double maxKey)
+{
+ TreeNode *keyNode = search_for_node(root, maxKey);
+
+ if (keyNode == NIL) {
+ /*fprintf(stderr, "key node not found. error occured!\n");
+ //there is no point in the structure with key < maxKey */
+ return SMALLEST_GRADIENT;
+ exit(1);
+ }
+
+ double max = find_max_value(keyNode->left);
+ double tmpMax;
+
+ 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;
+ }
+ keyNode = keyNode->parent;
+ }
+
+ return max;
+}
+
+
+void left_rotate(TreeNode ** root, TreeNode * x)
+{
+ TreeNode *y;
+
+ y = x->right;
+
+ /*maintain augmentation */
+ double tmpMax;
+
+ /*fix x */
+ tmpMax = x->left->value.maxGradient > y->left->value.maxGradient ?
+ x->left->value.maxGradient : y->left->value.maxGradient;
+
+ if (tmpMax > x->value.gradient)
+ x->value.maxGradient = tmpMax;
+ else
+ x->value.maxGradient = x->value.gradient;
+
+
+ /*fix y */
+ tmpMax = x->value.maxGradient > y->right->value.maxGradient ?
+ x->value.maxGradient : y->right->value.maxGradient;
+
+ if (tmpMax > y->value.gradient)
+ y->value.maxGradient = tmpMax;
+ else
+ y->value.maxGradient = y->value.gradient;
+
+ /*left rotation
+ //see pseudocode on page 278 in CLRS */
+
+ x->right = y->left; /*turn y's left subtree into x's right subtree */
+ y->left->parent = x;
+
+ y->parent = x->parent; /*link x's parent to y */
+
+ if (x->parent == NIL) {
+ *root = y;
+ }
+ else {
+ if (x == x->parent->left)
+ x->parent->left = y;
+ else
+ x->parent->right = y;
+ }
+
+ y->left = x;
+ x->parent = y;
+
+ return;
+}
+
+void right_rotate(TreeNode ** root, TreeNode * y)
+{
+ TreeNode *x;
+
+ x = y->left;
+
+ /*maintain augmentation
+ //fix y */
+ double tmpMax;
+
+ tmpMax = x->right->value.maxGradient > y->right->value.maxGradient ?
+ x->right->value.maxGradient : y->right->value.maxGradient;
+
+ if (tmpMax > y->value.gradient)
+ y->value.maxGradient = tmpMax;
+ else
+ y->value.maxGradient = y->value.gradient;
+
+ /*fix x */
+ tmpMax = x->left->value.maxGradient > y->value.maxGradient ?
+ x->left->value.maxGradient : y->value.maxGradient;
+
+ if (tmpMax > x->value.gradient)
+ x->value.maxGradient = tmpMax;
+ else
+ x->value.maxGradient = x->value.gradient;
+
+ /*ratation */
+ y->left = x->right;
+ x->right->parent = y;
+
+ x->parent = y->parent;
+
+ if (y->parent == NIL) {
+ *root = x;
+ }
+ else {
+ if (y->parent->left == y)
+ y->parent->left = x;
+ else
+ y->parent->right = x;
+ }
+
+ x->right = y;
+ y->parent = x;
+
+ return;
+}
Property changes on: grass/trunk/raster/r.viewshed/rbbst.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/rbbst.h
===================================================================
--- grass/trunk/raster/r.viewshed/rbbst.h (rev 0)
+++ grass/trunk/raster/r.viewshed/rbbst.h 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,156 @@
+/****************************************************************************
+ *
+ * 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 __RB_BINARY_SEARCH_TREE__
+#define __RB_BINARY_SEARCH_TREE__
+
+#define SMALLEST_GRADIENT (- 9999999999999999999999.0)
+/*this value is returned by findMaxValueWithinDist() is there is no
+ key within that distance. The largest double value is 1.7 E 308*/
+
+
+
+#define RB_RED (0)
+#define RB_BLACK (1)
+
+/*<===========================================>
+ //public:
+ //The value that's stored in the tree
+ //Change this structure to avoid type casting at run time */
+typedef struct tree_value_
+{
+ /*this field is mandatory and cannot be removed.
+ //the tree is indexed by this "key". */
+ double key;
+
+ /*anything else below this line is optional */
+ double gradient;
+ double maxGradient;
+} TreeValue;
+
+
+/*The node of a tree */
+typedef struct tree_node_
+{
+ TreeValue value;
+
+ char color;
+
+ struct tree_node_ *left;
+ struct tree_node_ *right;
+ struct tree_node_ *parent;
+
+} TreeNode;
+
+typedef struct rbtree_
+{
+ TreeNode *root;
+} RBTree;
+
+
+
+RBTree *create_tree(TreeValue tv);
+void delete_tree(RBTree * t);
+void destroy_sub_tree(TreeNode * node);
+void insert_into(RBTree * rbt, TreeValue value);
+void delete_from(RBTree * rbt, double key);
+TreeNode *search_for_node_with_key(RBTree * rbt, double key);
+
+
+/*------------The following is designed for kreveld's algorithm-------*/
+double find_max_gradient_within_key(RBTree * rbt, double key);
+
+/*LT: not sure if this is correct */
+int is_empty(RBTree * t);
+
+
+
+
+
+/*<================================================>
+ //private:
+ //The below are private functions you should not
+ //call directly when using the Tree
+
+ //<--------------------------------->
+ //for RB tree only
+
+ //in RB TREE, used to replace NULL */
+void init_nil_node();
+
+
+/*Left and Right Rotation
+ //the root of the tree may be modified during the rotations
+ //so TreeNode** is passed into the functions */
+void left_rotate(TreeNode ** root, TreeNode * x);
+void right_rotate(TreeNode ** root, TreeNode * y);
+void rb_insert_fixup(TreeNode ** root, TreeNode * z);
+void rb_delete_fixup(TreeNode ** root, TreeNode * x);
+
+/*<------------------------------------> */
+
+
+/*compare function used by findMaxValue
+ //-1: v1 < v2
+ //0: v1 = v2
+ //2: v1 > v2 */
+char compare_values(TreeValue * v1, TreeValue * v2);
+
+/*a function used to compare two doubles */
+char compare_double(double a, double b);
+
+/*create a tree node */
+TreeNode *create_tree_node(TreeValue value);
+
+/*create node with its value set to the value given
+ //and insert the node into the tree */
+void insert_into_tree(TreeNode ** root, TreeValue value);
+
+/*delete the node out of the tree */
+void delete_from_tree(TreeNode ** root, double key);
+
+/*search for a node with the given key */
+TreeNode *search_for_node(TreeNode * root, double key);
+
+/*find the max value in the given tree
+ //you need to provide a compare function to compare the nodes */
+double find_max_value(TreeNode * root);
+
+/*find max within the max key */
+double find_max_value_within_key(TreeNode * root, double maxKey);
+
+#endif
Property changes on: grass/trunk/raster/r.viewshed/rbbst.h
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/statusstructure.cc
===================================================================
--- grass/trunk/raster/r.viewshed/statusstructure.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/statusstructure.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,214 @@
+/****************************************************************************
+ *
+ * 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 <math.h>
+#include <stdio.h>
+#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"
+
+
+/*SMALLEST_GRADIENT is defined in rbbst.h */
+
+
+/* ------------------------------------------------------------ */
+/*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
+ earth */
+float get_vertical_angle(Viewpoint vp, StatusNode sn, int doCurv) {
+
+ /*determine the difference in elevation, based on the curvature*/
+ float diffElev;
+ diffElev = vp.elev - sn.elev;
+
+ /*calculate and return the angle in degrees*/
+ assert(fabs(sn.dist2vp) > 0.001);
+ if(diffElev >= 0.0)
+ return (atan(sn.dist2vp / diffElev) * (180/PI));
+ else
+ return ((atan(fabs(diffElev) / sn.dist2vp) * (180/PI)) + 90);
+}
+
+
+
+/* ------------------------------------------------------------ */
+/*return an estimate of the size of active structure */
+long long get_active_str_size_bytes(GridHeader * hd)
+{
+
+ long long sizeBytes;
+
+ printf("Estimated size active structure:");
+ printf(" (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);
+ return sizeBytes;
+}
+
+
+
+
+/* ------------------------------------------------------------ */
+/*given a StatusNode, fill in its dist2vp and gradient */
+void calculate_dist_n_gradient(StatusNode * sn, Viewpoint * vp)
+{
+
+ 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); */
+ sn->dist2vp = (sn->row - vp->row) * (sn->row - vp->row) +
+ (sn->col - vp->col) * (sn->col - vp->col);
+ sn->gradient =
+ (sn->elev - vp->elev) * (sn->elev - vp->elev) / (sn->dist2vp);
+ /*maintain sign */
+ if (sn->elev < vp->elev)
+ sn->gradient = -sn->gradient;
+ 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.key = 0;
+ tv.maxGradient = SMALLEST_GRADIENT;
+
+
+ sl->rbt = create_tree(tv);
+ return sl;
+}
+
+
+/* ------------------------------------------------------------ */
+/*delete a status structure */
+void delete_status_structure(StatusList * sl)
+{
+ assert(sl);
+ delete_tree(sl->rbt);
+#ifdef __GRASS__
+ G_free(sl);
+#else
+ free(sl);
+#endif
+ return;
+}
+
+
+/* ------------------------------------------------------------ */
+/*delete the statusNode with the given key */
+void delete_from_status_struct(StatusList * sl, double dist2vp)
+{
+ assert(sl);
+ delete_from(sl->rbt, dist2vp);
+ return;
+}
+
+
+
+
+/* ------------------------------------------------------------ */
+/*insert the element into the status structure */
+void insert_into_status_struct(StatusNode sn, StatusList * sl)
+{
+ assert(sl);
+ TreeValue tv;
+
+ tv.key = sn.dist2vp;
+ tv.gradient = sn.gradient;
+ tv.maxGradient = SMALLEST_GRADIENT;
+ insert_into(sl->rbt, tv);
+ return;
+}
+
+
+/* ------------------------------------------------------------ */
+/*find the node with max Gradient within the distance (from viewpoint)
+ //given */
+double find_max_gradient_in_status_struct(StatusList * sl, double dist)
+{
+ assert(sl);
+ /*note: if there is nothing in the status struccture, it means this
+ cell is VISIBLE */
+ if (is_empty(sl))
+ return SMALLEST_GRADIENT;
+ /*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);
+}
+
+/*returns true is it is empty */
+int is_empty(StatusList * sl)
+{
+ assert(sl);
+ return (is_empty(sl->rbt) ||
+ sl->rbt->root->value.maxGradient == SMALLEST_GRADIENT);
+}
+
Property changes on: grass/trunk/raster/r.viewshed/statusstructure.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/statusstructure.h
===================================================================
--- grass/trunk/raster/r.viewshed/statusstructure.h (rev 0)
+++ grass/trunk/raster/r.viewshed/statusstructure.h 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,107 @@
+/****************************************************************************
+ *
+ * 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 _YZ__STATUSSTRUCTURE_H
+#define _YZ__STATUSSTRUCTURE_H
+
+/*
+ This header file defines the status structure and related functions.
+ */
+
+#ifdef __GRASS__
+#include <grass/config.h>
+#endif
+
+#include "grid.h"
+#include "rbbst.h"
+#include "visibility.h"
+
+#define PI 3.1416
+
+typedef struct statusnode_
+{
+ dimensionType row, col; /*position of the cell */
+ float elev; /*elevation of cell */
+ double dist2vp; /*distance to the viewpoint */
+ double gradient; /*gradient of the Line of Sight */
+} StatusNode;
+
+
+typedef struct statuslist_
+{
+ RBTree *rbt; /*pointer to the root of the bst */
+} StatusList;
+
+
+
+/* ------------------------------------------------------------ */
+
+/*return an estimate of the size of active structure */
+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);
+
+
+/*create an empty status list. */
+StatusList *create_status_struct();
+
+void delete_status_structure(StatusList * sl);
+
+/*returns true is it is empty */
+int is_empty(StatusList * sl);
+
+
+/*delete the statusNode with the given key */
+void delete_from_status_struct(StatusList * sl, double dist2vp);
+
+/*insert the element into the status structure */
+void insert_into_status_struct(StatusNode sn, StatusList * sl);
+
+/*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);
+
+/*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);
+
+
+#endif
Property changes on: grass/trunk/raster/r.viewshed/statusstructure.h
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/sweep1.png
===================================================================
(Binary files differ)
Property changes on: grass/trunk/raster/r.viewshed/sweep1.png
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: grass/trunk/raster/r.viewshed/sweep2.png
===================================================================
(Binary files differ)
Property changes on: grass/trunk/raster/r.viewshed/sweep2.png
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: grass/trunk/raster/r.viewshed/viewshed.cc
===================================================================
--- grass/trunk/raster/r.viewshed/viewshed.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/viewshed.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,542 @@
+/****************************************************************************
+ *
+ * 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 <stdlib.h>
+#include <stdio.h>
+
+#include "viewshed.h"
+#include "visibility.h"
+#include "eventlist.h"
+#include "statusstructure.h"
+
+#ifdef __GRASS
+#include "grass.h"
+#endif
+
+
+#define VIEWSHEDDEBUG if(0)
+
+
+
+
+/* ------------------------------------------------------------ */
+/* return the memory usage (in bytes) of viewshed */
+long long get_viewshed_memory_usage(GridHeader* hd) {
+
+
+ assert(hd);
+ /* the output visibility grid */
+ long long totalcells = (long long)hd->nrows * (long long)hd->ncols;
+ VIEWSHEDDEBUG {printf("rows=%d, cols=%d, total = %lld\n", hd->nrows, hd->ncols, totalcells);}
+ long long gridMemUsage = totalcells * sizeof(float);
+ VIEWSHEDDEBUG {printf("grid usage=%lld\n", gridMemUsage);}
+
+ /* the event array */
+ long long eventListMemUsage = totalcells * 3 * sizeof(AEvent);
+ VIEWSHEDDEBUG {printf("memory_usage: eventList=%lld\n", 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("get_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));
+
+ return (gridMemUsage + eventListMemUsage + dataMemUsage);
+
+}
+
+
+/* ------------------------------------------------------------ */
+void
+print_viewshed_timings(Rtimer initEventTime,
+ Rtimer sortEventTime, Rtimer sweepTime)
+{
+
+ char timeused[1000];
+
+ printf("sweep timings:\n");
+ rt_sprint_safe(timeused, initEventTime);
+ printf("\t%30s", "init events: ");
+ printf(timeused);
+ printf("\n");
+
+ rt_sprint_safe(timeused, sortEventTime);
+ printf("\t%30s", "sort events: ");
+ printf(timeused);
+ printf("\n");
+
+ rt_sprint_safe(timeused, sweepTime);
+ printf("\t%30s", "process events: ");
+ printf(timeused);
+ printf("\n");
+ fflush(stdout);
+ return;
+}
+
+
+/* ------------------------------------------------------------ */
+static void print_event(StatusNode sn)
+{
+ printf("processing (row=%d, col=%d, elev=%f, dist=%f, grad=%f)",
+ sn.row, sn.col, sn.elev, sn.dist2vp, sn.gradient);
+ return;
+}
+
+
+
+/* ------------------------------------------------------------ */
+/* allocates the eventlist array used by kreveled_in_memory */
+AEvent*
+allocate_eventlist(GridHeader* hd) {
+
+ AEvent* eventList;
+
+ long long totalsize = hd->ncols * hd->nrows * 3;
+ totalsize *= sizeof(AEvent);
+ printf("total size of eventlist is %lld B (%d MB); ",
+ totalsize, (int)(totalsize>>20));
+ printf("size_t is %lu B\n", sizeof(size_t));
+
+ /* checking whether allocating totalsize causes an overflow */
+ long long maxsizet = ((long long)1<<(sizeof(size_t)*8)) -1;
+ printf("max size_t is %lld\n", maxsizet);
+ 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);
+ }
+
+ printf("allocating..");
+#ifdef __GRASS__
+ eventList = (AEvent *) G_malloc(totalsize);
+#else
+ eventList = (AEvent *) malloc(totalsize*sizeof(char));
+#endif
+ assert(eventList);
+ printf("..ok\n");
+
+ return eventList;
+}
+
+
+
+
+/*///////////////////////////////////////////////////////////
+ ------------------------------------------------------------ run
+ Viewshed's sweep algorithm on the grid stored in the given file, and
+ with the given viewpoint. Create a visibility grid and return
+ it. The computation runs in memory, which means the input grid, the
+ status structure and the output grid are stored in arrays in
+ memory.
+
+
+ The output: A cell x in the visibility grid is recorded as follows:
+
+ if it is NODATA, then x is set to NODATA
+ if it is invisible, then x is set to INVISIBLE
+ if it is visible, then x is set to the vertical angle wrt to viewpoint
+
+*/
+MemoryVisibilityGrid *viewshed_in_memory(char* inputfname, GridHeader * hd,
+ Viewpoint*vp,ViewOptions viewOptions){
+
+ assert(inputfname && hd && vp);
+ printf("Start sweeping.\n");
+ fflush(stdout);
+
+ /* ------------------------------ */
+ /* create the visibility grid */
+ MemoryVisibilityGrid *visgrid;
+ visgrid = create_inmem_visibilitygrid(*hd, *vp);
+ /* 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",
+ 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;
+ size_t nevents;
+
+ Rtimer initEventTime;
+ rt_start(initEventTime);
+
+ AEvent *eventList = allocate_eventlist(hd);
+#ifdef __GRASS__
+ nevents = grass_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);
+
+ /* ------------------------------ */
+ /*sort the events radially by angle */
+ Rtimer sortEventTime;
+ rt_start(sortEventTime);
+ printf("sorting events..");
+ fflush(stdout);
+
+ /*this is recursive and seg faults for large arrays
+ //qsort(eventList, nevents, sizeof(AEvent), radial_compare_events);
+
+ //this is too slow...
+ //heapsort(eventList, nevents, sizeof(AEvent), radial_compare_events);
+
+ //iostream quicksort */
+ RadialCompare cmpObj;
+ quicksort(eventList, nevents, cmpObj);
+ printf("done\n"); fflush(stdout);
+ rt_stop(sortEventTime);
+
+
+
+
+ /* ------------------------------ */
+ /*create the status structure */
+ StatusList *status_struct = create_status_struct();
+
+ /*Put cells that are initially on the sweepline into status structure */
+ Rtimer sweepTime;
+ StatusNode sn;
+
+ rt_start(sweepTime);
+ for (dimensionType i = vp->col + 1; i < hd->ncols; i++) {
+ sn.col = i;
+ sn.row = vp->row;
+ sn.elev = data[i];
+ if (!is_nodata(visgrid->grid->hd, sn.elev) &&
+ !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_event(sn);
+ printf("\n");
+ }
+ /*insert sn into the status structure */
+ insert_into_status_struct(sn, status_struct);
+ }
+ }
+#ifdef __GRASS__
+ G_free(data);
+#else
+ free(data);
+#endif
+
+
+
+ /* ------------------------------ */
+ /*sweep the event list */
+ long nvis = 0; /*number of visible cells */
+ AEvent *e;
+
+ for (size_t i = 0; i < nevents; i++) {
+
+ /*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;
+
+ /*calculate Distance to VP and Gradient */
+ calculate_dist_n_gradient(&sn, vp);
+ VIEWSHEDDEBUG {
+ printf("next event: ");
+ print_event(sn);
+ }
+
+ switch (e->eventType) {
+ case ENTERING_EVENT:
+ /*insert node into structure */
+ VIEWSHEDDEBUG {
+ printf("..ENTER-EVENT: insert\n");
+ }
+ insert_into_status_struct(sn, status_struct);
+ break;
+
+ case EXITING_EVENT:
+ /*delete node out of status structure */
+ VIEWSHEDDEBUG {
+ printf("..EXIT-EVENT: delete\n");
+ }
+ delete_from_status_struct(status_struct, sn.dist2vp);
+ break;
+
+ case CENTER_EVENT:
+ VIEWSHEDDEBUG {
+ printf("..QUERY-EVENT: query\n");
+ }
+ /*calculate visibility */
+ double max;
+ max=find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
+
+ /*the point is visible: store its vertical angle */
+ if (max <= sn.gradient) {
+ 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);
+ /* when you write the visibility grid you assume that
+ visible values are positive */
+ nvis++;
+ }
+ //else {
+ /* cell is invisible */
+ /* the visibility grid is initialized all invisible */
+ //visgrid->grid->grid_data[sn.row][sn.col] = INVISIBLE;
+ //}
+ break;
+ }
+ }
+ rt_stop(sweepTime);
+
+ printf("Sweeping done.\n");
+ printf("Total cells %ld, visible cells %ld (%.1f percent).\n",
+ (long)visgrid->grid->hd->nrows * visgrid->grid->hd->ncols,
+ nvis,
+ (float)((float)nvis * 100 /
+ (float)(visgrid->grid->hd->nrows *
+ visgrid->grid->hd->ncols)));
+
+ print_viewshed_timings(initEventTime, sortEventTime, sweepTime);
+
+ /*cleanup */
+#ifdef __GRASS__
+ G_free(eventList);
+#else
+ free(eventList);
+#endif
+
+ return visgrid;
+
+}
+
+
+
+
+
+
+
+
+/*///////////////////////////////////////////////////////////
+ ------------------------------------------------------------
+ run Viewshed's algorithm on the grid stored in the given file, and
+ with the given viewpoint. Create a visibility grid and return it. It
+ runs in external memory, i.e. the input grid and the outpt grid are
+ stored as streams
+ */
+
+IOVisibilityGrid *viewshed_external(char* inputfname, GridHeader * hd,
+ Viewpoint* vp, ViewOptions viewOptions){
+
+ assert(inputfname && hd && vp);
+ printf("Start sweeping.\n");
+ fflush(stdout);
+
+
+ /* ------------------------------ */
+ /*initialize the visibility grid */
+ IOVisibilityGrid *visgrid;
+ visgrid = init_io_visibilitygrid(*hd, *vp);
+
+
+ /* ------------------------------ */
+ /* construct the event list corresponding to the give input file and
+ viewpoint; this creates an array of all the cells on
+ the same row as the viewpoint */
+ Rtimer initEventTime, sortEventTime, sweepTime;
+ AMI_STREAM < AEvent > *eventList;
+ double *data;
+ rt_start(initEventTime);
+#ifdef __GRASS__
+ eventList = grass_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);
+ /*printf("Event stream length: %lu\n", (unsigned long)eventList->stream_len()); */
+
+
+ /* ------------------------------ */
+ /*sort the events radially by angle */
+ rt_start(sortEventTime);
+ sort_event_list(&eventList);
+ eventList->seek(0); /*this does not seem to be ensured by sort?? */
+ rt_stop(sortEventTime);
+
+
+ /* ------------------------------ */
+ /*create the status structure */
+ StatusList *status_struct = create_status_struct();
+
+ /* Put cells that are initially on the sweepline into status
+ structure */
+ StatusNode sn;
+ rt_start(sweepTime);
+ for (dimensionType i = vp->col + 1; i < hd->ncols; i++) {
+ sn.col = i;
+ sn.row = vp->row;
+ sn.elev = data[i];
+ if (!is_nodata(visgrid->hd, sn.elev) &&
+ !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_event(sn);
+ printf("\n");
+ }
+ /*insert sn into the status structure */
+ insert_into_status_struct(sn, status_struct);
+ }
+ }
+#ifdef __GRASS__
+ G_free(data);
+#else
+ free(data);
+#endif
+
+
+ /* ------------------------------ */
+ /*sweep the event list */
+ long nvis = 0; /*number of visible cells */
+ VisCell viscell;
+ AEvent *e;
+ AMI_err ae;
+ off_t nbEvents = eventList->stream_len();
+
+ /*printf("nbEvents = %ld\n", (long) nbEvents); */
+
+ 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);
+
+ sn.col = e->col;
+ sn.row = e->row;
+ sn.elev = e->elev;
+ /*calculate Distance to VP and Gradient */
+ calculate_dist_n_gradient(&sn, vp);
+ VIEWSHEDDEBUG {
+ printf("next event: ");
+ print_event(sn);
+ }
+
+ switch (e->eventType) {
+ case ENTERING_EVENT:
+ /*insert node into structure */
+ VIEWSHEDDEBUG {
+ printf("..ENTER-EVENT: insert\n");
+ }
+ insert_into_status_struct(sn, status_struct);
+ break;
+
+ case EXITING_EVENT:
+ /*delete node out of status structure */
+ VIEWSHEDDEBUG {
+ printf("..EXIT-EVENT: delete\n");
+ }
+ delete_from_status_struct(status_struct, sn.dist2vp);
+ break;
+
+ case CENTER_EVENT:
+ VIEWSHEDDEBUG {
+ printf("..QUERY-EVENT: query\n");
+ }
+ /*calculate visibility */
+ viscell.row = sn.row;
+ viscell.col = sn.col;
+ double max;
+ max=find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
+
+ /*the point is visible */
+ if (max <= sn.gradient) {
+ viscell.angle = get_vertical_angle(*vp, sn, viewOptions.doCurv);
+ assert(viscell.angle >= 0);
+ /* viscell.vis = VISIBLE; */
+ add_result_to_io_visibilitygrid(visgrid, &viscell);
+ nvis++;
+ }
+ else {
+ /* else the cell is invisible; we do not write it to the
+ visibility stream because we only record visible cells, and
+ nodata cells; */
+ /* viscell.vis = INVISIBLE; */
+ /* add_result_to_io_visibilitygrid(visgrid, &viscell); */
+ }
+ break;
+ }
+ } /* for each event */
+ rt_stop(sweepTime);
+
+ printf("Sweeping done.\n");
+ printf("Total cells %ld, visible cells %ld (%.1f percent).\n",
+ (long)visgrid->hd->nrows * visgrid->hd->ncols,
+ nvis,
+ (float)((float)nvis * 100 /
+ (float)(visgrid->hd->nrows * visgrid->hd->ncols)));
+
+ print_viewshed_timings(initEventTime, sortEventTime, sweepTime);
+
+ /*cleanup */
+ delete eventList;
+
+ return visgrid;
+}
Property changes on: grass/trunk/raster/r.viewshed/viewshed.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/viewshed.h
===================================================================
--- grass/trunk/raster/r.viewshed/viewshed.h (rev 0)
+++ grass/trunk/raster/r.viewshed/viewshed.h 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,110 @@
+
+/****************************************************************************
+ *
+ * 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 _KREVELD_H
+#define _KREVELD_H
+
+#include "visibility.h"
+#include "grid.h"
+#include "eventlist.h"
+
+#ifdef __GRASS__
+#include "grass.h"
+#include <grass/iostream/ami.h>
+#else
+#include <ami.h>
+#endif
+
+
+/* ------------------------------------------------------------ */
+/*return the memory usage in bytes of the algorithm when running in
+ memory*/
+long long get_viewshed_memory_usage(GridHeader* hd);
+
+
+
+
+/* ------------------------------------------------------------ */
+/* run Kreveld's algorithm on the grid stored in the given file, and
+ with the given viewpoint. Create a visibility grid and return
+ it. It runs in memory, i.e. the input grid and output grid are
+ stored in 2D arrays in memory. The computation runs in memory,
+ which means the input grid, the status structure and the output
+ grid are stored in arrays in memory.
+
+ The output: A cell x in the visibility grid is recorded as follows:
+
+ if it is NODATA, then x is set to NODATA if it is invisible, then x
+ is set to INVISIBLE if it is visible, then x is set to the vertical
+ angle wrt to viewpoint
+
+*/
+MemoryVisibilityGrid *viewshed_in_memory(char* inputfname,
+ GridHeader * hd,
+ Viewpoint * vp,
+ ViewOptions viewOptions);
+
+
+
+
+/* ------------------------------------------------------------ */
+/* compute viewshed on the grid stored in the given file, and with the
+ given viewpoint. Create a visibility grid and return it. The
+ program runs in external memory.
+
+ The output: A cell x in the visibility grid is recorded as follows:
+
+ if it is NODATA, then x is set to NODATA if it is invisible, then x
+ is set to INVISIBLE if it is visible, then x is set to the vertical
+ angle wrt to viewpoint.
+
+*/
+IOVisibilityGrid *viewshed_external(char* inputfname,
+ GridHeader * hd,
+ Viewpoint * vp,
+ ViewOptions viewOptions);
+
+
+
+void print_viewshed_timings(Rtimer initEventTime, Rtimer sortEventTime,
+ Rtimer sweepTime);
+
+
+
+#endif
Property changes on: grass/trunk/raster/r.viewshed/viewshed.h
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/visibility.cc
===================================================================
--- grass/trunk/raster/r.viewshed/visibility.cc (rev 0)
+++ grass/trunk/raster/r.viewshed/visibility.cc 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,564 @@
+/****************************************************************************
+ *
+ * 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 <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"
+#include "grass.h"
+
+
+
+/* ------------------------------------------------------------ */
+/* viewpoint functions */
+void print_viewpoint(Viewpoint vp)
+{
+ printf("vp=(%d, %d, %.1f) ", vp.row, vp.col, vp.elev);
+ return;
+}
+
+/* ------------------------------------------------------------ */
+void set_viewpoint_coord(Viewpoint * vp, dimensionType row, dimensionType col)
+{
+ assert(vp);
+ vp->row = row;
+ vp->col = col;
+ return;
+}
+
+/* ------------------------------------------------------------ */
+void set_viewpoint_elev(Viewpoint * vp, float elev)
+{
+ assert(vp);
+ vp->elev = elev;
+ return;
+}
+
+/* ------------------------------------------------------------ */
+/*copy from b to a */
+void copy_viewpoint(Viewpoint * a, Viewpoint b)
+{
+ assert(a);
+ a->row = b.row;
+ a->col = b.col;
+ a->elev = b.elev;
+ return;
+}
+
+
+/* ------------------------------------------------------------ */
+/* MemoryVisibilityGrid functions */
+
+/* create and return a grid of the sizes specified in the header */
+MemoryVisibilityGrid *create_inmem_visibilitygrid(GridHeader hd, Viewpoint vp){
+
+ MemoryVisibilityGrid *visgrid;
+#ifdef __GRASS__
+ visgrid = (MemoryVisibilityGrid *) G_malloc(sizeof(MemoryVisibilityGrid));
+#else
+ visgrid = (MemoryVisibilityGrid *) malloc(sizeof(MemoryVisibilityGrid));
+#endif
+ assert(visgrid);
+
+
+ /* create the grid */
+ visgrid->grid = create_empty_grid();
+ 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 */
+ copy_header(visgrid->grid->hd, hd);
+
+ /* allocate the Grid data */
+ 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);
+
+ return visgrid;
+}
+
+
+
+
+/* ------------------------------------------------------------ */
+void free_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid)
+{
+
+ assert(visgrid);
+
+ 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;
+}
+
+
+
+/* ------------------------------------------------------------ */
+/*set all values of visgrid's Grid to the given value*/
+void set_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid, float val)
+{
+
+ assert(visgrid && visgrid->grid && visgrid->grid->hd &&
+ visgrid->grid->grid_data);
+
+ dimensionType i, j;
+
+ for (i = 0; i < visgrid->grid->hd->nrows; i++) {
+ assert(visgrid->grid->grid_data[i]);
+ for (j = 0; j < visgrid->grid->hd->ncols; j++) {
+ visgrid->grid->grid_data[i][j] = val;
+ }
+ }
+ return;
+}
+
+
+
+/* ------------------------------------------------------------ */
+/*set the (i,j) value of visgrid's Grid to the given value*/
+void add_result_to_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid,
+ dimensionType i, dimensionType j,
+ float val)
+{
+
+ assert(visgrid && visgrid->grid && visgrid->grid->hd &&
+ visgrid->grid->grid_data);
+ assert(i < visgrid->grid->hd->nrows);
+ assert(j < visgrid->grid->hd->ncols);
+ assert(visgrid->grid->grid_data[i]);
+
+ visgrid->grid->grid_data[i][j] = val;
+
+ return;
+}
+
+
+
+/* ------------------------------------------------------------ */
+/* The following functions are used to convert the visibility results
+ recorded during the viewshed computation into the output grid into
+ tehe output required by the user.
+
+ x is assumed to be the visibility value computed for a cell during the
+ viewshed computation.
+
+ The value computed during the viewshed is the following:
+
+ x is NODATA if the cell is NODATA;
+
+
+ x is INVISIBLE if the cell is invisible;
+
+ x is the vertical angle of the cell wrt the viewpoint if the cell is
+ visible---the angle is a value in (0,180).
+*/
+int is_visible(float x) {
+ /* 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);
+ if (isnull) return 0;
+ else return (x >= 0);
+#else
+ return (x >=0);
+#endif
+}
+int is_invisible_not_nodata(float x) {
+
+ return ( (int)x == (int)INVISIBLE);
+}
+
+int is_invisible_nodata(float x) {
+
+ return (!is_visible(x)) && (!is_invisible_not_nodata(x));
+}
+
+/* ------------------------------------------------------------ */
+/* This function is called when the program runs in
+ viewOptions.outputMode == OUTPUT_BOOL. */
+float booleanVisibilityOutput(float x) {
+ /* NODATA and INVISIBLE are both negative values */
+ if (is_visible(x))
+ return BOOL_VISIBLE;
+ else
+ return BOOL_INVISIBLE;
+}
+/* ------------------------------------------------------------ */
+/* This function is called when the program runs in
+ viewOptions.outputMode == OUTPUT_ANGLE. In this case x represents
+ the right value. */
+float angleVisibilityOutput(float x) {
+
+ return x;
+}
+
+
+
+
+
+
+/* ------------------------------------------------------------ */
+/* visgrid is the structure that records the visibility information
+ after the sweep is done. Use it to write the visibility output
+ grid and then distroy it.
+*/
+void save_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid,
+ ViewOptions viewOptions, Viewpoint vp) {
+
+#ifdef __GRASS__
+ if (viewOptions.outputMode == OUTPUT_BOOL)
+ save_grid_to_GRASS(visgrid->grid, viewOptions.outputfname, CELL_TYPE,
+ booleanVisibilityOutput);
+ else if (viewOptions.outputMode == OUTPUT_ANGLE)
+ save_grid_to_GRASS(visgrid->grid, viewOptions.outputfname, FCELL_TYPE,
+ angleVisibilityOutput);
+ else
+ /* elevation output */
+ 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);
+
+ return;
+}
+
+
+
+
+
+/* ------------------------------------------------------------ */
+/* 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);
+
+ /*stream */
+ visgrid->visStr = new AMI_STREAM < VisCell > ();
+ assert(visgrid->visStr);
+
+ return visgrid;
+}
+
+
+
+/* ------------------------------------------------------------ */
+/*free the grid */
+void free_io_visibilitygrid(IOVisibilityGrid * grid)
+{
+ assert(grid);
+#ifdef __GRASS__
+ if (grid->hd)
+ G_free(grid->hd);
+ if (grid->vp)
+ G_free(grid->vp);
+ if (grid->visStr)
+ 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, VisCell * cell)
+{
+
+ assert(visgrid && cell);
+
+ AMI_err ae;
+
+ assert(visgrid->visStr);
+ ae = visgrid->visStr->write_item(*cell);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ return;
+}
+
+
+
+/* ------------------------------------------------------------ */
+/*compare function, (i,j) grid order */
+int IJCompare::compare(const VisCell & a, const VisCell & b)
+{
+ if (a.row > b.row)
+ return 1;
+ if (a.row < b.row)
+ return -1;
+
+ /*a.row==b.row */
+ if (a.col > b.col)
+ return 1;
+ if (a.col < b.col)
+ return -1;
+ /*all equal */
+ return 0;
+}
+
+
+/* ------------------------------------------------------------ */
+/*sort stream in grid order */
+void sort_io_visibilitygrid(IOVisibilityGrid * visGrid)
+{
+
+ assert(visGrid);
+ assert(visGrid->visStr);
+ if (visGrid->visStr->stream_len() == 0)
+ return;
+
+ AMI_STREAM < VisCell > *sortedStr;
+ AMI_err ae;
+ IJCompare cmpObj;
+
+ ae = AMI_sort(visGrid->visStr, &sortedStr, &cmpObj, 1);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ assert(sortedStr);
+ sortedStr->seek(0);
+
+ visGrid->visStr = sortedStr;
+ return;
+}
+
+
+
+/* ------------------------------------------------------------ */
+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);
+
+ else if (viewOptions.outputMode == OUTPUT_ANGLE)
+ save_io_visibilitygrid_to_GRASS(visgrid, viewOptions.outputfname,
+ FCELL_TYPE, angleVisibilityOutput);
+ else
+ /* elevation output */
+ 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;
+}
+
+
Property changes on: grass/trunk/raster/r.viewshed/visibility.cc
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/raster/r.viewshed/visibility.h
===================================================================
--- grass/trunk/raster/r.viewshed/visibility.h (rev 0)
+++ grass/trunk/raster/r.viewshed/visibility.h 2008-08-13 19:57:41 UTC (rev 32752)
@@ -0,0 +1,255 @@
+/****************************************************************************
+ *
+ * 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 visibility_h
+#define visibility_h
+
+#ifdef __GRASS__
+#include <grass/config.h>
+#include <grass/iostream/ami.h>
+
+#else
+#include <ami.h>
+#endif
+
+#include "grid.h"
+
+
+
+/* default max distance */
+#define INFINITY_DISTANCE -1
+
+
+
+typedef struct viewpoint_ {
+ dimensionType row, col;
+ float elev;
+} Viewpoint;
+
+
+typedef enum {
+ VISIBLE = 1,
+ INVISIBLE = -1,
+
+ /*boolean values for output */
+ BOOL_VISIBLE = 1,
+ BOOL_INVISIBLE = 0
+} VisMode;
+
+
+typedef struct visCell_ {
+ dimensionType row;
+ dimensionType col;
+ /* VisMode vis; */
+ float angle;
+} VisCell;
+
+
+
+typedef enum outputMode_{
+ OUTPUT_ANGLE = 0,
+ OUTPUT_BOOL = 1,
+ OUTPUT_ELEV = 2
+} OutputMode;
+
+
+typedef struct viewOptions_ {
+
+ /* the name of the input raster */
+ char inputfname [100];
+
+ /* the name of the output raster */
+ char outputfname[100];
+
+ float obsElev;
+ /* observer elevation above the terrain */
+
+ float maxDist;
+ /* points that are farther than this distance from the viewpoint are
+ not visible */
+
+ OutputMode outputMode;
+ /* The mode the viewshed is output;
+
+ in angle mode, the values recorded are {NODATA, INVISIBLE, angle}
+
+ in boolean mode, the values recorded are {BOOL_INVISIBLE, BOOL_VISIBLE}
+
+ in elev mode, the values recorded are {NODATA, INVISIBLE, elevation}
+ */
+
+ int doCurv;
+ /*determines if the curvature of the earth should be considered
+ when calculating. Only implemented for GRASS version. */
+
+ double ellps_a; /* the parameter of teh ellipsoid */
+ float cellsize; /* the cell resolution */
+} ViewOptions;
+
+
+
+
+/*memory visibility grid */
+typedef struct memory_visibility_grid_ {
+ Grid *grid;
+ Viewpoint *vp;
+} MemoryVisibilityGrid;
+
+
+/*io-efficient visibility grid */
+typedef struct IOvisibility_grid_ {
+ GridHeader *hd;
+ Viewpoint *vp;
+ AMI_STREAM < VisCell > *visStr;
+} IOVisibilityGrid;
+
+
+
+
+/* ------------------------------------------------------------ */
+/* visibility output functions */
+
+/* The following functions are used to convert the visibility results
+ recorded during the viewshed computation into the output grid into
+ the format required by the user. x is assumed to be the
+ visibility angle computed for a cell during the viewshed
+ computation.
+
+ The value passed to this function is the following: x is NODATA if the
+ cell is NODATA; x is INVISIBLE if the cell is invisible; x is the
+ vertical angle of the cell wrt the viewpoint if the cell is
+ visible---the angle is a value in (0,180).
+*/
+/* these functions assume that x is a value computed during the
+viewshed computation; right now x represents the vertical angle of a
+visible point wrt to the viewpoint; INVISIBLE if invisible; NODATA if
+nodata. They return true if x is visible, invisible but nodata,
+andnodata, respectively */
+int is_visible(float x);
+int is_invisible_not_nodata(float x);
+int is_invisible_nodata(float x);
+
+/* This function is called when the program runs in
+ viewOptions.outputMode == OUTPUT_BOOL. */
+float booleanVisibilityOutput(float x);
+
+/* This function is called when the program runs in
+ viewOptions.outputMode == OUTPUT_ANGLE. */
+float angleVisibilityOutput(float x);
+
+
+
+
+
+/* ------------------------------------------------------------ */
+/* viewpoint functions */
+
+void print_viewpoint(Viewpoint vp);
+
+/*copy from b to a */
+void copy_viewpoint(Viewpoint * a, Viewpoint b);
+
+void
+set_viewpoint_coord(Viewpoint * vp, dimensionType row, dimensionType col);
+
+void set_viewpoint_elev(Viewpoint * vp, float elev);
+
+
+
+/* ------------------------------------------------------------ */
+/* MemoryVisibilityGrid functions */
+
+MemoryVisibilityGrid *create_inmem_visibilitygrid(GridHeader hd, Viewpoint vp);
+
+void free_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid);
+
+void set_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid, float val);
+
+void add_result_to_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid,
+ dimensionType i, dimensionType j,
+ float val);
+
+void save_inmem_visibilitygrid(MemoryVisibilityGrid * vigrid,
+ ViewOptions viewopt, Viewpoint vp);
+
+
+/* ------------------------------------------------------------ */
+/* IOVisibilityGrid functions */
+
+/*create grid from given header and viewpoint */
+IOVisibilityGrid *init_io_visibilitygrid(GridHeader hd, Viewpoint vp);
+
+/*frees a visibility grid */
+void free_io_visibilitygrid(IOVisibilityGrid * grid);
+
+/*write cell to stream */
+void add_result_to_io_visibilitygrid(IOVisibilityGrid * visgrid,
+ VisCell * cell);
+
+/*void
+ addResult(IOVisibilityGrid* visgrid, DimensionType row, DimensionType col,
+ VisMode vis);
+*/
+
+
+/* write visibility grid. assume all cells that are not in stream are
+ NOT visible. assume stream is sorted. */
+void
+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);
+
+class IJCompare
+{
+ public:
+ int compare(const VisCell &, const VisCell &);
+};
+
+
+
+#endif
Property changes on: grass/trunk/raster/r.viewshed/visibility.h
___________________________________________________________________
Name: svn:executable
+ *
More information about the grass-commit
mailing list