GDALCircleLoS for checking line of sight of a circular path over a height map

[Ryanf55]

Feature description

Add an algorithm similar to GDALIsLineOfSightVisible, this time for circles, perhaps called GDALIsCircleLineOfSightVisible

Proposed name and function signature

bool CPL_DLL GDALIsCircleLineOfSightVisible(const GDALRasterBandH, const int x,
  const int y, const double z,
  const int radius,
  CSLConstList papszOptions);

Given a prescribed raster X-Y coordinate, a heightmap raster band, and a "radius" (number of pixels), and an elevation, return true if the path traced by the circle is safe above terrain.

Additional context

This is used in aerial robotics to complete basic planners. Now that we have line of sight checks for straight line distances, this will provide the other primitive - fixed-altitude turns!

As before with LoS, I'm happy to do the implementation in GDAL, but would like to confirm this feature would be accepted and do a bit of planning before I start coding.

Further extensions

After circles, a potential next step would be 3D Dubins path with would leverage straight lines and circles primitives.

[rouault]

That seems fairly specialized to me, at least at the API level. What would be more acceptable (but more complicated to implement) would be a method to check if an arbitrary geometry can be seen from a point:

bool CPL_DLL GDALIsGeometryLineOfSightVisible(const GDALRasterBandH, const int x,
  const int y, const double z,
  OGRGeometryH hGeometry,
  CSLConstList papszOptions);

We could potentially restrict the geometry to be a curve, and potentially a circle (using a OGRCircularString) initially

[Ryanf55]

OGRGeometryH

Neat, sounds reasonable and less specific to circles. If you had any strong recommendations on implementation, let me know, otherwise I can probably get a minimal prototype running for some thoughts.

[rouault]

If you had any strong recommendations on implementation

not really. I'm wondering how you are going to discretize things. Bresenham? For circular strings, you can use the getLinearGeometry() method that will return a linestring approximation of it (discretizing by default every 4 degrees), but I guess there might be a Bresenham-like implementation for circular arcs that leads to a better result.

[Ryanf55]

I don't know anything off the top of my head but can do some research to find out.

[Ryanf55]

There is an equivalent to bresenham's line algorithm for circles - https://en.wikipedia.org/wiki/Midpoint_circle_algorithm

What I'm actually looking to implement is whether the path tracing a circle at a specific altitude is above or below terrain.

The application, combined with the existing LoS for straight lines, is for path planning over DEM for aerial vehicles that follow dubin paths (straight lines and constant-curvature turns). https://en.wikipedia.org/wiki/Dubins_path

[Ryanf55]

Jesko's method mentioned in wikipedia is nice and simple to discretize the circle. Here's a demo, adapted from here.

// sudo apt install libgnuplot-iostream-dev
// g++ main.cpp -o main -L/usr/local/lib/ -lboost_iostreams

#include <vector>
#include <utility>
#include <stdio.h>
#include "gnuplot-iostream.h"

void
circle( int mx, int my, int r )
{
    Gnuplot gp;
    gp << "set xrange [9:31]\nset yrange [19:41]\n";
    gp << "plot '-' title 'circle'\n";

    std::vector<std::pair<int, int>> pts;
    int t1 = r / 16;
    int t2 = 0;
    int x  = r;
    int y  = 0;
    while ( x >= y )
    {
        pts.emplace_back( mx + x, my + y );
        pts.emplace_back(  mx + x, my - y );
        pts.emplace_back(  mx - x, my + y );
        pts.emplace_back(  mx - x, my - y );
        pts.emplace_back(  mx + y, my + x );
        pts.emplace_back(  mx + y, my - x );
        pts.emplace_back(  mx - y, my + x );
        pts.emplace_back(  mx - y, my - x );
        y  = y + 1;
        t1 = t1 + y;
        t2 = t1 - x;
        if ( t2 >= 0 )
        {
            t1 = t2;
            x  = x - 1;
        }
    }

    for (auto const& p: pts) {
        std::cout << p.first << "," << p.second << "\n";
    }
     gp.send1d(pts);
}

int
main()
{
    circle( 20, 30, 9 );
    return 0;
}