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
|
