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 79893 : double CalcHeightLine(int nDistance, double Za)
53 : {
54 79893 : nDistance = std::abs(nDistance);
55 79893 : assert(nDistance != 1);
56 79893 : 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 2802530 : double CalcHeightEdge(int i, int j, double Za, double Zb)
73 : {
74 2802530 : assert(i != j);
75 2802530 : 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 2796330 : double doEdge(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
86 : double dfLastPrev)
87 : {
88 2796330 : if (nXOffset >= nYOffset)
89 1164280 : return CalcHeightEdge(nYOffset, nXOffset, dfLastPrev, dfThisPrev);
90 : else
91 1632060 : 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_adfTransform.data());
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 2859900 : void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
180 : double dfZ)
181 : {
182 2859900 : if (oOpts.outputMode != OutputMode::Normal)
183 : {
184 28435 : dfResult += (dfZ - dfCellVal);
185 28435 : dfResult = std::max(0.0, dfResult);
186 : }
187 : else
188 2831460 : dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
189 : : oOpts.visibleVal;
190 2859820 : dfCellVal = std::max(dfCellVal, dfZ);
191 2872260 : }
192 :
193 : /// Read a line of raster data.
194 : ///
195 : /// @param nLine Line number to read.
196 : /// @param data Pointer to location in which to store data.
197 : /// @return Success or failure.
198 28952 : bool ViewshedExecutor::readLine(int nLine, double *data)
199 : {
200 57851 : std::lock_guard g(iMutex);
201 :
202 28954 : if (GDALRasterIO(&m_srcBand, GF_Read, oOutExtent.xStart, nLine,
203 : oOutExtent.xSize(), 1, data, oOutExtent.xSize(), 1,
204 28899 : GDT_Float64, 0, 0))
205 : {
206 0 : CPLError(CE_Failure, CPLE_AppDefined,
207 : "RasterIO error when reading DEM at position (%d,%d), "
208 : "size (%d,%d)",
209 0 : oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
210 0 : return false;
211 : }
212 28899 : return true;
213 : }
214 :
215 : /// Write an output line of either visibility or height data.
216 : ///
217 : /// @param nLine Line number being written.
218 : /// @param vResult Result line to write.
219 : /// @return True on success, false otherwise.
220 28805 : bool ViewshedExecutor::writeLine(int nLine, std::vector<double> &vResult)
221 : {
222 : // GDALRasterIO isn't thread-safe.
223 57706 : std::lock_guard g(oMutex);
224 :
225 57597 : if (GDALRasterIO(&m_dstBand, GF_Write, 0, nLine - oOutExtent.yStart,
226 28830 : oOutExtent.xSize(), 1, vResult.data(), oOutExtent.xSize(),
227 28901 : 1, GDT_Float64, 0, 0))
228 : {
229 0 : CPLError(CE_Failure, CPLE_AppDefined,
230 : "RasterIO error when writing target raster at position "
231 : "(%d,%d), size (%d,%d)",
232 0 : 0, nLine - oOutExtent.yStart, oOutExtent.xSize(), 1);
233 0 : return false;
234 : }
235 28901 : return true;
236 : }
237 :
238 : /// Adjust the height of the line of data by the observer height and the curvature of the
239 : /// earth.
240 : ///
241 : /// @param nYOffset Y offset of the line being adjusted.
242 : /// @param vThisLineVal Line height data.
243 : /// @return Processing limits of the line based on min/max distance.
244 28899 : LineLimits ViewshedExecutor::adjustHeight(int nYOffset,
245 : std::vector<double> &vThisLineVal)
246 : {
247 28899 : 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 28820 : int nXStart = oCurExtent.clampX(m_nX);
251 :
252 5782650 : const auto CheckNoData = [this](double val)
253 : {
254 5770360 : if (!m_hasFoundNoData &&
255 2894010 : ((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 2905200 : };
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 28845 : if (static_cast<bool>(m_dfHeightAdjFactor) || m_dfMaxDistance2 > 0 ||
271 0 : m_dfMinDistance2 > 0)
272 : {
273 : // Hoist invariants from the loops.
274 28845 : const double dfLineX = m_adfTransform[2] * nYOffset;
275 28787 : const double dfLineY = m_adfTransform[5] * nYOffset;
276 :
277 : // Go left
278 28809 : double *pdfHeight = vThisLineVal.data() + nXStart;
279 1469590 : for (int nXOffset = nXStart - m_nX; nXOffset >= -m_nX;
280 : nXOffset--, pdfHeight--)
281 : {
282 1441610 : double dfX = m_adfTransform[1] * nXOffset + dfLineX;
283 1441560 : double dfY = m_adfTransform[4] * nXOffset + dfLineY;
284 1445840 : double dfR2 = dfX * dfX + dfY * dfY;
285 :
286 1445840 : if (dfR2 < m_dfMinDistance2)
287 6 : ll.leftMin--;
288 1445840 : else if (dfR2 > m_dfMaxDistance2)
289 : {
290 26 : ll.left = nXOffset + m_nX + 1;
291 26 : break;
292 : }
293 :
294 1445820 : CheckNoData(*pdfHeight);
295 1440780 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
296 : }
297 :
298 : // Go right.
299 28008 : pdfHeight = vThisLineVal.data() + nXStart + 1;
300 28911 : for (int nXOffset = nXStart - m_nX + 1;
301 1500350 : nXOffset < oCurExtent.xSize() - m_nX; nXOffset++, pdfHeight++)
302 : {
303 1473710 : double dfX = m_adfTransform[1] * nXOffset + dfLineX;
304 1475650 : double dfY = m_adfTransform[4] * nXOffset + dfLineY;
305 1475580 : double dfR2 = dfX * dfX + dfY * dfY;
306 1475580 : if (dfR2 < m_dfMinDistance2)
307 3 : ll.rightMin++;
308 1475580 : else if (dfR2 > m_dfMaxDistance2)
309 : {
310 26 : ll.right = nXOffset + m_nX;
311 26 : break;
312 : }
313 :
314 1475560 : CheckNoData(*pdfHeight);
315 1471440 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
316 28945 : }
317 : }
318 : else
319 : {
320 : // No curvature adjustment. Just normalize for the observer height.
321 0 : double *pdfHeight = vThisLineVal.data();
322 0 : for (int i = 0; i < oCurExtent.xSize(); ++i)
323 : {
324 4 : CheckNoData(*pdfHeight);
325 0 : *pdfHeight -= m_dfZObserver;
326 0 : pdfHeight++;
327 : }
328 : }
329 28941 : return ll;
330 : }
331 :
332 : /// Process the first line (the one with the Y coordinate the same as the observer).
333 : ///
334 : /// @param vLastLineVal Vector in which to store the read line. Becomes the last line
335 : /// in further processing.
336 : /// @return True on success, false otherwise.
337 235 : bool ViewshedExecutor::processFirstLine(std::vector<double> &vLastLineVal)
338 : {
339 235 : int nLine = oOutExtent.clampY(m_nY);
340 235 : int nYOffset = nLine - m_nY;
341 :
342 470 : std::vector<double> vResult(oOutExtent.xSize());
343 470 : std::vector<double> vThisLineVal(oOutExtent.xSize());
344 :
345 235 : if (!readLine(nLine, vThisLineVal.data()))
346 0 : return false;
347 :
348 : // If the observer is outside of the raster, take the specified value as the Z height,
349 : // otherwise, take it as an offset from the raster height at that location.
350 235 : m_dfZObserver = oOpts.observer.z;
351 235 : if (oCurExtent.containsX(m_nX))
352 : {
353 229 : m_dfZObserver += vThisLineVal[m_nX];
354 229 : if (oOpts.outputMode == OutputMode::Normal)
355 212 : vResult[m_nX] = oOpts.visibleVal;
356 : }
357 235 : m_dfHeightAdjFactor = calcHeightAdjFactor();
358 :
359 : // In DEM mode the base is the pre-adjustment value. In ground mode the base is zero.
360 235 : if (oOpts.outputMode == OutputMode::DEM)
361 16 : vResult = vThisLineVal;
362 :
363 235 : LineLimits ll = adjustHeight(nYOffset, vThisLineVal);
364 235 : if (oCurExtent.containsX(m_nX) && ll.leftMin != ll.rightMin)
365 1 : vResult[m_nX] = oOpts.outOfRangeVal;
366 :
367 235 : if (!oCurExtent.containsY(m_nY))
368 4 : processFirstLineTopOrBottom(ll, vResult, vThisLineVal);
369 : else
370 : {
371 462 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
372 231 : pQueue->SubmitJob([&]()
373 231 : { processFirstLineLeft(ll, vResult, vThisLineVal); });
374 :
375 231 : pQueue->SubmitJob(
376 231 : [&]() { processFirstLineRight(ll, vResult, vThisLineVal); });
377 231 : pQueue->WaitCompletion();
378 : }
379 :
380 : // Make the current line the last line.
381 235 : vLastLineVal = std::move(vThisLineVal);
382 :
383 : // Create the output writer.
384 235 : if (!writeLine(nLine, vResult))
385 0 : return false;
386 :
387 235 : return oProgress.lineComplete();
388 : }
389 :
390 : // If the observer is above or below the raster, set all cells in the first line near the
391 : // observer as observable provided they're in range. Mark cells out of range as such.
392 : /// @param ll Line limits for processing.
393 : /// @param vResult Result line.
394 : /// @param vThisLineVal Heights of the cells in the target line
395 4 : void ViewshedExecutor::processFirstLineTopOrBottom(
396 : const LineLimits &ll, std::vector<double> &vResult,
397 : std::vector<double> &vThisLineVal)
398 : {
399 4 : double *pResult = vResult.data() + ll.left;
400 4 : double *pThis = vThisLineVal.data() + ll.left;
401 24 : for (int iPixel = ll.left; iPixel < ll.right; ++iPixel, ++pResult, pThis++)
402 : {
403 20 : if (oOpts.outputMode == OutputMode::Normal)
404 0 : *pResult = oOpts.visibleVal;
405 : else
406 20 : setOutput(*pResult, *pThis, *pThis);
407 : }
408 :
409 4 : std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
410 4 : std::fill(vResult.begin() + ll.right, vResult.begin() + oCurExtent.xStop,
411 4 : oOpts.outOfRangeVal);
412 4 : }
413 :
414 : /// Process the part of the first line to the left of the observer.
415 : ///
416 : /// @param ll Line limits for masking.
417 : /// @param vResult Vector in which to store the visibility/height results.
418 : /// @param vThisLineVal Height of each cell in the line being processed.
419 231 : void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll,
420 : std::vector<double> &vResult,
421 : std::vector<double> &vThisLineVal)
422 : {
423 231 : int iEnd = ll.left - 1;
424 231 : int iStart = m_nX - 1; // One left of the observer.
425 :
426 : // If end is to the right of start, everything is taken care of by right processing.
427 231 : if (iEnd >= iStart)
428 30 : return;
429 :
430 201 : iStart = oCurExtent.clampX(iStart);
431 :
432 201 : double *pThis = vThisLineVal.data() + iStart;
433 :
434 : // If the start cell is next to the observer, just mark it visible.
435 201 : if (iStart + 1 == m_nX || iStart + 1 == oCurExtent.xStop)
436 : {
437 201 : double dfZ = *pThis;
438 201 : if (oOpts.outputMode == OutputMode::Normal)
439 190 : vResult[iStart] = oOpts.visibleVal;
440 : else
441 : {
442 11 : maskLowPitch(dfZ, 1, 0);
443 11 : setOutput(vResult[iStart], *pThis, dfZ);
444 : }
445 201 : maskHighPitch(vResult[iStart], dfZ, 1, 0);
446 201 : iStart--;
447 201 : pThis--;
448 : }
449 :
450 : // Go from the observer to the left, calculating Z as we go.
451 10357 : for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--)
452 : {
453 10156 : int nXOffset = std::abs(iPixel - m_nX);
454 10156 : double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
455 10156 : maskLowPitch(dfZ, nXOffset, 0);
456 10156 : setOutput(vResult[iPixel], *pThis, dfZ);
457 10156 : maskHighPitch(vResult[iPixel], dfZ, nXOffset, 0);
458 : }
459 :
460 201 : maskLineLeft(vResult, ll, m_nY);
461 : }
462 :
463 : /// Mask cells based on angle intersection to the left of the observer.
464 : ///
465 : /// @param vResult Result raaster line.
466 : /// @param nLine Line number.
467 : /// @return True when all cells have been masked.
468 25958 : bool ViewshedExecutor::maskAngleLeft(std::vector<double> &vResult, int nLine)
469 : {
470 38 : auto clamp = [this](int x)
471 20 : { return (x < 0 || x >= m_nX) ? INVALID_ISECT : x; };
472 :
473 25958 : if (!oOpts.angleMasking())
474 25893 : return false;
475 :
476 13 : if (nLine != m_nY)
477 : {
478 : int startAngleX =
479 10 : clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
480 10 : int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
481 : // If neither X intersect is in the quadrant and a ray in the quadrant isn't
482 : // between start and stop, fill it all and return true. If it is in between
483 : // start and stop, we're done.
484 10 : if (invalid(startAngleX) && invalid(endAngleX))
485 : {
486 : // Choose a test angle in quadrant II or III depending on the line.
487 7 : double testAngle = nLine < m_nY ? m_testAngle[2] : m_testAngle[3];
488 7 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
489 : {
490 2 : std::fill(vResult.begin(), vResult.begin() + m_nX,
491 2 : oOpts.outOfRangeVal);
492 7 : return true;
493 : }
494 5 : return false;
495 : }
496 3 : if (nLine > m_nY)
497 0 : std::swap(startAngleX, endAngleX);
498 3 : if (invalid(startAngleX))
499 3 : startAngleX = 0;
500 3 : if (invalid(endAngleX))
501 0 : endAngleX = m_nX - 1;
502 3 : if (startAngleX <= endAngleX)
503 : {
504 3 : std::fill(vResult.begin(), vResult.begin() + startAngleX,
505 3 : oOpts.outOfRangeVal);
506 3 : std::fill(vResult.begin() + endAngleX + 1, vResult.begin() + m_nX,
507 3 : oOpts.outOfRangeVal);
508 : }
509 : else
510 : {
511 0 : std::fill(vResult.begin() + endAngleX + 1,
512 0 : vResult.begin() + startAngleX, oOpts.outOfRangeVal);
513 : }
514 : }
515 : // nLine == m_nY
516 3 : else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI))
517 : {
518 0 : std::fill(vResult.begin(), vResult.begin() + m_nX, oOpts.outOfRangeVal);
519 0 : return true;
520 : }
521 4 : return false;
522 : }
523 :
524 : /// Mask cells based on angle intersection to the right 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 28907 : bool ViewshedExecutor::maskAngleRight(std::vector<double> &vResult, int nLine)
530 : {
531 28907 : int lineLength = static_cast<int>(vResult.size());
532 :
533 54 : auto clamp = [this, lineLength](int x)
534 36 : { return (x <= m_nX || x >= lineLength) ? INVALID_ISECT : x; };
535 :
536 28877 : if (oOpts.startAngle == oOpts.endAngle)
537 28845 : return false;
538 :
539 32 : if (nLine != m_nY)
540 : {
541 : int startAngleX =
542 18 : clamp(hIntersect(oOpts.startAngle, m_nX, m_nY, nLine));
543 18 : int endAngleX = clamp(hIntersect(oOpts.endAngle, m_nX, m_nY, nLine));
544 :
545 : // If neither X intersect is in the quadrant and a ray in the quadrant isn't
546 : // between start and stop, fill it all and return true. If it is in between
547 : // start and stop, we're done.
548 18 : if (invalid(startAngleX) && invalid(endAngleX))
549 : {
550 : // Choose a test angle in quadrant I or IV depending on the line.
551 10 : double testAngle = nLine < m_nY ? m_testAngle[1] : m_testAngle[4];
552 10 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, testAngle))
553 : {
554 0 : std::fill(vResult.begin() + m_nX + 1, vResult.end(),
555 0 : oOpts.outOfRangeVal);
556 10 : return true;
557 : }
558 10 : return false;
559 : }
560 :
561 8 : if (nLine > m_nY)
562 0 : std::swap(startAngleX, endAngleX);
563 8 : if (invalid(endAngleX))
564 0 : endAngleX = lineLength - 1;
565 8 : if (invalid(startAngleX))
566 8 : startAngleX = m_nX + 1;
567 8 : if (startAngleX <= endAngleX)
568 : {
569 8 : std::fill(vResult.begin() + m_nX + 1, vResult.begin() + startAngleX,
570 8 : oOpts.outOfRangeVal);
571 8 : std::fill(vResult.begin() + endAngleX + 1, vResult.end(),
572 8 : oOpts.outOfRangeVal);
573 : }
574 : else
575 : {
576 0 : std::fill(vResult.begin() + endAngleX + 1,
577 0 : vResult.begin() + startAngleX, oOpts.outOfRangeVal);
578 : }
579 : }
580 : // nLine == m_nY
581 14 : else if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 0))
582 : {
583 1 : std::fill(vResult.begin() + m_nX + 1, vResult.end(),
584 1 : oOpts.outOfRangeVal);
585 1 : return true;
586 : }
587 9 : return false;
588 : }
589 :
590 : /// Perform angle and min/max masking to the left of the observer.
591 : ///
592 : /// @param vResult Raster line to mask.
593 : /// @param ll Min/max line limits.
594 : /// @param nLine Line number.
595 25947 : void ViewshedExecutor::maskLineLeft(std::vector<double> &vResult,
596 : const LineLimits &ll, int nLine)
597 : {
598 : // If we've already masked with angles everything, just return.
599 25947 : if (maskAngleLeft(vResult, nLine))
600 2 : return;
601 :
602 : // Mask cells from the left edge to the left limit.
603 25923 : std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
604 : // Mask cells from the left min to the observer.
605 25874 : if (ll.leftMin < m_nX)
606 3 : std::fill(vResult.begin() + ll.leftMin, vResult.begin() + m_nX,
607 3 : oOpts.outOfRangeVal);
608 : }
609 :
610 : /// Perform angle and min/max masking to the right of the observer.
611 : ///
612 : /// @param vResult Raster line to mask.
613 : /// @param ll Min/max line limits.
614 : /// @param nLine Line number.
615 28898 : void ViewshedExecutor::maskLineRight(std::vector<double> &vResult,
616 : const LineLimits &ll, int nLine)
617 : {
618 : // If we've already masked with angles everything, just return.
619 28898 : if (maskAngleRight(vResult, nLine))
620 1 : return;
621 :
622 : // Mask cells from the observer to right min.
623 28879 : std::fill(vResult.begin() + m_nX + 1, vResult.begin() + ll.rightMin,
624 28892 : oOpts.outOfRangeVal);
625 : // Mask cells from the right limit to the right edge.
626 28888 : if (ll.right + 1 < static_cast<int>(vResult.size()))
627 8 : std::fill(vResult.begin() + ll.right + 1, vResult.end(),
628 8 : oOpts.outOfRangeVal);
629 : }
630 :
631 : /// Process the part of the first line to the right of the observer.
632 : ///
633 : /// @param ll Line limits
634 : /// @param vResult Vector in which to store the visibility/height results.
635 : /// @param vThisLineVal Height of each cell in the line being processed.
636 231 : void ViewshedExecutor::processFirstLineRight(const LineLimits &ll,
637 : std::vector<double> &vResult,
638 : std::vector<double> &vThisLineVal)
639 : {
640 231 : int iStart = m_nX + 1;
641 231 : int iEnd = ll.right;
642 :
643 : // If start is to the right of end, everything is taken care of by left processing.
644 231 : if (iStart >= iEnd)
645 2 : return;
646 :
647 229 : iStart = oCurExtent.clampX(iStart);
648 :
649 229 : double *pThis = vThisLineVal.data() + iStart;
650 :
651 : // If the start cell is next to the observer, just mark it visible.
652 229 : if (iStart - 1 == m_nX || iStart == oCurExtent.xStart)
653 : {
654 229 : double dfZ = *pThis;
655 229 : if (oOpts.outputMode == OutputMode::Normal)
656 212 : vResult[iStart] = oOpts.visibleVal;
657 : else
658 : {
659 17 : maskLowPitch(dfZ, 1, 0);
660 17 : setOutput(vResult[iStart], *pThis, dfZ);
661 : }
662 229 : maskHighPitch(vResult[iStart], dfZ, 1, 0);
663 229 : iStart++;
664 229 : pThis++;
665 : }
666 :
667 : // Go from the observer to the right, calculating Z as we go.
668 10812 : for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++)
669 : {
670 10583 : int nXOffset = std::abs(iPixel - m_nX);
671 10583 : double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
672 10583 : maskLowPitch(dfZ, nXOffset, 0);
673 10583 : setOutput(vResult[iPixel], *pThis, dfZ);
674 10583 : maskHighPitch(vResult[iPixel], dfZ, nXOffset, 0);
675 : }
676 :
677 229 : maskLineRight(vResult, ll, m_nY);
678 : }
679 :
680 : /// Process a line to the left of the observer.
681 : ///
682 : /// @param nYOffset Offset of the line being processed from the observer
683 : /// @param ll Line limits
684 : /// @param vResult Vector in which to store the visibility/height results.
685 : /// @param vThisLineVal Height of each cell in the line being processed.
686 : /// @param vLastLineVal Observable height of each cell in the previous line processed.
687 28619 : void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
688 : std::vector<double> &vResult,
689 : std::vector<double> &vThisLineVal,
690 : std::vector<double> &vLastLineVal)
691 : {
692 28619 : int iStart = m_nX - 1;
693 28619 : int iEnd = ll.left - 1;
694 28619 : int nLine = m_nY + nYOffset;
695 : // If start to the left of end, everything is taken care of by processing right.
696 28619 : if (iStart <= iEnd)
697 2926 : return;
698 25693 : iStart = oCurExtent.clampX(iStart);
699 :
700 25638 : nYOffset = std::abs(nYOffset);
701 25638 : double *pThis = vThisLineVal.data() + iStart;
702 25570 : double *pLast = vLastLineVal.data() + iStart;
703 :
704 : // If the observer is to the right of the raster, mark the first cell to the left as
705 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
706 : // with the out of range assignment at the end.
707 :
708 25629 : if (iStart == oCurExtent.xStop - 1)
709 : {
710 6 : if (oOpts.outputMode == OutputMode::Normal)
711 0 : vResult[iStart] = oOpts.visibleVal;
712 : else
713 : {
714 6 : maskLowPitch(*pThis, m_nX - iStart, nYOffset);
715 6 : setOutput(vResult[iStart], *pThis, *pThis);
716 6 : maskHighPitch(vResult[iStart], *pThis, m_nX - iStart, nYOffset);
717 : }
718 6 : iStart--;
719 6 : pThis--;
720 6 : pLast--;
721 : }
722 :
723 : // Go from the observer to the left, calculating Z as we go.
724 1416400 : for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--, pLast--)
725 : {
726 1390640 : int nXOffset = std::abs(iPixel - m_nX);
727 : double dfZ;
728 1390640 : if (nXOffset == nYOffset)
729 : {
730 15644 : if (nXOffset == 1)
731 375 : dfZ = *pThis;
732 : else
733 15269 : dfZ = CalcHeightLine(nXOffset, *(pLast + 1));
734 : }
735 : else
736 1391850 : dfZ =
737 1375000 : oZcalc(nXOffset, nYOffset, *(pThis + 1), *pLast, *(pLast + 1));
738 1407500 : maskLowPitch(dfZ, nXOffset, nYOffset);
739 1399520 : setOutput(vResult[iPixel], *pThis, dfZ);
740 1393600 : maskHighPitch(vResult[iPixel], dfZ, nXOffset, nYOffset);
741 : }
742 :
743 25755 : maskLineLeft(vResult, ll, nLine);
744 : }
745 :
746 : /// Process a line to the right of the observer.
747 : ///
748 : /// @param nYOffset Offset of the line being processed from the observer
749 : /// @param ll Line limits
750 : /// @param vResult Vector in which to store the visibility/height results.
751 : /// @param vThisLineVal Height of each cell in the line being processed.
752 : /// @param vLastLineVal Observable height of each cell in the previous line processed.
753 28644 : void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
754 : std::vector<double> &vResult,
755 : std::vector<double> &vThisLineVal,
756 : std::vector<double> &vLastLineVal)
757 : {
758 28644 : int iStart = m_nX + 1;
759 28644 : int iEnd = ll.right;
760 28644 : int nLine = m_nY + nYOffset;
761 :
762 : // If start is to the right of end, everything is taken care of by processing left.
763 28644 : if (iStart >= iEnd)
764 12 : return;
765 28632 : iStart = oCurExtent.clampX(iStart);
766 :
767 28626 : nYOffset = std::abs(nYOffset);
768 28626 : double *pThis = vThisLineVal.data() + iStart;
769 28569 : double *pLast = vLastLineVal.data() + iStart;
770 :
771 : // If the observer is to the left of the raster, mark the first cell to the right as
772 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
773 : // with the out of range assignment at the end.
774 28634 : if (iStart == 0)
775 : {
776 6 : if (oOpts.outputMode == OutputMode::Normal)
777 0 : vResult[iStart] = oOpts.visibleVal;
778 : else
779 : {
780 6 : maskLowPitch(*pThis, m_nX, nYOffset);
781 6 : setOutput(vResult[0], *pThis, *pThis);
782 6 : maskHighPitch(vResult[0], *pThis, m_nX, nYOffset);
783 : }
784 6 : iStart++;
785 6 : pThis++;
786 6 : pLast++;
787 : }
788 :
789 : // Go from the observer to the right, calculating Z as we go.
790 1475280 : for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++, pLast++)
791 : {
792 1446600 : int nXOffset = std::abs(iPixel - m_nX);
793 : double dfZ;
794 1446600 : if (nXOffset == nYOffset)
795 : {
796 16127 : if (nXOffset == 1)
797 416 : dfZ = *pThis;
798 : else
799 15711 : dfZ = CalcHeightLine(nXOffset, *(pLast - 1));
800 : }
801 : else
802 1444900 : dfZ =
803 1430480 : oZcalc(nXOffset, nYOffset, *(pThis - 1), *pLast, *(pLast - 1));
804 1461020 : maskLowPitch(dfZ, nXOffset, nYOffset);
805 1456480 : setOutput(vResult[iPixel], *pThis, dfZ);
806 1457750 : maskHighPitch(vResult[iPixel], dfZ, nXOffset, nYOffset);
807 : }
808 :
809 28680 : maskLineRight(vResult, ll, nLine);
810 : }
811 :
812 : /// Apply angular mask to the initial X position. Assumes m_nX is in the raster.
813 : /// @param vResult Raster line on which to apply mask.
814 : /// @param nLine Line number.
815 28667 : void ViewshedExecutor::maskInitial(std::vector<double> &vResult, int nLine)
816 : {
817 28667 : if (!oOpts.angleMasking())
818 28561 : return;
819 :
820 19 : if (nLine < m_nY)
821 : {
822 13 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI / 2))
823 0 : vResult[m_nX] = oOpts.outOfRangeVal;
824 : }
825 6 : else if (nLine > m_nY)
826 : {
827 5 : if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 3 * M_PI / 2))
828 0 : vResult[m_nX] = oOpts.outOfRangeVal;
829 : }
830 : }
831 :
832 : /// Process a line above or below the observer.
833 : ///
834 : /// @param nLine Line number being processed.
835 : /// @param vLastLineVal Vector in which to store the read line. Becomes the last line
836 : /// in further processing.
837 : /// @return True on success, false otherwise.
838 28724 : bool ViewshedExecutor::processLine(int nLine, std::vector<double> &vLastLineVal)
839 : {
840 28724 : int nYOffset = nLine - m_nY;
841 57442 : std::vector<double> vResult(oOutExtent.xSize());
842 57434 : std::vector<double> vThisLineVal(oOutExtent.xSize());
843 :
844 28723 : if (!readLine(nLine, vThisLineVal.data()))
845 0 : return false;
846 :
847 : // In DEM mode the base is the input DEM value.
848 : // In ground mode the base is zero.
849 28691 : if (oOpts.outputMode == OutputMode::DEM)
850 163 : vResult = vThisLineVal;
851 :
852 : // Adjust height of the read line.
853 28691 : LineLimits ll = adjustHeight(nYOffset, vThisLineVal);
854 :
855 : // Handle the initial position on the line.
856 28706 : if (oCurExtent.containsX(m_nX))
857 : {
858 28672 : if (ll.left < ll.right && ll.leftMin == ll.rightMin)
859 : {
860 : double dfZ;
861 28663 : if (std::abs(nYOffset) == 1)
862 414 : dfZ = vThisLineVal[m_nX];
863 : else
864 28249 : dfZ = CalcHeightLine(nYOffset, vLastLineVal[m_nX]);
865 28628 : maskLowPitch(dfZ, 0, nYOffset);
866 28613 : setOutput(vResult[m_nX], vThisLineVal[m_nX], dfZ);
867 28609 : maskHighPitch(vResult[m_nX], dfZ, 0, nYOffset);
868 : }
869 : else
870 9 : vResult[m_nX] = oOpts.outOfRangeVal;
871 :
872 28603 : maskInitial(vResult, nLine);
873 : }
874 :
875 : // process left half then right half of line
876 57317 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
877 28571 : pQueue->SubmitJob(
878 28615 : [&]() {
879 28615 : processLineLeft(nYOffset, ll, vResult, vThisLineVal, vLastLineVal);
880 28556 : });
881 28634 : pQueue->SubmitJob(
882 28642 : [&]() {
883 28642 : processLineRight(nYOffset, ll, vResult, vThisLineVal, vLastLineVal);
884 28676 : });
885 28641 : pQueue->WaitCompletion();
886 : // Make the current line the last line.
887 28650 : vLastLineVal = std::move(vThisLineVal);
888 :
889 28606 : if (!writeLine(nLine, vResult))
890 0 : return false;
891 :
892 28661 : return oProgress.lineComplete();
893 : }
894 :
895 : // Calculate the ray angle from the origin to middle of the top or bottom
896 : // of each quadrant.
897 2 : void ViewshedExecutor::calcTestAngles()
898 : {
899 : // Quadrant 1.
900 : {
901 2 : int ysize = m_nY + 1;
902 2 : int xsize = oCurExtent.xStop - m_nX;
903 2 : m_testAngle[1] = atan2(ysize, xsize / 2.0);
904 : }
905 :
906 : // Quadrant 2.
907 : {
908 2 : int ysize = m_nY + 1;
909 2 : int xsize = m_nX + 1;
910 2 : m_testAngle[2] = atan2(ysize, -xsize / 2.0);
911 : }
912 :
913 : // Quadrant 3.
914 : {
915 2 : int ysize = oCurExtent.yStop - m_nY;
916 2 : int xsize = m_nX + 1;
917 2 : m_testAngle[3] = atan2(-ysize, -xsize / 2.0);
918 : }
919 :
920 : // Quadrant 4.
921 : {
922 2 : int ysize = oCurExtent.yStop - m_nY;
923 2 : int xsize = oCurExtent.xStop - m_nX;
924 2 : m_testAngle[4] = atan2(-ysize, xsize / 2.0);
925 : }
926 :
927 : // Adjust range to [0, 2 * M_PI)
928 10 : for (int i = 1; i <= 4; ++i)
929 8 : if (m_testAngle[i] < 0)
930 4 : m_testAngle[i] += (2 * M_PI);
931 2 : }
932 :
933 : /// Run the viewshed computation
934 : /// @return Success as true or false.
935 235 : bool ViewshedExecutor::run()
936 : {
937 : // If we're doing angular masking, calculate the test angles used later.
938 235 : if (oOpts.angleMasking())
939 2 : calcTestAngles();
940 :
941 470 : std::vector<double> vFirstLineVal(oCurExtent.xSize());
942 235 : if (!processFirstLine(vFirstLineVal))
943 0 : return false;
944 :
945 235 : if (oOpts.cellMode == CellMode::Edge)
946 235 : oZcalc = doEdge;
947 0 : else if (oOpts.cellMode == CellMode::Diagonal)
948 0 : oZcalc = doDiagonal;
949 0 : else if (oOpts.cellMode == CellMode::Min)
950 0 : oZcalc = doMin;
951 0 : else if (oOpts.cellMode == CellMode::Max)
952 0 : oZcalc = doMax;
953 :
954 : // scan upwards
955 235 : int yStart = oCurExtent.clampY(m_nY);
956 235 : std::atomic<bool> err(false);
957 235 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
958 235 : pQueue->SubmitJob(
959 235 : [&]()
960 : {
961 470 : std::vector<double> vLastLineVal = vFirstLineVal;
962 :
963 13515 : for (int nLine = yStart - 1; nLine >= oCurExtent.yStart && !err;
964 : nLine--)
965 13281 : if (!processLine(nLine, vLastLineVal))
966 0 : err = true;
967 235 : });
968 :
969 : // scan downwards
970 235 : pQueue->SubmitJob(
971 235 : [&]()
972 : {
973 470 : std::vector<double> vLastLineVal = vFirstLineVal;
974 :
975 15677 : for (int nLine = yStart + 1; nLine < oCurExtent.yStop && !err;
976 : nLine++)
977 15442 : if (!processLine(nLine, vLastLineVal))
978 0 : err = true;
979 235 : });
980 235 : return true;
981 : }
982 :
983 : /// Mask cells lower than the low pitch angle of intersection by setting the value
984 : /// to the intersection value.
985 : ///
986 : /// @param dfZ Initial/modified cell height.
987 : /// @param nXOffset Cell X offset from observer.
988 : /// @param nYOffset Cell Y offset from observer.
989 2893160 : void ViewshedExecutor::maskLowPitch(double &dfZ, int nXOffset, int nYOffset)
990 : {
991 2893160 : if (std::isnan(m_lowTanPitch))
992 2883890 : return;
993 :
994 0 : double dfX = m_adfTransform[1] * nXOffset + m_adfTransform[2] * nYOffset;
995 24 : double dfY = m_adfTransform[4] * nXOffset + m_adfTransform[5] * nYOffset;
996 24 : double dfR2 = dfX * dfX + dfY * dfY;
997 24 : double dfDist = std::sqrt(dfR2);
998 24 : double dfZmask = dfDist * m_lowTanPitch;
999 24 : if (dfZmask > dfZ)
1000 16 : dfZ = dfZmask;
1001 : }
1002 :
1003 : /// Mask cells higher than the high pitch angle of intersection by setting the value
1004 : /// to out-of-range.
1005 : ///
1006 : /// @param dfResult Result value.
1007 : /// @param dfZ Cell height.
1008 : /// @param nXOffset Cell X offset from observer.
1009 : /// @param nYOffset Cell Y offset from observer.
1010 2860390 : void ViewshedExecutor::maskHighPitch(double &dfResult, double dfZ, int nXOffset,
1011 : int nYOffset)
1012 : {
1013 2860390 : if (std::isnan(m_highTanPitch))
1014 2852720 : return;
1015 :
1016 0 : double dfX = m_adfTransform[1] * nXOffset + m_adfTransform[2] * nYOffset;
1017 224 : double dfY = m_adfTransform[4] * nXOffset + m_adfTransform[5] * nYOffset;
1018 223 : double dfR2 = dfX * dfX + dfY * dfY;
1019 223 : double dfDist = std::sqrt(dfR2);
1020 223 : double dfZmask = dfDist * m_highTanPitch;
1021 223 : if (dfZmask < dfZ)
1022 3 : dfResult = oOpts.outOfRangeVal;
1023 : }
1024 :
1025 : } // namespace viewshed
1026 : } // namespace gdal
|
