Make standard processing work and tests pass.

Author: abellgithub

alg/viewshed/viewshed_executor.cpp

@@ -26,6 +26,11 @@ namespace gdal
 namespace viewshed
 {
 
+CPLErr DummyBand::IReadBlock(int, int, void *)
+{
+    return static_cast<CPLErr>(CPLE_NotSupported);
+}
+
 namespace
 {
 
@@ -164,23 +169,9 @@ ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
                                    const Window &outExtent,
                                    const Window &curExtent, const Options &opts,
                                    Progress &progress, bool emitWarningIfNoData)
-    : m_pool(4), m_srcBand(srcBand), m_sdBand(srcBand), m_dstBand(dstBand),
-      m_emitWarningIfNoData(emitWarningIfNoData), oOutExtent(outExtent),
-      oCurExtent(curExtent), m_nX(nX - oOutExtent.xStart), m_nY(nY),
-      oOpts(opts), oProgress(progress),
-      m_dfMinDistance2(opts.minDistance * opts.minDistance),
-      m_dfMaxDistance2(opts.maxDistance * opts.maxDistance)
+    : ViewshedExecutor(srcBand, m_dummyBand, dstBand, nX, nY, outExtent,
+                       curExtent, opts, progress, emitWarningIfNoData)
 {
-    if (m_dfMaxDistance2 == 0)
-        m_dfMaxDistance2 = std::numeric_limits<double>::max();
-    if (opts.lowPitch != -90.0)
-        m_lowTanPitch = std::tan(oOpts.lowPitch * (2 * M_PI / 360.0));
-    if (opts.highPitch != 90.0)
-        m_highTanPitch = std::tan(oOpts.highPitch * (2 * M_PI / 360.0));
-    m_srcBand.GetDataset()->GetGeoTransform(m_gt);
-    int hasNoData = false;
-    m_noDataValue = m_srcBand.GetNoDataValue(&hasNoData);
-    m_hasNoData = hasNoData;
 }
 
 // calculate the height adjustment factor.
@@ -212,9 +203,41 @@ double ViewshedExecutor::calcHeightAdjFactor()
 ///
 /// dfResult  Reference to the result cell
 /// dfCellVal  Reference to the current cell height. Replace with observable height.
-/// dfZ  Minimum observable height at cell.
+/// dfSdCellVal  Reference to the SD value of the cell.
+/// dfZ    Minimum observable height at cell.
+/// dfSdZ  Minimum possible observable height at cell.
 void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
                                  double &dfSdCellVal, double dfZ, double dfSdZ)
+{
+    if (!sdMode())
+        setOutputNormal(dfResult, dfCellVal, dfZ);
+    else
+        setOutputSd(dfResult, dfCellVal, dfSdCellVal, dfZ, dfSdZ);
+}
+
+// Set the output Z value depending on the observable height and computation mode.
+///
+/// dfResult  Reference to the result cell
+/// dfCellVal  Reference to the current cell height. Replace with observable height.
+/// dfZ  Minimum observable height at cell.
+void ViewshedExecutor::setOutputNormal(double &dfResult, double &dfCellVal,
+                                       double dfZ)
+{
+    if (oOpts.outputMode != OutputMode::Normal)
+    {
+        double adjustment = dfZ - dfCellVal;
+        if (adjustment > 0)
+            dfResult += adjustment;
+    }
+    else
+        dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
+                                                          : oOpts.visibleVal;
+    dfCellVal = std::max(dfCellVal, dfZ);
+}
+
+void ViewshedExecutor::setOutputSd(double &dfResult, double &dfCellVal,
+                                   double &dfSdCellVal, double dfZ,
+                                   double dfSdZ)
 {
     if (oOpts.outputMode != OutputMode::Normal)
     {
@@ -225,38 +248,28 @@ void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
     }
     else
     {
-        /**
-        std::cerr << "cell height/dfz/dfsdz/sd = " <<
-            (dfCellVal + oOpts.targetHeight) << "/" << dfZ << "/" << dfSdZ << "/" <<
-                dfSdCellVal;
-        **/
-        if (dfZ <= dfCellVal + oOpts.targetHeight)
-        {
-            //            std::cerr << " - visible!\n";
-            dfResult = oOpts.visibleVal;
-            if (dfSdCellVal > 1)      // Changing SD to height.
-                dfSdCellVal = dfSdZ;  // See-through case.
-            else
-                dfSdCellVal = dfCellVal;  // Blocking case
-        }
-        else if (dfSdZ <= dfCellVal + oOpts.targetHeight)
+        double cellHeight = dfCellVal + oOpts.targetHeight;
+        if (cellHeight < dfSdZ)
         {
-            //            std::cerr << " - maybe visible!\n";
-            dfResult = oOpts.maybeVisibleVal;
-            if (dfSdCellVal > 1)      // Changing SD to height.
-                dfSdCellVal = dfSdZ;  // See-through case.
-            else
-                dfSdCellVal = dfCellVal;  // Blocking case
-            dfCellVal = dfZ;
+            dfResult = oOpts.invisibleVal;
+            dfSdCellVal =
+                dfSdZ;  // SD height at cell is projected plane height.
         }
         else
         {
-            //            std::cerr << " - not visible!\n";
-            dfResult = oOpts.invisibleVal;
-            dfCellVal = dfZ;
-            dfSdCellVal = dfSdZ;
+            dfResult =
+                (cellHeight < dfZ ? oOpts.maybeVisibleVal : oOpts.visibleVal);
+
+            // Changing SD to height.
+            if (dfSdCellVal > 1)
+                dfSdCellVal =
+                    dfSdZ;  // SD height at cell is projected plane height (maybe vis).
+            else
+                dfSdCellVal =
+                    dfCellVal;  // SD height is cell height (blocking).
         }
     }
+    dfCellVal = std::max(dfCellVal, dfZ);
 }
 
 /// Read a line of raster data.
