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 : ******************************************************************************
8 : *
9 : * Permission is hereby granted, free of charge, to any person obtaining a
10 : * copy of this software and associated documentation files (the "Software"),
11 : * to deal in the Software without restriction, including without limitation
12 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
13 : * and/or sell copies of the Software, and to permit persons to whom the
14 : * Software is furnished to do so, subject to the following conditions:
15 : *
16 : * The above copyright notice and this permission notice shall be included
17 : * in all copies or substantial portions of the Software.
18 : *
19 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
20 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
22 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
24 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
25 : * DEALINGS IN THE SOFTWARE.
26 : ****************************************************************************/
27 :
28 : #include <algorithm>
29 : #include <array>
30 : #include <limits>
31 :
32 : #include "gdal_alg.h"
33 : #include "gdal_priv.h"
34 : #include "gdal_priv_templates.hpp"
35 :
36 : #include "viewshed.h"
37 :
38 : /************************************************************************/
39 : /* GDALViewshedGenerate() */
40 : /************************************************************************/
41 :
42 : /**
43 : * Create viewshed from raster DEM.
44 : *
45 : * This algorithm will generate a viewshed raster from an input DEM raster
46 : * by using a modified algorithm of "Generating Viewsheds without Using
47 : * Sightlines" published at
48 : * https://www.asprs.org/wp-content/uploads/pers/2000journal/january/2000_jan_87-90.pdf
49 : * This appoach provides a relatively fast calculation, since the output raster
50 : * is generated in a single scan. The gdal/apps/gdal_viewshed.cpp mainline can
51 : * be used as an example of how to use this function. The output raster will be
52 : * of type Byte or Float64.
53 : *
54 : * \note The algorithm as implemented currently will only output meaningful
55 : * results if the georeferencing is in a projected coordinate reference system.
56 : *
57 : * @param hBand The band to read the DEM data from. Only the part of the raster
58 : * within the specified maxdistance around the observer point is processed.
59 : *
60 : * @param pszDriverName Driver name (GTiff if set to NULL)
61 : *
62 : * @param pszTargetRasterName The name of the target raster to be generated.
63 : * Must not be NULL
64 : *
65 : * @param papszCreationOptions creation options.
66 : *
67 : * @param dfObserverX observer X value (in SRS units)
68 : *
69 : * @param dfObserverY observer Y value (in SRS units)
70 : *
71 : * @param dfObserverHeight The height of the observer above the DEM surface.
72 : *
73 : * @param dfTargetHeight The height of the target above the DEM surface.
74 : * (default 0)
75 : *
76 : * @param dfVisibleVal pixel value for visibility (default 255)
77 : *
78 : * @param dfInvisibleVal pixel value for invisibility (default 0)
79 : *
80 : * @param dfOutOfRangeVal The value to be set for the cells that fall outside of
81 : * the range specified by dfMaxDistance.
82 : *
83 : * @param dfNoDataVal The value to be set for the cells that have no data.
84 : * If set to a negative value, nodata is not set.
85 : * Note: currently, no special processing of input cells at a
86 : * nodata value is done (which may result in erroneous results).
87 : *
88 : * @param dfCurvCoeff Coefficient to consider the effect of the curvature and
89 : * refraction. The height of the DEM is corrected according to the following
90 : * formula: [Height] -= dfCurvCoeff * [Target Distance]^2 / [Earth Diameter] For
91 : * the effect of the atmospheric refraction we can use 0.85714.
92 : *
93 : * @param eMode The mode of the viewshed calculation.
94 : * Possible values GVM_Diagonal = 1, GVM_Edge = 2 (default), GVM_Max = 3,
95 : * GVM_Min = 4.
96 : *
97 : * @param dfMaxDistance maximum distance range to compute viewshed.
98 : * It is also used to clamp the extent of the output
99 : * raster. If set to 0, then unlimited range is assumed, that is to say the
100 : * computation is performed on the extent of the whole
101 : * raster.
102 : *
103 : * @param pfnProgress A GDALProgressFunc that may be used to report progress
104 : * to the user, or to interrupt the algorithm. May be NULL if not required.
105 : *
106 : * @param pProgressArg The callback data for the pfnProgress function.
107 : *
108 : * @param heightMode Type of information contained in output raster. Possible
109 : * values GVOT_NORMAL = 1 (default), GVOT_MIN_TARGET_HEIGHT_FROM_DEM = 2,
110 : * GVOT_MIN_TARGET_HEIGHT_FROM_GROUND = 3
111 : *
112 : * GVOT_NORMAL returns a raster of type Byte containing
113 : * visible locations.
114 : *
115 : * GVOT_MIN_TARGET_HEIGHT_FROM_DEM and
116 : * GVOT_MIN_TARGET_HEIGHT_FROM_GROUND will return a raster of type Float64
117 : * containing the minimum target height for target to be visible from the DEM
118 : * surface or ground level respectively. Parameters dfTargetHeight, dfVisibleVal
119 : * and dfInvisibleVal will be ignored.
120 : *
121 : *
122 : * @param papszExtraOptions Future extra options. Must be set to NULL currently.
123 : *
124 : * @return not NULL output dataset on success (to be closed with GDALClose()) or
125 : * NULL if an error occurs.
126 : *
127 : * @since GDAL 3.1
128 : */
129 1 : GDALDatasetH GDALViewshedGenerate(
130 : GDALRasterBandH hBand, const char *pszDriverName,
131 : const char *pszTargetRasterName, CSLConstList papszCreationOptions,
132 : double dfObserverX, double dfObserverY, double dfObserverHeight,
133 : double dfTargetHeight, double dfVisibleVal, double dfInvisibleVal,
134 : double dfOutOfRangeVal, double dfNoDataVal, double dfCurvCoeff,
135 : GDALViewshedMode eMode, double dfMaxDistance, GDALProgressFunc pfnProgress,
136 : void *pProgressArg, GDALViewshedOutputType heightMode,
137 : [[maybe_unused]] CSLConstList papszExtraOptions)
138 : {
139 : using namespace gdal;
140 :
141 2 : Viewshed::Options oOpts;
142 1 : oOpts.outputFormat = pszDriverName;
143 1 : oOpts.outputFilename = pszTargetRasterName;
144 1 : oOpts.creationOpts = papszCreationOptions;
145 1 : oOpts.observer.x = dfObserverX;
146 1 : oOpts.observer.y = dfObserverY;
147 1 : oOpts.observer.z = dfObserverHeight;
148 1 : oOpts.targetHeight = dfTargetHeight;
149 1 : oOpts.curveCoeff = dfCurvCoeff;
150 1 : oOpts.maxDistance = dfMaxDistance;
151 1 : oOpts.nodataVal = dfNoDataVal;
152 :
153 1 : switch (eMode)
154 : {
155 1 : case GVM_Edge:
156 1 : oOpts.cellMode = Viewshed::CellMode::Edge;
157 1 : break;
158 0 : case GVM_Diagonal:
159 0 : oOpts.cellMode = Viewshed::CellMode::Diagonal;
160 0 : break;
161 0 : case GVM_Min:
162 0 : oOpts.cellMode = Viewshed::CellMode::Min;
163 0 : break;
164 0 : case GVM_Max:
165 0 : oOpts.cellMode = Viewshed::CellMode::Max;
166 0 : break;
167 : }
168 :
169 1 : switch (heightMode)
170 : {
171 0 : case GVOT_MIN_TARGET_HEIGHT_FROM_DEM:
172 0 : oOpts.outputMode = Viewshed::OutputMode::DEM;
173 0 : break;
174 1 : case GVOT_MIN_TARGET_HEIGHT_FROM_GROUND:
175 1 : oOpts.outputMode = Viewshed::OutputMode::Ground;
176 1 : break;
177 0 : case GVOT_NORMAL:
178 0 : oOpts.outputMode = Viewshed::OutputMode::Normal;
179 0 : break;
180 : }
181 :
182 1 : if (!GDALIsValueInRange<uint8_t>(dfVisibleVal))
183 : {
184 0 : CPLError(CE_Failure, CPLE_AppDefined,
185 : "dfVisibleVal out of range. Must be [0, 255].");
186 0 : return nullptr;
187 : }
188 1 : if (!GDALIsValueInRange<uint8_t>(dfInvisibleVal))
189 : {
190 0 : CPLError(CE_Failure, CPLE_AppDefined,
191 : "dfInvisibleVal out of range. Must be [0, 255].");
192 0 : return nullptr;
193 : }
194 1 : if (!GDALIsValueInRange<uint8_t>(dfOutOfRangeVal))
195 : {
196 0 : CPLError(CE_Failure, CPLE_AppDefined,
197 : "dfOutOfRangeVal out of range. Must be [0, 255].");
198 0 : return nullptr;
199 : }
200 1 : oOpts.visibleVal = static_cast<uint8_t>(dfVisibleVal);
201 1 : oOpts.invisibleVal = static_cast<uint8_t>(dfInvisibleVal);
202 1 : oOpts.outOfRangeVal = static_cast<uint8_t>(dfOutOfRangeVal);
203 :
204 1 : gdal::Viewshed v(oOpts);
205 :
206 : //ABELL - Make a function for progress that captures the progress argument.
207 1 : v.run(hBand, pfnProgress, pProgressArg);
208 :
209 1 : return GDALDataset::FromHandle(v.output().release());
210 : }
211 :
212 : namespace gdal
213 : {
214 :
215 : namespace
216 : {
217 :
218 88480 : void SetVisibility(int iPixel, double dfZ, double dfZTarget, double *padfZVal,
219 : std::vector<GByte> &vResult, GByte byVisibleVal,
220 : GByte byInvisibleVal)
221 : {
222 88480 : if (padfZVal[iPixel] + dfZTarget < dfZ)
223 57873 : vResult[iPixel] = byInvisibleVal;
224 : else
225 30607 : vResult[iPixel] = byVisibleVal;
226 :
227 88480 : if (padfZVal[iPixel] < dfZ)
228 57968 : padfZVal[iPixel] = dfZ;
229 88480 : }
230 :
231 88602 : bool AdjustHeightInRange(const double *adfGeoTransform, int iPixel, int iLine,
232 : double &dfHeight, double dfDistance2,
233 : double dfCurvCoeff, double dfSphereDiameter)
234 : {
235 88602 : if (dfDistance2 <= 0 && dfCurvCoeff == 0)
236 0 : return true;
237 :
238 88602 : double dfX = adfGeoTransform[1] * iPixel + adfGeoTransform[2] * iLine;
239 88602 : double dfY = adfGeoTransform[4] * iPixel + adfGeoTransform[5] * iLine;
240 88602 : double dfR2 = dfX * dfX + dfY * dfY;
241 :
242 : /* calc adjustment */
243 175184 : if (dfCurvCoeff != 0 &&
244 86582 : dfSphereDiameter != std::numeric_limits<double>::infinity())
245 86514 : dfHeight -= dfCurvCoeff * dfR2 / dfSphereDiameter;
246 :
247 88602 : if (dfDistance2 > 0 && dfR2 > dfDistance2)
248 104 : return false;
249 :
250 88498 : return true;
251 : }
252 :
253 2852 : double CalcHeightLine(int i, double Za, double Zo)
254 : {
255 2852 : if (i == 1)
256 54 : return Za;
257 : else
258 2798 : return (Za - Zo) / (i - 1) + Za;
259 : }
260 :
261 0 : double CalcHeightDiagonal(int i, int j, double Za, double Zb, double Zo)
262 : {
263 0 : return ((Za - Zo) * i + (Zb - Zo) * j) / (i + j - 1) + Zo;
264 : }
265 :
266 86936 : double CalcHeightEdge(int i, int j, double Za, double Zb, double Zo)
267 : {
268 86936 : if (i == j)
269 1308 : return CalcHeightLine(i, Za, Zo);
270 : else
271 85628 : return ((Za - Zo) * i + (Zb - Zo) * (j - i)) / (j - 1) + Zo;
272 : }
273 :
274 86936 : double CalcHeight(double dfZ, double dfZ2, Viewshed::CellMode eMode)
275 : {
276 86936 : double dfHeight = dfZ;
277 :
278 86936 : switch (eMode)
279 : {
280 86936 : case Viewshed::CellMode::Edge:
281 86936 : dfHeight = dfZ2;
282 86936 : break;
283 0 : case Viewshed::CellMode::Max:
284 0 : dfHeight = std::max(dfZ, dfZ2);
285 0 : break;
286 0 : case Viewshed::CellMode::Min:
287 0 : dfHeight = std::min(dfZ, dfZ2);
288 0 : break;
289 0 : case Viewshed::CellMode::Diagonal:
290 0 : dfHeight = dfZ;
291 0 : break;
292 : }
293 86936 : return dfHeight;
294 : }
295 :
296 : } // unnamed namespace
297 :
298 12 : bool Viewshed::run(GDALRasterBandH hBand, GDALProgressFunc pfnProgress,
299 : void *pProgressArg)
300 : {
301 12 : if (!pfnProgress)
302 1 : pfnProgress = GDALDummyProgress;
303 :
304 12 : if (!pfnProgress(0.0, "", pProgressArg))
305 : {
306 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
307 0 : return false;
308 : }
309 :
310 : /* set up geotransformation */
311 12 : std::array<double, 6> adfGeoTransform{{0.0, 1.0, 0.0, 0.0, 0.0, 1.0}};
312 12 : GDALDatasetH hSrcDS = GDALGetBandDataset(hBand);
313 12 : if (hSrcDS != nullptr)
314 12 : GDALGetGeoTransform(hSrcDS, adfGeoTransform.data());
315 :
316 : double adfInvGeoTransform[6];
317 12 : if (!GDALInvGeoTransform(adfGeoTransform.data(), adfInvGeoTransform))
318 : {
319 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
320 0 : return false;
321 : }
322 :
323 : /* calculate observer position */
324 : double dfX, dfY;
325 12 : GDALApplyGeoTransform(adfInvGeoTransform, oOpts.observer.x,
326 : oOpts.observer.y, &dfX, &dfY);
327 12 : int nX = static_cast<int>(dfX);
328 12 : int nY = static_cast<int>(dfY);
329 :
330 12 : int nXSize = GDALGetRasterBandXSize(hBand);
331 12 : int nYSize = GDALGetRasterBandYSize(hBand);
332 :
333 12 : if (nX < 0 || nX > nXSize || nY < 0 || nY > nYSize)
334 : {
335 1 : CPLError(CE_Failure, CPLE_AppDefined,
336 : "The observer location falls outside of the DEM area");
337 1 : return false;
338 : }
339 :
340 : /* calculate the area of interest */
341 11 : constexpr double EPSILON = 1e-8;
342 :
343 11 : int nXStart = 0;
344 11 : int nYStart = 0;
345 11 : int nXStop = nXSize;
346 11 : int nYStop = nYSize;
347 11 : if (oOpts.maxDistance > 0)
348 : {
349 1 : nXStart = static_cast<int>(std::floor(
350 1 : nX - adfInvGeoTransform[1] * oOpts.maxDistance + EPSILON));
351 1 : nXStop = static_cast<int>(
352 1 : std::ceil(nX + adfInvGeoTransform[1] * oOpts.maxDistance -
353 1 : EPSILON) +
354 : 1);
355 1 : nYStart =
356 1 : static_cast<int>(std::floor(
357 1 : nY - std::fabs(adfInvGeoTransform[5]) * oOpts.maxDistance +
358 1 : EPSILON)) -
359 1 : (adfInvGeoTransform[5] > 0 ? 1 : 0);
360 1 : nYStop = static_cast<int>(
361 1 : std::ceil(nY +
362 1 : std::fabs(adfInvGeoTransform[5]) * oOpts.maxDistance -
363 1 : EPSILON) +
364 1 : (adfInvGeoTransform[5] < 0 ? 1 : 0));
365 : }
366 11 : nXStart = std::max(nXStart, 0);
367 11 : nYStart = std::max(nYStart, 0);
368 11 : nXStop = std::min(nXStop, nXSize);
369 11 : nYStop = std::min(nYStop, nYSize);
370 :
371 : /* normalize horizontal index (0 - nXSize) */
372 11 : nXSize = nXStop - nXStart;
373 11 : nX -= nXStart;
374 :
375 11 : nYSize = nYStop - nYStart;
376 :
377 11 : if (nXSize == 0 || nYSize == 0)
378 : {
379 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid target raster size");
380 0 : return false;
381 : }
382 :
383 22 : std::vector<double> vFirstLineVal;
384 22 : std::vector<double> vLastLineVal;
385 22 : std::vector<double> vThisLineVal;
386 22 : std::vector<GByte> vResult;
387 22 : std::vector<double> vHeightResult;
388 :
389 : try
390 : {
391 11 : vFirstLineVal.resize(nXSize);
392 11 : vLastLineVal.resize(nXSize);
393 11 : vThisLineVal.resize(nXSize);
394 11 : vResult.resize(nXSize);
395 :
396 11 : if (oOpts.outputMode != OutputMode::Normal)
397 3 : vHeightResult.resize(nXSize);
398 : }
399 0 : catch (...)
400 : {
401 0 : CPLError(CE_Failure, CPLE_AppDefined,
402 : "Cannot allocate vectors for viewshed");
403 0 : return false;
404 : }
405 :
406 11 : double *padfFirstLineVal = vFirstLineVal.data();
407 11 : double *padfLastLineVal = vLastLineVal.data();
408 11 : double *padfThisLineVal = vThisLineVal.data();
409 11 : GByte *pabyResult = vResult.data();
410 11 : double *dfHeightResult = vHeightResult.data();
411 :
412 11 : GDALDriverManager *hMgr = GetGDALDriverManager();
413 11 : GDALDriver *hDriver = hMgr->GetDriverByName(oOpts.outputFormat.c_str());
414 11 : if (!hDriver)
415 : {
416 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot get driver");
417 1 : return false;
418 : }
419 :
420 : /* create output raster */
421 10 : poDstDS.reset(hDriver->Create(
422 : oOpts.outputFilename.c_str(), nXSize, nYSize, 1,
423 10 : oOpts.outputMode == OutputMode::Normal ? GDT_Byte : GDT_Float64,
424 10 : const_cast<char **>(oOpts.creationOpts.List())));
425 10 : if (!poDstDS)
426 : {
427 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create dataset for %s",
428 : oOpts.outputFilename.c_str());
429 1 : return false;
430 : }
431 : /* copy srs */
432 9 : if (hSrcDS)
433 18 : poDstDS->SetSpatialRef(
434 9 : GDALDataset::FromHandle(hSrcDS)->GetSpatialRef());
435 :
436 : std::array<double, 6> adfDstGeoTransform;
437 9 : adfDstGeoTransform[0] = adfGeoTransform[0] + adfGeoTransform[1] * nXStart +
438 9 : adfGeoTransform[2] * nYStart;
439 9 : adfDstGeoTransform[1] = adfGeoTransform[1];
440 9 : adfDstGeoTransform[2] = adfGeoTransform[2];
441 9 : adfDstGeoTransform[3] = adfGeoTransform[3] + adfGeoTransform[4] * nXStart +
442 9 : adfGeoTransform[5] * nYStart;
443 9 : adfDstGeoTransform[4] = adfGeoTransform[4];
444 9 : adfDstGeoTransform[5] = adfGeoTransform[5];
445 9 : poDstDS->SetGeoTransform(adfDstGeoTransform.data());
446 :
447 9 : auto hTargetBand = poDstDS->GetRasterBand(1);
448 9 : if (hTargetBand == nullptr)
449 : {
450 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot get band for %s",
451 : oOpts.outputFilename.c_str());
452 0 : return false;
453 : }
454 :
455 9 : if (oOpts.nodataVal >= 0)
456 1 : GDALSetRasterNoDataValue(hTargetBand, oOpts.nodataVal);
457 :
458 : /* process first line */
459 9 : if (GDALRasterIO(hBand, GF_Read, nXStart, nY, nXSize, 1, padfFirstLineVal,
460 9 : nXSize, 1, GDT_Float64, 0, 0))
461 : {
462 0 : CPLError(
463 : CE_Failure, CPLE_AppDefined,
464 : "RasterIO error when reading DEM at position(%d, %d), size(%d, %d)",
465 : nXStart, nY, nXSize, 1);
466 0 : return false;
467 : }
468 :
469 9 : const double dfZObserver = oOpts.observer.z + padfFirstLineVal[nX];
470 9 : double dfZ = 0.0;
471 9 : const double dfDistance2 = oOpts.maxDistance * oOpts.maxDistance;
472 :
473 : /* If we can't get a SemiMajor axis from the SRS, it will be
474 : * SRS_WGS84_SEMIMAJOR
475 : */
476 9 : double dfSphereDiameter(std::numeric_limits<double>::infinity());
477 9 : const OGRSpatialReference *poDstSRS = poDstDS->GetSpatialRef();
478 9 : if (poDstSRS)
479 : {
480 : OGRErr eSRSerr;
481 7 : double dfSemiMajor = poDstSRS->GetSemiMajor(&eSRSerr);
482 :
483 : /* If we fetched the axis from the SRS, use it */
484 7 : if (eSRSerr != OGRERR_FAILURE)
485 7 : dfSphereDiameter = dfSemiMajor * 2.0;
486 : else
487 0 : CPLDebug("GDALViewshedGenerate",
488 : "Unable to fetch SemiMajor axis from spatial reference");
489 : }
490 :
491 : /* mark the observer point as visible */
492 9 : double dfGroundLevel = 0;
493 9 : if (oOpts.outputMode == OutputMode::DEM)
494 1 : dfGroundLevel = padfFirstLineVal[nX];
495 :
496 9 : pabyResult[nX] = oOpts.visibleVal;
497 :
498 : //ABELL - Do we care about this conditional?
499 9 : if (oOpts.outputMode != OutputMode::Normal)
500 3 : dfHeightResult[nX] = dfGroundLevel;
501 :
502 9 : dfGroundLevel = 0;
503 9 : if (nX > 0)
504 : {
505 9 : if (oOpts.outputMode == OutputMode::DEM)
506 1 : dfGroundLevel = padfFirstLineVal[nX - 1];
507 9 : CPL_IGNORE_RET_VAL(AdjustHeightInRange(
508 9 : adfGeoTransform.data(), 1, 0, padfFirstLineVal[nX - 1], dfDistance2,
509 : oOpts.curveCoeff, dfSphereDiameter));
510 9 : pabyResult[nX - 1] = oOpts.visibleVal;
511 9 : if (oOpts.outputMode != OutputMode::Normal)
512 3 : dfHeightResult[nX - 1] = dfGroundLevel;
513 : }
514 9 : if (nX < nXSize - 1)
515 : {
516 9 : if (oOpts.outputMode == OutputMode::DEM)
517 1 : dfGroundLevel = padfFirstLineVal[nX + 1];
518 9 : CPL_IGNORE_RET_VAL(AdjustHeightInRange(
519 9 : adfGeoTransform.data(), 1, 0, padfFirstLineVal[nX + 1], dfDistance2,
520 : oOpts.curveCoeff, dfSphereDiameter));
521 9 : pabyResult[nX + 1] = oOpts.visibleVal;
522 9 : if (oOpts.outputMode != OutputMode::Normal)
523 3 : dfHeightResult[nX + 1] = dfGroundLevel;
524 : }
525 :
526 : /* process left direction */
527 337 : for (int iPixel = nX - 2; iPixel >= 0; iPixel--)
528 : {
529 328 : dfGroundLevel = 0;
530 328 : if (oOpts.outputMode == OutputMode::DEM)
531 50 : dfGroundLevel = padfFirstLineVal[iPixel];
532 :
533 328 : if (AdjustHeightInRange(adfGeoTransform.data(), nX - iPixel, 0,
534 328 : padfFirstLineVal[iPixel], dfDistance2,
535 : oOpts.curveCoeff, dfSphereDiameter))
536 : {
537 327 : dfZ = CalcHeightLine(nX - iPixel, padfFirstLineVal[iPixel + 1],
538 : dfZObserver);
539 :
540 327 : if (oOpts.outputMode != OutputMode::Normal)
541 150 : dfHeightResult[iPixel] = std::max(
542 150 : 0.0, (dfZ - padfFirstLineVal[iPixel] + dfGroundLevel));
543 :
544 327 : SetVisibility(iPixel, dfZ, oOpts.targetHeight, padfFirstLineVal,
545 327 : vResult, oOpts.visibleVal, oOpts.invisibleVal);
546 : }
547 : else
548 : {
549 2 : for (; iPixel >= 0; iPixel--)
550 : {
551 1 : pabyResult[iPixel] = oOpts.outOfRangeVal;
552 1 : if (oOpts.outputMode != OutputMode::Normal)
553 0 : dfHeightResult[iPixel] = oOpts.outOfRangeVal;
554 : }
555 : }
556 : }
557 : /* process right direction */
558 337 : for (int iPixel = nX + 2; iPixel < nXSize; iPixel++)
559 : {
560 328 : dfGroundLevel = 0;
561 328 : if (oOpts.outputMode == OutputMode::DEM)
562 50 : dfGroundLevel = padfFirstLineVal[iPixel];
563 328 : if (AdjustHeightInRange(adfGeoTransform.data(), iPixel - nX, 0,
564 328 : padfFirstLineVal[iPixel], dfDistance2,
565 : oOpts.curveCoeff, dfSphereDiameter))
566 : {
567 327 : dfZ = CalcHeightLine(iPixel - nX, padfFirstLineVal[iPixel - 1],
568 : dfZObserver);
569 :
570 327 : if (oOpts.outputMode != OutputMode::Normal)
571 150 : dfHeightResult[iPixel] = std::max(
572 150 : 0.0, (dfZ - padfFirstLineVal[iPixel] + dfGroundLevel));
573 :
574 327 : SetVisibility(iPixel, dfZ, oOpts.targetHeight, padfFirstLineVal,
575 327 : vResult, oOpts.visibleVal, oOpts.invisibleVal);
576 : }
577 : else
578 : {
579 2 : for (; iPixel < nXSize; iPixel++)
580 : {
581 1 : pabyResult[iPixel] = oOpts.outOfRangeVal;
582 1 : if (oOpts.outputMode != OutputMode::Normal)
583 0 : dfHeightResult[iPixel] = oOpts.outOfRangeVal;
584 : }
585 : }
586 : }
587 : /* write result line */
588 :
589 : void *data;
590 : GDALDataType dataType;
591 9 : if (oOpts.outputMode == OutputMode::Normal)
592 : {
593 6 : data = static_cast<void *>(pabyResult);
594 6 : dataType = GDT_Byte;
595 : }
596 : else
597 : {
598 3 : data = static_cast<void *>(dfHeightResult);
599 3 : dataType = GDT_Float64;
600 : }
601 9 : if (GDALRasterIO(hTargetBand, GF_Write, 0, nY - nYStart, nXSize, 1, data,
602 9 : nXSize, 1, dataType, 0, 0))
603 : {
604 0 : CPLError(CE_Failure, CPLE_AppDefined,
605 : "RasterIO error when writing target raster at position "
606 : "(%d,%d), size (%d,%d)",
607 : 0, nY - nYStart, nXSize, 1);
608 0 : return false;
609 : }
610 :
611 : /* scan upwards */
612 9 : std::copy(vFirstLineVal.begin(), vFirstLineVal.end(), vLastLineVal.begin());
613 458 : for (int iLine = nY - 1; iLine >= nYStart; iLine--)
614 : {
615 449 : if (GDALRasterIO(hBand, GF_Read, nXStart, iLine, nXSize, 1,
616 449 : padfThisLineVal, nXSize, 1, GDT_Float64, 0, 0))
617 : {
618 0 : CPLError(CE_Failure, CPLE_AppDefined,
619 : "RasterIO error when reading DEM at position (%d,%d), "
620 : "size (%d,%d)",
621 : nXStart, iLine, nXSize, 1);
622 0 : return false;
623 : }
624 :
625 : /* set up initial point on the scanline */
626 449 : dfGroundLevel = 0;
627 449 : if (oOpts.outputMode == OutputMode::DEM)
628 70 : dfGroundLevel = padfThisLineVal[nX];
629 449 : if (AdjustHeightInRange(adfGeoTransform.data(), 0, nY - iLine,
630 449 : padfThisLineVal[nX], dfDistance2,
631 : oOpts.curveCoeff, dfSphereDiameter))
632 : {
633 448 : dfZ = CalcHeightLine(nY - iLine, padfLastLineVal[nX], dfZObserver);
634 :
635 448 : if (oOpts.outputMode != OutputMode::Normal)
636 210 : dfHeightResult[nX] =
637 210 : std::max(0.0, (dfZ - padfThisLineVal[nX] + dfGroundLevel));
638 :
639 448 : SetVisibility(nX, dfZ, oOpts.targetHeight, padfThisLineVal, vResult,
640 448 : oOpts.visibleVal, oOpts.invisibleVal);
641 : }
642 : else
643 : {
644 1 : pabyResult[nX] = oOpts.outOfRangeVal;
645 1 : if (oOpts.outputMode != OutputMode::Normal)
646 0 : dfHeightResult[nX] = oOpts.outOfRangeVal;
647 : }
648 :
649 : /* process left direction */
650 22361 : for (int iPixel = nX - 1; iPixel >= 0; iPixel--)
651 : {
652 21912 : dfGroundLevel = 0;
653 21912 : if (oOpts.outputMode == OutputMode::DEM)
654 3570 : dfGroundLevel = padfThisLineVal[iPixel];
655 21912 : if (AdjustHeightInRange(adfGeoTransform.data(), nX - iPixel,
656 21912 : nY - iLine, padfThisLineVal[iPixel],
657 : dfDistance2, oOpts.curveCoeff,
658 : dfSphereDiameter))
659 : {
660 21887 : if (oOpts.cellMode != CellMode::Edge)
661 0 : dfZ = CalcHeightDiagonal(
662 0 : nX - iPixel, nY - iLine, padfThisLineVal[iPixel + 1],
663 0 : padfLastLineVal[iPixel], dfZObserver);
664 :
665 21887 : if (oOpts.cellMode != CellMode::Diagonal)
666 : {
667 : double dfZ2 =
668 21887 : nX - iPixel >= nY - iLine
669 21887 : ? CalcHeightEdge(nY - iLine, nX - iPixel,
670 8202 : padfLastLineVal[iPixel + 1],
671 8202 : padfThisLineVal[iPixel + 1],
672 : dfZObserver)
673 13685 : : CalcHeightEdge(nX - iPixel, nY - iLine,
674 13685 : padfLastLineVal[iPixel + 1],
675 13685 : padfLastLineVal[iPixel],
676 21887 : dfZObserver);
677 21887 : dfZ = CalcHeight(dfZ, dfZ2, oOpts.cellMode);
678 : }
679 :
680 21887 : if (oOpts.outputMode != OutputMode::Normal)
681 10710 : dfHeightResult[iPixel] = std::max(
682 10710 : 0.0, (dfZ - padfThisLineVal[iPixel] + dfGroundLevel));
683 :
684 21887 : SetVisibility(iPixel, dfZ, oOpts.targetHeight, padfThisLineVal,
685 21887 : vResult, oOpts.visibleVal, oOpts.invisibleVal);
686 : }
687 : else
688 : {
689 195 : for (; iPixel >= 0; iPixel--)
690 : {
691 170 : pabyResult[iPixel] = oOpts.outOfRangeVal;
692 170 : if (oOpts.outputMode != OutputMode::Normal)
693 0 : dfHeightResult[iPixel] = oOpts.outOfRangeVal;
694 : }
695 : }
696 : }
697 : /* process right direction */
698 22361 : for (int iPixel = nX + 1; iPixel < nXSize; iPixel++)
699 : {
700 21912 : dfGroundLevel = 0;
701 21912 : if (oOpts.outputMode == OutputMode::DEM)
702 3570 : dfGroundLevel = padfThisLineVal[iPixel];
703 :
704 21912 : if (AdjustHeightInRange(adfGeoTransform.data(), iPixel - nX,
705 21912 : nY - iLine, padfThisLineVal[iPixel],
706 : dfDistance2, oOpts.curveCoeff,
707 : dfSphereDiameter))
708 : {
709 21887 : if (oOpts.cellMode != CellMode::Edge)
710 0 : dfZ = CalcHeightDiagonal(
711 0 : iPixel - nX, nY - iLine, padfThisLineVal[iPixel - 1],
712 0 : padfLastLineVal[iPixel], dfZObserver);
713 21887 : if (oOpts.cellMode != CellMode::Diagonal)
714 : {
715 : double dfZ2 =
716 21887 : iPixel - nX >= nY - iLine
717 21887 : ? CalcHeightEdge(nY - iLine, iPixel - nX,
718 8202 : padfLastLineVal[iPixel - 1],
719 8202 : padfThisLineVal[iPixel - 1],
720 : dfZObserver)
721 13685 : : CalcHeightEdge(iPixel - nX, nY - iLine,
722 13685 : padfLastLineVal[iPixel - 1],
723 13685 : padfLastLineVal[iPixel],
724 21887 : dfZObserver);
725 21887 : dfZ = CalcHeight(dfZ, dfZ2, oOpts.cellMode);
726 : }
727 :
728 21887 : if (oOpts.outputMode != OutputMode::Normal)
729 10710 : dfHeightResult[iPixel] = std::max(
730 10710 : 0.0, (dfZ - padfThisLineVal[iPixel] + dfGroundLevel));
731 :
732 21887 : SetVisibility(iPixel, dfZ, oOpts.targetHeight, padfThisLineVal,
733 21887 : vResult, oOpts.visibleVal, oOpts.invisibleVal);
734 : }
735 : else
736 : {
737 195 : for (; iPixel < nXSize; iPixel++)
738 : {
739 170 : pabyResult[iPixel] = oOpts.outOfRangeVal;
740 170 : if (oOpts.outputMode != OutputMode::Normal)
741 0 : dfHeightResult[iPixel] = oOpts.outOfRangeVal;
742 : }
743 : }
744 : }
745 :
746 : /* write result line */
747 449 : if (GDALRasterIO(hTargetBand, GF_Write, 0, iLine - nYStart, nXSize, 1,
748 449 : data, nXSize, 1, dataType, 0, 0))
749 : {
750 0 : CPLError(CE_Failure, CPLE_AppDefined,
751 : "RasterIO error when writing target raster at position "
752 : "(%d,%d), size (%d,%d)",
753 : 0, iLine - nYStart, nXSize, 1);
754 0 : return false;
755 : }
756 :
757 449 : std::swap(padfLastLineVal, padfThisLineVal);
758 :
759 449 : if (!pfnProgress((nY - iLine) / static_cast<double>(nYSize), "",
760 : pProgressArg))
761 : {
762 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
763 0 : return false;
764 : }
765 : }
766 :
767 : /* scan downwards */
768 9 : memcpy(padfLastLineVal, padfFirstLineVal, nXSize * sizeof(double));
769 452 : for (int iLine = nY + 1; iLine < nYStop; iLine++)
770 : {
771 443 : if (GDALRasterIO(hBand, GF_Read, nXStart, iLine, nXSize, 1,
772 443 : padfThisLineVal, nXSize, 1, GDT_Float64, 0, 0))
773 : {
774 0 : CPLError(CE_Failure, CPLE_AppDefined,
775 : "RasterIO error when reading DEM at position (%d,%d), "
776 : "size (%d,%d)",
777 : nXStart, iLine, nXStop - nXStart, 1);
778 0 : return false;
779 : }
780 :
781 : /* set up initial point on the scanline */
782 443 : dfGroundLevel = 0;
783 443 : if (oOpts.outputMode == OutputMode::DEM)
784 69 : dfGroundLevel = padfThisLineVal[nX];
785 :
786 443 : if (AdjustHeightInRange(adfGeoTransform.data(), 0, iLine - nY,
787 443 : padfThisLineVal[nX], dfDistance2,
788 : oOpts.curveCoeff, dfSphereDiameter))
789 : {
790 442 : dfZ = CalcHeightLine(iLine - nY, padfLastLineVal[nX], dfZObserver);
791 :
792 442 : if (oOpts.outputMode != OutputMode::Normal)
793 207 : dfHeightResult[nX] =
794 207 : std::max(0.0, (dfZ - padfThisLineVal[nX] + dfGroundLevel));
795 :
796 442 : SetVisibility(nX, dfZ, oOpts.targetHeight, padfThisLineVal, vResult,
797 442 : oOpts.visibleVal, oOpts.invisibleVal);
798 : }
799 : else
800 : {
801 1 : pabyResult[nX] = oOpts.outOfRangeVal;
802 1 : if (oOpts.outputMode != OutputMode::Normal)
803 0 : dfHeightResult[nX] = oOpts.outOfRangeVal;
804 : }
805 :
806 : /* process left direction */
807 22049 : for (int iPixel = nX - 1; iPixel >= 0; iPixel--)
808 : {
809 21606 : dfGroundLevel = 0;
810 21606 : if (oOpts.outputMode == OutputMode::DEM)
811 3519 : dfGroundLevel = padfThisLineVal[iPixel];
812 :
813 21606 : if (AdjustHeightInRange(adfGeoTransform.data(), nX - iPixel,
814 21606 : iLine - nY, padfThisLineVal[iPixel],
815 : dfDistance2, oOpts.curveCoeff,
816 : dfSphereDiameter))
817 : {
818 21581 : if (oOpts.cellMode != CellMode::Edge)
819 0 : dfZ = CalcHeightDiagonal(
820 0 : nX - iPixel, iLine - nY, padfThisLineVal[iPixel + 1],
821 0 : padfLastLineVal[iPixel], dfZObserver);
822 :
823 21581 : if (oOpts.cellMode != CellMode::Diagonal)
824 : {
825 : double dfZ2 =
826 21581 : nX - iPixel >= iLine - nY
827 21581 : ? CalcHeightEdge(iLine - nY, nX - iPixel,
828 8202 : padfLastLineVal[iPixel + 1],
829 8202 : padfThisLineVal[iPixel + 1],
830 : dfZObserver)
831 13379 : : CalcHeightEdge(nX - iPixel, iLine - nY,
832 13379 : padfLastLineVal[iPixel + 1],
833 13379 : padfLastLineVal[iPixel],
834 21581 : dfZObserver);
835 21581 : dfZ = CalcHeight(dfZ, dfZ2, oOpts.cellMode);
836 : }
837 :
838 21581 : if (oOpts.outputMode != OutputMode::Normal)
839 10557 : dfHeightResult[iPixel] = std::max(
840 10557 : 0.0, (dfZ - padfThisLineVal[iPixel] + dfGroundLevel));
841 :
842 21581 : SetVisibility(iPixel, dfZ, oOpts.targetHeight, padfThisLineVal,
843 21581 : vResult, oOpts.visibleVal, oOpts.invisibleVal);
844 : }
845 : else
846 : {
847 195 : for (; iPixel >= 0; iPixel--)
848 : {
849 170 : pabyResult[iPixel] = oOpts.outOfRangeVal;
850 170 : if (oOpts.outputMode != OutputMode::Normal)
851 0 : dfHeightResult[iPixel] = oOpts.outOfRangeVal;
852 : }
853 : }
854 : }
855 : /* process right direction */
856 22049 : for (int iPixel = nX + 1; iPixel < nXSize; iPixel++)
857 : {
858 21606 : dfGroundLevel = 0;
859 21606 : if (oOpts.outputMode == OutputMode::DEM)
860 3519 : dfGroundLevel = padfThisLineVal[iPixel];
861 :
862 21606 : if (AdjustHeightInRange(adfGeoTransform.data(), iPixel - nX,
863 21606 : iLine - nY, padfThisLineVal[iPixel],
864 : dfDistance2, oOpts.curveCoeff,
865 : dfSphereDiameter))
866 : {
867 21581 : if (oOpts.cellMode != CellMode::Edge)
868 0 : dfZ = CalcHeightDiagonal(
869 0 : iPixel - nX, iLine - nY, padfThisLineVal[iPixel - 1],
870 0 : padfLastLineVal[iPixel], dfZObserver);
871 :
872 21581 : if (oOpts.cellMode != CellMode::Diagonal)
873 : {
874 : double dfZ2 =
875 21581 : iPixel - nX >= iLine - nY
876 21581 : ? CalcHeightEdge(iLine - nY, iPixel - nX,
877 8202 : padfLastLineVal[iPixel - 1],
878 8202 : padfThisLineVal[iPixel - 1],
879 : dfZObserver)
880 13379 : : CalcHeightEdge(iPixel - nX, iLine - nY,
881 13379 : padfLastLineVal[iPixel - 1],
882 13379 : padfLastLineVal[iPixel],
883 21581 : dfZObserver);
884 21581 : dfZ = CalcHeight(dfZ, dfZ2, oOpts.cellMode);
885 : }
886 :
887 21581 : if (oOpts.outputMode != OutputMode::Normal)
888 10557 : dfHeightResult[iPixel] = std::max(
889 10557 : 0.0, (dfZ - padfThisLineVal[iPixel] + dfGroundLevel));
890 :
891 21581 : SetVisibility(iPixel, dfZ, oOpts.targetHeight, padfThisLineVal,
892 21581 : vResult, oOpts.visibleVal, oOpts.invisibleVal);
893 : }
894 : else
895 : {
896 195 : for (; iPixel < nXSize; iPixel++)
897 : {
898 170 : pabyResult[iPixel] = oOpts.outOfRangeVal;
899 170 : if (oOpts.outputMode != OutputMode::Normal)
900 0 : dfHeightResult[iPixel] = oOpts.outOfRangeVal;
901 : }
902 : }
903 : }
904 :
905 : /* write result line */
906 443 : if (GDALRasterIO(hTargetBand, GF_Write, 0, iLine - nYStart, nXSize, 1,
907 443 : data, nXSize, 1, dataType, 0, 0))
908 : {
909 0 : CPLError(CE_Failure, CPLE_AppDefined,
910 : "RasterIO error when writing target raster at position "
911 : "(%d,%d), size (%d,%d)",
912 : 0, iLine - nYStart, nXSize, 1);
913 0 : return false;
914 : }
915 :
916 443 : std::swap(padfLastLineVal, padfThisLineVal);
917 :
918 443 : if (!pfnProgress((iLine - nYStart) / static_cast<double>(nYSize), "",
919 : pProgressArg))
920 : {
921 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
922 0 : return false;
923 : }
924 : }
925 :
926 9 : if (!pfnProgress(1.0, "", pProgressArg))
927 : {
928 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
929 0 : return false;
930 : }
931 :
932 9 : return true;
933 : }
934 :
935 : } // namespace gdal
|
