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 : // cppcheck-suppress-begin knownConditionTrueFalse
25 : namespace gdal
26 : {
27 : namespace viewshed
28 : {
29 :
30 : //! @cond Doxygen_Suppress
31 0 : CPLErr DummyBand::IReadBlock(int, int, void *)
32 : {
33 0 : return CE_Failure;
34 : }
35 :
36 : //! @endcond
37 :
38 : namespace
39 : {
40 :
41 : /// Determines whether a value is a valid intersection coordinate.
42 : /// @param i Value to test.
43 : /// @return True if the value doesn't represent an invalid intersection.
44 94 : bool valid(int i)
45 : {
46 94 : return i != INVALID_ISECT;
47 : }
48 :
49 : /// Determines whether a value is an invalid intersection coordinate.
50 : /// @param i Value to test.
51 : /// @return True if the value represents an invalid intersection.
52 94 : bool invalid(int i)
53 : {
54 94 : return !valid(i);
55 : }
56 :
57 : /// Calculate the height at nDistance units along a line through the origin given the height
58 : /// at nDistance - 1 units along the line.
59 : /// \param nDistance Distance along the line for the target point.
60 : /// \param Za Height at the line one unit previous to the target point.
61 80988 : double CalcHeightLine(int nDistance, double Za)
62 : {
63 80988 : assert(nDistance > 1);
64 80988 : return Za * nDistance / (nDistance - 1);
65 : }
66 :
67 : /// Calculate the height at nDistance units along a line through the origin given the height
68 : /// at nDistance - 1 units along the line.
69 : /// \param nDistance Distance along the line for the target point.
70 : /// \param Zcur Height at the line at the target point.
71 : /// \param Za Height at the line one unit previous to the target point.
72 61202 : double CalcHeightLine(int nDistance, double Zcur, double Za)
73 : {
74 61202 : nDistance = std::abs(nDistance);
75 61202 : assert(nDistance > 0);
76 61202 : if (nDistance == 1)
77 1225 : return Zcur;
78 59977 : return CalcHeightLine(nDistance, Za);
79 : }
80 :
81 : // Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
82 : // and passing through the line connecting (i - 1, j, Za) and (i, j - 1, Zb).
83 : // In other words, the origin and the two points form a plane and we're calculating Zc
84 : // of the point (i, j, Zc), also on the plane.
85 0 : double CalcHeightDiagonal(int i, int j, double Za, double Zb)
86 : {
87 0 : return (Za * i + Zb * j) / (i + j - 1);
88 : }
89 :
90 : // Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
91 : // and through the line connecting (i -1, j - 1, Za) and (i - 1, j, Zb). In other words,
92 : // the origin and the other two points form a plane and we're calculating Zc of the
93 : // point (i, j, Zc), also on the plane.
94 2917460 : double CalcHeightEdge(int i, int j, double Za, double Zb)
95 : {
96 2917460 : assert(i != j);
97 2917460 : return (Za * i + Zb * (j - i)) / (j - 1);
98 : }
99 :
100 0 : double doDiagonal(int nXOffset, [[maybe_unused]] int nYOffset,
101 : double dfThisPrev, double dfLast,
102 : [[maybe_unused]] double dfLastPrev)
103 : {
104 0 : return CalcHeightDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast);
105 : }
106 :
107 2917460 : double doEdge(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
108 : double dfLastPrev)
109 : {
110 2917460 : if (nXOffset >= nYOffset)
111 1192650 : return CalcHeightEdge(nYOffset, nXOffset, dfLastPrev, dfThisPrev);
112 : else
113 1724820 : return CalcHeightEdge(nXOffset, nYOffset, dfLastPrev, dfLast);
114 : }
115 :
116 0 : double doMin(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
117 : double dfLastPrev)
118 : {
119 0 : double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
120 : double dfDiagonal =
121 0 : doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
122 0 : return std::min(dfEdge, dfDiagonal);
123 : }
124 :
125 0 : double doMax(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
126 : double dfLastPrev)
127 : {
128 0 : double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
129 : double dfDiagonal =
130 0 : doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
131 0 : return std::max(dfEdge, dfDiagonal);
132 : }
133 :
134 : } // unnamed namespace
135 :
136 : /// Constructor - the viewshed algorithm executor
137 : /// @param srcBand Source raster band
138 : /// @param sdBand Standard-deviation raster band
139 : /// @param dstBand Destination raster band
140 : /// @param nX X position of observer
141 : /// @param nY Y position of observer
142 : /// @param outExtent Extent of output raster (relative to input)
143 : /// @param curExtent Extent of active raster.
144 : /// @param opts Configuration options.
145 : /// @param progress Reference to the progress tracker.
146 : /// @param emitWarningIfNoData Whether a warning must be emitted if an input
147 : /// pixel is at the nodata value.
148 244 : ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
149 : GDALRasterBand &sdBand,
150 : GDALRasterBand &dstBand, int nX, int nY,
151 : const Window &outExtent,
152 : const Window &curExtent, const Options &opts,
153 244 : Progress &progress, bool emitWarningIfNoData)
154 : : m_pool(4), m_dummyBand(), m_srcBand(srcBand), m_sdBand(sdBand),
155 : m_dstBand(dstBand),
156 : // If the standard deviation band isn't a dummy band, we're in SD mode.
157 244 : m_hasSdBand(dynamic_cast<DummyBand *>(&m_sdBand) == nullptr),
158 : m_emitWarningIfNoData(emitWarningIfNoData), oOutExtent(outExtent),
159 244 : oCurExtent(curExtent), m_nX(nX - oOutExtent.xStart), m_nY(nY),
160 : oOpts(opts), oProgress(progress),
161 244 : m_dfMinDistance2(opts.minDistance * opts.minDistance),
162 488 : m_dfMaxDistance2(opts.maxDistance * opts.maxDistance)
163 : {
164 244 : if (m_dfMaxDistance2 == 0)
165 241 : m_dfMaxDistance2 = std::numeric_limits<double>::max();
166 244 : if (opts.lowPitch != -90.0)
167 1 : m_lowTanPitch = std::tan(oOpts.lowPitch * (2 * M_PI / 360.0));
168 244 : if (opts.highPitch != 90.0)
169 1 : m_highTanPitch = std::tan(oOpts.highPitch * (2 * M_PI / 360.0));
170 244 : m_srcBand.GetDataset()->GetGeoTransform(m_gt);
171 244 : int hasNoData = false;
172 244 : m_noDataValue = m_srcBand.GetNoDataValue(&hasNoData);
173 244 : m_hasNoData = hasNoData;
174 244 : }
175 :
176 : /// Constructor - the viewshed algorithm executor
177 : /// @param srcBand Source raster band
178 : /// @param dstBand Destination raster band
179 : /// @param nX X position of observer
180 : /// @param nY Y position of observer
181 : /// @param outExtent Extent of output raster (relative to input)
182 : /// @param curExtent Extent of active raster.
183 : /// @param opts Configuration options.
184 : /// @param progress Reference to the progress tracker.
185 : /// @param emitWarningIfNoData Whether a warning must be emitted if an input
186 : /// pixel is at the nodata value.
187 235 : ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
188 : GDALRasterBand &dstBand, int nX, int nY,
189 : const Window &outExtent,
190 : const Window &curExtent, const Options &opts,
191 235 : Progress &progress, bool emitWarningIfNoData)
192 : : ViewshedExecutor(srcBand, m_dummyBand, dstBand, nX, nY, outExtent,
193 235 : curExtent, opts, progress, emitWarningIfNoData)
194 : {
195 235 : }
196 :
197 : // calculate the height adjustment factor.
198 244 : double ViewshedExecutor::calcHeightAdjFactor()
199 : {
200 488 : std::lock_guard g(oMutex);
201 :
202 : const OGRSpatialReference *poDstSRS =
203 244 : m_dstBand.GetDataset()->GetSpatialRef();
204 :
205 244 : if (poDstSRS)
206 : {
207 : OGRErr eSRSerr;
208 :
209 : // If we can't get a SemiMajor axis from the SRS, it will be SRS_WGS84_SEMIMAJOR
210 13 : double dfSemiMajor = poDstSRS->GetSemiMajor(&eSRSerr);
211 :
212 : /* If we fetched the axis from the SRS, use it */
213 13 : if (eSRSerr != OGRERR_FAILURE)
214 13 : return oOpts.curveCoeff / (dfSemiMajor * 2.0);
215 :
216 0 : CPLDebug("GDALViewshedGenerate",
217 : "Unable to fetch SemiMajor axis from spatial reference");
218 : }
219 231 : return 0;
220 : }
221 :
222 : /// Set the output Z value depending on the observable height and computation mode
223 : /// in normal mode.
224 : ///
225 : /// dfResult Reference to the result cell
226 : /// dfCellVal Reference to the current cell height. Replace with observable height.
227 : /// dfZ Minimum observable height at cell.
228 2985250 : void ViewshedExecutor::setOutputNormal(Lines &lines, int pos, double dfZ)
229 : {
230 2985250 : double &cur = lines.cur[pos];
231 2985250 : double &result = lines.result[pos];
232 :
233 2985250 : if (oOpts.outputMode != OutputMode::Normal)
234 : {
235 29100 : double adjustment = dfZ - cur;
236 29100 : if (adjustment > 0)
237 15777 : result += adjustment;
238 : }
239 : else
240 : {
241 2956150 : double cellHeight = cur + oOpts.targetHeight;
242 2956150 : result = (cellHeight < dfZ) ? oOpts.invisibleVal : oOpts.visibleVal;
243 : }
244 2985250 : cur = std::max(cur, dfZ);
245 2985250 : }
246 :
247 : /// Set the output Z value depending on the observable height and computation when
248 : /// making an standard deviation pass.
249 : ///
250 : /// dfResult Reference to the result cell
251 : /// dfCellVal Reference to the current cell height. Replace with observable height.
252 : /// dfZ Minimum observable height at cell.
253 14487 : void ViewshedExecutor::setOutputSd(Lines &lines, int pos, double dfZ)
254 : {
255 14487 : double &cur = lines.cur[pos];
256 14487 : double &result = lines.result[pos];
257 14487 : double &sd = lines.sd[pos];
258 :
259 14487 : assert(oOpts.outputMode == OutputMode::Normal);
260 14487 : if (result == oOpts.invisibleVal)
261 : {
262 7887 : double cellHeight = cur + oOpts.targetHeight;
263 7887 : if (cellHeight > dfZ)
264 2056 : result = oOpts.maybeVisibleVal;
265 : }
266 :
267 14487 : if (sd <= 1)
268 2323 : cur = std::max(dfZ, cur);
269 : else
270 12164 : cur = dfZ;
271 14487 : }
272 :
273 : /// Read a line of raster data.
274 : ///
275 : /// @param nLine Line number to read.
276 : /// @param lines Raster line to fill.
277 : /// @return Success or failure.
278 29125 : bool ViewshedExecutor::readLine(int nLine, Lines &lines)
279 : {
280 58250 : std::lock_guard g(iMutex);
281 :
282 58250 : if (GDALRasterIO(&m_srcBand, GF_Read, oOutExtent.xStart, nLine,
283 29125 : oOutExtent.xSize(), 1, lines.cur.data(),
284 29125 : oOutExtent.xSize(), 1, GDT_Float64, 0, 0))
285 : {
286 0 : CPLError(CE_Failure, CPLE_AppDefined,
287 : "RasterIO error when reading DEM at position (%d,%d), "
288 : "size (%d,%d)",
289 0 : oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
290 0 : return false;
291 : }
292 :
293 29125 : if (sdMode())
294 : {
295 166 : double nodata = m_sdBand.GetNoDataValue();
296 332 : CPLErr sdStatus = m_sdBand.RasterIO(
297 166 : GF_Read, oOutExtent.xStart, nLine, oOutExtent.xSize(), 1,
298 166 : lines.sd.data(), oOutExtent.xSize(), 1, GDT_Float64, 0, 0, nullptr);
299 166 : if (sdStatus != CE_None)
300 : {
301 0 : CPLError(CE_Failure, CPLE_AppDefined,
302 : "RasterIO error when reading standard deviation band at "
303 : "position (%d,%d), "
304 : "size (%d,%d)",
305 0 : oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
306 0 : return false;
307 : }
308 : // Set the standard deviation to 1000 if nodata is found.
309 14682 : for (size_t i = 0; i < lines.sd.size(); ++i)
310 14516 : if (lines.sd[i] == nodata)
311 476 : lines.sd[i] = 1000.0;
312 : }
313 :
314 : // Initialize the result line.
315 : // In DEM mode the base is the pre-adjustment value. In ground mode the base is zero.
316 29125 : if (oOpts.outputMode == OutputMode::DEM)
317 179 : lines.result = lines.cur;
318 28946 : else if (oOpts.outputMode == OutputMode::Ground)
319 150 : std::fill(lines.result.begin(), lines.result.end(), 0);
320 :
321 29125 : return true;
322 : }
323 :
324 : /// Write an output line of either visibility or height data.
325 : ///
326 : /// @param nLine Line number being written.
327 : /// @param vResult Result line to write.
328 : /// @return True on success, false otherwise.
329 29125 : bool ViewshedExecutor::writeLine(int nLine, std::vector<double> &vResult)
330 : {
331 : // GDALRasterIO isn't thread-safe.
332 58250 : std::lock_guard g(oMutex);
333 :
334 58250 : if (GDALRasterIO(&m_dstBand, GF_Write, 0, nLine - oOutExtent.yStart,
335 29125 : oOutExtent.xSize(), 1, vResult.data(), oOutExtent.xSize(),
336 29125 : 1, GDT_Float64, 0, 0))
337 : {
338 0 : CPLError(CE_Failure, CPLE_AppDefined,
339 : "RasterIO error when writing target raster at position "
340 : "(%d,%d), size (%d,%d)",
341 0 : 0, nLine - oOutExtent.yStart, oOutExtent.xSize(), 1);
342 0 : return false;
343 : }
344 29125 : return true;
345 : }
346 :
347 : /// Adjust the height of the line of data by the observer height and the curvature of the
348 : /// earth.
349 : ///
350 : /// @param nYOffset Y offset of the line being adjusted.
351 : /// @param lines Raster lines to adjust.
352 : /// @return Processing limits of the line based on min/max distance.
353 29125 : LineLimits ViewshedExecutor::adjustHeight(int nYOffset, Lines &lines)
354 : {
355 29125 : LineLimits ll(0, m_nX + 1, m_nX + 1, oCurExtent.xSize());
356 :
357 : // Find the starting point in the raster (m_nX may be outside)
358 29125 : int nXStart = oCurExtent.clampX(m_nX);
359 :
360 5971800 : const auto CheckNoData = [this](double val)
361 : {
362 5971800 : if (!m_hasFoundNoData &&
363 2985900 : ((m_hasNoData && val == m_noDataValue) || std::isnan(val)))
364 : {
365 0 : m_hasFoundNoData = true;
366 0 : if (m_emitWarningIfNoData)
367 : {
368 0 : CPLError(CE_Warning, CPLE_AppDefined,
369 : "Nodata value found in input DEM. Output will be "
370 : "likely incorrect");
371 : }
372 : }
373 3015020 : };
374 :
375 : // If there is a height adjustment factor other than zero or a max distance,
376 : // calculate the adjusted height of the cell, stopping if we've exceeded the max
377 : // distance.
378 27571 : if (static_cast<bool>(m_dfHeightAdjFactor) || oOpts.pitchMasking() ||
379 56696 : m_dfMaxDistance2 > 0 || m_dfMinDistance2 > 0)
380 : {
381 : // Hoist invariants from the loops.
382 29125 : const double dfLineX = m_gt[2] * nYOffset;
383 29125 : const double dfLineY = m_gt[5] * nYOffset;
384 :
385 : // Go left
386 29125 : double *pdfHeight = lines.cur.data() + nXStart;
387 1509090 : for (int nXOffset = nXStart - m_nX; nXOffset >= -m_nX;
388 : nXOffset--, pdfHeight--)
389 : {
390 1479990 : double dfX = m_gt[1] * nXOffset + dfLineX;
391 1479990 : double dfY = m_gt[4] * nXOffset + dfLineY;
392 1479990 : double dfR2 = dfX * dfX + dfY * dfY;
393 :
394 1479990 : if (dfR2 < m_dfMinDistance2)
395 6 : ll.leftMin--;
396 1479980 : else if (dfR2 > m_dfMaxDistance2)
397 : {
398 26 : ll.left = nXOffset + m_nX + 1;
399 26 : break;
400 : }
401 :
402 1479960 : CheckNoData(*pdfHeight);
403 1479960 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
404 1479960 : if (oOpts.pitchMasking())
405 75 : calcPitchMask(*pdfHeight, std::sqrt(dfR2),
406 75 : lines.result[m_nX + nXOffset],
407 75 : lines.pitchMask[m_nX + nXOffset]);
408 : }
409 :
410 : // Go right.
411 29125 : pdfHeight = lines.cur.data() + nXStart + 1;
412 1535060 : for (int nXOffset = nXStart - m_nX + 1;
413 1535060 : nXOffset < oCurExtent.xSize() - m_nX; nXOffset++, pdfHeight++)
414 : {
415 1505960 : double dfX = m_gt[1] * nXOffset + dfLineX;
416 1505960 : double dfY = m_gt[4] * nXOffset + dfLineY;
417 1505960 : double dfR2 = dfX * dfX + dfY * dfY;
418 :
419 1505960 : if (dfR2 < m_dfMinDistance2)
420 3 : ll.rightMin++;
421 1505960 : else if (dfR2 > m_dfMaxDistance2)
422 : {
423 26 : ll.right = nXOffset + m_nX;
424 26 : break;
425 : }
426 :
427 1505940 : CheckNoData(*pdfHeight);
428 1505940 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
429 1505940 : if (oOpts.pitchMasking())
430 175 : calcPitchMask(*pdfHeight, std::sqrt(dfR2),
431 175 : lines.result[m_nX + nXOffset],
432 175 : lines.pitchMask[m_nX + nXOffset]);
433 : }
434 : }
435 : else
436 : {
437 : // No curvature adjustment. Just normalize for the observer height.
438 0 : double *pdfHeight = lines.cur.data();
439 0 : for (int i = 0; i < oCurExtent.xSize(); ++i)
440 : {
441 0 : CheckNoData(*pdfHeight);
442 0 : *pdfHeight -= m_dfZObserver;
443 0 : pdfHeight++;
444 : }
445 : }
446 29125 : return ll;
447 : }
448 :
449 : /// Calculate the pitch masking value to apply after running the viewshed algorithm.
450 : ///
451 : /// @param dfZ Adjusted input height.
452 : /// @param dfDist Distance from observer to cell.
453 : /// @param dfResult Result value to which adjustment may be added.
454 : /// @param maskVal Output mask value.
455 250 : void ViewshedExecutor::calcPitchMask(double dfZ, double dfDist, double dfResult,
456 : double &maskVal)
457 : {
458 250 : if (oOpts.lowPitchMasking())
459 : {
460 25 : double dfZMask = dfDist * m_lowTanPitch;
461 25 : double adjustment = dfZMask - dfZ;
462 25 : if (adjustment > 0)
463 : {
464 48 : maskVal = (oOpts.outputMode == OutputMode::Normal
465 24 : ? std::numeric_limits<double>::infinity()
466 : : adjustment + dfResult);
467 24 : return;
468 : }
469 : }
470 226 : if (oOpts.highPitchMasking())
471 : {
472 225 : double dfZMask = dfDist * m_highTanPitch;
473 225 : if (dfZ > dfZMask)
474 7 : maskVal = std::numeric_limits<double>::infinity();
475 : }
476 : }
477 :
478 : /// Process the first line (the one with the Y coordinate the same as the observer).
479 : ///
480 : /// @param lines Raster lines to process.
481 : /// @return True on success, false otherwise.
482 244 : bool ViewshedExecutor::processFirstLine(Lines &lines)
483 : {
484 244 : int nLine = oOutExtent.clampY(m_nY);
485 244 : int nYOffset = nLine - m_nY;
486 :
487 244 : if (!readLine(nLine, lines))
488 0 : return false;
489 :
490 : // If the observer is outside of the raster, take the specified value as the Z height,
491 : // otherwise, take it as an offset from the raster height at that location.
492 244 : m_dfZObserver = oOpts.observer.z;
493 244 : if (oCurExtent.containsX(m_nX))
494 238 : m_dfZObserver += lines.cur[m_nX];
495 :
496 244 : LineLimits ll = adjustHeight(nYOffset, lines);
497 :
498 488 : std::vector<double> savedInput;
499 244 : if (sdMode())
500 9 : savedInput = lines.cur;
501 :
502 244 : if (oCurExtent.containsX(m_nX))
503 : {
504 238 : if (ll.leftMin != ll.rightMin)
505 1 : lines.result[m_nX] = oOpts.outOfRangeVal;
506 237 : else if (oOpts.outputMode == OutputMode::Normal)
507 220 : lines.result[m_nX] = oOpts.visibleVal;
508 : }
509 :
510 1004 : auto process = [this, &ll, &lines](bool sdCalc)
511 : {
512 253 : if (!oCurExtent.containsY(m_nY))
513 4 : processFirstLineTopOrBottom(ll, lines);
514 : else
515 : {
516 498 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
517 249 : pQueue->SubmitJob([&]()
518 249 : { processFirstLineLeft(ll, lines, sdCalc); });
519 249 : pQueue->SubmitJob([&]()
520 249 : { processFirstLineRight(ll, lines, sdCalc); });
521 249 : pQueue->WaitCompletion();
522 : }
523 253 : };
524 :
525 244 : process(false);
526 244 : lines.prev = lines.cur;
527 244 : if (sdMode())
528 : {
529 9 : lines.cur = std::move(savedInput);
530 9 : process(true);
531 9 : lines.prevTmp = lines.cur;
532 : }
533 :
534 244 : if (oOpts.pitchMasking())
535 2 : applyPitchMask(lines.result, lines.pitchMask);
536 244 : if (!writeLine(nLine, lines.result))
537 0 : return false;
538 :
539 244 : return oProgress.lineComplete();
540 : }
541 :
542 : /// Set the pitch masked value into the result vector when applicable.
543 : ///
544 : /// @param vResult Result vector.
545 : /// @param vPitchMaskVal Pitch mask values (nan is no masking, inf is out-of-range, else
546 : /// actual value).
547 20 : void ViewshedExecutor::applyPitchMask(std::vector<double> &vResult,
548 : const std::vector<double> &vPitchMaskVal)
549 : {
550 270 : for (size_t i = 0; i < vResult.size(); ++i)
551 : {
552 250 : if (std::isnan(vPitchMaskVal[i]))
553 219 : continue;
554 31 : if (std::isinf(vPitchMaskVal[i]))
555 7 : vResult[i] = oOpts.outOfRangeVal;
556 : else
557 24 : vResult[i] = vPitchMaskVal[i];
558 : }
559 20 : }
560 :
561 : /// If the observer is above or below the raster, set all cells in the first line near the
562 : /// observer as observable provided they're in range. Mark cells out of range as such.
563 : /// @param ll Line limits for processing.
564 : /// @param lines Raster lines to process.
565 4 : void ViewshedExecutor::processFirstLineTopOrBottom(const LineLimits &ll,
566 : Lines &lines)
567 : {
568 24 : for (int iPixel = ll.left; iPixel < ll.right; ++iPixel)
569 : {
570 20 : if (oOpts.outputMode == OutputMode::Normal)
571 0 : lines.result[iPixel] = oOpts.visibleVal;
572 : else
573 20 : setOutputNormal(lines, iPixel, lines.cur[iPixel]);
574 : }
575 :
576 4 : std::fill(lines.result.begin(), lines.result.begin() + ll.left,
577 4 : oOpts.outOfRangeVal);
578 4 : std::fill(lines.result.begin() + ll.right,
579 4 : lines.result.begin() + oCurExtent.xStop, oOpts.outOfRangeVal);
580 4 : }
581 :
582 : /// Process the part of the first line to the left of the observer.
583 : ///
584 : /// @param ll Line limits for masking.
585 : /// @param sdCalc True when doing standard deviation calculation.
586 : /// @param lines Raster lines to process.
587 249 : void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll, Lines &lines,
588 : bool sdCalc)
589 : {
590 249 : int iEnd = ll.left - 1;
591 249 : int iStart = m_nX - 1; // One left of the observer.
592 :
593 : // If end is to the right of start, everything is taken care of by right processing.
594 249 : if (iEnd >= iStart)
595 : {
596 40 : maskLineLeft(lines.result, ll, m_nY);
597 40 : return;
598 : }
599 :
600 209 : iStart = oCurExtent.clampX(iStart);
601 :
602 : // If the start cell is next to the observer, just mark it visible.
603 209 : if (iStart + 1 == m_nX || iStart + 1 == oCurExtent.xStop)
604 : {
605 209 : double dfZ = lines.cur[iStart];
606 209 : if (oOpts.outputMode == OutputMode::Normal)
607 : {
608 198 : lines.result[iStart] = oOpts.visibleVal;
609 198 : if (sdCalc)
610 4 : if (lines.sd[iStart] > 1)
611 3 : lines.cur[iStart] =
612 3 : m_dfZObserver; // Should this be a minimum value?
613 : }
614 : else
615 11 : setOutputNormal(lines, iStart, dfZ);
616 209 : iStart--;
617 : }
618 :
619 : // Go from the observer to the left, calculating Z as we go.
620 10501 : for (int iPixel = iStart; iPixel > iEnd; iPixel--)
621 : {
622 10292 : int nXOffset = std::abs(iPixel - m_nX);
623 10292 : double dfZ = CalcHeightLine(nXOffset, lines.cur[iPixel + 1]);
624 10292 : if (!sdCalc)
625 10224 : setOutputNormal(lines, iPixel, dfZ);
626 : else
627 68 : setOutputSd(lines, iPixel, dfZ);
628 : }
629 :
630 209 : maskLineLeft(lines.result, ll, m_nY);
631 : }
632 :
633 : /// Mask cells based on angle intersection to the left of the observer.
634 : ///
635 : /// @param vResult Result raaster line.
636 : /// @param nLine Line number.
637 : /// @return True when all cells have been masked.
638 29287 : bool ViewshedExecutor::maskAngleLeft(std::vector<double> &vResult, int nLine)
639 : {
640 70 : auto clamp = [this](int x)
641 36 : { return (x < 0 || x >= m_nX) ? INVALID_ISECT : x; };
642 :
643 29287 : if (!oOpts.angleMasking())
644 29267 : return false;
645 :
646 20 : if (nLine != m_nY)
647 : {
648 : int startAngleX =
649 18 : clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
650 18 : int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
651 : // If neither X intersect is in the quadrant and a ray in the quadrant isn't
652 : // between start and stop, fill it all and return true. If it is in between
653 : // start and stop, we're done.
654 18 : if (invalid(startAngleX) && invalid(endAngleX))
655 : {
656 : // Choose a test angle in quadrant II or III depending on the line.
657 15 : double testAngle = nLine < m_nY ? m_testAngle[2] : m_testAngle[3];
658 15 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
659 : {
660 10 : std::fill(vResult.begin(), vResult.begin() + m_nX,
661 10 : oOpts.outOfRangeVal);
662 15 : return true;
663 : }
664 5 : return false;
665 : }
666 3 : if (nLine > m_nY)
667 0 : std::swap(startAngleX, endAngleX);
668 3 : if (invalid(startAngleX))
669 3 : startAngleX = 0;
670 3 : if (invalid(endAngleX))
671 0 : endAngleX = m_nX - 1;
672 3 : if (startAngleX <= endAngleX)
673 : {
674 3 : std::fill(vResult.begin(), vResult.begin() + startAngleX,
675 3 : oOpts.outOfRangeVal);
676 3 : std::fill(vResult.begin() + endAngleX + 1, vResult.begin() + m_nX,
677 3 : oOpts.outOfRangeVal);
678 : }
679 : else
680 : {
681 0 : std::fill(vResult.begin() + endAngleX + 1,
682 0 : vResult.begin() + startAngleX, oOpts.outOfRangeVal);
683 : }
684 : }
685 : // nLine == m_nY
686 2 : else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI))
687 : {
688 1 : std::fill(vResult.begin(), vResult.begin() + m_nX, oOpts.outOfRangeVal);
689 1 : return true;
690 : }
691 4 : return false;
692 : }
693 :
694 : /// Mask cells based on angle intersection to the right of the observer.
695 : ///
696 : /// @param vResult Result raaster line.
697 : /// @param nLine Line number.
698 : /// @return True when all cells have been masked.
699 29287 : bool ViewshedExecutor::maskAngleRight(std::vector<double> &vResult, int nLine)
700 : {
701 29287 : int lineLength = static_cast<int>(vResult.size());
702 :
703 54 : auto clamp = [this, lineLength](int x)
704 36 : { return (x <= m_nX || x >= lineLength) ? INVALID_ISECT : x; };
705 :
706 29287 : if (oOpts.startAngle == oOpts.endAngle)
707 29267 : return false;
708 :
709 20 : if (nLine != m_nY)
710 : {
711 : int startAngleX =
712 18 : clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
713 18 : int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
714 :
715 : // If neither X intersect is in the quadrant and a ray in the quadrant isn't
716 : // between start and stop, fill it all and return true. If it is in between
717 : // start and stop, we're done.
718 18 : if (invalid(startAngleX) && invalid(endAngleX))
719 : {
720 : // Choose a test angle in quadrant I or IV depending on the line.
721 10 : double testAngle = nLine < m_nY ? m_testAngle[1] : m_testAngle[4];
722 10 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
723 : {
724 0 : std::fill(vResult.begin() + m_nX + 1, vResult.end(),
725 0 : oOpts.outOfRangeVal);
726 10 : return true;
727 : }
728 10 : return false;
729 : }
730 :
731 8 : if (nLine > m_nY)
732 0 : std::swap(startAngleX, endAngleX);
733 8 : if (invalid(endAngleX))
734 0 : endAngleX = lineLength - 1;
735 8 : if (invalid(startAngleX))
736 8 : startAngleX = m_nX + 1;
737 8 : if (startAngleX <= endAngleX)
738 : {
739 8 : std::fill(vResult.begin() + m_nX + 1, vResult.begin() + startAngleX,
740 8 : oOpts.outOfRangeVal);
741 8 : std::fill(vResult.begin() + endAngleX + 1, vResult.end(),
742 8 : oOpts.outOfRangeVal);
743 : }
744 : else
745 : {
746 0 : std::fill(vResult.begin() + endAngleX + 1,
747 0 : vResult.begin() + startAngleX, oOpts.outOfRangeVal);
748 : }
749 : }
750 : // nLine == m_nY
751 2 : else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 0))
752 : {
753 1 : std::fill(vResult.begin() + m_nX + 1, vResult.end(),
754 1 : oOpts.outOfRangeVal);
755 1 : return true;
756 : }
757 9 : return false;
758 : }
759 :
760 : /// Perform angle and min/max masking to the left of the observer.
761 : ///
762 : /// @param vResult Raster line to mask.
763 : /// @param ll Min/max line limits.
764 : /// @param nLine Line number.
765 29287 : void ViewshedExecutor::maskLineLeft(std::vector<double> &vResult,
766 : const LineLimits &ll, int nLine)
767 : {
768 : // If we've already masked with angles everything, just return.
769 29287 : if (maskAngleLeft(vResult, nLine))
770 11 : return;
771 :
772 : // Mask cells from the left edge to the left limit.
773 29276 : std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
774 : // Mask cells from the left min to the observer.
775 29276 : if (ll.leftMin < m_nX)
776 3 : std::fill(vResult.begin() + ll.leftMin, vResult.begin() + m_nX,
777 3 : oOpts.outOfRangeVal);
778 : }
779 :
780 : /// Perform angle and min/max masking to the right of the observer.
781 : ///
782 : /// @param vResult Raster line to mask.
783 : /// @param ll Min/max line limits.
784 : /// @param nLine Line number.
785 29287 : void ViewshedExecutor::maskLineRight(std::vector<double> &vResult,
786 : const LineLimits &ll, int nLine)
787 : {
788 : // If we've already masked with angles everything, just return.
789 29287 : if (maskAngleRight(vResult, nLine))
790 1 : return;
791 :
792 : // Mask cells from the observer to right min.
793 29286 : std::fill(vResult.begin() + m_nX + 1, vResult.begin() + ll.rightMin,
794 29286 : oOpts.outOfRangeVal);
795 : // Mask cells from the right limit to the right edge.
796 :
797 : //
798 : //ABELL - Changed from ll.right + 1
799 : //
800 29286 : if (ll.right <= static_cast<int>(vResult.size()))
801 29286 : std::fill(vResult.begin() + ll.right, vResult.end(),
802 29286 : oOpts.outOfRangeVal);
803 : }
804 :
805 : /// Process the part of the first line to the right of the observer.
806 : ///
807 : /// @param ll Line limits
808 : /// @param sdCalc True when doing standard deviation calcuation.
809 : /// @param lines Raster lines to process.
810 249 : void ViewshedExecutor::processFirstLineRight(const LineLimits &ll, Lines &lines,
811 : bool sdCalc)
812 : {
813 249 : int iStart = m_nX + 1;
814 249 : int iEnd = ll.right;
815 :
816 : // If start is to the right of end, everything is taken care of by left processing.
817 249 : if (iStart >= iEnd)
818 : {
819 12 : maskLineRight(lines.result, ll, m_nY);
820 12 : return;
821 : }
822 :
823 237 : iStart = oCurExtent.clampX(iStart);
824 :
825 : // If the start cell is next to the observer, just mark it visible.
826 237 : if (iStart - 1 == m_nX || iStart == oCurExtent.xStart)
827 : {
828 237 : double dfZ = lines.cur[iStart];
829 237 : if (oOpts.outputMode == OutputMode::Normal)
830 : {
831 220 : lines.result[iStart] = oOpts.visibleVal;
832 220 : if (sdCalc)
833 4 : if (lines.sd[iStart] > 1)
834 4 : lines.cur[iStart] =
835 4 : m_dfZObserver; // Use some minimum value instead?
836 : }
837 : else
838 17 : setOutputNormal(lines, iStart, dfZ);
839 237 : iStart++;
840 : }
841 :
842 : // Go from the observer to the right, calculating Z as we go.
843 10956 : for (int iPixel = iStart; iPixel < iEnd; iPixel++)
844 : {
845 10719 : int nXOffset = std::abs(iPixel - m_nX);
846 10719 : double dfZ = CalcHeightLine(nXOffset, lines.cur[iPixel - 1]);
847 10719 : if (!sdCalc)
848 10651 : setOutputNormal(lines, iPixel, dfZ);
849 : else
850 68 : setOutputSd(lines, iPixel, dfZ);
851 : }
852 :
853 237 : maskLineRight(lines.result, ll, m_nY);
854 : }
855 :
856 : /// Process a line to the left of the observer.
857 : ///
858 : /// @param nYOffset Offset of the line being processed from the observer
859 : /// @param ll Line limits
860 : /// @param lines Raster lines to process.
861 : /// @param sdCalc standard deviation calculation indicator.
862 29038 : void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
863 : Lines &lines, bool sdCalc)
864 : {
865 29038 : int iStart = m_nX - 1;
866 29038 : int iEnd = ll.left - 1;
867 29038 : int nLine = m_nY + nYOffset;
868 :
869 : // If start to the left of end, everything is taken care of by processing right.
870 29038 : if (iStart <= iEnd)
871 : {
872 2971 : maskLineLeft(lines.result, ll, nLine);
873 2971 : return;
874 : }
875 26067 : iStart = oCurExtent.clampX(iStart);
876 :
877 : // If the observer is to the right of the raster, mark the first cell to the left as
878 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
879 : // with the out of range assignment at the end.
880 26067 : if (iStart == oCurExtent.xStop - 1)
881 : {
882 6 : if (oOpts.outputMode == OutputMode::Normal)
883 0 : lines.result[iStart] = oOpts.visibleVal;
884 : else
885 6 : setOutputNormal(lines, iStart, lines.cur[iStart]);
886 6 : iStart--;
887 : }
888 :
889 : // Go from the observer to the left, calculating Z as we go.
890 26067 : nYOffset = std::abs(nYOffset);
891 1473580 : for (int iPixel = iStart; iPixel > iEnd; iPixel--)
892 : {
893 1447510 : int nXOffset = std::abs(iPixel - m_nX);
894 : double dfZ;
895 1447510 : if (nXOffset == nYOffset)
896 15856 : dfZ = CalcHeightLine(nYOffset, lines.cur[iPixel],
897 15856 : lines.prev[iPixel + 1]);
898 : else
899 1431650 : dfZ = oZcalc(nXOffset, nYOffset, lines.cur[iPixel + 1],
900 1431650 : lines.prev[iPixel], lines.prev[iPixel + 1]);
901 1447510 : if (!sdCalc)
902 1440410 : setOutputNormal(lines, iPixel, dfZ);
903 : else
904 7103 : setOutputSd(lines, iPixel, dfZ);
905 : }
906 :
907 26067 : maskLineLeft(lines.result, ll, nLine);
908 : }
909 :
910 : /// Process a line to the right of the observer.
911 : ///
912 : /// @param nYOffset Offset of the line being processed from the observer
913 : /// @param ll Line limits
914 : /// @param lines Raster lines to process.
915 : /// @param sdCalc standard deviation calculation indicator.
916 29038 : void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
917 : Lines &lines, bool sdCalc)
918 : {
919 29038 : int iStart = m_nX + 1;
920 29038 : int iEnd = ll.right;
921 29038 : int nLine = m_nY + nYOffset;
922 :
923 : // If start is to the right of end, everything is taken care of by processing left.
924 29038 : if (iStart >= iEnd)
925 : {
926 44 : maskLineRight(lines.result, ll, nLine);
927 44 : return;
928 : }
929 28994 : iStart = oCurExtent.clampX(iStart);
930 :
931 : // If the observer is to the left of the raster, mark the first cell to the right as
932 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
933 : // with the out of range assignment at the end.
934 28994 : if (iStart == 0)
935 : {
936 6 : if (oOpts.outputMode == OutputMode::Normal)
937 0 : lines.result[iStart] = oOpts.visibleVal;
938 : else
939 6 : setOutputNormal(lines, 0, lines.cur[0]);
940 6 : iStart++;
941 : }
942 :
943 : // Go from the observer to the right, calculating Z as we go.
944 28994 : nYOffset = std::abs(nYOffset);
945 1531140 : for (int iPixel = iStart; iPixel < iEnd; iPixel++)
946 : {
947 1502150 : int nXOffset = std::abs(iPixel - m_nX);
948 : double dfZ;
949 1502150 : if (nXOffset == nYOffset)
950 : {
951 16339 : if (sdCalc && nXOffset == 1)
952 : {
953 4 : lines.result[iPixel] = oOpts.visibleVal;
954 4 : if (lines.sd[iPixel] > 1)
955 4 : lines.cur[iPixel] = m_dfZObserver;
956 4 : continue;
957 : }
958 16335 : dfZ = CalcHeightLine(nYOffset, lines.cur[iPixel],
959 16335 : lines.prev[iPixel - 1]);
960 : }
961 : else
962 1485810 : dfZ = oZcalc(nXOffset, nYOffset, lines.cur[iPixel - 1],
963 1485810 : lines.prev[iPixel], lines.prev[iPixel - 1]);
964 1502150 : if (!sdCalc)
965 1495050 : setOutputNormal(lines, iPixel, dfZ);
966 : else
967 7099 : setOutputSd(lines, iPixel, dfZ);
968 : }
969 :
970 28994 : maskLineRight(lines.result, ll, nLine);
971 : }
972 :
973 : /// Apply angular/distance mask to the initial X position. Assumes m_nX is in the raster.
974 : /// @param vResult Raster line on which to apply mask.
975 : /// @param ll Line limits.
976 : /// @param nLine Line number.
977 : /// @return True if the initial X position was masked.
978 28869 : bool ViewshedExecutor::maskInitial(std::vector<double> &vResult,
979 : const LineLimits &ll, int nLine)
980 : {
981 : // Mask min/max.
982 28869 : if (ll.left >= ll.right || ll.leftMin != ll.rightMin)
983 : {
984 7 : vResult[m_nX] = oOpts.outOfRangeVal;
985 7 : return true;
986 : }
987 :
988 28862 : if (!oOpts.angleMasking())
989 28844 : return false;
990 :
991 18 : if (nLine < m_nY)
992 : {
993 13 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI / 2))
994 : {
995 0 : vResult[m_nX] = oOpts.outOfRangeVal;
996 0 : return true;
997 : }
998 : }
999 5 : else if (nLine > m_nY)
1000 : {
1001 5 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 3 * M_PI / 2))
1002 : {
1003 0 : vResult[m_nX] = oOpts.outOfRangeVal;
1004 0 : return true;
1005 : }
1006 : }
1007 18 : return false;
1008 : }
1009 :
1010 : /// Process a line above or below the observer.
1011 : ///
1012 : /// @param nLine Line number being processed.
1013 : /// @param lines Raster lines to process.
1014 : /// @return True on success, false otherwise.
1015 28881 : bool ViewshedExecutor::processLine(int nLine, Lines &lines)
1016 : {
1017 28881 : int nYOffset = nLine - m_nY;
1018 :
1019 28881 : if (!readLine(nLine, lines))
1020 0 : return false;
1021 :
1022 : // Adjust height of the read line.
1023 28881 : LineLimits ll = adjustHeight(nYOffset, lines);
1024 :
1025 57762 : std::vector<double> savedLine;
1026 28881 : if (sdMode())
1027 157 : savedLine = lines.cur;
1028 :
1029 87114 : auto process = [this, nYOffset, &ll, &lines](bool sdCalc)
1030 : {
1031 58076 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
1032 29038 : pQueue->SubmitJob([&]()
1033 29038 : { processLineLeft(nYOffset, ll, lines, sdCalc); });
1034 29038 : pQueue->SubmitJob([&]()
1035 29038 : { processLineRight(nYOffset, ll, lines, sdCalc); });
1036 29038 : pQueue->WaitCompletion();
1037 29038 : };
1038 :
1039 28881 : bool masked = false;
1040 : // Handle initial position on the line.
1041 28881 : if (oCurExtent.containsX(m_nX))
1042 : {
1043 28869 : masked = maskInitial(lines.result, ll, nLine);
1044 28869 : if (!masked)
1045 : {
1046 28862 : double dfZ = CalcHeightLine(std::abs(nYOffset), lines.cur[m_nX],
1047 28862 : lines.prev[m_nX]);
1048 28862 : setOutputNormal(lines, m_nX, dfZ);
1049 : }
1050 : }
1051 :
1052 28881 : process(false);
1053 :
1054 : // Process standard deviation mode
1055 28881 : if (sdMode())
1056 : {
1057 157 : lines.prev = std::move(lines.prevTmp);
1058 157 : lines.prevTmp = std::move(lines.cur);
1059 157 : lines.cur = std::move(savedLine);
1060 : // Handle initial position on the line.
1061 157 : if (!masked && oCurExtent.containsX(m_nX))
1062 : {
1063 157 : if (std::abs(nYOffset) == 1)
1064 : {
1065 8 : lines.result[m_nX] = oOpts.visibleVal;
1066 8 : if (lines.sd[m_nX] > 1)
1067 3 : lines.cur[m_nX] = m_dfZObserver;
1068 : }
1069 : else
1070 : {
1071 :
1072 149 : double dfZ = CalcHeightLine(std::abs(nYOffset), lines.cur[m_nX],
1073 149 : lines.prev[m_nX]);
1074 149 : setOutputSd(lines, m_nX, dfZ);
1075 : }
1076 : }
1077 157 : process(true);
1078 157 : lines.prev = std::move(lines.prevTmp);
1079 157 : lines.prevTmp = lines.cur;
1080 : }
1081 : else
1082 28724 : lines.prev = lines.cur;
1083 :
1084 28881 : if (oOpts.pitchMasking())
1085 18 : applyPitchMask(lines.result, lines.pitchMask);
1086 28881 : if (!writeLine(nLine, lines.result))
1087 0 : return false;
1088 :
1089 28881 : return oProgress.lineComplete();
1090 : }
1091 :
1092 : // Calculate the ray angle from the origin to middle of the top or bottom
1093 : // of each quadrant.
1094 2 : void ViewshedExecutor::calcTestAngles()
1095 : {
1096 : // Quadrant 1.
1097 : {
1098 2 : int ysize = m_nY + 1;
1099 2 : int xsize = oCurExtent.xStop - m_nX;
1100 2 : m_testAngle[1] = atan2(ysize, xsize / 2.0);
1101 : }
1102 :
1103 : // Quadrant 2.
1104 : {
1105 2 : int ysize = m_nY + 1;
1106 2 : int xsize = m_nX + 1;
1107 2 : m_testAngle[2] = atan2(ysize, -xsize / 2.0);
1108 : }
1109 :
1110 : // Quadrant 3.
1111 : {
1112 2 : int ysize = oCurExtent.yStop - m_nY;
1113 2 : int xsize = m_nX + 1;
1114 2 : m_testAngle[3] = atan2(-ysize, -xsize / 2.0);
1115 : }
1116 :
1117 : // Quadrant 4.
1118 : {
1119 2 : int ysize = oCurExtent.yStop - m_nY;
1120 2 : int xsize = oCurExtent.xStop - m_nX;
1121 2 : m_testAngle[4] = atan2(-ysize, xsize / 2.0);
1122 : }
1123 :
1124 : // Adjust range to [0, 2 * M_PI)
1125 10 : for (int i = 1; i <= 4; ++i)
1126 8 : if (m_testAngle[i] < 0)
1127 4 : m_testAngle[i] += (2 * M_PI);
1128 2 : }
1129 :
1130 : /// Run the viewshed computation
1131 : /// @return Success as true or false.
1132 244 : bool ViewshedExecutor::run()
1133 : {
1134 : // If we're doing angular masking, calculate the test angles used later.
1135 244 : if (oOpts.angleMasking())
1136 2 : calcTestAngles();
1137 :
1138 488 : Lines firstLine(oCurExtent.xSize());
1139 244 : if (oOpts.pitchMasking())
1140 2 : firstLine.pitchMask.resize(oOutExtent.xSize(),
1141 4 : std::numeric_limits<double>::quiet_NaN());
1142 244 : if (sdMode())
1143 9 : firstLine.sd.resize(oOutExtent.xSize());
1144 :
1145 244 : m_dfHeightAdjFactor = calcHeightAdjFactor();
1146 :
1147 244 : if (!processFirstLine(firstLine))
1148 0 : return false;
1149 :
1150 244 : if (oOpts.cellMode == CellMode::Edge)
1151 244 : oZcalc = doEdge;
1152 0 : else if (oOpts.cellMode == CellMode::Diagonal)
1153 0 : oZcalc = doDiagonal;
1154 0 : else if (oOpts.cellMode == CellMode::Min)
1155 0 : oZcalc = doMin;
1156 0 : else if (oOpts.cellMode == CellMode::Max)
1157 0 : oZcalc = doMax;
1158 :
1159 : // scan upwards
1160 244 : int yStart = oCurExtent.clampY(m_nY);
1161 244 : std::atomic<bool> err(false);
1162 244 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
1163 244 : pQueue->SubmitJob(
1164 244 : [&]()
1165 : {
1166 488 : Lines lines(oCurExtent.xSize());
1167 244 : lines.prev = firstLine.prev;
1168 244 : lines.prevTmp = firstLine.prevTmp;
1169 244 : if (oOpts.pitchMasking())
1170 2 : lines.pitchMask.resize(
1171 2 : oOutExtent.xSize(),
1172 4 : std::numeric_limits<double>::quiet_NaN());
1173 244 : if (sdMode())
1174 9 : lines.sd.resize(oOutExtent.xSize());
1175 :
1176 13605 : for (int nLine = yStart - 1; nLine >= oCurExtent.yStart && !err;
1177 : nLine--)
1178 : {
1179 13361 : if (!processLine(nLine, lines))
1180 0 : err = true;
1181 13361 : if (oOpts.pitchMasking())
1182 9 : std::fill(lines.pitchMask.begin(), lines.pitchMask.end(),
1183 18 : std::numeric_limits<double>::quiet_NaN());
1184 : }
1185 244 : });
1186 :
1187 : // scan downwards
1188 244 : pQueue->SubmitJob(
1189 244 : [&]()
1190 : {
1191 488 : Lines lines(oCurExtent.xSize());
1192 244 : lines.prev = firstLine.prev;
1193 244 : lines.prevTmp = firstLine.prevTmp;
1194 244 : if (oOpts.pitchMasking())
1195 2 : lines.pitchMask.resize(
1196 2 : oOutExtent.xSize(),
1197 4 : std::numeric_limits<double>::quiet_NaN());
1198 244 : if (sdMode())
1199 9 : lines.sd.resize(oOutExtent.xSize());
1200 :
1201 15764 : for (int nLine = yStart + 1; nLine < oCurExtent.yStop && !err;
1202 : nLine++)
1203 : {
1204 15520 : if (!processLine(nLine, lines))
1205 0 : err = true;
1206 15520 : if (oOpts.pitchMasking())
1207 9 : std::fill(lines.pitchMask.begin(), lines.pitchMask.end(),
1208 18 : std::numeric_limits<double>::quiet_NaN());
1209 : }
1210 244 : });
1211 244 : return true;
1212 : }
1213 :
1214 : } // namespace viewshed
1215 : } // namespace gdal
1216 :
1217 : // cppcheck-suppress-end knownConditionTrueFalse
|
