[GRASS-SVN] r34825 - grass-addons/raster/r.viewshed
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Dec 11 03:10:58 EST 2008
Author: hamish
Date: 2008-12-11 03:10:58 -0500 (Thu, 11 Dec 2008)
New Revision: 34825
Modified:
grass-addons/raster/r.viewshed/distribute.cc
grass-addons/raster/r.viewshed/distribute.h
grass-addons/raster/r.viewshed/eventlist.cc
grass-addons/raster/r.viewshed/eventlist.h
grass-addons/raster/r.viewshed/grass.cc
grass-addons/raster/r.viewshed/grass.h
grass-addons/raster/r.viewshed/grid.cc
grass-addons/raster/r.viewshed/grid.h
grass-addons/raster/r.viewshed/main.cc
grass-addons/raster/r.viewshed/rbbst.h
grass-addons/raster/r.viewshed/statusstructure.cc
grass-addons/raster/r.viewshed/statusstructure.h
grass-addons/raster/r.viewshed/viewshed.cc
grass-addons/raster/r.viewshed/viewshed.h
grass-addons/raster/r.viewshed/visibility.cc
grass-addons/raster/r.viewshed/visibility.h
Log:
run grass_indent.sh
Modified: grass-addons/raster/r.viewshed/distribute.cc
===================================================================
--- grass-addons/raster/r.viewshed/distribute.cc 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/distribute.cc 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -52,7 +53,7 @@
/* include IOSTREAM header */
#include <grass/iostream/ami.h>
-#else
+#else
/* if GRASS is not defined */
#include <ami.h>
#endif
@@ -89,89 +90,95 @@
////////////////////////////////////////////////////////////////////////
/* distribution sweep: write results to visgrid.
*/
-IOVisibilityGrid *distribute_and_sweep(char* inputfname,
- GridHeader * hd,
- Viewpoint *vp,
- ViewOptions viewOptions) {
-
- assert(inputfname && hd && vp);
- print_message("Start distributed sweeping.\n");
- fflush(stdout);
-
- /* ------------------------------ */
- /*initialize the visibility grid */
- IOVisibilityGrid *visgrid;
- visgrid = init_io_visibilitygrid(*hd, *vp);
- printf("distribute_and_sweep: visgrid=%s\n", visgrid->visStr->name());
-
-
- /* ------------------------------ */
- /*construct event list corresp to the input file and viewpoint */
- Rtimer initEventTime;
- rt_start(initEventTime);
- AMI_STREAM < AEvent > *eventList;
+IOVisibilityGrid *distribute_and_sweep(char *inputfname,
+ GridHeader * hd,
+ Viewpoint * vp,
+ ViewOptions viewOptions)
+{
+
+ assert(inputfname && hd && vp);
+ print_message("Start distributed sweeping.\n");
+ fflush(stdout);
+
+ /* ------------------------------ */
+ /*initialize the visibility grid */
+ IOVisibilityGrid *visgrid;
+
+ visgrid = init_io_visibilitygrid(*hd, *vp);
+ printf("distribute_and_sweep: visgrid=%s\n", visgrid->visStr->name());
+
+
+ /* ------------------------------ */
+ /*construct event list corresp to the input file and viewpoint */
+ Rtimer initEventTime;
+
+ rt_start(initEventTime);
+ AMI_STREAM < AEvent > *eventList;
#ifdef __GRASS__
- eventList = grass_init_event_list(inputfname, vp, hd,
- viewOptions, NULL, visgrid);
+ eventList = grass_init_event_list(inputfname, vp, hd,
+ viewOptions, NULL, visgrid);
#else
- eventList = init_event_list(inputfname, vp, hd, viewOptions,NULL, visgrid);
+ eventList =
+ init_event_list(inputfname, vp, hd, viewOptions, NULL, visgrid);
#endif
- assert(eventList);
- eventList->seek(0);
- rt_stop(initEventTime);
- printf("distribute_and_sweep: eventlist=%s\n", eventList->sprint());
-
-
- /* ------------------------------ */
- /*sort the events concentrically */
- Rtimer sortEventTime;
- rt_start(sortEventTime);
- PRINT_DISTRIBUTE {
+ assert(eventList);
+ eventList->seek(0);
+ rt_stop(initEventTime);
+ printf("distribute_and_sweep: eventlist=%s\n", eventList->sprint());
+
+
+ /* ------------------------------ */
+ /*sort the events concentrically */
+ Rtimer sortEventTime;
+
+ rt_start(sortEventTime);
+ PRINT_DISTRIBUTE {
print_message("sorting events by distance from viewpoint..");
fflush(stdout);
- }
+ }
- sort_event_list_by_distance(&eventList, *vp);
- PRINT_DISTRIBUTE {
+ sort_event_list_by_distance(&eventList, *vp);
+ PRINT_DISTRIBUTE {
print_message("..sorting done.\n");
fflush(stdout);
- }
+ }
- /* debugging */
- /*sortCheck(eventList, *vp); */
- eventList->seek(0); /*this does not seem to be ensured by sort */
- rt_stop(sortEventTime);
- printf("distribute_and_sweep: eventlist=%s\n", eventList->sprint());
-
-
-
- /* ------------------------------ */
- /* start distribution */
- long nvis; /*number of cells visible. Returned by distribute_sector */
-
- Rtimer sweepTime;
- rt_start(sweepTime);
- /*distribute recursively the events and write results to visgrid.
- invariant: distribute_sector deletes its eventlist */
- nvis = distribute_sector(eventList, NULL, 0, ANGLE_FACTOR * 2 * M_PI,
- visgrid, vp, viewOptions);
- rt_stop(sweepTime);
-
-
- /* ------------------------------ */
- /*cleanup */
- print_message("Distribution sweeping done.\n");
- fflush(stdout);
-
- printf("Total cells %ld, visible cells %ld (%.1f percent).\n",
- (long)visgrid->hd->nrows * visgrid->hd->ncols,
- nvis,
- (float)((float)nvis * 100 /
- (float)(visgrid->hd->nrows * visgrid->hd->ncols)));
-
- print_viewshed_timings(initEventTime, sortEventTime, sweepTime);
-
- return visgrid;
+ /* debugging */
+ /*sortCheck(eventList, *vp); */
+ eventList->seek(0); /*this does not seem to be ensured by sort */
+ rt_stop(sortEventTime);
+ printf("distribute_and_sweep: eventlist=%s\n", eventList->sprint());
+
+
+
+ /* ------------------------------ */
+ /* start distribution */
+ long nvis; /*number of cells visible. Returned by distribute_sector */
+
+ Rtimer sweepTime;
+
+ rt_start(sweepTime);
+ /*distribute recursively the events and write results to visgrid.
+ invariant: distribute_sector deletes its eventlist */
+ nvis = distribute_sector(eventList, NULL, 0, ANGLE_FACTOR * 2 * M_PI,
+ visgrid, vp, viewOptions);
+ rt_stop(sweepTime);
+
+
+ /* ------------------------------ */
+ /*cleanup */
+ print_message("Distribution sweeping done.\n");
+ fflush(stdout);
+
+ printf("Total cells %ld, visible cells %ld (%.1f percent).\n",
+ (long)visgrid->hd->nrows * visgrid->hd->ncols,
+ nvis,
+ (float)((float)nvis * 100 /
+ (float)(visgrid->hd->nrows * visgrid->hd->ncols)));
+
+ print_viewshed_timings(initEventTime, sortEventTime, sweepTime);
+
+ return visgrid;
}
@@ -192,12 +199,13 @@
invariant: distribute_sector deletes eventList and enterBndEvents
*/
unsigned long distribute_sector(AMI_STREAM < AEvent > *eventList,
- AMI_STREAM < AEvent > *enterBndEvents,
- double start_angle, double end_angle,
- IOVisibilityGrid * visgrid, Viewpoint * vp,
- ViewOptions viewOptions) {
-
+ AMI_STREAM < AEvent > *enterBndEvents,
+ double start_angle, double end_angle,
+ IOVisibilityGrid * visgrid, Viewpoint * vp,
+ ViewOptions viewOptions)
+{
+
assert(eventList && visgrid && vp);
/*enterBndEvents may be NULL first time */
@@ -206,64 +214,66 @@
printf("initial_gradient: %lf\n", SMALLEST_GRADIENT);
printf("eventlist: %s\n", eventList->sprint());
if (enterBndEvents)
- printf("BndEvents: %s\n", enterBndEvents->sprint());
+ printf("BndEvents: %s\n", enterBndEvents->sprint());
PRINT_DISTRIBUTE LOG_avail_memo();
-
+
unsigned long nvis = 0;
-
+
/*******************************************************
BASE CASE
*******************************************************/
if (eventList->stream_len() * sizeof(AEvent) <
- MM_manager.memory_available()) {
- if (enterBndEvents) {
+ MM_manager.memory_available()) {
+ if (enterBndEvents) {
nvis += solve_in_memory(eventList, enterBndEvents,
- start_angle,end_angle,visgrid,vp,viewOptions);
+ start_angle, end_angle, visgrid, vp,
+ viewOptions);
return nvis;
- }
- else {
+ }
+ else {
/*we are here if the problem fits in memory, and
- //enterBNdEvents==NULL---this only happens the very first time;
- //technically at this point we should do one pass though the
- //data and collect the events that cross the 1st and 4th
- //quadrant; instead we force recursion, nect= 2 */
- }
+ //enterBNdEvents==NULL---this only happens the very first time;
+ //technically at this point we should do one pass though the
+ //data and collect the events that cross the 1st and 4th
+ //quadrant; instead we force recursion, nect= 2 */
+ }
}
-
+
/*else, must recurse */
PRINT_DISTRIBUTE print_message("in EXTERNAL memory\n");
-
+
/*compute number of sectors */
int nsect = compute_n_sectors();
- assert(nsect > 1);
-
+
+ assert(nsect > 1);
+
/*an array of streams, one for each sector; sector[i] will keep all
- the cells that are completely inside sector i */
+ the cells that are completely inside sector i */
AMI_STREAM < AEvent > *sector = new AMI_STREAM < AEvent >[nsect];
-
+
/*the array of gradient values, one for each sector; the gradient is
- the gradient of the center of a cell that spans the sector
- completely */
+ the gradient of the center of a cell that spans the sector
+ completely */
double *high = new double[nsect];
-
+
for (int i = 0; i < nsect; i++)
- high[i] = SMALLEST_GRADIENT;
-
+ high[i] = SMALLEST_GRADIENT;
+
/*an array of streams, one for each stream boundary; sectorBnd[0]
- will keep all the cells crossing into sector 0 from below; and so on. */
+ will keep all the cells crossing into sector 0 from below; and so on. */
AMI_STREAM < AEvent > *sectorBnd = new AMI_STREAM < AEvent >[nsect];
-
+
/*keeps stats for each sector */
long *total = new long[nsect];
long *insert = new long[nsect];
long *drop = new long[nsect];
long *bndInsert = new long[nsect];
long *bndDrop = new long[nsect];
-
+
for (int i = 0; i < nsect; i++)
- total[i] = insert[i] = drop[i] = bndInsert[i] = bndDrop[i] = 0;
+ total[i] = insert[i] = drop[i] = bndInsert[i] = bndDrop[i] = 0;
long longEvents = 0;
-
+
/*******************************************************
CONCENTRIC SWEEP
*******************************************************/
@@ -273,142 +283,142 @@
int exit_s, s, enter_s;
long boundaryEvents = 0;
off_t nbEvents = eventList->stream_len();
-
+
eventList->seek(0);
for (off_t i = 0; i < nbEvents; i++) {
-
- /*get out one event at a time and process it according to its type */
- ae = eventList->read_item(&e);
- assert(ae == AMI_ERROR_NO_ERROR);
- assert(is_inside(e, start_angle, end_angle));
-
- /*compute its sector */
- s = get_event_sector(e->angle, start_angle, end_angle, nsect);
- /*detect boundary cases ==> precision issues */
- if (is_almost_on_boundary(e->angle, s, start_angle, end_angle, nsect)) {
+
+ /*get out one event at a time and process it according to its type */
+ ae = eventList->read_item(&e);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ assert(is_inside(e, start_angle, end_angle));
+
+ /*compute its sector */
+ s = get_event_sector(e->angle, start_angle, end_angle, nsect);
+ /*detect boundary cases ==> precision issues */
+ if (is_almost_on_boundary(e->angle, s, start_angle, end_angle, nsect)) {
double ssize = (end_angle - start_angle) / nsect;
-
+
boundaryEvents++;
PRINTWARNING {
- print_message("WARNING!event ");
- print_event(*e);
- print_message("CLOSE TO BOUNDARY\n");
- printf("angle=%f close to sector boundaries=[%f, %f]\n",
- e->angle, s * ssize, (s + 1) * ssize);
+ print_message("WARNING!event ");
+ print_event(*e);
+ print_message("CLOSE TO BOUNDARY\n");
+ printf("angle=%f close to sector boundaries=[%f, %f]\n",
+ e->angle, s * ssize, (s + 1) * ssize);
}
- }
-
- DISTRIBDEBUG printf("event %7lu: ", (unsigned long)i);
- DISTRIBDEBUG print_event(*e);
- DISTRIBDEBUG printf("d=%8.1f, ",
- get_square_distance_from_viewpoint(*e, *vp));
- DISTRIBDEBUG printf("s=%3d ", s);
-
- assert(is_inside(s, nsect));
- total[s]++;
-
- /*insert event in sector if not occluded */
- insert_event_in_sector(e, s, §or[s], high[s], vp, insert, drop);
-
- switch (e->eventType) {
- case CENTER_EVENT:
+ }
+
+ DISTRIBDEBUG printf("event %7lu: ", (unsigned long)i);
+ DISTRIBDEBUG print_event(*e);
+ DISTRIBDEBUG printf("d=%8.1f, ",
+ get_square_distance_from_viewpoint(*e, *vp));
+ DISTRIBDEBUG printf("s=%3d ", s);
+
+ assert(is_inside(s, nsect));
+ total[s]++;
+
+ /*insert event in sector if not occluded */
+ insert_event_in_sector(e, s, §or[s], high[s], vp, insert, drop);
+
+ switch (e->eventType) {
+ case CENTER_EVENT:
break;
-
- case ENTERING_EVENT:
+
+ case ENTERING_EVENT:
/*find its corresponding exit event and its sector */
exit_angle = calculate_exit_angle(e->row, e->col, vp);
exit_s =
- get_event_sector(exit_angle, start_angle, end_angle, nsect);
+ get_event_sector(exit_angle, start_angle, end_angle, nsect);
DISTRIBDEBUG
- printf(" ENTER (a=%.2f,s=%3d)---> EXIT (a=%.2f,s=%3d) ",
- e->angle, s, exit_angle, exit_s);
+ printf(" ENTER (a=%.2f,s=%3d)---> EXIT (a=%.2f,s=%3d) ",
+ e->angle, s, exit_angle, exit_s);
/*note: exit_s can be -1 (outside) */
if (exit_s == s) {
- /*short event, fit in sector s */
-
+ /*short event, fit in sector s */
+
}
else if (exit_s == (s + 1) % nsect || (exit_s + 1) % nsect == s) {
- /*semi-short event; insert in sector s, and in sector boundary s+1
- NOTE: to avoid precision issues, the events are inserted
- when processing the EXIT_EVENT
- insertEventInSector(e, (s+1)%nsect, §orBnd[(s+1)%nsect],
- high[(s+1)%nsect], vp,bndInsert, bndDrop); */
-
+ /*semi-short event; insert in sector s, and in sector boundary s+1
+ NOTE: to avoid precision issues, the events are inserted
+ when processing the EXIT_EVENT
+ insertEventInSector(e, (s+1)%nsect, §orBnd[(s+1)%nsect],
+ high[(s+1)%nsect], vp,bndInsert, bndDrop); */
+
}
else {
- /*long event; insert in sector s, and in sector boundary exit_s */
- process_long_cell(s, exit_s, nsect, vp, e, high);
- longEvents++;
- /*insertEventInSector(e, exit_s, §orBnd[exit_s], high[exit_s],
- // vp, bndInsert, bndDrop); */
+ /*long event; insert in sector s, and in sector boundary exit_s */
+ process_long_cell(s, exit_s, nsect, vp, e, high);
+ longEvents++;
+ /*insertEventInSector(e, exit_s, §orBnd[exit_s], high[exit_s],
+ // vp, bndInsert, bndDrop); */
}
break;
- case EXITING_EVENT:
+ case EXITING_EVENT:
/*find its corresponding enter event and its sector */
enter_angle = calculate_enter_angle(e->row, e->col, vp);
enter_s =
- get_event_sector(enter_angle, start_angle, end_angle, nsect);
+ get_event_sector(enter_angle, start_angle, end_angle, nsect);
DISTRIBDEBUG
- printf(" EXIT (a=%.2f,s=%3d)--->ENTER (a=%.2f,s=%3d) ",
- e->angle, s, enter_angle, enter_s);
-
+ printf(" EXIT (a=%.2f,s=%3d)--->ENTER (a=%.2f,s=%3d) ",
+ e->angle, s, enter_angle, enter_s);
+
/*don't need to check spanned sectors because it is done on its
- //ENTER event; actually...you do, because its enter event may
- //not be in this sector==> enter_s = -1 (outside) */
+ //ENTER event; actually...you do, because its enter event may
+ //not be in this sector==> enter_s = -1 (outside) */
if (enter_s == s) {
- /*short event, fit in sector */
+ /*short event, fit in sector */
}
else if (enter_s == (s + 1) % nsect || (enter_s + 1) % nsect == s) {
- /*semi-short event
- //the corresponding ENTER event must insert itself in sectorBnd[s] */
- e->eventType = ENTERING_EVENT;
- BND_DEBUG {
+ /*semi-short event
+ //the corresponding ENTER event must insert itself in sectorBnd[s] */
+ e->eventType = ENTERING_EVENT;
+ BND_DEBUG {
print_message("BND event ");
print_event(*e);
printf("in bndSector %d\n", s);
fflush(stdout);
- }
- insert_event_in_sector(e, s, §orBnd[s], high[s],
- vp, bndInsert, bndDrop);
-
+ }
+ insert_event_in_sector(e, s, §orBnd[s], high[s],
+ vp, bndInsert, bndDrop);
+
}
else {
- /*long event */
- process_long_cell(enter_s, s, nsect, vp, e, high);
- longEvents++;
- /*the corresponding ENTER event must insert itself in sectorBnd[s] */
- e->eventType = ENTERING_EVENT;
- BND_DEBUG {
+ /*long event */
+ process_long_cell(enter_s, s, nsect, vp, e, high);
+ longEvents++;
+ /*the corresponding ENTER event must insert itself in sectorBnd[s] */
+ e->eventType = ENTERING_EVENT;
+ BND_DEBUG {
print_message("BND event ");
print_event(*e);
printf("in bndSector %d\n", s);
fflush(stdout);
- }
- insert_event_in_sector(e, s, §orBnd[s], high[s],
- vp, bndInsert, bndDrop);
+ }
+ insert_event_in_sector(e, s, §orBnd[s], high[s],
+ vp, bndInsert, bndDrop);
}
break;
-
- } /*switch event-type */
-
- DISTRIBDEBUG print_message("\n");
+
+ } /*switch event-type */
+
+ DISTRIBDEBUG print_message("\n");
} /*for event i */
-
-
+
+
/*distribute the enterBnd events to the boundary streams of the
- //sectors; note: the boundary streams are not sorted by distance. */
+ //sectors; note: the boundary streams are not sorted by distance. */
if (enterBndEvents)
- distribute_bnd_events(enterBndEvents, sectorBnd, nsect, vp,
- start_angle, end_angle, high, bndInsert,
- bndDrop);
-
+ distribute_bnd_events(enterBndEvents, sectorBnd, nsect, vp,
+ start_angle, end_angle, high, bndInsert,
+ bndDrop);
+
/*sanity checks */
PRINT_DISTRIBUTE printf("boundary events in distribution: %ld\n",
- boundaryEvents);
+ boundaryEvents);
print_sector_stats(nbEvents, nsect, high, total, insert, drop, sector,
- sectorBnd, bndInsert,
- longEvents, start_angle, end_angle);
+ sectorBnd, bndInsert,
+ longEvents, start_angle, end_angle);
/*cleanup after sector stats */
delete[]total;
delete[]insert;
@@ -416,13 +426,13 @@
delete[]high;
delete[]bndInsert;
delete[]bndDrop;
-
+
/*we finished processing this sector, delete the event list */
delete eventList;
if (enterBndEvents)
- delete enterBndEvents;
-
+ delete enterBndEvents;
+
/*save stream names of new sectors */
#ifdef __GRASS__
char **sectorName = (char **)G_malloc(nsect * sizeof(char *));
@@ -435,47 +445,50 @@
for (int i = 0; i < nsect; i++) {
sector[i].name(§orName[i]);
PRINT_DISTRIBUTE printf("saving stream %d: %s\t", i, sectorName[i]);
-
+
sector[i].persist(PERSIST_PERSISTENT);
sectorBnd[i].name(§orBndName[i]);
PRINT_DISTRIBUTE printf("saving BndStr %d: %s\n", i,
- sectorBndName[i]);
+ sectorBndName[i]);
sectorBnd[i].persist(PERSIST_PERSISTENT);
}
-
+
/*delete [] sector;
- //does this call delete on every single stream?
- //for (int i=0; i< nsect; i++ ) delete sector[i]; */
+ //does this call delete on every single stream?
+ //for (int i=0; i< nsect; i++ ) delete sector[i]; */
delete[]sector;
delete[]sectorBnd;
/*LOG_avail_memo(); */
-
+
/*solve recursively each sector */
for (int i = 0; i < nsect; i++) {
-
- /*recover stream */
- PRINT_DISTRIBUTE printf("\nopening sector stream %s ", sectorName[i]);
-
- AMI_STREAM < AEvent > *str =
+
+ /*recover stream */
+ PRINT_DISTRIBUTE printf("\nopening sector stream %s ", sectorName[i]);
+
+ AMI_STREAM < AEvent > *str =
new AMI_STREAM < AEvent > (sectorName[i]);
- assert(str);
- PRINT_DISTRIBUTE printf(" len=%lu\n",
- (unsigned long)str->stream_len());
- /*recover boundary stream */
- PRINT_DISTRIBUTE printf("opening boundary sector stream %s ",
- sectorBndName[i]);
- AMI_STREAM < AEvent > *bndStr =
+ assert(str);
+ PRINT_DISTRIBUTE printf(" len=%lu\n",
+ (unsigned long)str->stream_len());
+ /*recover boundary stream */
+ PRINT_DISTRIBUTE printf("opening boundary sector stream %s ",
+ sectorBndName[i]);
+ AMI_STREAM < AEvent > *bndStr =
new AMI_STREAM < AEvent > (sectorBndName[i]);
- assert(str);
- PRINT_DISTRIBUTE printf(" len=%lu\n",
- (unsigned long)bndStr->stream_len());
-
-
- nvis += distribute_sector(str, bndStr,
- start_angle+i*((end_angle-start_angle)/nsect),
- start_angle+(i+1)*((end_angle-start_angle)/nsect),
- visgrid, vp, viewOptions);
+ assert(str);
+ PRINT_DISTRIBUTE printf(" len=%lu\n",
+ (unsigned long)bndStr->stream_len());
+
+
+ nvis += distribute_sector(str, bndStr,
+ start_angle +
+ i * ((end_angle - start_angle) / nsect),
+ start_angle + (i +
+ 1) * ((end_angle -
+ start_angle) / nsect),
+ visgrid, vp, viewOptions);
}
/*cleanup */
@@ -486,10 +499,10 @@
free(sectorName);
free(sectorBndName);
#endif
-
+
PRINT_DISTRIBUTE printf("Distribute sector [ %.4f, %.4f] done.\n",
- start_angle, end_angle);
-
+ start_angle, end_angle);
+
return nvis;
}
@@ -512,45 +525,45 @@
{
PRINT_DISTRIBUTE printf("Distribute boundary of sector [ %.4f, %.4f] ",
- start_angle, end_angle);
+ start_angle, end_angle);
assert(bndEvents && sectorBnd && vp && high && insert && drop);
AEvent *e;
AMI_err ae;
double exit_angle;
int exit_s;
off_t nbEvents = bndEvents->stream_len();
-
+
bndEvents->seek(0);
for (off_t i = 0; i < nbEvents; i++) {
-
- /*get out one event at a time */
- ae = bndEvents->read_item(&e);
- assert(ae == AMI_ERROR_NO_ERROR);
- /*must be outside, but better not check to avoid precision issues
- //assert(!is_inside(e, start_angle, end_angle));
-
- //each event must be an ENTER event that falls in a diff sector
- //than its EXIT */
- assert(e->eventType == ENTERING_EVENT);
-
- /*find its corresponding exit event and its sector */
+
+ /*get out one event at a time */
+ ae = bndEvents->read_item(&e);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ /*must be outside, but better not check to avoid precision issues
+ //assert(!is_inside(e, start_angle, end_angle));
+
+ //each event must be an ENTER event that falls in a diff sector
+ //than its EXIT */
+ assert(e->eventType == ENTERING_EVENT);
+
+ /*find its corresponding exit event and its sector */
exit_angle = calculate_exit_angle(e->row, e->col, vp);
exit_s = get_event_sector(exit_angle, start_angle, end_angle, nsect);
-
+
/*exit_s cannot be outside sector; though we have to be careful
- //with precision */
+ //with precision */
assert(is_inside(exit_s, nsect));
-
+
/*insert this event in the boundary stream of this sub-sector */
insert_event_in_sector(e, exit_s, §orBnd[exit_s], high[exit_s],
- vp, insert, drop);
-
+ vp, insert, drop);
+
} /*for i */
-
+
PRINT_DISTRIBUTE
- printf("Distribute boundary of sector [ %.4f, %.4f] done.\n",
- start_angle, end_angle);
-
+ printf("Distribute boundary of sector [ %.4f, %.4f] done.\n",
+ start_angle, end_angle);
+
return;
}
@@ -580,27 +593,30 @@
printf("solve_in_memory: eventlist: %s\n", eventList->sprint());
if (enterBndEvents)
- printf("BndEvents: %s\n", enterBndEvents->sprint());
-
+ printf("BndEvents: %s\n", enterBndEvents->sprint());
+
if (eventList->stream_len() == 0) {
- delete eventList;
- if (enterBndEvents) delete enterBndEvents;
- return nvis;
+ delete eventList;
+
+ if (enterBndEvents)
+ delete enterBndEvents;
+
+ return nvis;
}
/*sort the events radially */
sort_event_list(&eventList);
-
+
/*create the status structure */
StatusList *status_struct = create_status_struct();
-
+
/*initialize status structure with all ENTER events whose EXIT
//events is inside the sector */
AEvent *e;
AMI_err ae;
StatusNode sn;
int inevents = 0;
-
+
/*
eventList->seek(0);
double enter_angle;
@@ -630,116 +646,118 @@
assert(ae == AMI_ERROR_END_OF_STREAM);
*/
if (enterBndEvents) {
- enterBndEvents->seek(0);
- inevents = enterBndEvents->stream_len();
- for (off_t i = 0; i < inevents; i++) {
+ enterBndEvents->seek(0);
+ inevents = enterBndEvents->stream_len();
+ for (off_t i = 0; i < inevents; i++) {
ae = enterBndEvents->read_item(&e);
assert(ae == AMI_ERROR_NO_ERROR);
DEBUGINIT {
- print_message("INMEM init: initializing boundary ");
- print_event(*e);
- print_message("\n");
+ print_message("INMEM init: initializing boundary ");
+ print_event(*e);
+ print_message("\n");
}
/*this must span the first boundary of this sector; insert it
- //in status structure */
+ //in status structure */
sn.col = e->col;
sn.row = e->row;
sn.elev = e->elev;
calculate_dist_n_gradient(&sn, vp);
insert_into_status_struct(sn, status_struct);
- }
+ }
}
PRINT_DISTRIBUTE {
- printf("initialized active structure with %d events\n", inevents);
- fflush(stdout);
+ printf("initialized active structure with %d events\n", inevents);
+ fflush(stdout);
}
-
-
+
+
/*sweep the event list */
VisCell viscell;
off_t nbEvents = eventList->stream_len();
-
+
/*printf("nbEvents = %ld\n", (long) nbEvents); */
eventList->seek(0);
for (off_t i = 0; i < nbEvents; i++) {
-
- /*get out one event at a time and process it according to its type */
- ae = eventList->read_item(&e);
- assert(ae == AMI_ERROR_NO_ERROR);
- SOLVEINMEMDEBUG {
+
+ /*get out one event at a time and process it according to its type */
+ ae = eventList->read_item(&e);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ SOLVEINMEMDEBUG {
print_message("INMEM sweep: next event: ");
print_event(*e);
- }
-
- sn.col = e->col;
- sn.row = e->row;
- sn.elev = e->elev;
- /*calculate Distance to VP and Gradient */
- calculate_dist_n_gradient(&sn, vp);
-
- switch (e->eventType) {
- case ENTERING_EVENT:
+ }
+
+ sn.col = e->col;
+ sn.row = e->row;
+ sn.elev = e->elev;
+ /*calculate Distance to VP and Gradient */
+ calculate_dist_n_gradient(&sn, vp);
+
+ switch (e->eventType) {
+ case ENTERING_EVENT:
/*insert node into structure */
SOLVEINMEMDEBUG {
- print_message("..ENTER-EVENT: insert\n");
+ print_message("..ENTER-EVENT: insert\n");
}
/*don't insert if its close to the boundary---the segment was
- //already inserted in initialization above
- //if (!is_almost_on_boundary(e->angle, start_angle))
- //insertIntoStatusStruct(sn,status_struct); */
+ //already inserted in initialization above
+ //if (!is_almost_on_boundary(e->angle, start_angle))
+ //insertIntoStatusStruct(sn,status_struct); */
insert_into_status_struct(sn, status_struct);
break;
-
- case EXITING_EVENT:
+
+ case EXITING_EVENT:
/*delete node out of status structure */
SOLVEINMEMDEBUG {
- print_message("..EXIT-EVENT: delete\n");
- /*find its corresponding enter event and its sector */
- double enter_angle =
+ print_message("..EXIT-EVENT: delete\n");
+ /*find its corresponding enter event and its sector */
+ double enter_angle =
calculate_enter_angle(e->row, e->col, vp);
- printf(" EXIT (a=%f)--->ENTER (a=%f) ", e->angle,
- enter_angle);
+ printf(" EXIT (a=%f)--->ENTER (a=%f) ", e->angle,
+ enter_angle);
}
delete_from_status_struct(status_struct, sn.dist2vp);
break;
-
- case CENTER_EVENT:
+
+ case CENTER_EVENT:
SOLVEINMEMDEBUG {
- print_message("..QUERY-EVENT: query\n");
+ print_message("..QUERY-EVENT: query\n");
}
/*calculate visibility
-
- //note: if there is nothing in the status structure, it means
- //there is no prior event to occlude it, so this cell is
- //VISIBLE. this is taken care of in the status structure --- if
- //a query event comes when the structure is empty, the query
- //returns minimum gradient */
+
+ //note: if there is nothing in the status structure, it means
+ //there is no prior event to occlude it, so this cell is
+ //VISIBLE. this is taken care of in the status structure --- if
+ //a query event comes when the structure is empty, the query
+ //returns minimum gradient */
double max;
-
- max = find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
-
+
+ max =
+ find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
+
viscell.row = sn.row;
viscell.col = sn.col;
-
- if (max <= sn.gradient) {
- /*the point is visible */
- viscell.angle = get_vertical_angle(*vp, sn, viewOptions.doCurv);
- assert(viscell.angle >0);
- /* viscell.vis = VISIBLE; */
- add_result_to_io_visibilitygrid(visgrid, &viscell);
- /*make sure nvis is correct*/
- nvis++;
+
+ if (max <= sn.gradient) {
+ /*the point is visible */
+ viscell.angle =
+ get_vertical_angle(*vp, sn, viewOptions.doCurv);
+ assert(viscell.angle > 0);
+ /* viscell.vis = VISIBLE; */
+ add_result_to_io_visibilitygrid(visgrid, &viscell);
+ /*make sure nvis is correct */
+ nvis++;
}
- else {
+ else {
/* else the point is invisible; we do not write it to the
- //visibility stream, because we only record visible and nodata
- //values to the stream */
- /* viscell.vis = INVISIBLE; */
- /* add_result_to_io_visibilitygrid(visgrid, &viscell); */
- }
+ //visibility stream, because we only record visible and nodata
+ //values to the stream */
+ /* viscell.vis = INVISIBLE; */
+ /* add_result_to_io_visibilitygrid(visgrid, &viscell); */
+ }
break;
- }
- } /* for each event */
+ }
+ } /* for each event */
PRINT_DISTRIBUTE print_message("in memory sweeping done.\n");
@@ -928,55 +946,55 @@
assert(insert[i] == sector[i].stream_len());
assert(bndInsert[i] == bndSector[i].stream_len());
-
+
totalSector += total[i];
totalDrop += drop[i];
totalInsert += insert[i];
-
+
}
assert(totalSector == nevents);
-
+
#ifdef __GRASS__
PRINT_DISTRIBUTE {
- G_message("-----nsectors=%d\n", nsect);
- for (int i = 0; i < nsect; i++) {
+ G_message("-----nsectors=%d\n", nsect);
+ for (int i = 0; i < nsect; i++) {
G_message("\ts=%3d ", i);
G_message("[%.4f, %.4f] ",
- get_sector_start(i, start_angle, end_angle, nsect),
- get_sector_end(i, start_angle, end_angle, nsect));
+ get_sector_start(i, start_angle, end_angle, nsect),
+ get_sector_end(i, start_angle, end_angle, nsect));
G_message("high = %9.1f, ", high[i]);
G_message("total = %10ld, ", total[i]);
G_message("inserted = %10ld, ", insert[i]);
G_message("dropped = %10ld, ", drop[i]);
G_message("BOUNDARY = %5ld", bndInsert[i]);
G_message("\n");
- }
+ }
}
G_message("Distribute [%.4f, %.4f]: nsect=%d, ", start_angle, end_angle,
- nsect);
+ nsect);
G_message
- ("total events %lu, inserted %lu, dropped %lu, long events=%ld\n",
- totalSector, totalInsert, totalDrop, longEvents);
+ ("total events %lu, inserted %lu, dropped %lu, long events=%ld\n",
+ totalSector, totalInsert, totalDrop, longEvents);
#else
PRINT_DISTRIBUTE {
- printf("-----nsectors=%d\n", nsect);
- for (int i = 0; i < nsect; i++) {
+ printf("-----nsectors=%d\n", nsect);
+ for (int i = 0; i < nsect; i++) {
printf("\ts=%3d ", i);
printf("[%.4f, %.4f] ",
- get_sector_start(i, start_angle, end_angle, nsect),
- get_sector_end(i, start_angle, end_angle, nsect));
+ get_sector_start(i, start_angle, end_angle, nsect),
+ get_sector_end(i, start_angle, end_angle, nsect));
printf("high = %9.1f, ", high[i]);
printf("total = %10ld, ", total[i]);
printf("inserted = %10ld, ", insert[i]);
printf("dropped = %10ld, ", drop[i]);
printf("BOUNDARY = %5ld", bndInsert[i]);
printf("\n");
- }
+ }
}
printf("Distribute [%.4f, %.4f]: nsect=%d, ", start_angle, end_angle,
- nsect);
+ nsect);
printf("total events %lu, inserted %lu, dropped %lu, long events=%ld\n",
- totalSector, totalInsert, totalDrop, longEvents);
+ totalSector, totalInsert, totalDrop, longEvents);
fflush(stdout);
#endif
return;
@@ -992,7 +1010,7 @@
long memSizeBytes = MM_manager.memory_available();
unsigned int blockSizeBytes = UntypedStream::get_block_length();
-
+
/*printf("computeNSect: block=%d, mem=%d\n", blockSizeBytes, (int)memSizeBytes); */
int nsect = (int)(memSizeBytes / (2 * blockSizeBytes));
@@ -1011,8 +1029,8 @@
if (nsect == 0 || nsect == 1)
nsect = 2;
else {
- /*we'll have 2 streams for each sector open; subtract 10 to be safe */
- if (2 * nsect > MAX_STREAMS_OPEN - 10)
+ /*we'll have 2 streams for each sector open; subtract 10 to be safe */
+ if (2 * nsect > MAX_STREAMS_OPEN - 10)
nsect = (MAX_STREAMS_OPEN - 10) / 2;
}
printf("nsectors set to %d\n", nsect);
Modified: grass-addons/raster/r.viewshed/distribute.h
===================================================================
--- grass-addons/raster/r.viewshed/distribute.h 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/distribute.h 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -48,10 +49,10 @@
/* distribution sweep: write results to visgrid.
*/
-IOVisibilityGrid *distribute_and_sweep(char* inputfname,
- GridHeader * hd,
- Viewpoint * vp,
- ViewOptions viewOptions);
+IOVisibilityGrid *distribute_and_sweep(char *inputfname,
+ GridHeader * hd,
+ Viewpoint * vp,
+ ViewOptions viewOptions);
@@ -70,10 +71,10 @@
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);
+ 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
@@ -98,10 +99,10 @@
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);
+ 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 */
@@ -109,7 +110,7 @@
/* 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);
+ double end_angle, int nsect);
/* computes the number of sector for the distribution sweep;
technically M/B */
@@ -121,7 +122,7 @@
/* 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);
+ int nsect);
/* insert event in this sector */
@@ -129,13 +130,13 @@
/* 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);
+ 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);
+ AMI_STREAM < AEvent > *str, long *insert);
@@ -150,24 +151,24 @@
/* 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);
+ 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);
+ 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);
+ int nsect);
/* return the start angle of the i-th sector. Assuming that
Modified: grass-addons/raster/r.viewshed/eventlist.cc
===================================================================
--- grass-addons/raster/r.viewshed/eventlist.cc 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/eventlist.cc 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -403,38 +404,43 @@
/* ------------------------------------------------------------ */
/*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. */
+ 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)
+ dimensionType row, dimensionType col,
+ float maxDist)
{
- /* it is not too smart to compare floats */
- if ((int)maxDist == INFINITY_DISTANCE)
+ /* 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) {
+ 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 {
+ }
+ else {
return 0;
- }
+ }
}
-void testeventlist(AEvent* elist, size_t n) {
+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);
+ 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);
}
/*///////////////////////////////////////////////////////////
@@ -447,104 +453,109 @@
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);
+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++) {
+
+ /* 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);
- }
- }
+
+ 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;
+ /*don't insert the viewpoint events in the list */
+ if (i == vp->row && j == vp->col)
+ 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;
+ /*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;
}
@@ -563,33 +574,34 @@
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)
+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);
- }
+ 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);
+ 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 = */
+ //the data
+ //GridHeader* foo = */
read_header_from_arcascii_file(grid_fp);
-
+
/*scan throught the arcascii file */
dimensionType i, j;
double ax, ay;
@@ -598,45 +610,46 @@
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) {
+ 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;
+ (*data)[j] = e.elev;
}
- }
+ }
- if (i == vp-> row && j == vp->col) {
+ 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);
+ 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)
+ }
+
+ /*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)) {
+
+ /*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};
+ 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))
+ /* 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;
@@ -654,10 +667,10 @@
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 */
+ } /* for j */
+ } /* for i */
- fclose(grid_fp);
+ fclose(grid_fp);
printf("..done\n");
printf("nbEvents = %lu\n", (unsigned long)eventList->stream_len());
@@ -667,7 +680,7 @@
long)(((long long)(eventList->stream_len() *
sizeof(AEvent))) >> 20));
fflush(stdout);
-
+
return eventList;
}
@@ -688,8 +701,8 @@
// */
// AMI_STREAM < AEvent > *
// init_event_list(FILE * grid_fp, Viewpoint * vp,
-// GridHeader * hd, ViewOptions viewOptions,
-// IOVisibilityGrid *visgrid) {
+// GridHeader * hd, ViewOptions viewOptions,
+// IOVisibilityGrid *visgrid) {
// printf("computing events..");
// fflush(stdout);
@@ -701,82 +714,82 @@
// /*we do this just to position the pointer after header for reading
-// //the data GridHeader* foo = */
+// //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));
+// for (j = 0; j < hd->ncols; j++) {
-// 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;
+// e.row = i;
+// e.col = j;
+// fscanf(grid_fp, "%f", &(e.elev));
-
-// /*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;
+// 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;
-// 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); */
-// }
+// /*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));
+// (unsigned long)eventList->stream_len(), (int)sizeof(AEvent),
+// (unsigned
+// long)(((long long)(eventList->stream_len() *
+// sizeof(AEvent))) >> 20));
// fflush(stdout);
// return eventList;
// }
@@ -807,19 +820,19 @@
calculate_event_position(a, globalVP.row, globalVP.col, &eventy, &eventx);
da = (eventx - globalVP.col) * (eventx - globalVP.col) +
- (eventy - globalVP.row) * (eventy - globalVP.row);
+ (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);
+ (eventy - globalVP.row) * (eventy - globalVP.row);
if (da > db) {
- return 1;
+ return 1;
}
else if (da < db) {
- return -1;
+ return -1;
}
else {
- return 0;
+ return 0;
}
return 0;
}
Modified: grass-addons/raster/r.viewshed/eventlist.h
===================================================================
--- grass-addons/raster/r.viewshed/eventlist.h 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/eventlist.h 2008-12-11 08:10:58 UTC (rev 34825)
@@ -45,7 +45,7 @@
#ifdef __GRASS__
#include <grass/iostream/ami.h>
-#else
+#else
#include <ami.h>
#endif
@@ -68,12 +68,12 @@
/* ------------------------------------------------------------ */
/*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);
+ 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);
@@ -83,12 +83,12 @@
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);
+init_event_list_in_memory(AEvent * eventList, char *inputfname,
+ Viewpoint * vp,
+ GridHeader * hd, ViewOptions viewOptions,
+ double **data, MemoryVisibilityGrid * visgrid);
@@ -100,11 +100,11 @@
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);
+ */
+AMI_STREAM < AEvent > *init_event_list(char *inputfname, Viewpoint * vp,
+ GridHeader * hd,
+ ViewOptions viewOptions, double **data,
+ IOVisibilityGrid * visgrid);
Modified: grass-addons/raster/r.viewshed/grass.cc
===================================================================
--- grass-addons/raster/r.viewshed/grass.cc 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/grass.cc 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -59,20 +60,25 @@
/* 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;
+ */
+float adjust_for_curvature(Viewpoint vp, dimensionType row,
+ dimensionType col, float h,
+ ViewOptions viewOptions)
+{
- 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));
+ 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));
}
@@ -81,23 +87,26 @@
/* ************************************************************ */
/*return a GridHeader that has all the relevant data filled in from
- GRASS */
-GridHeader *read_header_from_GRASS(char *rastName, Cell_head *region) {
+ 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;
+ hd->nrows = (dimensionType) nrows;
+ hd->ncols = (dimensionType) ncols;
}
else
G_fatal_error(_("Grid dimension too big for current precision"));
@@ -107,17 +116,17 @@
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);
- G_message("Nodata value set to %f", hd->nodata_value);
- return hd;
+ //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);
+ G_message("Nodata value set to %f", hd->nodata_value);
+ return hd;
}
@@ -126,75 +135,82 @@
/* ************************************************************ */
/* 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. */
+ 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
+ Viewpoint * vp, GridHeader * hd,
+ ViewOptions viewOptions, double **data,
+ MemoryVisibilityGrid * visgrid)
+{
- /*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_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);
-
- /*open map */
- int infd;
- if ((infd = G_open_cell_old(rastName, mapset)) < 0)
+
+ /*open map */
+ int infd;
+
+ if ((infd = G_open_cell_old(rastName, mapset)) < 0)
G_fatal_error(_("Cannot open raster file [%s]"), rastName);
-
- /*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++) {
+ /*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);
-
+ 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 < 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:
+ 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;
@@ -206,69 +222,72 @@
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) {
+ /* 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 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_warning(_("Viewpoint is NODATA."));
- G_message(_("Will assume its elevation is = %f"), vp->elev);
+ /*what to do when viewpoint is NODATA ? */
+ G_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) {
+ 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);
+ 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))
+ }
+
+ /* 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++;
-
+
+ /* 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;
+ G_message(_("...done creating event list"));
+ G_close_cell(infd);
+ return nevents;
}
@@ -284,146 +303,155 @@
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);
+AMI_STREAM < AEvent > *grass_init_event_list(char *rastName, Viewpoint * vp,
+ GridHeader * hd,
+ ViewOptions viewOptions,
+ double **data,
+ IOVisibilityGrid * visgrid)
+{
- if (data != NULL) {
+ 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. */
+ //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)
+ }
+
+ /*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);
-
- /*open map */
- int infd;
- if ((infd = G_open_cell_old(rastName, mapset)) < 0)
+
+ /*open map */
+ int infd;
+
+ if ((infd = G_open_cell_old(rastName, mapset)) < 0)
G_fatal_error(_("Cannot open raster file [%s]"), rastName);
- RASTER_MAP_TYPE data_type;
- data_type = G_raster_map_type(rastName, mapset);
- void *inrast;
- inrast = G_allocate_raster_buf(data_type);
- assert(inrast);
+ RASTER_MAP_TYPE data_type;
- /*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++) {
-
+ 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);
+ 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 < 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:
+ 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:
+ case FCELL_TYPE:
isnull = G_is_f_null_value(&(((FCELL *) inrast)[j]));
e.elev = (float)(((FCELL *) inrast)[j]);
break;
- case DCELL_TYPE:
+ 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) {
+ /* 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;
+ (*data)[j] = e.elev;
}
- }
-
- /* set the viewpoint */
- if (i == vp-> row && j == vp->col) {
+ }
+
+ /* 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)) {
- G_warning("Viewpoint is NODATA.");
- G_message("Will assume its elevation is %.f", e.elev);
+ G_warning("Viewpoint is NODATA.");
+ G_message("Will assume its elevation is %.f", e.elev);
};
- }
-
- /*don't insert viewpoint into eventlist */
- if (i == vp->row && j == vp->col)
+ }
+
+ /*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)) {
+
+ /*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};
+ 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))
+ }
+
+ /* 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);
- G_message("nbEvents = %lu", (unsigned long)eventList->stream_len());
- G_message("Event stream length: %lu x %dB (%lu MB)",
- (unsigned long)eventList->stream_len(), (int)sizeof(AEvent),
- (unsigned long)(((long long)(eventList->stream_len() *
- sizeof(AEvent))) >> 20));
- fflush(stdout);
- return eventList;
+
+ /*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);
+ G_message("nbEvents = %lu", (unsigned long)eventList->stream_len());
+ G_message("Event stream length: %lu x %dB (%lu MB)",
+ (unsigned long)eventList->stream_len(), (int)sizeof(AEvent),
+ (unsigned long)(((long long)(eventList->stream_len() *
+ sizeof(AEvent))) >> 20));
+ fflush(stdout);
+ return eventList;
}
@@ -433,44 +461,48 @@
/* ************************************************************ */
/* saves the grid into a GRASS raster. Loops through all elements x
- in row-column order and writes fun(x) to file.*/
+ 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);
+save_grid_to_GRASS(Grid * grid, char *filename, RASTER_MAP_TYPE type,
+ float (*fun) (float))
+{
- /*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);
+ G_message(_("Saving grid to <%s>"), filename);
+ assert(grid && filename);
- dimensionType i, j;
- for (i = 0; i < G_window_rows(); i++) {
+ /*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]);
+
+ 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]);
+ 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]);
+ case DCELL_TYPE:
+ ((DCELL *) outrast)[j] = (DCELL) fun(grid->grid_data[i][j]);
break;
- }
- } /* for j */
+ }
+ } /* for j */
G_put_raster_row(outfd, outrast, type);
- } /* for i */
+ } /* for i */
- G_close_cell(outfd);
- return;
+ G_close_cell(outfd);
+ return;
}
@@ -479,136 +511,150 @@
/* ************************************************************ */
/* 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 */
+ 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)
+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);
- /*open elevation map */
- int elevfd;
- if ((elevfd = G_open_cell_old(elevfname, mapset)) < 0)
+ /*open elevation map */
+ int elevfd;
+
+ if ((elevfd = G_open_cell_old(elevfname, mapset)) < 0)
G_fatal_error(_("Cannot open raster file [%s]"), elevfname);
-
- /*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 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);
+ 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:
+
+ /* 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:
+ case FCELL_TYPE:
isNull = G_is_f_null_value(&((FCELL *) elevrast)[j]);
elev = (double)(((FCELL *) elevrast)[j]);
break;
- case DCELL_TYPE:
+ 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])) {
+ }
+
+ 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])) {
+ 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);
+ assert(!isNull);
/* write INVISIBLE */
- viewshed_value = INVISIBLE;
- writeValue(visrast,j,viewshed_value, elev_data_type);
- } else {
- /* nodata */
+ 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*/
+
+ } /* for j */
G_put_raster_row(visfd, visrast, elev_data_type);
- } /* for i */
-
- G_close_cell(elevfd);
- G_close_cell(visfd);
- return;
+ } /* 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) {
-
+/* 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:
+ case CELL_TYPE:
((CELL *) bufrast)[j] = (CELL) x;
break;
- case FCELL_TYPE:
- ((FCELL *) bufrast)[j] = (FCELL) x;
+ case FCELL_TYPE:
+ ((FCELL *) bufrast)[j] = (FCELL) x;
break;
- case DCELL_TYPE:
- ((DCELL *) bufrast)[j] = (DCELL) x;
+ case DCELL_TYPE:
+ ((DCELL *) bufrast)[j] = (DCELL) x;
break;
- default:
- G_fatal_error(_("Unknown data type"));
- }
+ default:
+ G_fatal_error(_("Unknown data 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);
+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);
+ 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);
+ case DCELL_TYPE:
+ G_set_d_null_value(&((DCELL *) bufrast)[j], 1);
break;
- default:
- G_fatal_error(_("Unknown data type"));
- }
+ default:
+ G_fatal_error(_("Unknown data type"));
+ }
}
@@ -617,181 +663,194 @@
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))
+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) {
+ 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) {
+ 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);
-
+ 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++;
+ ae = vstr->read_item(&curResult);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ counter++;
}
- }
- else {
+ }
+ else {
/* this cell is not in stream, so it is invisible */
writeValue(visrast, j, fun(INVISIBLE), type);
- }
- } /* for j */
-
+ }
+ } /* for j */
+
G_put_raster_row(visfd, visrast, type);
- } /* for i */
-
- G_close_cell(visfd);
-}
+ } /* 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. */
+ 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)
+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);
- /*open elevation map */
- int elevfd;
- if ((elevfd = G_open_cell_old(elevfname, mapset)) < 0)
+ 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);
+
+ /*open elevation map */
+ int elevfd;
+
+ if ((elevfd = G_open_cell_old(elevfname, mapset)) < 0)
G_fatal_error(_("Cannot open raster file [%s]"), elevfname);
-
- /*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) {
+ /*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++) {
+ }
+ 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);
+ G_fatal_error(_("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:
+
+ /* 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:
+ case FCELL_TYPE:
isNull = G_is_f_null_value(&((FCELL *) elevrast)[j]);
elev = (double)(((FCELL *) elevrast)[j]);
break;
- case DCELL_TYPE:
+ 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) {
+
+ 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);
+ 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++;
+ ae = vstr->read_item(&curResult);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ counter++;
}
- }
- else {
+ }
+ 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;
+ }
+ } /* for j */
+
+ G_put_raster_row(visfd, visrast, elev_data_type);
+ } /* for i */
+
+ G_close_cell(elevfd);
+ G_close_cell(visfd);
+ return;
}
Modified: grass-addons/raster/r.viewshed/grass.h
===================================================================
--- grass-addons/raster/r.viewshed/grass.h 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/grass.h 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -58,34 +59,34 @@
/* 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);
+ */
+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);
-
+/* 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. */
+ 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 );
+ Viewpoint * vp, GridHeader * hd,
+ ViewOptions viewOptions, double **data,
+ MemoryVisibilityGrid * visgrid);
@@ -98,30 +99,31 @@
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);
+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.*/
+ 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));
+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 */
+ 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);
+save_vis_elev_to_GRASS(Grid * visgrid, char *elevfname, char *visfname,
+ float vp_elev);
/* ************************************************************ */
@@ -129,22 +131,22 @@
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));
+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. */
+ 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);
+save_io_vis_and_elev_to_GRASS(IOVisibilityGrid * visgrid, char *elevfname,
+ char *visfname, float vp_elev);
#endif/*_GRASS_H*/
Modified: grass-addons/raster/r.viewshed/grid.cc
===================================================================
--- grass-addons/raster/r.viewshed/grid.cc 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/grid.cc 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -56,28 +57,31 @@
/* ------------------------------------------------------------ */
/*read header from file and return it; */
-GridHeader *read_header_from_arcascii_file(char* fname)
+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;
+ 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)
+GridHeader *read_header_from_arcascii_file(FILE * fp)
{
assert(fp);
- rewind(fp);
+ rewind(fp);
GridHeader *hd = (GridHeader *) malloc(sizeof(GridHeader));
+
assert(hd);
int nrows, ncols;
@@ -166,7 +170,7 @@
int is_nodata(Grid * grid, float value)
{
assert(grid);
- return is_nodata(grid->hd, value);
+ return is_nodata(grid->hd, value);
}
@@ -219,16 +223,16 @@
for (i = 0; i < pgrid->hd->nrows; i++) {
#ifdef __GRASS__
- pgrid->grid_data[i] =
+ pgrid->grid_data[i] =
(float *)G_malloc(pgrid->hd->ncols * sizeof(float));
#else
- pgrid->grid_data[i] =
+ pgrid->grid_data[i] =
(float *)malloc(pgrid->hd->ncols * sizeof(float));
#endif
- assert(pgrid->grid_data[i]);
+ assert(pgrid->grid_data[i]);
}
-
+
#ifdef _DEBUG_ON
printf("**DEBUG: allocGridData\n");
fflush(stdout);
@@ -349,49 +353,51 @@
/* ------------------------------------------------------------ */
/*save the grid into an arcascii file. Loops through all elements x
- in row-column order and writes fun(x) to file */
+ in row-column order and writes fun(x) to file */
void
-save_grid_to_arcascii_file(Grid * grid, char *filename,
- float(*fun)(float)) {
+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 */
+ 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);
+ 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++) {
+ }
+ 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]) );
-
+
+ /* 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);
+ printf("**DEBUG: saveGridToArcasciiFile: saved to %s\n", filename);
+ fflush(stdout);
#endif
-
- return;
+
+ return;
}
Modified: grass-addons/raster/r.viewshed/grid.h
===================================================================
--- grass-addons/raster/r.viewshed/grid.h 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/grid.h 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -54,18 +55,20 @@
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 */
+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_ {
+typedef struct grid_
+{
GridHeader *hd;
/*two dimensional array holding all the values in the grid */
@@ -78,11 +81,11 @@
-/* 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(char *fname);
-/* create and return the header of the grid stored in this file;*/
-GridHeader *read_header_from_arcascii_file(FILE* fp);
+/* 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);
@@ -110,9 +113,9 @@
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));
+ row-column order and writes fun(x) to file */
+void save_grid_to_arcascii_file(Grid * grid, char *filename,
+ float (*fun) (float));
Modified: grass-addons/raster/r.viewshed/main.cc
===================================================================
--- grass-addons/raster/r.viewshed/main.cc 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/main.cc 2008-12-11 08:10:58 UTC (rev 34825)
@@ -74,7 +74,7 @@
/* observer elevation above the terrain */
-#define DEFAULT_OBS_ELEVATION 0
+#define DEFAULT_OBS_ELEVATION 0
/* All these flags are used for debugging */
@@ -98,30 +98,35 @@
/* ------------------------------------------------------------ */
/* forward declarations */
/* ------------------------------------------------------------ */
-void print_timings_internal(Rtimer sweepTime, Rtimer outputTime, Rtimer totalTime);
+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();
+ 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);
+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);
+void parse_args(int argc, char *argv[], int *vpRow, int *vpCol,
+ ViewOptions * viewOptions, long long *memSizeBytes);
#endif
-
+
/* ------------------------------------------------------------ */
-int main(int argc, char *argv[]) {
+int main(int argc, char *argv[])
+{
#ifdef __GRASS__
- /* GRASS initialization stuff */
+ /* GRASS initialization stuff */
struct GModule *module;
/*initialize GIS environment */
@@ -133,101 +138,110 @@
module->description = _("IO-efficient viewshed algorithm");
struct Cell_head region;
+
if (G_get_set_window(®ion) == -1)
G_fatal_error("Error getting current region");
#endif
- /* ************************************************************ */
- /* parameters set up */
+ /* ************************************************************ */
+ /* 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 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 */
+ 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;
+ ViewOptions viewOptions;
+ //viewOptions.inputfname = (char*)malloc(500);
+ //viewOptions.outputfname = (char*)malloc(500);
+ //assert(inputfname && outputfname);
+ viewOptions.obsElev = DEFAULT_OBS_ELEVATION;
+ viewOptions.maxDist = INFINITY_DISTANCE;
+ viewOptions.outputMode = OUTPUT_ANGLE;
+ viewOptions.doCurv = 0;
+
#ifdef __GRASS__
- parse_args(argc,argv, &vpRow, &vpCol,&viewOptions, &memSizeBytes,®ion);
+ parse_args(argc, argv, &vpRow, &vpCol, &viewOptions, &memSizeBytes,
+ ®ion);
#else
- parse_args(argc,argv, &vpRow, &vpCol, &viewOptions, &memSizeBytes);
+ 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 */
+ 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);
-
-
+ 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;
+ sure the requested viewpoint is on the map */
+ GridHeader *hd;
+
#ifdef __GRASS__
hd = read_header_from_GRASS(viewOptions.inputfname, ®ion);
assert(hd);
- /* LT: there is no need to exit if viewpoint is outside grid,
- the algorithm will work correctly in theory. But this
- requires some changes. To do.*/
+ /* 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_warning(_("Viewpoint outside grid"));
- G_warning(_("viewpont: (row=%d, col=%d)"), vp.row, vp.col);
- G_fatal_error(_("grid: (rows=%d, cols=%d)"), hd->nrows, hd->ncols);
+ G_warning(_("Viewpoint outside grid"));
+ G_warning(_("viewpont: (row=%d, col=%d)"), vp.row, vp.col);
+ G_fatal_error(_("grid: (rows=%d, cols=%d)"), hd->nrows, hd->ncols);
}
#else
- /*open file input file and read grid header from grid ascii file */
+ /*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.*/
- }
+ 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 */
+
+ /* 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, using sphere (6370997.0)"));
- viewOptions.ellps_a = 6370997.00;
- }
+ 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, using 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;
+ /* in standalone mode we do not know how to adjust for curvature */
+ assert(viewOptions.doCurv == 0);
+ viewOptions.ellps_a = 0;
#endif
@@ -237,18 +251,17 @@
/* 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);
+ 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));
+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");
+ IN_MEMORY = 1;
+ print_message("*************\nIN_MEMORY MODE\n*************\n");
}
else {
- print_message("*************\nEXTERNAL_MEMORY MODE\n**********\n");
- IN_MEMORY = 0;
+ 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
@@ -268,27 +281,28 @@
/* 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);
+ /*//////////////////////////////////////////////////// */
+ /*/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);
}
@@ -298,123 +312,130 @@
/* 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)
+
+ /* ************************************************************ */
+ /* 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 */
+
+ /*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;
+ BASE_CASE = 0;
#endif
-
-
-
-
- /* ************************************************************ */
- /* external memory, base case */
- /* ************************************************************ */
- if (BASE_CASE) {
- print_message
- ("---Active structure small, starting base case---\n");
-
+
+
+
+
+ /* ************************************************************ */
+ /* 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);
+ 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);
- }
-
-
-
-
+ outputTime, sortOutputTime);
+ }
+
+
+
+
/************************************************************/
- /* external memory, recursive distribution sweeping recursion */
+ /* external memory, recursive distribution sweeping recursion */
+
/************************************************************ */
- else { /* if not BASE_CASE */
+ 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;
-
+
+ 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);
-
+ 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 */
+ outputted */
rt_start(sortOutputTime);
sort_io_visibilitygrid(visgrid);
rt_stop(sortOutputTime);
@@ -422,36 +443,37 @@
rt_start(outputTime);
save_io_visibilitygrid(visgrid, viewOptions, vp);
rt_stop(outputTime);
-
-
+
+
rt_stop(totalTime);
print_timings_external_memory(totalTime, sweepTime,
- outputTime, sortOutputTime);
+ 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 */
+ /*following GRASS's coding standards for history and exiting */
struct History history;
- G_short_history(viewOptions.outputfname, "raster", &history);
+
+ G_short_history(viewOptions.outputfname, "raster", &history);
G_command_history(&history);
G_write_history(viewOptions.outputfname, &history);
- exit(EXIT_SUCCESS);
+ exit(EXIT_SUCCESS);
#else
- free(hd);
+ free(hd);
#endif
}
@@ -466,101 +488,115 @@
#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);
+ ViewOptions * viewOptions, long long *memSizeBytes,
+ Cell_head * window)
+{
- /* the input */
- struct Option *inputOpt;
- inputOpt = G_define_standard_option(G_OPT_R_ELEV);
- inputOpt->key = "input";
- inputOpt->guisection = _("Input_options");
+ assert(vpRow && vpCol && memSizeBytes && window);
- /* the output */
- struct Option *outputOpt;
- outputOpt = G_define_standard_option(G_OPT_R_OUTPUT);
- outputOpt->description = _("Name of output viewshed raster map\n\t\tdefault format: {NODATA, -1 (invisible), vertical angle wrt viewpoint (visible)}");
- outputOpt->guisection = _("Output_options");
+ /* the input */
+ struct Option *inputOpt;
- /* row-column flag */
- struct Flag *row_col;
- row_col = G_define_flag();
- row_col->key = 'r';
- row_col->description =
+ inputOpt = G_define_standard_option(G_OPT_R_ELEV);
+ inputOpt->key = "input";
+ inputOpt->guisection = _("Input_options");
+
+ /* the output */
+ struct Option *outputOpt;
+
+ outputOpt = G_define_standard_option(G_OPT_R_OUTPUT);
+ outputOpt->description =
+ _("Name of output viewshed raster map\n\t\tdefault format: {NODATA, -1 (invisible), vertical angle wrt viewpoint (visible)}");
+ outputOpt->guisection = _("Output_options");
+
+ /* 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 =
+
+ /* 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 =
+
+
+ /* 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 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 = "coordinate";
- viewLocOpt->type = TYPE_STRING;
- viewLocOpt->required = YES;
- viewLocOpt->key_desc = "lat,long";
- viewLocOpt->description =
+ /* viewpoint coordinates */
+ struct Option *viewLocOpt;
+
+ viewLocOpt = G_define_option();
+ viewLocOpt->key = "coordinate";
+ 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)");
- viewLocOpt->guisection = _("Input_options");
-
- /* observer elevation */
- struct Option *obsElevOpt;
- obsElevOpt = G_define_option();
- obsElevOpt->key = "obs_elev";
- obsElevOpt->type = TYPE_DOUBLE;
- obsElevOpt->required = NO;
- obsElevOpt->key_desc = "value";
- obsElevOpt->description = _("Viewing elevation above the ground");
- obsElevOpt->answer = "1.75";
- obsElevOpt->guisection = _("Input_options");
-
- /* 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 =
+ viewLocOpt->guisection = _("Input_options");
+
+ /* observer elevation */
+ struct Option *obsElevOpt;
+
+ obsElevOpt = G_define_option();
+ obsElevOpt->key = "obs_elev";
+ obsElevOpt->type = TYPE_DOUBLE;
+ obsElevOpt->required = NO;
+ obsElevOpt->key_desc = "value";
+ obsElevOpt->description = _("Viewing elevation above the ground");
+ obsElevOpt->answer = "1.75";
+ obsElevOpt->guisection = _("Input_options");
+
+ /* 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;
- maxDistOpt->guisection = _("Input_options");
+ char infdist[10];
- /* memory size */
+ sprintf(infdist, "%d", INFINITY_DISTANCE);
+ maxDistOpt->answer = infdist;
+ maxDistOpt->guisection = _("Input_options");
+
+ /* 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->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);
+ exit(EXIT_FAILURE);
/* store the parameters into a structure to be used along the way */
@@ -570,46 +606,47 @@
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"));
+ if (viewOptions->maxDist < 0 && viewOptions->maxDist != INFINITY_DISTANCE) {
+ G_fatal_error(_("Negative max distance value is not valid"));
}
-
- viewOptions->doCurv = curvature->answer;
- if (booleanOutput->answer)
- viewOptions->outputMode = OUTPUT_BOOL;
- else if (elevationFlag->answer)
- viewOptions->outputMode = OUTPUT_ELEV;
- else viewOptions->outputMode = OUTPUT_ANGLE;
+ 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);
- }
+
+ 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 */
+ 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);
+ *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);
-
- }
+ *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;
}
@@ -619,93 +656,96 @@
/* ------------------------------------------------------------ */
#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);
+parse_args(int argc, char *argv[], int *vpRow, int *vpCol,
+ ViewOptions * viewOptions, long long *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) {
+ 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;
+ /* inputfile name */
+ //*inputfname = optarg;
+ strcpy(viewOptions->inputfname, optarg);
+ gotinput = 1;
+ break;
case 'o':
- //*outputfname = optarg;
- strcpy(viewOptions->outputfname, optarg);
- gotoutput= 1;
- break;
+ //*outputfname = optarg;
+ strcpy(viewOptions->outputfname, optarg);
+ gotoutput = 1;
+ break;
case 'r':
*vpRow = atoi(optarg);
- gotrow=1;
+ gotrow = 1;
break;
case 'c':
*vpCol = atoi(optarg);
- gotcol = 1;
+ 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;
+ /* 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;
+ viewOptions->obsElev = atof(optarg);
+ break;
case 'd':
- viewOptions->maxDist = atof(optarg);
- if (viewOptions->maxDist < 0) {
+ viewOptions->maxDist = atof(optarg);
+ if (viewOptions->maxDist < 0) {
printf("max distance needs to be positive\n");
- exit(EXIT_FAILURE);
- }
- break;
+ exit(EXIT_FAILURE);
+ }
+ break;
case 'm':
- int memSizeMB;
- memSizeMB = atoi(optarg);
- if (memSizeMB <0) {
+ int memSizeMB;
+
+ memSizeMB = atoi(optarg);
+ if (memSizeMB < 0) {
printf("Memory cannot be negative.\n");
exit(1);
- }
- *memSizeBytes = (long long)memSizeMB;
- *memSizeBytes = (*memSizeBytes) << 20;
- break;
+ }
+ *memSizeBytes = (long long)memSizeMB;
+ *memSizeBytes = (*memSizeBytes) << 20;
+ break;
case '?':
- if (optopt == 'i' || optopt == 'o' || optopt == 'r' ||
- optopt == 'c' || optopt == 'e' || optopt == 'd' ||
- optopt == 'm')
+ 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))
+ else if (isprint(optopt))
fprintf(stderr, "Unknown option '-%c'.\n", optopt);
- else
+ else
fprintf(stderr, "Unknown option character '\\x%x.\n", optopt);
- print_usage();
- exit(1);
- //default:
- //exit(1);
+ print_usage();
+ exit(1);
+ //default:
+ //exit(1);
}
- } /* while getopt */
-
- /*check to make sure the required options are set.*/
- if(!(gotinput && gotoutput && gotrow &&gotcol)) {
+ } /* 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;
+ }
+ printf("viewpoint: (%d, %d)\n", *vpRow, *vpCol);
+ return;
}
#endif
@@ -714,24 +754,28 @@
/* ------------------------------------------------------------ */
/*print the timings for the internal memory method of computing the
- viewshed */
+ 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");
+print_timings_internal(Rtimer sweepTime, Rtimer outputTime, Rtimer totalTime)
+{
- rt_sprint_safe(timeused, totalTime);
- printf("\t%30s", "total:");
- printf(timeused);
- printf("\n");
+ 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");
}
@@ -766,91 +810,100 @@
}
/* ------------------------------------------------------------ */
-void print_status(Viewpoint vp,ViewOptions viewOptions, long long memSizeBytes)
+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"));
+ 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("---------------------------------\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("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");
-
+ 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");
+/*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");
+ 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");
}
Modified: grass-addons/raster/r.viewshed/rbbst.h
===================================================================
--- grass-addons/raster/r.viewshed/rbbst.h 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/rbbst.h 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -40,7 +41,7 @@
#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*/
+ key within that distance. The largest double value is 1.7 E 308 */
Modified: grass-addons/raster/r.viewshed/statusstructure.cc
===================================================================
--- grass-addons/raster/r.viewshed/statusstructure.cc 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/statusstructure.cc 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -59,22 +60,24 @@
/* ------------------------------------------------------------ */
/*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);
+ 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);
}
@@ -84,15 +87,15 @@
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;
+ 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;
}
@@ -133,14 +136,14 @@
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;
}
@@ -211,4 +214,3 @@
return (is_empty(sl->rbt) ||
sl->rbt->root->value.maxGradient == SMALLEST_GRADIENT);
}
-
Modified: grass-addons/raster/r.viewshed/statusstructure.h
===================================================================
--- grass-addons/raster/r.viewshed/statusstructure.h 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/statusstructure.h 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -99,8 +100,8 @@
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.*/
+ 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);
Modified: grass-addons/raster/r.viewshed/viewshed.cc
===================================================================
--- grass-addons/raster/r.viewshed/viewshed.cc 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/viewshed.cc 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -56,39 +57,45 @@
/* ------------------------------------------------------------ */
/* return the memory usage (in bytes) of viewshed */
-long long get_viewshed_memory_usage(GridHeader* hd) {
+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;
- 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("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);
-
+ assert(hd);
+ /* the output visibility grid */
+ long long totalcells = (long long)hd->nrows * (long long)hd->ncols;
+
+ 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("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)
+ Rtimer sortEventTime, Rtimer sweepTime)
{
char timeused[1000];
@@ -128,45 +135,50 @@
possible that the amount of memory required is more than that
supported by the platform; for e.g. on a 32-bt platform cannot
allocate more than 4GB. Try to detect this situation. */
-AEvent*
-allocate_eventlist(GridHeader* hd) {
-
- AEvent* eventList;
+AEvent *allocate_eventlist(GridHeader * hd)
+{
- 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));
+ AEvent *eventList;
- /* what's the size of size_t on this machine? */
- int sizet_size;
- sizet_size = (int) sizeof(size_t);
- printf("size_t is %lu B\n", sizeof(size_t));
-
- if (sizet_size >= 8) {
- printf("64-bit platform, great.\n");
- } else {
- /* this is the max value of size_t */
- long long maxsizet = ((long long)1<<(sizeof(size_t)*8)) -1;
- printf("max size_t is %lld\n", maxsizet);
-
- /* checking whether allocating totalsize causes an overflow */
- 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);
+ 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));
+
+ /* what's the size of size_t on this machine? */
+ int sizet_size;
+
+ sizet_size = (int)sizeof(size_t);
+ printf("size_t is %lu B\n", sizeof(size_t));
+
+ if (sizet_size >= 8) {
+ printf("64-bit platform, great.\n");
}
- }
-
- printf("allocating..");
+ else {
+ /* this is the max value of size_t */
+ long long maxsizet = ((long long)1 << (sizeof(size_t) * 8)) - 1;
+
+ printf("max size_t is %lld\n", maxsizet);
+
+ /* checking whether allocating totalsize causes an overflow */
+ 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);
+ eventList = (AEvent *) G_malloc(totalsize);
#else
- eventList = (AEvent *) malloc(totalsize*sizeof(char));
+ eventList = (AEvent *) malloc(totalsize * sizeof(char));
#endif
- assert(eventList);
- printf("..ok\n");
+ assert(eventList);
+ printf("..ok\n");
- return eventList;
+ return eventList;
}
@@ -186,31 +198,35 @@
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){
-
+
+ */
+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);
+ 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));
- }
+ 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 */
@@ -218,43 +234,48 @@
size_t nevents;
Rtimer initEventTime;
+
rt_start(initEventTime);
- AEvent *eventList = allocate_eventlist(hd);
+ AEvent *eventList = allocate_eventlist(hd);
+
#ifdef __GRASS__
- nevents = grass_init_event_list_in_memory(eventList, inputfname, vp, hd,
- viewOptions, &data, visgrid);
+ 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);
+ 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);
+ 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 */
+ //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);
+ printf("done\n");
+ fflush(stdout);
rt_stop(sortEventTime);
-
- /* ------------------------------ */
+
+ /* ------------------------------ */
/*create the status structure */
StatusList *status_struct = create_status_struct();
@@ -264,22 +285,22 @@
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)) {
+ 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_statusnode(sn);
- printf("\n");
+ printf("inserting: ");
+ print_statusnode(sn);
+ printf("\n");
}
/*insert sn into the status structure */
insert_into_status_struct(sn, status_struct);
- }
+ }
}
#ifdef __GRASS__
G_free(data);
@@ -289,78 +310,82 @@
- /* ------------------------------ */
+ /* ------------------------------ */
/*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);
- INMEMORY_DEBUG {
+ /*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);
+ INMEMORY_DEBUG {
printf("event: ");
print_event(*e);
- printf("sn.dist=%f, sn.gradient=%f\n", sn.dist2vp, sn.gradient);
- }
-
- switch (e->eventType) {
- case ENTERING_EVENT:
+ printf("sn.dist=%f, sn.gradient=%f\n", sn.dist2vp, sn.gradient);
+ }
+
+ switch (e->eventType) {
+ case ENTERING_EVENT:
/*insert node into structure */
VIEWSHEDDEBUG {
- printf("..ENTER-EVENT: insert\n");
+ printf("..ENTER-EVENT: insert\n");
}
insert_into_status_struct(sn, status_struct);
break;
-
- case EXITING_EVENT:
+
+ case EXITING_EVENT:
/*delete node out of status structure */
VIEWSHEDDEBUG {
- printf("..EXIT-EVENT: delete\n");
+ printf("..EXIT-EVENT: delete\n");
}
delete_from_status_struct(status_struct, sn.dist2vp);
break;
-
- case CENTER_EVENT:
+
+ case CENTER_EVENT:
VIEWSHEDDEBUG {
- printf("..QUERY-EVENT: query\n");
+ printf("..QUERY-EVENT: query\n");
}
/*calculate visibility */
double max;
- max=find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
-
+
+ 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++;
+ 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;
- //}
+ /* 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 *
+ (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);
@@ -371,7 +396,7 @@
#else
free(eventList);
#endif
-
+
return visgrid;
}
@@ -391,82 +416,88 @@
stored as streams
*/
-IOVisibilityGrid *viewshed_external(char* inputfname, GridHeader * hd,
- Viewpoint* vp, ViewOptions viewOptions){
+IOVisibilityGrid *viewshed_external(char *inputfname, GridHeader * hd,
+ Viewpoint * vp, ViewOptions viewOptions)
+{
- assert(inputfname && hd && vp);
- printf("Start sweeping.\n");
- fflush(stdout);
+ 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);
+ /* ------------------------------ */
+ /*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);
+ eventList = grass_init_event_list(inputfname, vp, hd, viewOptions,
+ &data, visgrid);
#else
- eventList = init_event_list(inputfname, vp,hd,viewOptions,&data,visgrid);
+ 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++) {
+ 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 {
+ 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_statusnode(sn);
printf("\n");
- }
- /*insert sn into the status structure */
- insert_into_status_struct(sn, status_struct);
+ }
+ /*insert sn into the status structure */
+ insert_into_status_struct(sn, status_struct);
}
- }
+ }
#ifdef __GRASS__
- G_free(data);
+ G_free(data);
#else
- free(data);
+ free(data);
#endif
-
-
- /* ------------------------------ */
+
+
+ /* ------------------------------ */
/*sweep the event list */
long nvis = 0; /*number of visible cells */
VisCell viscell;
@@ -481,75 +512,78 @@
/*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_statusnode(sn);
+ printf("next event: ");
+ print_statusnode(sn);
}
-
+
switch (e->eventType) {
case ENTERING_EVENT:
- /*insert node into structure */
- VIEWSHEDDEBUG {
+ /*insert node into structure */
+ VIEWSHEDDEBUG {
printf("..ENTER-EVENT: insert\n");
- }
- insert_into_status_struct(sn, status_struct);
- break;
-
+ }
+ insert_into_status_struct(sn, status_struct);
+ break;
+
case EXITING_EVENT:
- /*delete node out of status structure */
- VIEWSHEDDEBUG {
+ /*delete node out of status structure */
+ VIEWSHEDDEBUG {
printf("..EXIT-EVENT: delete\n");
- }
- delete_from_status_struct(status_struct, sn.dist2vp);
- break;
-
+ }
+ delete_from_status_struct(status_struct, sn.dist2vp);
+ break;
+
case CENTER_EVENT:
- VIEWSHEDDEBUG {
+ 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);
+ }
+ /*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 {
/* 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;
+ }
+ break;
}
- } /* for each event */
+ } /* 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)));
-
+ (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;
}
Modified: grass-addons/raster/r.viewshed/viewshed.h
===================================================================
--- grass-addons/raster/r.viewshed/viewshed.h 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/viewshed.h 2008-12-11 08:10:58 UTC (rev 34825)
@@ -54,8 +54,8 @@
/* ------------------------------------------------------------ */
/*return the memory usage in bytes of the algorithm when running in
- memory*/
-long long get_viewshed_memory_usage(GridHeader* hd);
+ memory */
+long long get_viewshed_memory_usage(GridHeader * hd);
@@ -69,41 +69,40 @@
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);
+ */
+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);
+ */
+IOVisibilityGrid *viewshed_external(char *inputfname,
+ GridHeader * hd,
+ Viewpoint * vp, ViewOptions viewOptions);
+
void print_viewshed_timings(Rtimer initEventTime, Rtimer sortEventTime,
- Rtimer sweepTime);
+ Rtimer sweepTime);
Modified: grass-addons/raster/r.viewshed/visibility.cc
===================================================================
--- grass-addons/raster/r.viewshed/visibility.cc 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/visibility.cc 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -92,22 +93,24 @@
/* MemoryVisibilityGrid functions */
/* create and return a grid of the sizes specified in the header */
-MemoryVisibilityGrid *create_inmem_visibilitygrid(GridHeader hd, Viewpoint vp){
-
- MemoryVisibilityGrid *visgrid;
+MemoryVisibilityGrid *create_inmem_visibilitygrid(GridHeader hd, Viewpoint vp)
+{
+
+ MemoryVisibilityGrid *visgrid;
+
#ifdef __GRASS__
- visgrid = (MemoryVisibilityGrid *) G_malloc(sizeof(MemoryVisibilityGrid));
+ visgrid = (MemoryVisibilityGrid *) G_malloc(sizeof(MemoryVisibilityGrid));
#else
- visgrid = (MemoryVisibilityGrid *) malloc(sizeof(MemoryVisibilityGrid));
+ visgrid = (MemoryVisibilityGrid *) malloc(sizeof(MemoryVisibilityGrid));
#endif
- assert(visgrid);
+ assert(visgrid);
-
- /* create the grid */
- visgrid->grid = create_empty_grid();
- assert(visgrid->grid);
-
- /* create the header */
+
+ /* 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
@@ -115,12 +118,12 @@
#endif
assert(visgrid->grid->hd);
- /* set the header */
+ /* 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));
@@ -142,20 +145,20 @@
assert(visgrid);
- if (visgrid->grid) {
- destroy_grid(visgrid->grid);
+ if (visgrid->grid) {
+ destroy_grid(visgrid->grid);
}
#ifdef __GRASS__
- if (visgrid->vp) {
- G_free(visgrid->vp);
+ if (visgrid->vp) {
+ G_free(visgrid->vp);
}
- G_free(visgrid);
-
+ G_free(visgrid);
+
#else
- if (visgrid->vp) {
- free(visgrid->vp);
+ if (visgrid->vp) {
+ free(visgrid->vp);
}
- free(visgrid);
+ free(visgrid);
#endif
return;
}
@@ -163,102 +166,111 @@
/* ------------------------------------------------------------ */
-/*set all values of visgrid's Grid to the given value*/
+/*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 && 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;
+ visgrid->grid->grid_data[i][j] = val;
}
- }
- return;
+ }
+ 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)
+/*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;
+ 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]);
- return;
+ 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.
+ 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.
+ 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:
+ The value computed during the viewshed is the following:
- x is NODATA if the cell is NODATA;
+ x is NODATA if the cell is NODATA;
- x is INVISIBLE if the cell is invisible;
+ 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 */
+ 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);
+ int isnull = G_is_f_null_value(&x);
+
+ if (isnull)
+ return 0;
+ else
+ return (x >= 0);
#else
- return (x >=0);
+ return (x >= 0);
#endif
}
-int is_invisible_not_nodata(float x) {
+int is_invisible_not_nodata(float x)
+{
- return ( (int)x == (int)INVISIBLE);
+ return ((int)x == (int)INVISIBLE);
}
-int is_invisible_nodata(float x) {
+int is_invisible_nodata(float x)
+{
- return (!is_visible(x)) && (!is_invisible_not_nodata(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;
+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;
+float angleVisibilityOutput(float x)
+{
+
+ return x;
}
@@ -270,42 +282,43 @@
/* 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)
+ 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
+ angleVisibilityOutput);
+ else
/* elevation output */
- save_vis_elev_to_GRASS(visgrid->grid, viewOptions.inputfname,
- viewOptions.outputfname,
- vp.elev + viewOptions.obsElev);
+ 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 {
+ 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);
+ 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;
+ free_inmem_visibilitygrid(visgrid);
+
+ return;
}
@@ -392,7 +405,8 @@
/* ------------------------------------------------------------ */
/*write cell to stream */
-void add_result_to_io_visibilitygrid(IOVisibilityGrid * visgrid, VisCell * cell)
+void add_result_to_io_visibilitygrid(IOVisibilityGrid * visgrid,
+ VisCell * cell)
{
assert(visgrid && cell);
@@ -453,43 +467,44 @@
/* ------------------------------------------------------------ */
void
-save_io_visibilitygrid(IOVisibilityGrid * visgrid,
- ViewOptions viewOptions, Viewpoint vp) {
+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
+#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);
+ 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 {
+ 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);
+ 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;
+ /*free visibiliyty grid */
+ free_io_visibilitygrid(visgrid);
+
+ return;
}
@@ -497,68 +512,69 @@
/* ------------------------------------------------------------ */
/*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) */
+ 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) {
+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 i = 0; i < visgrid->hd->nrows; i++) {
for (dimensionType j = 0; j < visgrid->hd->ncols; j++) {
-
- if (curResult->row == i && curResult->col == 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));
-
+ 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++;
+ ae = vstr->read_item(&curResult);
+ assert(ae == AMI_ERROR_NO_ERROR);
+ counter++;
}
- }
- else {
+ }
+ else {
/* this cell is not in stream, then it is invisible */
- fprintf(fp, "%.1f ", fun(INVISIBLE));
- }
- } /* for j */
+ 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;
+ } /* for i */
+
+ /* assert that we read teh entire file, otherwise something is
+ wrong */
+ assert(counter == streamLen);
+ fclose(fp);
+ return;
}
-
-
Modified: grass-addons/raster/r.viewshed/visibility.h
===================================================================
--- grass-addons/raster/r.viewshed/visibility.h 2008-12-11 08:07:02 UTC (rev 34824)
+++ grass-addons/raster/r.viewshed/visibility.h 2008-12-11 08:10:58 UTC (rev 34825)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.viewshed
@@ -42,7 +43,7 @@
#include <grass/config.h>
#include <grass/iostream/ami.h>
-#else
+#else
#include <ami.h>
#endif
@@ -55,121 +56,128 @@
-typedef struct viewpoint_ {
- dimensionType row, col;
- float elev;
+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
+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;
+typedef struct visCell_
+{
+ dimensionType row;
+ dimensionType col;
+ /* VisMode vis; */
+ float angle;
} VisCell;
-typedef enum outputMode_{
- OUTPUT_ANGLE = 0,
- OUTPUT_BOOL = 1,
- OUTPUT_ELEV = 2
+typedef enum outputMode_
+{
+ OUTPUT_ANGLE = 0,
+ OUTPUT_BOOL = 1,
+ OUTPUT_ELEV = 2
} OutputMode;
-typedef struct viewOptions_ {
+typedef struct viewOptions_
+{
- /* the name of the input raster */
- char inputfname [100];
+ /* the name of the input raster */
+ char inputfname[100];
- /* the name of the output raster */
- char outputfname[100];
+ /* the name of the output raster */
+ char outputfname[100];
- float obsElev;
- /* observer elevation above the terrain */
+ 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;
+ float maxDist;
+ /* points that are farther than this distance from the viewpoint are
+ not visible */
- in angle mode, the values recorded are {NODATA, INVISIBLE, angle}
+ OutputMode outputMode;
+ /* The mode the viewshed is output;
- in boolean mode, the values recorded are {BOOL_INVISIBLE, BOOL_VISIBLE}
+ in angle mode, the values recorded are {NODATA, INVISIBLE, angle}
- in elev mode, the values recorded are {NODATA, INVISIBLE, elevation}
- */
+ in boolean mode, the values recorded are {BOOL_INVISIBLE, BOOL_VISIBLE}
- int doCurv;
- /*determines if the curvature of the earth should be considered
- when calculating. Only implemented for GRASS version. */
+ in elev mode, the values recorded are {NODATA, INVISIBLE, elevation}
+ */
- double ellps_a; /* the parameter of the ellipsoid */
- float cellsize; /* the cell resolution */
-} ViewOptions;
+ 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 the ellipsoid */
+ float cellsize; /* the cell resolution */
+} ViewOptions;
+
/*memory visibility grid */
-typedef struct memory_visibility_grid_ {
- Grid *grid;
- Viewpoint *vp;
+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;
+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.
+ 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).
-*/
+ 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);
+ 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);
+float booleanVisibilityOutput(float x);
/* This function is called when the program runs in
viewOptions.outputMode == OUTPUT_ANGLE. */
-float angleVisibilityOutput(float x);
+float angleVisibilityOutput(float x);
@@ -193,18 +201,19 @@
/* ------------------------------------------------------------ */
/* MemoryVisibilityGrid functions */
-MemoryVisibilityGrid *create_inmem_visibilitygrid(GridHeader hd, Viewpoint vp);
+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 add_result_to_inmem_visibilitygrid(MemoryVisibilityGrid * visgrid,
+ dimensionType i, dimensionType j,
+ float val);
void save_inmem_visibilitygrid(MemoryVisibilityGrid * vigrid,
- ViewOptions viewopt, Viewpoint vp);
+ ViewOptions viewopt, Viewpoint vp);
/* ------------------------------------------------------------ */
@@ -217,27 +226,27 @@
void free_io_visibilitygrid(IOVisibilityGrid * grid);
/*write cell to stream */
-void add_result_to_io_visibilitygrid(IOVisibilityGrid * visgrid,
- VisCell * cell);
+void add_result_to_io_visibilitygrid(IOVisibilityGrid * visgrid,
+ VisCell * cell);
/*void
- addResult(IOVisibilityGrid* visgrid, DimensionType row, DimensionType col,
- VisMode vis);
-*/
+ 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);
+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));
+save_io_visibilitygrid_to_arcascii(IOVisibilityGrid * visgrid,
+ char *outfname, float (*fun) (float));
More information about the grass-commit
mailing list