Add special cases for LOS checks

* This adds optimizations for vertical/horizontal, where bresenham
  shouldn't be used.
* This avoids a bunch of sqrt options in GetZValueFromXY and should
  speed things up

Signed-off-by: Ryan Friedman <[email protected]>

alg/los.cpp

@@ -205,19 +205,113 @@ bool GDALIsLineOfSightVisible(const GDALRasterBandH hBand, const int xA,
         return true;
     }
 
-    // TODO if both X's or Y's are the same, it could be optimized for vertical/horizontal lines.
+    // Lambda for Linear interpolate like C++20 std::lerp.
+    auto lerp = [](const double a, const double b, const double t)
+    { return a + t * (b - a); };
+
+    // Lambda for getting Z test height given y input along the LOS line.
+    // Only to be used for vertical line checks.
+    auto GetZValueFromY = [&](const int y) -> double
+    {
+        // A ratio of 0.0 corresponds to being at yA.
+        const auto ratio =
+            static_cast<double>(y - yA) / static_cast<double>(yB - yA);
+        return lerp(zA, zB, ratio);
+    };
+
+    // Lambda for getting Z test height given x input along the LOS line.
+    // Only to be used for horizontal line checks.
+    auto GetZValueFromX = [&](const int x) -> double
+    {
+        // A ratio of 0.0 corresponds to being at xA.
+        const auto ratio =
+            static_cast<double>(x - xA) / static_cast<double>(xB - xA);
+        return lerp(zA, zB, ratio);
+    };
+
+    // Lambda for checking path safety of a vertical line.
+    // Returns true if the path has clear LOS.
+    auto CheckVerticalLine = [&]() -> bool
+    {
+        if (xA < xB)
+        {
+            for (int x = xA; x <= xB; ++x)
+            {
+                const auto zTest = GetZValueFromX(x);
+                if (!IsAboveTerrain(hBand, x, yA, zTest))
+                {
+                    return false;
+                }
+            }
+            return true;
+        }
+        if (xA > xB)
+        {
+            for (int x = xA; x >= xB; --x)
+            {
+                const auto zTest = GetZValueFromX(x);
+                if (!IsAboveTerrain(hBand, x, yA, zTest))
+                {
+                    return false;
+                }
+            }
+            return true;
+        }
+
+        // If we get here, it's a coding error.
+        return false;
+    };
+
+    // Lambda for checking path safety of a horizontal line.
+    // Returns true if the path has clear LOS.
+    auto CheckHorizontalLine = [&]() -> bool
+    {
+        if (yA < yB)
+        {
+            for (int y = yA; y <= yB; ++y)
+            {
+                const auto zTest = GetZValueFromY(y);
+                if (!IsAboveTerrain(hBand, xA, y, zTest))
+                {
+                    return false;
+                }
+            }
+            return true;
+        }
+        if (yA > yB)
+        {
+            for (int y = yA; y >= yB; --y)
+            {
+                const auto zTest = GetZValueFromX(y);
+                if (!IsAboveTerrain(hBand, xA, y, zTest))
+                {
+                    return false;
+                }
+            }
+            return true;
+        }
+
+        // If we get here, it's a coding error.
+        return false;
+    };
+
+    // Handle special cases if it's a vertical or horizontal line (don't use bresenham).
+    if (xA == xB)
+    {
+        return CheckVerticalLine();
+    }
+    if (yA == yB)
+    {
+        return CheckHorizontalLine();
+    }
 
     // Use an interpolated Z height with 2D bresenham for the remaining cases.
 
     // Lambda for computing the square of a number
     auto SQUARE = [](const double d) -> double { return d * d; };
 
-    // Lambda for Linear interpolate like C++20 std::lerp.
-    auto lerp = [](const double a, const double b, const double t)
-    { return a + t * (b - a); };
-
     // Lambda for getting Z test height given x-y input along the bresenham line.
-    auto GetZValue = [&](const int x, const int y) -> double
+    auto GetZValueFromXY = [&](const int x, const int y) -> double
     {
         const auto rNum = SQUARE(static_cast<double>(x - xA)) +
                           SQUARE(static_cast<double>(y - yA));
@@ -235,7 +329,7 @@ bool GDALIsLineOfSightVisible(const GDALRasterBandH hBand, const int xA,
     // Lambda to get elevation at a bresenham-computed location.
     auto OnBresenhamPoint = [&](const int x, const int y) -> bool
     {
-        const auto z = GetZValue(x, y);
+        const auto z = GetZValueFromXY(x, y);
         return IsAboveTerrain(hBand, x, y, z);
     };
 

autotest/cpp/test_alg.cpp

@@ -429,13 +429,25 @@ TEST_F(test_alg, GDALIsLineOfSightVisible_through_mountain)
     EXPECT_FALSE(
         GDALIsLineOfSightVisible(pBand, 0, 120, 203, 120, 0, 247, nullptr));
 
-    // Vertical line test with hill between two points.
+    // Vertical line tests with hill between two points, in both directions.
     EXPECT_FALSE(
         GDALIsLineOfSightVisible(pBand, 83, 111, 154, 83, 117, 198, nullptr));
+    EXPECT_FALSE(
+        GDALIsLineOfSightVisible(pBand, 83, 117, 198, 83, 111, 154, nullptr));
+    EXPECT_TRUE(
+        GDALIsLineOfSightVisible(pBand, 83, 111, 460, 83, 117, 460, nullptr));
+    EXPECT_TRUE(
+        GDALIsLineOfSightVisible(pBand, 83, 117, 460, 83, 111, 460, nullptr));
 
-    // Horizonal line test with hill between two points.
+    // Horizonal line tests with hill between two points, in both directions.
     EXPECT_FALSE(
         GDALIsLineOfSightVisible(pBand, 75, 115, 192, 89, 115, 191, nullptr));
+    EXPECT_FALSE(
+        GDALIsLineOfSightVisible(pBand, 89, 115, 191, 75, 115, 192, nullptr));
+    EXPECT_TRUE(
+        GDALIsLineOfSightVisible(pBand, 75, 115, 460, 89, 115, 460, nullptr));
+    EXPECT_TRUE(
+        GDALIsLineOfSightVisible(pBand, 89, 115, 460, 75, 115, 460, nullptr));
 }
 
 }  // namespace