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 <limits>
18 :
19 : #include "viewshed_executor.h"
20 : #include "progress.h"
21 :
22 : namespace gdal
23 : {
24 : namespace viewshed
25 : {
26 :
27 : namespace
28 : {
29 :
30 : /// Calculate the height at nDistance units along a line through the origin given the height
31 : /// at nDistance - 1 units along the line.
32 : /// \param nDistance Distance along the line for the target point.
33 : /// \param Za Height at the line one unit previous to the target point.
34 231037 : double CalcHeightLine(int nDistance, double Za)
35 : {
36 231037 : nDistance = std::abs(nDistance);
37 231037 : assert(nDistance != 1);
38 231037 : return Za * nDistance / (nDistance - 1);
39 : }
40 :
41 : // Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
42 : // and passing through the line connecting (i - 1, j, Za) and (i, j - 1, Zb).
43 : // In other words, the origin and the two points form a plane and we're calculating Zc
44 : // of the point (i, j, Zc), also on the plane.
45 0 : double CalcHeightDiagonal(int i, int j, double Za, double Zb)
46 : {
47 0 : return (Za * i + Zb * j) / (i + j - 1);
48 : }
49 :
50 : // Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
51 : // and through the line connecting (i -1, j - 1, Za) and (i - 1, j, Zb). In other words,
52 : // the origin and the other two points form a plane and we're calculating Zc of the
53 : // point (i, j, Zc), also on the plane.
54 8292380 : double CalcHeightEdge(int i, int j, double Za, double Zb)
55 : {
56 8292380 : assert(i != j);
57 8292380 : return (Za * i + Zb * (j - i)) / (j - 1);
58 : }
59 :
60 0 : double doDiagonal(int nXOffset, [[maybe_unused]] int nYOffset,
61 : double dfThisPrev, double dfLast,
62 : [[maybe_unused]] double dfLastPrev)
63 : {
64 0 : return CalcHeightDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast);
65 : }
66 :
67 8274870 : double doEdge(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
68 : double dfLastPrev)
69 : {
70 8274870 : if (nXOffset >= nYOffset)
71 3332280 : return CalcHeightEdge(nYOffset, nXOffset, dfLastPrev, dfThisPrev);
72 : else
73 4942590 : return CalcHeightEdge(nXOffset, nYOffset, dfLastPrev, dfLast);
74 : }
75 :
76 0 : double doMin(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
77 : double dfLastPrev)
78 : {
79 0 : double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
80 : double dfDiagonal =
81 0 : doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
82 0 : return std::min(dfEdge, dfDiagonal);
83 : }
84 :
85 0 : double doMax(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
86 : double dfLastPrev)
87 : {
88 0 : double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
89 : double dfDiagonal =
90 0 : doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
91 0 : return std::max(dfEdge, dfDiagonal);
92 : }
93 :
94 : } // unnamed namespace
95 :
96 : /// Constructor -- the viewshed algorithm executor
97 : /// @param srcBand Source raster band
98 : /// @param dstBand Destination raster band
99 : /// @param nX X position of observer
100 : /// @param nY Y position of observer
101 : /// @param outExtent Extent of output raster (relative to input)
102 : /// @param curExtent Extent of active raster.
103 : /// @param opts Configuration options.
104 : /// @param progress Reference to the progress tracker.
105 623 : ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
106 : GDALRasterBand &dstBand, int nX, int nY,
107 : const Window &outExtent,
108 : const Window &curExtent, const Options &opts,
109 623 : Progress &progress)
110 : : m_pool(4), m_srcBand(srcBand), m_dstBand(dstBand), oOutExtent(outExtent),
111 623 : oCurExtent(curExtent), m_nX(nX - oOutExtent.xStart), m_nY(nY),
112 : oOpts(opts), oProgress(progress),
113 623 : m_dfMaxDistance2(opts.maxDistance * opts.maxDistance)
114 : {
115 623 : if (m_dfMaxDistance2 == 0)
116 621 : m_dfMaxDistance2 = std::numeric_limits<double>::max();
117 623 : m_srcBand.GetDataset()->GetGeoTransform(m_adfTransform.data());
118 623 : }
119 :
120 : // calculate the height adjustment factor.
121 623 : double ViewshedExecutor::calcHeightAdjFactor()
122 : {
123 1246 : std::lock_guard g(oMutex);
124 :
125 : const OGRSpatialReference *poDstSRS =
126 623 : m_dstBand.GetDataset()->GetSpatialRef();
127 :
128 623 : if (poDstSRS)
129 : {
130 : OGRErr eSRSerr;
131 :
132 : // If we can't get a SemiMajor axis from the SRS, it will be SRS_WGS84_SEMIMAJOR
133 11 : double dfSemiMajor = poDstSRS->GetSemiMajor(&eSRSerr);
134 :
135 : /* If we fetched the axis from the SRS, use it */
136 11 : if (eSRSerr != OGRERR_FAILURE)
137 11 : return oOpts.curveCoeff / (dfSemiMajor * 2.0);
138 :
139 0 : CPLDebug("GDALViewshedGenerate",
140 : "Unable to fetch SemiMajor axis from spatial reference");
141 : }
142 612 : return 0;
143 : }
144 :
145 : /// Set the output Z value depending on the observable height and computation mode.
146 : ///
147 : /// dfResult Reference to the result cell
148 : /// dfCellVal Reference to the current cell height. Replace with observable height.
149 : /// dfZ Minimum observable height at cell.
150 8465900 : void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
151 : double dfZ)
152 : {
153 8465900 : if (oOpts.outputMode != OutputMode::Normal)
154 : {
155 43346 : dfResult += (dfZ - dfCellVal);
156 43346 : dfResult = std::max(0.0, dfResult);
157 : }
158 : else
159 8422550 : dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
160 : : oOpts.visibleVal;
161 8465610 : dfCellVal = std::max(dfCellVal, dfZ);
162 8508830 : }
163 :
164 : /// Read a line of raster data.
165 : ///
166 : /// @param nLine Line number to read.
167 : /// @param data Pointer to location in which to store data.
168 : /// @return Success or failure.
169 83721 : bool ViewshedExecutor::readLine(int nLine, double *data)
170 : {
171 167371 : std::lock_guard g(iMutex);
172 :
173 83728 : if (GDALRasterIO(&m_srcBand, GF_Read, oOutExtent.xStart, nLine,
174 : oOutExtent.xSize(), 1, data, oOutExtent.xSize(), 1,
175 83650 : GDT_Float64, 0, 0))
176 : {
177 0 : CPLError(CE_Failure, CPLE_AppDefined,
178 : "RasterIO error when reading DEM at position (%d,%d), "
179 : "size (%d,%d)",
180 0 : oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
181 0 : return false;
182 : }
183 83650 : return true;
184 : }
185 :
186 : /// Write an output line of either visibility or height data.
187 : ///
188 : /// @param nLine Line number being written.
189 : /// @param vResult Result line to write.
190 : /// @return True on success, false otherwise.
191 83517 : bool ViewshedExecutor::writeLine(int nLine, std::vector<double> &vResult)
192 : {
193 : // GDALRasterIO isn't thread-safe.
194 167171 : std::lock_guard g(oMutex);
195 :
196 166946 : if (GDALRasterIO(&m_dstBand, GF_Write, 0, nLine - oOutExtent.yStart,
197 83473 : oOutExtent.xSize(), 1, vResult.data(), oOutExtent.xSize(),
198 83654 : 1, GDT_Float64, 0, 0))
199 : {
200 0 : CPLError(CE_Failure, CPLE_AppDefined,
201 : "RasterIO error when writing target raster at position "
202 : "(%d,%d), size (%d,%d)",
203 0 : 0, nLine - oOutExtent.yStart, oOutExtent.xSize(), 1);
204 0 : return false;
205 : }
206 83654 : return true;
207 : }
208 :
209 : /// Adjust the height of the line of data by the observer height and the curvature of the
210 : /// earth.
211 : ///
212 : /// @param nYOffset Y offset of the line being adjusted.
213 : /// @param vThisLineVal Line height data.
214 : /// @return [left, right) Leftmost and one past the rightmost cell in the line within
215 : /// the max distance. Indices are limited to the raster extent (right may be just
216 : /// outside the raster).
217 : std::pair<int, int>
218 83676 : ViewshedExecutor::adjustHeight(int nYOffset, std::vector<double> &vThisLineVal)
219 : {
220 83676 : int nLeft = 0;
221 83676 : int nRight = oCurExtent.xSize();
222 :
223 : // Find the starting point in the raster (m_nX may be outside)
224 83586 : int nXStart = oCurExtent.clampX(m_nX);
225 :
226 : // If there is a height adjustment factor other than zero or a max distance,
227 : // calculate the adjusted height of the cell, stopping if we've exceeded the max
228 : // distance.
229 83631 : if (static_cast<bool>(m_dfHeightAdjFactor) || m_dfMaxDistance2 > 0)
230 : {
231 : // Hoist invariants from the loops.
232 83631 : const double dfLineX = m_adfTransform[2] * nYOffset;
233 83594 : const double dfLineY = m_adfTransform[5] * nYOffset;
234 :
235 : // Go left
236 83613 : double *pdfHeight = vThisLineVal.data() + nXStart;
237 4271430 : for (int nXOffset = nXStart - m_nX; nXOffset >= -m_nX;
238 : nXOffset--, pdfHeight--)
239 : {
240 4193890 : double dfX = m_adfTransform[1] * nXOffset + dfLineX;
241 4182150 : double dfY = m_adfTransform[4] * nXOffset + dfLineY;
242 4187970 : double dfR2 = dfX * dfX + dfY * dfY;
243 4187970 : if (dfR2 > m_dfMaxDistance2)
244 : {
245 102 : nLeft = nXOffset + m_nX + 1;
246 102 : break;
247 : }
248 4187860 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
249 : }
250 :
251 : // Go right.
252 77637 : pdfHeight = vThisLineVal.data() + nXStart + 1;
253 83679 : for (int nXOffset = nXStart - m_nX + 1;
254 4377040 : nXOffset < oCurExtent.xSize() - m_nX; nXOffset++, pdfHeight++)
255 : {
256 4290210 : double dfX = m_adfTransform[1] * nXOffset + dfLineX;
257 4284530 : double dfY = m_adfTransform[4] * nXOffset + dfLineY;
258 4293460 : double dfR2 = dfX * dfX + dfY * dfY;
259 4293460 : if (dfR2 > m_dfMaxDistance2)
260 : {
261 102 : nRight = nXOffset + m_nX;
262 102 : break;
263 : }
264 4293360 : *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
265 83702 : }
266 : }
267 : else
268 : {
269 : // No curvature adjustment. Just normalize for the observer height.
270 0 : double *pdfHeight = vThisLineVal.data();
271 0 : for (int i = 0; i < oCurExtent.xSize(); ++i)
272 : {
273 0 : *pdfHeight -= m_dfZObserver;
274 0 : pdfHeight++;
275 : }
276 : }
277 83702 : return {nLeft, nRight};
278 : }
279 :
280 : /// Process the first line (the one with the Y coordinate the same as the observer).
281 : ///
282 : /// @param vLastLineVal Vector in which to store the read line. Becomes the last line
283 : /// in further processing.
284 : /// @return True on success, false otherwise.
285 623 : bool ViewshedExecutor::processFirstLine(std::vector<double> &vLastLineVal)
286 : {
287 623 : int nLine = oOutExtent.clampY(m_nY);
288 623 : int nYOffset = nLine - m_nY;
289 :
290 1246 : std::vector<double> vResult(oOutExtent.xSize());
291 1246 : std::vector<double> vThisLineVal(oOutExtent.xSize());
292 :
293 623 : if (!readLine(nLine, vThisLineVal.data()))
294 0 : return false;
295 :
296 : // If the observer is outside of the raster, take the specified value as the Z height,
297 : // otherwise, take it as an offset from the raster height at that location.
298 623 : m_dfZObserver = oOpts.observer.z;
299 623 : if (oCurExtent.containsX(m_nX))
300 : {
301 617 : m_dfZObserver += vThisLineVal[m_nX];
302 617 : if (oOpts.outputMode == OutputMode::Normal)
303 600 : vResult[m_nX] = oOpts.visibleVal;
304 : }
305 623 : m_dfHeightAdjFactor = calcHeightAdjFactor();
306 :
307 : // In DEM mode the base is the pre-adjustment value. In ground mode the base is zero.
308 622 : if (oOpts.outputMode == OutputMode::DEM)
309 15 : vResult = vThisLineVal;
310 :
311 : // iLeft and iRight are the processing limits for the line.
312 622 : const auto [iLeft, iRight] = adjustHeight(nYOffset, vThisLineVal);
313 :
314 622 : if (!oCurExtent.containsY(m_nY))
315 4 : processFirstLineTopOrBottom(iLeft, iRight, vResult, vThisLineVal);
316 : else
317 : {
318 1237 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
319 1237 : pQueue->SubmitJob(
320 618 : [&, left = iLeft]() {
321 618 : processFirstLineLeft(m_nX - 1, left - 1, vResult, vThisLineVal);
322 618 : });
323 :
324 1237 : pQueue->SubmitJob(
325 618 : [&, right = iRight]()
326 619 : { processFirstLineRight(m_nX + 1, right, vResult, vThisLineVal); });
327 618 : pQueue->WaitCompletion();
328 : }
329 :
330 : // Make the current line the last line.
331 623 : vLastLineVal = std::move(vThisLineVal);
332 :
333 : // Create the output writer.
334 623 : if (!writeLine(nLine, vResult))
335 0 : return false;
336 :
337 623 : return oProgress.lineComplete();
338 : }
339 :
340 : // If the observer is above or below the raster, set all cells in the first line near the
341 : // observer as observable provided they're in range. Mark cells out of range as such.
342 : /// @param iLeft Leftmost observable raster position in range of the target line.
343 : /// @param iRight One past the rightmost observable raster position of the target line.
344 : /// @param vResult Result line.
345 : /// @param vThisLineVal Heights of the cells in the target line
346 4 : void ViewshedExecutor::processFirstLineTopOrBottom(
347 : int iLeft, int iRight, std::vector<double> &vResult,
348 : std::vector<double> &vThisLineVal)
349 : {
350 4 : double *pResult = vResult.data() + iLeft;
351 4 : double *pThis = vThisLineVal.data() + iLeft;
352 24 : for (int iPixel = iLeft; iPixel < iRight; ++iPixel, ++pResult, pThis++)
353 : {
354 20 : if (oOpts.outputMode == OutputMode::Normal)
355 0 : *pResult = oOpts.visibleVal;
356 : else
357 20 : setOutput(*pResult, *pThis, *pThis);
358 : }
359 4 : std::fill(vResult.begin(), vResult.begin() + iLeft, oOpts.outOfRangeVal);
360 4 : std::fill(vResult.begin() + iRight, vResult.begin() + oCurExtent.xStop,
361 4 : oOpts.outOfRangeVal);
362 4 : }
363 :
364 : /// Process the part of the first line to the left of the observer.
365 : ///
366 : /// @param iStart X coordinate of the first cell to the left of the observer to be procssed.
367 : /// @param iEnd X coordinate one past the last cell to be processed.
368 : /// @param vResult Vector in which to store the visibility/height results.
369 : /// @param vThisLineVal Height of each cell in the line being processed.
370 619 : void ViewshedExecutor::processFirstLineLeft(int iStart, int iEnd,
371 : std::vector<double> &vResult,
372 : std::vector<double> &vThisLineVal)
373 : {
374 : // If end is to the right of start, everything is taken care of by right processing.
375 619 : if (iEnd >= iStart)
376 36 : return;
377 :
378 583 : iStart = oCurExtent.clampX(iStart);
379 :
380 583 : double *pThis = vThisLineVal.data() + iStart;
381 :
382 : // If the start cell is next to the observer, just mark it visible.
383 583 : if (iStart + 1 == m_nX || iStart + 1 == oCurExtent.xStop)
384 : {
385 583 : if (oOpts.outputMode == OutputMode::Normal)
386 572 : vResult[iStart] = oOpts.visibleVal;
387 : else
388 11 : setOutput(vResult[iStart], *pThis, *pThis);
389 583 : iStart--;
390 583 : pThis--;
391 : }
392 :
393 : // Go from the observer to the left, calculating Z as we go.
394 29940 : for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--)
395 : {
396 29357 : int nXOffset = std::abs(iPixel - m_nX);
397 29357 : double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
398 29357 : setOutput(vResult[iPixel], *pThis, dfZ);
399 : }
400 : // For cells outside of the [start, end) range, set the outOfRange value.
401 583 : std::fill(vResult.begin(), vResult.begin() + iEnd + 1, oOpts.outOfRangeVal);
402 : }
403 :
404 : /// Process the part of the first line to the right of the observer.
405 : ///
406 : /// @param iStart X coordinate of the first cell to the right of the observer to be processed.
407 : /// @param iEnd X coordinate one past the last cell to be processed.
408 : /// @param vResult Vector in which to store the visibility/height results.
409 : /// @param vThisLineVal Height of each cell in the line being processed.
410 618 : void ViewshedExecutor::processFirstLineRight(int iStart, int iEnd,
411 : std::vector<double> &vResult,
412 : std::vector<double> &vThisLineVal)
413 : {
414 : // If start is to the right of end, everything is taken care of by left processing.
415 618 : if (iStart >= iEnd)
416 2 : return;
417 :
418 616 : iStart = oCurExtent.clampX(iStart);
419 :
420 616 : double *pThis = vThisLineVal.data() + iStart;
421 :
422 : // If the start cell is next to the observer, just mark it visible.
423 616 : if (iStart - 1 == m_nX || iStart == oCurExtent.xStart)
424 : {
425 616 : if (oOpts.outputMode == OutputMode::Normal)
426 599 : vResult[iStart] = oOpts.visibleVal;
427 : else
428 17 : setOutput(vResult[iStart], *pThis, *pThis);
429 616 : iStart++;
430 616 : pThis++;
431 : }
432 :
433 : // Go from the observer to the right, calculating Z as we go.
434 31141 : for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++)
435 : {
436 30524 : int nXOffset = std::abs(iPixel - m_nX);
437 30524 : double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
438 30521 : setOutput(vResult[iPixel], *pThis, dfZ);
439 : }
440 : // For cells outside of the [start, end) range, set the outOfRange value.
441 617 : std::fill(vResult.begin() + iEnd, vResult.end(), oOpts.outOfRangeVal);
442 : }
443 :
444 : /// Process a line to the left of the observer.
445 : ///
446 : /// @param nYOffset Offset of the line being processed from the observer
447 : /// @param iStart X coordinate of the first cell to the left of the observer to be processed.
448 : /// @param iEnd X coordinate one past the last cell to be processed.
449 : /// @param vResult Vector in which to store the visibility/height results.
450 : /// @param vThisLineVal Height of each cell in the line being processed.
451 : /// @param vLastLineVal Observable height of each cell in the previous line processed.
452 82936 : void ViewshedExecutor::processLineLeft(int nYOffset, int iStart, int iEnd,
453 : std::vector<double> &vResult,
454 : std::vector<double> &vThisLineVal,
455 : std::vector<double> &vLastLineVal)
456 : {
457 : // If start to the left of end, everything is taken care of by processing right.
458 82936 : if (iStart <= iEnd)
459 3891 : return;
460 79045 : iStart = oCurExtent.clampX(iStart);
461 :
462 79046 : nYOffset = std::abs(nYOffset);
463 79046 : double *pThis = vThisLineVal.data() + iStart;
464 78975 : double *pLast = vLastLineVal.data() + iStart;
465 :
466 : // If the observer is to the right of the raster, mark the first cell to the left as
467 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
468 : // with the out of range assignment at the end.
469 78969 : if (iStart == oCurExtent.xStop - 1)
470 : {
471 6 : if (oOpts.outputMode == OutputMode::Normal)
472 0 : vResult[iStart] = oOpts.visibleVal;
473 : else
474 6 : setOutput(vResult[iStart], *pThis, *pThis);
475 6 : iStart--;
476 6 : pThis--;
477 6 : pLast--;
478 : }
479 :
480 : // Go from the observer to the left, calculating Z as we go.
481 4190130 : for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--, pLast--)
482 : {
483 4110930 : int nXOffset = std::abs(iPixel - m_nX);
484 : double dfZ;
485 4110930 : if (nXOffset == nYOffset)
486 : {
487 45244 : if (nXOffset == 1)
488 1134 : dfZ = *pThis;
489 : else
490 44110 : dfZ = CalcHeightLine(nXOffset, *(pLast + 1));
491 : }
492 : else
493 : dfZ =
494 4065690 : oZcalc(nXOffset, nYOffset, *(pThis + 1), *pLast, *(pLast + 1));
495 4113890 : setOutput(vResult[iPixel], *pThis, dfZ);
496 : }
497 :
498 : // For cells outside of the [start, end) range, set the outOfRange value.
499 79199 : std::fill(vResult.begin(), vResult.begin() + iEnd + 1, oOpts.outOfRangeVal);
500 : }
501 :
502 : /// Process a line to the right of the observer.
503 : ///
504 : /// @param nYOffset Offset of the line being processed from the observer
505 : /// @param iStart X coordinate of the first cell to the right of the observer to be processed.
506 : /// @param iEnd X coordinate one past the last cell to be processed.
507 : /// @param vResult Vector in which to store the visibility/height results.
508 : /// @param vThisLineVal Height of each cell in the line being processed.
509 : /// @param vLastLineVal Observable height of each cell in the previous line processed.
510 82873 : void ViewshedExecutor::processLineRight(int nYOffset, int iStart, int iEnd,
511 : std::vector<double> &vResult,
512 : std::vector<double> &vThisLineVal,
513 : std::vector<double> &vLastLineVal)
514 : {
515 : // If start is to the right of end, everything is taken care of by processing left.
516 82873 : if (iStart >= iEnd)
517 10 : return;
518 82863 : iStart = oCurExtent.clampX(iStart);
519 :
520 82936 : nYOffset = std::abs(nYOffset);
521 82936 : double *pThis = vThisLineVal.data() + iStart;
522 82870 : double *pLast = vLastLineVal.data() + iStart;
523 :
524 : // If the observer is to the left of the raster, mark the first cell to the right as
525 : // visible. This may mark an out-of-range cell with a value, but this will be fixed
526 : // with the out of range assignment at the end.
527 82896 : if (iStart == 0)
528 : {
529 6 : if (oOpts.outputMode == OutputMode::Normal)
530 0 : vResult[iStart] = oOpts.visibleVal;
531 : else
532 6 : setOutput(vResult[0], *pThis, *pThis);
533 6 : iStart++;
534 6 : pThis++;
535 6 : pLast++;
536 : }
537 :
538 : // Go from the observer to the right, calculating Z as we go.
539 4371950 : for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++, pLast++)
540 : {
541 4288880 : int nXOffset = std::abs(iPixel - m_nX);
542 : double dfZ;
543 4288880 : if (nXOffset == nYOffset)
544 : {
545 46714 : if (nXOffset == 1)
546 1189 : dfZ = *pThis;
547 : else
548 45525 : dfZ = CalcHeightLine(nXOffset, *(pLast - 1));
549 : }
550 : else
551 : dfZ =
552 4242170 : oZcalc(nXOffset, nYOffset, *(pThis - 1), *pLast, *(pLast - 1));
553 4282150 : setOutput(vResult[iPixel], *pThis, dfZ);
554 : }
555 :
556 : // For cells outside of the [start, end) range, set the outOfRange value.
557 83069 : std::fill(vResult.begin() + iEnd, vResult.end(), oOpts.outOfRangeVal);
558 : }
559 :
560 : /// Process a line above or below the observer.
561 : ///
562 : /// @param nLine Line number being processed.
563 : /// @param vLastLineVal Vector in which to store the read line. Becomes the last line
564 : /// in further processing.
565 : /// @return True on success, false otherwise.
566 83104 : bool ViewshedExecutor::processLine(int nLine, std::vector<double> &vLastLineVal)
567 : {
568 83104 : int nYOffset = nLine - m_nY;
569 166210 : std::vector<double> vResult(oOutExtent.xSize());
570 166220 : std::vector<double> vThisLineVal(oOutExtent.xSize());
571 :
572 83114 : if (!readLine(nLine, vThisLineVal.data()))
573 0 : return false;
574 :
575 : // In DEM mode the base is the input DEM value.
576 : // In ground mode the base is zero.
577 83038 : if (oOpts.outputMode == OutputMode::DEM)
578 159 : vResult = vThisLineVal;
579 :
580 : // Adjust height of the read line.
581 83038 : const auto [iLeft, iRight] = adjustHeight(nYOffset, vThisLineVal);
582 :
583 : // Handle the initial position on the line.
584 82940 : if (oCurExtent.containsX(m_nX))
585 : {
586 82849 : if (iLeft < iRight)
587 : {
588 : double dfZ;
589 82845 : if (std::abs(nYOffset) == 1)
590 1189 : dfZ = vThisLineVal[m_nX];
591 : else
592 81656 : dfZ = CalcHeightLine(nYOffset, vLastLineVal[m_nX]);
593 82919 : setOutput(vResult[m_nX], vThisLineVal[m_nX], dfZ);
594 : }
595 : else
596 4 : vResult[m_nX] = oOpts.outOfRangeVal;
597 : }
598 :
599 : // process left half then right half of line
600 166025 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
601 165696 : pQueue->SubmitJob(
602 82872 : [&, left = iLeft]()
603 : {
604 82983 : processLineLeft(nYOffset, m_nX - 1, left - 1, vResult, vThisLineVal,
605 82983 : vLastLineVal);
606 82959 : });
607 165874 : pQueue->SubmitJob(
608 82905 : [&, right = iRight]()
609 : {
610 82922 : processLineRight(nYOffset, m_nX + 1, right, vResult, vThisLineVal,
611 82922 : vLastLineVal);
612 82999 : });
613 83015 : pQueue->WaitCompletion();
614 :
615 : // Make the current line the last line.
616 82987 : vLastLineVal = std::move(vThisLineVal);
617 :
618 82891 : if (!writeLine(nLine, vResult))
619 0 : return false;
620 :
621 83009 : return oProgress.lineComplete();
622 : }
623 :
624 : /// Run the viewshed computation
625 : /// @return Success as true or false.
626 623 : bool ViewshedExecutor::run()
627 : {
628 1246 : std::vector<double> vFirstLineVal(oCurExtent.xSize());
629 623 : if (!processFirstLine(vFirstLineVal))
630 0 : return false;
631 :
632 623 : if (oOpts.cellMode == CellMode::Edge)
633 623 : oZcalc = doEdge;
634 0 : else if (oOpts.cellMode == CellMode::Diagonal)
635 0 : oZcalc = doDiagonal;
636 0 : else if (oOpts.cellMode == CellMode::Min)
637 0 : oZcalc = doMin;
638 0 : else if (oOpts.cellMode == CellMode::Max)
639 0 : oZcalc = doMax;
640 :
641 : // scan upwards
642 623 : int yStart = oCurExtent.clampY(m_nY);
643 623 : std::atomic<bool> err(false);
644 623 : CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
645 623 : pQueue->SubmitJob(
646 623 : [&]()
647 : {
648 1245 : std::vector<double> vLastLineVal = vFirstLineVal;
649 :
650 41008 : for (int nLine = yStart - 1; nLine >= oCurExtent.yStart && !err;
651 : nLine--)
652 40385 : if (!processLine(nLine, vLastLineVal))
653 0 : err = true;
654 623 : });
655 :
656 : // scan downwards
657 623 : pQueue->SubmitJob(
658 622 : [&]()
659 : {
660 1253 : std::vector<double> vLastLineVal = vFirstLineVal;
661 :
662 43346 : for (int nLine = yStart + 1; nLine < oCurExtent.yStop && !err;
663 : nLine++)
664 42714 : if (!processLine(nLine, vLastLineVal))
665 0 : err = true;
666 623 : });
667 623 : return true;
668 : }
669 :
670 : } // namespace viewshed
671 : } // namespace gdal
|
