/**************************************************************************** * * MODULE: r.viewshed * * AUTHOR(S): Laura Toma, Bowdoin College - ltoma@bowdoin.edu * Yi Zhuang - yzhuang@bowdoin.edu * * Ported to GRASS by William Richard - * wkrichar@bowdoin.edu or willster3021@gmail.com * Markus Metz: surface interpolation * * Date: april 2011 * * PURPOSE: To calculate the viewshed (the visible cells in the * raster) for the given viewpoint (observer) location. The * visibility model is the following: Two points in the raster are * considered visible to each other if the cells where they belong are * visible to each other. Two cells are visible to each other if the * line-of-sight that connects their centers does not intersect the * terrain. The terrain is NOT viewed as a tesselation of flat cells, * i.e. if the line-of-sight does not pass through the cell center, * elevation is determined using bilinear interpolation. * The viewshed algorithm is efficient both in * terms of CPU operations and I/O operations. It has worst-case * complexity O(n lg n) in the RAM model and O(sort(n)) in the * I/O-model. For the algorithm and all the other details see the * paper: "Computing Visibility on * Terrains in External Memory" by * Herman Haverkort, Laura Toma and Yi Zhuang. * * COPYRIGHT: (C) 2008 by the GRASS Development Team * * This program is free software under the GNU General Public License * (>=v2). Read the file COPYING that comes with GRASS for details. * ****************************************************************************/ #include #include #include extern "C" { #include "grass/gis.h" #include "grass/glocale.h" } #include "viewshed.h" #include "visibility.h" #include "eventlist.h" #include "statusstructure.h" #include "grass.h" #define VIEWSHEDDEBUG if(0) #define INMEMORY_DEBUG if(0) /* ------------------------------------------------------------ */ /* return the memory usage (in bytes) of viewshed */ long long get_viewshed_memory_usage(GridHeader * hd) { assert(hd); /* the output visibility grid */ long long totalcells = (long long)hd->nrows * (long long)hd->ncols; G_verbose_message(_("rows=%d, cols=%d, total = %lld"), hd->nrows, hd->ncols, totalcells); long long gridMemUsage = totalcells * sizeof(float); G_debug(1, "grid usage=%lld", gridMemUsage); /* the event array */ long long eventListMemUsage = totalcells * 3 * sizeof(AEvent); G_debug(1, "memory_usage: eventList=%lld", eventListMemUsage); /* the double array that stores all the cells in the same row as the viewpoint */ long long dataMemUsage = (long long)(hd->ncols * sizeof(double)); G_debug(1, "viewshed memory usage: size AEvent=%dB, nevents=%lld, \ total=%lld B (%d MB)", (int)sizeof(AEvent), totalcells * 3, gridMemUsage + eventListMemUsage + dataMemUsage, (int)((gridMemUsage + eventListMemUsage + dataMemUsage) >> 20)); return (gridMemUsage + eventListMemUsage + dataMemUsage); } /* ------------------------------------------------------------ */ void print_viewshed_timings(Rtimer initEventTime, Rtimer sortEventTime, Rtimer sweepTime) { char timeused[1000]; G_verbose_message(_("Sweep timings:")); rt_sprint_safe(timeused, initEventTime); G_verbose_message("Init events: %s", timeused); rt_sprint_safe(timeused, sortEventTime); G_verbose_message("Sort events: %s", timeused); rt_sprint_safe(timeused, sweepTime); G_verbose_message("Process events: %s", timeused); return; } /* ------------------------------------------------------------ */ static void print_statusnode(StatusNode sn) { G_debug(3, "processing (row=%d, col=%d, dist=%f, grad=%f)", sn.row, sn.col, sn.dist2vp, sn.gradient[1]); return; } /* ------------------------------------------------------------ */ /* allocates the eventlist array used by kreveled_in_memory; it is 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; long long totalsize = hd->ncols * hd->nrows * 3; totalsize *= sizeof(AEvent); G_debug(1, "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); G_debug(1, "size_t is %d B", sizet_size); if (sizet_size >= 8) { G_debug(1, "64-bit platform, great."); } else { /* this is the max value of size_t */ long long maxsizet = ((long long)1 << (sizeof(size_t) * 8)) - 1; G_debug(1, "max size_t is %lld", maxsizet); /* checking whether allocating totalsize causes an overflow */ if (totalsize > maxsizet) { G_fatal_error(_("Running the program in-memory mode requires " \ "memory beyond the capability of the platform. " \ "Use external mode, or a 64-bit platform.")); } } G_debug(1, "allocating eventList..."); eventList = (AEvent *) G_malloc(totalsize); assert(eventList); G_debug(1, "...ok"); return eventList; } /*/////////////////////////////////////////////////////////// ------------------------------------------------------------ run Viewshed's sweep algorithm on the grid stored in the given file, and with the given viewpoint. Create a visibility grid and return it. The computation runs in memory, which means the input grid, the status structure and the output grid are stored in arrays in memory. The output: A cell x in the visibility grid is recorded as follows: if it is NODATA, then x is set to NODATA if it is invisible, then x is set to INVISIBLE if it is visible, then x is set to the vertical angle wrt to viewpoint */ MemoryVisibilityGrid *viewshed_in_memory(char *inputfname, GridHeader * hd, Viewpoint * vp, ViewOptions viewOptions) { assert(inputfname && hd && vp); G_verbose_message(_("Start sweeping.")); /* ------------------------------ */ /* create the visibility grid */ MemoryVisibilityGrid *visgrid; visgrid = create_inmem_visibilitygrid(*hd, *vp); /* set everything initially invisible */ set_inmem_visibilitygrid(visgrid, INVISIBLE); assert(visgrid); G_debug(1, "visibility grid size: %d x %d x %d B (%d MB)", 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 */ surface_type **data; size_t nevents; Rtimer initEventTime; rt_start(initEventTime); AEvent *eventList = allocate_eventlist(hd); nevents = init_event_list_in_memory(eventList, inputfname, vp, hd, viewOptions, &data, visgrid); assert(data); rt_stop(initEventTime); G_debug(1, "actual nb events is %lu", (long unsigned int)nevents); /* ------------------------------ */ /*sort the events radially by angle */ Rtimer sortEventTime; rt_start(sortEventTime); G_verbose_message(_("Sorting events...")); fflush(stdout); /*this is recursive and seg faults for large arrays //qsort(eventList, nevents, sizeof(AEvent), radial_compare_events); //this is too slow... //heapsort(eventList, nevents, sizeof(AEvent), radial_compare_events); //iostream quicksort */ RadialCompare cmpObj; quicksort(eventList, nevents, cmpObj); G_verbose_message(_("Done.")); fflush(stdout); rt_stop(sortEventTime); /* ------------------------------ */ /*create the status structure */ StatusList *status_struct = create_status_struct(); /*Put cells that are initially on the sweepline into status structure */ Rtimer sweepTime; StatusNode sn; rt_start(sweepTime); for (dimensionType i = vp->col + 1; i < hd->ncols; i++) { AEvent e; double ax, ay; sn.col = i; sn.row = vp->row; e.col = i; e.row = vp->row; e.elev[0] = data[0][i]; e.elev[1] = data[1][i]; e.elev[2] = data[2][i]; if (!is_nodata(visgrid->grid->hd, data[1][i]) && !is_point_outside_max_dist(*vp, *hd, sn.row, sn.col, viewOptions.maxDist)) { /*calculate Distance to VP and Gradient, store them into sn */ /* need either 3 elevation values or * 3 gradients calculated from 3 elevation values */ /* need also 3 angles */ e.eventType = ENTERING_EVENT; calculate_event_position(e, vp->row, vp->col, &ay, &ax); sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row); calculate_event_gradient(&sn, 0, ay, ax, e.elev[0], vp, *hd); e.eventType = CENTER_EVENT; calculate_event_position(e, vp->row, vp->col, &ay, &ax); sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row); calculate_dist_n_gradient(&sn, e.elev[1], vp, *hd); e.eventType = EXITING_EVENT; calculate_event_position(e, vp->row, vp->col, &ay, &ax); sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row); calculate_event_gradient(&sn, 2, ay, ax, e.elev[2], vp, *hd); assert(sn.angle[1] == 0); if (sn.angle[0] > sn.angle[1]) sn.angle[0] -= 2 * M_PI; G_debug(2, "inserting: "); print_statusnode(sn); /*insert sn into the status structure */ insert_into_status_struct(sn, status_struct); } } G_free(data[0]); G_free(data); /* ------------------------------ */ /*sweep the event list */ long nvis = 0; /*number of visible cells */ AEvent *e; G_important_message(_("Computing visibility...")); G_percent(0, 100, 2); for (size_t i = 0; i < nevents; i++) { int perc = (int)(1000000 * i / nevents); if (perc > 0 && perc < 1000000) G_percent(perc, 1000000, 1); /*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, e->elev[1] + vp->target_offset, vp, *hd); G_debug(3, "event: "); print_event(*e, 3); G_debug(3, "sn.dist=%f, sn.gradient=%f", sn.dist2vp, sn.gradient[1]); switch (e->eventType) { case ENTERING_EVENT: double ax, ay; /*insert node into structure */ G_debug(3, "..ENTER-EVENT: insert"); /* need either 3 elevation values or * 3 gradients calculated from 3 elevation values */ /* need also 3 angles */ calculate_event_position(*e, vp->row, vp->col, &ay, &ax); //sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row); sn.angle[0] = e->angle; calculate_event_gradient(&sn, 0, ay, ax, e->elev[0], vp, *hd); e->eventType = CENTER_EVENT; calculate_event_position(*e, vp->row, vp->col, &ay, &ax); sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row); calculate_dist_n_gradient(&sn, e->elev[1], vp, *hd); e->eventType = EXITING_EVENT; calculate_event_position(*e, vp->row, vp->col, &ay, &ax); sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row); calculate_event_gradient(&sn, 2, ay, ax, e->elev[2], vp, *hd); e->eventType = ENTERING_EVENT; if (e->angle < M_PI) { if (sn.angle[0] > sn.angle[1]) sn.angle[0] -= 2 * M_PI; } else { if (sn.angle[0] > sn.angle[1]) { sn.angle[1] += 2 * M_PI; sn.angle[2] += 2 * M_PI; } } insert_into_status_struct(sn, status_struct); break; case EXITING_EVENT: /*delete node out of status structure */ G_debug(3, "..EXIT-EVENT: delete"); /* need only distance */ delete_from_status_struct(status_struct, sn.dist2vp); break; case CENTER_EVENT: G_debug(3, "..QUERY-EVENT: query"); /*calculate visibility */ double max; /* consider current angle and gradient */ max = find_max_gradient_in_status_struct(status_struct, sn.dist2vp, e->angle, sn.gradient[1]); /*the point is visible: store its vertical angle */ if (max <= sn.gradient[1]) { float vert_angle = get_vertical_angle(*vp, sn, e->elev[1] + vp->target_offset, viewOptions.doCurv); add_result_to_inmem_visibilitygrid(visgrid, sn.row, sn.col, vert_angle); assert(vert_angle >= 0); /* when you write the visibility grid you assume that visible values are positive */ nvis++; } //else { /* cell is invisible */ /* the visibility grid is initialized all invisible */ //visgrid->grid->grid_data[sn.row][sn.col] = INVISIBLE; //} break; } } rt_stop(sweepTime); G_percent(1, 1, 1); G_verbose_message(_("Sweeping done.")); G_verbose_message(_("Total cells %ld, visible cells %ld (%.1f percent)."), (long)visgrid->grid->hd->nrows * visgrid->grid->hd->ncols, nvis, (float)((float)nvis * 100 / (float)(visgrid->grid->hd->nrows * visgrid->grid->hd->ncols))); print_viewshed_timings(initEventTime, sortEventTime, sweepTime); /*cleanup */ G_free(eventList); return visgrid; } /*/////////////////////////////////////////////////////////// ------------------------------------------------------------ run Viewshed's algorithm on the grid stored in the given file, and with the given viewpoint. Create a visibility grid and return it. It runs in external memory, i.e. the input grid and the outpt grid are stored as streams */ IOVisibilityGrid *viewshed_external(char *inputfname, GridHeader * hd, Viewpoint * vp, ViewOptions viewOptions) { assert(inputfname && hd && vp); G_message(_("Start sweeping.")); /* ------------------------------ */ /*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; surface_type **data; rt_start(initEventTime); eventList = init_event_list(inputfname, vp, hd, viewOptions, &data, visgrid); 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 */ G_verbose_message(_("Sorting events...")); 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; G_message(_("Initialize sweepline...")); rt_start(sweepTime); for (dimensionType i = vp->col + 1; i < hd->ncols; i++) { AEvent e; double ax, ay; G_percent(i, hd->ncols, 2); sn.col = i; sn.row = vp->row; e.col = i; e.row = vp->row; e.elev[0] = data[0][i]; e.elev[1] = data[1][i]; e.elev[2] = data[2][i]; if (!is_nodata(visgrid->hd, data[1][i]) && !is_point_outside_max_dist(*vp, *hd, sn.row, sn.col, viewOptions.maxDist)) { /*calculate Distance to VP and Gradient, store them into sn */ /* need either 3 elevation values or * 3 gradients calculated from 3 elevation values */ /* need also 3 angles */ e.eventType = ENTERING_EVENT; calculate_event_position(e, vp->row, vp->col, &ay, &ax); sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row); calculate_event_gradient(&sn, 0, ay, ax, e.elev[0], vp, *hd); e.eventType = CENTER_EVENT; calculate_event_position(e, vp->row, vp->col, &ay, &ax); sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row); calculate_dist_n_gradient(&sn, e.elev[1], vp, *hd); e.eventType = EXITING_EVENT; calculate_event_position(e, vp->row, vp->col, &ay, &ax); sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row); calculate_event_gradient(&sn, 2, ay, ax, e.elev[2], vp, *hd); assert(sn.angle[1] == 0); if (sn.angle[0] > sn.angle[1]) sn.angle[0] -= 2 * M_PI; G_debug(3, "inserting: "); print_statusnode(sn); /*insert sn into the status structure */ insert_into_status_struct(sn, status_struct); } } G_percent(hd->ncols, hd->ncols, 2); G_free(data[0]); G_free(data); /* ------------------------------ */ /*sweep the event list */ long nvis = 0; /*number of visible cells */ VisCell viscell; AEvent *e; AMI_err ae; off_t nbEvents = eventList->stream_len(); /*printf("nbEvents = %ld\n", (long) nbEvents); */ G_message(_("Determine visibility...")); G_percent(0, 100, 2); for (off_t i = 0; i < nbEvents; i++) { int perc = (int)(1000000 * i / nbEvents); if (perc > 0) G_percent(perc, 1000000, 1); /*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, e->elev[1] + vp->target_offset, vp, *hd); G_debug(3, "next event: "); print_statusnode(sn); switch (e->eventType) { case ENTERING_EVENT: double ax, ay; /*insert node into structure */ /* need either 3 elevation values or * 3 gradients calculated from 3 elevation values */ /* need also 3 angles */ calculate_event_position(*e, vp->row, vp->col, &ay, &ax); //sn.angle[0] = calculate_angle(ax, ay, vp->col, vp->row); sn.angle[0] = e->angle; calculate_event_gradient(&sn, 0, ay, ax, e->elev[0], vp, *hd); e->eventType = CENTER_EVENT; calculate_event_position(*e, vp->row, vp->col, &ay, &ax); sn.angle[1] = calculate_angle(ax, ay, vp->col, vp->row); calculate_dist_n_gradient(&sn, e->elev[1], vp, *hd); e->eventType = EXITING_EVENT; calculate_event_position(*e, vp->row, vp->col, &ay, &ax); sn.angle[2] = calculate_angle(ax, ay, vp->col, vp->row); calculate_event_gradient(&sn, 2, ay, ax, e->elev[2], vp, *hd); e->eventType = ENTERING_EVENT; if (e->angle < M_PI) { if (sn.angle[0] > sn.angle[1]) sn.angle[0] -= 2 * M_PI; } else { if (sn.angle[0] > sn.angle[1]) { sn.angle[1] += 2 * M_PI; sn.angle[2] += 2 * M_PI; } } G_debug(3, "..ENTER-EVENT: insert"); insert_into_status_struct(sn, status_struct); break; case EXITING_EVENT: /*delete node out of status structure */ G_debug(3, "..EXIT-EVENT: delete"); delete_from_status_struct(status_struct, sn.dist2vp); break; case CENTER_EVENT: G_debug(3, "..QUERY-EVENT: query"); /*calculate visibility */ viscell.row = sn.row; viscell.col = sn.col; double max; max = find_max_gradient_in_status_struct(status_struct, sn.dist2vp, e->angle, sn.gradient[1]); /*the point is visible */ if (max <= sn.gradient[1]) { viscell.angle = get_vertical_angle(*vp, sn, e->elev[1] + vp->target_offset, viewOptions.doCurv); assert(viscell.angle >= 0); /* viscell.vis = VISIBLE; */ add_result_to_io_visibilitygrid(visgrid, &viscell); nvis++; } else { /* else the cell is invisible; we do not write it to the visibility stream because we only record visible cells, and nodata cells; */ /* viscell.vis = INVISIBLE; */ /* add_result_to_io_visibilitygrid(visgrid, &viscell); */ } break; } } /* for each event */ rt_stop(sweepTime); G_percent(1, 1, 1); G_message(_("Sweeping done.")); G_verbose_message(_("Total cells %ld, visible cells %ld (%.1f percent)."), (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; }