[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