[GRASS-SVN] r62081 - grass-addons/grass6/raster/r.viewshed

svn_grass at osgeo.org svn_grass at osgeo.org
Thu Sep 25 03:44:26 PDT 2014


Author: neteler
Date: 2014-09-25 03:44:26 -0700 (Thu, 25 Sep 2014)
New Revision: 62081

Added:
   grass-addons/grass6/raster/r.viewshed/description.html
Removed:
   grass-addons/grass6/raster/r.viewshed/r.viewshed.html
Modified:
   grass-addons/grass6/raster/r.viewshed/distribute.cpp
   grass-addons/grass6/raster/r.viewshed/grass.cpp
   grass-addons/grass6/raster/r.viewshed/main.cpp
   grass-addons/grass6/raster/r.viewshed/rbbst.cpp
   grass-addons/grass6/raster/r.viewshed/rbbst.h
   grass-addons/grass6/raster/r.viewshed/statusstructure.cpp
   grass-addons/grass6/raster/r.viewshed/viewshed.cpp
Log:
r.viewshed addon: backported fixes from GRASS 7 trunk

Copied: grass-addons/grass6/raster/r.viewshed/description.html (from rev 62070, grass-addons/grass6/raster/r.viewshed/r.viewshed.html)
===================================================================
--- grass-addons/grass6/raster/r.viewshed/description.html	                        (rev 0)
+++ grass-addons/grass6/raster/r.viewshed/description.html	2014-09-25 10:44:26 UTC (rev 62081)
@@ -0,0 +1,240 @@
+<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 elevation
+map name, an output raster map name, and the location of the
+viewpoint.
+
+<p>Viewpoint (<em>coordinate</em> parameter): For the time being the viewpoint
+is assumed to be located inside the terrain.  The viewpoint location is
+given in map coordinates.
+
+<p>Output: The output map may have one of three possible formats,
+based on which flags are set.
+
+<p>
+By default, if no flag is set, the output is in angle-mode, and
+each point in the output map is marked as NULL if the point is not
+visible or the respective point in the elevation map is NULL.
+
+Otherwise, a value in [0, 180] representing the vertical angle with
+regard to the viewpoint, in degrees, if the point is visible.
+A value of 0 is directly below the specified viewing position,
+90 is due horizontal. The angle to the cell containing the viewing
+position is undefined and set to 180. 
+
+<p>
+If the <b>-b</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 no-data/null or not visible
+  <li> 1 if the point is visible.
+</ul>
+
+
+<p>
+If the <b>-e</b> flag is set, the output is in elevation-mode, and each point
+in the output map is marked as:
+<!-- Check & FIXME -->
+<ul>
+  <li> no-data (null), if the respective point in the elevation map is no-data (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>
+
+<p>
+If you wish to identify the area of the map which is within the search
+radius but not visible, a combination of <em>r.buffer</em> and
+<em>r.mapcalc</em> can be used to create a negative of the viewshed map.
+
+
+
+<p>
+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 of 1.75 map units above the terrain.  The user can change this using option
+<em>obs_elev</em>. The value entered is in the same units
+as the elevation.
+
+<p>
+Target elevation: By default the target is assumed to have
+height of 0 map units above the terrain.  The user can change this using option
+<em>tgt_elev</em> to determine if objects of a given height would be 
+visible. The value entered is in the same units as the elevation.
+
+<p>
+Maximum visibility distance: By default there is no restriction on
+the maximum distance to which the observer can see.  The user can set
+a maximum distance of visibility using option <em>max_dist</em>.
+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>
+
+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>
+
+<em>r.viewshed</em> uses the following model for determining
+visibility: The height of a cell is assumed to be variable, and the 
+actual height of a point falling into a cell, but not identical the cell 
+center, is interpolated. Thus the terrain is viewed as a smooth surface. 
+Two points are visible to each other if their line-of-sight does not
+intersect the terrain. The height for an arbitrary point x in the terrain 
+is interpolated from the 4 surrounding neighbours. This means that this 
+model does a bilinear interpolation of heights.
+
+This model is suitable for both low and high resolution rasters as well 
+as terrain with flat and steep slopes.
+
+<p>The core of the algorithm is determining, for each cell, the
+line-of-sight and its intersections with the cells in the terrain. For
+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>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.png" width="200" alt="[SDF]" border="0"></th>
+      <th><img src="sweep2.png" width="200" alt="[SDF]" border="0"></th>
+  </tr>
+  <tr>	
+    <th>The sweep-line.</th>
+    <th>The active cells.</th>
+  </tr> 
+</table>
+
+
+<h2>EXAMPLES</h2>
+
+Using the North Carolina dataset:  Compute viewshed from a observation 
+point (coordinates: 638728.087167, 220609.261501) which is 5 meters 
+above ground:
+
+<div class="code"><pre>
+g.region rast=elev_lid792_1m -p
+r.viewshed input=elev_lid792_1m output=elev_lid792_1m_viewshed \
+  coordinate=638728.087167,220609.261501 obs_elev=5.0
+</pre></div>
+
+<div align="center" style="margin: 10px">
+<img src="r_viewshed.png" border=0><br>
+<i>Viewshed</i>
+</div>
+<!--
+# image generated using
+d.mon x0
+d.rast elev_lid792_1m
+d.rast elev_lid792_1m_viewshed -o
+# save image to files
+# crop the background using Gimp or ImageMagic
+mogrify -trim *.png
+# some bounding box problems noticed when opening mogrify result in Gimp
+-->
+
+Using the Spearfish dataset:  calculating the viewpoint from top
+of a mountain:
+
+<div class="code"><pre>
+g.region rast=elevation.10m
+r.viewshed input=elevation.10m output=viewshed \
+  coordinate=598869,4916642 mem=800
+</pre></div>
+
+
+<h2>REFERENCES</h2>
+
+<ul>
+
+<li>Computing Visibility on Terrains in External Memory. Herman 
+Haverkort, Laura Toma and Yi Zhuang. In <em>ACM Journal on Experimental 
+Algorithmics (JEA)</em> 13 (2009).</li>
+	 
+<li><a href="http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.76.4282&rep=rep1&type=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>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.los.html">r.los</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>
+</em>
+
+<h2>AUTHORS</h2>
+
+<p>Laura Toma (Bowdoin College): <tt>ltoma at bowdoin.edu</tt>
+<p> Yi Zhuang (Carnegie-Mellon University): <tt>yzhuang at andrew.cmu.edu</tt>
+<p>William Richard (Bowdoin College): <tt>willster3021 at gmail.com </tt>
+<p>Markus Metz
+
+<p>
+<i>Last changed: $Date$</i>

Modified: grass-addons/grass6/raster/r.viewshed/distribute.cpp
===================================================================
--- grass-addons/grass6/raster/r.viewshed/distribute.cpp	2014-09-25 08:12:03 UTC (rev 62080)
+++ grass-addons/grass6/raster/r.viewshed/distribute.cpp	2014-09-25 10:44:26 UTC (rev 62081)
@@ -52,8 +52,8 @@
 #include <grass/iostream/ami.h>
 
 
-/*note: MAX_STREAM_OPEN defined in IOStream/include/ami_stream.h, which
-   is included by ami.h  */
+/*note: MAX_STREAMS_OPEN defined in include/iostream/ami_stream.h
+  which is included by ami.h  */
 
 #include "distribute.h"
 #include "viewshed.h"

Modified: grass-addons/grass6/raster/r.viewshed/grass.cpp
===================================================================
--- grass-addons/grass6/raster/r.viewshed/grass.cpp	2014-09-25 08:12:03 UTC (rev 62080)
+++ grass-addons/grass6/raster/r.viewshed/grass.cpp	2014-09-25 10:44:26 UTC (rev 62081)
@@ -108,8 +108,11 @@
 	hd->nrows = (dimensionType) nrows;
 	hd->ncols = (dimensionType) ncols;
     }
