Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Viewshed Generation
4 : * Purpose: Core algorithm implementation for viewshed generation.
5 : * Author: Tamas Szekeres, szekerest@gmail.com
6 : *
7 : * (c) 2024 info@hobu.co
8 : *
9 : ******************************************************************************
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include <algorithm>
15 : #include <atomic>
16 : #include <cassert>
17 : #include <cmath>
18 : #include <limits>
19 :
20 : #include "viewshed_executor.h"
21 : #include "progress.h"
22 : #include "util.h"
23 :
24 : namespace gdal
25 : {
26 : namespace viewshed
27 : {
28 :
29 : namespace
30 : {
31 :
32 : /// Determines whether a value is a valid intersection coordinate.
33 : /// @param i Value to test.
34 : /// @return True if the value doesn't represent an invalid intersection.
35 78 : bool valid(int i)
36 : {
37 78 : return i != INVALID_ISECT;
38 : }
39 :
40 : /// Determines whether a value is an invalid intersection coordinate.
41 : /// @param i Value to test.
42 : /// @return True if the value represents an invalid intersection.
43 78 : bool invalid(int i)
44 : {
45 78 : return !valid(i);
46 : }
47 :
48 : /// Calculate the height at nDistance units along a line through the origin given the height
49 : /// at nDistance - 1 units along the line.
50 : /// \param nDistance Distance along the line for the target point.
51 : /// \param Za Height at the line one unit previous to the target point.
52 79803 : double CalcHeightLine(int nDistance, double Za)
53 : {
54 79803 : nDistance = std::abs(nDistance);
55 79803 : assert(nDistance != 1);
56 79803 : 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 2842140 : double CalcHeightEdge(int i, int j, double Za, double Zb)
73 : {
74 2842140 : assert(i != j);
75 2842140 : 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 2842300 : double doEdge(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
86 : double dfLastPrev)
87 : {
88 2842300 : if (nXOffset >= nYOffset)
89 1173250 : return CalcHeightEdge(nYOffset, nXOffset, dfLastPrev, dfThisPrev);
90 : else
91 1669050 : 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 235 : m_hasNoData = hasNoData;
147 235 : }
148 :
149 : // calculate the height adjustment factor.
150 235 : double ViewshedExecutor::calcHeightAdjFactor()
151 : {
152 470 : std::lock_guard g(oMutex);
153 :
154 : const OGRSpatialReference *poDstSRS =
155 235 : m_dstBand.GetDataset()->GetSpatialRef();
156 :
157 235 : 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 223 : 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 2914590 : void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
180 : double dfZ)
181 : {
182 2914590 : if (oOpts.outputMode != OutputMode::Normal)
183 : {
184 28786 : double adjustment = dfZ - dfCellVal;
185 28786 : if (adjustment > 0)
186 15686 : dfResult += adjustment;
187 : }
188 : else
189 2885800 : dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
190 : : oOpts.visibleVal;
191 2914590 : dfCellVal = std::max(dfCellVal, dfZ);
192 2915300 : }
193 :
194 : /// Read a line of raster data.
195 : ///
196 : /// @param nLine Line number to read.
197 : /// @param data Pointer to location in which to store data.
198 : /// @return Success or failure.
199 28912 : bool ViewshedExecutor::readLine(int nLine, double *data)
200 : {
201 57786 : std::lock_guard g(iMutex);
202 :
203 28939 : if (GDALRasterIO(&m_srcBand, GF_Read, oOutExtent.xStart, nLine,
204 : oOutExtent.xSize(), 1, data, oOutExtent.xSize(), 1,
205 28874 : 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 28874 : 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 28800 : bool ViewshedExecutor::writeLine(int nLine, std::vector<double> &vResult)
222 : {
223 : // GDALRasterIO isn't thread-safe.
224 57622 : std::lock_guard g(oMutex);
225 :
226 57586 : if (GDALRasterIO(&m_dstBand, GF_Write, 0, nLine - oOutExtent.yStart,
227 28794 : oOutExtent.xSize(), 1, vResult.data(), oOutExtent.xSize(),
228 28822 : 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 28822 : return true;
237 : }
238 :
239 : /// Adjust the height of the line of data by the observer height and the curvature of the
240 : /// earth.
241 : ///
242 : /// @param nYOffset Y offset of the line being adjusted.
243 : /// @param vThisLineVal Line height data.
244 : /// @param vPitchMaskVal Pitch masking line.
245 : /// @param vResult Initialized Result vector.
246 : /// @return Processing limits of the line based on min/max distance.
247 28870 : LineLimits ViewshedExecutor::adjustHeight(int nYOffset,
248 : std::vector<double> &vThisLineVal,
249 : const std::vector<double> &vResult,
250 : std::vector<double> &vPitchMaskVal)
251 : {
252 28870 : LineLimits ll(0, m_nX + 1, m_nX + 1, oCurExtent.xSize());
253 :
254 : // Find the starting point in the raster (m_nX may be outside)
255 28784 : int nXStart = oCurExtent.clampX(m_nX);
256 :
257 5834030 : const auto CheckNoData = [this](double val)
258 : {
259 5828730 : if (!m_hasFoundNoData &&
260 2918520 : ((m_hasNoData && val == m_noDataValue) || std::isnan(val)))
261 : {
262 0 : m_hasFoundNoData = true;
263 0 : if (m_emitWarningIfNoData)
264 : {
265 0 : CPLError(CE_Warning, CPLE_AppDefined,
266 : "Nodata value found in input DEM. Output will be "
267 : "likely incorrect");
268 : }
269 : }
270 2939000 : };
271 :
272 : // If there is a height adjustment factor other than zero or a max distance,
273 : // calculate the adjusted height of the cell, stopping if we've exceeded the max
274 : // distance.
275 27405 : if (static_cast<bool>(m_dfHeightAdjFactor) || oOpts.pitchMasking() ||
276 56240 : m_dfMaxDistance2 > 0 || m_dfMinDistance2 > 0)
277 : {
278 : // Hoist invariants from the loops.
279 28796 : const double dfLineX = m_gt[2] * nYOffset;
280 28798 : const double dfLineY = m_gt[5] * nYOffset;
281 :
282 : // Go left
283 28792 : double *pdfHeight = vThisLineVal.data() + nXStart;
284 1474460 : for (int nXOffset = nXStart - m_nX; nXOffset >= -m_nX;
285 : nXOffset--, pdfHeight--)
286 : {
287 1446220 : double dfX = m_gt[1] * nXOffset + dfLineX;
288 1446480 : double dfY = m_gt[4] * nXOffset + dfLineY;
289 1451410 : double dfR2 = dfX * dfX + dfY * dfY;
290 :
291 1451410 : if (dfR2 < m_dfMinDistance2)
292 6 : ll.leftMin--;
293 1451410 : else if (dfR2 > m_dfMaxDistance2)
294 : {
295 26 : ll.left = nXOffset + m_nX + 1;
296 26 : break;
297 : }
298 :
299 1451390 : CheckNoData(*pdfHeight);
300 1443030 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
301 1443030 : if (oOpts.pitchMasking())
302 75 : calcPitchMask(*pdfHeight, std::sqrt(dfR2),
303 75 : vResult[m_nX + nXOffset],
304 75 : vPitchMaskVal[m_nX + nXOffset]);
305 : }
306 :
307 : // Go right.
308 28261 : pdfHeight = vThisLineVal.data() + nXStart + 1;
309 1516310 : for (int nXOffset = nXStart - m_nX + 1;
310 1516310 : nXOffset < oCurExtent.xSize() - m_nX; nXOffset++, pdfHeight++)
311 : {
312 1487780 : double dfX = m_gt[1] * nXOffset + dfLineX;
313 1485240 : double dfY = m_gt[4] * nXOffset + dfLineY;
314 1488300 : double dfR2 = dfX * dfX + dfY * dfY;
315 :
316 1488300 : if (dfR2 < m_dfMinDistance2)
317 3 : ll.rightMin++;
318 1488300 : else if (dfR2 > m_dfMaxDistance2)
319 : {
320 26 : ll.right = nXOffset + m_nX;
321 26 : break;
322 : }
323 :
324 1488270 : CheckNoData(*pdfHeight);
325 1486570 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
326 1486570 : if (oOpts.pitchMasking())
327 175 : calcPitchMask(*pdfHeight, std::sqrt(dfR2),
328 175 : vResult[m_nX + nXOffset],
329 175 : vPitchMaskVal[m_nX + nXOffset]);
330 : }
331 : }
332 : else
333 : {
334 : // No curvature adjustment. Just normalize for the observer height.
335 39 : double *pdfHeight = vThisLineVal.data();
336 0 : for (int i = 0; i < oCurExtent.xSize(); ++i)
337 : {
338 4 : CheckNoData(*pdfHeight);
339 0 : *pdfHeight -= m_dfZObserver;
340 0 : pdfHeight++;
341 : }
342 : }
343 28944 : return ll;
344 : }
345 :
346 : /// Calculate the pitch masking value to apply after running the viewshed algorithm.
347 : ///
348 : /// @param dfZ Adjusted input height.
349 : /// @param dfDist Distance from observer to cell.
350 : /// @param dfResult Result value to which adjustment may be added.
351 : /// @param maskVal Output mask value.
352 250 : void ViewshedExecutor::calcPitchMask(double dfZ, double dfDist, double dfResult,
353 : double &maskVal)
354 : {
355 250 : if (oOpts.lowPitchMasking())
356 : {
357 25 : double dfZMask = dfDist * m_lowTanPitch;
358 25 : double adjustment = dfZMask - dfZ;
359 25 : if (adjustment > 0)
360 : {
361 48 : maskVal = (oOpts.outputMode == OutputMode::Normal
362 24 : ? std::numeric_limits<double>::infinity()
363 : : adjustment + dfResult);
364 24 : return;
365 : }
366 : }
367 226 : if (oOpts.highPitchMasking())
368 : {
369 225 : double dfZMask = dfDist * m_highTanPitch;
370 225 : if (dfZ > dfZMask)
371 7 : maskVal = std::numeric_limits<double>::infinity();
372 : }
373 : }
374 :
375 : /// Process the first line (the one with the Y coordinate the same as the observer).
376 : ///
377 : /// @param vLastLineVal Vector in which to store the read line. Becomes the last line
378 : /// in further processing.
379 : /// @return True on success, false otherwise.
380 235 : bool ViewshedExecutor::processFirstLine(std::vector<double> &vLastLineVal)
381 : {
382 235 : int nLine = oOutExtent.clampY(m_nY);
383 235 : int nYOffset = nLine - m_nY;
384 :
385 470 : std::vector<double> vResult(oOutExtent.xSize());
386 470 : std::vector<double> vThisLineVal(oOutExtent.xSize());
387 470 : std::vector<double> vPitchMaskVal;
388 235 : if (oOpts.pitchMasking())
389 2 : vPitchMaskVal.resize(oOutExtent.xSize(),
390 4 : std::numeric_limits<double>::quiet_NaN());
391 :
392 235 : if (!readLine(nLine, vThisLineVal.data()))
393 0 : return false;
394 :
395 : // If the observer is outside of the raster, take the specified value as the Z height,
396 : // otherwise, take it as an offset from the raster height at that location.
397 235 : m_dfZObserver = oOpts.observer.z;
398 235 : if (oCurExtent.containsX(m_nX))
399 : {
400 229 : m_dfZObserver += vThisLineVal[m_nX];
401 229 : if (oOpts.outputMode == OutputMode::Normal)
402 212 : vResult[m_nX] = oOpts.visibleVal;
403 : }
404 235 : m_dfHeightAdjFactor = calcHeightAdjFactor();
405 :
406 : // In DEM mode the base is the pre-adjustment value. In ground mode the base is zero.
407 235 : if (oOpts.outputMode == OutputMode::DEM)
408 16 : vResult = vThisLineVal;
409 :
410 : LineLimits ll =
411 235 : adjustHeight(nYOffset, vThisLineVal, vResult, vPitchMaskVal);
412 235 : if (oCurExtent.containsX(m_nX) && ll.leftMin != ll.rightMin)
413 1 : vResult[m_nX] = oOpts.outOfRangeVal;
414 :
415 235 : if (!oCurExtent.containsY(m_nY))
416 4 : processFirstLineTopOrBottom(ll, vResult, vThisLineVal);
417 : else
418 : {
419 462 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
420 231 : pQueue->SubmitJob([&]()
421 231 : { processFirstLineLeft(ll, vResult, vThisLineVal); });
422 231 : pQueue->SubmitJob(
423 231 : [&]() { processFirstLineRight(ll, vResult, vThisLineVal); });
424 231 : pQueue->WaitCompletion();
425 : }
426 :
427 : // Make the current line the last line.
428 235 : vLastLineVal = std::move(vThisLineVal);
429 :
430 235 : if (oOpts.pitchMasking())
431 2 : applyPitchMask(vResult, vPitchMaskVal);
432 235 : if (!writeLine(nLine, vResult))
433 0 : return false;
434 :
435 235 : return oProgress.lineComplete();
436 : }
437 :
438 : /// Set the pitch masked value into the result vector when applicable.
439 : ///
440 : /// @param vResult Result vector.
441 : /// @param vPitchMaskVal Pitch mask values (nan is no masking, inf is out-of-range, else
442 : /// actual value).
443 20 : void ViewshedExecutor::applyPitchMask(std::vector<double> &vResult,
444 : const std::vector<double> &vPitchMaskVal)
445 : {
446 270 : for (size_t i = 0; i < vResult.size(); ++i)
447 : {
448 250 : if (std::isnan(vPitchMaskVal[i]))
449 219 : continue;
450 31 : if (std::isinf(vPitchMaskVal[i]))
451 7 : vResult[i] = oOpts.outOfRangeVal;
452 : else
453 24 : vResult[i] = vPitchMaskVal[i];
454 : }
455 20 : }
456 :
457 : /// If the observer is above or below the raster, set all cells in the first line near the
458 : /// observer as observable provided they're in range. Mark cells out of range as such.
459 : /// @param ll Line limits for processing.
460 : /// @param vResult Result line.
461 : /// @param vThisLineVal Heights of the cells in the target line
462 4 : void ViewshedExecutor::processFirstLineTopOrBottom(
463 : const LineLimits &ll, std::vector<double> &vResult,
464 : std::vector<double> &vThisLineVal)
465 : {
466 4 : double *pResult = vResult.data() + ll.left;
467 4 : double *pThis = vThisLineVal.data() + ll.left;
468 24 : for (int iPixel = ll.left; iPixel < ll.right; ++iPixel, ++pResult, pThis++)
469 : {
470 20 : if (oOpts.outputMode == OutputMode::Normal)
471 0 : *pResult = oOpts.visibleVal;
472 : else
473 20 : setOutput(*pResult, *pThis, *pThis);
474 : }
475 :
476 4 : std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
477 4 : std::fill(vResult.begin() + ll.right, vResult.begin() + oCurExtent.xStop,
478 4 : oOpts.outOfRangeVal);
479 4 : }
480 :
481 : /// Process the part of the first line to the left of the observer.
482 : ///
483 : /// @param ll Line limits for masking.
484 : /// @param vResult Vector in which to store the visibility/height results.
485 : /// @param vThisLineVal Height of each cell in the line being processed.
486 231 : void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll,
487 : std::vector<double> &vResult,
488 : std::vector<double> &vThisLineVal)
489 : {
490 231 : int iEnd = ll.left - 1;
491 231 : int iStart = m_nX - 1; // One left of the observer.
492 :
493 : // If end is to the right of start, everything is taken care of by right processing.
494 231 : if (iEnd >= iStart)
495 30 : return;
496 :
497 201 : iStart = oCurExtent.clampX(iStart);
498 :
499 201 : double *pThis = vThisLineVal.data() + iStart;
500 :
501 : // If the start cell is next to the observer, just mark it visible.
502 201 : if (iStart + 1 == m_nX || iStart + 1 == oCurExtent.xStop)
503 : {
504 201 : double dfZ = *pThis;
505 201 : if (oOpts.outputMode == OutputMode::Normal)
506 190 : vResult[iStart] = oOpts.visibleVal;
507 : else
508 11 : setOutput(vResult[iStart], *pThis, dfZ);
509 201 : iStart--;
510 201 : pThis--;
511 : }
512 :
513 : // Go from the observer to the left, calculating Z as we go.
514 10357 : for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--)
515 : {
516 10156 : int nXOffset = std::abs(iPixel - m_nX);
517 10156 : double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
518 10156 : setOutput(vResult[iPixel], *pThis, dfZ);
519 : }
520 :
521 201 : maskLineLeft(vResult, ll, m_nY);
522 : }
523 :
524 : /// Mask cells based on angle intersection to the left of the observer.
525 : ///
526 : /// @param vResult Result raaster line.
527 : /// @param nLine Line number.
528 : /// @return True when all cells have been masked.
529 25935 : bool ViewshedExecutor::maskAngleLeft(std::vector<double> &vResult, int nLine)
530 : {
531 38 : auto clamp = [this](int x)
532 20 : { return (x < 0 || x >= m_nX) ? INVALID_ISECT : x; };
533 :
534 25935 : if (!oOpts.angleMasking())
535 25859 : return false;
536 :
537 20 : if (nLine != m_nY)
538 : {
539 : int startAngleX =
540 10 : clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
541 10 : int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
542 : // If neither X intersect is in the quadrant and a ray in the quadrant isn't
543 : // between start and stop, fill it all and return true. If it is in between
544 : // start and stop, we're done.
545 10 : if (invalid(startAngleX) && invalid(endAngleX))
546 : {
547 : // Choose a test angle in quadrant II or III depending on the line.
548 7 : double testAngle = nLine < m_nY ? m_testAngle[2] : m_testAngle[3];
549 7 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
550 : {
551 2 : std::fill(vResult.begin(), vResult.begin() + m_nX,
552 2 : oOpts.outOfRangeVal);
553 7 : return true;
554 : }
555 5 : return false;
556 : }
557 3 : if (nLine > m_nY)
558 0 : std::swap(startAngleX, endAngleX);
559 3 : if (invalid(startAngleX))
560 3 : startAngleX = 0;
561 3 : if (invalid(endAngleX))
562 0 : endAngleX = m_nX - 1;
563 3 : if (startAngleX <= endAngleX)
564 : {
565 3 : std::fill(vResult.begin(), vResult.begin() + startAngleX,
566 3 : oOpts.outOfRangeVal);
567 3 : std::fill(vResult.begin() + endAngleX + 1, vResult.begin() + m_nX,
568 3 : oOpts.outOfRangeVal);
569 : }
570 : else
571 : {
572 0 : std::fill(vResult.begin() + endAngleX + 1,
573 0 : vResult.begin() + startAngleX, oOpts.outOfRangeVal);
574 : }
575 : }
576 : // nLine == m_nY
577 10 : else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI))
578 : {
579 0 : std::fill(vResult.begin(), vResult.begin() + m_nX, oOpts.outOfRangeVal);
580 0 : return true;
581 : }
582 4 : return false;
583 : }
584 :
585 : /// Mask cells based on angle intersection to the right of the observer.
586 : ///
587 : /// @param vResult Result raaster line.
588 : /// @param nLine Line number.
589 : /// @return True when all cells have been masked.
590 28900 : bool ViewshedExecutor::maskAngleRight(std::vector<double> &vResult, int nLine)
591 : {
592 28900 : int lineLength = static_cast<int>(vResult.size());
593 :
594 54 : auto clamp = [this, lineLength](int x)
595 36 : { return (x <= m_nX || x >= lineLength) ? INVALID_ISECT : x; };
596 :
597 28885 : if (oOpts.startAngle == oOpts.endAngle)
598 28865 : return false;
599 :
600 20 : if (nLine != m_nY)
601 : {
602 : int startAngleX =
603 18 : clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
604 18 : int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
605 :
606 : // If neither X intersect is in the quadrant and a ray in the quadrant isn't
607 : // between start and stop, fill it all and return true. If it is in between
608 : // start and stop, we're done.
609 18 : if (invalid(startAngleX) && invalid(endAngleX))
610 : {
611 : // Choose a test angle in quadrant I or IV depending on the line.
612 10 : double testAngle = nLine < m_nY ? m_testAngle[1] : m_testAngle[4];
613 10 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
614 : {
615 0 : std::fill(vResult.begin() + m_nX + 1, vResult.end(),
616 0 : oOpts.outOfRangeVal);
617 10 : return true;
618 : }
619 10 : return false;
620 : }
621 :
622 8 : if (nLine > m_nY)
623 0 : std::swap(startAngleX, endAngleX);
624 8 : if (invalid(endAngleX))
625 0 : endAngleX = lineLength - 1;
626 8 : if (invalid(startAngleX))
627 8 : startAngleX = m_nX + 1;
628 8 : if (startAngleX <= endAngleX)
629 : {
630 8 : std::fill(vResult.begin() + m_nX + 1, vResult.begin() + startAngleX,
631 8 : oOpts.outOfRangeVal);
632 8 : std::fill(vResult.begin() + endAngleX + 1, vResult.end(),
633 8 : oOpts.outOfRangeVal);
634 : }
635 : else
636 : {
637 0 : std::fill(vResult.begin() + endAngleX + 1,
638 0 : vResult.begin() + startAngleX, oOpts.outOfRangeVal);
639 : }
640 : }
641 : // nLine == m_nY
642 2 : else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 0))
643 : {
644 1 : std::fill(vResult.begin() + m_nX + 1, vResult.end(),
645 1 : oOpts.outOfRangeVal);
646 1 : return true;
647 : }
648 9 : return false;
649 : }
650 :
651 : /// Perform angle and min/max masking to the left of the observer.
652 : ///
653 : /// @param vResult Raster line to mask.
654 : /// @param ll Min/max line limits.
655 : /// @param nLine Line number.
656 25950 : void ViewshedExecutor::maskLineLeft(std::vector<double> &vResult,
657 : const LineLimits &ll, int nLine)
658 : {
659 : // If we've already masked with angles everything, just return.
660 25950 : if (maskAngleLeft(vResult, nLine))
661 2 : return;
662 :
663 : // Mask cells from the left edge to the left limit.
664 25911 : std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
665 : // Mask cells from the left min to the observer.
666 25781 : if (ll.leftMin < m_nX)
667 3 : std::fill(vResult.begin() + ll.leftMin, vResult.begin() + m_nX,
668 3 : oOpts.outOfRangeVal);
669 : }
670 :
671 : /// Perform angle and min/max masking to the right of the observer.
672 : ///
673 : /// @param vResult Raster line to mask.
674 : /// @param ll Min/max line limits.
675 : /// @param nLine Line number.
676 28926 : void ViewshedExecutor::maskLineRight(std::vector<double> &vResult,
677 : const LineLimits &ll, int nLine)
678 : {
679 : // If we've already masked with angles everything, just return.
680 28926 : if (maskAngleRight(vResult, nLine))
681 1 : return;
682 :
683 : // Mask cells from the observer to right min.
684 28870 : std::fill(vResult.begin() + m_nX + 1, vResult.begin() + ll.rightMin,
685 28894 : oOpts.outOfRangeVal);
686 : // Mask cells from the right limit to the right edge.
687 28821 : if (ll.right + 1 < static_cast<int>(vResult.size()))
688 8 : std::fill(vResult.begin() + ll.right + 1, vResult.end(),
689 8 : oOpts.outOfRangeVal);
690 : }
691 :
692 : /// Process the part of the first line to the right of the observer.
693 : ///
694 : /// @param ll Line limits
695 : /// @param vResult Vector in which to store the visibility/height results.
696 : /// @param vThisLineVal Height of each cell in the line being processed.
697 231 : void ViewshedExecutor::processFirstLineRight(const LineLimits &ll,
698 : std::vector<double> &vResult,
699 : std::vector<double> &vThisLineVal)
700 : {
701 231 : int iStart = m_nX + 1;
702 231 : int iEnd = ll.right;
703 :
704 : // If start is to the right of end, everything is taken care of by left processing.
705 231 : if (iStart >= iEnd)
706 2 : return;
707 :
708 229 : iStart = oCurExtent.clampX(iStart);
709 :
710 229 : double *pThis = vThisLineVal.data() + iStart;
711 :
712 : // If the start cell is next to the observer, just mark it visible.
713 229 : if (iStart - 1 == m_nX || iStart == oCurExtent.xStart)
714 : {
715 229 : double dfZ = *pThis;
716 229 : if (oOpts.outputMode == OutputMode::Normal)
717 212 : vResult[iStart] = oOpts.visibleVal;
718 : else
719 17 : setOutput(vResult[iStart], *pThis, dfZ);
720 229 : iStart++;
721 229 : pThis++;
722 : }
723 :
724 : // Go from the observer to the right, calculating Z as we go.
725 10770 : for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++)
726 : {
727 10541 : int nXOffset = std::abs(iPixel - m_nX);
728 10541 : double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
729 10535 : setOutput(vResult[iPixel], *pThis, dfZ);
730 : }
731 :
732 229 : maskLineRight(vResult, ll, m_nY);
733 : }
734 :
735 : /// Process a line to the left of the observer.
736 : ///
737 : /// @param nYOffset Offset of the line being processed from the observer
738 : /// @param ll Line limits
739 : /// @param vResult Vector in which to store the visibility/height results.
740 : /// @param vThisLineVal Height of each cell in the line being processed.
741 : /// @param vLastLineVal Observable height of each cell in the previous line processed.
742 28630 : void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
743 : std::vector<double> &vResult,
744 : std::vector<double> &vThisLineVal,
745 : std::vector<double> &vLastLineVal)
746 : {
747 28630 : int iStart = m_nX - 1;
748 28630 : int iEnd = ll.left - 1;
749 28630 : int nLine = m_nY + nYOffset;
750 : // If start to the left of end, everything is taken care of by processing right.
751 28630 : if (iStart <= iEnd)
752 2922 : return;
753 25708 : iStart = oCurExtent.clampX(iStart);
754 :
755 25679 : nYOffset = std::abs(nYOffset);
756 25679 : double *pThis = vThisLineVal.data() + iStart;
757 25595 : double *pLast = vLastLineVal.data() + iStart;
758 :
759 : // If the observer is to the right of the raster, mark the first cell to the left as
760 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
761 : // with the out of range assignment at the end.
762 :
763 25632 : if (iStart == oCurExtent.xStop - 1)
764 : {
765 6 : if (oOpts.outputMode == OutputMode::Normal)
766 0 : vResult[iStart] = oOpts.visibleVal;
767 : else
768 6 : setOutput(vResult[iStart], *pThis, *pThis);
769 6 : iStart--;
770 6 : pThis--;
771 6 : pLast--;
772 : }
773 :
774 : // Go from the observer to the left, calculating Z as we go.
775 1445030 : for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--, pLast--)
776 : {
777 1419260 : int nXOffset = std::abs(iPixel - m_nX);
778 : double dfZ;
779 1419260 : if (nXOffset == nYOffset)
780 : {
781 15647 : if (nXOffset == 1)
782 375 : dfZ = *pThis;
783 : else
784 15272 : dfZ = CalcHeightLine(nXOffset, *(pLast + 1));
785 : }
786 : else
787 : dfZ =
788 1403610 : oZcalc(nXOffset, nYOffset, *(pThis + 1), *pLast, *(pLast + 1));
789 1421680 : setOutput(vResult[iPixel], *pThis, dfZ);
790 : }
791 :
792 25770 : maskLineLeft(vResult, ll, nLine);
793 : }
794 :
795 : /// Process a line to the right of the observer.
796 : ///
797 : /// @param nYOffset Offset of the line being processed from the observer
798 : /// @param ll Line limits
799 : /// @param vResult Vector in which to store the visibility/height results.
800 : /// @param vThisLineVal Height of each cell in the line being processed.
801 : /// @param vLastLineVal Observable height of each cell in the previous line processed.
802 28631 : void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
803 : std::vector<double> &vResult,
804 : std::vector<double> &vThisLineVal,
805 : std::vector<double> &vLastLineVal)
806 : {
807 28631 : int iStart = m_nX + 1;
808 28631 : int iEnd = ll.right;
809 28631 : int nLine = m_nY + nYOffset;
810 :
811 : // If start is to the right of end, everything is taken care of by processing left.
812 28631 : if (iStart >= iEnd)
813 12 : return;
814 28619 : iStart = oCurExtent.clampX(iStart);
815 :
816 28631 : nYOffset = std::abs(nYOffset);
817 28631 : double *pThis = vThisLineVal.data() + iStart;
818 28584 : double *pLast = vLastLineVal.data() + iStart;
819 :
820 : // If the observer is to the left of the raster, mark the first cell to the right as
821 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
822 : // with the out of range assignment at the end.
823 28593 : if (iStart == 0)
824 : {
825 6 : if (oOpts.outputMode == OutputMode::Normal)
826 0 : vResult[iStart] = oOpts.visibleVal;
827 : else
828 6 : setOutput(vResult[0], *pThis, *pThis);
829 6 : iStart++;
830 6 : pThis++;
831 6 : pLast++;
832 : }
833 :
834 : // Go from the observer to the right, calculating Z as we go.
835 1492220 : for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++, pLast++)
836 : {
837 1463530 : int nXOffset = std::abs(iPixel - m_nX);
838 : double dfZ;
839 1463530 : if (nXOffset == nYOffset)
840 : {
841 16128 : if (nXOffset == 1)
842 415 : dfZ = *pThis;
843 : else
844 15713 : dfZ = CalcHeightLine(nXOffset, *(pLast - 1));
845 : }
846 : else
847 : dfZ =
848 1447400 : oZcalc(nXOffset, nYOffset, *(pThis - 1), *pLast, *(pLast - 1));
849 1466730 : setOutput(vResult[iPixel], *pThis, dfZ);
850 : }
851 :
852 28686 : maskLineRight(vResult, ll, nLine);
853 : }
854 :
855 : /// Apply angular mask to the initial X position. Assumes m_nX is in the raster.
856 : /// @param vResult Raster line on which to apply mask.
857 : /// @param nLine Line number.
858 28610 : void ViewshedExecutor::maskInitial(std::vector<double> &vResult, int nLine)
859 : {
860 28610 : if (!oOpts.angleMasking())
861 28538 : return;
862 :
863 41 : if (nLine < m_nY)
864 : {
865 13 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI / 2))
866 0 : vResult[m_nX] = oOpts.outOfRangeVal;
867 : }
868 28 : else if (nLine > m_nY)
869 : {
870 5 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 3 * M_PI / 2))
871 0 : vResult[m_nX] = oOpts.outOfRangeVal;
872 : }
873 : }
874 :
875 : /// Process a line above or below the observer.
876 : ///
877 : /// @param nLine Line number being processed.
878 : /// @param vLastLineVal Vector in which to store the read line. Becomes the last line
879 : /// in further processing.
880 : /// @return True on success, false otherwise.
881 28721 : bool ViewshedExecutor::processLine(int nLine, std::vector<double> &vLastLineVal)
882 : {
883 28721 : int nYOffset = nLine - m_nY;
884 57439 : std::vector<double> vResult(oOutExtent.xSize());
885 57413 : std::vector<double> vThisLineVal(oOutExtent.xSize());
886 57436 : std::vector<double> vPitchMaskVal;
887 28703 : if (oOpts.pitchMasking())
888 18 : vPitchMaskVal.resize(oOutExtent.xSize(),
889 36 : std::numeric_limits<double>::quiet_NaN());
890 :
891 28709 : if (!readLine(nLine, vThisLineVal.data()))
892 0 : return false;
893 :
894 : // In DEM mode the base is the input DEM value.
895 : // In ground mode the base is zero.
896 28658 : if (oOpts.outputMode == OutputMode::DEM)
897 163 : vResult = vThisLineVal;
898 :
899 : // Adjust height of the read line.
900 : LineLimits ll =
901 28658 : adjustHeight(nYOffset, vThisLineVal, vResult, vPitchMaskVal);
902 :
903 : // Handle the initial position on the line.
904 28690 : if (oCurExtent.containsX(m_nX))
905 : {
906 28626 : if (ll.left < ll.right && ll.leftMin == ll.rightMin)
907 : {
908 : double dfZ;
909 28613 : if (std::abs(nYOffset) == 1)
910 414 : dfZ = vThisLineVal[m_nX];
911 : else
912 28199 : dfZ = CalcHeightLine(nYOffset, vLastLineVal[m_nX]);
913 28593 : setOutput(vResult[m_nX], vThisLineVal[m_nX], dfZ);
914 : }
915 : else
916 13 : vResult[m_nX] = oOpts.outOfRangeVal;
917 :
918 28567 : maskInitial(vResult, nLine);
919 : }
920 :
921 : // process left half then right half of line
922 57346 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
923 28533 : pQueue->SubmitJob(
924 28596 : [&]() {
925 28596 : processLineLeft(nYOffset, ll, vResult, vThisLineVal, vLastLineVal);
926 28524 : });
927 28642 : pQueue->SubmitJob(
928 28654 : [&]() {
929 28654 : processLineRight(nYOffset, ll, vResult, vThisLineVal, vLastLineVal);
930 28623 : });
931 28642 : pQueue->WaitCompletion();
932 : // Make the current line the last line.
933 28672 : vLastLineVal = std::move(vThisLineVal);
934 :
935 28573 : if (oOpts.pitchMasking())
936 18 : applyPitchMask(vResult, vPitchMaskVal);
937 28579 : if (!writeLine(nLine, vResult))
938 0 : return false;
939 :
940 28597 : return oProgress.lineComplete();
941 : }
942 :
943 : // Calculate the ray angle from the origin to middle of the top or bottom
944 : // of each quadrant.
945 2 : void ViewshedExecutor::calcTestAngles()
946 : {
947 : // Quadrant 1.
948 : {
949 2 : int ysize = m_nY + 1;
950 2 : int xsize = oCurExtent.xStop - m_nX;
951 2 : m_testAngle[1] = atan2(ysize, xsize / 2.0);
952 : }
953 :
954 : // Quadrant 2.
955 : {
956 2 : int ysize = m_nY + 1;
957 2 : int xsize = m_nX + 1;
958 2 : m_testAngle[2] = atan2(ysize, -xsize / 2.0);
959 : }
960 :
961 : // Quadrant 3.
962 : {
963 2 : int ysize = oCurExtent.yStop - m_nY;
964 2 : int xsize = m_nX + 1;
965 2 : m_testAngle[3] = atan2(-ysize, -xsize / 2.0);
966 : }
967 :
968 : // Quadrant 4.
969 : {
970 2 : int ysize = oCurExtent.yStop - m_nY;
971 2 : int xsize = oCurExtent.xStop - m_nX;
972 2 : m_testAngle[4] = atan2(-ysize, xsize / 2.0);
973 : }
974 :
975 : // Adjust range to [0, 2 * M_PI)
976 10 : for (int i = 1; i <= 4; ++i)
977 8 : if (m_testAngle[i] < 0)
978 4 : m_testAngle[i] += (2 * M_PI);
979 2 : }
980 :
981 : /// Run the viewshed computation
982 : /// @return Success as true or false.
983 235 : bool ViewshedExecutor::run()
984 : {
985 : // If we're doing angular masking, calculate the test angles used later.
986 235 : if (oOpts.angleMasking())
987 2 : calcTestAngles();
988 :
989 470 : std::vector<double> vFirstLineVal(oCurExtent.xSize());
990 235 : if (!processFirstLine(vFirstLineVal))
991 0 : return false;
992 :
993 235 : if (oOpts.cellMode == CellMode::Edge)
994 235 : oZcalc = doEdge;
995 0 : else if (oOpts.cellMode == CellMode::Diagonal)
996 0 : oZcalc = doDiagonal;
997 0 : else if (oOpts.cellMode == CellMode::Min)
998 0 : oZcalc = doMin;
999 0 : else if (oOpts.cellMode == CellMode::Max)
1000 0 : oZcalc = doMax;
1001 :
1002 : // scan upwards
1003 235 : int yStart = oCurExtent.clampY(m_nY);
1004 235 : std::atomic<bool> err(false);
1005 235 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
1006 235 : pQueue->SubmitJob(
1007 235 : [&]()
1008 : {
1009 469 : std::vector<double> vLastLineVal = vFirstLineVal;
1010 :
1011 13514 : for (int nLine = yStart - 1; nLine >= oCurExtent.yStart && !err;
1012 : nLine--)
1013 13281 : if (!processLine(nLine, vLastLineVal))
1014 0 : err = true;
1015 235 : });
1016 :
1017 : // scan downwards
1018 235 : pQueue->SubmitJob(
1019 235 : [&]()
1020 : {
1021 470 : std::vector<double> vLastLineVal = vFirstLineVal;
1022 :
1023 15674 : for (int nLine = yStart + 1; nLine < oCurExtent.yStop && !err;
1024 : nLine++)
1025 15437 : if (!processLine(nLine, vLastLineVal))
1026 0 : err = true;
1027 235 : });
1028 235 : return true;
1029 : }
1030 :
1031 : } // namespace viewshed
1032 : } // namespace gdal
|
