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 79714 : double CalcHeightLine(int nDistance, double Za)
53 : {
54 79714 : nDistance = std::abs(nDistance);
55 79714 : assert(nDistance != 1);
56 79714 : 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 2822740 : double CalcHeightEdge(int i, int j, double Za, double Zb)
73 : {
74 2822740 : assert(i != j);
75 2822740 : 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 2818240 : double doEdge(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
86 : double dfLastPrev)
87 : {
88 2818240 : if (nXOffset >= nYOffset)
89 1167790 : return CalcHeightEdge(nYOffset, nXOffset, dfLastPrev, dfThisPrev);
90 : else
91 1650450 : return CalcHeightEdge(nXOffset, nYOffset, dfLastPrev, dfLast);
92 : }
93 :
94 0 : double doMin(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
95 : double dfLastPrev)
96 : {
97 0 : double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
98 : double dfDiagonal =
99 0 : doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
100 0 : return std::min(dfEdge, dfDiagonal);
101 : }
102 :
103 0 : double doMax(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
104 : double dfLastPrev)
105 : {
106 0 : double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
107 : double dfDiagonal =
108 0 : doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
109 0 : return std::max(dfEdge, dfDiagonal);
110 : }
111 :
112 : } // unnamed namespace
113 :
114 : /// Constructor - the viewshed algorithm executor
115 : /// @param srcBand Source raster band
116 : /// @param dstBand Destination raster band
117 : /// @param nX X position of observer
118 : /// @param nY Y position of observer
119 : /// @param outExtent Extent of output raster (relative to input)
120 : /// @param curExtent Extent of active raster.
121 : /// @param opts Configuration options.
122 : /// @param progress Reference to the progress tracker.
123 : /// @param emitWarningIfNoData Whether a warning must be emitted if an input
124 : /// pixel is at the nodata value.
125 235 : ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
126 : GDALRasterBand &dstBand, int nX, int nY,
127 : const Window &outExtent,
128 : const Window &curExtent, const Options &opts,
129 235 : Progress &progress, bool emitWarningIfNoData)
130 : : m_pool(4), m_srcBand(srcBand), m_dstBand(dstBand),
131 : m_emitWarningIfNoData(emitWarningIfNoData), oOutExtent(outExtent),
132 235 : oCurExtent(curExtent), m_nX(nX - oOutExtent.xStart), m_nY(nY),
133 : oOpts(opts), oProgress(progress),
134 235 : m_dfMinDistance2(opts.minDistance * opts.minDistance),
135 235 : m_dfMaxDistance2(opts.maxDistance * opts.maxDistance)
136 : {
137 235 : if (m_dfMaxDistance2 == 0)
138 232 : m_dfMaxDistance2 = std::numeric_limits<double>::max();
139 235 : if (opts.lowPitch != -90.0)
140 1 : m_lowTanPitch = std::tan(oOpts.lowPitch * (2 * M_PI / 360.0));
141 235 : if (opts.highPitch != 90.0)
142 1 : m_highTanPitch = std::tan(oOpts.highPitch * (2 * M_PI / 360.0));
143 235 : m_srcBand.GetDataset()->GetGeoTransform(m_gt);
144 235 : int hasNoData = false;
145 235 : m_noDataValue = m_srcBand.GetNoDataValue(&hasNoData);
146 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 2888150 : void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
180 : double dfZ)
181 : {
182 2888150 : if (oOpts.outputMode != OutputMode::Normal)
183 : {
184 28621 : double adjustment = dfZ - dfCellVal;
185 28621 : if (adjustment > 0)
186 15582 : dfResult += adjustment;
187 : }
188 : else
189 2859530 : dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
190 : : oOpts.visibleVal;
191 2888150 : dfCellVal = std::max(dfCellVal, dfZ);
192 2889660 : }
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 28911 : bool ViewshedExecutor::readLine(int nLine, double *data)
200 : {
201 57768 : std::lock_guard g(iMutex);
202 :
203 28932 : if (GDALRasterIO(&m_srcBand, GF_Read, oOutExtent.xStart, nLine,
204 : oOutExtent.xSize(), 1, data, oOutExtent.xSize(), 1,
205 28857 : 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 28857 : 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 28859 : bool ViewshedExecutor::writeLine(int nLine, std::vector<double> &vResult)
222 : {
223 : // GDALRasterIO isn't thread-safe.
224 57673 : std::lock_guard g(oMutex);
225 :
226 57577 : if (GDALRasterIO(&m_dstBand, GF_Write, 0, nLine - oOutExtent.yStart,
227 28804 : oOutExtent.xSize(), 1, vResult.data(), oOutExtent.xSize(),
228 28814 : 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 28814 : 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 : /// @return Processing limits of the line based on min/max distance.
246 28845 : LineLimits ViewshedExecutor::adjustHeight(int nYOffset,
247 : std::vector<double> &vThisLineVal,
248 : std::vector<double> &vPitchMaskVal)
249 : {
250 28845 : LineLimits ll(0, m_nX + 1, m_nX + 1, oCurExtent.xSize());
251 :
252 : // Find the starting point in the raster (m_nX may be outside)
253 28748 : int nXStart = oCurExtent.clampX(m_nX);
254 :
255 5803540 : const auto CheckNoData = [this](double val)
256 : {
257 5800360 : if (!m_hasFoundNoData &&
258 2904350 : ((m_hasNoData && val == m_noDataValue) || std::isnan(val)))
259 : {
260 0 : m_hasFoundNoData = true;
261 0 : if (m_emitWarningIfNoData)
262 : {
263 0 : CPLError(CE_Warning, CPLE_AppDefined,
264 : "Nodata value found in input DEM. Output will be "
265 : "likely incorrect");
266 : }
267 : }
268 2924730 : };
269 :
270 : // If there is a height adjustment factor other than zero or a max distance,
271 : // calculate the adjusted height of the cell, stopping if we've exceeded the max
272 : // distance.
273 27339 : if (static_cast<bool>(m_dfHeightAdjFactor) || oOpts.pitchMasking() ||
274 56120 : m_dfMaxDistance2 > 0 || m_dfMinDistance2 > 0)
275 : {
276 : // Hoist invariants from the loops.
277 28733 : const double dfLineX = m_gt[2] * nYOffset;
278 28733 : const double dfLineY = m_gt[5] * nYOffset;
279 :
280 : // Go left
281 28731 : double *pdfHeight = vThisLineVal.data() + nXStart;
282 1483710 : for (int nXOffset = nXStart - m_nX; nXOffset >= -m_nX;
283 : nXOffset--, pdfHeight--)
284 : {
285 1454740 : double dfX = m_gt[1] * nXOffset + dfLineX;
286 1453690 : double dfY = m_gt[4] * nXOffset + dfLineY;
287 1458340 : double dfR2 = dfX * dfX + dfY * dfY;
288 :
289 1458340 : if (dfR2 < m_dfMinDistance2)
290 6 : ll.leftMin--;
291 1458330 : else if (dfR2 > m_dfMaxDistance2)
292 : {
293 26 : ll.left = nXOffset + m_nX + 1;
294 26 : break;
295 : }
296 :
297 1458310 : CheckNoData(*pdfHeight);
298 1454880 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
299 1454880 : if (oOpts.pitchMasking())
300 75 : calcPitchMask(*pdfHeight, std::sqrt(dfR2),
301 75 : vPitchMaskVal[m_nX + nXOffset]);
302 : }
303 :
304 : // Go right.
305 29005 : pdfHeight = vThisLineVal.data() + nXStart + 1;
306 1495350 : for (int nXOffset = nXStart - m_nX + 1;
307 1495350 : nXOffset < oCurExtent.xSize() - m_nX; nXOffset++, pdfHeight++)
308 : {
309 1464740 : double dfX = m_gt[1] * nXOffset + dfLineX;
310 1463440 : double dfY = m_gt[4] * nXOffset + dfLineY;
311 1470180 : double dfR2 = dfX * dfX + dfY * dfY;
312 :
313 1470180 : if (dfR2 < m_dfMinDistance2)
314 3 : ll.rightMin++;
315 1470180 : else if (dfR2 > m_dfMaxDistance2)
316 : {
317 26 : ll.right = nXOffset + m_nX;
318 26 : break;
319 : }
320 :
321 1470160 : CheckNoData(*pdfHeight);
322 1463420 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
323 1463420 : if (oOpts.pitchMasking())
324 175 : calcPitchMask(*pdfHeight, std::sqrt(dfR2),
325 175 : vPitchMaskVal[m_nX + nXOffset]);
326 : }
327 : }
328 : else
329 : {
330 : // No curvature adjustment. Just normalize for the observer height.
331 48 : double *pdfHeight = vThisLineVal.data();
332 0 : for (int i = 0; i < oCurExtent.xSize(); ++i)
333 : {
334 1 : CheckNoData(*pdfHeight);
335 0 : *pdfHeight -= m_dfZObserver;
336 0 : pdfHeight++;
337 : }
338 : }
339 28935 : return ll;
340 : }
341 :
342 250 : void ViewshedExecutor::calcPitchMask(double dfZ, double dfDist, double &maskVal)
343 : {
344 250 : if (oOpts.lowPitchMasking())
345 : {
346 25 : double dfZMask = dfDist * m_lowTanPitch;
347 25 : double adjustment = dfZMask - dfZ;
348 25 : if (adjustment > 0)
349 : {
350 24 : maskVal =
351 24 : (oOpts.outputMode == OutputMode::Normal ? oOpts.outOfRangeVal
352 : : adjustment);
353 24 : return;
354 : }
355 : }
356 226 : if (oOpts.highPitchMasking())
357 : {
358 225 : double dfZMask = dfDist * m_highTanPitch;
359 225 : if (dfZ > dfZMask)
360 7 : maskVal = oOpts.outOfRangeVal;
361 : }
362 : }
363 :
364 : /// Process the first line (the one with the Y coordinate the same as the observer).
365 : ///
366 : /// @param vLastLineVal Vector in which to store the read line. Becomes the last line
367 : /// in further processing.
368 : /// @return True on success, false otherwise.
369 235 : bool ViewshedExecutor::processFirstLine(std::vector<double> &vLastLineVal)
370 : {
371 235 : int nLine = oOutExtent.clampY(m_nY);
372 235 : int nYOffset = nLine - m_nY;
373 :
374 470 : std::vector<double> vResult(oOutExtent.xSize());
375 470 : std::vector<double> vThisLineVal(oOutExtent.xSize());
376 470 : std::vector<double> vPitchMaskVal;
377 235 : if (oOpts.pitchMasking())
378 2 : vPitchMaskVal.resize(oOutExtent.xSize(),
379 4 : std::numeric_limits<double>::quiet_NaN());
380 :
381 235 : if (!readLine(nLine, vThisLineVal.data()))
382 0 : return false;
383 :
384 : // If the observer is outside of the raster, take the specified value as the Z height,
385 : // otherwise, take it as an offset from the raster height at that location.
386 235 : m_dfZObserver = oOpts.observer.z;
387 235 : if (oCurExtent.containsX(m_nX))
388 : {
389 229 : m_dfZObserver += vThisLineVal[m_nX];
390 229 : if (oOpts.outputMode == OutputMode::Normal)
391 212 : vResult[m_nX] = oOpts.visibleVal;
392 : }
393 235 : m_dfHeightAdjFactor = calcHeightAdjFactor();
394 :
395 : // In DEM mode the base is the pre-adjustment value. In ground mode the base is zero.
396 235 : if (oOpts.outputMode == OutputMode::DEM)
397 16 : vResult = vThisLineVal;
398 :
399 235 : LineLimits ll = adjustHeight(nYOffset, vThisLineVal, vPitchMaskVal);
400 235 : if (oCurExtent.containsX(m_nX) && ll.leftMin != ll.rightMin)
401 1 : vResult[m_nX] = oOpts.outOfRangeVal;
402 :
403 235 : if (!oCurExtent.containsY(m_nY))
404 4 : processFirstLineTopOrBottom(ll, vResult, vThisLineVal);
405 : else
406 : {
407 462 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
408 231 : pQueue->SubmitJob([&]()
409 231 : { processFirstLineLeft(ll, vResult, vThisLineVal); });
410 231 : pQueue->SubmitJob(
411 231 : [&]() { processFirstLineRight(ll, vResult, vThisLineVal); });
412 231 : pQueue->WaitCompletion();
413 : }
414 :
415 : // Make the current line the last line.
416 235 : vLastLineVal = std::move(vThisLineVal);
417 :
418 235 : if (oOpts.pitchMasking())
419 2 : applyPitchMask(vResult, vPitchMaskVal);
420 235 : if (!writeLine(nLine, vResult))
421 0 : return false;
422 :
423 235 : return oProgress.lineComplete();
424 : }
425 :
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 (vPitchMaskVal[i] == oOpts.outOfRangeVal)
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 vResult Result line.
444 : /// @param vThisLineVal Heights of the cells in the target line
445 4 : void ViewshedExecutor::processFirstLineTopOrBottom(
446 : const LineLimits &ll, std::vector<double> &vResult,
447 : std::vector<double> &vThisLineVal)
448 : {
449 4 : double *pResult = vResult.data() + ll.left;
450 4 : double *pThis = vThisLineVal.data() + ll.left;
451 24 : for (int iPixel = ll.left; iPixel < ll.right; ++iPixel, ++pResult, pThis++)
452 : {
453 20 : if (oOpts.outputMode == OutputMode::Normal)
454 0 : *pResult = oOpts.visibleVal;
455 : else
456 20 : setOutput(*pResult, *pThis, *pThis);
457 : }
458 :
459 4 : std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
460 4 : std::fill(vResult.begin() + ll.right, vResult.begin() + oCurExtent.xStop,
461 4 : oOpts.outOfRangeVal);
462 4 : }
463 :
464 : /// Process the part of the first line to the left of the observer.
465 : ///
466 : /// @param ll Line limits for masking.
467 : /// @param vResult Vector in which to store the visibility/height results.
468 : /// @param vThisLineVal Height of each cell in the line being processed.
469 231 : void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll,
470 : std::vector<double> &vResult,
471 : std::vector<double> &vThisLineVal)
472 : {
473 231 : int iEnd = ll.left - 1;
474 231 : int iStart = m_nX - 1; // One left of the observer.
475 :
476 : // If end is to the right of start, everything is taken care of by right processing.
477 231 : if (iEnd >= iStart)
478 30 : return;
479 :
480 201 : iStart = oCurExtent.clampX(iStart);
481 :
482 201 : double *pThis = vThisLineVal.data() + iStart;
483 :
484 : // If the start cell is next to the observer, just mark it visible.
485 201 : if (iStart + 1 == m_nX || iStart + 1 == oCurExtent.xStop)
486 : {
487 201 : double dfZ = *pThis;
488 201 : if (oOpts.outputMode == OutputMode::Normal)
489 190 : vResult[iStart] = oOpts.visibleVal;
490 : else
491 11 : setOutput(vResult[iStart], *pThis, dfZ);
492 201 : iStart--;
493 201 : pThis--;
494 : }
495 :
496 : // Go from the observer to the left, calculating Z as we go.
497 10357 : for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--)
498 : {
499 10156 : int nXOffset = std::abs(iPixel - m_nX);
500 10156 : double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
501 10156 : setOutput(vResult[iPixel], *pThis, dfZ);
502 : }
503 :
504 201 : maskLineLeft(vResult, ll, m_nY);
505 : }
506 :
507 : /// Mask cells based on angle intersection to the left of the observer.
508 : ///
509 : /// @param vResult Result raaster line.
510 : /// @param nLine Line number.
511 : /// @return True when all cells have been masked.
512 25946 : bool ViewshedExecutor::maskAngleLeft(std::vector<double> &vResult, int nLine)
513 : {
514 38 : auto clamp = [this](int x)
515 20 : { return (x < 0 || x >= m_nX) ? INVALID_ISECT : x; };
516 :
517 25946 : if (!oOpts.angleMasking())
518 25906 : return false;
519 :
520 11 : if (nLine != m_nY)
521 : {
522 : int startAngleX =
523 10 : clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
524 10 : int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
525 : // If neither X intersect is in the quadrant and a ray in the quadrant isn't
526 : // between start and stop, fill it all and return true. If it is in between
527 : // start and stop, we're done.
528 10 : if (invalid(startAngleX) && invalid(endAngleX))
529 : {
530 : // Choose a test angle in quadrant II or III depending on the line.
531 7 : double testAngle = nLine < m_nY ? m_testAngle[2] : m_testAngle[3];
532 7 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
533 : {
534 2 : std::fill(vResult.begin(), vResult.begin() + m_nX,
535 2 : oOpts.outOfRangeVal);
536 7 : return true;
537 : }
538 5 : return false;
539 : }
540 3 : if (nLine > m_nY)
541 0 : std::swap(startAngleX, endAngleX);
542 3 : if (invalid(startAngleX))
543 3 : startAngleX = 0;
544 3 : if (invalid(endAngleX))
545 0 : endAngleX = m_nX - 1;
546 3 : if (startAngleX <= endAngleX)
547 : {
548 3 : std::fill(vResult.begin(), vResult.begin() + startAngleX,
549 3 : oOpts.outOfRangeVal);
550 3 : std::fill(vResult.begin() + endAngleX + 1, vResult.begin() + m_nX,
551 3 : oOpts.outOfRangeVal);
552 : }
553 : else
554 : {
555 0 : std::fill(vResult.begin() + endAngleX + 1,
556 0 : vResult.begin() + startAngleX, oOpts.outOfRangeVal);
557 : }
558 : }
559 : // nLine == m_nY
560 1 : else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI))
561 : {
562 0 : std::fill(vResult.begin(), vResult.begin() + m_nX, oOpts.outOfRangeVal);
563 0 : return true;
564 : }
565 4 : return false;
566 : }
567 :
568 : /// Mask cells based on angle intersection to the right of the observer.
569 : ///
570 : /// @param vResult Result raaster line.
571 : /// @param nLine Line number.
572 : /// @return True when all cells have been masked.
573 28911 : bool ViewshedExecutor::maskAngleRight(std::vector<double> &vResult, int nLine)
574 : {
575 28911 : int lineLength = static_cast<int>(vResult.size());
576 :
577 54 : auto clamp = [this, lineLength](int x)
578 36 : { return (x <= m_nX || x >= lineLength) ? INVALID_ISECT : x; };
579 :
580 28912 : if (oOpts.startAngle == oOpts.endAngle)
581 28894 : return false;
582 :
583 18 : if (nLine != m_nY)
584 : {
585 : int startAngleX =
586 18 : clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
587 18 : int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
588 :
589 : // If neither X intersect is in the quadrant and a ray in the quadrant isn't
590 : // between start and stop, fill it all and return true. If it is in between
591 : // start and stop, we're done.
592 18 : if (invalid(startAngleX) && invalid(endAngleX))
593 : {
594 : // Choose a test angle in quadrant I or IV depending on the line.
595 10 : double testAngle = nLine < m_nY ? m_testAngle[1] : m_testAngle[4];
596 10 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
597 : {
598 0 : std::fill(vResult.begin() + m_nX + 1, vResult.end(),
599 0 : oOpts.outOfRangeVal);
600 10 : return true;
601 : }
602 10 : return false;
603 : }
604 :
605 8 : if (nLine > m_nY)
606 0 : std::swap(startAngleX, endAngleX);
607 8 : if (invalid(endAngleX))
608 0 : endAngleX = lineLength - 1;
609 8 : if (invalid(startAngleX))
610 8 : startAngleX = m_nX + 1;
611 8 : if (startAngleX <= endAngleX)
612 : {
613 8 : std::fill(vResult.begin() + m_nX + 1, vResult.begin() + startAngleX,
614 8 : oOpts.outOfRangeVal);
615 8 : std::fill(vResult.begin() + endAngleX + 1, vResult.end(),
616 8 : oOpts.outOfRangeVal);
617 : }
618 : else
619 : {
620 0 : std::fill(vResult.begin() + endAngleX + 1,
621 0 : vResult.begin() + startAngleX, oOpts.outOfRangeVal);
622 : }
623 : }
624 : // nLine == m_nY
625 0 : else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 0))
626 : {
627 1 : std::fill(vResult.begin() + m_nX + 1, vResult.end(),
628 1 : oOpts.outOfRangeVal);
629 1 : return true;
630 : }
631 9 : return false;
632 : }
633 :
634 : /// Perform angle and min/max masking to the left of the observer.
635 : ///
636 : /// @param vResult Raster line to mask.
637 : /// @param ll Min/max line limits.
638 : /// @param nLine Line number.
639 25959 : void ViewshedExecutor::maskLineLeft(std::vector<double> &vResult,
640 : const LineLimits &ll, int nLine)
641 : {
642 : // If we've already masked with angles everything, just return.
643 25959 : if (maskAngleLeft(vResult, nLine))
644 2 : return;
645 :
646 : // Mask cells from the left edge to the left limit.
647 25927 : std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
648 : // Mask cells from the left min to the observer.
649 25874 : if (ll.leftMin < m_nX)
650 3 : std::fill(vResult.begin() + ll.leftMin, vResult.begin() + m_nX,
651 3 : oOpts.outOfRangeVal);
652 : }
653 :
654 : /// Perform angle and min/max masking to the right of the observer.
655 : ///
656 : /// @param vResult Raster line to mask.
657 : /// @param ll Min/max line limits.
658 : /// @param nLine Line number.
659 28922 : void ViewshedExecutor::maskLineRight(std::vector<double> &vResult,
660 : const LineLimits &ll, int nLine)
661 : {
662 : // If we've already masked with angles everything, just return.
663 28922 : if (maskAngleRight(vResult, nLine))
664 1 : return;
665 :
666 : // Mask cells from the observer to right min.
667 28887 : std::fill(vResult.begin() + m_nX + 1, vResult.begin() + ll.rightMin,
668 28919 : oOpts.outOfRangeVal);
669 : // Mask cells from the right limit to the right edge.
670 28864 : if (ll.right + 1 < static_cast<int>(vResult.size()))
671 8 : std::fill(vResult.begin() + ll.right + 1, vResult.end(),
672 8 : oOpts.outOfRangeVal);
673 : }
674 :
675 : /// Process the part of the first line to the right of the observer.
676 : ///
677 : /// @param ll Line limits
678 : /// @param vResult Vector in which to store the visibility/height results.
679 : /// @param vThisLineVal Height of each cell in the line being processed.
680 231 : void ViewshedExecutor::processFirstLineRight(const LineLimits &ll,
681 : std::vector<double> &vResult,
682 : std::vector<double> &vThisLineVal)
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 2 : return;
690 :
691 229 : iStart = oCurExtent.clampX(iStart);
692 :
693 229 : double *pThis = vThisLineVal.data() + iStart;
694 :
695 : // If the start cell is next to the observer, just mark it visible.
696 229 : if (iStart - 1 == m_nX || iStart == oCurExtent.xStart)
697 : {
698 229 : double dfZ = *pThis;
699 229 : if (oOpts.outputMode == OutputMode::Normal)
700 212 : vResult[iStart] = oOpts.visibleVal;
701 : else
702 17 : setOutput(vResult[iStart], *pThis, dfZ);
703 229 : iStart++;
704 229 : pThis++;
705 : }
706 :
707 : // Go from the observer to the right, calculating Z as we go.
708 10812 : for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++)
709 : {
710 10583 : int nXOffset = std::abs(iPixel - m_nX);
711 10583 : double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
712 10583 : setOutput(vResult[iPixel], *pThis, dfZ);
713 : }
714 :
715 229 : maskLineRight(vResult, ll, m_nY);
716 : }
717 :
718 : /// Process a line to the left of the observer.
719 : ///
720 : /// @param nYOffset Offset of the line being processed from the observer
721 : /// @param ll Line limits
722 : /// @param vResult Vector in which to store the visibility/height results.
723 : /// @param vThisLineVal Height of each cell in the line being processed.
724 : /// @param vLastLineVal Observable height of each cell in the previous line processed.
725 28633 : void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
726 : std::vector<double> &vResult,
727 : std::vector<double> &vThisLineVal,
728 : std::vector<double> &vLastLineVal)
729 : {
730 28633 : int iStart = m_nX - 1;
731 28633 : int iEnd = ll.left - 1;
732 28633 : int nLine = m_nY + nYOffset;
733 : // If start to the left of end, everything is taken care of by processing right.
734 28633 : if (iStart <= iEnd)
735 2924 : return;
736 25709 : iStart = oCurExtent.clampX(iStart);
737 :
738 25664 : nYOffset = std::abs(nYOffset);
739 25664 : double *pThis = vThisLineVal.data() + iStart;
740 25647 : double *pLast = vLastLineVal.data() + iStart;
741 :
742 : // If the observer is to the right of the raster, mark the first cell to the left as
743 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
744 : // with the out of range assignment at the end.
745 :
746 25674 : if (iStart == oCurExtent.xStop - 1)
747 : {
748 6 : if (oOpts.outputMode == OutputMode::Normal)
749 0 : vResult[iStart] = oOpts.visibleVal;
750 : else
751 6 : setOutput(vResult[iStart], *pThis, *pThis);
752 6 : iStart--;
753 6 : pThis--;
754 6 : pLast--;
755 : }
756 :
757 : // Go from the observer to the left, calculating Z as we go.
758 1435330 : for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--, pLast--)
759 : {
760 1409560 : int nXOffset = std::abs(iPixel - m_nX);
761 : double dfZ;
762 1409560 : if (nXOffset == nYOffset)
763 : {
764 15647 : if (nXOffset == 1)
765 375 : dfZ = *pThis;
766 : else
767 15272 : dfZ = CalcHeightLine(nXOffset, *(pLast + 1));
768 : }
769 : else
770 : dfZ =
771 1393920 : oZcalc(nXOffset, nYOffset, *(pThis + 1), *pLast, *(pLast + 1));
772 1411470 : setOutput(vResult[iPixel], *pThis, dfZ);
773 : }
774 :
775 25769 : maskLineLeft(vResult, ll, nLine);
776 : }
777 :
778 : /// Process a line to the right of the observer.
779 : ///
780 : /// @param nYOffset Offset of the line being processed from the observer
781 : /// @param ll Line limits
782 : /// @param vResult Vector in which to store the visibility/height results.
783 : /// @param vThisLineVal Height of each cell in the line being processed.
784 : /// @param vLastLineVal Observable height of each cell in the previous line processed.
785 28661 : void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
786 : std::vector<double> &vResult,
787 : std::vector<double> &vThisLineVal,
788 : std::vector<double> &vLastLineVal)
789 : {
790 28661 : int iStart = m_nX + 1;
791 28661 : int iEnd = ll.right;
792 28661 : int nLine = m_nY + nYOffset;
793 :
794 : // If start is to the right of end, everything is taken care of by processing left.
795 28661 : if (iStart >= iEnd)
796 12 : return;
797 28649 : iStart = oCurExtent.clampX(iStart);
798 :
799 28647 : nYOffset = std::abs(nYOffset);
800 28647 : double *pThis = vThisLineVal.data() + iStart;
801 28625 : double *pLast = vLastLineVal.data() + iStart;
802 :
803 : // If the observer is to the left of the raster, mark the first cell to the right as
804 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
805 : // with the out of range assignment at the end.
806 28631 : if (iStart == 0)
807 : {
808 6 : if (oOpts.outputMode == OutputMode::Normal)
809 0 : vResult[iStart] = oOpts.visibleVal;
810 : else
811 6 : setOutput(vResult[0], *pThis, *pThis);
812 6 : iStart++;
813 6 : pThis++;
814 6 : pLast++;
815 : }
816 :
817 : // Go from the observer to the right, calculating Z as we go.
818 1486320 : for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++, pLast++)
819 : {
820 1457640 : int nXOffset = std::abs(iPixel - m_nX);
821 : double dfZ;
822 1457640 : if (nXOffset == nYOffset)
823 : {
824 16122 : if (nXOffset == 1)
825 416 : dfZ = *pThis;
826 : else
827 15706 : dfZ = CalcHeightLine(nXOffset, *(pLast - 1));
828 : }
829 : else
830 : dfZ =
831 1441520 : oZcalc(nXOffset, nYOffset, *(pThis - 1), *pLast, *(pLast - 1));
832 1459210 : setOutput(vResult[iPixel], *pThis, dfZ);
833 : }
834 :
835 28684 : maskLineRight(vResult, ll, nLine);
836 : }
837 :
838 : /// Apply angular mask to the initial X position. Assumes m_nX is in the raster.
839 : /// @param vResult Raster line on which to apply mask.
840 : /// @param nLine Line number.
841 28546 : void ViewshedExecutor::maskInitial(std::vector<double> &vResult, int nLine)
842 : {
843 28546 : if (!oOpts.angleMasking())
844 28482 : return;
845 :
846 21 : if (nLine < m_nY)
847 : {
848 13 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI / 2))
849 0 : vResult[m_nX] = oOpts.outOfRangeVal;
850 : }
851 8 : else if (nLine > m_nY)
852 : {
853 5 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 3 * M_PI / 2))
854 0 : vResult[m_nX] = oOpts.outOfRangeVal;
855 : }
856 : }
857 :
858 : /// Process a line above or below the observer.
859 : ///
860 : /// @param nLine Line number being processed.
861 : /// @param vLastLineVal Vector in which to store the read line. Becomes the last line
862 : /// in further processing.
863 : /// @return True on success, false otherwise.
864 28719 : bool ViewshedExecutor::processLine(int nLine, std::vector<double> &vLastLineVal)
865 : {
866 28719 : int nYOffset = nLine - m_nY;
867 57431 : std::vector<double> vResult(oOutExtent.xSize());
868 57418 : std::vector<double> vThisLineVal(oOutExtent.xSize());
869 57432 : std::vector<double> vPitchMaskVal;
870 28685 : if (oOpts.pitchMasking())
871 18 : vPitchMaskVal.resize(oOutExtent.xSize(),
872 36 : std::numeric_limits<double>::quiet_NaN());
873 :
874 28684 : if (!readLine(nLine, vThisLineVal.data()))
875 0 : return false;
876 :
877 : // In DEM mode the base is the input DEM value.
878 : // In ground mode the base is zero.
879 28632 : if (oOpts.outputMode == OutputMode::DEM)
880 163 : vResult = vThisLineVal;
881 :
882 : // Adjust height of the read line.
883 28632 : LineLimits ll = adjustHeight(nYOffset, vThisLineVal, vPitchMaskVal);
884 :
885 : // Handle the initial position on the line.
886 28702 : if (oCurExtent.containsX(m_nX))
887 : {
888 28562 : if (ll.left < ll.right && ll.leftMin == ll.rightMin)
889 : {
890 : double dfZ;
891 28553 : if (std::abs(nYOffset) == 1)
892 414 : dfZ = vThisLineVal[m_nX];
893 : else
894 28139 : dfZ = CalcHeightLine(nYOffset, vLastLineVal[m_nX]);
895 28549 : setOutput(vResult[m_nX], vThisLineVal[m_nX], dfZ);
896 : }
897 : else
898 9 : vResult[m_nX] = oOpts.outOfRangeVal;
899 :
900 28503 : maskInitial(vResult, nLine);
901 : }
902 :
903 : // process left half then right half of line
904 57296 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
905 28513 : pQueue->SubmitJob(
906 28629 : [&]() {
907 28629 : processLineLeft(nYOffset, ll, vResult, vThisLineVal, vLastLineVal);
908 28605 : });
909 28615 : pQueue->SubmitJob(
910 28674 : [&]() {
911 28674 : processLineRight(nYOffset, ll, vResult, vThisLineVal, vLastLineVal);
912 28659 : });
913 28634 : pQueue->WaitCompletion();
914 : // Make the current line the last line.
915 28674 : vLastLineVal = std::move(vThisLineVal);
916 :
917 28600 : if (oOpts.pitchMasking())
918 18 : applyPitchMask(vResult, vPitchMaskVal);
919 28510 : if (!writeLine(nLine, vResult))
920 0 : return false;
921 :
922 28575 : return oProgress.lineComplete();
923 : }
924 :
925 : // Calculate the ray angle from the origin to middle of the top or bottom
926 : // of each quadrant.
927 2 : void ViewshedExecutor::calcTestAngles()
928 : {
929 : // Quadrant 1.
930 : {
931 2 : int ysize = m_nY + 1;
932 2 : int xsize = oCurExtent.xStop - m_nX;
933 2 : m_testAngle[1] = atan2(ysize, xsize / 2.0);
934 : }
935 :
936 : // Quadrant 2.
937 : {
938 2 : int ysize = m_nY + 1;
939 2 : int xsize = m_nX + 1;
940 2 : m_testAngle[2] = atan2(ysize, -xsize / 2.0);
941 : }
942 :
943 : // Quadrant 3.
944 : {
945 2 : int ysize = oCurExtent.yStop - m_nY;
946 2 : int xsize = m_nX + 1;
947 2 : m_testAngle[3] = atan2(-ysize, -xsize / 2.0);
948 : }
949 :
950 : // Quadrant 4.
951 : {
952 2 : int ysize = oCurExtent.yStop - m_nY;
953 2 : int xsize = oCurExtent.xStop - m_nX;
954 2 : m_testAngle[4] = atan2(-ysize, xsize / 2.0);
955 : }
956 :
957 : // Adjust range to [0, 2 * M_PI)
958 10 : for (int i = 1; i <= 4; ++i)
959 8 : if (m_testAngle[i] < 0)
960 4 : m_testAngle[i] += (2 * M_PI);
961 2 : }
962 :
963 : /// Run the viewshed computation
964 : /// @return Success as true or false.
965 235 : bool ViewshedExecutor::run()
966 : {
967 : // If we're doing angular masking, calculate the test angles used later.
968 235 : if (oOpts.angleMasking())
969 2 : calcTestAngles();
970 :
971 470 : std::vector<double> vFirstLineVal(oCurExtent.xSize());
972 235 : if (!processFirstLine(vFirstLineVal))
973 0 : return false;
974 :
975 235 : if (oOpts.cellMode == CellMode::Edge)
976 235 : oZcalc = doEdge;
977 0 : else if (oOpts.cellMode == CellMode::Diagonal)
978 0 : oZcalc = doDiagonal;
979 0 : else if (oOpts.cellMode == CellMode::Min)
980 0 : oZcalc = doMin;
981 0 : else if (oOpts.cellMode == CellMode::Max)
982 0 : oZcalc = doMax;
983 :
984 : // scan upwards
985 235 : int yStart = oCurExtent.clampY(m_nY);
986 235 : std::atomic<bool> err(false);
987 235 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
988 235 : pQueue->SubmitJob(
989 235 : [&]()
990 : {
991 471 : std::vector<double> vLastLineVal = vFirstLineVal;
992 :
993 13517 : for (int nLine = yStart - 1; nLine >= oCurExtent.yStart && !err;
994 : nLine--)
995 13280 : if (!processLine(nLine, vLastLineVal))
996 0 : err = true;
997 235 : });
998 :
999 : // scan downwards
1000 235 : pQueue->SubmitJob(
1001 235 : [&]()
1002 : {
1003 469 : std::vector<double> vLastLineVal = vFirstLineVal;
1004 :
1005 15672 : for (int nLine = yStart + 1; nLine < oCurExtent.yStop && !err;
1006 : nLine++)
1007 15437 : if (!processLine(nLine, vLastLineVal))
1008 0 : err = true;
1009 235 : });
1010 235 : return true;
1011 : }
1012 :
1013 : } // namespace viewshed
1014 : } // namespace gdal
|