-    else
-	G_fatal_error(_("Grid dimension too big for current precision"));
+    else {
+	G_warning("ERROR: nrows (%d) > maxDimension (%d) AND/OR ncols (%d) > maxDimension (%d)", 
+	           nrows, maxDimension, ncols, maxDimension);
+	G_fatal_error(_("Computational region too large. Use smaller area or lower raster resolution"));
+    }
 
 
     /*fill in rest of header */
@@ -127,7 +130,7 @@
     hd->ns_res = region->ns_res;
     //store the null value of the map
     G_set_null_value(&(hd->nodata_value), 1, G_SURFACE_TYPE);
-    G_message("Nodata value set to %f", hd->nodata_value);
+    G_verbose_message("Nodata value set to %f", hd->nodata_value);
     
     
     
@@ -228,19 +231,14 @@
     G_set_null_value(inrast[1], ncols, data_type);
     G_set_null_value(inrast[2], ncols, data_type);
 
-    /*DCELL to check for loss of prescion- haven't gotten that to work
-       yet though */
-    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;
-    int row1, col1;
+    dimensionType i, j;
     double ax, ay;
-    G_SURFACE_T *elev1, *elev2, *elev3, *elev4;
     AEvent e;
     
     /* read first row */
@@ -249,7 +247,6 @@
     e.angle = -1;
     for (i = 0; i < nrows; i++) {
 	/*read in the raster row */
-	int rasterRowResult = 1;
 	
 	G_SURFACE_T *tmprast = inrast[0];
 	inrast[0] = inrast[1];
@@ -257,14 +254,10 @@
 	inrast[2] = tmprast;
 
 	if (i < nrows - 1)
-	    rasterRowResult = G_get_raster_row(infd, inrast[2], i + 1, data_type);
+	    G_get_raster_row(infd, inrast[2], i + 1, data_type);
 	else
 	    G_set_null_value(inrast[2], ncols, data_type);
 
-	if (rasterRowResult <= 0)
-	    G_fatal_error(_("Coord not read from row %d of <%s>"), i,
-			  rastName);
-
 	G_percent(i, nrows, 2);
 
 	/*fill event list with events from this row */
@@ -405,7 +398,7 @@
 					     IOVisibilityGrid * visgrid)
 {
 
-    G_message(_("Computing events ..."));
+    G_message(_("Computing events..."));
     assert(rastName && vp && hd && visgrid);
 
     if (data != NULL) {
@@ -458,12 +451,9 @@
     G_set_null_value(inrast[2], ncols, data_type);
 
     /*scan through the raster data */
-    DCELL d;
     int isnull = 0;
-    dimensionType i, j, k;
-    int row1, col1;
+    dimensionType i, j;
     double ax, ay;
-    G_SURFACE_T *elev1, *elev2, *elev3, *elev4;
     AEvent e;
 
     /* read first row */
@@ -477,7 +467,6 @@
 	G_percent(i, nrows, 2);
 
 	/*read in the raster row */
-	int rasterRowResult = 1;
 	
 	G_SURFACE_T *tmprast = inrast[0];
 	inrast[0] = inrast[1];
@@ -485,14 +474,10 @@
 	inrast[2] = tmprast;
 
 	if (i < nrows - 1)
-	    rasterRowResult = G_get_raster_row(infd, inrast[2], i + 1, data_type);
+	    G_get_raster_row(infd, inrast[2], i + 1, data_type);
 	else
 	    G_set_null_value(inrast[2], ncols, data_type);
 
-	if (rasterRowResult <= 0)
-	    G_fatal_error(_("Coord not read from row %d of <%s>"), i,
-			  rastName);
-
 	/*fill event list with events from this row */
 	for (j = 0; j < ncols; j++) {
 
@@ -632,7 +617,7 @@
 		   float (*fun) (float))
 {
 
-    G_message(_("Saving grid to <%s>"), filename);
+    G_important_message(_("Writing output raster map..."));
     assert(grid && filename);
 
     /*open the new raster  */

Modified: grass-addons/grass6/raster/r.viewshed/main.cpp
===================================================================
--- grass-addons/grass6/raster/r.viewshed/main.cpp	2014-09-25 08:12:03 UTC (rev 62080)
+++ grass-addons/grass6/raster/r.viewshed/main.cpp	2014-09-25 10:44:26 UTC (rev 62081)
@@ -113,7 +113,7 @@
     /*initialize module */
     module = G_define_module();
     module->keywords = _("raster, viewshed, line of sight");
-    module->description = _("IO-efficient viewshed algorithm");
+    module->description = _("Computes the viewshed of a point on an elevation raster map.");
 
     struct Cell_head region;
     char *optstreamdir;
@@ -292,12 +292,13 @@
 	}
 	G_important_message(_("Intermediate files will not be deleted \
 		              in case of abnormal termination."));
+	G_important_message(_("Intermediate location: %s"), viewOptions.streamdir);
 	G_important_message(_("To save space delete these files manually!"));
 
 
 	/* initialize IOSTREAM memory manager */
 	MM_manager.set_memory_limit(memSizeBytes);
-	MM_manager.warn_memory_limit();
+	MM_manager.ignore_memory_limit();
 	MM_manager.print_limit_mode();
 
 

Deleted: grass-addons/grass6/raster/r.viewshed/r.viewshed.html
===================================================================
--- grass-addons/grass6/raster/r.viewshed/r.viewshed.html	2014-09-25 08:12:03 UTC (rev 62080)
+++ grass-addons/grass6/raster/r.viewshed/r.viewshed.html	2014-09-25 10:44:26 UTC (rev 62081)
@@ -1,240 +0,0 @@
-<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 elevation
-map name, an output raster map name, and the location of the
-viewpoint.
-
-<p>Viewpoint (<em>coordinate</em> parameter): For the time being the viewpoint
-is assumed to be located inside the terrain.  The viewpoint location is
-given in map coordinates.
-
-<p>Output: The output map may have one of three possible formats,
-based on which flags are set.
-
-<p>
-By default, if no flag is set, the output is in angle-mode, and
-each point in the output map is marked as NULL if the point is not
-visible or the respective point in the elevation map is NULL.
-
-Otherwise, a value in [0, 180] representing the vertical angle with
-regard to the viewpoint, in degrees, if the point is visible.
-A value of 0 is directly below the specified viewing position,
-90 is due horizontal. The angle to the cell containing the viewing
-position is undefined and set to 180. 
-
-<p>
-If the <b>-b</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 no-data/null or not visible
-  <li> 1 if the point is visible.
-</ul>
-
-
-<p>
-If the <b>-e</b> flag is set, the output is in elevation-mode, and each point
-in the output map is marked as:
-<!-- Check & FIXME -->
-<ul>
-  <li> no-data (null), if the respective point in the elevation map is no-data (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>
-
-<p>
-If you wish to identify the area of the map which is within the search
-radius but not visible, a combination of <em>r.buffer</em> and
-<em>r.mapcalc</em> can be used to create a negative of the viewshed map.
-
-
-
-<p>
-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 of 1.75 map units above the terrain.  The user can change this using option
-<em>obs_elev</em>. The value entered is in the same units
-as the elevation.
-
-<p>
-Target elevation: By default the target is assumed to have
-height of 0 map units above the terrain.  The user can change this using option
-<em>tgt_elev</em> to determine if objects of a given height would be 
-visible. The value entered is in the same units as the elevation.
-
-<p>
-Maximum visibility distance: By default there is no restriction on
-the maximum distance to which the observer can see.  The user can set
-a maximum distance of visibility using option <em>max_dist</em>.
-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>
-
-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>
-
-<em>r.viewshed</em> uses the following model for determining
-visibility: The height of a cell is assumed to be variable, and the 
-actual height of a point falling into a cell, but not identical the cell 
-center, is interpolated. Thus the terrain is viewed as a smooth surface. 
-Two points are visible to each other if their line-of-sight does not
-intersect the terrain. The height for an arbitrary point x in the terrain 
-is interpolated from the 4 surrounding neighbours. This means that this 
-model does a bilinear interpolation of heights.
-
-This model is suitable for both low and high resolution rasters as well 
-as terrain with flat and steep slopes.
-
-<p>The core of the algorithm is determining, for each cell, the
-line-of-sight and its intersections with the cells in the terrain. For
-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>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.png" width="200" alt="[SDF]" border="0"></th>
-      <th><img src="sweep2.png" width="200" alt="[SDF]" border="0"></th>
-  </tr>
-  <tr>	
-    <th>The sweep-line.</th>
-    <th>The active cells.</th>
-  </tr> 
-</table>
-
-
-<h2>EXAMPLES</h2>
-
-Using the North Carolina dataset:  Compute viewshed from a observation 
-point (coordinates: 638728.087167, 220609.261501) which is 5 meters 
-above ground:
-
-<div class="code"><pre>
-g.region rast=elev_lid792_1m -p
-r.viewshed input=elev_lid792_1m output=elev_lid792_1m_viewshed \
-  coordinate=638728.087167,220609.261501 obs_elev=5.0
-</pre></div>
-
-<div align="center" style="margin: 10px">
-<img src="r_viewshed.png" border=0><br>
-<i>Viewshed</i>
-</div>
-<!--
-# image generated using
-d.mon x0
-d.rast elev_lid792_1m
-d.rast elev_lid792_1m_viewshed -o
-# save image to files
-# crop the background using Gimp or ImageMagic
-mogrify -trim *.png
-# some bounding box problems noticed when opening mogrify result in Gimp
--->
-
-Using the Spearfish dataset:  calculating the viewpoint from top
-of a mountain:
-
-<div class="code"><pre>
-g.region rast=elevation.10m
-r.viewshed input=elevation.10m output=viewshed \
-  coordinate=598869,4916642 mem=800
-</pre></div>
-
-
-<h2>REFERENCES</h2>
-
-<ul>
-
-<li>Computing Visibility on Terrains in External Memory. Herman 
-Haverkort, Laura Toma and Yi Zhuang. In <em>ACM Journal on Experimental 
-Algorithmics (JEA)</em> 13 (2009).</li>
-	 
-<li><a href="http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.76.4282&rep=rep1&type=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>
-
-<h2>SEE ALSO</h2>
-
-<em>
-<a href="r.los.html">r.los</a>,
-<a href="r.mapcalc.html">r.mapcalc</a>
-</em>
-
-<h2>AUTHORS</h2>
-
-<p>Laura Toma (Bowdoin College): <tt>ltoma at bowdoin.edu</tt>
-<p> Yi Zhuang (Carnegie-Mellon University): <tt>yzhuang at andrew.cmu.edu</tt>
-<p>William Richard (Bowdoin College): <tt>willster3021 at gmail.com </tt>
-<p>Markus Metz
-
-<p>
-<i>Last changed: $Date$</i>

Modified: grass-addons/grass6/raster/r.viewshed/rbbst.cpp
===================================================================
--- grass-addons/grass6/raster/r.viewshed/rbbst.cpp	2014-09-25 08:12:03 UTC (rev 62080)
+++ grass-addons/grass6/raster/r.viewshed/rbbst.cpp	2014-09-25 10:44:26 UTC (rev 62081)
@@ -188,14 +188,12 @@
 
 
 /*a function used to compare two doubles */
+/* a < b : -1
+ * a > b : 1
+ * a == b : 0 */
 char compare_double(double a, double b)
 {
-    if (fabs(a - b) < EPSILON)
-	return 0;
-    if (a - b < 0)
-	return -1;
-
-    return 1;
+    return (a < b ? -1 : (a > b));
 }
 
 
@@ -264,7 +262,7 @@
     TreeNode *inserted = nextNode;
 
     /*update augmented maxGradient */
-    nextNode->value.maxGradient = nextNode->value.gradient[1];
+    nextNode->value.maxGradient = find_value_min_value(&nextNode->value);
     while (nextNode->parent != NIL) {
 	if (nextNode->parent->value.maxGradient < nextNode->value.maxGradient)
 	    nextNode->parent->value.maxGradient = nextNode->value.maxGradient;
@@ -432,7 +430,7 @@
     double left, right;
 
     while (curNode->parent != NIL) {
-	if (curNode->parent->value.maxGradient == y->value.gradient[1]) {
+	if (curNode->parent->value.maxGradient == find_value_min_value(&y->value)) {
 	    left = find_max_value(curNode->parent->left);
 	    right = find_max_value(curNode->parent->right);
 
@@ -441,10 +439,10 @@
 	    else
 		curNode->parent->value.maxGradient = right;
 
-	    if (curNode->parent->value.gradient[1] >
+	    if (find_value_min_value(&curNode->parent->value) >
 		curNode->parent->value.maxGradient)
 		curNode->parent->value.maxGradient =
-		    curNode->parent->value.gradient[1];
+		    find_value_min_value(&curNode->parent->value);
 	}
 	else {
 	    break;
@@ -458,14 +456,14 @@
 	toFix->left->value.maxGradient >
 	toFix->right->value.maxGradient ? toFix->left->value.
 	maxGradient : toFix->right->value.maxGradient;
-    if (tmpMax > toFix->value.gradient[1])
+    if (tmpMax > find_value_min_value(&toFix->value))
 	toFix->value.maxGradient = tmpMax;
     else
-	toFix->value.maxGradient = toFix->value.gradient[1];
+	toFix->value.maxGradient = find_value_min_value(&toFix->value);
 
     /*13-15 */
     if (y != NIL && y != z) {
-	double zGradient = z->value.gradient[1];
+	double zGradient = find_value_min_value(&z->value);
 
 	z->value.key = y->value.key;
 	z->value.gradient[0] = y->value.gradient[0];
@@ -482,14 +480,14 @@
 	    toFix->left->value.maxGradient >
 	    toFix->right->value.maxGradient ? toFix->left->value.
 	    maxGradient : toFix->right->value.maxGradient;
-	if (tmpMax > toFix->value.gradient[1])
+	if (tmpMax > find_value_min_value(&toFix->value))
 	    toFix->value.maxGradient = tmpMax;
 	else
-	    toFix->value.maxGradient = toFix->value.gradient[1];
+	    toFix->value.maxGradient = find_value_min_value(&toFix->value);
 
 	while (z->parent != NIL) {
 	    if (z->parent->value.maxGradient == zGradient) {
-		if (z->parent->value.gradient[1] != zGradient &&
+		if (find_value_min_value(&z->parent->value) != zGradient &&
 		    (!(z->parent->left->value.maxGradient == zGradient &&
 		       z->parent->right->value.maxGradient == zGradient))) {
 
@@ -501,10 +499,10 @@
 		    else
 			z->parent->value.maxGradient = right;
 
-		    if (z->parent->value.gradient[1] >
+		    if (find_value_min_value(&z->parent->value) >
 			z->parent->value.maxGradient)
 			z->parent->value.maxGradient =
-			    z->parent->value.gradient[1];
+			    find_value_min_value(&z->parent->value);
 
 		}
 
@@ -608,6 +606,24 @@
     return;
 }
 
+/*find the min value in the given tree value */
+double find_value_min_value(TreeValue * v)
+{
+    if (v->gradient[0] < v->gradient[1]) {
+	if (v->gradient[0] < v->gradient[2])
+	    return v->gradient[0];
+	else
+	    return v->gradient[2];
+    }
+    else {
+	if (v->gradient[1] < v->gradient[2])
+	    return v->gradient[1];
+	else
+	    return v->gradient[2];
+    }
+    return v->gradient[0];
+}
+
 /*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)
@@ -629,10 +645,11 @@
     TreeNode *keyNode = search_for_node(root, maxKey);
 
     if (keyNode == NIL) {
-	/*fprintf(stderr, "key node not found. error occured!\n");
+	/*fprintf(stderr, "key node not found. error occurred!\n");
 	   //there is no point in the structure with key < maxKey */
+	/*node not found */
+	G_fatal_error(_("Attempt to find node with key=%f failed"), maxKey);
 	return SMALLEST_GRADIENT;
-	exit(1);
     }
 
     TreeNode *currNode = keyNode;
@@ -640,17 +657,33 @@
     double tmpMax;
     double curr_gradient;
 
+    while (currNode->parent != NIL) {
+	if (currNode == currNode->parent->right) {	/*its the right node of its parent; */
+	    tmpMax = find_max_value(currNode->parent->left);
+	    if (tmpMax > max)
+		max = tmpMax;
+	    if (find_value_min_value(&currNode->parent->value) > max)
+		max = find_value_min_value(&currNode->parent->value);
+	}
+	currNode = currNode->parent;
+    }
+
+    if (max > gradient)
+	return max;
+
     /* traverse all nodes with smaller distance */
+    max = SMALLEST_GRADIENT;
+    currNode = keyNode;
     while (currNode != NIL) {
 	int checkme = (currNode->value.angle[0] <= angle &&
 	              currNode->value.angle[2] >= angle);
 		      
 	if (!checkme && currNode->value.key > 0) {
-	    G_warning(_("\nangles outside angle %.4f"), angle);
+	    G_warning(_("Angles outside angle %.4f"), angle);
 	    G_warning(_("ENTER angle %.4f"), currNode->value.angle[0]);
 	    G_warning(_("CENTER angle %.4f"), currNode->value.angle[1]);
 	    G_warning(_("EXIT angle %.4f"), currNode->value.angle[2]);
-	    G_warning(_("\nENTER gradient %.4f"), currNode->value.gradient[0]);
+	    G_warning(_("ENTER gradient %.4f"), currNode->value.gradient[0]);
 	    G_warning(_("CENTER gradient %.4f"), currNode->value.gradient[1]);
 	    G_warning(_("EXIT gradient %.4f"), currNode->value.gradient[2]);
 	}
@@ -729,20 +762,20 @@
     tmpMax = x->left->value.maxGradient > y->left->value.maxGradient ?
 	x->left->value.maxGradient : y->left->value.maxGradient;
 
-    if (tmpMax > x->value.gradient[1])
+    if (tmpMax > find_value_min_value(&x->value))
 	x->value.maxGradient = tmpMax;
     else
-	x->value.maxGradient = x->value.gradient[1];
+	x->value.maxGradient = find_value_min_value(&x->value);
 
 
     /*fix y */
     tmpMax = x->value.maxGradient > y->right->value.maxGradient ?
 	x->value.maxGradient : y->right->value.maxGradient;
 
-    if (tmpMax > y->value.gradient[1])
+    if (tmpMax > find_value_min_value(&y->value))
 	y->value.maxGradient = tmpMax;
     else
-	y->value.maxGradient = y->value.gradient[1];
+	y->value.maxGradient = find_value_min_value(&y->value);
 
     /*left rotation
        //see pseudocode on page 278 in CLRS */
@@ -781,19 +814,19 @@
     tmpMax = x->right->value.maxGradient > y->right->value.maxGradient ?
 	x->right->value.maxGradient : y->right->value.maxGradient;
 
-    if (tmpMax > y->value.gradient[1])
+    if (tmpMax > find_value_min_value(&y->value))
 	y->value.maxGradient = tmpMax;
     else
-	y->value.maxGradient = y->value.gradient[1];
+	y->value.maxGradient = find_value_min_value(&y->value);
 
     /*fix x */
     tmpMax = x->left->value.maxGradient > y->value.maxGradient ?
 	x->left->value.maxGradient : y->value.maxGradient;
 
-    if (tmpMax > x->value.gradient[1])
+    if (tmpMax > find_value_min_value(&x->value))
 	x->value.maxGradient = tmpMax;
     else
-	x->value.maxGradient = x->value.gradient[1];
+	x->value.maxGradient = find_value_min_value(&x->value);
 
     /*ratation */
     y->left = x->right;

Modified: grass-addons/grass6/raster/r.viewshed/rbbst.h
===================================================================
--- grass-addons/grass6/raster/r.viewshed/rbbst.h	2014-09-25 08:12:03 UTC (rev 62080)
+++ grass-addons/grass6/raster/r.viewshed/rbbst.h	2014-09-25 10:44:26 UTC (rev 62081)
@@ -147,6 +147,9 @@
 /*search for a node with the given key */
 TreeNode *search_for_node(TreeNode * root, double key);
 
+/*find the min value in the given tree value */
+double find_value_min_value(TreeValue * v);
+
 /*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);

Modified: grass-addons/grass6/raster/r.viewshed/statusstructure.cpp
===================================================================
--- grass-addons/grass6/raster/r.viewshed/statusstructure.cpp	2014-09-25 08:12:03 UTC (rev 62080)
+++ grass-addons/grass6/raster/r.viewshed/statusstructure.cpp	2014-09-25 10:44:26 UTC (rev 62081)
@@ -98,7 +98,7 @@
     G_verbose_message(_(" (key=%d, ptr=%d, total node=%d B)"),
 	   (int)sizeof(TreeValue),
 	   (int)sizeof(TreeNode *), (int)sizeof(TreeNode));
-    sizeBytes = sizeof(TreeNode) * max(hd->ncols, hd->nrows);
+    sizeBytes = sizeof(TreeNode) * std::max(hd->ncols, hd->nrows);
     G_verbose_message(_(" Total= %lld B"), sizeBytes);
     return sizeBytes;
 }

Modified: grass-addons/grass6/raster/r.viewshed/viewshed.cpp
===================================================================
--- grass-addons/grass6/raster/r.viewshed/viewshed.cpp	2014-09-25 08:12:03 UTC (rev 62080)
+++ grass-addons/grass6/raster/r.viewshed/viewshed.cpp	2014-09-25 10:44:26 UTC (rev 62081)
@@ -320,14 +320,14 @@
     long nvis = 0;		/*number of visible cells */
     AEvent *e;
 
-    G_message(_("Determine visibility..."));
+    G_important_message(_("Computing visibility..."));
     G_percent(0, 100, 2);
 
     for (size_t i = 0; i < nevents; i++) {
 
-	int perc = (int)((double)i / nevents * 1000000.);
+	int perc = (int)(1000000 * i / nevents);
 	if (perc > 0 && perc < 1000000)
-	    G_percent(perc, 1000000, 2);
+	    G_percent(perc, 1000000, 1);
 
 	/*get out one event at a time and process it according to its type */
 	e = &(eventList[i]);
@@ -506,7 +506,7 @@
        structure */
     StatusNode sn;
 
-    G_message(_("Calculating distances..."));
+    G_message(_("Initialize sweepline..."));
 
     rt_start(sweepTime);
     for (dimensionType i = vp->col + 1; i < hd->ncols; i++) {
@@ -576,9 +576,9 @@
 
     for (off_t i = 0; i < nbEvents; i++) {
 
-	int perc = (int)((double)i / nbEvents * 1000000.);
+	int perc = (int)(1000000 * i / nbEvents);
 	if (perc > 0)
-	    G_percent(perc, 1000000, 2);
+	    G_percent(perc, 1000000, 1);
 
 	/*get out one event at a time and process it according to its type */
 	ae = eventList->read_item(&e);



More information about the grass-commit mailing list