LCOV - code coverage report
Current view: top level - alg - los.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 132 132 100.0 %
Date: 2024-11-21 22:18:42 Functions: 13 13 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Line of Sight
       4             :  * Purpose:  Core algorithm implementation for line of sight algorithms.
       5             :  * Author:   Ryan Friedman, ryanfriedman5410+gdal@gmail.com
       6             :  *
       7             :  ******************************************************************************
       8             :  *
       9             :  * SPDX-License-Identifier: MIT
      10             :  ****************************************************************************/
      11             : 
      12             : #include <functional>
      13             : #include <cmath>
      14             : 
      15             : #include "cpl_port.h"
      16             : #include "gdal_alg.h"
      17             : 
      18             : // There's a plethora of bresenham implementations, all questionable production quality.
      19             : // Bresenham optimizes for integer math, which makes sense for raster datasets in 2D.
      20             : // For 3D, a 3D bresenham could be used if the altitude is also integer resolution.
      21             : // 2D:
      22             : // https://codereview.stackexchange.com/questions/77460/bresenhams-line-algorithm-optimization
      23             : // https://gist.github.com/ssavi-ict/092501c69e2ffec65e96a8865470ad2f
      24             : // https://blog.demofox.org/2015/01/17/bresenhams-drawing-algorithms/
      25             : // https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm
      26             : // https://www.cs.helsinki.fi/group/goa/mallinnus/lines/bresenh.html
      27             : // https://stackoverflow.com/questions/10060046/drawing-lines-with-bresenhams-line-algorithm
      28             : // http://www.edepot.com/linebresenham.html
      29             : // 3D:
      30             : // https://gist.github.com/yamamushi/5823518
      31             : 
      32             : // Run bresenham terrain checking from (x1, y1) to (x2, y2).
      33             : // The callback is run at every point along the line,
      34             : // which should return True if the point is above terrain.
      35             : // Bresenham2D will return true if all points have LOS between the start and end.
      36             : static bool
      37           8 : Bresenham2D(const int x1, const int y1, const int x2, const int y2,
      38             :             std::function<auto(const int, const int)->bool> OnBresenhamPoint)
      39             : {
      40           8 :     bool isAboveTerrain = true;
      41             :     int dx, dy;
      42             :     int incx, incy;
      43             : 
      44           8 :     if (x2 >= x1)
      45             :     {
      46           7 :         dx = x2 - x1;
      47           7 :         incx = 1;
      48             :     }
      49             :     else
      50             :     {
      51           1 :         dx = x1 - x2;
      52           1 :         incx = -1;
      53             :     }
      54             : 
      55           8 :     if (y2 >= y1)
      56             :     {
      57           6 :         dy = y2 - y1;
      58           6 :         incy = 1;
      59             :     }
      60             :     else
      61             :     {
      62           2 :         dy = y1 - y2;
      63           2 :         incy = -1;
      64             :     }
      65             : 
      66           8 :     auto x = x1;
      67           8 :     auto y = y1;
      68             :     int balance;
      69             : 
      70           8 :     if (dx >= dy)
      71             :     {
      72           7 :         dy <<= 1;
      73           7 :         balance = dy - dx;
      74           7 :         dx *= 2;
      75             : 
      76         254 :         while (x != x2 && isAboveTerrain)
      77             :         {
      78         247 :             isAboveTerrain &= OnBresenhamPoint(x, y);
      79         247 :             if (balance >= 0)
      80             :             {
      81         247 :                 y += incy;
      82         247 :                 balance -= dx;
      83             :             }
      84         247 :             balance += dy;
      85         247 :             x += incx;
      86             :         }
      87           7 :         isAboveTerrain &= OnBresenhamPoint(x, y);
      88             :     }
      89             :     else
      90             :     {
      91           1 :         dx *= 2;
      92           1 :         balance = dx - dy;
      93           1 :         dy *= 2;
      94             : 
      95           7 :         while (y != y2 && isAboveTerrain)
      96             :         {
      97           6 :             isAboveTerrain &= OnBresenhamPoint(x, y);
      98           6 :             if (balance >= 0)
      99             :             {
     100           4 :                 x += incx;
     101           4 :                 balance -= dy;
     102             :             }
     103           6 :             balance += dx;
     104           6 :             y += incy;
     105             :         }
     106           1 :         isAboveTerrain &= OnBresenhamPoint(x, y);
     107             :     }
     108           8 :     return isAboveTerrain;
     109             : }
     110             : 
     111             : // Get the elevation of a single point.
     112         637 : static bool GetElevation(const GDALRasterBandH hBand, const int x, const int y,
     113             :                          double &val)
     114             : {
     115             :     /// @todo GDALCachedPixelAccessor may give increased performance.
     116         637 :     return GDALRasterIO(hBand, GF_Read, x, y, 1, 1, &val, 1, 1, GDT_Float64, 0,
     117         637 :                         0) == CE_None;
     118             : }
     119             : 
     120             : // Check a single location is above terrain.
     121         637 : static bool IsAboveTerrain(const GDALRasterBandH hBand, const int x,
     122             :                            const int y, const double z)
     123             : {
     124             :     double terrainHeight;
     125         637 :     if (GetElevation(hBand, x, y, terrainHeight))
     126             :     {
     127         636 :         return z > terrainHeight;
     128             :     }
     129             :     else
     130             :     {
     131           1 :         return false;
     132             :     }
     133             : }
     134             : 
     135             : /************************************************************************/
     136             : /*                        GDALIsLineOfSightVisible()                    */
     137             : /************************************************************************/
     138             : 
     139             : /**
     140             :  * Check Line of Sight between two points.
     141             :  * Both input coordinates must be within the raster coordinate bounds.
     142             :  *
     143             :  * This algorithm will check line of sight using a Bresenham algorithm.
     144             :  * https://www.researchgate.net/publication/2411280_Efficient_Line-of-Sight_Algorithms_for_Real_Terrain_Data
     145             :  * Line of sight is computed in raster coordinate space, and thus may not be appropriate.
     146             :  * For example, datasets referenced against geographic coordinate at high latitudes may have issues.
     147             :  *
     148             :  * @param hBand The band to read the DEM data from. This must NOT be null.
     149             :  *
     150             :  * @param xA The X location (raster column) of the first point to check on the raster.
     151             :  *
     152             :  * @param yA The Y location (raster row) of the first point to check on the raster.
     153             :  *
     154             :  * @param zA The Z location (height) of the first point to check.
     155             :  *
     156             :  * @param xB The X location (raster column) of the second point to check on the raster.
     157             :  *
     158             :  * @param yB The Y location (raster row) of the second point to check on the raster.
     159             :  *
     160             :  * @param zB The Z location (height) of the second point to check.
     161             :  *
     162             :  * @param[out] pnxTerrainIntersection The X location where the LOS line
     163             :  *             intersects with terrain, or nullptr if it does not intersect
     164             :  *             terrain.
     165             :  *
     166             :  * @param[out] pnyTerrainIntersection The Y location where the LOS line
     167             :  *             intersects with terrain, or nullptr if it does not intersect
     168             :  *             terrain.
     169             :  *
     170             :  * @param papszOptions Options for the line of sight algorithm (currently ignored).
     171             :  *
     172             :  * @return True if the two points are within Line of Sight.
     173             :  *
     174             :  * @since GDAL 3.9
     175             :  */
     176             : 
     177          36 : bool GDALIsLineOfSightVisible(const GDALRasterBandH hBand, const int xA,
     178             :                               const int yA, const double zA, const int xB,
     179             :                               const int yB, const double zB,
     180             :                               int *pnxTerrainIntersection,
     181             :                               int *pnyTerrainIntersection,
     182             :                               CPL_UNUSED CSLConstList papszOptions)
     183             : {
     184          36 :     VALIDATE_POINTER1(hBand, "GDALIsLineOfSightVisible", false);
     185             : 
     186             :     // A lambda to set the X-Y intersection if it's not null
     187          24 :     auto SetXYIntersection = [&](const int x, const int y)
     188             :     {
     189          24 :         if (pnxTerrainIntersection != nullptr)
     190             :         {
     191           9 :             *pnxTerrainIntersection = x;
     192             :         }
     193          24 :         if (pnyTerrainIntersection != nullptr)
     194             :         {
     195           9 :             *pnyTerrainIntersection = y;
     196             :         }
     197          60 :     };
     198             : 
     199          36 :     if (pnxTerrainIntersection)
     200          12 :         *pnxTerrainIntersection = -1;
     201             : 
     202          36 :     if (pnyTerrainIntersection)
     203          12 :         *pnyTerrainIntersection = -1;
     204             : 
     205             :     // Perform a preliminary check of the start and end points.
     206          36 :     if (!IsAboveTerrain(hBand, xA, yA, zA))
     207             :     {
     208          11 :         SetXYIntersection(xA, yA);
     209          11 :         return false;
     210             :     }
     211          25 :     if (!IsAboveTerrain(hBand, xB, yB, zB))
     212             :     {
     213           5 :         SetXYIntersection(xB, yB);
     214           5 :         return false;
     215             :     }
     216             : 
     217             :     // If both X and Y are the same, no further checks are needed.
     218          20 :     if (xA == xB && yA == yB)
     219             :     {
     220           3 :         return true;
     221             :     }
     222             : 
     223             :     // Lambda for Linear interpolate like C++20 std::lerp.
     224         315 :     auto lerp = [](const double a, const double b, const double t)
     225         315 :     { return a + t * (b - a); };
     226             : 
     227             :     // Lambda for getting Z test height given y input along the LOS line.
     228             :     // Only to be used for vertical line checks.
     229          18 :     auto GetZValueFromY = [&](const int y) -> double
     230             :     {
     231             :         // A ratio of 0.0 corresponds to being at yA.
     232          18 :         const auto ratio =
     233          18 :             static_cast<double>(y - yA) / static_cast<double>(yB - yA);
     234          18 :         return lerp(zA, zB, ratio);
     235          17 :     };
     236             : 
     237             :     // Lambda for getting Z test height given x input along the LOS line.
     238             :     // Only to be used for horizontal line checks.
     239          36 :     auto GetZValueFromX = [&](const int x) -> double
     240             :     {
     241             :         // A ratio of 0.0 corresponds to being at xA.
     242          36 :         const auto ratio =
     243          36 :             static_cast<double>(x - xA) / static_cast<double>(xB - xA);
     244          36 :         return lerp(zA, zB, ratio);
     245          17 :     };
     246             : 
     247             :     // Lambda for checking path safety of a vertical line.
     248             :     // Returns true if the path has clear LOS.
     249           4 :     auto CheckVerticalLine = [&]() -> bool
     250             :     {
     251           4 :         CPLAssert(xA == xB);
     252           4 :         CPLAssert(yA != yB);
     253             : 
     254           4 :         if (yA < yB)
     255             :         {
     256          10 :             for (int y = yA; y <= yB; ++y)
     257             :             {
     258           9 :                 const auto zTest = GetZValueFromY(y);
     259           9 :                 if (!IsAboveTerrain(hBand, xA, y, zTest))
     260             :                 {
     261           1 :                     SetXYIntersection(xA, y);
     262           1 :                     return false;
     263             :                 }
     264             :             }
     265           1 :             return true;
     266             :         }
     267             :         else
     268             :         {
     269          10 :             for (int y = yA; y >= yB; --y)
     270             :             {
     271           9 :                 const auto zTest = GetZValueFromY(y);
     272           9 :                 if (!IsAboveTerrain(hBand, xA, y, zTest))
     273             :                 {
     274           1 :                     SetXYIntersection(xA, y);
     275           1 :                     return false;
     276             :                 }
     277             :             }
     278           1 :             return true;
     279             :         }
     280          17 :     };
     281             : 
     282             :     // Lambda for checking path safety of a horizontal line.
     283             :     // Returns true if the path has clear LOS.
     284           5 :     auto CheckHorizontalLine = [&]() -> bool
     285             :     {
     286           5 :         CPLAssert(yA == yB);
     287           5 :         CPLAssert(xA != xB);
     288             : 
     289           5 :         if (xA < xB)
     290             :         {
     291          21 :             for (int x = xA; x <= xB; ++x)
     292             :             {
     293          19 :                 const auto zTest = GetZValueFromX(x);
     294          19 :                 if (!IsAboveTerrain(hBand, x, yA, zTest))
     295             :                 {
     296           1 :                     SetXYIntersection(x, yA);
     297           1 :                     return false;
     298             :                 }
     299             :             }
     300           2 :             return true;
     301             :         }
     302             :         else
     303             :         {
     304          18 :             for (int x = xA; x >= xB; --x)
     305             :             {
     306          17 :                 const auto zTest = GetZValueFromX(x);
     307          17 :                 if (!IsAboveTerrain(hBand, x, yA, zTest))
     308             :                 {
     309           1 :                     SetXYIntersection(x, yA);
     310           1 :                     return false;
     311             :                 }
     312             :             }
     313           1 :             return true;
     314             :         }
     315          17 :     };
     316             : 
     317             :     // Handle special cases if it's a vertical or horizontal line (don't use bresenham).
     318          17 :     if (xA == xB)
     319             :     {
     320           4 :         return CheckVerticalLine();
     321             :     }
     322          13 :     if (yA == yB)
     323             :     {
     324           5 :         return CheckHorizontalLine();
     325             :     }
     326             : 
     327             :     // Use an interpolated Z height with 2D bresenham for the remaining cases.
     328             : 
     329             :     // Lambda for computing the square of a number
     330        1044 :     auto SQUARE = [](const double d) -> double { return d * d; };
     331             : 
     332             :     // Lambda for getting Z test height given x-y input along the bresenham line.
     333         261 :     auto GetZValueFromXY = [&](const int x, const int y) -> double
     334             :     {
     335         261 :         const auto rNum = SQUARE(static_cast<double>(x - xA)) +
     336         261 :                           SQUARE(static_cast<double>(y - yA));
     337         261 :         const auto rDenom = SQUARE(static_cast<double>(xB - xA)) +
     338         261 :                             SQUARE(static_cast<double>(yB - yA));
     339             :         /// @todo In order to reduce CPU cost and avoid a sqrt operation, consider
     340             :         /// the approach to just the ratio along x or y depending on whether
     341             :         /// the line is steep or shallow.
     342             :         /// See https://github.com/OSGeo/gdal/pull/9506#discussion_r1532459689.
     343             :         const double ratio =
     344         261 :             sqrt(static_cast<double>(rNum) / static_cast<double>(rDenom));
     345         261 :         return lerp(zA, zB, ratio);
     346           8 :     };
     347             : 
     348             :     // Lambda to get elevation at a bresenham-computed location.
     349         261 :     auto OnBresenhamPoint = [&](const int x, const int y) -> bool
     350             :     {
     351         261 :         const auto z = GetZValueFromXY(x, y);
     352         261 :         const auto isAbove = IsAboveTerrain(hBand, x, y, z);
     353         261 :         if (!isAbove)
     354             :         {
     355           4 :             SetXYIntersection(x, y);
     356             :         }
     357         261 :         return IsAboveTerrain(hBand, x, y, z);
     358           8 :     };
     359             : 
     360           8 :     return Bresenham2D(xA, yA, xB, yB, OnBresenhamPoint);
     361             : }

Generated by: LCOV version 1.14