First working version. Uses grass/lib/iostrean, Makefile and include changes from Paul kelly

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@32752 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.viewshed/Makefile

@@ -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)

raster/r.viewshed/distribute.h

@@ -0,0 +1,184 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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

raster/r.viewshed/eventlist.cc

@@ -0,0 +1,961 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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;
+}

raster/r.viewshed/eventlist.h

@@ -0,0 +1,170 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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

raster/r.viewshed/grass.cc

@@ -0,0 +1,822 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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

raster/r.viewshed/grass.h

@@ -0,0 +1,151 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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__*/

raster/r.viewshed/grid.cc

@@ -0,0 +1,397 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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;
+}

raster/r.viewshed/grid.h

@@ -0,0 +1,119 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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

raster/r.viewshed/main.cc

@@ -0,0 +1,863 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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(&region) == -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,&region);
+#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, &region);
+    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");
+}

raster/r.viewshed/print_message.cc

@@ -0,0 +1,64 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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
+}

raster/r.viewshed/print_message.h

@@ -0,0 +1,46 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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

raster/r.viewshed/r.viewshed.html

@@ -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@bowdoin.edu</tt>
+
+<p> Yi Zhuang (Carnegie-Mellon University): <tt>yzhuang@andrew.cmu.edu</tt>
+
+<p>William Richard (Bowdoin College): <tt>willster3021@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>&copy; 2003-2008 <a href="http://grass.osgeo.org">GRASS Development Team</a></p>
+</body>
+</html>

raster/r.viewshed/rbbst.cc

@@ -0,0 +1,764 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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;
+}

raster/r.viewshed/rbbst.h

@@ -0,0 +1,156 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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

raster/r.viewshed/statusstructure.cc

@@ -0,0 +1,214 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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);
+}
+

raster/r.viewshed/statusstructure.h

@@ -0,0 +1,107 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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

raster/r.viewshed/viewshed.cc

@@ -0,0 +1,542 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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;
+}

raster/r.viewshed/viewshed.h

@@ -0,0 +1,110 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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

raster/r.viewshed/visibility.cc

@@ -0,0 +1,564 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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;
+}
+
+

raster/r.viewshed/visibility.h

@@ -0,0 +1,255 @@
+/****************************************************************************
+ *
+ * MODULE:       r.viewshed
+ *
+ * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
+ *               Yi Zhuang - yzhuang@bowdoin.edu
+
+ *               Ported to GRASS by William Richard -
+ *               wkrichar@bowdoin.edu or willster3021@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