@@ -287,25 +300,28 @@ bool ViewshedExecutor::readLine(int nLine, std::vector<Cell> &vThisLine)
     for (size_t i = 0; i < line.size(); ++i)
         vThisLine[i].val = line[i];
 
-    int sdStatus = 0;
-    double nodata = std::numeric_limits<double>::max();
+    if (sdMode())
     {
-        std::lock_guard g(iMutex);
-        nodata = m_sdBand.GetNoDataValue();
-        sdStatus = m_sdBand.RasterIO(
-            GF_Read, oOutExtent.xStart, nLine, oOutExtent.xSize(), 1,
-            line.data(), oOutExtent.xSize(), 1, GDT_Float64, 0, 0, nullptr);
-    }
-    if (sdStatus)
-    {
-        CPLError(CE_Failure, CPLE_AppDefined,
-                 "RasterIO error when reading SD band at position (%d,%d), "
-                 "size (%d,%d)",
-                 oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
-        return false;
+        int sdStatus = 0;
+        double nodata = std::numeric_limits<double>::max();
+        {
+            std::lock_guard g(iMutex);
+            nodata = m_sdBand.GetNoDataValue();
+            sdStatus = m_sdBand.RasterIO(
+                GF_Read, oOutExtent.xStart, nLine, oOutExtent.xSize(), 1,
+                line.data(), oOutExtent.xSize(), 1, GDT_Float64, 0, 0, nullptr);
+        }
+        if (sdStatus)
+        {
+            CPLError(CE_Failure, CPLE_AppDefined,
+                     "RasterIO error when reading SD band at position (%d,%d), "
+                     "size (%d,%d)",
+                     oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
+            return false;
+        }
+        for (size_t i = 0; i < line.size(); ++i)
+            vThisLine[i].sd = (line[i] == nodata ? 1000.0 : line[i]);
     }
-    for (size_t i = 0; i < line.size(); ++i)
-        vThisLine[i].sd = (line[i] == nodata ? 1000.0 : line[i]);
 
     return true;
 }
@@ -1162,5 +1178,11 @@ bool ViewshedExecutor::run()
     return true;
 }
 
+bool ViewshedExecutor::sdMode() const
+{
+    // If the SD band isn't a dummy band, we're in SD mode.
+    return dynamic_cast<DummyBand *>(&m_sdBand) == nullptr;
+}
+
 }  // namespace viewshed
 }  // namespace gdal

alg/viewshed/viewshed_executor.h

@@ -27,6 +27,11 @@ namespace gdal
 namespace viewshed
 {
 
+class DummyBand : public GDALRasterBand
+{
+    virtual CPLErr IReadBlock(int, int, void *) override;
+};
+
 class Progress;
 
 /// Executes a viewshed computation on a source band, placing the result
@@ -60,6 +65,7 @@ class ViewshedExecutor
 
   private:
     CPLWorkerThreadPool m_pool;
+    DummyBand m_dummyBand;  // Something to point at when there's no SD band.
     GDALRasterBand &m_srcBand;
     GDALRasterBand &m_sdBand;
     GDALRasterBand &m_dstBand;
@@ -85,6 +91,9 @@ class ViewshedExecutor
     double calcHeightAdjFactor();
     void setOutput(double &dfResult, double &dfCellVal, double &sdCellVal,
                    double dfZ, double dfSdZ);
+    void setOutputSd(double &dfResult, double &dfCellVal, double &sdCellVal,
+                     double dfZ, double dfSdZ);
+    void setOutputNormal(double &dfResult, double &dfCellVal, double dfZ);
     bool readLine(int nLine, std::vector<Cell> &vThisLine);
     bool writeLine(int nLine, std::vector<Cell> &vThisLine);
     bool processLine(int nLine, std::vector<Cell> &vLastLine);
@@ -121,6 +130,7 @@ class ViewshedExecutor
     void applyPitchMask(std::vector<Cell> &vThisLine,
                         const std::vector<double> &vPitchMaskVal);
     void calcTestAngles();
+    bool sdMode() const;
 };
 
 }  // namespace viewshed

autotest/cpp/test_viewshed.cpp

@@ -304,7 +304,8 @@ TEST(Viewshed, high_mask)
     {
         for (int x = 0; x < 15; ++x)
         {
-            EXPECT_EQ(*e, *o);
+            EXPECT_EQ(*e, *o) << "X/Y exp/act = " << x << "/" << y << "/"
+                              << (int)*e << "/" << (int)*o << "!\n";
             e++;
             o++;
         }