r.viewshed deleted; moved to addons

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

raster/r.viewshed/Makefile

@@ -1,52 +0,0 @@
-# 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

@@ -1,184 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,961 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,170 +0,0 @@
-
-/****************************************************************************
- *
- * 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

@@ -1,822 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,151 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,397 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,119 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,863 +0,0 @@
-
-/****************************************************************************
- *
- * 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

@@ -1,64 +0,0 @@
-
-/****************************************************************************
- *
- * 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

@@ -1,46 +0,0 @@
-
-/****************************************************************************
- *
- * 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

@@ -1,225 +0,0 @@
-<!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

@@ -1,764 +0,0 @@
-
-/****************************************************************************
- *
- * 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

@@ -1,156 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,214 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,107 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,542 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,110 +0,0 @@
-
-/****************************************************************************
- *
- * 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

@@ -1,564 +0,0 @@
-/****************************************************************************
- *
- * 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

@@ -1,255 +0,0 @@
-/****************************************************************************
- *
- * 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