LCOV - code coverage report
Current view: top level - alg/viewshed - viewshed_executor.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 416 472 88.1 %
Date: 2025-10-22 21:46:17 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          94 : bool valid(int i)
      36             : {
      37          94 :     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          94 : bool invalid(int i)
      44             : {
      45          94 :     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       79945 : double CalcHeightLine(int nDistance, double Za)
      53             : {
      54       79945 :     nDistance = std::abs(nDistance);
      55       79945 :     assert(nDistance != 1);
      56       79945 :     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     2834170 : double CalcHeightEdge(int i, int j, double Za, double Zb)
      73             : {
      74     2834170 :     assert(i != j);
      75     2834170 :     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     2833230 : double doEdge(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
      86             :               double dfLastPrev)
      87             : {
      88     2833230 :     if (nXOffset >= nYOffset)
      89     1172420 :         return CalcHeightEdge(nYOffset, nXOffset, dfLastPrev, dfThisPrev);
      90             :     else
      91     1660810 :         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         234 : ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
     126             :                                    GDALRasterBand &dstBand, int nX, int nY,
     127             :                                    const Window &outExtent,
     128             :                                    const Window &curExtent, const Options &opts,
     129         234 :                                    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         234 :       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     2909840 : void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
     180             :                                  double dfZ)
     181             : {
     182     2909840 :     if (oOpts.outputMode != OutputMode::Normal)
     183             :     {
     184       29041 :         double adjustment = dfZ - dfCellVal;
     185       29041 :         if (adjustment > 0)
     186       15760 :             dfResult += adjustment;
     187             :     }
     188             :     else
     189     2880790 :         dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
     190             :                                                           : oOpts.visibleVal;
     191     2909840 :     dfCellVal = std::max(dfCellVal, dfZ);
     192     2913950 : }
     193             : 
     194             : /// Read a line of raster data.
     195             : ///
     196             : /// @param  nLine  Line number to read.
     197             : /// @param  line   Raster line to fill.
     198             : /// @return  Success or failure.
     199       28942 : bool ViewshedExecutor::readLine(int nLine, std::vector<double> &line)
     200             : {
     201       57867 :     std::lock_guard g(iMutex);
     202             : 
     203       57889 :     if (GDALRasterIO(&m_srcBand, GF_Read, oOutExtent.xStart, nLine,
     204       28938 :                      oOutExtent.xSize(), 1, line.data(), oOutExtent.xSize(), 1,
     205       28925 :                      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       28925 :     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       28905 : bool ViewshedExecutor::writeLine(int nLine, std::vector<double> &vResult)
     222             : {
     223             :     // GDALRasterIO isn't thread-safe.
     224       57784 :     std::lock_guard g(oMutex);
     225             : 
     226       57744 :     if (GDALRasterIO(&m_dstBand, GF_Write, 0, nLine - oOutExtent.yStart,
     227       28871 :                      oOutExtent.xSize(), 1, vResult.data(), oOutExtent.xSize(),
     228       28879 :                      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       28879 :     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  lines  Raster lines to adjust.
     244             : /// @return  Processing limits of the line based on min/max distance.
     245       28925 : LineLimits ViewshedExecutor::adjustHeight(int nYOffset, Lines &lines)
     246             : {
     247       28925 :     LineLimits ll(0, m_nX + 1, m_nX + 1, oCurExtent.xSize());
     248             : 
     249             :     // Find the starting point in the raster (m_nX may be outside)
     250       28859 :     int nXStart = oCurExtent.clampX(m_nX);
     251             : 
     252     5844540 :     const auto CheckNoData = [this](double val)
     253             :     {
     254     5843130 :         if (!m_hasFoundNoData &&
     255     2922550 :             ((m_hasNoData && val == m_noDataValue) || std::isnan(val)))
     256             :         {
     257           0 :             m_hasFoundNoData = true;
     258           0 :             if (m_emitWarningIfNoData)
     259             :             {
     260           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     261             :                          "Nodata value found in input DEM. Output will be "
     262             :                          "likely incorrect");
     263             :             }
     264             :         }
     265     2949470 :     };
     266             : 
     267             :     // If there is a height adjustment factor other than zero or a max distance,
     268             :     // calculate the adjusted height of the cell, stopping if we've exceeded the max
     269             :     // distance.
     270       27459 :     if (static_cast<bool>(m_dfHeightAdjFactor) || oOpts.pitchMasking() ||
     271       56357 :         m_dfMaxDistance2 > 0 || m_dfMinDistance2 > 0)
     272             :     {
     273             :         // Hoist invariants from the loops.
     274       28901 :         const double dfLineX = m_gt[2] * nYOffset;
     275       28870 :         const double dfLineY = m_gt[5] * nYOffset;
     276             : 
     277             :         // Go left
     278       28855 :         double *pdfHeight = lines.cur.data() + nXStart;
     279     1479220 :         for (int nXOffset = nXStart - m_nX; nXOffset >= -m_nX;
     280             :              nXOffset--, pdfHeight--)
     281             :         {
     282     1450070 :             double dfX = m_gt[1] * nXOffset + dfLineX;
     283     1450820 :             double dfY = m_gt[4] * nXOffset + dfLineY;
     284     1455480 :             double dfR2 = dfX * dfX + dfY * dfY;
     285             : 
     286     1455480 :             if (dfR2 < m_dfMinDistance2)
     287           6 :                 ll.leftMin--;
     288     1455470 :             else if (dfR2 > m_dfMaxDistance2)
     289             :             {
     290          26 :                 ll.left = nXOffset + m_nX + 1;
     291          26 :                 break;
     292             :             }
     293             : 
     294     1455450 :             CheckNoData(*pdfHeight);
     295     1449600 :             *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
     296     1449600 :             if (oOpts.pitchMasking())
     297          75 :                 calcPitchMask(*pdfHeight, std::sqrt(dfR2),
     298          75 :                               lines.result[m_nX + nXOffset],
     299          75 :                               lines.pitchMask[m_nX + nXOffset]);
     300             :         }
     301             : 
     302             :         // Go right.
     303       29168 :         pdfHeight = lines.cur.data() + nXStart + 1;
     304     1512460 :         for (int nXOffset = nXStart - m_nX + 1;
     305     1512460 :              nXOffset < oCurExtent.xSize() - m_nX; nXOffset++, pdfHeight++)
     306             :         {
     307     1487180 :             double dfX = m_gt[1] * nXOffset + dfLineX;
     308     1486890 :             double dfY = m_gt[4] * nXOffset + dfLineY;
     309     1489470 :             double dfR2 = dfX * dfX + dfY * dfY;
     310             : 
     311     1489470 :             if (dfR2 < m_dfMinDistance2)
     312           3 :                 ll.rightMin++;
     313     1489470 :             else if (dfR2 > m_dfMaxDistance2)
     314             :             {
     315          26 :                 ll.right = nXOffset + m_nX;
     316          26 :                 break;
     317             :             }
     318             : 
     319     1489450 :             CheckNoData(*pdfHeight);
     320     1488470 :             *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
     321     1488470 :             if (oOpts.pitchMasking())
     322         175 :                 calcPitchMask(*pdfHeight, std::sqrt(dfR2),
     323         175 :                               lines.result[m_nX + nXOffset],
     324         175 :                               lines.pitchMask[m_nX + nXOffset]);
     325             :         }
     326             :     }
     327             :     else
     328             :     {
     329             :         // No curvature adjustment. Just normalize for the observer height.
     330           0 :         double *pdfHeight = lines.cur.data();
     331           0 :         for (int i = 0; i < oCurExtent.xSize(); ++i)
     332             :         {
     333           0 :             CheckNoData(*pdfHeight);
     334           0 :             *pdfHeight -= m_dfZObserver;
     335           0 :             pdfHeight++;
     336             :         }
     337             :     }
     338       28951 :     return ll;
     339             : }
     340             : 
     341             : /// Calculate the pitch masking value to apply after running the viewshed algorithm.
     342             : ///
     343             : /// @param  dfZ  Adjusted input height.
     344             : /// @param  dfDist  Distance from observer to cell.
     345             : /// @param  dfResult  Result value to which adjustment may be added.
     346             : /// @param  maskVal  Output mask value.
     347         250 : void ViewshedExecutor::calcPitchMask(double dfZ, double dfDist, double dfResult,
     348             :                                      double &maskVal)
     349             : {
     350         250 :     if (oOpts.lowPitchMasking())
     351             :     {
     352          25 :         double dfZMask = dfDist * m_lowTanPitch;
     353          25 :         double adjustment = dfZMask - dfZ;
     354          25 :         if (adjustment > 0)
     355             :         {
     356          48 :             maskVal = (oOpts.outputMode == OutputMode::Normal
     357          24 :                            ? std::numeric_limits<double>::infinity()
     358             :                            : adjustment + dfResult);
     359          24 :             return;
     360             :         }
     361             :     }
     362         226 :     if (oOpts.highPitchMasking())
     363             :     {
     364         225 :         double dfZMask = dfDist * m_highTanPitch;
     365         225 :         if (dfZ > dfZMask)
     366           7 :             maskVal = std::numeric_limits<double>::infinity();
     367             :     }
     368             : }
     369             : 
     370             : /// Process the first line (the one with the Y coordinate the same as the observer).
     371             : ///
     372             : /// @param lines  Raster lines to process.
     373             : /// @return True on success, false otherwise.
     374         234 : bool ViewshedExecutor::processFirstLine(Lines &lines)
     375             : {
     376         234 :     int nLine = oOutExtent.clampY(m_nY);
     377         234 :     int nYOffset = nLine - m_nY;
     378             : 
     379         234 :     if (!readLine(nLine, lines.cur))
     380           0 :         return false;
     381             : 
     382             :     // If the observer is outside of the raster, take the specified value as the Z height,
     383             :     // otherwise, take it as an offset from the raster height at that location.
     384         234 :     m_dfZObserver = oOpts.observer.z;
     385         234 :     if (oCurExtent.containsX(m_nX))
     386             :     {
     387         229 :         m_dfZObserver += lines.cur[m_nX];
     388         229 :         if (oOpts.outputMode == OutputMode::Normal)
     389         212 :             lines.result[m_nX] = oOpts.visibleVal;
     390             :     }
     391         234 :     m_dfHeightAdjFactor = calcHeightAdjFactor();
     392             : 
     393             :     // In DEM mode the base is the pre-adjustment value.  In ground mode the base is zero.
     394         234 :     if (oOpts.outputMode == OutputMode::DEM)
     395          16 :         lines.result = lines.cur;
     396         218 :     else if (oOpts.outputMode == OutputMode::Ground)
     397           7 :         std::fill(lines.result.begin(), lines.result.end(), 0);
     398             : 
     399         234 :     LineLimits ll = adjustHeight(nYOffset, lines);
     400         235 :     if (oCurExtent.containsX(m_nX) && ll.leftMin != ll.rightMin)
     401           1 :         lines.result[m_nX] = oOpts.outOfRangeVal;
     402             : 
     403         235 :     if (!oCurExtent.containsY(m_nY))
     404           4 :         processFirstLineTopOrBottom(ll, lines);
     405             :     else
     406             :     {
     407         461 :         CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
     408         461 :         pQueue->SubmitJob([&]() { processFirstLineLeft(ll, lines); });
     409         462 :         pQueue->SubmitJob([&]() { processFirstLineRight(ll, lines); });
     410         231 :         pQueue->WaitCompletion();
     411             :     }
     412             : 
     413         235 :     if (oOpts.pitchMasking())
     414           2 :         applyPitchMask(lines.result, lines.pitchMask);
     415         235 :     if (!writeLine(nLine, lines.result))
     416           0 :         return false;
     417             : 
     418         235 :     return oProgress.lineComplete();
     419             : }
     420             : 
     421             : /// Set the pitch masked value into the result vector when applicable.
     422             : ///
     423             : /// @param  vResult  Result vector.
     424             : /// @param  vPitchMaskVal  Pitch mask values (nan is no masking, inf is out-of-range, else
     425             : ///                        actual value).
     426          20 : void ViewshedExecutor::applyPitchMask(std::vector<double> &vResult,
     427             :                                       const std::vector<double> &vPitchMaskVal)
     428             : {
     429         270 :     for (size_t i = 0; i < vResult.size(); ++i)
     430             :     {
     431         250 :         if (std::isnan(vPitchMaskVal[i]))
     432         219 :             continue;
     433          31 :         if (std::isinf(vPitchMaskVal[i]))
     434           7 :             vResult[i] = oOpts.outOfRangeVal;
     435             :         else
     436          24 :             vResult[i] = vPitchMaskVal[i];
     437             :     }
     438          20 : }
     439             : 
     440             : /// If the observer is above or below the raster, set all cells in the first line near the
     441             : /// observer as observable provided they're in range. Mark cells out of range as such.
     442             : /// @param  ll  Line limits for processing.
     443             : /// @param  lines  Raster lines to process.
     444           4 : void ViewshedExecutor::processFirstLineTopOrBottom(const LineLimits &ll,
     445             :                                                    Lines &lines)
     446             : {
     447           4 :     double *pResult = lines.result.data() + ll.left;
     448           4 :     double *pThis = lines.cur.data() + ll.left;
     449          24 :     for (int iPixel = ll.left; iPixel < ll.right; ++iPixel, ++pResult, pThis++)
     450             :     {
     451          20 :         if (oOpts.outputMode == OutputMode::Normal)
     452           0 :             *pResult = oOpts.visibleVal;
     453             :         else
     454          20 :             setOutput(*pResult, *pThis, *pThis);
     455             :     }
     456             : 
     457           4 :     std::fill(lines.result.begin(), lines.result.begin() + ll.left,
     458           4 :               oOpts.outOfRangeVal);
     459           4 :     std::fill(lines.result.begin() + ll.right,
     460           4 :               lines.result.begin() + oCurExtent.xStop, oOpts.outOfRangeVal);
     461           4 : }
     462             : 
     463             : /// Process the part of the first line to the left of the observer.
     464             : ///
     465             : /// @param ll  Line limits for masking.
     466             : /// @param lines  Raster lines to process.
     467         231 : void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll, Lines &lines)
     468             : {
     469         231 :     int iEnd = ll.left - 1;
     470         231 :     int iStart = m_nX - 1;  // One left of the observer.
     471             : 
     472             :     // If end is to the right of start, everything is taken care of by right processing.
     473         231 :     if (iEnd >= iStart)
     474             :     {
     475          30 :         maskLineLeft(lines.result, ll, m_nY);
     476          30 :         return;
     477             :     }
     478             : 
     479         201 :     iStart = oCurExtent.clampX(iStart);
     480             : 
     481         201 :     double *pThis = lines.cur.data() + iStart;
     482             : 
     483             :     // If the start cell is next to the observer, just mark it visible.
     484         201 :     if (iStart + 1 == m_nX || iStart + 1 == oCurExtent.xStop)
     485             :     {
     486         201 :         double dfZ = *pThis;
     487         201 :         if (oOpts.outputMode == OutputMode::Normal)
     488         190 :             lines.result[iStart] = oOpts.visibleVal;
     489             :         else
     490          11 :             setOutput(lines.result[iStart], *pThis, dfZ);
     491         201 :         iStart--;
     492         201 :         pThis--;
     493             :     }
     494             : 
     495             :     // Go from the observer to the left, calculating Z as we go.
     496       10357 :     for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--)
     497             :     {
     498       10156 :         int nXOffset = std::abs(iPixel - m_nX);
     499       10156 :         double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
     500       10156 :         setOutput(lines.result[iPixel], *pThis, dfZ);
     501             :     }
     502             : 
     503         201 :     maskLineLeft(lines.result, ll, m_nY);
     504             : }
     505             : 
     506             : /// Mask cells based on angle intersection to the left of the observer.
     507             : ///
     508             : /// @param vResult  Result raaster line.
     509             : /// @param nLine  Line number.
     510             : /// @return  True when all cells have been masked.
     511       28916 : bool ViewshedExecutor::maskAngleLeft(std::vector<double> &vResult, int nLine)
     512             : {
     513          70 :     auto clamp = [this](int x)
     514          36 :     { return (x < 0 || x >= m_nX) ? INVALID_ISECT : x; };
     515             : 
     516       28916 :     if (!oOpts.angleMasking())
     517       28845 :         return false;
     518             : 
     519          31 :     if (nLine != m_nY)
     520             :     {
     521             :         int startAngleX =
     522          18 :             clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
     523          18 :         int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
     524             :         // If neither X intersect is in the quadrant and a ray in the quadrant isn't
     525             :         // between start and stop, fill it all and return true.  If it is in between
     526             :         // start and stop, we're done.
     527          18 :         if (invalid(startAngleX) && invalid(endAngleX))
     528             :         {
     529             :             // Choose a test angle in quadrant II or III depending on the line.
     530          15 :             double testAngle = nLine < m_nY ? m_testAngle[2] : m_testAngle[3];
     531          15 :             if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
     532             :             {
     533          10 :                 std::fill(vResult.begin(), vResult.begin() + m_nX,
     534          10 :                           oOpts.outOfRangeVal);
     535          15 :                 return true;
     536             :             }
     537           5 :             return false;
     538             :         }
     539           3 :         if (nLine > m_nY)
     540           0 :             std::swap(startAngleX, endAngleX);
     541           3 :         if (invalid(startAngleX))
     542           3 :             startAngleX = 0;
     543           3 :         if (invalid(endAngleX))
     544           0 :             endAngleX = m_nX - 1;
     545           3 :         if (startAngleX <= endAngleX)
     546             :         {
     547           3 :             std::fill(vResult.begin(), vResult.begin() + startAngleX,
     548           3 :                       oOpts.outOfRangeVal);
     549           3 :             std::fill(vResult.begin() + endAngleX + 1, vResult.begin() + m_nX,
     550           3 :                       oOpts.outOfRangeVal);
     551             :         }
     552             :         else
     553             :         {
     554           0 :             std::fill(vResult.begin() + endAngleX + 1,
     555           0 :                       vResult.begin() + startAngleX, oOpts.outOfRangeVal);
     556             :         }
     557             :     }
     558             :     // nLine == m_nY
     559          13 :     else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI))
     560             :     {
     561           1 :         std::fill(vResult.begin(), vResult.begin() + m_nX, oOpts.outOfRangeVal);
     562           1 :         return true;
     563             :     }
     564           4 :     return false;
     565             : }
     566             : 
     567             : /// Mask cells based on angle intersection to the right of the observer.
     568             : ///
     569             : /// @param vResult  Result raaster line.
     570             : /// @param nLine  Line number.
     571             : /// @return  True when all cells have been masked.
     572       28928 : bool ViewshedExecutor::maskAngleRight(std::vector<double> &vResult, int nLine)
     573             : {
     574       28928 :     int lineLength = static_cast<int>(vResult.size());
     575             : 
     576          54 :     auto clamp = [this, lineLength](int x)
     577          36 :     { return (x <= m_nX || x >= lineLength) ? INVALID_ISECT : x; };
     578             : 
     579       28909 :     if (oOpts.startAngle == oOpts.endAngle)
     580       28884 :         return false;
     581             : 
     582          25 :     if (nLine != m_nY)
     583             :     {
     584             :         int startAngleX =
     585          18 :             clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
     586          18 :         int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
     587             : 
     588             :         // If neither X intersect is in the quadrant and a ray in the quadrant isn't
     589             :         // between start and stop, fill it all and return true.  If it is in between
     590             :         // start and stop, we're done.
     591          18 :         if (invalid(startAngleX) && invalid(endAngleX))
     592             :         {
     593             :             // Choose a test angle in quadrant I or IV depending on the line.
     594          10 :             double testAngle = nLine < m_nY ? m_testAngle[1] : m_testAngle[4];
     595          10 :             if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
     596             :             {
     597           0 :                 std::fill(vResult.begin() + m_nX + 1, vResult.end(),
     598           0 :                           oOpts.outOfRangeVal);
     599          10 :                 return true;
     600             :             }
     601          10 :             return false;
     602             :         }
     603             : 
     604           8 :         if (nLine > m_nY)
     605           0 :             std::swap(startAngleX, endAngleX);
     606           8 :         if (invalid(endAngleX))
     607           0 :             endAngleX = lineLength - 1;
     608           8 :         if (invalid(startAngleX))
     609           8 :             startAngleX = m_nX + 1;
     610           8 :         if (startAngleX <= endAngleX)
     611             :         {
     612           8 :             std::fill(vResult.begin() + m_nX + 1, vResult.begin() + startAngleX,
     613           8 :                       oOpts.outOfRangeVal);
     614           8 :             std::fill(vResult.begin() + endAngleX + 1, vResult.end(),
     615           8 :                       oOpts.outOfRangeVal);
     616             :         }
     617             :         else
     618             :         {
     619           0 :             std::fill(vResult.begin() + endAngleX + 1,
     620           0 :                       vResult.begin() + startAngleX, oOpts.outOfRangeVal);
     621             :         }
     622             :     }
     623             :     // nLine == m_nY
     624           7 :     else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 0))
     625             :     {
     626           1 :         std::fill(vResult.begin() + m_nX + 1, vResult.end(),
     627           1 :                   oOpts.outOfRangeVal);
     628           1 :         return true;
     629             :     }
     630           9 :     return false;
     631             : }
     632             : 
     633             : /// Perform angle and min/max masking to the left of the observer.
     634             : ///
     635             : /// @param vResult  Raster line to mask.
     636             : /// @param ll  Min/max line limits.
     637             : /// @param nLine  Line number.
     638       28915 : void ViewshedExecutor::maskLineLeft(std::vector<double> &vResult,
     639             :                                     const LineLimits &ll, int nLine)
     640             : {
     641             :     // If we've already masked with angles everything, just return.
     642       28915 :     if (maskAngleLeft(vResult, nLine))
     643          11 :         return;
     644             : 
     645             :     // Mask cells from the left edge to the left limit.
     646       28864 :     std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
     647             :     // Mask cells from the left min to the observer.
     648       28789 :     if (ll.leftMin < m_nX)
     649           3 :         std::fill(vResult.begin() + ll.leftMin, vResult.begin() + m_nX,
     650           3 :                   oOpts.outOfRangeVal);
     651             : }
     652             : 
     653             : /// Perform angle and min/max masking to the right of the observer.
     654             : ///
     655             : /// @param vResult  Raster line to mask.
     656             : /// @param ll  Min/max line limits.
     657             : /// @param nLine  Line number.
     658       28934 : void ViewshedExecutor::maskLineRight(std::vector<double> &vResult,
     659             :                                      const LineLimits &ll, int nLine)
     660             : {
     661             :     // If we've already masked with angles everything, just return.
     662       28934 :     if (maskAngleRight(vResult, nLine))
     663           1 :         return;
     664             : 
     665             :     // Mask cells from the observer to right min.
     666       28887 :     std::fill(vResult.begin() + m_nX + 1, vResult.begin() + ll.rightMin,
     667       28899 :               oOpts.outOfRangeVal);
     668             :     // Mask cells from the right limit to the right edge.
     669             : 
     670             :     //
     671             :     //ABELL - Changed from ll.right + 1
     672             :     //
     673       28859 :     if (ll.right <= static_cast<int>(vResult.size()))
     674       28880 :         std::fill(vResult.begin() + ll.right, vResult.end(),
     675       28860 :                   oOpts.outOfRangeVal);
     676             : }
     677             : 
     678             : /// Process the part of the first line to the right of the observer.
     679             : ///
     680             : /// @param ll  Line limits
     681             : /// @param lines  Raster lines to process.
     682         231 : void ViewshedExecutor::processFirstLineRight(const LineLimits &ll, Lines &lines)
     683             : {
     684         231 :     int iStart = m_nX + 1;
     685         231 :     int iEnd = ll.right;
     686             : 
     687             :     // If start is to the right of end, everything is taken care of by left processing.
     688         231 :     if (iStart >= iEnd)
     689             :     {
     690           2 :         maskLineRight(lines.result, ll, m_nY);
     691           2 :         return;
     692             :     }
     693             : 
     694         229 :     iStart = oCurExtent.clampX(iStart);
     695             : 
     696         229 :     double *pThis = lines.cur.data() + iStart;
     697             : 
     698             :     // If the start cell is next to the observer, just mark it visible.
     699         229 :     if (iStart - 1 == m_nX || iStart == oCurExtent.xStart)
     700             :     {
     701         229 :         double dfZ = *pThis;
     702         229 :         if (oOpts.outputMode == OutputMode::Normal)
     703         212 :             lines.result[iStart] = oOpts.visibleVal;
     704             :         else
     705          17 :             setOutput(lines.result[iStart], *pThis, dfZ);
     706         229 :         iStart++;
     707         229 :         pThis++;
     708             :     }
     709             : 
     710             :     // Go from the observer to the right, calculating Z as we go.
     711       10812 :     for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++)
     712             :     {
     713       10583 :         int nXOffset = std::abs(iPixel - m_nX);
     714       10583 :         double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
     715       10583 :         setOutput(lines.result[iPixel], *pThis, dfZ);
     716             :     }
     717             : 
     718         229 :     maskLineRight(lines.result, ll, m_nY);
     719             : }
     720             : 
     721             : /// Process a line to the left of the observer.
     722             : ///
     723             : /// @param nYOffset  Offset of the line being processed from the observer
     724             : /// @param ll  Line limits
     725             : /// @param lines  Raster lines to process.
     726       28596 : void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
     727             :                                        Lines &lines)
     728             : {
     729       28596 :     int iStart = m_nX - 1;
     730       28596 :     int iEnd = ll.left - 1;
     731       28596 :     int nLine = m_nY + nYOffset;
     732             : 
     733             :     // If start to the left of end, everything is taken care of by processing right.
     734       28596 :     if (iStart <= iEnd)
     735             :     {
     736        2921 :         maskLineLeft(lines.result, ll, nLine);
     737        2919 :         return;
     738             :     }
     739       25675 :     iStart = oCurExtent.clampX(iStart);
     740             : 
     741       25684 :     nYOffset = std::abs(nYOffset);
     742       25684 :     double *pThis = lines.cur.data() + iStart;
     743       25644 :     double *pLast = lines.prev.data() + iStart;
     744             : 
     745             :     // If the observer is to the right of the raster, mark the first cell to the left as
     746             :     // visible. This may mark an out-of-range cell with a value, but this will be fixed
     747             :     // with the out of range assignment at the end.
     748             : 
     749       25668 :     if (iStart == oCurExtent.xStop - 1)
     750             :     {
     751           6 :         if (oOpts.outputMode == OutputMode::Normal)
     752           0 :             lines.result[iStart] = oOpts.visibleVal;
     753             :         else
     754           6 :             setOutput(lines.result[iStart], *pThis, *pThis);
     755           6 :         iStart--;
     756           6 :         pThis--;
     757           6 :         pLast--;
     758             :     }
     759             : 
     760             :     // Go from the observer to the left, calculating Z as we go.
     761     1430290 :     for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--, pLast--)
     762             :     {
     763     1404530 :         int nXOffset = std::abs(iPixel - m_nX);
     764             :         double dfZ;
     765     1404530 :         if (nXOffset == nYOffset)
     766             :         {
     767       15643 :             if (nXOffset == 1)
     768         375 :                 dfZ = *pThis;
     769             :             else
     770       15268 :                 dfZ = CalcHeightLine(nXOffset, *(pLast + 1));
     771             :         }
     772             :         else
     773             :             dfZ =
     774     1388880 :                 oZcalc(nXOffset, nYOffset, *(pThis + 1), *pLast, *(pLast + 1));
     775     1407350 :         setOutput(lines.result[iPixel], *pThis, dfZ);
     776             :     }
     777             : 
     778       25765 :     maskLineLeft(lines.result, ll, nLine);
     779             : }
     780             : 
     781             : /// Process a line to the right of the observer.
     782             : ///
     783             : /// @param nYOffset  Offset of the line being processed from the observer
     784             : /// @param ll  Line limits
     785             : /// @param lines  Raster lines to process.
     786       28633 : void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
     787             :                                         Lines &lines)
     788             : {
     789       28633 :     int iStart = m_nX + 1;
     790       28633 :     int iEnd = ll.right;
     791       28633 :     int nLine = m_nY + nYOffset;
     792             : 
     793             :     // If start is to the right of end, everything is taken care of by processing left.
     794       28633 :     if (iStart >= iEnd)
     795             :     {
     796          12 :         maskLineRight(lines.result, ll, nLine);
     797          12 :         return;
     798             :     }
     799       28621 :     iStart = oCurExtent.clampX(iStart);
     800             : 
     801       28649 :     nYOffset = std::abs(nYOffset);
     802       28649 :     double *pThis = lines.cur.data() + iStart;
     803       28587 :     double *pLast = lines.prev.data() + iStart;
     804             : 
     805             :     // If the observer is to the left of the raster, mark the first cell to the right as
     806             :     // visible. This may mark an out-of-range cell with a value, but this will be fixed
     807             :     // with the out of range assignment at the end.
     808       28626 :     if (iStart == 0)
     809             :     {
     810           6 :         if (oOpts.outputMode == OutputMode::Normal)
     811           0 :             lines.result[iStart] = oOpts.visibleVal;
     812             :         else
     813           6 :             setOutput(lines.result[0], *pThis, *pThis);
     814           6 :         iStart++;
     815           6 :         pThis++;
     816           6 :         pLast++;
     817             :     }
     818             : 
     819             :     // Go from the observer to the right, calculating Z as we go.
     820     1500920 :     for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++, pLast++)
     821             :     {
     822     1472240 :         int nXOffset = std::abs(iPixel - m_nX);
     823             :         double dfZ;
     824     1472240 :         if (nXOffset == nYOffset)
     825             :         {
     826       16131 :             if (nXOffset == 1)
     827         416 :                 dfZ = *pThis;
     828             :             else
     829       15715 :                 dfZ = CalcHeightLine(nXOffset, *(pLast - 1));
     830             :         }
     831             :         else
     832             :             dfZ =
     833     1456110 :                 oZcalc(nXOffset, nYOffset, *(pThis - 1), *pLast, *(pLast - 1));
     834     1473470 :         setOutput(lines.result[iPixel], *pThis, dfZ);
     835             :     }
     836             : 
     837       28674 :     maskLineRight(lines.result, ll, nLine);
     838             : }
     839             : 
     840             : /// Apply angular mask to the initial X position.  Assumes m_nX is in the raster.
     841             : /// @param vResult  Raster line on which to apply mask.
     842             : /// @param nLine  Line number.
     843       28642 : void ViewshedExecutor::maskInitial(std::vector<double> &vResult, int nLine)
     844             : {
     845       28642 :     if (!oOpts.angleMasking())
     846       28600 :         return;
     847             : 
     848          17 :     if (nLine < m_nY)
     849             :     {
     850          13 :         if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI / 2))
     851           0 :             vResult[m_nX] = oOpts.outOfRangeVal;
     852             :     }
     853           4 :     else if (nLine > m_nY)
     854             :     {
     855           5 :         if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 3 * M_PI / 2))
     856           0 :             vResult[m_nX] = oOpts.outOfRangeVal;
     857             :     }
     858             : }
     859             : 
     860             : /// Process a line above or below the observer.
     861             : ///
     862             : /// @param nLine  Line number being processed.
     863             : /// @param lines  Raster lines to process.
     864             : /// @return True on success, false otherwise.
     865       28715 : bool ViewshedExecutor::processLine(int nLine, Lines &lines)
     866             : {
     867       28715 :     int nYOffset = nLine - m_nY;
     868             : 
     869       28715 :     if (!readLine(nLine, lines.cur))
     870           0 :         return false;
     871             : 
     872             :     // In DEM mode the base is the input DEM value.
     873             :     // In ground mode the base is zero.
     874       28664 :     if (oOpts.outputMode == OutputMode::DEM)
     875         163 :         lines.result = lines.cur;
     876       28501 :     else if (oOpts.outputMode == OutputMode::Ground)
     877         143 :         std::fill(lines.result.begin(), lines.result.end(), 0);
     878             : 
     879             :     // Adjust height of the read line.
     880       28664 :     LineLimits ll = adjustHeight(nYOffset, lines);
     881             : 
     882             :     // Handle the initial position on the line.
     883       28704 :     if (oCurExtent.containsX(m_nX))
     884             :     {
     885       28641 :         if (ll.left < ll.right && ll.leftMin == ll.rightMin)
     886             :         {
     887             :             double dfZ;
     888       28648 :             if (std::abs(nYOffset) == 1)
     889         414 :                 dfZ = lines.cur[m_nX];
     890             :             else
     891       28234 :                 dfZ = CalcHeightLine(nYOffset, lines.prev[m_nX]);
     892       28640 :             setOutput(lines.result[m_nX], lines.cur[m_nX], dfZ);
     893             :         }
     894             :         else
     895           0 :             lines.result[m_nX] = oOpts.outOfRangeVal;
     896             : 
     897       28623 :         maskInitial(lines.result, nLine);
     898             :     }
     899             : 
     900             :     // process left half then right half of line
     901       57372 :     CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
     902       57238 :     pQueue->SubmitJob([&]() { processLineLeft(nYOffset, ll, lines); });
     903       57311 :     pQueue->SubmitJob([&]() { processLineRight(nYOffset, ll, lines); });
     904       28658 :     pQueue->WaitCompletion();
     905             : 
     906       28690 :     if (oOpts.pitchMasking())
     907          18 :         applyPitchMask(lines.result, lines.pitchMask);
     908       28655 :     if (!writeLine(nLine, lines.result))
     909           0 :         return false;
     910             : 
     911       28636 :     return oProgress.lineComplete();
     912             : }
     913             : 
     914             : // Calculate the ray angle from the origin to middle of the top or bottom
     915             : // of each quadrant.
     916           2 : void ViewshedExecutor::calcTestAngles()
     917             : {
     918             :     // Quadrant 1.
     919             :     {
     920           2 :         int ysize = m_nY + 1;
     921           2 :         int xsize = oCurExtent.xStop - m_nX;
     922           2 :         m_testAngle[1] = atan2(ysize, xsize / 2.0);
     923             :     }
     924             : 
     925             :     // Quadrant 2.
     926             :     {
     927           2 :         int ysize = m_nY + 1;
     928           2 :         int xsize = m_nX + 1;
     929           2 :         m_testAngle[2] = atan2(ysize, -xsize / 2.0);
     930             :     }
     931             : 
     932             :     // Quadrant 3.
     933             :     {
     934           2 :         int ysize = oCurExtent.yStop - m_nY;
     935           2 :         int xsize = m_nX + 1;
     936           2 :         m_testAngle[3] = atan2(-ysize, -xsize / 2.0);
     937             :     }
     938             : 
     939             :     // Quadrant 4.
     940             :     {
     941           2 :         int ysize = oCurExtent.yStop - m_nY;
     942           2 :         int xsize = oCurExtent.xStop - m_nX;
     943           2 :         m_testAngle[4] = atan2(-ysize, xsize / 2.0);
     944             :     }
     945             : 
     946             :     // Adjust range to [0, 2 * M_PI)
     947          10 :     for (int i = 1; i <= 4; ++i)
     948           8 :         if (m_testAngle[i] < 0)
     949           4 :             m_testAngle[i] += (2 * M_PI);
     950           2 : }
     951             : 
     952             : /// Run the viewshed computation
     953             : /// @return  Success as true or false.
     954         235 : bool ViewshedExecutor::run()
     955             : {
     956             :     // If we're doing angular masking, calculate the test angles used later.
     957         235 :     if (oOpts.angleMasking())
     958           2 :         calcTestAngles();
     959             : 
     960         469 :     Lines firstLine(oCurExtent.xSize());
     961         235 :     if (oOpts.pitchMasking())
     962           2 :         firstLine.pitchMask.resize(oOutExtent.xSize(),
     963           4 :                                    std::numeric_limits<double>::quiet_NaN());
     964             : 
     965         234 :     if (!processFirstLine(firstLine))
     966           0 :         return false;
     967             : 
     968         235 :     if (oOpts.cellMode == CellMode::Edge)
     969         235 :         oZcalc = doEdge;
     970           0 :     else if (oOpts.cellMode == CellMode::Diagonal)
     971           0 :         oZcalc = doDiagonal;
     972           0 :     else if (oOpts.cellMode == CellMode::Min)
     973           0 :         oZcalc = doMin;
     974           0 :     else if (oOpts.cellMode == CellMode::Max)
     975           0 :         oZcalc = doMax;
     976             : 
     977             :     // scan upwards
     978         235 :     int yStart = oCurExtent.clampY(m_nY);
     979         235 :     std::atomic<bool> err(false);
     980         235 :     CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
     981         235 :     pQueue->SubmitJob(
     982         235 :         [&]()
     983             :         {
     984         471 :             Lines lines(oCurExtent.xSize());
     985         235 :             lines.prev = firstLine.cur;
     986         235 :             if (oOpts.pitchMasking())
     987           2 :                 lines.pitchMask.resize(
     988           2 :                     oOutExtent.xSize(),
     989           4 :                     std::numeric_limits<double>::quiet_NaN());
     990             : 
     991       13517 :             for (int nLine = yStart - 1; nLine >= oCurExtent.yStart && !err;
     992             :                  nLine--)
     993             :             {
     994       13277 :                 if (!processLine(nLine, lines))
     995           0 :                     err = true;
     996       13281 :                 lines.prev = lines.cur;
     997       13276 :                 if (oOpts.pitchMasking())
     998           9 :                     std::fill(lines.pitchMask.begin(), lines.pitchMask.end(),
     999          14 :                               std::numeric_limits<double>::quiet_NaN());
    1000             :             }
    1001         235 :         });
    1002             : 
    1003             :     // scan downwards
    1004         235 :     pQueue->SubmitJob(
    1005         235 :         [&]()
    1006             :         {
    1007         469 :             Lines lines(oCurExtent.xSize());
    1008         235 :             lines.prev = firstLine.cur;
    1009         235 :             if (oOpts.pitchMasking())
    1010           2 :                 lines.pitchMask.resize(
    1011           2 :                     oOutExtent.xSize(),
    1012           4 :                     std::numeric_limits<double>::quiet_NaN());
    1013             : 
    1014       15676 :             for (int nLine = yStart + 1; nLine < oCurExtent.yStop && !err;
    1015             :                  nLine++)
    1016             :             {
    1017       15439 :                 if (!processLine(nLine, lines))
    1018           0 :                     err = true;
    1019       15442 :                 lines.prev = lines.cur;
    1020       15440 :                 if (oOpts.pitchMasking())
    1021           9 :                     std::fill(lines.pitchMask.begin(), lines.pitchMask.end(),
    1022          19 :                               std::numeric_limits<double>::quiet_NaN());
    1023             :             }
    1024         235 :         });
    1025         235 :     return true;
    1026             : }
    1027             : 
    1028             : }  // namespace viewshed
    1029             : }  // namespace gdal

Generated by: LCOV version 1.14