LCOV - code coverage report
Current view: top level - alg/viewshed - viewshed_executor.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 407 463 87.9 %
Date: 2025-09-10 17:48:50 Functions: 36 40 90.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Viewshed Generation
       4             :  * Purpose:  Core algorithm implementation for viewshed generation.
       5             :  * Author:   Tamas Szekeres, szekerest@gmail.com
       6             :  *
       7             :  * (c) 2024 info@hobu.co
       8             :  *
       9             :  ******************************************************************************
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include <algorithm>
      15             : #include <atomic>
      16             : #include <cassert>
      17             : #include <cmath>
      18             : #include <limits>
      19             : 
      20             : #include "viewshed_executor.h"
      21             : #include "progress.h"
      22             : #include "util.h"
      23             : 
      24             : namespace gdal
      25             : {
      26             : namespace viewshed
      27             : {
      28             : 
      29             : namespace
      30             : {
      31             : 
      32             : /// Determines whether a value is a valid intersection coordinate.
      33             : /// @param  i  Value to test.
      34             : /// @return  True if the value doesn't represent an invalid intersection.
      35          78 : bool valid(int i)
      36             : {
      37          78 :     return i != INVALID_ISECT;
      38             : }
      39             : 
      40             : /// Determines whether a value is an invalid intersection coordinate.
      41             : /// @param  i  Value to test.
      42             : /// @return  True if the value represents an invalid intersection.
      43          78 : bool invalid(int i)
      44             : {
      45          78 :     return !valid(i);
      46             : }
      47             : 
      48             : /// Calculate the height at nDistance units along a line through the origin given the height
      49             : /// at nDistance - 1 units along the line.
      50             : /// \param nDistance  Distance along the line for the target point.
      51             : /// \param Za  Height at the line one unit previous to the target point.
      52       79856 : double CalcHeightLine(int nDistance, double Za)
      53             : {
      54       79856 :     nDistance = std::abs(nDistance);
      55       79856 :     assert(nDistance != 1);
      56       79856 :     return Za * nDistance / (nDistance - 1);
      57             : }
      58             : 
      59             : // Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
      60             : // and passing through the line connecting (i - 1, j, Za) and (i, j - 1, Zb).
      61             : // In other words, the origin and the two points form a plane and we're calculating Zc
      62             : // of the point (i, j, Zc), also on the plane.
      63           0 : double CalcHeightDiagonal(int i, int j, double Za, double Zb)
      64             : {
      65           0 :     return (Za * i + Zb * j) / (i + j - 1);
      66             : }
      67             : 
      68             : // Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
      69             : // and through the line connecting (i -1, j - 1, Za) and (i - 1, j, Zb). In other words,
      70             : // the origin and the other two points form a plane and we're calculating Zc of the
      71             : // point (i, j, Zc), also on the plane.
      72     2843830 : double CalcHeightEdge(int i, int j, double Za, double Zb)
      73             : {
      74     2843830 :     assert(i != j);
      75     2843830 :     return (Za * i + Zb * (j - i)) / (j - 1);
      76             : }
      77             : 
      78           0 : double doDiagonal(int nXOffset, [[maybe_unused]] int nYOffset,
      79             :                   double dfThisPrev, double dfLast,
      80             :                   [[maybe_unused]] double dfLastPrev)
      81             : {
      82           0 :     return CalcHeightDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast);
      83             : }
      84             : 
      85     2845610 : double doEdge(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
      86             :               double dfLastPrev)
      87             : {
      88     2845610 :     if (nXOffset >= nYOffset)
      89     1173890 :         return CalcHeightEdge(nYOffset, nXOffset, dfLastPrev, dfThisPrev);
      90             :     else
      91     1671710 :         return CalcHeightEdge(nXOffset, nYOffset, dfLastPrev, dfLast);
      92             : }
      93             : 
      94           0 : double doMin(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
      95             :              double dfLastPrev)
      96             : {
      97           0 :     double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
      98             :     double dfDiagonal =
      99           0 :         doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
     100           0 :     return std::min(dfEdge, dfDiagonal);
     101             : }
     102             : 
     103           0 : double doMax(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
     104             :              double dfLastPrev)
     105             : {
     106           0 :     double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
     107             :     double dfDiagonal =
     108           0 :         doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
     109           0 :     return std::max(dfEdge, dfDiagonal);
     110             : }
     111             : 
     112             : }  // unnamed namespace
     113             : 
     114             : /// Constructor - the viewshed algorithm executor
     115             : /// @param srcBand  Source raster band
     116             : /// @param dstBand  Destination raster band
     117             : /// @param nX  X position of observer
     118             : /// @param nY  Y position of observer
     119             : /// @param outExtent  Extent of output raster (relative to input)
     120             : /// @param curExtent  Extent of active raster.
     121             : /// @param opts  Configuration options.
     122             : /// @param progress  Reference to the progress tracker.
     123             : /// @param emitWarningIfNoData  Whether a warning must be emitted if an input
     124             : ///                             pixel is at the nodata value.
     125         235 : ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
     126             :                                    GDALRasterBand &dstBand, int nX, int nY,
     127             :                                    const Window &outExtent,
     128             :                                    const Window &curExtent, const Options &opts,
     129         235 :                                    Progress &progress, bool emitWarningIfNoData)
     130             :     : m_pool(4), m_srcBand(srcBand), m_dstBand(dstBand),
     131             :       m_emitWarningIfNoData(emitWarningIfNoData), oOutExtent(outExtent),
     132         235 :       oCurExtent(curExtent), m_nX(nX - oOutExtent.xStart), m_nY(nY),
     133             :       oOpts(opts), oProgress(progress),
     134         235 :       m_dfMinDistance2(opts.minDistance * opts.minDistance),
     135         235 :       m_dfMaxDistance2(opts.maxDistance * opts.maxDistance)
     136             : {
     137         235 :     if (m_dfMaxDistance2 == 0)
     138         232 :         m_dfMaxDistance2 = std::numeric_limits<double>::max();
     139         235 :     if (opts.lowPitch != -90.0)
     140           1 :         m_lowTanPitch = std::tan(oOpts.lowPitch * (2 * M_PI / 360.0));
     141         235 :     if (opts.highPitch != 90.0)
     142           1 :         m_highTanPitch = std::tan(oOpts.highPitch * (2 * M_PI / 360.0));
     143         235 :     m_srcBand.GetDataset()->GetGeoTransform(m_gt);
     144         235 :     int hasNoData = false;
     145         235 :     m_noDataValue = m_srcBand.GetNoDataValue(&hasNoData);
     146         234 :     m_hasNoData = hasNoData;
     147         234 : }
     148             : 
     149             : // calculate the height adjustment factor.
     150         235 : double ViewshedExecutor::calcHeightAdjFactor()
     151             : {
     152         469 :     std::lock_guard g(oMutex);
     153             : 
     154             :     const OGRSpatialReference *poDstSRS =
     155         235 :         m_dstBand.GetDataset()->GetSpatialRef();
     156             : 
     157         234 :     if (poDstSRS)
     158             :     {
     159             :         OGRErr eSRSerr;
     160             : 
     161             :         // If we can't get a SemiMajor axis from the SRS, it will be SRS_WGS84_SEMIMAJOR
     162          12 :         double dfSemiMajor = poDstSRS->GetSemiMajor(&eSRSerr);
     163             : 
     164             :         /* If we fetched the axis from the SRS, use it */
     165          12 :         if (eSRSerr != OGRERR_FAILURE)
     166          12 :             return oOpts.curveCoeff / (dfSemiMajor * 2.0);
     167             : 
     168           0 :         CPLDebug("GDALViewshedGenerate",
     169             :                  "Unable to fetch SemiMajor axis from spatial reference");
     170             :     }
     171         222 :     return 0;
     172             : }
     173             : 
     174             : /// Set the output Z value depending on the observable height and computation mode.
     175             : ///
     176             : /// dfResult  Reference to the result cell
     177             : /// dfCellVal  Reference to the current cell height. Replace with observable height.
     178             : /// dfZ  Minimum observable height at cell.
     179     2916860 : void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
     180             :                                  double dfZ)
     181             : {
     182     2916860 :     if (oOpts.outputMode != OutputMode::Normal)
     183             :     {
     184       28555 :         double adjustment = dfZ - dfCellVal;
     185       28555 :         if (adjustment > 0)
     186       15673 :             dfResult += adjustment;
     187             :     }
     188             :     else
     189     2888310 :         dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
     190             :                                                           : oOpts.visibleVal;
     191     2916860 :     dfCellVal = std::max(dfCellVal, dfZ);
     192     2915270 : }
     193             : 
     194             : /// Read a line of raster data.
     195             : ///
     196             : /// @param  nLine  Line number to read.
     197             : /// @param  data  Pointer to location in which to store data.
     198             : /// @return  Success or failure.
     199       28933 : bool ViewshedExecutor::readLine(int nLine, double *data)
     200             : {
     201       57824 :     std::lock_guard g(iMutex);
     202             : 
     203       28943 :     if (GDALRasterIO(&m_srcBand, GF_Read, oOutExtent.xStart, nLine,
     204             :                      oOutExtent.xSize(), 1, data, oOutExtent.xSize(), 1,
     205       28891 :                      GDT_Float64, 0, 0))
     206             :     {
     207           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     208             :                  "RasterIO error when reading DEM at position (%d,%d), "
     209             :                  "size (%d,%d)",
     210           0 :                  oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
     211           0 :         return false;
     212             :     }
     213       28891 :     return true;
     214             : }
     215             : 
     216             : /// Write an output line of either visibility or height data.
     217             : ///
     218             : /// @param  nLine  Line number being written.
     219             : /// @param vResult  Result line to write.
     220             : /// @return  True on success, false otherwise.
     221       28835 : bool ViewshedExecutor::writeLine(int nLine, std::vector<double> &vResult)
     222             : {
     223             :     // GDALRasterIO isn't thread-safe.
     224       57690 :     std::lock_guard g(oMutex);
     225             : 
     226       57605 :     if (GDALRasterIO(&m_dstBand, GF_Write, 0, nLine - oOutExtent.yStart,
     227       28820 :                      oOutExtent.xSize(), 1, vResult.data(), oOutExtent.xSize(),
     228       28855 :                      1, GDT_Float64, 0, 0))
     229             :     {
     230           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     231             :                  "RasterIO error when writing target raster at position "
     232             :                  "(%d,%d), size (%d,%d)",
     233           0 :                  0, nLine - oOutExtent.yStart, oOutExtent.xSize(), 1);
     234           0 :         return false;
     235             :     }
     236       28855 :     return true;
     237             : }
     238             : 
     239             : /// Adjust the height of the line of data by the observer height and the curvature of the
     240             : /// earth.
     241             : ///
     242             : /// @param  nYOffset  Y offset of the line being adjusted.
     243             : /// @param  vThisLineVal  Line height data.
     244             : /// @param  vPitchMaskVal  Pitch masking line.
     245             : /// @param  vResult  Initialized Result vector.
     246             : /// @return  Processing limits of the line based on min/max distance.
     247       28911 : LineLimits ViewshedExecutor::adjustHeight(int nYOffset,
     248             :                                           std::vector<double> &vThisLineVal,
     249             :                                           const std::vector<double> &vResult,
     250             :                                           std::vector<double> &vPitchMaskVal)
     251             : {
     252       28911 :     LineLimits ll(0, m_nX + 1, m_nX + 1, oCurExtent.xSize());
     253             : 
     254             :     // Find the starting point in the raster (m_nX may be outside)
     255       28848 :     int nXStart = oCurExtent.clampX(m_nX);
     256             : 
     257     5867940 :     const auto CheckNoData = [this](double val)
     258             :     {
     259     5856580 :         if (!m_hasFoundNoData &&
     260     2936070 :             ((m_hasNoData && val == m_noDataValue) || std::isnan(val)))
     261             :         {
     262           0 :             m_hasFoundNoData = true;
     263           0 :             if (m_emitWarningIfNoData)
     264             :             {
     265           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     266             :                          "Nodata value found in input DEM. Output will be "
     267             :                          "likely incorrect");
     268             :             }
     269             :         }
     270     2949340 :     };
     271             : 
     272             :     // If there is a height adjustment factor other than zero or a max distance,
     273             :     // calculate the adjusted height of the cell, stopping if we've exceeded the max
     274             :     // distance.
     275       27418 :     if (static_cast<bool>(m_dfHeightAdjFactor) || oOpts.pitchMasking() ||
     276       56274 :         m_dfMaxDistance2 > 0 || m_dfMinDistance2 > 0)
     277             :     {
     278             :         // Hoist invariants from the loops.
     279       28842 :         const double dfLineX = m_gt[2] * nYOffset;
     280       28820 :         const double dfLineY = m_gt[5] * nYOffset;
     281             : 
     282             :         // Go left
     283       28809 :         double *pdfHeight = vThisLineVal.data() + nXStart;
     284     1478830 :         for (int nXOffset = nXStart - m_nX; nXOffset >= -m_nX;
     285             :              nXOffset--, pdfHeight--)
     286             :         {
     287     1450350 :             double dfX = m_gt[1] * nXOffset + dfLineX;
     288     1450320 :             double dfY = m_gt[4] * nXOffset + dfLineY;
     289     1454060 :             double dfR2 = dfX * dfX + dfY * dfY;
     290             : 
     291     1454060 :             if (dfR2 < m_dfMinDistance2)
     292           6 :                 ll.leftMin--;
     293     1454050 :             else if (dfR2 > m_dfMaxDistance2)
     294             :             {
     295          26 :                 ll.left = nXOffset + m_nX + 1;
     296          26 :                 break;
     297             :             }
     298             : 
     299     1454030 :             CheckNoData(*pdfHeight);
     300     1444660 :             *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
     301     1444660 :             if (oOpts.pitchMasking())
     302          75 :                 calcPitchMask(*pdfHeight, std::sqrt(dfR2),
     303          75 :                               vResult[m_nX + nXOffset],
     304          75 :                               vPitchMaskVal[m_nX + nXOffset]);
     305             :         }
     306             : 
     307             :         // Go right.
     308       28511 :         pdfHeight = vThisLineVal.data() + nXStart + 1;
     309     1521950 :         for (int nXOffset = nXStart - m_nX + 1;
     310     1521950 :              nXOffset < oCurExtent.xSize() - m_nX; nXOffset++, pdfHeight++)
     311             :         {
     312     1491380 :             double dfX = m_gt[1] * nXOffset + dfLineX;
     313     1490940 :             double dfY = m_gt[4] * nXOffset + dfLineY;
     314     1492830 :             double dfR2 = dfX * dfX + dfY * dfY;
     315             : 
     316     1492830 :             if (dfR2 < m_dfMinDistance2)
     317           3 :                 ll.rightMin++;
     318     1492830 :             else if (dfR2 > m_dfMaxDistance2)
     319             :             {
     320          26 :                 ll.right = nXOffset + m_nX;
     321          26 :                 break;
     322             :             }
     323             : 
     324     1492810 :             CheckNoData(*pdfHeight);
     325     1492280 :             *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
     326     1492280 :             if (oOpts.pitchMasking())
     327         175 :                 calcPitchMask(*pdfHeight, std::sqrt(dfR2),
     328         175 :                               vResult[m_nX + nXOffset],
     329         175 :                               vPitchMaskVal[m_nX + nXOffset]);
     330             :         }
     331             :     }
     332             :     else
     333             :     {
     334             :         // No curvature adjustment. Just normalize for the observer height.
     335          14 :         double *pdfHeight = vThisLineVal.data();
     336           0 :         for (int i = 0; i < oCurExtent.xSize(); ++i)
     337             :         {
     338           0 :             CheckNoData(*pdfHeight);
     339           0 :             *pdfHeight -= m_dfZObserver;
     340           0 :             pdfHeight++;
     341             :         }
     342             :     }
     343       28951 :     return ll;
     344             : }
     345             : 
     346             : /// Calculate the pitch masking value to apply after running the viewshed algorithm.
     347             : ///
     348             : /// @param  dfZ  Adjusted input height.
     349             : /// @param  dfDist  Distance from observer to cell.
     350             : /// @param  dfResult  Result value to which adjustment may be added.
     351             : /// @param  maskVal  Output mask value.
     352         250 : void ViewshedExecutor::calcPitchMask(double dfZ, double dfDist, double dfResult,
     353             :                                      double &maskVal)
     354             : {
     355         250 :     if (oOpts.lowPitchMasking())
     356             :     {
     357          25 :         double dfZMask = dfDist * m_lowTanPitch;
     358          25 :         double adjustment = dfZMask - dfZ;
     359          25 :         if (adjustment > 0)
     360             :         {
     361          48 :             maskVal = (oOpts.outputMode == OutputMode::Normal
     362          24 :                            ? std::numeric_limits<double>::infinity()
     363             :                            : adjustment + dfResult);
     364          24 :             return;
     365             :         }
     366             :     }
     367         226 :     if (oOpts.highPitchMasking())
     368             :     {
     369         225 :         double dfZMask = dfDist * m_highTanPitch;
     370         225 :         if (dfZ > dfZMask)
     371           7 :             maskVal = std::numeric_limits<double>::infinity();
     372             :     }
     373             : }
     374             : 
     375             : /// Process the first line (the one with the Y coordinate the same as the observer).
     376             : ///
     377             : /// @param vLastLineVal  Vector in which to store the read line. Becomes the last line
     378             : ///    in further processing.
     379             : /// @return True on success, false otherwise.
     380         235 : bool ViewshedExecutor::processFirstLine(std::vector<double> &vLastLineVal)
     381             : {
     382         235 :     int nLine = oOutExtent.clampY(m_nY);
     383         234 :     int nYOffset = nLine - m_nY;
     384             : 
     385         469 :     std::vector<double> vResult(oOutExtent.xSize());
     386         470 :     std::vector<double> vThisLineVal(oOutExtent.xSize());
     387         470 :     std::vector<double> vPitchMaskVal;
     388         235 :     if (oOpts.pitchMasking())
     389           2 :         vPitchMaskVal.resize(oOutExtent.xSize(),
     390           4 :                              std::numeric_limits<double>::quiet_NaN());
     391             : 
     392         235 :     if (!readLine(nLine, vThisLineVal.data()))
     393           0 :         return false;
     394             : 
     395             :     // If the observer is outside of the raster, take the specified value as the Z height,
     396             :     // otherwise, take it as an offset from the raster height at that location.
     397         234 :     m_dfZObserver = oOpts.observer.z;
     398         234 :     if (oCurExtent.containsX(m_nX))
     399             :     {
     400         229 :         m_dfZObserver += vThisLineVal[m_nX];
     401         228 :         if (oOpts.outputMode == OutputMode::Normal)
     402         212 :             vResult[m_nX] = oOpts.visibleVal;
     403             :     }
     404         233 :     m_dfHeightAdjFactor = calcHeightAdjFactor();
     405             : 
     406             :     // In DEM mode the base is the pre-adjustment value.  In ground mode the base is zero.
     407         234 :     if (oOpts.outputMode == OutputMode::DEM)
     408          16 :         vResult = vThisLineVal;
     409             : 
     410             :     LineLimits ll =
     411         234 :         adjustHeight(nYOffset, vThisLineVal, vResult, vPitchMaskVal);
     412         235 :     if (oCurExtent.containsX(m_nX) && ll.leftMin != ll.rightMin)
     413           1 :         vResult[m_nX] = oOpts.outOfRangeVal;
     414             : 
     415         235 :     if (!oCurExtent.containsY(m_nY))
     416           4 :         processFirstLineTopOrBottom(ll, vResult, vThisLineVal);
     417             :     else
     418             :     {
     419         462 :         CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
     420         231 :         pQueue->SubmitJob([&]()
     421         231 :                           { processFirstLineLeft(ll, vResult, vThisLineVal); });
     422         231 :         pQueue->SubmitJob(
     423         231 :             [&]() { processFirstLineRight(ll, vResult, vThisLineVal); });
     424         231 :         pQueue->WaitCompletion();
     425             :     }
     426             : 
     427             :     // Make the current line the last line.
     428         235 :     vLastLineVal = std::move(vThisLineVal);
     429             : 
     430         235 :     if (oOpts.pitchMasking())
     431           2 :         applyPitchMask(vResult, vPitchMaskVal);
     432         235 :     if (!writeLine(nLine, vResult))
     433           0 :         return false;
     434             : 
     435         235 :     return oProgress.lineComplete();
     436             : }
     437             : 
     438             : /// Set the pitch masked value into the result vector when applicable.
     439             : ///
     440             : /// @param  vResult  Result vector.
     441             : /// @param  vPitchMaskVal  Pitch mask values (nan is no masking, inf is out-of-range, else
     442             : ///                        actual value).
     443          20 : void ViewshedExecutor::applyPitchMask(std::vector<double> &vResult,
     444             :                                       const std::vector<double> &vPitchMaskVal)
     445             : {
     446         270 :     for (size_t i = 0; i < vResult.size(); ++i)
     447             :     {
     448         250 :         if (std::isnan(vPitchMaskVal[i]))
     449         219 :             continue;
     450          31 :         if (std::isinf(vPitchMaskVal[i]))
     451           7 :             vResult[i] = oOpts.outOfRangeVal;
     452             :         else
     453          24 :             vResult[i] = vPitchMaskVal[i];
     454             :     }
     455          20 : }
     456             : 
     457             : /// If the observer is above or below the raster, set all cells in the first line near the
     458             : /// observer as observable provided they're in range. Mark cells out of range as such.
     459             : /// @param  ll  Line limits for processing.
     460             : /// @param  vResult  Result line.
     461             : /// @param  vThisLineVal  Heights of the cells in the target line
     462           4 : void ViewshedExecutor::processFirstLineTopOrBottom(
     463             :     const LineLimits &ll, std::vector<double> &vResult,
     464             :     std::vector<double> &vThisLineVal)
     465             : {
     466           4 :     double *pResult = vResult.data() + ll.left;
     467           4 :     double *pThis = vThisLineVal.data() + ll.left;
     468          24 :     for (int iPixel = ll.left; iPixel < ll.right; ++iPixel, ++pResult, pThis++)
     469             :     {
     470          20 :         if (oOpts.outputMode == OutputMode::Normal)
     471           0 :             *pResult = oOpts.visibleVal;
     472             :         else
     473          20 :             setOutput(*pResult, *pThis, *pThis);
     474             :     }
     475             : 
     476           4 :     std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
     477           4 :     std::fill(vResult.begin() + ll.right, vResult.begin() + oCurExtent.xStop,
     478           4 :               oOpts.outOfRangeVal);
     479           4 : }
     480             : 
     481             : /// Process the part of the first line to the left of the observer.
     482             : ///
     483             : /// @param ll  Line limits for masking.
     484             : /// @param vResult  Vector in which to store the visibility/height results.
     485             : /// @param vThisLineVal  Height of each cell in the line being processed.
     486         231 : void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll,
     487             :                                             std::vector<double> &vResult,
     488             :                                             std::vector<double> &vThisLineVal)
     489             : {
     490         231 :     int iEnd = ll.left - 1;
     491         231 :     int iStart = m_nX - 1;  // One left of the observer.
     492             : 
     493             :     // If end is to the right of start, everything is taken care of by right processing.
     494         231 :     if (iEnd >= iStart)
     495          29 :         return;
     496             : 
     497         202 :     iStart = oCurExtent.clampX(iStart);
     498             : 
     499         201 :     double *pThis = vThisLineVal.data() + iStart;
     500             : 
     501             :     // If the start cell is next to the observer, just mark it visible.
     502         201 :     if (iStart + 1 == m_nX || iStart + 1 == oCurExtent.xStop)
     503             :     {
     504         201 :         double dfZ = *pThis;
     505         201 :         if (oOpts.outputMode == OutputMode::Normal)
     506         190 :             vResult[iStart] = oOpts.visibleVal;
     507             :         else
     508          11 :             setOutput(vResult[iStart], *pThis, dfZ);
     509         201 :         iStart--;
     510         201 :         pThis--;
     511             :     }
     512             : 
     513             :     // Go from the observer to the left, calculating Z as we go.
     514       10336 :     for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--)
     515             :     {
     516       10135 :         int nXOffset = std::abs(iPixel - m_nX);
     517       10135 :         double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
     518       10129 :         setOutput(vResult[iPixel], *pThis, dfZ);
     519             :     }
     520             : 
     521         201 :     maskLineLeft(vResult, ll, m_nY);
     522             : }
     523             : 
     524             : /// Mask cells based on angle intersection to the left of the observer.
     525             : ///
     526             : /// @param vResult  Result raaster line.
     527             : /// @param nLine  Line number.
     528             : /// @return  True when all cells have been masked.
     529       25939 : bool ViewshedExecutor::maskAngleLeft(std::vector<double> &vResult, int nLine)
     530             : {
     531          38 :     auto clamp = [this](int x)
     532          20 :     { return (x < 0 || x >= m_nX) ? INVALID_ISECT : x; };
     533             : 
     534       25939 :     if (!oOpts.angleMasking())
     535       25879 :         return false;
     536             : 
     537          24 :     if (nLine != m_nY)
     538             :     {
     539             :         int startAngleX =
     540          10 :             clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
     541          10 :         int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
     542             :         // If neither X intersect is in the quadrant and a ray in the quadrant isn't
     543             :         // between start and stop, fill it all and return true.  If it is in between
     544             :         // start and stop, we're done.
     545          10 :         if (invalid(startAngleX) && invalid(endAngleX))
     546             :         {
     547             :             // Choose a test angle in quadrant II or III depending on the line.
     548           7 :             double testAngle = nLine < m_nY ? m_testAngle[2] : m_testAngle[3];
     549           7 :             if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
     550             :             {
     551           2 :                 std::fill(vResult.begin(), vResult.begin() + m_nX,
     552           2 :                           oOpts.outOfRangeVal);
     553           7 :                 return true;
     554             :             }
     555           5 :             return false;
     556             :         }
     557           3 :         if (nLine > m_nY)
     558           0 :             std::swap(startAngleX, endAngleX);
     559           3 :         if (invalid(startAngleX))
     560           3 :             startAngleX = 0;
     561           3 :         if (invalid(endAngleX))
     562           0 :             endAngleX = m_nX - 1;
     563           3 :         if (startAngleX <= endAngleX)
     564             :         {
     565           3 :             std::fill(vResult.begin(), vResult.begin() + startAngleX,
     566           3 :                       oOpts.outOfRangeVal);
     567           3 :             std::fill(vResult.begin() + endAngleX + 1, vResult.begin() + m_nX,
     568           3 :                       oOpts.outOfRangeVal);
     569             :         }
     570             :         else
     571             :         {
     572           0 :             std::fill(vResult.begin() + endAngleX + 1,
     573           0 :                       vResult.begin() + startAngleX, oOpts.outOfRangeVal);
     574             :         }
     575             :     }
     576             :     // nLine == m_nY
     577          14 :     else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI))
     578             :     {
     579           0 :         std::fill(vResult.begin(), vResult.begin() + m_nX, oOpts.outOfRangeVal);
     580           0 :         return true;
     581             :     }
     582           4 :     return false;
     583             : }
     584             : 
     585             : /// Mask cells based on angle intersection to the right of the observer.
     586             : ///
     587             : /// @param vResult  Result raaster line.
     588             : /// @param nLine  Line number.
     589             : /// @return  True when all cells have been masked.
     590       28899 : bool ViewshedExecutor::maskAngleRight(std::vector<double> &vResult, int nLine)
     591             : {
     592       28899 :     int lineLength = static_cast<int>(vResult.size());
     593             : 
     594          54 :     auto clamp = [this, lineLength](int x)
     595          36 :     { return (x <= m_nX || x >= lineLength) ? INVALID_ISECT : x; };
     596             : 
     597       28892 :     if (oOpts.startAngle == oOpts.endAngle)
     598       28856 :         return false;
     599             : 
     600          36 :     if (nLine != m_nY)
     601             :     {
     602             :         int startAngleX =
     603          18 :             clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
     604          18 :         int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
     605             : 
     606             :         // If neither X intersect is in the quadrant and a ray in the quadrant isn't
     607             :         // between start and stop, fill it all and return true.  If it is in between
     608             :         // start and stop, we're done.
     609          18 :         if (invalid(startAngleX) && invalid(endAngleX))
     610             :         {
     611             :             // Choose a test angle in quadrant I or IV depending on the line.
     612          10 :             double testAngle = nLine < m_nY ? m_testAngle[1] : m_testAngle[4];
     613          10 :             if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
     614             :             {
     615           0 :                 std::fill(vResult.begin() + m_nX + 1, vResult.end(),
     616           0 :                           oOpts.outOfRangeVal);
     617          10 :                 return true;
     618             :             }
     619          10 :             return false;
     620             :         }
     621             : 
     622           8 :         if (nLine > m_nY)
     623           0 :             std::swap(startAngleX, endAngleX);
     624           8 :         if (invalid(endAngleX))
     625           0 :             endAngleX = lineLength - 1;
     626           8 :         if (invalid(startAngleX))
     627           8 :             startAngleX = m_nX + 1;
     628           8 :         if (startAngleX <= endAngleX)
     629             :         {
     630           8 :             std::fill(vResult.begin() + m_nX + 1, vResult.begin() + startAngleX,
     631           8 :                       oOpts.outOfRangeVal);
     632           8 :             std::fill(vResult.begin() + endAngleX + 1, vResult.end(),
     633           8 :                       oOpts.outOfRangeVal);
     634             :         }
     635             :         else
     636             :         {
     637           0 :             std::fill(vResult.begin() + endAngleX + 1,
     638           0 :                       vResult.begin() + startAngleX, oOpts.outOfRangeVal);
     639             :         }
     640             :     }
     641             :     // nLine == m_nY
     642          18 :     else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 0))
     643             :     {
     644           1 :         std::fill(vResult.begin() + m_nX + 1, vResult.end(),
     645           1 :                   oOpts.outOfRangeVal);
     646           1 :         return true;
     647             :     }
     648           9 :     return false;
     649             : }
     650             : 
     651             : /// Perform angle and min/max masking to the left of the observer.
     652             : ///
     653             : /// @param vResult  Raster line to mask.
     654             : /// @param ll  Min/max line limits.
     655             : /// @param nLine  Line number.
     656       25961 : void ViewshedExecutor::maskLineLeft(std::vector<double> &vResult,
     657             :                                     const LineLimits &ll, int nLine)
     658             : {
     659             :     // If we've already masked with angles everything, just return.
     660       25961 :     if (maskAngleLeft(vResult, nLine))
     661           2 :         return;
     662             : 
     663             :     // Mask cells from the left edge to the left limit.
     664       25911 :     std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
     665             :     // Mask cells from the left min to the observer.
     666       25818 :     if (ll.leftMin < m_nX)
     667           3 :         std::fill(vResult.begin() + ll.leftMin, vResult.begin() + m_nX,
     668           3 :                   oOpts.outOfRangeVal);
     669             : }
     670             : 
     671             : /// Perform angle and min/max masking to the right of the observer.
     672             : ///
     673             : /// @param vResult  Raster line to mask.
     674             : /// @param ll  Min/max line limits.
     675             : /// @param nLine  Line number.
     676       28914 : void ViewshedExecutor::maskLineRight(std::vector<double> &vResult,
     677             :                                      const LineLimits &ll, int nLine)
     678             : {
     679             :     // If we've already masked with angles everything, just return.
     680       28914 :     if (maskAngleRight(vResult, nLine))
     681           1 :         return;
     682             : 
     683             :     // Mask cells from the observer to right min.
     684       28857 :     std::fill(vResult.begin() + m_nX + 1, vResult.begin() + ll.rightMin,
     685       28874 :               oOpts.outOfRangeVal);
     686             :     // Mask cells from the right limit to the right edge.
     687       28833 :     if (ll.right + 1 < static_cast<int>(vResult.size()))
     688           8 :         std::fill(vResult.begin() + ll.right + 1, vResult.end(),
     689           8 :                   oOpts.outOfRangeVal);
     690             : }
     691             : 
     692             : /// Process the part of the first line to the right of the observer.
     693             : ///
     694             : /// @param ll  Line limits
     695             : /// @param vResult  Vector in which to store the visibility/height results.
     696             : /// @param vThisLineVal  Height of each cell in the line being processed.
     697         230 : void ViewshedExecutor::processFirstLineRight(const LineLimits &ll,
     698             :                                              std::vector<double> &vResult,
     699             :                                              std::vector<double> &vThisLineVal)
     700             : {
     701         230 :     int iStart = m_nX + 1;
     702         230 :     int iEnd = ll.right;
     703             : 
     704             :     // If start is to the right of end, everything is taken care of by left processing.
     705         230 :     if (iStart >= iEnd)
     706           2 :         return;
     707             : 
     708         228 :     iStart = oCurExtent.clampX(iStart);
     709             : 
     710         229 :     double *pThis = vThisLineVal.data() + iStart;
     711             : 
     712             :     // If the start cell is next to the observer, just mark it visible.
     713         228 :     if (iStart - 1 == m_nX || iStart == oCurExtent.xStart)
     714             :     {
     715         228 :         double dfZ = *pThis;
     716         228 :         if (oOpts.outputMode == OutputMode::Normal)
     717         211 :             vResult[iStart] = oOpts.visibleVal;
     718             :         else
     719          17 :             setOutput(vResult[iStart], *pThis, dfZ);
     720         228 :         iStart++;
     721         228 :         pThis++;
     722             :     }
     723             : 
     724             :     // Go from the observer to the right, calculating Z as we go.
     725       10783 :     for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++)
     726             :     {
     727       10554 :         int nXOffset = std::abs(iPixel - m_nX);
     728       10554 :         double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
     729       10555 :         setOutput(vResult[iPixel], *pThis, dfZ);
     730             :     }
     731             : 
     732         229 :     maskLineRight(vResult, ll, m_nY);
     733             : }
     734             : 
     735             : /// Process a line to the left of the observer.
     736             : ///
     737             : /// @param nYOffset  Offset of the line being processed from the observer
     738             : /// @param ll  Line limits
     739             : /// @param vResult  Vector in which to store the visibility/height results.
     740             : /// @param vThisLineVal  Height of each cell in the line being processed.
     741             : /// @param vLastLineVal  Observable height of each cell in the previous line processed.
     742       28624 : void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
     743             :                                        std::vector<double> &vResult,
     744             :                                        std::vector<double> &vThisLineVal,
     745             :                                        std::vector<double> &vLastLineVal)
     746             : {
     747       28624 :     int iStart = m_nX - 1;
     748       28624 :     int iEnd = ll.left - 1;
     749       28624 :     int nLine = m_nY + nYOffset;
     750             :     // If start to the left of end, everything is taken care of by processing right.
     751       28624 :     if (iStart <= iEnd)
     752        2918 :         return;
     753       25706 :     iStart = oCurExtent.clampX(iStart);
     754             : 
     755       25640 :     nYOffset = std::abs(nYOffset);
     756       25640 :     double *pThis = vThisLineVal.data() + iStart;
     757       25560 :     double *pLast = vLastLineVal.data() + iStart;
     758             : 
     759             :     // If the observer is to the right of the raster, mark the first cell to the left as
     760             :     // visible. This may mark an out-of-range cell with a value, but this will be fixed
     761             :     // with the out of range assignment at the end.
     762             : 
     763       25588 :     if (iStart == oCurExtent.xStop - 1)
     764             :     {
     765           6 :         if (oOpts.outputMode == OutputMode::Normal)
     766           0 :             vResult[iStart] = oOpts.visibleVal;
     767             :         else
     768           6 :             setOutput(vResult[iStart], *pThis, *pThis);
     769           6 :         iStart--;
     770           6 :         pThis--;
     771           6 :         pLast--;
     772             :     }
     773             : 
     774             :     // Go from the observer to the left, calculating Z as we go.
     775     1442780 :     for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--, pLast--)
     776             :     {
     777     1417010 :         int nXOffset = std::abs(iPixel - m_nX);
     778             :         double dfZ;
     779     1417010 :         if (nXOffset == nYOffset)
     780             :         {
     781       15647 :             if (nXOffset == 1)
     782         375 :                 dfZ = *pThis;
     783             :             else
     784       15272 :                 dfZ = CalcHeightLine(nXOffset, *(pLast + 1));
     785             :         }
     786             :         else
     787             :             dfZ =
     788     1401360 :                 oZcalc(nXOffset, nYOffset, *(pThis + 1), *pLast, *(pLast + 1));
     789     1422160 :         setOutput(vResult[iPixel], *pThis, dfZ);
     790             :     }
     791             : 
     792       25774 :     maskLineLeft(vResult, ll, nLine);
     793             : }
     794             : 
     795             : /// Process a line to the right of the observer.
     796             : ///
     797             : /// @param nYOffset  Offset of the line being processed from the observer
     798             : /// @param ll  Line limits
     799             : /// @param vResult  Vector in which to store the visibility/height results.
     800             : /// @param vThisLineVal  Height of each cell in the line being processed.
     801             : /// @param vLastLineVal  Observable height of each cell in the previous line processed.
     802       28640 : void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
     803             :                                         std::vector<double> &vResult,
     804             :                                         std::vector<double> &vThisLineVal,
     805             :                                         std::vector<double> &vLastLineVal)
     806             : {
     807       28640 :     int iStart = m_nX + 1;
     808       28640 :     int iEnd = ll.right;
     809       28640 :     int nLine = m_nY + nYOffset;
     810             : 
     811             :     // If start is to the right of end, everything is taken care of by processing left.
     812       28640 :     if (iStart >= iEnd)
     813          12 :         return;
     814       28628 :     iStart = oCurExtent.clampX(iStart);
     815             : 
     816       28612 :     nYOffset = std::abs(nYOffset);
     817       28612 :     double *pThis = vThisLineVal.data() + iStart;
     818       28587 :     double *pLast = vLastLineVal.data() + iStart;
     819             : 
     820             :     // If the observer is to the left of the raster, mark the first cell to the right as
     821             :     // visible. This may mark an out-of-range cell with a value, but this will be fixed
     822             :     // with the out of range assignment at the end.
     823       28593 :     if (iStart == 0)
     824             :     {
     825           6 :         if (oOpts.outputMode == OutputMode::Normal)
     826           0 :             vResult[iStart] = oOpts.visibleVal;
     827             :         else
     828           6 :             setOutput(vResult[0], *pThis, *pThis);
     829           6 :         iStart++;
     830           6 :         pThis++;
     831           6 :         pLast++;
     832             :     }
     833             : 
     834             :     // Go from the observer to the right, calculating Z as we go.
     835     1493160 :     for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++, pLast++)
     836             :     {
     837     1464490 :         int nXOffset = std::abs(iPixel - m_nX);
     838             :         double dfZ;
     839     1464490 :         if (nXOffset == nYOffset)
     840             :         {
     841       16127 :             if (nXOffset == 1)
     842         415 :                 dfZ = *pThis;
     843             :             else
     844       15712 :                 dfZ = CalcHeightLine(nXOffset, *(pLast - 1));
     845             :         }
     846             :         else
     847             :             dfZ =
     848     1448360 :                 oZcalc(nXOffset, nYOffset, *(pThis - 1), *pLast, *(pLast - 1));
     849     1460310 :         setOutput(vResult[iPixel], *pThis, dfZ);
     850             :     }
     851             : 
     852       28672 :     maskLineRight(vResult, ll, nLine);
     853             : }
     854             : 
     855             : /// Apply angular mask to the initial X position.  Assumes m_nX is in the raster.
     856             : /// @param vResult  Raster line on which to apply mask.
     857             : /// @param nLine  Line number.
     858       28578 : void ViewshedExecutor::maskInitial(std::vector<double> &vResult, int nLine)
     859             : {
     860       28578 :     if (!oOpts.angleMasking())
     861       28543 :         return;
     862             : 
     863          27 :     if (nLine < m_nY)
     864             :     {
     865          13 :         if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI / 2))
     866           0 :             vResult[m_nX] = oOpts.outOfRangeVal;
     867             :     }
     868          14 :     else if (nLine > m_nY)
     869             :     {
     870           5 :         if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 3 * M_PI / 2))
     871           0 :             vResult[m_nX] = oOpts.outOfRangeVal;
     872             :     }
     873             : }
     874             : 
     875             : /// Process a line above or below the observer.
     876             : ///
     877             : /// @param nLine  Line number being processed.
     878             : /// @param vLastLineVal  Vector in which to store the read line. Becomes the last line
     879             : ///    in further processing.
     880             : /// @return True on success, false otherwise.
     881       28722 : bool ViewshedExecutor::processLine(int nLine, std::vector<double> &vLastLineVal)
     882             : {
     883       28722 :     int nYOffset = nLine - m_nY;
     884       57439 :     std::vector<double> vResult(oOutExtent.xSize());
     885       57435 :     std::vector<double> vThisLineVal(oOutExtent.xSize());
     886       57438 :     std::vector<double> vPitchMaskVal;
     887       28710 :     if (oOpts.pitchMasking())
     888          18 :         vPitchMaskVal.resize(oOutExtent.xSize(),
     889          36 :                              std::numeric_limits<double>::quiet_NaN());
     890             : 
     891       28704 :     if (!readLine(nLine, vThisLineVal.data()))
     892           0 :         return false;
     893             : 
     894             :     // In DEM mode the base is the input DEM value.
     895             :     // In ground mode the base is zero.
     896       28670 :     if (oOpts.outputMode == OutputMode::DEM)
     897         163 :         vResult = vThisLineVal;
     898             : 
     899             :     // Adjust height of the read line.
     900             :     LineLimits ll =
     901       28670 :         adjustHeight(nYOffset, vThisLineVal, vResult, vPitchMaskVal);
     902             : 
     903             :     // Handle the initial position on the line.
     904       28714 :     if (oCurExtent.containsX(m_nX))
     905             :     {
     906       28617 :         if (ll.left < ll.right && ll.leftMin == ll.rightMin)
     907             :         {
     908             :             double dfZ;
     909       28603 :             if (std::abs(nYOffset) == 1)
     910         414 :                 dfZ = vThisLineVal[m_nX];
     911             :             else
     912       28189 :                 dfZ = CalcHeightLine(nYOffset, vLastLineVal[m_nX]);
     913       28594 :             setOutput(vResult[m_nX], vThisLineVal[m_nX], dfZ);
     914             :         }
     915             :         else
     916          14 :             vResult[m_nX] = oOpts.outOfRangeVal;
     917             : 
     918       28562 :         maskInitial(vResult, nLine);
     919             :     }
     920             : 
     921             :     // process left half then right half of line
     922       57376 :     CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
     923       28537 :     pQueue->SubmitJob(
     924       28589 :         [&]() {
     925       28589 :             processLineLeft(nYOffset, ll, vResult, vThisLineVal, vLastLineVal);
     926       28555 :         });
     927       28643 :     pQueue->SubmitJob(
     928       28645 :         [&]() {
     929       28645 :             processLineRight(nYOffset, ll, vResult, vThisLineVal, vLastLineVal);
     930       28635 :         });
     931       28649 :     pQueue->WaitCompletion();
     932             :     // Make the current line the last line.
     933       28664 :     vLastLineVal = std::move(vThisLineVal);
     934             : 
     935       28583 :     if (oOpts.pitchMasking())
     936          18 :         applyPitchMask(vResult, vPitchMaskVal);
     937       28556 :     if (!writeLine(nLine, vResult))
     938           0 :         return false;
     939             : 
     940       28565 :     return oProgress.lineComplete();
     941             : }
     942             : 
     943             : // Calculate the ray angle from the origin to middle of the top or bottom
     944             : // of each quadrant.
     945           2 : void ViewshedExecutor::calcTestAngles()
     946             : {
     947             :     // Quadrant 1.
     948             :     {
     949           2 :         int ysize = m_nY + 1;
     950           2 :         int xsize = oCurExtent.xStop - m_nX;
     951           2 :         m_testAngle[1] = atan2(ysize, xsize / 2.0);
     952             :     }
     953             : 
     954             :     // Quadrant 2.
     955             :     {
     956           2 :         int ysize = m_nY + 1;
     957           2 :         int xsize = m_nX + 1;
     958           2 :         m_testAngle[2] = atan2(ysize, -xsize / 2.0);
     959             :     }
     960             : 
     961             :     // Quadrant 3.
     962             :     {
     963           2 :         int ysize = oCurExtent.yStop - m_nY;
     964           2 :         int xsize = m_nX + 1;
     965           2 :         m_testAngle[3] = atan2(-ysize, -xsize / 2.0);
     966             :     }
     967             : 
     968             :     // Quadrant 4.
     969             :     {
     970           2 :         int ysize = oCurExtent.yStop - m_nY;
     971           2 :         int xsize = oCurExtent.xStop - m_nX;
     972           2 :         m_testAngle[4] = atan2(-ysize, xsize / 2.0);
     973             :     }
     974             : 
     975             :     // Adjust range to [0, 2 * M_PI)
     976          10 :     for (int i = 1; i <= 4; ++i)
     977           8 :         if (m_testAngle[i] < 0)
     978           4 :             m_testAngle[i] += (2 * M_PI);
     979           2 : }
     980             : 
     981             : /// Run the viewshed computation
     982             : /// @return  Success as true or false.
     983         235 : bool ViewshedExecutor::run()
     984             : {
     985             :     // If we're doing angular masking, calculate the test angles used later.
     986         235 :     if (oOpts.angleMasking())
     987           2 :         calcTestAngles();
     988             : 
     989         470 :     std::vector<double> vFirstLineVal(oCurExtent.xSize());
     990         234 :     if (!processFirstLine(vFirstLineVal))
     991           0 :         return false;
     992             : 
     993         235 :     if (oOpts.cellMode == CellMode::Edge)
     994         235 :         oZcalc = doEdge;
     995           0 :     else if (oOpts.cellMode == CellMode::Diagonal)
     996           0 :         oZcalc = doDiagonal;
     997           0 :     else if (oOpts.cellMode == CellMode::Min)
     998           0 :         oZcalc = doMin;
     999           0 :     else if (oOpts.cellMode == CellMode::Max)
    1000           0 :         oZcalc = doMax;
    1001             : 
    1002             :     // scan upwards
    1003         235 :     int yStart = oCurExtent.clampY(m_nY);
    1004         235 :     std::atomic<bool> err(false);
    1005         235 :     CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
    1006         235 :     pQueue->SubmitJob(
    1007         235 :         [&]()
    1008             :         {
    1009         469 :             std::vector<double> vLastLineVal = vFirstLineVal;
    1010             : 
    1011       13516 :             for (int nLine = yStart - 1; nLine >= oCurExtent.yStart && !err;
    1012             :                  nLine--)
    1013       13282 :                 if (!processLine(nLine, vLastLineVal))
    1014           0 :                     err = true;
    1015         235 :         });
    1016             : 
    1017             :     // scan downwards
    1018         235 :     pQueue->SubmitJob(
    1019         235 :         [&]()
    1020             :         {
    1021         470 :             std::vector<double> vLastLineVal = vFirstLineVal;
    1022             : 
    1023       15674 :             for (int nLine = yStart + 1; nLine < oCurExtent.yStop && !err;
    1024             :                  nLine++)
    1025       15439 :                 if (!processLine(nLine, vLastLineVal))
    1026           0 :                     err = true;
    1027         235 :         });
    1028         235 :     return true;
    1029             : }
    1030             : 
    1031             : }  // namespace viewshed
    1032             : }  // namespace gdal

Generated by: LCOV version 1.14