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