Checkpoint.

Author: abellgithub

alg/viewshed/util.cpp

@@ -5,6 +5,7 @@
  ****************************************************************************/
 
 #include <array>
+#include <iostream>
 
 #include "gdal_priv.h"
 #include "util.h"
@@ -15,6 +16,81 @@ namespace gdal
 namespace viewshed
 {
 
+/// Normalize a masking angle. Change from clockwise with 0 north (up) to counterclockwise
+/// with 0 to the east (right) and change to radians.
+double normalizeAngle(double maskAngle)
+{
+    maskAngle = 90 - maskAngle;
+    if (maskAngle < 0)
+        maskAngle += 360;
+    return maskAngle * (M_PI / 180);
+}
+
+/// Compute the X intersect position on the line Y = y with a ray extending
+/// from (nX, nY) along `angle`.
+///ABELL doc args
+double horizontalIntersect(double angle, int nX, int nY, int y)
+{
+    double x = std::numeric_limits<double>::quiet_NaN();
+
+    if (nY == y)
+        x = nX;
+    else if (nY > y)
+    {
+        if (ARE_REAL_EQUAL(angle, M_PI / 2))
+            x = nX;
+        else if (angle > 0 && angle < M_PI)
+            x = nX + (nY - y) / std::tan(angle);
+    }
+    else  // nY < y
+    {
+        if (ARE_REAL_EQUAL(angle, 3 * M_PI / 2))
+            x = nX;
+        else if (angle > M_PI)
+            x = nX - (y - nY) / std::tan(angle);
+    }
+    return x;
+}
+
+/// Compute the Y intersect position on the line X = x with a ray extending
+/// from (nX, nY) along `angle`.
+///ABELL doc args
+double verticalIntersect(double angle, int nX, int nY, int x)
+{
+    double y = std::numeric_limits<double>::quiet_NaN();
+
+    if (nX == x)
+        y = nY;
+    else if (nX < x)
+    {
+        if (ARE_REAL_EQUAL(angle, 0))
+            y = nY;
+        else if (angle < M_PI / 2 || angle > M_PI * 3 / 2)
+            y = nY + (nX - x) * std::tan(angle);
+    }
+    else  // nX > x
+    {
+        if (ARE_REAL_EQUAL(angle, M_PI))
+            y = nY;
+        else if (angle > M_PI / 2 && angle < M_PI * 3 / 2)
+            y = nY - (x - nX) * std::tan(angle);
+    }
+    return y;
+}
+
+// Determine if ray is in the slice between two rays starting at `start` and
+// going clockwise to `end`.
+bool rayBetween(double start, double end, double test)
+{
+    // Our angles go counterclockwise, so swap start and end
+    std::swap(start, end);
+    if (start < end)
+        return (test > start && test < end);
+    else if (start > end)
+        return (test > start || test < end);
+    return false;
+}
+
 /// Get the band size
 ///
 /// @param  band Raster band

alg/viewshed/util.h

@@ -14,6 +14,10 @@ namespace gdal
 namespace viewshed
 {
 
+double normalizeAngle(double maskAngle);
+double horizontalIntersect(double angle, int nX, int nY, int y);
+double verticalIntersect(double angle, int nX, int nY, int x);
+bool rayBetween(double start, double end, double test);
 size_t bandSize(GDALRasterBand &band);
 
 DatasetPtr createOutputDataset(GDALRasterBand &srcBand, const Options &opts,

alg/viewshed/viewshed.cpp

@@ -12,6 +12,9 @@
  ****************************************************************************/
 
 #include <algorithm>
+//ABELL
+#include <iostream>
+#include <iomanip>
 
 #include "gdal_alg.h"
 #include "gdal_priv_templates.hpp"
@@ -216,6 +219,77 @@ bool getTransforms(GDALRasterBand &band, double *pFwdTransform,
     return true;
 }
 
+void shrinkWindowForAngles(Window &oOutExtent, int nX, int nY,
+                           double startAngle, double endAngle)
+{
+    if (startAngle == endAngle)
+        return;
+
+    Window win = oOutExtent;
+
+    // Set the X boundaries for the angles
+    //ABELL - Verify for out-of-raster.
+    double xStart = horizontalIntersect(startAngle, nX, nY, win.yStart);
+    if (isnan(xStart))
+        xStart = horizontalIntersect(startAngle, nX, nY, win.yStop);
+
+    double xStop = horizontalIntersect(endAngle, nX, nY, win.yStart);
+    if (isnan(xStop))
+        xStop = horizontalIntersect(endAngle, nX, nY, win.yStop);
+
+    double xmax = nX;
+    if (!rayBetween(startAngle, endAngle, 0))
+    {
+        if (!isnan(xStart))
+            xmax = std::max(xmax, static_cast<double>(xStart));
+        if (!isnan(xStop))
+            xmax = std::max(xmax, static_cast<double>(xStop));
+        oOutExtent.xStop =
+            std::min(oOutExtent.xStop, static_cast<int>(std::round(xmax)));
+    }
+    double xmin = nX;
+    if (!rayBetween(startAngle, endAngle, M_PI))
+    {
+        if (!isnan(xStart))
+            xmin = std::min(xmin, xStart);
+        if (!isnan(xStop))
+            xmin = std::min(xmin, xStop);
+        oOutExtent.xStart =
+            std::max(oOutExtent.xStart, static_cast<int>(std::round(xmin)));
+    }
+
+    // Set the Y boundaries for the angles
+    //ABELL - Verify for out-of-raster.
+    double yStart = verticalIntersect(startAngle, nX, nY, win.xStart);
+    if (isnan(yStart))
+        yStart = verticalIntersect(startAngle, nX, nY, win.xStop);
+
+    double yStop = verticalIntersect(endAngle, nX, nY, win.xStart);
+    if (isnan(yStop))
+        yStop = verticalIntersect(endAngle, nX, nY, win.xStop);
+
+    double ymin = nY;
+    if (!rayBetween(startAngle, endAngle, M_PI / 2))
+    {
+        if (!isnan(yStart))
+            ymin = std::min(ymin, static_cast<double>(yStart));
+        if (!isnan(yStop))
+            ymin = std::min(ymin, static_cast<double>(yStop));
+        oOutExtent.yStart =
+            std::max(oOutExtent.yStart, static_cast<int>(std::round(ymin)));
+    }
+    double ymax = nY;
+    if (!rayBetween(startAngle, endAngle, 3 * M_PI / 2))
+    {
+        if (!isnan(yStart))
+            ymax = std::max(ymax, static_cast<double>(yStart));
+        if (!isnan(yStop))
+            ymax = std::max(ymax, static_cast<double>(yStop));
+        oOutExtent.yStop =
+            std::min(oOutExtent.yStop, static_cast<int>(std::round(ymax)));
+    }
+}
+
 }  // unnamed namespace
 
 Viewshed::Viewshed(const Options &opts) : oOpts{opts}
@@ -284,6 +358,8 @@ bool Viewshed::calcExtents(int nX, int nY,
         return false;
     }
 
+    shrinkWindowForAngles(oOutExtent, nX, nY, oOpts.startAngle, oOpts.endAngle);
+
     // normalize horizontal index to [ 0, oOutExtent.xSize() )
     oCurExtent = oOutExtent;
     oCurExtent.shiftX(-oOutExtent.xStart);
@@ -325,6 +401,17 @@ bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
     int nX = static_cast<int>(dfX);
     int nY = static_cast<int>(dfY);
 
+    if (oOpts.startAngle < 0 || oOpts.startAngle >= 360)
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "Start angle out of range. Must be [0, 360).");
+    if (oOpts.endAngle < 0 || oOpts.endAngle >= 360)
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "End angle out of range. Must be [0, 360).");
+
+    // Normalize angle to radians and standard math arrangement.
+    oOpts.startAngle = normalizeAngle(oOpts.startAngle);
+    oOpts.endAngle = normalizeAngle(oOpts.endAngle);
+
     // Must calculate extents in order to make the output dataset.
     if (!calcExtents(nX, nY, adfInvTransform))
         return false;
@@ -370,5 +457,11 @@ double adjustCurveCoeff(double curveCoeff, GDALDatasetH hSrcDS)
     return curveCoeff;
 }
 
+void testShrinkWindowForAngles(Window &oOutExtent, int nX, int nY,
+                               double startAngle, double endAngle)
+{
+    shrinkWindowForAngles(oOutExtent, nX, nY, startAngle, endAngle);
+}
+
 }  // namespace viewshed
 }  // namespace gdal