Feature Request: Add Line of Sight function

[Ryanf55]

Expected behavior and actual behavior.

I would like to be able to call a C++ function in GDAL to compute the Line of Sight visibility between two points over a DEM raster band. Currently, there only seems to be viewshed support, which is a lot of overhead if I just need to check line of sight between two points.

There was a request for an algorithm on the forums; it never made it into a supported feature: https://gis.stackexchange.com/questions/228726/algorithm-or-gdal-module-to-calculate-los

Desired behavior

GDAL could add a new function like so:

struct Loc {
    double x;
    double y;
    double z;
}

static bool IsLineOfSightVisible(const GDALRasterBand& band, const Loc& locA, const Loc& locB); 

If LocA and LocB, and the geodetic line between them are all above terrain, it returns true. Otherwise false.

If either of the locations are outside the raster band, or the raster band reports no data at any of the locations along the line, then it would also return false.

From my perspective, it would be even better if there were an API that I could pass in the locations as LLH and have GDAL deal with the projections/conversions internally.

Papers

https://www.researchgate.net/publication/2411280_Efficient_Line-of-Sight_Algorithms_for_Real_Terrain_Data

This recommends use of Bresenham's line of sight algorithm in 3D.

Operating system

Ubuntu 22.04

GDAL version and provenance

master

Other notes

Pending the support of the GDAL maintainers to add this feature, and they agree on the algorithm I found, I am happy to try to implement this and contribute it.

[rouault]

I'd recommend a C compatible signature like

bool CPL_DLL GDALIsLineOfSightVisible(GDALRasterBandH, double xA, double yA, double zA, double xB, double yB, double zB, char** papszOptions);

implemented in a alg/los.cpp file

It would be good if we could be consistent of the results of GDALViewshedGenerate(), in particular the curvature coefficient (https://gdal.org/programs/gdal_viewshed.html#cmdoption-gdal_viewshed-cc), although that could be a later added option.

CC'ing @szekerest the author of gdal_viewshed

[Ryanf55]

Awesome. Thanks for the tips. I appreciate any further recommendations.

The use case driving this is for aerial drone path planning. The line of sight check would be used to check a planned path does not request the drone to fly through terrain and crash. Thus, accounting for refraction of light is not strictly necessary for the feature to still be useful. Current open-source implementations of path planners over terrain use a more reactive approach to adjust the path on the fly, but this is less robust than what could be done with a line of sight algorithm.

[sa2329]

One approach to checking line-of sight for two points is to translate / warp the elevation data to extract elevations only between the two points of interest, then compute the viewshed on that. GDAL will not calculate the viewshed for points outside the elevation data, so that can reduce the amount of computation needed substantially.

[Ryanf55]

Would that be more efficient you think than bresenham?

[sa2329]

I don't know about efficiency, particularly if you're doing a large number of line of sight computations. I'm not familiar with bresenham. However, it's fairly straightforward to implement and test. The warp operation will only compute points for the extracted elevations, so it should run fairly fast.

[Ryanf55]

I'd recommend a C compatible signature like

bool CPL_DLL GDALIsLineOfSightVisible(GDALRasterBandH, double xA, double yA, double zA, double xB, double yB, double zB, char** papszOptions);

implemented in a alg/los.cpp file

It would be good if we could be consistent of the results of GDALViewshedGenerate(), in particular the curvature coefficient (https://gdal.org/programs/gdal_viewshed.html#cmdoption-gdal_viewshed-cc), although that could be a later added option.

CC'ing @szekerest the author of gdal_viewshed

I started implementation of this. I was curious why you suggest GDALIsLineOfSightVisible take in location coordinates as a double, but the raster coordinates for RasterIO are int. Do you intend that people would use projection-space for their LOS points, and expect the function to convert to raster-space coordinates inside? Or, should that be a user-selectable behavior depending on the value of papszOptions?

If the first implementation could be restricted to raster-coordinate LOS computation, the PR would be more straightforward.

Thanks!

[rouault]

I was curious why you suggest GDALIsLineOfSightVisible take in location coordinates as a double, but the raster coordinates for RasterIO are int

int for x and y would be indeed more in line with RasterIO() API (unless sub-pixel positioning would be needed). I was mostly reacting to introducing a struct Loc which would be a small complication for little benefit, in particular for bindings to other languages.

[Ryanf55]

Hello. I've pushed a draft PR to garner some feedback. It meets the minimum requirements of what was discussed here and comes with some tests. Please take a look. It should be open to extension without breaking ABI for other algorithm approaches than bresenham.

[Ryanf55]

Awesome, glad to get this in.

Future plans (contributions welcome)

• Add reporting of where the LOS check fails (may break ABI - so I will try to do this before release).
• Support special case optimization for exact vertical and exact horizontal checks
• Add benchmarking and optimize performance of LOS
• GPU accelerate it

[rouault]

• GPU accelerate it

Tough subject. I suspect benchmarking will show that the bottleneck is in the pixel acquisition part, not the computation. But GPU acceleration is a subject where we might welcome expertise.

[Ryanf55]

• GPU accelerate it

Tough subject. I suspect benchmarking will show that the bottleneck is in the pixel acquisition part, not the computation. But GPU acceleration is a subject where we might welcome expertise.

Yea. I don't have much either, but I know there's a lot of work on GPU ray casting which what essentially what LOS checking is doing here. I won't personally be spending time on it until my team shows that we need to make it faster.