Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL DEM Utilities
4 : * Purpose:
5 : * Authors: Matthew Perry, perrygeo at gmail.com
6 : * Even Rouault, even dot rouault at spatialys.com
7 : * Howard Butler, hobu.inc at gmail.com
8 : * Chris Yesson, chris dot yesson at ioz dot ac dot uk
9 : *
10 : ******************************************************************************
11 : * Copyright (c) 2006, 2009 Matthew Perry
12 : * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>
13 : * Portions derived from GRASS 4.1 (public domain) See
14 : * http://trac.osgeo.org/gdal/ticket/2975 for more information regarding
15 : * history of this code
16 : *
17 : * SPDX-License-Identifier: MIT
18 : ****************************************************************************
19 : *
20 : * Slope and aspect calculations based on original method for GRASS GIS 4.1
21 : * by Michael Shapiro, U.S.Army Construction Engineering Research Laboratory
22 : * Olga Waupotitsch, U.S.Army Construction Engineering Research Laboratory
23 : * Marjorie Larson, U.S.Army Construction Engineering Research Laboratory
24 : * as found in GRASS's r.slope.aspect module.
25 : *
26 : * Horn's formula is used to find the first order derivatives in x and y
27 : *directions for slope and aspect calculations: Horn, B. K. P. (1981). "Hill
28 : *Shading and the Reflectance Map", Proceedings of the IEEE, 69(1):14-47.
29 : *
30 : * Other reference :
31 : * Burrough, P.A. and McDonell, R.A., 1998. Principles of Geographical
32 : *Information Systems. p. 190.
33 : *
34 : * Shaded relief based on original method for GRASS GIS 4.1 by Jim Westervelt,
35 : * U.S. Army Construction Engineering Research Laboratory
36 : * as found in GRASS's r.shaded.relief (formerly shade.rel.sh) module.
37 : * ref: "r.mapcalc: An Algebra for GIS and Image Processing",
38 : * by Michael Shapiro and Jim Westervelt, U.S. Army Construction Engineering
39 : * Research Laboratory (March/1991)
40 : *
41 : * Color table of named colors and lookup code derived from
42 : *src/libes/gis/named_colr.c of GRASS 4.1
43 : *
44 : * TRI -
45 : * For bathymetric use cases, implements
46 : * Terrain Ruggedness Index is as described in Wilson et al. (2007)
47 : * this is based on the method of Valentine et al. (2004)
48 : *
49 : * For terrestrial use cases, implements
50 : * Riley, S.J., De Gloria, S.D., Elliot, R. (1999): A Terrain Ruggedness
51 : * that Quantifies Topographic Heterogeneity. Intermountain Journal of Science,
52 : *Vol.5, No.1-4, pp.23-27
53 : *
54 : *
55 : * TPI - Topographic Position Index follows the description in
56 : * Wilson et al. (2007), following Weiss (2001). The radius is fixed
57 : * at 1 cell width/height
58 : *
59 : * Roughness - follows the definition in Wilson et al. (2007), which follows
60 : * Dartnell (2000).
61 : *
62 : * References for TRI/TPI/Roughness:
63 : * Dartnell, P. 2000. Applying Remote Sensing Techniques to map Seafloor
64 : * Geology/Habitat Relationships. Masters Thesis, San Francisco State
65 : * University, pp. 108.
66 : * Valentine, P. C., S. J. Fuller, L. A. Scully. 2004. Terrain Ruggedness
67 : * Analysis and Distribution of Boulder Ridges in the Stellwagen Bank National
68 : * Marine Sanctuary Region (poster). Galway, Ireland: 5th International
69 : * Symposium on Marine Geological and Biological Habitat Mapping (GeoHAB),
70 : * May 2004.
71 : * Weiss, A. D. 2001. Topographic Positions and Landforms Analysis (poster),
72 : * ESRI International User Conference, July 2001. San Diego, CA: ESRI.
73 : * Wilson, M. F. J.; O'Connell, B.; Brown, C.; Guinan, J. C. & Grehan, A. J.
74 : * Multiscale terrain analysis of multibeam bathymetry data for habitat mapping
75 : * on the continental slope Marine Geodesy, 2007, 30, 3-35
76 : ****************************************************************************/
77 :
78 : // Include before others for mingw for VSIStatBufL
79 : #include "cpl_conv.h"
80 :
81 : #include "cpl_port.h"
82 : #include "gdal_utils.h"
83 : #include "gdal_utils_priv.h"
84 : #include "commonutils.h"
85 : #include "gdalargumentparser.h"
86 :
87 : #include <cassert>
88 : #include <cfloat>
89 : #include <cmath>
90 : #include <cstdio>
91 : #include <cstdlib>
92 : #include <cstring>
93 : #if HAVE_SYS_STAT_H
94 : #include <sys/stat.h>
95 : #endif
96 :
97 : #include <algorithm>
98 : #include <array>
99 : #include <cmath>
100 : #include <limits>
101 :
102 : #include "cpl_error.h"
103 : #include "cpl_float.h"
104 : #include "cpl_progress.h"
105 : #include "cpl_string.h"
106 : #include "cpl_vsi.h"
107 : #include "cpl_vsi_virtual.h"
108 : #include "gdal.h"
109 : #include "gdal_priv.h"
110 :
111 : #if defined(__x86_64__) || defined(_M_X64)
112 : #define HAVE_16_SSE_REG
113 : #include "emmintrin.h"
114 : #include "gdalsse_priv.h"
115 : #endif
116 :
117 : static const double kdfDegreesToRadians = M_PI / 180.0;
118 : static const double kdfRadiansToDegrees = 180.0 / M_PI;
119 :
120 : typedef enum
121 : {
122 : COLOR_SELECTION_INTERPOLATE,
123 : COLOR_SELECTION_NEAREST_ENTRY,
124 : COLOR_SELECTION_EXACT_ENTRY
125 : } ColorSelectionMode;
126 :
127 : namespace gdal::GDALDEM
128 : {
129 : enum class GradientAlg
130 : {
131 : HORN,
132 : ZEVENBERGEN_THORNE,
133 : };
134 :
135 : enum class TRIAlg
136 : {
137 : WILSON,
138 : RILEY,
139 : };
140 : } // namespace gdal::GDALDEM
141 :
142 : using namespace gdal::GDALDEM;
143 :
144 : struct GDALDEMProcessingOptions
145 : {
146 : /*! output format. Use the short format name. */
147 : std::string osFormat{};
148 :
149 : /*! the progress function to use */
150 : GDALProgressFunc pfnProgress = nullptr;
151 :
152 : /*! pointer to the progress data variable */
153 : void *pProgressData = nullptr;
154 :
155 : double z = 1.0;
156 : double globalScale = std::numeric_limits<
157 : double>::quiet_NaN(); // when set, copied to xscale and yscale
158 : double xscale = std::numeric_limits<double>::quiet_NaN();
159 : double yscale = std::numeric_limits<double>::quiet_NaN();
160 : double az = 315.0;
161 : double alt = 45.0;
162 : bool bSlopeFormatUseDegrees =
163 : true; // false = 'percent' or true = 'degrees'
164 : bool bAddAlpha = false;
165 : bool bZeroForFlat = false;
166 : bool bAngleAsAzimuth = true;
167 : ColorSelectionMode eColorSelectionMode = COLOR_SELECTION_INTERPOLATE;
168 : bool bComputeAtEdges = false;
169 : bool bGradientAlgSpecified = false;
170 : GradientAlg eGradientAlg = GradientAlg::HORN;
171 : bool bTRIAlgSpecified = false;
172 : TRIAlg eTRIAlg = TRIAlg::RILEY;
173 : bool bCombined = false;
174 : bool bIgor = false;
175 : bool bMultiDirectional = false;
176 : CPLStringList aosCreationOptions{};
177 : int nBand = 1;
178 : };
179 :
180 : /************************************************************************/
181 : /* AlgorithmParameters */
182 : /************************************************************************/
183 :
184 106 : struct AlgorithmParameters
185 : {
186 100 : AlgorithmParameters() = default;
187 : virtual ~AlgorithmParameters();
188 6 : AlgorithmParameters(const AlgorithmParameters &) = default;
189 : AlgorithmParameters &operator=(const AlgorithmParameters &) = delete;
190 : AlgorithmParameters(AlgorithmParameters &&) = delete;
191 : AlgorithmParameters &operator=(AlgorithmParameters &&) = delete;
192 :
193 : virtual std::unique_ptr<AlgorithmParameters>
194 : CreateScaledParameters(double dfXRatio, double dfYRatio) = 0;
195 : };
196 :
197 : AlgorithmParameters::~AlgorithmParameters() = default;
198 :
199 : /************************************************************************/
200 : /* ComputeVal() */
201 : /************************************************************************/
202 :
203 : template <class T> struct GDALGeneric3x3ProcessingAlg
204 : {
205 : typedef float (*type)(const T *pafWindow, float fDstNoDataValue,
206 : const AlgorithmParameters *pData);
207 : };
208 :
209 : template <class T> struct GDALGeneric3x3ProcessingAlg_multisample
210 : {
211 : typedef int (*type)(const T *pafFirstLine, const T *pafSecondLine,
212 : const T *pafThirdLine, int nXSize,
213 : const AlgorithmParameters *pData, float *pafOutputBuf);
214 : };
215 :
216 : template <class T>
217 : static float ComputeVal(bool bSrcHasNoData, T fSrcNoDataValue,
218 : bool bIsSrcNoDataNan, T *afWin, float fDstNoDataValue,
219 : typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlg,
220 : const AlgorithmParameters *pData, bool bComputeAtEdges);
221 :
222 : template <>
223 111005 : float ComputeVal(bool bSrcHasNoData, float fSrcNoDataValue,
224 : bool bIsSrcNoDataNan, float *afWin, float fDstNoDataValue,
225 : GDALGeneric3x3ProcessingAlg<float>::type pfnAlg,
226 : const AlgorithmParameters *pData, bool bComputeAtEdges)
227 : {
228 164414 : if (bSrcHasNoData &&
229 53409 : ((!bIsSrcNoDataNan && ARE_REAL_EQUAL(afWin[4], fSrcNoDataValue)) ||
230 0 : (bIsSrcNoDataNan && std::isnan(afWin[4]))))
231 : {
232 8 : return fDstNoDataValue;
233 : }
234 110997 : else if (bSrcHasNoData)
235 : {
236 527225 : for (int k = 0; k < 9; k++)
237 : {
238 474949 : if ((!bIsSrcNoDataNan &&
239 949898 : ARE_REAL_EQUAL(afWin[k], fSrcNoDataValue)) ||
240 0 : (bIsSrcNoDataNan && std::isnan(afWin[k])))
241 : {
242 1133 : if (bComputeAtEdges)
243 8 : afWin[k] = afWin[4];
244 : else
245 1125 : return fDstNoDataValue;
246 : }
247 : }
248 : }
249 :
250 109872 : return pfnAlg(afWin, fDstNoDataValue, pData);
251 : }
252 :
253 : template <>
254 1142590 : float ComputeVal(bool bSrcHasNoData, GInt32 fSrcNoDataValue,
255 : bool /* bIsSrcNoDataNan */, GInt32 *afWin,
256 : float fDstNoDataValue,
257 : GDALGeneric3x3ProcessingAlg<GInt32>::type pfnAlg,
258 : const AlgorithmParameters *pData, bool bComputeAtEdges)
259 : {
260 1142590 : if (bSrcHasNoData && afWin[4] == fSrcNoDataValue)
261 : {
262 584 : return fDstNoDataValue;
263 : }
264 1142010 : else if (bSrcHasNoData)
265 : {
266 5543190 : for (int k = 0; k < 9; k++)
267 : {
268 4989200 : if (afWin[k] == fSrcNoDataValue)
269 : {
270 852 : if (bComputeAtEdges)
271 8 : afWin[k] = afWin[4];
272 : else
273 844 : return fDstNoDataValue;
274 : }
275 : }
276 : }
277 :
278 1141160 : return pfnAlg(afWin, fDstNoDataValue, pData);
279 : }
280 :
281 : /************************************************************************/
282 : /* INTERPOL() */
283 : /************************************************************************/
284 :
285 : template <class T>
286 : static T INTERPOL(T a, T b, int bSrcHasNodata, T fSrcNoDataValue);
287 :
288 : template <>
289 1464 : float INTERPOL(float a, float b, int bSrcHasNoData, float fSrcNoDataValue)
290 : {
291 2904 : if (bSrcHasNoData && (ARE_REAL_EQUAL(a, fSrcNoDataValue) ||
292 1440 : ARE_REAL_EQUAL(b, fSrcNoDataValue)))
293 24 : return fSrcNoDataValue;
294 1440 : const float fVal = 2 * a - b;
295 1440 : if (bSrcHasNoData && ARE_REAL_EQUAL(fVal, fSrcNoDataValue))
296 : return fSrcNoDataValue *
297 0 : (1 + 3 * std::numeric_limits<float>::epsilon());
298 1440 : return fVal;
299 : }
300 :
301 : template <>
302 69144 : GInt32 INTERPOL(GInt32 a, GInt32 b, int bSrcHasNoData, GInt32 fSrcNoDataValue)
303 : {
304 69144 : if (bSrcHasNoData && ((a == fSrcNoDataValue) || (b == fSrcNoDataValue)))
305 24 : return fSrcNoDataValue;
306 : const int nVal = static_cast<int>(
307 69120 : std::clamp<int64_t>(2 * static_cast<int64_t>(a) - b, INT_MIN, INT_MAX));
308 69120 : if (bSrcHasNoData && fSrcNoDataValue == nVal)
309 0 : return nVal == INT_MAX ? INT_MAX - 1 : nVal + 1;
310 69120 : return nVal;
311 : }
312 :
313 : /************************************************************************/
314 : /* GDALGeneric3x3Processing() */
315 : /************************************************************************/
316 :
317 : template <class T>
318 75 : static CPLErr GDALGeneric3x3Processing(
319 : GDALRasterBandH hSrcBand, GDALRasterBandH hDstBand,
320 : typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlg,
321 : typename GDALGeneric3x3ProcessingAlg_multisample<T>::type
322 : pfnAlg_multisample,
323 : std::unique_ptr<AlgorithmParameters> pData, bool bComputeAtEdges,
324 : GDALProgressFunc pfnProgress, void *pProgressData)
325 : {
326 75 : if (pfnProgress == nullptr)
327 69 : pfnProgress = GDALDummyProgress;
328 :
329 : /* -------------------------------------------------------------------- */
330 : /* Initialize progress counter. */
331 : /* -------------------------------------------------------------------- */
332 75 : if (!pfnProgress(0.0, nullptr, pProgressData))
333 : {
334 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
335 0 : return CE_Failure;
336 : }
337 :
338 75 : const int nXSize = GDALGetRasterBandXSize(hSrcBand);
339 75 : const int nYSize = GDALGetRasterBandYSize(hSrcBand);
340 :
341 : // 1 line destination buffer.
342 : float *pafOutputBuf =
343 75 : static_cast<float *>(VSI_MALLOC2_VERBOSE(sizeof(float), nXSize));
344 : // 3 line rotating source buffer.
345 : T *pafThreeLineWin =
346 75 : static_cast<T *>(VSI_MALLOC2_VERBOSE(3 * sizeof(T), nXSize));
347 75 : if (pafOutputBuf == nullptr || pafThreeLineWin == nullptr)
348 : {
349 0 : VSIFree(pafOutputBuf);
350 0 : VSIFree(pafThreeLineWin);
351 0 : return CE_Failure;
352 : }
353 :
354 : GDALDataType eReadDT;
355 75 : int bSrcHasNoData = FALSE;
356 : const double dfNoDataValue =
357 75 : GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
358 :
359 75 : bool bIsSrcNoDataNan = false;
360 75 : T fSrcNoDataValue = 0;
361 : if constexpr (std::numeric_limits<T>::is_integer)
362 : {
363 63 : eReadDT = GDT_Int32;
364 63 : if (bSrcHasNoData)
365 : {
366 61 : GDALDataType eSrcDT = GDALGetRasterDataType(hSrcBand);
367 61 : CPLAssert(eSrcDT == GDT_Byte || eSrcDT == GDT_UInt16 ||
368 : eSrcDT == GDT_Int16);
369 61 : const int nMinVal = (eSrcDT == GDT_Byte) ? 0
370 : : (eSrcDT == GDT_UInt16) ? 0
371 : : -32768;
372 61 : const int nMaxVal = (eSrcDT == GDT_Byte) ? 255
373 : : (eSrcDT == GDT_UInt16) ? 65535
374 : : 32767;
375 :
376 61 : if (fabs(dfNoDataValue - floor(dfNoDataValue + 0.5)) < 1e-2 &&
377 61 : dfNoDataValue >= nMinVal && dfNoDataValue <= nMaxVal)
378 : {
379 61 : fSrcNoDataValue = static_cast<T>(floor(dfNoDataValue + 0.5));
380 : }
381 : else
382 : {
383 0 : bSrcHasNoData = FALSE;
384 : }
385 : }
386 : }
387 : else
388 : {
389 12 : eReadDT = GDT_Float32;
390 12 : fSrcNoDataValue = static_cast<T>(dfNoDataValue);
391 12 : bIsSrcNoDataNan = bSrcHasNoData && std::isnan(dfNoDataValue);
392 : }
393 :
394 75 : int bDstHasNoData = FALSE;
395 75 : float fDstNoDataValue =
396 75 : static_cast<float>(GDALGetRasterNoDataValue(hDstBand, &bDstHasNoData));
397 75 : if (!bDstHasNoData)
398 0 : fDstNoDataValue = 0.0;
399 :
400 75 : int nLine1Off = 0;
401 75 : int nLine2Off = nXSize;
402 75 : int nLine3Off = 2 * nXSize;
403 :
404 : // Move a 3x3 pafWindow over each cell
405 : // (where the cell in question is #4)
406 : //
407 : // 0 1 2
408 : // 3 4 5
409 : // 6 7 8
410 :
411 : /* Preload the first 2 lines */
412 :
413 75 : bool abLineHasNoDataValue[3] = {CPL_TO_BOOL(bSrcHasNoData),
414 75 : CPL_TO_BOOL(bSrcHasNoData),
415 75 : CPL_TO_BOOL(bSrcHasNoData)};
416 :
417 225 : for (int i = 0; i < 2 && i < nYSize; i++)
418 : {
419 300 : if (GDALRasterIO(hSrcBand, GF_Read, 0, i, nXSize, 1,
420 150 : pafThreeLineWin + i * nXSize, nXSize, 1, eReadDT, 0,
421 150 : 0) != CE_None)
422 : {
423 0 : CPLFree(pafOutputBuf);
424 0 : CPLFree(pafThreeLineWin);
425 :
426 0 : return CE_Failure;
427 : }
428 150 : if (bSrcHasNoData)
429 : {
430 146 : abLineHasNoDataValue[i] = false;
431 : if constexpr (std::numeric_limits<T>::is_integer)
432 : {
433 11496 : for (int iX = 0; iX < nXSize; iX++)
434 : {
435 11402 : if (pafThreeLineWin[i * nXSize + iX] == fSrcNoDataValue)
436 : {
437 28 : abLineHasNoDataValue[i] = true;
438 28 : break;
439 : }
440 : }
441 : }
442 : else
443 : {
444 1476 : for (int iX = 0; iX < nXSize; iX++)
445 : {
446 2916 : if (pafThreeLineWin[i * nXSize + iX] == fSrcNoDataValue ||
447 1452 : std::isnan(pafThreeLineWin[i * nXSize + iX]))
448 : {
449 12 : abLineHasNoDataValue[i] = true;
450 12 : break;
451 : }
452 : }
453 : }
454 : }
455 : }
456 :
457 75 : CPLErr eErr = CE_None;
458 75 : if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
459 : {
460 2204 : for (int j = 0; j < nXSize; j++)
461 : {
462 2184 : int jmin = (j == 0) ? j : j - 1;
463 2184 : int jmax = (j == nXSize - 1) ? j : j + 1;
464 :
465 6552 : T afWin[9] = {
466 2184 : INTERPOL(pafThreeLineWin[jmin], pafThreeLineWin[nXSize + jmin],
467 : bSrcHasNoData, fSrcNoDataValue),
468 2184 : INTERPOL(pafThreeLineWin[j], pafThreeLineWin[nXSize + j],
469 : bSrcHasNoData, fSrcNoDataValue),
470 2184 : INTERPOL(pafThreeLineWin[jmax], pafThreeLineWin[nXSize + jmax],
471 : bSrcHasNoData, fSrcNoDataValue),
472 2184 : pafThreeLineWin[jmin],
473 2184 : pafThreeLineWin[j],
474 2184 : pafThreeLineWin[jmax],
475 2184 : pafThreeLineWin[nXSize + jmin],
476 2184 : pafThreeLineWin[nXSize + j],
477 2184 : pafThreeLineWin[nXSize + jmax]};
478 2184 : pafOutputBuf[j] = ComputeVal(
479 2184 : CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, bIsSrcNoDataNan,
480 2184 : afWin, fDstNoDataValue, pfnAlg, pData.get(), bComputeAtEdges);
481 : }
482 20 : eErr = GDALRasterIO(hDstBand, GF_Write, 0, 0, nXSize, 1, pafOutputBuf,
483 20 : nXSize, 1, GDT_Float32, 0, 0);
484 : }
485 : else
486 : {
487 : // Exclude the edges
488 5300 : for (int j = 0; j < nXSize; j++)
489 : {
490 5245 : pafOutputBuf[j] = fDstNoDataValue;
491 : }
492 55 : eErr = GDALRasterIO(hDstBand, GF_Write, 0, 0, nXSize, 1, pafOutputBuf,
493 : nXSize, 1, GDT_Float32, 0, 0);
494 :
495 55 : if (eErr == CE_None && nYSize > 1)
496 : {
497 55 : eErr = GDALRasterIO(hDstBand, GF_Write, 0, nYSize - 1, nXSize, 1,
498 : pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
499 : }
500 : }
501 75 : if (eErr != CE_None)
502 : {
503 0 : CPLFree(pafOutputBuf);
504 0 : CPLFree(pafThreeLineWin);
505 :
506 0 : return eErr;
507 : }
508 :
509 75 : int i = 1; // Used after for.
510 7613 : for (; i < nYSize - 1; i++)
511 : {
512 : /* Read third line of the line buffer */
513 : eErr =
514 15076 : GDALRasterIO(hSrcBand, GF_Read, 0, i + 1, nXSize, 1,
515 7538 : pafThreeLineWin + nLine3Off, nXSize, 1, eReadDT, 0, 0);
516 7538 : if (eErr != CE_None)
517 : {
518 0 : CPLFree(pafOutputBuf);
519 0 : CPLFree(pafThreeLineWin);
520 :
521 0 : return eErr;
522 : }
523 :
524 : // In case none of the 3 lines have nodata values, then no need to
525 : // check it in ComputeVal()
526 7538 : bool bOneOfThreeLinesHasNoData = CPL_TO_BOOL(bSrcHasNoData);
527 7538 : if (bSrcHasNoData)
528 : {
529 : if constexpr (std::numeric_limits<T>::is_integer)
530 : {
531 6078 : bool bLastLineHasNoDataValue = false;
532 6078 : int iX = 0;
533 180093 : for (; iX + 3 < nXSize; iX += 4)
534 : {
535 174249 : if (pafThreeLineWin[nLine3Off + iX] == fSrcNoDataValue ||
536 174015 : pafThreeLineWin[nLine3Off + iX + 1] ==
537 174015 : fSrcNoDataValue ||
538 174015 : pafThreeLineWin[nLine3Off + iX + 2] ==
539 174015 : fSrcNoDataValue ||
540 174015 : pafThreeLineWin[nLine3Off + iX + 3] == fSrcNoDataValue)
541 : {
542 234 : bLastLineHasNoDataValue = true;
543 234 : break;
544 : }
545 : }
546 6078 : if (!bLastLineHasNoDataValue)
547 : {
548 11941 : for (; iX < nXSize; iX++)
549 : {
550 6097 : if (pafThreeLineWin[nLine3Off + iX] == fSrcNoDataValue)
551 : {
552 234 : bLastLineHasNoDataValue = true;
553 : }
554 : }
555 : }
556 6078 : abLineHasNoDataValue[nLine3Off / nXSize] =
557 : bLastLineHasNoDataValue;
558 :
559 11689 : bOneOfThreeLinesHasNoData = abLineHasNoDataValue[0] ||
560 11689 : abLineHasNoDataValue[1] ||
561 : abLineHasNoDataValue[2];
562 : }
563 : else
564 : {
565 1264 : bool bLastLineHasNoDataValue = false;
566 1264 : int iX = 0;
567 30984 : for (; iX + 3 < nXSize; iX += 4)
568 : {
569 29720 : if (pafThreeLineWin[nLine3Off + iX] == fSrcNoDataValue ||
570 29720 : std::isnan(pafThreeLineWin[nLine3Off + iX]) ||
571 29720 : pafThreeLineWin[nLine3Off + iX + 1] ==
572 29720 : fSrcNoDataValue ||
573 29720 : std::isnan(pafThreeLineWin[nLine3Off + iX + 1]) ||
574 29720 : pafThreeLineWin[nLine3Off + iX + 2] ==
575 29720 : fSrcNoDataValue ||
576 29720 : std::isnan(pafThreeLineWin[nLine3Off + iX + 2]) ||
577 29720 : pafThreeLineWin[nLine3Off + iX + 3] ==
578 59656 : fSrcNoDataValue ||
579 29720 : std::isnan(pafThreeLineWin[nLine3Off + iX + 3]))
580 : {
581 216 : bLastLineHasNoDataValue = true;
582 216 : break;
583 : }
584 : }
585 1264 : if (!bLastLineHasNoDataValue)
586 : {
587 2432 : for (; iX < nXSize; iX++)
588 : {
589 2768 : if (pafThreeLineWin[nLine3Off + iX] ==
590 2458 : fSrcNoDataValue ||
591 1074 : std::isnan(pafThreeLineWin[nLine3Off + iX]))
592 : {
593 310 : bLastLineHasNoDataValue = true;
594 : }
595 : }
596 : }
597 1264 : abLineHasNoDataValue[nLine3Off / nXSize] =
598 : bLastLineHasNoDataValue;
599 :
600 2002 : bOneOfThreeLinesHasNoData = abLineHasNoDataValue[0] ||
601 2002 : abLineHasNoDataValue[1] ||
602 : abLineHasNoDataValue[2];
603 : }
604 : }
605 :
606 7538 : if (bComputeAtEdges && nXSize >= 2)
607 : {
608 2144 : int j = 0;
609 8576 : T afWin[9] = {INTERPOL(pafThreeLineWin[nLine1Off + j],
610 2144 : pafThreeLineWin[nLine1Off + j + 1],
611 : bSrcHasNoData, fSrcNoDataValue),
612 2144 : pafThreeLineWin[nLine1Off + j],
613 2144 : pafThreeLineWin[nLine1Off + j + 1],
614 4168 : INTERPOL(pafThreeLineWin[nLine2Off + j],
615 2144 : pafThreeLineWin[nLine2Off + j + 1],
616 : bSrcHasNoData, fSrcNoDataValue),
617 2144 : pafThreeLineWin[nLine2Off + j],
618 2144 : pafThreeLineWin[nLine2Off + j + 1],
619 4168 : INTERPOL(pafThreeLineWin[nLine3Off + j],
620 2144 : pafThreeLineWin[nLine3Off + j + 1],
621 : bSrcHasNoData, fSrcNoDataValue),
622 2144 : pafThreeLineWin[nLine3Off + j],
623 2144 : pafThreeLineWin[nLine3Off + j + 1]};
624 :
625 2144 : pafOutputBuf[j] = ComputeVal(
626 : bOneOfThreeLinesHasNoData, fSrcNoDataValue, bIsSrcNoDataNan,
627 4288 : afWin, fDstNoDataValue, pfnAlg, pData.get(), bComputeAtEdges);
628 : }
629 : else
630 : {
631 : // Exclude the edges
632 5394 : pafOutputBuf[0] = fDstNoDataValue;
633 : }
634 :
635 7538 : int j = 1;
636 7538 : if (pfnAlg_multisample && !bOneOfThreeLinesHasNoData)
637 : {
638 1183 : j = pfnAlg_multisample(
639 1183 : pafThreeLineWin + nLine1Off, pafThreeLineWin + nLine2Off,
640 1183 : pafThreeLineWin + nLine3Off, nXSize, pData.get(), pafOutputBuf);
641 : }
642 :
643 741011 : for (; j < nXSize - 1; j++)
644 : {
645 733473 : T afWin[9] = {pafThreeLineWin[nLine1Off + j - 1],
646 733473 : pafThreeLineWin[nLine1Off + j],
647 733473 : pafThreeLineWin[nLine1Off + j + 1],
648 733473 : pafThreeLineWin[nLine2Off + j - 1],
649 733473 : pafThreeLineWin[nLine2Off + j],
650 733473 : pafThreeLineWin[nLine2Off + j + 1],
651 733473 : pafThreeLineWin[nLine3Off + j - 1],
652 733473 : pafThreeLineWin[nLine3Off + j],
653 733473 : pafThreeLineWin[nLine3Off + j + 1]};
654 :
655 733473 : pafOutputBuf[j] = ComputeVal(
656 : bOneOfThreeLinesHasNoData, fSrcNoDataValue, bIsSrcNoDataNan,
657 733473 : afWin, fDstNoDataValue, pfnAlg, pData.get(), bComputeAtEdges);
658 : }
659 :
660 7538 : if (bComputeAtEdges && nXSize >= 2)
661 : {
662 2144 : j = nXSize - 1;
663 :
664 8576 : T afWin[9] = {pafThreeLineWin[nLine1Off + j - 1],
665 2144 : pafThreeLineWin[nLine1Off + j],
666 4168 : INTERPOL(pafThreeLineWin[nLine1Off + j],
667 2144 : pafThreeLineWin[nLine1Off + j - 1],
668 : bSrcHasNoData, fSrcNoDataValue),
669 2144 : pafThreeLineWin[nLine2Off + j - 1],
670 2144 : pafThreeLineWin[nLine2Off + j],
671 4168 : INTERPOL(pafThreeLineWin[nLine2Off + j],
672 2144 : pafThreeLineWin[nLine2Off + j - 1],
673 : bSrcHasNoData, fSrcNoDataValue),
674 2144 : pafThreeLineWin[nLine3Off + j - 1],
675 2144 : pafThreeLineWin[nLine3Off + j],
676 4168 : INTERPOL(pafThreeLineWin[nLine3Off + j],
677 2144 : pafThreeLineWin[nLine3Off + j - 1],
678 : bSrcHasNoData, fSrcNoDataValue)};
679 :
680 2144 : pafOutputBuf[j] = ComputeVal(
681 : bOneOfThreeLinesHasNoData, fSrcNoDataValue, bIsSrcNoDataNan,
682 4288 : afWin, fDstNoDataValue, pfnAlg, pData.get(), bComputeAtEdges);
683 : }
684 : else
685 : {
686 : // Exclude the edges
687 5394 : if (nXSize > 1)
688 5394 : pafOutputBuf[nXSize - 1] = fDstNoDataValue;
689 : }
690 :
691 : /* -----------------------------------------
692 : * Write Line to Raster
693 : */
694 7538 : eErr = GDALRasterIO(hDstBand, GF_Write, 0, i, nXSize, 1, pafOutputBuf,
695 : nXSize, 1, GDT_Float32, 0, 0);
696 7538 : if (eErr != CE_None)
697 : {
698 0 : CPLFree(pafOutputBuf);
699 0 : CPLFree(pafThreeLineWin);
700 :
701 0 : return eErr;
702 : }
703 :
704 7538 : if (!pfnProgress(1.0 * (i + 1) / nYSize, nullptr, pProgressData))
705 : {
706 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
707 0 : eErr = CE_Failure;
708 :
709 0 : CPLFree(pafOutputBuf);
710 0 : CPLFree(pafThreeLineWin);
711 :
712 0 : return eErr;
713 : }
714 :
715 7538 : const int nTemp = nLine1Off;
716 7538 : nLine1Off = nLine2Off;
717 7538 : nLine2Off = nLine3Off;
718 7538 : nLine3Off = nTemp;
719 : }
720 :
721 75 : if (bComputeAtEdges && nXSize >= 2 && nYSize >= 2)
722 : {
723 2204 : for (int j = 0; j < nXSize; j++)
724 : {
725 2184 : int jmin = (j == 0) ? j : j - 1;
726 2184 : int jmax = (j == nXSize - 1) ? j : j + 1;
727 :
728 8736 : T afWin[9] = {
729 2184 : pafThreeLineWin[nLine1Off + jmin],
730 2184 : pafThreeLineWin[nLine1Off + j],
731 2184 : pafThreeLineWin[nLine1Off + jmax],
732 2184 : pafThreeLineWin[nLine2Off + jmin],
733 2184 : pafThreeLineWin[nLine2Off + j],
734 2184 : pafThreeLineWin[nLine2Off + jmax],
735 4244 : INTERPOL(pafThreeLineWin[nLine2Off + jmin],
736 2184 : pafThreeLineWin[nLine1Off + jmin], bSrcHasNoData,
737 : fSrcNoDataValue),
738 4244 : INTERPOL(pafThreeLineWin[nLine2Off + j],
739 2184 : pafThreeLineWin[nLine1Off + j], bSrcHasNoData,
740 : fSrcNoDataValue),
741 4244 : INTERPOL(pafThreeLineWin[nLine2Off + jmax],
742 2184 : pafThreeLineWin[nLine1Off + jmax], bSrcHasNoData,
743 : fSrcNoDataValue),
744 : };
745 :
746 2184 : pafOutputBuf[j] = ComputeVal(
747 2184 : CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, bIsSrcNoDataNan,
748 2184 : afWin, fDstNoDataValue, pfnAlg, pData.get(), bComputeAtEdges);
749 : }
750 20 : eErr = GDALRasterIO(hDstBand, GF_Write, 0, i, nXSize, 1, pafOutputBuf,
751 : nXSize, 1, GDT_Float32, 0, 0);
752 20 : if (eErr != CE_None)
753 : {
754 0 : CPLFree(pafOutputBuf);
755 0 : CPLFree(pafThreeLineWin);
756 :
757 0 : return eErr;
758 : }
759 : }
760 :
761 75 : pfnProgress(1.0, nullptr, pProgressData);
762 75 : eErr = CE_None;
763 :
764 75 : CPLFree(pafOutputBuf);
765 75 : CPLFree(pafThreeLineWin);
766 :
767 75 : return eErr;
768 : }
769 :
770 : /************************************************************************/
771 : /* GradientAlg */
772 : /************************************************************************/
773 :
774 : template <class T, GradientAlg alg> struct Gradient
775 : {
776 : static void inline calc(const T *afWin, double inv_ewres, double inv_nsres,
777 : double &x, double &y);
778 : };
779 :
780 : template <class T> struct Gradient<T, GradientAlg::HORN>
781 : {
782 338425 : static void calc(const T *afWin, double inv_ewres, double inv_nsres,
783 : double &x, double &y)
784 : {
785 338425 : x = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
786 338425 : (afWin[2] + afWin[5] + afWin[5] + afWin[8])) *
787 : inv_ewres;
788 :
789 338425 : y = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
790 338425 : (afWin[0] + afWin[1] + afWin[1] + afWin[2])) *
791 : inv_nsres;
792 338425 : }
793 : };
794 :
795 : template <class T> struct Gradient<T, GradientAlg::ZEVENBERGEN_THORNE>
796 : {
797 142570 : static void calc(const T *afWin, double inv_ewres, double inv_nsres,
798 : double &x, double &y)
799 : {
800 142570 : x = (afWin[3] - afWin[5]) * inv_ewres;
801 142570 : y = (afWin[7] - afWin[1]) * inv_nsres;
802 142570 : }
803 : };
804 :
805 : /************************************************************************/
806 : /* GDALHillshade() */
807 : /************************************************************************/
808 :
809 : struct GDALHillshadeAlgData final : public AlgorithmParameters
810 : {
811 : double inv_nsres_yscale = 0;
812 : double inv_ewres_xscale = 0;
813 : double sin_altRadians = 0;
814 : double cos_alt_mul_z = 0;
815 : double azRadians = 0;
816 : double cos_az_mul_cos_alt_mul_z = 0;
817 : double sin_az_mul_cos_alt_mul_z = 0;
818 : double square_z = 0;
819 : double sin_altRadians_mul_254 = 0;
820 : double cos_az_mul_cos_alt_mul_z_mul_254 = 0;
821 : double sin_az_mul_cos_alt_mul_z_mul_254 = 0;
822 :
823 : double square_z_mul_square_inv_res = 0;
824 : double cos_az_mul_cos_alt_mul_z_mul_254_mul_inv_res = 0;
825 : double sin_az_mul_cos_alt_mul_z_mul_254_mul_inv_res = 0;
826 : double z_factor = 0;
827 :
828 : std::unique_ptr<AlgorithmParameters>
829 : CreateScaledParameters(double dfXRatio, double dfYRatio) override;
830 : };
831 :
832 : std::unique_ptr<AlgorithmParameters>
833 2 : GDALHillshadeAlgData::CreateScaledParameters(double dfXRatio, double dfYRatio)
834 : {
835 4 : auto newData = std::make_unique<GDALHillshadeAlgData>(*this);
836 2 : newData->inv_ewres_xscale /= dfXRatio;
837 2 : newData->inv_nsres_yscale /= dfYRatio;
838 :
839 2 : newData->square_z_mul_square_inv_res /= dfXRatio * dfXRatio;
840 2 : newData->cos_az_mul_cos_alt_mul_z_mul_254_mul_inv_res /= dfXRatio;
841 2 : newData->sin_az_mul_cos_alt_mul_z_mul_254_mul_inv_res /= dfXRatio;
842 :
843 4 : return newData;
844 : }
845 :
846 : /* Unoptimized formulas are :
847 : x = psData->z*((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
848 : (afWin[2] + afWin[5] + afWin[5] + afWin[8])) /
849 : (8.0 * psData->ewres * psData->xscale);
850 :
851 : y = psData->z*((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
852 : (afWin[0] + afWin[1] + afWin[1] + afWin[2])) /
853 : (8.0 * psData->nsres * psData->yscale);
854 :
855 : slope = atan(sqrt(x*x + y*y));
856 :
857 : aspect = atan2(y,x);
858 :
859 : cang = sin(alt) * cos(slope) +
860 : cos(alt) * sin(slope) *
861 : cos(az - M_PI/2 - aspect);
862 :
863 : We can avoid a lot of trigonometric computations:
864 :
865 : since cos(atan(x)) = 1 / sqrt(1+x^2)
866 : ==> cos(slope) = 1 / sqrt(1+ x*x+y*y)
867 :
868 : and sin(atan(x)) = x / sqrt(1+x^2)
869 : ==> sin(slope) = sqrt(x*x + y*y) / sqrt(1+ x*x+y*y)
870 :
871 : and cos(az - M_PI/2 - aspect)
872 : = cos(-az + M_PI/2 + aspect)
873 : = cos(M_PI/2 - (az - aspect))
874 : = sin(az - aspect)
875 : = -sin(aspect-az)
876 :
877 : ==> cang = (sin(alt) - cos(alt) * sqrt(x*x + y*y) * sin(aspect-az)) /
878 : sqrt(1+ x*x+y*y)
879 :
880 : But:
881 : sin(aspect - az) = sin(aspect)*cos(az) - cos(aspect)*sin(az))
882 :
883 : and as sin(aspect)=sin(atan2(y,x)) = y / sqrt(xx_plus_yy)
884 : and cos(aspect)=cos(atan2(y,x)) = x / sqrt(xx_plus_yy)
885 :
886 : sin(aspect - az) = (y * cos(az) - x * sin(az)) / sqrt(xx_plus_yy)
887 :
888 : so we get a final formula with just one transcendental function
889 : (reciprocal of square root):
890 :
891 : cang = (psData->sin_altRadians -
892 : (y * psData->cos_az_mul_cos_alt_mul_z -
893 : x * psData->sin_az_mul_cos_alt_mul_z)) /
894 : sqrt(1 + psData->square_z * xx_plus_yy);
895 : */
896 :
897 : #ifdef HAVE_SSE2
898 : inline double ApproxADivByInvSqrtB(double a, double b)
899 : {
900 : __m128d regB = _mm_load_sd(&b);
901 : __m128d regB_half = _mm_mul_sd(regB, _mm_set1_pd(0.5));
902 : // Compute rough approximation of 1 / sqrt(b) with _mm_rsqrt_ss
903 : regB =
904 : _mm_cvtss_sd(regB, _mm_rsqrt_ss(_mm_cvtsd_ss(_mm_setzero_ps(), regB)));
905 : // And perform one step of Newton-Raphson approximation to improve it
906 : // approx_inv_sqrt_x = approx_inv_sqrt_x*(1.5 -
907 : // 0.5*x*approx_inv_sqrt_x*approx_inv_sqrt_x);
908 : regB = _mm_mul_sd(
909 : regB, _mm_sub_sd(_mm_set1_pd(1.5),
910 : _mm_mul_sd(regB_half, _mm_mul_sd(regB, regB))));
911 : double dOut;
912 : _mm_store_sd(&dOut, regB);
913 : return a * dOut;
914 : }
915 : #else
916 565642 : inline double ApproxADivByInvSqrtB(double a, double b)
917 : {
918 565642 : return a / sqrt(b);
919 : }
920 : #endif
921 :
922 87846 : static double NormalizeAngle(double angle, double normalizer)
923 : {
924 87846 : angle = std::fmod(angle, normalizer);
925 87846 : if (angle < 0)
926 63029 : angle = normalizer + angle;
927 :
928 87846 : return angle;
929 : }
930 :
931 43923 : static double DifferenceBetweenAngles(double angle1, double angle2,
932 : double normalizer)
933 : {
934 : double diff =
935 43923 : NormalizeAngle(angle1, normalizer) - NormalizeAngle(angle2, normalizer);
936 43923 : diff = std::abs(diff);
937 43923 : if (diff > normalizer / 2)
938 10547 : diff = normalizer - diff;
939 43923 : return diff;
940 : }
941 :
942 : template <class T, GradientAlg alg>
943 43923 : static float GDALHillshadeIgorAlg(const T *afWin, float /*fDstNoDataValue*/,
944 : const AlgorithmParameters *pData)
945 : {
946 43923 : const GDALHillshadeAlgData *psData =
947 : static_cast<const GDALHillshadeAlgData *>(pData);
948 :
949 : double slopeDegrees;
950 : if (alg == GradientAlg::HORN)
951 : {
952 29282 : const double dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
953 29282 : (afWin[2] + afWin[5] + afWin[5] + afWin[8])) *
954 29282 : psData->inv_ewres_xscale;
955 :
956 29282 : const double dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
957 29282 : (afWin[0] + afWin[1] + afWin[1] + afWin[2])) *
958 29282 : psData->inv_nsres_yscale;
959 :
960 29282 : const double key = (dx * dx + dy * dy);
961 29282 : slopeDegrees = atan(sqrt(key) * psData->z_factor) * kdfRadiansToDegrees;
962 : }
963 : else // ZEVENBERGEN_THORNE
964 : {
965 14641 : const double dx = (afWin[3] - afWin[5]) * psData->inv_ewres_xscale;
966 14641 : const double dy = (afWin[7] - afWin[1]) * psData->inv_nsres_yscale;
967 14641 : const double key = dx * dx + dy * dy;
968 :
969 14641 : slopeDegrees = atan(sqrt(key) * psData->z_factor) * kdfRadiansToDegrees;
970 : }
971 :
972 : double aspect;
973 : if (alg == GradientAlg::HORN)
974 : {
975 29282 : const double dx = ((afWin[2] + afWin[5] + afWin[5] + afWin[8]) -
976 29282 : (afWin[0] + afWin[3] + afWin[3] + afWin[6]));
977 :
978 29282 : const double dy2 = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
979 29282 : (afWin[0] + afWin[1] + afWin[1] + afWin[2]));
980 :
981 29282 : aspect = atan2(dy2, -dx);
982 : }
983 : else // ZEVENBERGEN_THORNE
984 : {
985 14641 : const double dx = afWin[5] - afWin[3];
986 14641 : const double dy = afWin[7] - afWin[1];
987 14641 : aspect = atan2(dy, -dx);
988 : }
989 :
990 43923 : double slopeStrength = slopeDegrees / 90;
991 :
992 43923 : double aspectDiff = DifferenceBetweenAngles(
993 43923 : aspect, M_PI * 3 / 2 - psData->azRadians, M_PI * 2);
994 :
995 43923 : double aspectStrength = 1 - aspectDiff / M_PI;
996 :
997 43923 : double shadowness = 1.0 - slopeStrength * aspectStrength;
998 :
999 43923 : return static_cast<float>(255.0 * shadowness);
1000 : }
1001 :
1002 : template <class T, GradientAlg alg>
1003 294982 : static float GDALHillshadeAlg(const T *afWin, float /*fDstNoDataValue*/,
1004 : const AlgorithmParameters *pData)
1005 : {
1006 294982 : const GDALHillshadeAlgData *psData =
1007 : static_cast<const GDALHillshadeAlgData *>(pData);
1008 :
1009 : // First Slope ...
1010 : double x, y;
1011 294982 : Gradient<T, alg>::calc(afWin, psData->inv_ewres_xscale,
1012 294982 : psData->inv_nsres_yscale, x, y);
1013 :
1014 294982 : const double xx_plus_yy = x * x + y * y;
1015 :
1016 : // ... then the shade value
1017 : const double cang_mul_254 =
1018 294982 : ApproxADivByInvSqrtB(psData->sin_altRadians_mul_254 -
1019 294982 : (y * psData->cos_az_mul_cos_alt_mul_z_mul_254 -
1020 294982 : x * psData->sin_az_mul_cos_alt_mul_z_mul_254),
1021 294982 : 1 + psData->square_z * xx_plus_yy);
1022 :
1023 294982 : const double cang = cang_mul_254 <= 0.0 ? 1.0 : 1.0 + cang_mul_254;
1024 :
1025 294982 : return static_cast<float>(cang);
1026 : }
1027 :
1028 : template <class T>
1029 97511 : static float GDALHillshadeAlg_same_res(const T *afWin,
1030 : float /*fDstNoDataValue*/,
1031 : const AlgorithmParameters *pData)
1032 : {
1033 97511 : const GDALHillshadeAlgData *psData =
1034 : static_cast<const GDALHillshadeAlgData *>(pData);
1035 :
1036 : // First Slope ...
1037 : /*x = (afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
1038 : (afWin[2] + afWin[5] + afWin[5] + afWin[8]);
1039 :
1040 : y = (afWin[0] + afWin[1] + afWin[1] + afWin[2]) -
1041 : (afWin[6] + afWin[7] + afWin[7] + afWin[8]);*/
1042 :
1043 97511 : T accX = afWin[0] - afWin[8];
1044 97511 : const T six_minus_two = afWin[6] - afWin[2];
1045 97511 : T accY = accX;
1046 97511 : const T three_minus_five = afWin[3] - afWin[5];
1047 97511 : const T one_minus_seven = afWin[1] - afWin[7];
1048 97511 : accX += three_minus_five;
1049 97511 : accY += one_minus_seven;
1050 97511 : accX += three_minus_five;
1051 97511 : accY += one_minus_seven;
1052 97511 : accX += six_minus_two;
1053 97511 : accY -= six_minus_two;
1054 97511 : const double x = accX;
1055 97511 : const double y = accY;
1056 :
1057 97511 : const double xx_plus_yy = x * x + y * y;
1058 :
1059 : // ... then the shade value
1060 97511 : const double cang_mul_254 = ApproxADivByInvSqrtB(
1061 97511 : psData->sin_altRadians_mul_254 +
1062 97511 : (x * psData->sin_az_mul_cos_alt_mul_z_mul_254_mul_inv_res +
1063 97511 : y * psData->cos_az_mul_cos_alt_mul_z_mul_254_mul_inv_res),
1064 97511 : 1 + psData->square_z_mul_square_inv_res * xx_plus_yy);
1065 :
1066 97511 : const double cang = cang_mul_254 <= 0.0 ? 1.0 : 1.0 + cang_mul_254;
1067 :
1068 97511 : return static_cast<float>(cang);
1069 : }
1070 :
1071 : #ifdef HAVE_16_SSE_REG
1072 :
1073 : template <class T, class REG_T>
1074 1659 : static int GDALHillshadeAlg_same_res_multisample(
1075 : const T *pafFirstLine, const T *pafSecondLine, const T *pafThirdLine,
1076 : int nXSize, const AlgorithmParameters *pData, float *pafOutputBuf)
1077 : {
1078 : // Only valid for T == int
1079 1659 : const GDALHillshadeAlgData *psData =
1080 : static_cast<const GDALHillshadeAlgData *>(pData);
1081 1659 : const auto reg_fact_x = XMMReg4Double::Set1(
1082 1659 : psData->sin_az_mul_cos_alt_mul_z_mul_254_mul_inv_res);
1083 1659 : const auto reg_fact_y = XMMReg4Double::Set1(
1084 1659 : psData->cos_az_mul_cos_alt_mul_z_mul_254_mul_inv_res);
1085 1659 : const auto reg_constant_num =
1086 1659 : XMMReg4Double::Set1(psData->sin_altRadians_mul_254);
1087 1659 : const auto reg_constant_denom =
1088 1659 : XMMReg4Double::Set1(psData->square_z_mul_square_inv_res);
1089 1659 : const auto reg_half = XMMReg4Double::Set1(0.5);
1090 1659 : const auto reg_one = reg_half + reg_half;
1091 1659 : const auto reg_one_float = XMMReg4Float::Set1(1.0f);
1092 :
1093 1659 : int j = 1; // Used after for.
1094 48650 : for (; j < nXSize - 4; j += 4)
1095 : {
1096 46991 : const T *firstLine = pafFirstLine + j - 1;
1097 46991 : const T *secondLine = pafSecondLine + j - 1;
1098 46991 : const T *thirdLine = pafThirdLine + j - 1;
1099 :
1100 46991 : const auto firstLine0 = REG_T::Load4Val(firstLine);
1101 46991 : const auto firstLine1 = REG_T::Load4Val(firstLine + 1);
1102 46991 : const auto firstLine2 = REG_T::Load4Val(firstLine + 2);
1103 46991 : const auto thirdLine0 = REG_T::Load4Val(thirdLine);
1104 46991 : const auto thirdLine1 = REG_T::Load4Val(thirdLine + 1);
1105 46991 : const auto thirdLine2 = REG_T::Load4Val(thirdLine + 2);
1106 46991 : auto accX = firstLine0 - thirdLine2;
1107 46991 : const auto six_minus_two = thirdLine0 - firstLine2;
1108 46991 : auto accY = accX;
1109 46991 : const auto three_minus_five =
1110 : REG_T::Load4Val(secondLine) - REG_T::Load4Val(secondLine + 2);
1111 46991 : const auto one_minus_seven = firstLine1 - thirdLine1;
1112 46991 : accX += three_minus_five;
1113 46991 : accY += one_minus_seven;
1114 46991 : accX += three_minus_five;
1115 46991 : accY += one_minus_seven;
1116 46991 : accX += six_minus_two;
1117 46991 : accY -= six_minus_two;
1118 :
1119 46991 : const auto reg_x = accX.cast_to_double();
1120 46991 : const auto reg_y = accY.cast_to_double();
1121 46991 : const auto reg_xx_plus_yy = reg_x * reg_x + reg_y * reg_y;
1122 46991 : const auto reg_numerator =
1123 : reg_constant_num + reg_fact_x * reg_x + reg_fact_y * reg_y;
1124 46991 : const auto reg_denominator =
1125 : reg_one + reg_constant_denom * reg_xx_plus_yy;
1126 46991 : const auto num_div_sqrt_denom =
1127 : reg_numerator * reg_denominator.approx_inv_sqrt(reg_one, reg_half);
1128 :
1129 46991 : auto res = num_div_sqrt_denom.cast_to_float();
1130 46991 : res = XMMReg4Float::Max(reg_one_float, res + reg_one_float);
1131 46991 : res.Store4Val(pafOutputBuf + j);
1132 : }
1133 1659 : return j;
1134 : }
1135 : #endif
1136 :
1137 : static const double INV_SQUARE_OF_HALF_PI = 1.0 / ((M_PI * M_PI) / 4);
1138 :
1139 : template <class T, GradientAlg alg>
1140 142090 : static float GDALHillshadeCombinedAlg(const T *afWin, float /*fDstNoDataValue*/,
1141 : const AlgorithmParameters *pData)
1142 : {
1143 142090 : const GDALHillshadeAlgData *psData =
1144 : static_cast<const GDALHillshadeAlgData *>(pData);
1145 :
1146 : // First Slope ...
1147 : double x, y;
1148 142090 : Gradient<T, alg>::calc(afWin, psData->inv_ewres_xscale,
1149 142090 : psData->inv_nsres_yscale, x, y);
1150 :
1151 142090 : const double xx_plus_yy = x * x + y * y;
1152 :
1153 142090 : const double slope = xx_plus_yy * psData->square_z;
1154 :
1155 : // ... then the shade value
1156 284180 : double cang = acos(ApproxADivByInvSqrtB(
1157 142090 : psData->sin_altRadians - (y * psData->cos_az_mul_cos_alt_mul_z -
1158 142090 : x * psData->sin_az_mul_cos_alt_mul_z),
1159 : 1 + slope));
1160 :
1161 : // combined shading
1162 142090 : cang = 1 - cang * atan(sqrt(slope)) * INV_SQUARE_OF_HALF_PI;
1163 :
1164 142090 : const float fcang =
1165 142086 : cang <= 0.0 ? 1.0f : static_cast<float>(1.0 + (254.0 * cang));
1166 :
1167 142090 : return fcang;
1168 : }
1169 :
1170 : static std::unique_ptr<AlgorithmParameters>
1171 69 : GDALCreateHillshadeData(const double *adfGeoTransform, double z, double xscale,
1172 : double yscale, double alt, double az, GradientAlg eAlg)
1173 : {
1174 138 : auto pData = std::make_unique<GDALHillshadeAlgData>();
1175 :
1176 69 : pData->inv_nsres_yscale = 1.0 / (adfGeoTransform[5] * yscale);
1177 69 : pData->inv_ewres_xscale = 1.0 / (adfGeoTransform[1] * xscale);
1178 69 : pData->sin_altRadians = sin(alt * kdfDegreesToRadians);
1179 69 : pData->azRadians = az * kdfDegreesToRadians;
1180 69 : pData->z_factor = z / (eAlg == GradientAlg::ZEVENBERGEN_THORNE ? 2 : 8);
1181 69 : pData->cos_alt_mul_z = cos(alt * kdfDegreesToRadians) * pData->z_factor;
1182 138 : pData->cos_az_mul_cos_alt_mul_z =
1183 69 : cos(pData->azRadians) * pData->cos_alt_mul_z;
1184 138 : pData->sin_az_mul_cos_alt_mul_z =
1185 69 : sin(pData->azRadians) * pData->cos_alt_mul_z;
1186 69 : pData->square_z = pData->z_factor * pData->z_factor;
1187 :
1188 69 : pData->sin_altRadians_mul_254 = 254.0 * pData->sin_altRadians;
1189 138 : pData->cos_az_mul_cos_alt_mul_z_mul_254 =
1190 69 : 254.0 * pData->cos_az_mul_cos_alt_mul_z;
1191 138 : pData->sin_az_mul_cos_alt_mul_z_mul_254 =
1192 69 : 254.0 * pData->sin_az_mul_cos_alt_mul_z;
1193 :
1194 69 : if (adfGeoTransform[1] == -adfGeoTransform[5] && xscale == yscale)
1195 : {
1196 86 : pData->square_z_mul_square_inv_res =
1197 43 : pData->square_z * pData->inv_ewres_xscale * pData->inv_ewres_xscale;
1198 86 : pData->cos_az_mul_cos_alt_mul_z_mul_254_mul_inv_res =
1199 43 : pData->cos_az_mul_cos_alt_mul_z_mul_254 * -pData->inv_ewres_xscale;
1200 43 : pData->sin_az_mul_cos_alt_mul_z_mul_254_mul_inv_res =
1201 43 : pData->sin_az_mul_cos_alt_mul_z_mul_254 * pData->inv_ewres_xscale;
1202 : }
1203 :
1204 138 : return pData;
1205 : }
1206 :
1207 : /************************************************************************/
1208 : /* GDALHillshadeMultiDirectional() */
1209 : /************************************************************************/
1210 :
1211 : struct GDALHillshadeMultiDirectionalAlgData final : public AlgorithmParameters
1212 : {
1213 : double inv_nsres_yscale = 0;
1214 : double inv_ewres_xscale = 0;
1215 : double square_z = 0;
1216 : double sin_altRadians_mul_127 = 0;
1217 : double sin_altRadians_mul_254 = 0;
1218 :
1219 : double cos_alt_mul_z_mul_127 = 0;
1220 : double cos225_az_mul_cos_alt_mul_z_mul_127 = 0;
1221 :
1222 : std::unique_ptr<AlgorithmParameters>
1223 : CreateScaledParameters(double dfXRatio, double dfYRatio) override;
1224 : };
1225 :
1226 : std::unique_ptr<AlgorithmParameters>
1227 0 : GDALHillshadeMultiDirectionalAlgData::CreateScaledParameters(double dfXRatio,
1228 : double dfYRatio)
1229 : {
1230 : auto newData =
1231 0 : std::make_unique<GDALHillshadeMultiDirectionalAlgData>(*this);
1232 0 : newData->inv_ewres_xscale /= dfXRatio;
1233 0 : newData->inv_nsres_yscale /= dfYRatio;
1234 0 : return newData;
1235 : }
1236 :
1237 : template <class T, GradientAlg alg>
1238 43923 : static float GDALHillshadeMultiDirectionalAlg(const T *afWin,
1239 : float /*fDstNoDataValue*/,
1240 : const AlgorithmParameters *pData)
1241 : {
1242 43923 : const GDALHillshadeMultiDirectionalAlgData *psData =
1243 : static_cast<const GDALHillshadeMultiDirectionalAlgData *>(pData);
1244 :
1245 : // First Slope ...
1246 : double x, y;
1247 43923 : Gradient<T, alg>::calc(afWin, psData->inv_ewres_xscale,
1248 43923 : psData->inv_nsres_yscale, x, y);
1249 :
1250 : // See http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
1251 : // W225 = sin^2(aspect - 225) = 0.5 * (1 - 2 * sin(aspect) * cos(aspect))
1252 : // W270 = sin^2(aspect - 270) = cos^2(aspect)
1253 : // W315 = sin^2(aspect - 315) = 0.5 * (1 + 2 * sin(aspect) * cos(aspect))
1254 : // W360 = sin^2(aspect - 360) = sin^2(aspect)
1255 : // hillshade= 0.5 * (W225 * hillshade(az=225) +
1256 : // W270 * hillshade(az=270) +
1257 : // W315 * hillshade(az=315) +
1258 : // W360 * hillshade(az=360))
1259 :
1260 43923 : const double xx = x * x;
1261 43923 : const double yy = y * y;
1262 43923 : const double xx_plus_yy = xx + yy;
1263 43923 : if (xx_plus_yy == 0.0)
1264 12864 : return static_cast<float>(1.0 + psData->sin_altRadians_mul_254);
1265 :
1266 : // ... then the shade value from different azimuth
1267 31059 : double val225_mul_127 =
1268 31059 : psData->sin_altRadians_mul_127 +
1269 31059 : (x - y) * psData->cos225_az_mul_cos_alt_mul_z_mul_127;
1270 31059 : val225_mul_127 = (val225_mul_127 <= 0.0) ? 0.0 : val225_mul_127;
1271 31059 : double val270_mul_127 =
1272 31059 : psData->sin_altRadians_mul_127 - x * psData->cos_alt_mul_z_mul_127;
1273 31059 : val270_mul_127 = (val270_mul_127 <= 0.0) ? 0.0 : val270_mul_127;
1274 31059 : double val315_mul_127 =
1275 31059 : psData->sin_altRadians_mul_127 +
1276 31059 : (x + y) * psData->cos225_az_mul_cos_alt_mul_z_mul_127;
1277 31059 : val315_mul_127 = (val315_mul_127 <= 0.0) ? 0.0 : val315_mul_127;
1278 31059 : double val360_mul_127 =
1279 31059 : psData->sin_altRadians_mul_127 - y * psData->cos_alt_mul_z_mul_127;
1280 31059 : val360_mul_127 = (val360_mul_127 <= 0.0) ? 0.0 : val360_mul_127;
1281 :
1282 : // ... then the weighted shading
1283 31059 : const double weight_225 = 0.5 * xx_plus_yy - x * y;
1284 31059 : const double weight_270 = xx;
1285 31059 : const double weight_315 = xx_plus_yy - weight_225;
1286 31059 : const double weight_360 = yy;
1287 31059 : const double cang_mul_127 = ApproxADivByInvSqrtB(
1288 31059 : (weight_225 * val225_mul_127 + weight_270 * val270_mul_127 +
1289 31059 : weight_315 * val315_mul_127 + weight_360 * val360_mul_127) /
1290 : xx_plus_yy,
1291 31059 : 1 + psData->square_z * xx_plus_yy);
1292 :
1293 31059 : const double cang = 1.0 + cang_mul_127;
1294 :
1295 31059 : return static_cast<float>(cang);
1296 : }
1297 :
1298 : static std::unique_ptr<AlgorithmParameters>
1299 3 : GDALCreateHillshadeMultiDirectionalData(const double *adfGeoTransform, double z,
1300 : double xscale, double yscale,
1301 : double alt, GradientAlg eAlg)
1302 : {
1303 6 : auto pData = std::make_unique<GDALHillshadeMultiDirectionalAlgData>();
1304 :
1305 3 : pData->inv_nsres_yscale = 1.0 / (adfGeoTransform[5] * yscale);
1306 3 : pData->inv_ewres_xscale = 1.0 / (adfGeoTransform[1] * xscale);
1307 3 : const double z_factor =
1308 3 : z / (eAlg == GradientAlg::ZEVENBERGEN_THORNE ? 2 : 8);
1309 3 : const double cos_alt_mul_z = cos(alt * kdfDegreesToRadians) * z_factor;
1310 3 : pData->square_z = z_factor * z_factor;
1311 :
1312 3 : pData->sin_altRadians_mul_127 = 127.0 * sin(alt * kdfDegreesToRadians);
1313 3 : pData->sin_altRadians_mul_254 = 254.0 * sin(alt * kdfDegreesToRadians);
1314 3 : pData->cos_alt_mul_z_mul_127 = 127.0 * cos_alt_mul_z;
1315 6 : pData->cos225_az_mul_cos_alt_mul_z_mul_127 =
1316 3 : 127.0 * cos(225 * kdfDegreesToRadians) * cos_alt_mul_z;
1317 :
1318 6 : return pData;
1319 : }
1320 :
1321 : /************************************************************************/
1322 : /* GDALSlope() */
1323 : /************************************************************************/
1324 :
1325 : struct GDALSlopeAlgData final : public AlgorithmParameters
1326 : {
1327 : double nsres_yscale = 0;
1328 : double ewres_xscale = 0;
1329 : int slopeFormat = 0;
1330 :
1331 : std::unique_ptr<AlgorithmParameters>
1332 : CreateScaledParameters(double dfXRatio, double dfYRatio) override;
1333 : };
1334 :
1335 : std::unique_ptr<AlgorithmParameters>
1336 2 : GDALSlopeAlgData::CreateScaledParameters(double dfXRatio, double dfYRatio)
1337 : {
1338 4 : auto newData = std::make_unique<GDALSlopeAlgData>(*this);
1339 2 : newData->ewres_xscale *= dfXRatio;
1340 2 : newData->nsres_yscale *= dfYRatio;
1341 4 : return newData;
1342 : }
1343 :
1344 : template <class T>
1345 123610 : static float GDALSlopeHornAlg(const T *afWin, float /*fDstNoDataValue*/,
1346 : const AlgorithmParameters *pData)
1347 : {
1348 123610 : const GDALSlopeAlgData *psData =
1349 : static_cast<const GDALSlopeAlgData *>(pData);
1350 :
1351 123610 : const double dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
1352 123610 : (afWin[2] + afWin[5] + afWin[5] + afWin[8])) /
1353 123610 : psData->ewres_xscale;
1354 :
1355 123610 : const double dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
1356 123610 : (afWin[0] + afWin[1] + afWin[1] + afWin[2])) /
1357 123610 : psData->nsres_yscale;
1358 :
1359 123610 : const double key = (dx * dx + dy * dy);
1360 :
1361 123610 : if (psData->slopeFormat == 1)
1362 108969 : return static_cast<float>(atan(sqrt(key) / 8) * kdfRadiansToDegrees);
1363 :
1364 14641 : return static_cast<float>(100 * (sqrt(key) / 8));
1365 : }
1366 :
1367 : template <class T>
1368 85446 : static float GDALSlopeZevenbergenThorneAlg(const T *afWin,
1369 : float /*fDstNoDataValue*/,
1370 : const AlgorithmParameters *pData)
1371 : {
1372 85446 : const GDALSlopeAlgData *psData =
1373 : static_cast<const GDALSlopeAlgData *>(pData);
1374 :
1375 85446 : const double dx = (afWin[3] - afWin[5]) / psData->ewres_xscale;
1376 85446 : const double dy = (afWin[7] - afWin[1]) / psData->nsres_yscale;
1377 85446 : const double key = dx * dx + dy * dy;
1378 :
1379 85446 : if (psData->slopeFormat == 1)
1380 85446 : return static_cast<float>(atan(sqrt(key) / 2) * kdfRadiansToDegrees);
1381 :
1382 0 : return static_cast<float>(100 * (sqrt(key) / 2));
1383 : }
1384 :
1385 : static std::unique_ptr<AlgorithmParameters>
1386 16 : GDALCreateSlopeData(double *adfGeoTransform, double xscale, double yscale,
1387 : int slopeFormat)
1388 : {
1389 32 : auto pData = std::make_unique<GDALSlopeAlgData>();
1390 16 : pData->nsres_yscale = adfGeoTransform[5] * yscale;
1391 16 : pData->ewres_xscale = adfGeoTransform[1] * xscale;
1392 16 : pData->slopeFormat = slopeFormat;
1393 32 : return pData;
1394 : }
1395 :
1396 : /************************************************************************/
1397 : /* GDALAspect() */
1398 : /************************************************************************/
1399 :
1400 : struct GDALAspectAlgData final : public AlgorithmParameters
1401 : {
1402 : bool bAngleAsAzimuth = false;
1403 :
1404 : std::unique_ptr<AlgorithmParameters>
1405 : CreateScaledParameters(double, double) override;
1406 : };
1407 :
1408 : std::unique_ptr<AlgorithmParameters>
1409 2 : GDALAspectAlgData::CreateScaledParameters(double, double)
1410 : {
1411 2 : return std::make_unique<GDALAspectAlgData>(*this);
1412 : }
1413 :
1414 : template <class T>
1415 108969 : static float GDALAspectAlg(const T *afWin, float fDstNoDataValue,
1416 : const AlgorithmParameters *pData)
1417 : {
1418 108969 : const GDALAspectAlgData *psData =
1419 : static_cast<const GDALAspectAlgData *>(pData);
1420 :
1421 108969 : const double dx = ((afWin[2] + afWin[5] + afWin[5] + afWin[8]) -
1422 108969 : (afWin[0] + afWin[3] + afWin[3] + afWin[6]));
1423 :
1424 108969 : const double dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
1425 108969 : (afWin[0] + afWin[1] + afWin[1] + afWin[2]));
1426 :
1427 108969 : float aspect = static_cast<float>(atan2(dy, -dx) / kdfDegreesToRadians);
1428 :
1429 108969 : if (dx == 0 && dy == 0)
1430 : {
1431 : /* Flat area */
1432 31487 : aspect = fDstNoDataValue;
1433 : }
1434 77482 : else if (psData->bAngleAsAzimuth)
1435 : {
1436 67088 : if (aspect > 90.0f)
1437 8356 : aspect = 450.0f - aspect;
1438 : else
1439 58732 : aspect = 90.0f - aspect;
1440 : }
1441 : else
1442 : {
1443 10394 : if (aspect < 0)
1444 6464 : aspect += 360.0f;
1445 : }
1446 :
1447 108969 : if (aspect == 360.0f)
1448 0 : aspect = 0.0;
1449 :
1450 108969 : return aspect;
1451 : }
1452 :
1453 : template <class T>
1454 42963 : static float GDALAspectZevenbergenThorneAlg(const T *afWin,
1455 : float fDstNoDataValue,
1456 : const AlgorithmParameters *pData)
1457 : {
1458 42963 : const GDALAspectAlgData *psData =
1459 : static_cast<const GDALAspectAlgData *>(pData);
1460 :
1461 42963 : const double dx = afWin[5] - afWin[3];
1462 42963 : const double dy = afWin[7] - afWin[1];
1463 42963 : float aspect = static_cast<float>(atan2(dy, -dx) / kdfDegreesToRadians);
1464 42963 : if (dx == 0 && dy == 0)
1465 : {
1466 : /* Flat area */
1467 12974 : aspect = fDstNoDataValue;
1468 : }
1469 29989 : else if (psData->bAngleAsAzimuth)
1470 : {
1471 29989 : if (aspect > 90.0f)
1472 3780 : aspect = 450.0f - aspect;
1473 : else
1474 26209 : aspect = 90.0f - aspect;
1475 : }
1476 : else
1477 : {
1478 0 : if (aspect < 0)
1479 0 : aspect += 360.0f;
1480 : }
1481 :
1482 42963 : if (aspect == 360.0f)
1483 0 : aspect = 0.0;
1484 :
1485 42963 : return aspect;
1486 : }
1487 :
1488 : static std::unique_ptr<AlgorithmParameters>
1489 12 : GDALCreateAspectData(bool bAngleAsAzimuth)
1490 : {
1491 24 : auto pData = std::make_unique<GDALAspectAlgData>();
1492 12 : pData->bAngleAsAzimuth = bAngleAsAzimuth;
1493 24 : return pData;
1494 : }
1495 :
1496 : /************************************************************************/
1497 : /* GDALColorRelief() */
1498 : /************************************************************************/
1499 :
1500 410 : static int GDALColorReliefSortColors(const GDALColorAssociation &pA,
1501 : const GDALColorAssociation &pB)
1502 : {
1503 : /* Sort NaN in first position */
1504 819 : return (std::isnan(pA.dfVal) && !std::isnan(pB.dfVal)) ||
1505 819 : pA.dfVal < pB.dfVal;
1506 : }
1507 :
1508 39 : static void GDALColorReliefProcessColors(
1509 : std::vector<GDALColorAssociation> &asColorAssociation, int bSrcHasNoData,
1510 : double dfSrcNoDataValue, ColorSelectionMode eColorSelectionMode)
1511 : {
1512 39 : std::stable_sort(asColorAssociation.begin(), asColorAssociation.end(),
1513 : GDALColorReliefSortColors);
1514 :
1515 39 : size_t nRepeatedEntryIndex = 0;
1516 39 : const size_t nInitialSize = asColorAssociation.size();
1517 222 : for (size_t i = 1; i < nInitialSize; ++i)
1518 : {
1519 183 : const GDALColorAssociation *pPrevious = &asColorAssociation[i - 1];
1520 183 : const GDALColorAssociation *pCurrent = &asColorAssociation[i];
1521 :
1522 : // NaN comparison is always false, so it handles itself
1523 183 : if (eColorSelectionMode != COLOR_SELECTION_EXACT_ENTRY &&
1524 153 : bSrcHasNoData && pCurrent->dfVal == dfSrcNoDataValue)
1525 : {
1526 : // Check if there is enough distance between the nodata value and
1527 : // its predecessor.
1528 6 : const double dfNewValue = std::nextafter(
1529 6 : pCurrent->dfVal, -std::numeric_limits<double>::infinity());
1530 6 : if (dfNewValue > pPrevious->dfVal)
1531 : {
1532 : // add one just below the nodata value
1533 6 : GDALColorAssociation sNew = *pPrevious;
1534 6 : sNew.dfVal = dfNewValue;
1535 6 : asColorAssociation.push_back(std::move(sNew));
1536 6 : }
1537 : }
1538 177 : else if (eColorSelectionMode != COLOR_SELECTION_EXACT_ENTRY &&
1539 147 : bSrcHasNoData && pPrevious->dfVal == dfSrcNoDataValue)
1540 : {
1541 : // Check if there is enough distance between the nodata value and
1542 : // its successor.
1543 7 : const double dfNewValue = std::nextafter(
1544 7 : pPrevious->dfVal, std::numeric_limits<double>::infinity());
1545 7 : if (dfNewValue < pCurrent->dfVal)
1546 : {
1547 : // add one just above the nodata value
1548 7 : GDALColorAssociation sNew = *pCurrent;
1549 7 : sNew.dfVal = dfNewValue;
1550 7 : asColorAssociation.push_back(std::move(sNew));
1551 7 : }
1552 : }
1553 170 : else if (nRepeatedEntryIndex == 0 &&
1554 168 : pCurrent->dfVal == pPrevious->dfVal)
1555 : {
1556 : // second of a series of equivalent entries
1557 2 : nRepeatedEntryIndex = i;
1558 : }
1559 168 : else if (nRepeatedEntryIndex != 0 &&
1560 2 : pCurrent->dfVal != pPrevious->dfVal)
1561 : {
1562 : // Get the distance between the predecessor and successor of the
1563 : // equivalent entries.
1564 2 : double dfTotalDist = 0.0;
1565 2 : double dfLeftDist = 0.0;
1566 2 : if (nRepeatedEntryIndex >= 2)
1567 : {
1568 : const GDALColorAssociation *pLower =
1569 2 : &asColorAssociation[nRepeatedEntryIndex - 2];
1570 2 : dfTotalDist = pCurrent->dfVal - pLower->dfVal;
1571 2 : dfLeftDist = pPrevious->dfVal - pLower->dfVal;
1572 : }
1573 : else
1574 : {
1575 0 : dfTotalDist = pCurrent->dfVal - pPrevious->dfVal;
1576 : }
1577 :
1578 : // check if this distance is enough
1579 2 : const size_t nEquivalentCount = i - nRepeatedEntryIndex + 1;
1580 2 : if (dfTotalDist >
1581 2 : std::abs(pPrevious->dfVal) * nEquivalentCount * DBL_EPSILON)
1582 : {
1583 : // balance the alterations
1584 2 : double dfMultiplier =
1585 2 : 0.5 - double(nEquivalentCount) * dfLeftDist / dfTotalDist;
1586 6 : for (auto j = nRepeatedEntryIndex - 1; j < i; ++j)
1587 : {
1588 8 : asColorAssociation[j].dfVal +=
1589 4 : (std::abs(pPrevious->dfVal) * dfMultiplier) *
1590 : DBL_EPSILON;
1591 4 : dfMultiplier += 1.0;
1592 : }
1593 : }
1594 : else
1595 : {
1596 : // Fallback to the old behavior: keep equivalent entries as
1597 : // they are.
1598 : }
1599 :
1600 2 : nRepeatedEntryIndex = 0;
1601 : }
1602 : }
1603 :
1604 39 : if (nInitialSize != asColorAssociation.size())
1605 : {
1606 11 : std::stable_sort(asColorAssociation.begin(), asColorAssociation.end(),
1607 : GDALColorReliefSortColors);
1608 : }
1609 39 : }
1610 :
1611 294622 : static bool GDALColorReliefGetRGBA(
1612 : const std::vector<GDALColorAssociation> &asColorAssociation, double dfVal,
1613 : ColorSelectionMode eColorSelectionMode, int *pnR, int *pnG, int *pnB,
1614 : int *pnA)
1615 : {
1616 294622 : CPLAssert(!asColorAssociation.empty());
1617 :
1618 294622 : size_t lower = 0;
1619 :
1620 : // Special case for NaN
1621 294622 : if (std::isnan(asColorAssociation[0].dfVal))
1622 : {
1623 4 : if (std::isnan(dfVal))
1624 : {
1625 1 : *pnR = asColorAssociation[0].nR;
1626 1 : *pnG = asColorAssociation[0].nG;
1627 1 : *pnB = asColorAssociation[0].nB;
1628 1 : *pnA = asColorAssociation[0].nA;
1629 1 : return true;
1630 : }
1631 : else
1632 : {
1633 3 : lower = 1;
1634 : }
1635 : }
1636 :
1637 : // Find the index of the first element in the LUT input array that
1638 : // is not smaller than the dfVal value.
1639 294621 : size_t i = 0;
1640 294621 : size_t upper = asColorAssociation.size() - 1;
1641 : while (true)
1642 : {
1643 951279 : const size_t mid = (lower + upper) / 2;
1644 951279 : if (upper - lower <= 1)
1645 : {
1646 294621 : if (dfVal <= asColorAssociation[lower].dfVal)
1647 12 : i = lower;
1648 294609 : else if (dfVal <= asColorAssociation[upper].dfVal)
1649 293599 : i = upper;
1650 : else
1651 1010 : i = upper + 1;
1652 294621 : break;
1653 : }
1654 656658 : else if (asColorAssociation[mid].dfVal >= dfVal)
1655 : {
1656 387778 : upper = mid;
1657 : }
1658 : else
1659 : {
1660 268880 : lower = mid;
1661 : }
1662 656658 : }
1663 :
1664 294621 : if (i == 0)
1665 : {
1666 11 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
1667 2 : asColorAssociation[0].dfVal != dfVal)
1668 : {
1669 0 : *pnR = 0;
1670 0 : *pnG = 0;
1671 0 : *pnB = 0;
1672 0 : *pnA = 0;
1673 0 : return false;
1674 : }
1675 : else
1676 : {
1677 9 : *pnR = asColorAssociation[0].nR;
1678 9 : *pnG = asColorAssociation[0].nG;
1679 9 : *pnB = asColorAssociation[0].nB;
1680 9 : *pnA = asColorAssociation[0].nA;
1681 9 : return true;
1682 : }
1683 : }
1684 294612 : else if (i == asColorAssociation.size())
1685 : {
1686 1262 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
1687 252 : asColorAssociation[i - 1].dfVal != dfVal)
1688 : {
1689 252 : *pnR = 0;
1690 252 : *pnG = 0;
1691 252 : *pnB = 0;
1692 252 : *pnA = 0;
1693 252 : return false;
1694 : }
1695 : else
1696 : {
1697 758 : *pnR = asColorAssociation[i - 1].nR;
1698 758 : *pnG = asColorAssociation[i - 1].nG;
1699 758 : *pnB = asColorAssociation[i - 1].nB;
1700 758 : *pnA = asColorAssociation[i - 1].nA;
1701 758 : return true;
1702 : }
1703 : }
1704 : else
1705 : {
1706 293602 : if (asColorAssociation[i - 1].dfVal == dfVal)
1707 : {
1708 0 : *pnR = asColorAssociation[i - 1].nR;
1709 0 : *pnG = asColorAssociation[i - 1].nG;
1710 0 : *pnB = asColorAssociation[i - 1].nB;
1711 0 : *pnA = asColorAssociation[i - 1].nA;
1712 0 : return true;
1713 : }
1714 :
1715 293602 : if (asColorAssociation[i].dfVal == dfVal)
1716 : {
1717 93543 : *pnR = asColorAssociation[i].nR;
1718 93543 : *pnG = asColorAssociation[i].nG;
1719 93543 : *pnB = asColorAssociation[i].nB;
1720 93543 : *pnA = asColorAssociation[i].nA;
1721 93543 : return true;
1722 : }
1723 :
1724 200059 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
1725 : {
1726 20182 : *pnR = 0;
1727 20182 : *pnG = 0;
1728 20182 : *pnB = 0;
1729 20182 : *pnA = 0;
1730 20182 : return false;
1731 : }
1732 :
1733 210024 : if (eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY &&
1734 30147 : asColorAssociation[i - 1].dfVal != dfVal)
1735 : {
1736 30147 : const size_t index = (dfVal - asColorAssociation[i - 1].dfVal <
1737 30147 : asColorAssociation[i].dfVal - dfVal)
1738 30147 : ? i - 1
1739 30147 : : i;
1740 30147 : *pnR = asColorAssociation[index].nR;
1741 30147 : *pnG = asColorAssociation[index].nG;
1742 30147 : *pnB = asColorAssociation[index].nB;
1743 30147 : *pnA = asColorAssociation[index].nA;
1744 30147 : return true;
1745 : }
1746 :
1747 149730 : if (std::isnan(asColorAssociation[i - 1].dfVal))
1748 : {
1749 0 : *pnR = asColorAssociation[i].nR;
1750 0 : *pnG = asColorAssociation[i].nG;
1751 0 : *pnB = asColorAssociation[i].nB;
1752 0 : *pnA = asColorAssociation[i].nA;
1753 0 : return true;
1754 : }
1755 :
1756 : const double dfRatio =
1757 149730 : (dfVal - asColorAssociation[i - 1].dfVal) /
1758 149730 : (asColorAssociation[i].dfVal - asColorAssociation[i - 1].dfVal);
1759 598920 : const auto LinearInterpolation = [dfRatio](int nValBefore, int nVal)
1760 : {
1761 1197840 : return std::clamp(static_cast<int>(0.5 + nValBefore +
1762 598920 : dfRatio * (nVal - nValBefore)),
1763 598920 : 0, 255);
1764 149730 : };
1765 :
1766 149730 : *pnR = LinearInterpolation(asColorAssociation[i - 1].nR,
1767 149730 : asColorAssociation[i].nR);
1768 149730 : *pnG = LinearInterpolation(asColorAssociation[i - 1].nG,
1769 149730 : asColorAssociation[i].nG);
1770 149730 : *pnB = LinearInterpolation(asColorAssociation[i - 1].nB,
1771 149730 : asColorAssociation[i].nB);
1772 149730 : *pnA = LinearInterpolation(asColorAssociation[i - 1].nA,
1773 149730 : asColorAssociation[i].nA);
1774 :
1775 149730 : return true;
1776 : }
1777 : }
1778 :
1779 : static std::vector<GDALColorAssociation>
1780 41 : GDALColorReliefParseColorFile(GDALRasterBandH hSrcBand,
1781 : const char *pszColorFilename,
1782 : ColorSelectionMode eColorSelectionMode)
1783 : {
1784 : std::vector<GDALColorAssociation> asColorAssociation = GDALLoadTextColorMap(
1785 82 : pszColorFilename, GDALRasterBand::FromHandle(hSrcBand));
1786 41 : if (asColorAssociation.empty())
1787 : {
1788 2 : return {};
1789 : }
1790 :
1791 39 : int bSrcHasNoData = FALSE;
1792 : const double dfSrcNoDataValue =
1793 39 : GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
1794 :
1795 39 : GDALColorReliefProcessColors(asColorAssociation, bSrcHasNoData,
1796 : dfSrcNoDataValue, eColorSelectionMode);
1797 :
1798 39 : return asColorAssociation;
1799 : }
1800 :
1801 25 : static GByte *GDALColorReliefPrecompute(
1802 : GDALRasterBandH hSrcBand,
1803 : const std::vector<GDALColorAssociation> &asColorAssociation,
1804 : ColorSelectionMode eColorSelectionMode, int *pnIndexOffset)
1805 : {
1806 25 : const GDALDataType eDT = GDALGetRasterDataType(hSrcBand);
1807 25 : GByte *pabyPrecomputed = nullptr;
1808 25 : const int nIndexOffset = (eDT == GDT_Int16) ? 32768 : 0;
1809 25 : *pnIndexOffset = nIndexOffset;
1810 25 : const int nXSize = GDALGetRasterBandXSize(hSrcBand);
1811 25 : const int nYSize = GDALGetRasterBandYSize(hSrcBand);
1812 25 : if (eDT == GDT_Byte || ((eDT == GDT_Int16 || eDT == GDT_UInt16) &&
1813 14 : static_cast<GIntBig>(nXSize) * nYSize > 65536))
1814 : {
1815 7 : const int iMax = (eDT == GDT_Byte) ? 256 : 65536;
1816 7 : pabyPrecomputed = static_cast<GByte *>(VSI_MALLOC2_VERBOSE(4, iMax));
1817 7 : if (pabyPrecomputed)
1818 : {
1819 1799 : for (int i = 0; i < iMax; i++)
1820 : {
1821 1792 : int nR = 0;
1822 1792 : int nG = 0;
1823 1792 : int nB = 0;
1824 1792 : int nA = 0;
1825 1792 : GDALColorReliefGetRGBA(asColorAssociation, i - nIndexOffset,
1826 : eColorSelectionMode, &nR, &nG, &nB, &nA);
1827 1792 : pabyPrecomputed[4 * i] = static_cast<GByte>(nR);
1828 1792 : pabyPrecomputed[4 * i + 1] = static_cast<GByte>(nG);
1829 1792 : pabyPrecomputed[4 * i + 2] = static_cast<GByte>(nB);
1830 1792 : pabyPrecomputed[4 * i + 3] = static_cast<GByte>(nA);
1831 : }
1832 : }
1833 : }
1834 25 : return pabyPrecomputed;
1835 : }
1836 :
1837 : /************************************************************************/
1838 : /* ==================================================================== */
1839 : /* GDALColorReliefDataset */
1840 : /* ==================================================================== */
1841 : /************************************************************************/
1842 :
1843 : class GDALColorReliefRasterBand;
1844 :
1845 : class GDALColorReliefDataset : public GDALDataset
1846 : {
1847 : friend class GDALColorReliefRasterBand;
1848 :
1849 : GDALDatasetH hSrcDS;
1850 : GDALRasterBandH hSrcBand;
1851 : std::vector<GDALColorAssociation> asColorAssociation{};
1852 : ColorSelectionMode eColorSelectionMode;
1853 : GByte *pabyPrecomputed;
1854 : int nIndexOffset;
1855 : float *pafSourceBuf;
1856 : int *panSourceBuf;
1857 : int nCurBlockXOff;
1858 : int nCurBlockYOff;
1859 :
1860 : CPL_DISALLOW_COPY_ASSIGN(GDALColorReliefDataset)
1861 :
1862 : public:
1863 : GDALColorReliefDataset(GDALDatasetH hSrcDS, GDALRasterBandH hSrcBand,
1864 : const char *pszColorFilename,
1865 : ColorSelectionMode eColorSelectionMode, int bAlpha);
1866 : ~GDALColorReliefDataset();
1867 :
1868 2 : bool InitOK() const
1869 : {
1870 2 : return pafSourceBuf != nullptr || panSourceBuf != nullptr;
1871 : }
1872 :
1873 : CPLErr GetGeoTransform(double *padfGeoTransform) override;
1874 : const OGRSpatialReference *GetSpatialRef() const override;
1875 : };
1876 :
1877 : /************************************************************************/
1878 : /* ==================================================================== */
1879 : /* GDALColorReliefRasterBand */
1880 : /* ==================================================================== */
1881 : /************************************************************************/
1882 :
1883 : class GDALColorReliefRasterBand : public GDALRasterBand
1884 : {
1885 : friend class GDALColorReliefDataset;
1886 :
1887 : public:
1888 : GDALColorReliefRasterBand(GDALColorReliefDataset *, int);
1889 :
1890 : virtual CPLErr IReadBlock(int, int, void *) override;
1891 : virtual GDALColorInterp GetColorInterpretation() override;
1892 : };
1893 :
1894 2 : GDALColorReliefDataset::GDALColorReliefDataset(
1895 : GDALDatasetH hSrcDSIn, GDALRasterBandH hSrcBandIn,
1896 : const char *pszColorFilename, ColorSelectionMode eColorSelectionModeIn,
1897 2 : int bAlpha)
1898 : : hSrcDS(hSrcDSIn), hSrcBand(hSrcBandIn),
1899 : eColorSelectionMode(eColorSelectionModeIn), pabyPrecomputed(nullptr),
1900 : nIndexOffset(0), pafSourceBuf(nullptr), panSourceBuf(nullptr),
1901 2 : nCurBlockXOff(-1), nCurBlockYOff(-1)
1902 : {
1903 4 : asColorAssociation = GDALColorReliefParseColorFile(
1904 2 : hSrcBand, pszColorFilename, eColorSelectionMode);
1905 :
1906 2 : nRasterXSize = GDALGetRasterXSize(hSrcDS);
1907 2 : nRasterYSize = GDALGetRasterYSize(hSrcDS);
1908 :
1909 2 : int nBlockXSize = 0;
1910 2 : int nBlockYSize = 0;
1911 2 : GDALGetBlockSize(hSrcBand, &nBlockXSize, &nBlockYSize);
1912 :
1913 4 : pabyPrecomputed = GDALColorReliefPrecompute(
1914 2 : hSrcBand, asColorAssociation, eColorSelectionMode, &nIndexOffset);
1915 :
1916 8 : for (int i = 0; i < ((bAlpha) ? 4 : 3); i++)
1917 : {
1918 6 : SetBand(i + 1, new GDALColorReliefRasterBand(this, i + 1));
1919 : }
1920 :
1921 2 : if (pabyPrecomputed)
1922 0 : panSourceBuf = static_cast<int *>(
1923 0 : VSI_MALLOC3_VERBOSE(sizeof(int), nBlockXSize, nBlockYSize));
1924 : else
1925 2 : pafSourceBuf = static_cast<float *>(
1926 2 : VSI_MALLOC3_VERBOSE(sizeof(float), nBlockXSize, nBlockYSize));
1927 2 : }
1928 :
1929 4 : GDALColorReliefDataset::~GDALColorReliefDataset()
1930 : {
1931 2 : CPLFree(pabyPrecomputed);
1932 2 : CPLFree(panSourceBuf);
1933 2 : CPLFree(pafSourceBuf);
1934 4 : }
1935 :
1936 2 : CPLErr GDALColorReliefDataset::GetGeoTransform(double *padfGeoTransform)
1937 : {
1938 2 : return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
1939 : }
1940 :
1941 2 : const OGRSpatialReference *GDALColorReliefDataset::GetSpatialRef() const
1942 : {
1943 2 : return GDALDataset::FromHandle(hSrcDS)->GetSpatialRef();
1944 : }
1945 :
1946 6 : GDALColorReliefRasterBand::GDALColorReliefRasterBand(
1947 6 : GDALColorReliefDataset *poDSIn, int nBandIn)
1948 : {
1949 6 : poDS = poDSIn;
1950 6 : nBand = nBandIn;
1951 6 : eDataType = GDT_Byte;
1952 6 : GDALGetBlockSize(poDSIn->hSrcBand, &nBlockXSize, &nBlockYSize);
1953 6 : }
1954 :
1955 36 : CPLErr GDALColorReliefRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
1956 : void *pImage)
1957 : {
1958 : GDALColorReliefDataset *poGDS =
1959 36 : cpl::down_cast<GDALColorReliefDataset *>(poDS);
1960 72 : const int nReqXSize = (nBlockXOff + 1) * nBlockXSize >= nRasterXSize
1961 36 : ? nRasterXSize - nBlockXOff * nBlockXSize
1962 : : nBlockXSize;
1963 :
1964 72 : const int nReqYSize = (nBlockYOff + 1) * nBlockYSize >= nRasterYSize
1965 36 : ? nRasterYSize - nBlockYOff * nBlockYSize
1966 : : nBlockYSize;
1967 :
1968 36 : if (poGDS->nCurBlockXOff != nBlockXOff ||
1969 34 : poGDS->nCurBlockYOff != nBlockYOff)
1970 : {
1971 12 : poGDS->nCurBlockXOff = nBlockXOff;
1972 12 : poGDS->nCurBlockYOff = nBlockYOff;
1973 :
1974 24 : const CPLErr eErr = GDALRasterIO(
1975 12 : poGDS->hSrcBand, GF_Read, nBlockXOff * nBlockXSize,
1976 12 : nBlockYOff * nBlockYSize, nReqXSize, nReqYSize,
1977 12 : (poGDS->panSourceBuf) ? static_cast<void *>(poGDS->panSourceBuf)
1978 : : static_cast<void *>(poGDS->pafSourceBuf),
1979 : nReqXSize, nReqYSize,
1980 12 : (poGDS->panSourceBuf) ? GDT_Int32 : GDT_Float32, 0, 0);
1981 12 : if (eErr != CE_None)
1982 : {
1983 0 : memset(pImage, 0, static_cast<size_t>(nBlockXSize) * nBlockYSize);
1984 0 : return eErr;
1985 : }
1986 : }
1987 :
1988 36 : int j = 0;
1989 36 : if (poGDS->panSourceBuf)
1990 : {
1991 0 : for (int y = 0; y < nReqYSize; y++)
1992 : {
1993 0 : for (int x = 0; x < nReqXSize; x++)
1994 : {
1995 0 : const int nIndex = poGDS->panSourceBuf[j] + poGDS->nIndexOffset;
1996 0 : static_cast<GByte *>(pImage)[y * nBlockXSize + x] =
1997 0 : poGDS->pabyPrecomputed[4 * nIndex + nBand - 1];
1998 0 : j++;
1999 : }
2000 : }
2001 : }
2002 : else
2003 : {
2004 36 : int anComponents[4] = {0, 0, 0, 0};
2005 762 : for (int y = 0; y < nReqYSize; y++)
2006 : {
2007 88572 : for (int x = 0; x < nReqXSize; x++)
2008 : {
2009 87846 : GDALColorReliefGetRGBA(
2010 87846 : poGDS->asColorAssociation, poGDS->pafSourceBuf[j],
2011 : poGDS->eColorSelectionMode, &anComponents[0],
2012 : &anComponents[1], &anComponents[2], &anComponents[3]);
2013 87846 : static_cast<GByte *>(pImage)[y * nBlockXSize + x] =
2014 87846 : static_cast<GByte>(anComponents[nBand - 1]);
2015 87846 : j++;
2016 : }
2017 : }
2018 : }
2019 :
2020 36 : return CE_None;
2021 : }
2022 :
2023 12 : GDALColorInterp GDALColorReliefRasterBand::GetColorInterpretation()
2024 : {
2025 12 : return static_cast<GDALColorInterp>(GCI_RedBand + nBand - 1);
2026 : }
2027 :
2028 : static CPLErr
2029 25 : GDALColorRelief(GDALRasterBandH hSrcBand, GDALRasterBandH hDstBand1,
2030 : GDALRasterBandH hDstBand2, GDALRasterBandH hDstBand3,
2031 : GDALRasterBandH hDstBand4, const char *pszColorFilename,
2032 : ColorSelectionMode eColorSelectionMode,
2033 : GDALProgressFunc pfnProgress, void *pProgressData)
2034 : {
2035 25 : if (hSrcBand == nullptr || hDstBand1 == nullptr || hDstBand2 == nullptr ||
2036 : hDstBand3 == nullptr)
2037 0 : return CE_Failure;
2038 :
2039 : const auto asColorAssociation = GDALColorReliefParseColorFile(
2040 50 : hSrcBand, pszColorFilename, eColorSelectionMode);
2041 25 : if (asColorAssociation.empty())
2042 2 : return CE_Failure;
2043 :
2044 23 : if (pfnProgress == nullptr)
2045 16 : pfnProgress = GDALDummyProgress;
2046 :
2047 : /* -------------------------------------------------------------------- */
2048 : /* Precompute the map from values to RGBA quadruplets */
2049 : /* for GDT_Byte, GDT_Int16 or GDT_UInt16 */
2050 : /* -------------------------------------------------------------------- */
2051 23 : int nIndexOffset = 0;
2052 : std::unique_ptr<GByte, VSIFreeReleaser> pabyPrecomputed(
2053 : GDALColorReliefPrecompute(hSrcBand, asColorAssociation,
2054 46 : eColorSelectionMode, &nIndexOffset));
2055 :
2056 : /* -------------------------------------------------------------------- */
2057 : /* Initialize progress counter. */
2058 : /* -------------------------------------------------------------------- */
2059 :
2060 23 : const int nXSize = GDALGetRasterBandXSize(hSrcBand);
2061 23 : const int nYSize = GDALGetRasterBandYSize(hSrcBand);
2062 :
2063 23 : std::unique_ptr<float, VSIFreeReleaser> pafSourceBuf;
2064 23 : std::unique_ptr<int, VSIFreeReleaser> panSourceBuf;
2065 23 : if (pabyPrecomputed)
2066 7 : panSourceBuf.reset(
2067 7 : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nXSize)));
2068 : else
2069 16 : pafSourceBuf.reset(
2070 16 : static_cast<float *>(VSI_MALLOC2_VERBOSE(sizeof(float), nXSize)));
2071 : std::unique_ptr<GByte, VSIFreeReleaser> pabyDestBuf(
2072 46 : static_cast<GByte *>(VSI_MALLOC2_VERBOSE(4, nXSize)));
2073 23 : GByte *pabyDestBuf1 = pabyDestBuf.get();
2074 23 : GByte *pabyDestBuf2 = pabyDestBuf1 ? pabyDestBuf1 + nXSize : nullptr;
2075 23 : GByte *pabyDestBuf3 = pabyDestBuf2 ? pabyDestBuf2 + nXSize : nullptr;
2076 23 : GByte *pabyDestBuf4 = pabyDestBuf3 ? pabyDestBuf3 + nXSize : nullptr;
2077 :
2078 53 : if ((pabyPrecomputed != nullptr && panSourceBuf == nullptr) ||
2079 53 : (pabyPrecomputed == nullptr && pafSourceBuf == nullptr) ||
2080 : pabyDestBuf1 == nullptr)
2081 : {
2082 0 : return CE_Failure;
2083 : }
2084 :
2085 23 : if (!pfnProgress(0.0, nullptr, pProgressData))
2086 : {
2087 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
2088 0 : return CE_Failure;
2089 : }
2090 :
2091 23 : int nR = 0;
2092 23 : int nG = 0;
2093 23 : int nB = 0;
2094 23 : int nA = 0;
2095 :
2096 1729 : for (int i = 0; i < nYSize; i++)
2097 : {
2098 : /* Read source buffer */
2099 5118 : CPLErr eErr = GDALRasterIO(
2100 : hSrcBand, GF_Read, 0, i, nXSize, 1,
2101 7 : panSourceBuf ? static_cast<void *>(panSourceBuf.get())
2102 3405 : : static_cast<void *>(pafSourceBuf.get()),
2103 1706 : nXSize, 1, panSourceBuf ? GDT_Int32 : GDT_Float32, 0, 0);
2104 1706 : if (eErr != CE_None)
2105 : {
2106 0 : return eErr;
2107 : }
2108 :
2109 1706 : if (pabyPrecomputed)
2110 : {
2111 7 : const auto pabyPrecomputedRaw = pabyPrecomputed.get();
2112 7 : const auto panSourceBufRaw = panSourceBuf.get();
2113 32 : for (int j = 0; j < nXSize; j++)
2114 : {
2115 25 : int nIndex = panSourceBufRaw[j] + nIndexOffset;
2116 25 : pabyDestBuf1[j] = pabyPrecomputedRaw[4 * nIndex];
2117 25 : pabyDestBuf2[j] = pabyPrecomputedRaw[4 * nIndex + 1];
2118 25 : pabyDestBuf3[j] = pabyPrecomputedRaw[4 * nIndex + 2];
2119 25 : pabyDestBuf4[j] = pabyPrecomputedRaw[4 * nIndex + 3];
2120 : }
2121 : }
2122 : else
2123 : {
2124 1699 : const auto pafSourceBufRaw = pafSourceBuf.get();
2125 206683 : for (int j = 0; j < nXSize; j++)
2126 : {
2127 204984 : GDALColorReliefGetRGBA(asColorAssociation, pafSourceBufRaw[j],
2128 : eColorSelectionMode, &nR, &nG, &nB, &nA);
2129 204984 : pabyDestBuf1[j] = static_cast<GByte>(nR);
2130 204984 : pabyDestBuf2[j] = static_cast<GByte>(nG);
2131 204984 : pabyDestBuf3[j] = static_cast<GByte>(nB);
2132 204984 : pabyDestBuf4[j] = static_cast<GByte>(nA);
2133 : }
2134 : }
2135 :
2136 : /* -----------------------------------------
2137 : * Write Line to Raster
2138 : */
2139 1706 : eErr = GDALRasterIO(hDstBand1, GF_Write, 0, i, nXSize, 1, pabyDestBuf1,
2140 : nXSize, 1, GDT_Byte, 0, 0);
2141 1706 : if (eErr == CE_None)
2142 : {
2143 1706 : eErr = GDALRasterIO(hDstBand2, GF_Write, 0, i, nXSize, 1,
2144 : pabyDestBuf2, nXSize, 1, GDT_Byte, 0, 0);
2145 : }
2146 1706 : if (eErr == CE_None)
2147 : {
2148 1706 : eErr = GDALRasterIO(hDstBand3, GF_Write, 0, i, nXSize, 1,
2149 : pabyDestBuf3, nXSize, 1, GDT_Byte, 0, 0);
2150 : }
2151 1706 : if (eErr == CE_None && hDstBand4)
2152 : {
2153 242 : eErr = GDALRasterIO(hDstBand4, GF_Write, 0, i, nXSize, 1,
2154 : pabyDestBuf4, nXSize, 1, GDT_Byte, 0, 0);
2155 : }
2156 :
2157 3412 : if (eErr == CE_None &&
2158 1706 : !pfnProgress(1.0 * (i + 1) / nYSize, nullptr, pProgressData))
2159 : {
2160 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
2161 0 : eErr = CE_Failure;
2162 : }
2163 1706 : if (eErr != CE_None)
2164 : {
2165 0 : return eErr;
2166 : }
2167 : }
2168 :
2169 23 : pfnProgress(1.0, nullptr, pProgressData);
2170 :
2171 23 : return CE_None;
2172 : }
2173 :
2174 : /************************************************************************/
2175 : /* GDALGenerateVRTColorRelief() */
2176 : /************************************************************************/
2177 :
2178 14 : static bool GDALGenerateVRTColorRelief(const char *pszDstFilename,
2179 : GDALDatasetH hSrcDataset,
2180 : GDALRasterBandH hSrcBand,
2181 : const char *pszColorFilename,
2182 : ColorSelectionMode eColorSelectionMode,
2183 : bool bAddAlpha)
2184 : {
2185 : const auto asColorAssociation = GDALColorReliefParseColorFile(
2186 28 : hSrcBand, pszColorFilename, eColorSelectionMode);
2187 14 : if (asColorAssociation.empty())
2188 0 : return false;
2189 :
2190 14 : VSILFILE *fp = VSIFOpenL(pszDstFilename, "wt");
2191 14 : if (fp == nullptr)
2192 : {
2193 0 : return false;
2194 : }
2195 :
2196 14 : const int nXSize = GDALGetRasterBandXSize(hSrcBand);
2197 14 : const int nYSize = GDALGetRasterBandYSize(hSrcBand);
2198 :
2199 : bool bOK =
2200 14 : VSIFPrintfL(fp, "<VRTDataset rasterXSize=\"%d\" rasterYSize=\"%d\">\n",
2201 14 : nXSize, nYSize) > 0;
2202 14 : const char *pszProjectionRef = GDALGetProjectionRef(hSrcDataset);
2203 14 : if (pszProjectionRef && pszProjectionRef[0] != '\0')
2204 : {
2205 : char *pszEscapedString =
2206 7 : CPLEscapeString(pszProjectionRef, -1, CPLES_XML);
2207 7 : bOK &= VSIFPrintfL(fp, " <SRS>%s</SRS>\n", pszEscapedString) > 0;
2208 7 : VSIFree(pszEscapedString);
2209 : }
2210 14 : double adfGT[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2211 14 : if (GDALGetGeoTransform(hSrcDataset, adfGT) == CE_None)
2212 : {
2213 8 : bOK &= VSIFPrintfL(fp,
2214 : " <GeoTransform> %.16g, %.16g, %.16g, "
2215 : "%.16g, %.16g, %.16g</GeoTransform>\n",
2216 : adfGT[0], adfGT[1], adfGT[2], adfGT[3], adfGT[4],
2217 8 : adfGT[5]) > 0;
2218 : }
2219 :
2220 14 : int nBlockXSize = 0;
2221 14 : int nBlockYSize = 0;
2222 14 : GDALGetBlockSize(hSrcBand, &nBlockXSize, &nBlockYSize);
2223 :
2224 14 : int bRelativeToVRT = FALSE;
2225 14 : const CPLString osPath = CPLGetPathSafe(pszDstFilename);
2226 14 : char *pszSourceFilename = CPLStrdup(CPLExtractRelativePath(
2227 : osPath.c_str(), GDALGetDescription(hSrcDataset), &bRelativeToVRT));
2228 :
2229 14 : const int nBands = 3 + (bAddAlpha ? 1 : 0);
2230 :
2231 57 : for (int iBand = 0; iBand < nBands; iBand++)
2232 : {
2233 43 : bOK &=
2234 43 : VSIFPrintfL(fp, " <VRTRasterBand dataType=\"Byte\" band=\"%d\">\n",
2235 43 : iBand + 1) > 0;
2236 43 : bOK &= VSIFPrintfL(
2237 : fp, " <ColorInterp>%s</ColorInterp>\n",
2238 : GDALGetColorInterpretationName(
2239 43 : static_cast<GDALColorInterp>(GCI_RedBand + iBand))) > 0;
2240 43 : bOK &= VSIFPrintfL(fp, " <ComplexSource>\n") > 0;
2241 43 : bOK &= VSIFPrintfL(fp,
2242 : " <SourceFilename "
2243 : "relativeToVRT=\"%d\">%s</SourceFilename>\n",
2244 43 : bRelativeToVRT, pszSourceFilename) > 0;
2245 43 : bOK &= VSIFPrintfL(fp, " <SourceBand>%d</SourceBand>\n",
2246 43 : GDALGetBandNumber(hSrcBand)) > 0;
2247 43 : bOK &= VSIFPrintfL(fp,
2248 : " <SourceProperties RasterXSize=\"%d\" "
2249 : "RasterYSize=\"%d\" DataType=\"%s\" "
2250 : "BlockXSize=\"%d\" BlockYSize=\"%d\"/>\n",
2251 : nXSize, nYSize,
2252 : GDALGetDataTypeName(GDALGetRasterDataType(hSrcBand)),
2253 43 : nBlockXSize, nBlockYSize) > 0;
2254 43 : bOK &=
2255 43 : VSIFPrintfL(
2256 : fp,
2257 : " "
2258 : "<SrcRect xOff=\"0\" yOff=\"0\" xSize=\"%d\" ySize=\"%d\"/>\n",
2259 43 : nXSize, nYSize) > 0;
2260 43 : bOK &=
2261 43 : VSIFPrintfL(
2262 : fp,
2263 : " "
2264 : "<DstRect xOff=\"0\" yOff=\"0\" xSize=\"%d\" ySize=\"%d\"/>\n",
2265 43 : nXSize, nYSize) > 0;
2266 :
2267 43 : bOK &= VSIFPrintfL(fp, " <LUT>") > 0;
2268 :
2269 302 : for (size_t iColor = 0; iColor < asColorAssociation.size(); iColor++)
2270 : {
2271 259 : const double dfVal = asColorAssociation[iColor].dfVal;
2272 259 : if (iColor > 0)
2273 216 : bOK &= VSIFPrintfL(fp, ",") > 0;
2274 :
2275 216 : if (iColor > 0 &&
2276 475 : eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY &&
2277 : dfVal !=
2278 60 : std::nextafter(asColorAssociation[iColor - 1].dfVal,
2279 : std::numeric_limits<double>::infinity()))
2280 : {
2281 : const double dfMidVal =
2282 54 : (dfVal + asColorAssociation[iColor - 1].dfVal) / 2.0;
2283 54 : bOK &=
2284 162 : VSIFPrintfL(
2285 : fp, "%.18g:%d",
2286 : std::nextafter(
2287 54 : dfMidVal, -std::numeric_limits<double>::infinity()),
2288 90 : (iBand == 0) ? asColorAssociation[iColor - 1].nR
2289 18 : : (iBand == 1) ? asColorAssociation[iColor - 1].nG
2290 18 : : (iBand == 2) ? asColorAssociation[iColor - 1].nB
2291 0 : : asColorAssociation[iColor - 1].nA) > 0;
2292 108 : bOK &= VSIFPrintfL(
2293 : fp, ",%.18g:%d", dfMidVal,
2294 90 : (iBand == 0) ? asColorAssociation[iColor].nR
2295 18 : : (iBand == 1) ? asColorAssociation[iColor].nG
2296 18 : : (iBand == 2) ? asColorAssociation[iColor].nB
2297 54 : : asColorAssociation[iColor].nA) > 0;
2298 : }
2299 : else
2300 : {
2301 205 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
2302 : {
2303 45 : bOK &=
2304 45 : VSIFPrintfL(
2305 : fp, "%.18g:0,",
2306 : std::nextafter(
2307 : dfVal,
2308 90 : -std::numeric_limits<double>::infinity())) > 0;
2309 : }
2310 205 : if (dfVal != static_cast<double>(static_cast<int>(dfVal)))
2311 18 : bOK &= VSIFPrintfL(fp, "%.18g", dfVal) > 0;
2312 : else
2313 187 : bOK &= VSIFPrintfL(fp, "%d", static_cast<int>(dfVal)) > 0;
2314 410 : bOK &= VSIFPrintfL(
2315 : fp, ":%d",
2316 344 : (iBand == 0) ? asColorAssociation[iColor].nR
2317 66 : : (iBand == 1) ? asColorAssociation[iColor].nG
2318 66 : : (iBand == 2) ? asColorAssociation[iColor].nB
2319 212 : : asColorAssociation[iColor].nA) > 0;
2320 : }
2321 :
2322 259 : if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
2323 : {
2324 45 : bOK &= VSIFPrintfL(
2325 : fp, ",%.18g:0",
2326 : std::nextafter(
2327 : dfVal,
2328 45 : std::numeric_limits<double>::infinity())) > 0;
2329 : }
2330 : }
2331 43 : bOK &= VSIFPrintfL(fp, "</LUT>\n") > 0;
2332 :
2333 43 : bOK &= VSIFPrintfL(fp, " </ComplexSource>\n") > 0;
2334 43 : bOK &= VSIFPrintfL(fp, " </VRTRasterBand>\n") > 0;
2335 : }
2336 :
2337 14 : CPLFree(pszSourceFilename);
2338 :
2339 14 : bOK &= VSIFPrintfL(fp, "</VRTDataset>\n") > 0;
2340 :
2341 14 : VSIFCloseL(fp);
2342 :
2343 14 : return bOK;
2344 : }
2345 :
2346 : /************************************************************************/
2347 : /* GDALTRIAlg() */
2348 : /************************************************************************/
2349 :
2350 : // Implements Wilson et al. (2007), for bathymetric use cases
2351 : template <class T>
2352 28802 : static float GDALTRIAlgWilson(const T *afWin, float /*fDstNoDataValue*/,
2353 : const AlgorithmParameters * /*pData*/)
2354 : {
2355 : // Terrain Ruggedness is average difference in height
2356 28802 : return (std::abs(afWin[0] - afWin[4]) + std::abs(afWin[1] - afWin[4]) +
2357 28802 : std::abs(afWin[2] - afWin[4]) + std::abs(afWin[3] - afWin[4]) +
2358 28802 : std::abs(afWin[5] - afWin[4]) + std::abs(afWin[6] - afWin[4]) +
2359 28802 : std::abs(afWin[7] - afWin[4]) + std::abs(afWin[8] - afWin[4])) *
2360 28802 : 0.125f;
2361 : }
2362 :
2363 : // Implements Riley, S.J., De Gloria, S.D., Elliot, R. (1999): A Terrain
2364 : // Ruggedness that Quantifies Topographic Heterogeneity. Intermountain Journal
2365 : // of Science, Vol.5, No.1-4, pp.23-27 for terrestrial use cases
2366 : template <class T>
2367 86886 : static float GDALTRIAlgRiley(const T *afWin, float /*fDstNoDataValue*/,
2368 : const AlgorithmParameters * /*pData*/)
2369 : {
2370 695088 : const auto square = [](double x) { return x * x; };
2371 :
2372 : return static_cast<float>(
2373 86886 : std::sqrt(square(afWin[0] - afWin[4]) + square(afWin[1] - afWin[4]) +
2374 86886 : square(afWin[2] - afWin[4]) + square(afWin[3] - afWin[4]) +
2375 86886 : square(afWin[5] - afWin[4]) + square(afWin[6] - afWin[4]) +
2376 86886 : square(afWin[7] - afWin[4]) + square(afWin[8] - afWin[4])));
2377 : }
2378 :
2379 : /************************************************************************/
2380 : /* GDALTPIAlg() */
2381 : /************************************************************************/
2382 :
2383 : template <class T>
2384 72245 : static float GDALTPIAlg(const T *afWin, float /*fDstNoDataValue*/,
2385 : const AlgorithmParameters * /*pData*/)
2386 : {
2387 : // Terrain Position is the difference between
2388 : // The central cell and the mean of the surrounding cells
2389 72245 : return afWin[4] - ((afWin[0] + afWin[1] + afWin[2] + afWin[3] + afWin[5] +
2390 72245 : afWin[6] + afWin[7] + afWin[8]) *
2391 72245 : 0.125f);
2392 : }
2393 :
2394 : /************************************************************************/
2395 : /* GDALRoughnessAlg() */
2396 : /************************************************************************/
2397 :
2398 : template <class T>
2399 79687 : static float GDALRoughnessAlg(const T *afWin, float /*fDstNoDataValue*/,
2400 : const AlgorithmParameters * /*pData*/)
2401 : {
2402 : // Roughness is the largest difference between any two cells
2403 :
2404 79687 : T fRoughnessMin = afWin[0];
2405 79687 : T fRoughnessMax = afWin[0];
2406 :
2407 717183 : for (int k = 1; k < 9; k++)
2408 : {
2409 637496 : if (afWin[k] > fRoughnessMax)
2410 : {
2411 75885 : fRoughnessMax = afWin[k];
2412 : }
2413 637496 : if (afWin[k] < fRoughnessMin)
2414 : {
2415 141951 : fRoughnessMin = afWin[k];
2416 : }
2417 : }
2418 79687 : return static_cast<float>(fRoughnessMax - fRoughnessMin);
2419 : }
2420 :
2421 : /************************************************************************/
2422 : /* ==================================================================== */
2423 : /* GDALGeneric3x3Dataset */
2424 : /* ==================================================================== */
2425 : /************************************************************************/
2426 :
2427 : template <class T> class GDALGeneric3x3RasterBand;
2428 :
2429 : template <class T> class GDALGeneric3x3Dataset final : public GDALDataset
2430 : {
2431 : friend class GDALGeneric3x3RasterBand<T>;
2432 :
2433 : const typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlg;
2434 : const typename GDALGeneric3x3ProcessingAlg_multisample<T>::type
2435 : pfnAlg_multisample;
2436 : std::unique_ptr<AlgorithmParameters> pAlgData;
2437 : GDALDatasetH hSrcDS = nullptr;
2438 : GDALRasterBandH hSrcBand = nullptr;
2439 : std::array<T *, 3> apafSourceBuf = {nullptr, nullptr, nullptr};
2440 : std::array<bool, 3> abLineHasNoDataValue = {false, false, false};
2441 : std::unique_ptr<float, VSIFreeReleaser> pafOutputBuf{};
2442 : int bDstHasNoData = false;
2443 : double dfDstNoDataValue = 0;
2444 : int nCurLine = -1;
2445 : const bool bComputeAtEdges;
2446 : const bool bTakeReference;
2447 :
2448 : using GDALDatasetRefCountedPtr =
2449 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>;
2450 :
2451 : std::vector<GDALDatasetRefCountedPtr> m_apoOverviewDS{};
2452 :
2453 : CPL_DISALLOW_COPY_ASSIGN(GDALGeneric3x3Dataset)
2454 :
2455 : public:
2456 : GDALGeneric3x3Dataset(
2457 : GDALDatasetH hSrcDS, GDALRasterBandH hSrcBand,
2458 : GDALDataType eDstDataType, int bDstHasNoData, double dfDstNoDataValue,
2459 : typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlg,
2460 : typename GDALGeneric3x3ProcessingAlg_multisample<T>::type
2461 : pfnAlg_multisample,
2462 : std::unique_ptr<AlgorithmParameters> pAlgData, bool bComputeAtEdges,
2463 : bool bTakeReferenceIn);
2464 : ~GDALGeneric3x3Dataset();
2465 :
2466 53 : bool InitOK() const
2467 : {
2468 106 : return apafSourceBuf[0] != nullptr && apafSourceBuf[1] != nullptr &&
2469 106 : apafSourceBuf[2] != nullptr;
2470 : }
2471 :
2472 : CPLErr GetGeoTransform(double *padfGeoTransform) override;
2473 : const OGRSpatialReference *GetSpatialRef() const override;
2474 : };
2475 :
2476 : /************************************************************************/
2477 : /* ==================================================================== */
2478 : /* GDALGeneric3x3RasterBand */
2479 : /* ==================================================================== */
2480 : /************************************************************************/
2481 :
2482 : template <class T> class GDALGeneric3x3RasterBand final : public GDALRasterBand
2483 : {
2484 : friend class GDALGeneric3x3Dataset<T>;
2485 : int bSrcHasNoData = false;
2486 : T fSrcNoDataValue = 0;
2487 : bool bIsSrcNoDataNan = false;
2488 : GDALDataType eReadDT = GDT_Unknown;
2489 :
2490 : void InitWithNoData(void *pImage);
2491 :
2492 : public:
2493 : GDALGeneric3x3RasterBand(GDALGeneric3x3Dataset<T> *poDSIn,
2494 : GDALDataType eDstDataType);
2495 :
2496 : virtual CPLErr IReadBlock(int, int, void *) override;
2497 : virtual double GetNoDataValue(int *pbHasNoData) override;
2498 :
2499 58 : int GetOverviewCount() override
2500 : {
2501 58 : auto poGDS = cpl::down_cast<GDALGeneric3x3Dataset<T> *>(poDS);
2502 58 : return static_cast<int>(poGDS->m_apoOverviewDS.size());
2503 : }
2504 :
2505 32 : GDALRasterBand *GetOverview(int idx) override
2506 : {
2507 32 : auto poGDS = cpl::down_cast<GDALGeneric3x3Dataset<T> *>(poDS);
2508 28 : return idx >= 0 && idx < GetOverviewCount()
2509 60 : ? poGDS->m_apoOverviewDS[idx]->GetRasterBand(1)
2510 32 : : nullptr;
2511 : }
2512 : };
2513 :
2514 : template <class T>
2515 53 : GDALGeneric3x3Dataset<T>::GDALGeneric3x3Dataset(
2516 : GDALDatasetH hSrcDSIn, GDALRasterBandH hSrcBandIn,
2517 : GDALDataType eDstDataType, int bDstHasNoDataIn, double dfDstNoDataValueIn,
2518 : typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlgIn,
2519 : typename GDALGeneric3x3ProcessingAlg_multisample<T>::type
2520 : pfnAlg_multisampleIn,
2521 : std::unique_ptr<AlgorithmParameters> pAlgDataIn, bool bComputeAtEdgesIn,
2522 : bool bTakeReferenceIn)
2523 : : pfnAlg(pfnAlgIn), pfnAlg_multisample(pfnAlg_multisampleIn),
2524 53 : pAlgData(std::move(pAlgDataIn)), hSrcDS(hSrcDSIn), hSrcBand(hSrcBandIn),
2525 : bDstHasNoData(bDstHasNoDataIn), dfDstNoDataValue(dfDstNoDataValueIn),
2526 53 : bComputeAtEdges(bComputeAtEdgesIn), bTakeReference(bTakeReferenceIn)
2527 : {
2528 53 : CPLAssert(eDstDataType == GDT_Byte || eDstDataType == GDT_Float32);
2529 :
2530 53 : if (bTakeReference)
2531 45 : GDALReferenceDataset(hSrcDS);
2532 :
2533 53 : nRasterXSize = GDALGetRasterXSize(hSrcDS);
2534 53 : nRasterYSize = GDALGetRasterYSize(hSrcDS);
2535 :
2536 53 : SetBand(1, new GDALGeneric3x3RasterBand<T>(this, eDstDataType));
2537 :
2538 106 : apafSourceBuf[0] =
2539 53 : static_cast<T *>(VSI_MALLOC2_VERBOSE(sizeof(T), nRasterXSize));
2540 106 : apafSourceBuf[1] =
2541 53 : static_cast<T *>(VSI_MALLOC2_VERBOSE(sizeof(T), nRasterXSize));
2542 106 : apafSourceBuf[2] =
2543 53 : static_cast<T *>(VSI_MALLOC2_VERBOSE(sizeof(T), nRasterXSize));
2544 53 : if (pfnAlg_multisample && eDstDataType == GDT_Byte)
2545 : {
2546 8 : pafOutputBuf.reset(static_cast<float *>(
2547 4 : VSI_MALLOC2_VERBOSE(sizeof(float), nRasterXSize)));
2548 : }
2549 :
2550 53 : const int nOvrCount = GDALGetOverviewCount(hSrcBandIn);
2551 61 : for (int i = 0;
2552 61 : i < nOvrCount && m_apoOverviewDS.size() == static_cast<size_t>(i); ++i)
2553 : {
2554 8 : auto hOvrBand = GDALGetOverview(hSrcBandIn, i);
2555 8 : auto hOvrDS = GDALGetBandDataset(hOvrBand);
2556 8 : if (hOvrDS && hOvrDS != hSrcDSIn)
2557 : {
2558 16 : auto poOvrDS = std::make_unique<GDALGeneric3x3Dataset>(
2559 8 : hOvrDS, hOvrBand, eDstDataType, bDstHasNoData, dfDstNoDataValue,
2560 8 : pfnAlg, pfnAlg_multisampleIn,
2561 6 : pAlgData ? pAlgData->CreateScaledParameters(
2562 6 : static_cast<double>(nRasterXSize) /
2563 6 : GDALGetRasterXSize(hOvrDS),
2564 6 : static_cast<double>(nRasterYSize) /
2565 6 : GDALGetRasterYSize(hOvrDS))
2566 : : nullptr,
2567 8 : bComputeAtEdges, false);
2568 8 : if (poOvrDS->InitOK())
2569 : {
2570 8 : m_apoOverviewDS.emplace_back(poOvrDS.release());
2571 : }
2572 : }
2573 : }
2574 53 : }
2575 :
2576 106 : template <class T> GDALGeneric3x3Dataset<T>::~GDALGeneric3x3Dataset()
2577 : {
2578 53 : if (bTakeReference)
2579 45 : GDALReleaseDataset(hSrcDS);
2580 :
2581 53 : CPLFree(apafSourceBuf[0]);
2582 53 : CPLFree(apafSourceBuf[1]);
2583 53 : CPLFree(apafSourceBuf[2]);
2584 159 : }
2585 :
2586 : template <class T>
2587 32 : CPLErr GDALGeneric3x3Dataset<T>::GetGeoTransform(double *padfGeoTransform)
2588 : {
2589 32 : return GDALGetGeoTransform(hSrcDS, padfGeoTransform);
2590 : }
2591 :
2592 : template <class T>
2593 33 : const OGRSpatialReference *GDALGeneric3x3Dataset<T>::GetSpatialRef() const
2594 : {
2595 33 : return GDALDataset::FromHandle(hSrcDS)->GetSpatialRef();
2596 : }
2597 :
2598 : template <class T>
2599 53 : GDALGeneric3x3RasterBand<T>::GDALGeneric3x3RasterBand(
2600 53 : GDALGeneric3x3Dataset<T> *poDSIn, GDALDataType eDstDataType)
2601 : {
2602 53 : poDS = poDSIn;
2603 53 : nBand = 1;
2604 53 : eDataType = eDstDataType;
2605 53 : nBlockXSize = poDS->GetRasterXSize();
2606 53 : nBlockYSize = 1;
2607 :
2608 : const double dfNoDataValue =
2609 53 : GDALGetRasterNoDataValue(poDSIn->hSrcBand, &bSrcHasNoData);
2610 : if (std::numeric_limits<T>::is_integer)
2611 : {
2612 52 : eReadDT = GDT_Int32;
2613 52 : if (bSrcHasNoData)
2614 : {
2615 52 : GDALDataType eSrcDT = GDALGetRasterDataType(poDSIn->hSrcBand);
2616 52 : CPLAssert(eSrcDT == GDT_Byte || eSrcDT == GDT_UInt16 ||
2617 : eSrcDT == GDT_Int16);
2618 52 : const int nMinVal = (eSrcDT == GDT_Byte) ? 0
2619 : : (eSrcDT == GDT_UInt16) ? 0
2620 : : -32768;
2621 52 : const int nMaxVal = (eSrcDT == GDT_Byte) ? 255
2622 : : (eSrcDT == GDT_UInt16) ? 65535
2623 : : 32767;
2624 :
2625 52 : if (fabs(dfNoDataValue - floor(dfNoDataValue + 0.5)) < 1e-2 &&
2626 52 : dfNoDataValue >= nMinVal && dfNoDataValue <= nMaxVal)
2627 : {
2628 52 : fSrcNoDataValue = static_cast<T>(floor(dfNoDataValue + 0.5));
2629 : }
2630 : else
2631 : {
2632 0 : bSrcHasNoData = FALSE;
2633 : }
2634 : }
2635 : }
2636 : else
2637 : {
2638 1 : eReadDT = GDT_Float32;
2639 1 : fSrcNoDataValue = static_cast<T>(dfNoDataValue);
2640 1 : bIsSrcNoDataNan = bSrcHasNoData && std::isnan(dfNoDataValue);
2641 : }
2642 53 : }
2643 :
2644 : template <class T>
2645 20 : void GDALGeneric3x3RasterBand<T>::InitWithNoData(void *pImage)
2646 : {
2647 20 : auto poGDS = cpl::down_cast<GDALGeneric3x3Dataset<T> *>(poDS);
2648 20 : if (eDataType == GDT_Byte)
2649 : {
2650 732 : for (int j = 0; j < nBlockXSize; j++)
2651 726 : static_cast<GByte *>(pImage)[j] =
2652 726 : static_cast<GByte>(poGDS->dfDstNoDataValue);
2653 : }
2654 : else
2655 : {
2656 1708 : for (int j = 0; j < nBlockXSize; j++)
2657 1694 : static_cast<float *>(pImage)[j] =
2658 1694 : static_cast<float>(poGDS->dfDstNoDataValue);
2659 : }
2660 20 : }
2661 :
2662 : template <class T>
2663 4965 : CPLErr GDALGeneric3x3RasterBand<T>::IReadBlock(int /*nBlockXOff*/,
2664 : int nBlockYOff, void *pImage)
2665 : {
2666 4965 : auto poGDS = cpl::down_cast<GDALGeneric3x3Dataset<T> *>(poDS);
2667 :
2668 1753923 : const auto UpdateLineNoDataFlag = [this, poGDS](int iLine)
2669 : {
2670 4965 : if (bSrcHasNoData)
2671 : {
2672 4965 : poGDS->abLineHasNoDataValue[iLine] = false;
2673 576450 : for (int i = 0; i < nRasterXSize; ++i)
2674 : {
2675 : if constexpr (std::numeric_limits<T>::is_integer)
2676 : {
2677 556844 : if (poGDS->apafSourceBuf[iLine][i] == fSrcNoDataValue)
2678 : {
2679 0 : poGDS->abLineHasNoDataValue[iLine] = true;
2680 0 : break;
2681 : }
2682 : }
2683 : else
2684 : {
2685 29282 : if (poGDS->apafSourceBuf[iLine][i] == fSrcNoDataValue ||
2686 14641 : std::isnan(poGDS->apafSourceBuf[iLine][i]))
2687 : {
2688 0 : poGDS->abLineHasNoDataValue[iLine] = true;
2689 0 : break;
2690 : }
2691 : }
2692 : }
2693 : }
2694 : };
2695 :
2696 4965 : if (poGDS->bComputeAtEdges && nRasterXSize >= 2 && nRasterYSize >= 2)
2697 : {
2698 3755 : if (nBlockYOff == 0)
2699 : {
2700 105 : for (int i = 0; i < 2; i++)
2701 : {
2702 70 : CPLErr eErr = GDALRasterIO(
2703 : poGDS->hSrcBand, GF_Read, 0, i, nBlockXSize, 1,
2704 70 : poGDS->apafSourceBuf[i + 1], nBlockXSize, 1, eReadDT, 0, 0);
2705 70 : if (eErr != CE_None)
2706 : {
2707 0 : InitWithNoData(pImage);
2708 0 : return eErr;
2709 : }
2710 70 : UpdateLineNoDataFlag(i + 1);
2711 : }
2712 35 : poGDS->nCurLine = 0;
2713 :
2714 3790 : for (int j = 0; j < nRasterXSize; j++)
2715 : {
2716 3755 : int jmin = (j == 0) ? j : j - 1;
2717 3755 : int jmax = (j == nRasterXSize - 1) ? j : j + 1;
2718 :
2719 3755 : T afWin[9] = {INTERPOL(poGDS->apafSourceBuf[1][jmin],
2720 3755 : poGDS->apafSourceBuf[2][jmin],
2721 : bSrcHasNoData, fSrcNoDataValue),
2722 3755 : INTERPOL(poGDS->apafSourceBuf[1][j],
2723 3755 : poGDS->apafSourceBuf[2][j],
2724 : bSrcHasNoData, fSrcNoDataValue),
2725 3755 : INTERPOL(poGDS->apafSourceBuf[1][jmax],
2726 3755 : poGDS->apafSourceBuf[2][jmax],
2727 : bSrcHasNoData, fSrcNoDataValue),
2728 3755 : poGDS->apafSourceBuf[1][jmin],
2729 3755 : poGDS->apafSourceBuf[1][j],
2730 3755 : poGDS->apafSourceBuf[1][jmax],
2731 3755 : poGDS->apafSourceBuf[2][jmin],
2732 3755 : poGDS->apafSourceBuf[2][j],
2733 3755 : poGDS->apafSourceBuf[2][jmax]};
2734 :
2735 3755 : const float fVal = ComputeVal(
2736 3755 : CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue,
2737 3755 : bIsSrcNoDataNan, afWin,
2738 3755 : static_cast<float>(poGDS->dfDstNoDataValue), poGDS->pfnAlg,
2739 3755 : poGDS->pAlgData.get(), poGDS->bComputeAtEdges);
2740 :
2741 3755 : if (eDataType == GDT_Byte)
2742 485 : static_cast<GByte *>(pImage)[j] =
2743 485 : static_cast<GByte>(fVal + 0.5);
2744 : else
2745 3270 : static_cast<float *>(pImage)[j] = fVal;
2746 : }
2747 :
2748 35 : return CE_None;
2749 : }
2750 3720 : else if (nBlockYOff == nRasterYSize - 1)
2751 : {
2752 35 : if (poGDS->nCurLine != nRasterYSize - 2)
2753 : {
2754 0 : for (int i = 0; i < 2; i++)
2755 : {
2756 0 : CPLErr eErr = GDALRasterIO(
2757 0 : poGDS->hSrcBand, GF_Read, 0, nRasterYSize - 2 + i,
2758 0 : nBlockXSize, 1, poGDS->apafSourceBuf[i + 1],
2759 : nBlockXSize, 1, eReadDT, 0, 0);
2760 0 : if (eErr != CE_None)
2761 : {
2762 0 : InitWithNoData(pImage);
2763 0 : return eErr;
2764 : }
2765 0 : UpdateLineNoDataFlag(i + 1);
2766 : }
2767 : }
2768 :
2769 3790 : for (int j = 0; j < nRasterXSize; j++)
2770 : {
2771 3755 : int jmin = (j == 0) ? j : j - 1;
2772 3755 : int jmax = (j == nRasterXSize - 1) ? j : j + 1;
2773 :
2774 3755 : T afWin[9] = {poGDS->apafSourceBuf[1][jmin],
2775 3755 : poGDS->apafSourceBuf[1][j],
2776 3755 : poGDS->apafSourceBuf[1][jmax],
2777 3755 : poGDS->apafSourceBuf[2][jmin],
2778 3755 : poGDS->apafSourceBuf[2][j],
2779 3755 : poGDS->apafSourceBuf[2][jmax],
2780 3755 : INTERPOL(poGDS->apafSourceBuf[2][jmin],
2781 3755 : poGDS->apafSourceBuf[1][jmin],
2782 : bSrcHasNoData, fSrcNoDataValue),
2783 3755 : INTERPOL(poGDS->apafSourceBuf[2][j],
2784 3755 : poGDS->apafSourceBuf[1][j],
2785 : bSrcHasNoData, fSrcNoDataValue),
2786 3755 : INTERPOL(poGDS->apafSourceBuf[2][jmax],
2787 3755 : poGDS->apafSourceBuf[1][jmax],
2788 : bSrcHasNoData, fSrcNoDataValue)};
2789 :
2790 3755 : const float fVal = ComputeVal(
2791 3755 : CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue,
2792 3755 : bIsSrcNoDataNan, afWin,
2793 3755 : static_cast<float>(poGDS->dfDstNoDataValue), poGDS->pfnAlg,
2794 3755 : poGDS->pAlgData.get(), poGDS->bComputeAtEdges);
2795 :
2796 3755 : if (eDataType == GDT_Byte)
2797 485 : static_cast<GByte *>(pImage)[j] =
2798 485 : static_cast<GByte>(fVal + 0.5);
2799 : else
2800 3270 : static_cast<float *>(pImage)[j] = fVal;
2801 : }
2802 :
2803 35 : return CE_None;
2804 3685 : }
2805 : }
2806 1210 : else if (nBlockYOff == 0 || nBlockYOff == nRasterYSize - 1)
2807 : {
2808 20 : InitWithNoData(pImage);
2809 20 : return CE_None;
2810 : }
2811 :
2812 4875 : if (poGDS->nCurLine != nBlockYOff)
2813 : {
2814 4875 : if (poGDS->nCurLine + 1 == nBlockYOff)
2815 : {
2816 4865 : T *pafTmp = poGDS->apafSourceBuf[0];
2817 4865 : poGDS->apafSourceBuf[0] = poGDS->apafSourceBuf[1];
2818 4865 : poGDS->apafSourceBuf[1] = poGDS->apafSourceBuf[2];
2819 4865 : poGDS->apafSourceBuf[2] = pafTmp;
2820 :
2821 4865 : CPLErr eErr = GDALRasterIO(
2822 : poGDS->hSrcBand, GF_Read, 0, nBlockYOff + 1, nBlockXSize, 1,
2823 4865 : poGDS->apafSourceBuf[2], nBlockXSize, 1, eReadDT, 0, 0);
2824 :
2825 4865 : if (eErr != CE_None)
2826 : {
2827 0 : InitWithNoData(pImage);
2828 0 : return eErr;
2829 : }
2830 :
2831 4865 : UpdateLineNoDataFlag(2);
2832 : }
2833 : else
2834 : {
2835 40 : for (int i = 0; i < 3; i++)
2836 : {
2837 90 : const CPLErr eErr = GDALRasterIO(
2838 30 : poGDS->hSrcBand, GF_Read, 0, nBlockYOff + i - 1,
2839 30 : nBlockXSize, 1, poGDS->apafSourceBuf[i], nBlockXSize, 1,
2840 : eReadDT, 0, 0);
2841 30 : if (eErr != CE_None)
2842 : {
2843 0 : InitWithNoData(pImage);
2844 0 : return eErr;
2845 : }
2846 :
2847 30 : UpdateLineNoDataFlag(i);
2848 : }
2849 : }
2850 :
2851 4875 : poGDS->nCurLine = nBlockYOff;
2852 : }
2853 :
2854 4875 : if (poGDS->bComputeAtEdges && nRasterXSize >= 2)
2855 : {
2856 3685 : int j = 0;
2857 33165 : T afWin[9] = {
2858 3685 : INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j + 1],
2859 : bSrcHasNoData, fSrcNoDataValue),
2860 3685 : poGDS->apafSourceBuf[0][j],
2861 3685 : poGDS->apafSourceBuf[0][j + 1],
2862 3685 : INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j + 1],
2863 : bSrcHasNoData, fSrcNoDataValue),
2864 3685 : poGDS->apafSourceBuf[1][j],
2865 3685 : poGDS->apafSourceBuf[1][j + 1],
2866 3685 : INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j + 1],
2867 : bSrcHasNoData, fSrcNoDataValue),
2868 3685 : poGDS->apafSourceBuf[2][j],
2869 3685 : poGDS->apafSourceBuf[2][j + 1]};
2870 :
2871 : {
2872 3685 : const float fVal = ComputeVal(
2873 3685 : CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, bIsSrcNoDataNan,
2874 3685 : afWin, static_cast<float>(poGDS->dfDstNoDataValue),
2875 3685 : poGDS->pfnAlg, poGDS->pAlgData.get(), poGDS->bComputeAtEdges);
2876 :
2877 3685 : if (eDataType == GDT_Byte)
2878 475 : static_cast<GByte *>(pImage)[j] =
2879 475 : static_cast<GByte>(fVal + 0.5);
2880 : else
2881 3210 : static_cast<float *>(pImage)[j] = fVal;
2882 : }
2883 :
2884 3685 : j = nRasterXSize - 1;
2885 :
2886 3685 : afWin[0] = poGDS->apafSourceBuf[0][j - 1];
2887 3685 : afWin[1] = poGDS->apafSourceBuf[0][j];
2888 3685 : afWin[2] =
2889 3685 : INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j - 1],
2890 : bSrcHasNoData, fSrcNoDataValue);
2891 3685 : afWin[3] = poGDS->apafSourceBuf[1][j - 1];
2892 3685 : afWin[4] = poGDS->apafSourceBuf[1][j];
2893 3685 : afWin[5] =
2894 3685 : INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j - 1],
2895 : bSrcHasNoData, fSrcNoDataValue);
2896 3685 : afWin[6] = poGDS->apafSourceBuf[2][j - 1];
2897 3685 : afWin[7] = poGDS->apafSourceBuf[2][j];
2898 3685 : afWin[8] =
2899 3685 : INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j - 1],
2900 : bSrcHasNoData, fSrcNoDataValue);
2901 :
2902 3685 : const float fVal = ComputeVal(
2903 3685 : CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, bIsSrcNoDataNan, afWin,
2904 3685 : static_cast<float>(poGDS->dfDstNoDataValue), poGDS->pfnAlg,
2905 3685 : poGDS->pAlgData.get(), poGDS->bComputeAtEdges);
2906 :
2907 3685 : if (eDataType == GDT_Byte)
2908 475 : static_cast<GByte *>(pImage)[j] = static_cast<GByte>(fVal + 0.5);
2909 : else
2910 3685 : static_cast<float *>(pImage)[j] = fVal;
2911 : }
2912 : else
2913 : {
2914 1190 : if (eDataType == GDT_Byte)
2915 : {
2916 357 : static_cast<GByte *>(pImage)[0] =
2917 357 : static_cast<GByte>(poGDS->dfDstNoDataValue);
2918 357 : if (nBlockXSize > 1)
2919 357 : static_cast<GByte *>(pImage)[nBlockXSize - 1] =
2920 357 : static_cast<GByte>(poGDS->dfDstNoDataValue);
2921 : }
2922 : else
2923 : {
2924 833 : static_cast<float *>(pImage)[0] =
2925 833 : static_cast<float>(poGDS->dfDstNoDataValue);
2926 833 : if (nBlockXSize > 1)
2927 833 : static_cast<float *>(pImage)[nBlockXSize - 1] =
2928 833 : static_cast<float>(poGDS->dfDstNoDataValue);
2929 : }
2930 : }
2931 :
2932 4875 : int j = 1;
2933 10226 : if (poGDS->pfnAlg_multisample &&
2934 476 : (eDataType == GDT_Float32 || poGDS->pafOutputBuf) &&
2935 5827 : !poGDS->abLineHasNoDataValue[0] && !poGDS->abLineHasNoDataValue[1] &&
2936 476 : !poGDS->abLineHasNoDataValue[2])
2937 : {
2938 1428 : j = poGDS->pfnAlg_multisample(
2939 476 : poGDS->apafSourceBuf[0], poGDS->apafSourceBuf[1],
2940 476 : poGDS->apafSourceBuf[2], nRasterXSize, poGDS->pAlgData.get(),
2941 952 : poGDS->pafOutputBuf ? poGDS->pafOutputBuf.get()
2942 : : static_cast<float *>(pImage));
2943 :
2944 476 : if (poGDS->pafOutputBuf)
2945 : {
2946 476 : GDALCopyWords64(poGDS->pafOutputBuf.get() + 1, GDT_Float32,
2947 : static_cast<int>(sizeof(float)),
2948 : static_cast<GByte *>(pImage) + 1, GDT_Byte, 1,
2949 476 : j - 1);
2950 : }
2951 : }
2952 :
2953 501464 : for (; j < nBlockXSize - 1; j++)
2954 : {
2955 4469303 : T afWin[9] = {
2956 496589 : poGDS->apafSourceBuf[0][j - 1], poGDS->apafSourceBuf[0][j],
2957 496589 : poGDS->apafSourceBuf[0][j + 1], poGDS->apafSourceBuf[1][j - 1],
2958 496589 : poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j + 1],
2959 496589 : poGDS->apafSourceBuf[2][j - 1], poGDS->apafSourceBuf[2][j],
2960 496589 : poGDS->apafSourceBuf[2][j + 1]};
2961 :
2962 496589 : const float fVal = ComputeVal(
2963 496589 : CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, bIsSrcNoDataNan, afWin,
2964 496589 : static_cast<float>(poGDS->dfDstNoDataValue), poGDS->pfnAlg,
2965 496589 : poGDS->pAlgData.get(), poGDS->bComputeAtEdges);
2966 :
2967 496589 : if (eDataType == GDT_Byte)
2968 36712 : static_cast<GByte *>(pImage)[j] = static_cast<GByte>(fVal + 0.5);
2969 : else
2970 459877 : static_cast<float *>(pImage)[j] = fVal;
2971 : }
2972 :
2973 4875 : return CE_None;
2974 : }
2975 :
2976 : template <class T>
2977 113 : double GDALGeneric3x3RasterBand<T>::GetNoDataValue(int *pbHasNoData)
2978 : {
2979 113 : auto poGDS = cpl::down_cast<GDALGeneric3x3Dataset<T> *>(poDS);
2980 113 : if (pbHasNoData)
2981 78 : *pbHasNoData = poGDS->bDstHasNoData;
2982 113 : return poGDS->dfDstNoDataValue;
2983 : }
2984 :
2985 : /************************************************************************/
2986 : /* GetAlgorithm() */
2987 : /************************************************************************/
2988 :
2989 : typedef enum
2990 : {
2991 : INVALID,
2992 : HILL_SHADE,
2993 : SLOPE,
2994 : ASPECT,
2995 : COLOR_RELIEF,
2996 : TRI,
2997 : TPI,
2998 : ROUGHNESS
2999 : } Algorithm;
3000 :
3001 161 : static Algorithm GetAlgorithm(const char *pszProcessing)
3002 : {
3003 161 : if (EQUAL(pszProcessing, "shade") || EQUAL(pszProcessing, "hillshade"))
3004 : {
3005 72 : return HILL_SHADE;
3006 : }
3007 89 : else if (EQUAL(pszProcessing, "slope"))
3008 : {
3009 16 : return SLOPE;
3010 : }
3011 73 : else if (EQUAL(pszProcessing, "aspect"))
3012 : {
3013 12 : return ASPECT;
3014 : }
3015 61 : else if (EQUAL(pszProcessing, "color-relief"))
3016 : {
3017 41 : return COLOR_RELIEF;
3018 : }
3019 20 : else if (EQUAL(pszProcessing, "TRI"))
3020 : {
3021 8 : return TRI;
3022 : }
3023 12 : else if (EQUAL(pszProcessing, "TPI"))
3024 : {
3025 5 : return TPI;
3026 : }
3027 7 : else if (EQUAL(pszProcessing, "roughness"))
3028 : {
3029 7 : return ROUGHNESS;
3030 : }
3031 : else
3032 : {
3033 0 : return INVALID;
3034 : }
3035 : }
3036 :
3037 : /************************************************************************/
3038 : /* GDALDEMAppOptionsGetParser() */
3039 : /************************************************************************/
3040 :
3041 166 : static std::unique_ptr<GDALArgumentParser> GDALDEMAppOptionsGetParser(
3042 : GDALDEMProcessingOptions *psOptions,
3043 : GDALDEMProcessingOptionsForBinary *psOptionsForBinary)
3044 : {
3045 : auto argParser = std::make_unique<GDALArgumentParser>(
3046 166 : "gdaldem", /* bForBinary=*/psOptionsForBinary != nullptr);
3047 :
3048 166 : argParser->add_description(_("Tools to analyze and visualize DEMs."));
3049 :
3050 166 : argParser->add_epilog(_("For more details, consult "
3051 166 : "https://gdal.org/programs/gdaldem.html"));
3052 :
3053 : // Common options (for all modes)
3054 : auto addCommonOptions =
3055 5964 : [psOptions, psOptionsForBinary](GDALArgumentParser *subParser)
3056 : {
3057 1162 : subParser->add_output_format_argument(psOptions->osFormat);
3058 :
3059 1162 : subParser->add_argument("-compute_edges")
3060 1162 : .flag()
3061 1162 : .store_into(psOptions->bComputeAtEdges)
3062 : .help(_(
3063 1162 : "Do the computation at raster edges and near nodata values."));
3064 :
3065 1162 : auto &bandArg = subParser->add_argument("-b")
3066 2324 : .metavar("<value>")
3067 1162 : .scan<'i', int>()
3068 1162 : .store_into(psOptions->nBand)
3069 1162 : .help(_("Select an input band."));
3070 :
3071 1162 : subParser->add_hidden_alias_for(bandArg, "--b");
3072 :
3073 1162 : subParser->add_creation_options_argument(psOptions->aosCreationOptions);
3074 :
3075 1162 : if (psOptionsForBinary)
3076 : {
3077 154 : subParser->add_quiet_argument(&psOptionsForBinary->bQuiet);
3078 : }
3079 1162 : };
3080 :
3081 : // Hillshade
3082 :
3083 166 : auto subCommandHillshade = argParser->add_subparser(
3084 : "hillshade", /* bForBinary=*/psOptionsForBinary != nullptr);
3085 166 : subCommandHillshade->add_description(_("Compute hillshade."));
3086 :
3087 166 : if (psOptionsForBinary)
3088 : {
3089 22 : subCommandHillshade->add_argument("input_dem")
3090 22 : .store_into(psOptionsForBinary->osSrcFilename)
3091 22 : .help(_("The input DEM raster to be processed."));
3092 :
3093 22 : subCommandHillshade->add_argument("output_hillshade")
3094 22 : .store_into(psOptionsForBinary->osDstFilename)
3095 22 : .help(_("The output raster to be produced."));
3096 : }
3097 :
3098 166 : subCommandHillshade->add_argument("-alg")
3099 332 : .metavar("<Horn|ZevenbergenThorne>")
3100 : .action(
3101 105 : [psOptions](const std::string &s)
3102 : {
3103 56 : if (EQUAL(s.c_str(), "ZevenbergenThorne"))
3104 : {
3105 20 : psOptions->bGradientAlgSpecified = true;
3106 20 : psOptions->eGradientAlg = GradientAlg::ZEVENBERGEN_THORNE;
3107 : }
3108 36 : else if (EQUAL(s.c_str(), "Horn"))
3109 : {
3110 29 : psOptions->bGradientAlgSpecified = true;
3111 29 : psOptions->eGradientAlg = GradientAlg::HORN;
3112 : }
3113 : else
3114 : {
3115 : throw std::invalid_argument(
3116 7 : CPLSPrintf("Invalid value for -alg: %s.", s.c_str()));
3117 : }
3118 215 : })
3119 166 : .help(_("Choose between ZevenbergenThorne or Horn algorithms."));
3120 :
3121 166 : subCommandHillshade->add_argument("-z")
3122 332 : .metavar("<factor>")
3123 166 : .scan<'g', double>()
3124 166 : .store_into(psOptions->z)
3125 166 : .help(_("Vertical exaggeration."));
3126 :
3127 : auto &scaleHillshadeArg =
3128 166 : subCommandHillshade->add_argument("-s")
3129 332 : .metavar("<scale>")
3130 166 : .scan<'g', double>()
3131 166 : .store_into(psOptions->globalScale)
3132 166 : .help(_("Ratio of vertical units to horizontal units."));
3133 :
3134 166 : subCommandHillshade->add_hidden_alias_for(scaleHillshadeArg, "--s");
3135 166 : subCommandHillshade->add_hidden_alias_for(scaleHillshadeArg, "-scale");
3136 166 : subCommandHillshade->add_hidden_alias_for(scaleHillshadeArg, "--scale");
3137 :
3138 : auto &xscaleHillshadeArg =
3139 166 : subCommandHillshade->add_argument("-xscale")
3140 332 : .metavar("<scale>")
3141 166 : .scan<'g', double>()
3142 166 : .store_into(psOptions->xscale)
3143 166 : .help(_("Ratio of vertical units to horizontal X axis units."));
3144 166 : subCommandHillshade->add_hidden_alias_for(xscaleHillshadeArg, "--xscale");
3145 :
3146 : auto &yscaleHillshadeArg =
3147 166 : subCommandHillshade->add_argument("-yscale")
3148 332 : .metavar("<scale>")
3149 166 : .scan<'g', double>()
3150 166 : .store_into(psOptions->yscale)
3151 166 : .help(_("Ratio of vertical units to horizontal Y axis units."));
3152 166 : subCommandHillshade->add_hidden_alias_for(yscaleHillshadeArg, "--yscale");
3153 :
3154 : auto &azimuthHillshadeArg =
3155 166 : subCommandHillshade->add_argument("-az")
3156 332 : .metavar("<azimuth>")
3157 166 : .scan<'g', double>()
3158 166 : .store_into(psOptions->az)
3159 166 : .help(_("Azimuth of the light, in degrees."));
3160 :
3161 166 : subCommandHillshade->add_hidden_alias_for(azimuthHillshadeArg, "--az");
3162 166 : subCommandHillshade->add_hidden_alias_for(azimuthHillshadeArg, "-azimuth");
3163 166 : subCommandHillshade->add_hidden_alias_for(azimuthHillshadeArg, "--azimuth");
3164 :
3165 : auto &altitudeHillshadeArg =
3166 166 : subCommandHillshade->add_argument("-alt")
3167 332 : .metavar("<altitude>")
3168 166 : .scan<'g', double>()
3169 166 : .store_into(psOptions->alt)
3170 166 : .help(_("Altitude of the light, in degrees."));
3171 :
3172 166 : subCommandHillshade->add_hidden_alias_for(altitudeHillshadeArg, "--alt");
3173 : subCommandHillshade->add_hidden_alias_for(altitudeHillshadeArg,
3174 166 : "--altitude");
3175 : subCommandHillshade->add_hidden_alias_for(altitudeHillshadeArg,
3176 166 : "-altitude");
3177 :
3178 166 : auto &shadingGroup = subCommandHillshade->add_mutually_exclusive_group();
3179 :
3180 166 : auto &combinedHillshadeArg = shadingGroup.add_argument("-combined")
3181 166 : .flag()
3182 166 : .store_into(psOptions->bCombined)
3183 166 : .help(_("Use combined shading."));
3184 :
3185 : subCommandHillshade->add_hidden_alias_for(combinedHillshadeArg,
3186 166 : "--combined");
3187 :
3188 : auto &multidirectionalHillshadeArg =
3189 166 : shadingGroup.add_argument("-multidirectional")
3190 166 : .flag()
3191 166 : .store_into(psOptions->bMultiDirectional)
3192 166 : .help(_("Use multidirectional shading."));
3193 :
3194 : subCommandHillshade->add_hidden_alias_for(multidirectionalHillshadeArg,
3195 166 : "--multidirectional");
3196 :
3197 : auto &igorHillshadeArg =
3198 166 : shadingGroup.add_argument("-igor")
3199 166 : .flag()
3200 166 : .store_into(psOptions->bIgor)
3201 : .help(_("Shading which tries to minimize effects on other map "
3202 166 : "features beneath."));
3203 :
3204 166 : subCommandHillshade->add_hidden_alias_for(igorHillshadeArg, "--igor");
3205 :
3206 166 : addCommonOptions(subCommandHillshade);
3207 :
3208 : // Slope
3209 :
3210 166 : auto subCommandSlope = argParser->add_subparser(
3211 : "slope", /* bForBinary=*/psOptionsForBinary != nullptr);
3212 :
3213 166 : subCommandSlope->add_description(_("Compute slope."));
3214 :
3215 166 : if (psOptionsForBinary)
3216 : {
3217 22 : subCommandSlope->add_argument("input_dem")
3218 22 : .store_into(psOptionsForBinary->osSrcFilename)
3219 22 : .help(_("The input DEM raster to be processed."));
3220 :
3221 22 : subCommandSlope->add_argument("output_slope_map")
3222 22 : .store_into(psOptionsForBinary->osDstFilename)
3223 22 : .help(_("The output raster to be produced."));
3224 : }
3225 :
3226 166 : subCommandSlope->add_argument("-alg")
3227 332 : .metavar("Horn|ZevenbergenThorne")
3228 : .action(
3229 9 : [psOptions](const std::string &s)
3230 : {
3231 8 : if (EQUAL(s.c_str(), "ZevenbergenThorne"))
3232 : {
3233 0 : psOptions->bGradientAlgSpecified = true;
3234 0 : psOptions->eGradientAlg = GradientAlg::ZEVENBERGEN_THORNE;
3235 : }
3236 8 : else if (EQUAL(s.c_str(), "Horn"))
3237 : {
3238 1 : psOptions->bGradientAlgSpecified = true;
3239 1 : psOptions->eGradientAlg = GradientAlg::HORN;
3240 : }
3241 : else
3242 : {
3243 : throw std::invalid_argument(
3244 7 : CPLSPrintf("Invalid value for -alg: %s.", s.c_str()));
3245 : }
3246 167 : })
3247 166 : .help(_("Choose between ZevenbergenThorne or Horn algorithms."));
3248 :
3249 : subCommandSlope->add_inverted_logic_flag(
3250 : "-p", &psOptions->bSlopeFormatUseDegrees,
3251 166 : _("Express slope as a percentage."));
3252 :
3253 166 : subCommandSlope->add_argument("-s")
3254 332 : .metavar("<scale>")
3255 166 : .scan<'g', double>()
3256 166 : .store_into(psOptions->globalScale)
3257 166 : .help(_("Ratio of vertical units to horizontal."));
3258 :
3259 : auto &xscaleSlopeArg =
3260 166 : subCommandSlope->add_argument("-xscale")
3261 332 : .metavar("<scale>")
3262 166 : .scan<'g', double>()
3263 166 : .store_into(psOptions->xscale)
3264 166 : .help(_("Ratio of vertical units to horizontal X axis units."));
3265 166 : subCommandSlope->add_hidden_alias_for(xscaleSlopeArg, "--xscale");
3266 :
3267 : auto &yscaleSlopeArg =
3268 166 : subCommandSlope->add_argument("-yscale")
3269 332 : .metavar("<scale>")
3270 166 : .scan<'g', double>()
3271 166 : .store_into(psOptions->yscale)
3272 166 : .help(_("Ratio of vertical units to horizontal Y axis units."));
3273 166 : subCommandSlope->add_hidden_alias_for(yscaleSlopeArg, "--yscale");
3274 :
3275 166 : addCommonOptions(subCommandSlope);
3276 :
3277 : // Aspect
3278 :
3279 166 : auto subCommandAspect = argParser->add_subparser(
3280 : "aspect", /* bForBinary=*/psOptionsForBinary != nullptr);
3281 :
3282 166 : subCommandAspect->add_description(_("Compute aspect."));
3283 :
3284 166 : if (psOptionsForBinary)
3285 : {
3286 22 : subCommandAspect->add_argument("input_dem")
3287 22 : .store_into(psOptionsForBinary->osSrcFilename)
3288 22 : .help(_("The input DEM raster to be processed."));
3289 :
3290 22 : subCommandAspect->add_argument("output_aspect_map")
3291 22 : .store_into(psOptionsForBinary->osDstFilename)
3292 22 : .help(_("The output raster to be produced."));
3293 : }
3294 :
3295 166 : subCommandAspect->add_argument("-alg")
3296 332 : .metavar("Horn|ZevenbergenThorne")
3297 : .action(
3298 11 : [psOptions](const std::string &s)
3299 : {
3300 9 : if (EQUAL(s.c_str(), "ZevenbergenThorne"))
3301 : {
3302 0 : psOptions->bGradientAlgSpecified = true;
3303 0 : psOptions->eGradientAlg = GradientAlg::ZEVENBERGEN_THORNE;
3304 : }
3305 9 : else if (EQUAL(s.c_str(), "Horn"))
3306 : {
3307 2 : psOptions->bGradientAlgSpecified = true;
3308 2 : psOptions->eGradientAlg = GradientAlg::HORN;
3309 : }
3310 : else
3311 : {
3312 : throw std::invalid_argument(
3313 7 : CPLSPrintf("Invalid value for -alg: %s.", s.c_str()));
3314 : }
3315 168 : })
3316 166 : .help(_("Choose between ZevenbergenThorne or Horn algorithms."));
3317 :
3318 : subCommandAspect->add_inverted_logic_flag(
3319 : "-trigonometric", &psOptions->bAngleAsAzimuth,
3320 166 : _("Express aspect in trigonometric format."));
3321 :
3322 166 : subCommandAspect->add_argument("-zero_for_flat")
3323 166 : .flag()
3324 166 : .store_into(psOptions->bZeroForFlat)
3325 166 : .help(_("Return 0 for flat areas with slope=0, instead of -9999."));
3326 :
3327 166 : addCommonOptions(subCommandAspect);
3328 :
3329 : // Color-relief
3330 :
3331 166 : auto subCommandColorRelief = argParser->add_subparser(
3332 : "color-relief", /* bForBinary=*/psOptionsForBinary != nullptr);
3333 :
3334 : subCommandColorRelief->add_description(
3335 : _("Color relief computed from the elevation and a text-based color "
3336 166 : "configuration file."));
3337 :
3338 166 : if (psOptionsForBinary)
3339 : {
3340 22 : subCommandColorRelief->add_argument("input_dem")
3341 22 : .store_into(psOptionsForBinary->osSrcFilename)
3342 22 : .help(_("The input DEM raster to be processed."));
3343 :
3344 22 : subCommandColorRelief->add_argument("color_text_file")
3345 22 : .store_into(psOptionsForBinary->osColorFilename)
3346 22 : .help(_("Text-based color configuration file."));
3347 :
3348 22 : subCommandColorRelief->add_argument("output_color_relief_map")
3349 22 : .store_into(psOptionsForBinary->osDstFilename)
3350 22 : .help(_("The output raster to be produced."));
3351 : }
3352 :
3353 166 : subCommandColorRelief->add_argument("-alpha")
3354 166 : .flag()
3355 166 : .store_into(psOptions->bAddAlpha)
3356 166 : .help(_("Add an alpha channel to the output raster."));
3357 :
3358 166 : subCommandColorRelief->add_argument("-exact_color_entry")
3359 166 : .nargs(0)
3360 : .action(
3361 7 : [psOptions](const auto &)
3362 166 : { psOptions->eColorSelectionMode = COLOR_SELECTION_EXACT_ENTRY; })
3363 : .help(
3364 166 : _("Use strict matching when searching in the configuration file."));
3365 :
3366 166 : subCommandColorRelief->add_argument("-nearest_color_entry")
3367 166 : .nargs(0)
3368 : .action(
3369 9 : [psOptions](const auto &)
3370 166 : { psOptions->eColorSelectionMode = COLOR_SELECTION_NEAREST_ENTRY; })
3371 : .help(_("Use the RGBA corresponding to the closest entry in the "
3372 166 : "configuration file."));
3373 :
3374 166 : addCommonOptions(subCommandColorRelief);
3375 :
3376 : // TRI
3377 :
3378 166 : auto subCommandTRI = argParser->add_subparser(
3379 : "TRI", /* bForBinary=*/psOptionsForBinary != nullptr);
3380 :
3381 166 : subCommandTRI->add_description(_("Compute the Terrain Ruggedness Index."));
3382 :
3383 166 : if (psOptionsForBinary)
3384 : {
3385 :
3386 22 : subCommandTRI->add_argument("input_dem")
3387 22 : .store_into(psOptionsForBinary->osSrcFilename)
3388 22 : .help(_("The input DEM raster to be processed."));
3389 :
3390 22 : subCommandTRI->add_argument("output_TRI_map")
3391 22 : .store_into(psOptionsForBinary->osDstFilename)
3392 22 : .help(_("The output raster to be produced."));
3393 : }
3394 :
3395 166 : subCommandTRI->add_argument("-alg")
3396 332 : .metavar("Wilson|Riley")
3397 : .action(
3398 14 : [psOptions](const std::string &s)
3399 : {
3400 7 : if (EQUAL(s.c_str(), "Wilson"))
3401 : {
3402 2 : psOptions->bTRIAlgSpecified = true;
3403 2 : psOptions->eTRIAlg = TRIAlg::WILSON;
3404 : }
3405 5 : else if (EQUAL(s.c_str(), "Riley"))
3406 : {
3407 5 : psOptions->bTRIAlgSpecified = true;
3408 5 : psOptions->eTRIAlg = TRIAlg::RILEY;
3409 : }
3410 : else
3411 : {
3412 : throw std::invalid_argument(
3413 0 : CPLSPrintf("Invalid value for -alg: %s.", s.c_str()));
3414 : }
3415 173 : })
3416 166 : .help(_("Choose between Wilson or Riley algorithms."));
3417 :
3418 166 : addCommonOptions(subCommandTRI);
3419 :
3420 : // TPI
3421 :
3422 166 : auto subCommandTPI = argParser->add_subparser(
3423 : "TPI", /* bForBinary=*/psOptionsForBinary != nullptr);
3424 :
3425 : subCommandTPI->add_description(
3426 166 : _("Compute the Topographic Position Index."));
3427 :
3428 166 : if (psOptionsForBinary)
3429 : {
3430 22 : subCommandTPI->add_argument("input_dem")
3431 22 : .store_into(psOptionsForBinary->osSrcFilename)
3432 22 : .help(_("The input DEM raster to be processed."));
3433 :
3434 22 : subCommandTPI->add_argument("output_TPI_map")
3435 22 : .store_into(psOptionsForBinary->osDstFilename)
3436 22 : .help(_("The output raster to be produced."));
3437 : }
3438 :
3439 166 : addCommonOptions(subCommandTPI);
3440 :
3441 : // Roughness
3442 :
3443 166 : auto subCommandRoughness = argParser->add_subparser(
3444 : "roughness", /* bForBinary=*/psOptionsForBinary != nullptr);
3445 :
3446 : subCommandRoughness->add_description(
3447 166 : _("Compute the roughness of the DEM."));
3448 :
3449 166 : if (psOptionsForBinary)
3450 : {
3451 22 : subCommandRoughness->add_argument("input_dem")
3452 22 : .store_into(psOptionsForBinary->osSrcFilename)
3453 22 : .help(_("The input DEM raster to be processed."));
3454 :
3455 22 : subCommandRoughness->add_argument("output_roughness_map")
3456 22 : .store_into(psOptionsForBinary->osDstFilename)
3457 22 : .help(_("The output raster to be produced."));
3458 : }
3459 :
3460 166 : addCommonOptions(subCommandRoughness);
3461 :
3462 332 : return argParser;
3463 : }
3464 :
3465 : /************************************************************************/
3466 : /* GDALDEMAppGetParserUsage() */
3467 : /************************************************************************/
3468 :
3469 1 : std::string GDALDEMAppGetParserUsage(const std::string &osProcessingMode)
3470 : {
3471 : try
3472 : {
3473 2 : GDALDEMProcessingOptions sOptions;
3474 2 : GDALDEMProcessingOptionsForBinary sOptionsForBinary;
3475 : auto argParser =
3476 2 : GDALDEMAppOptionsGetParser(&sOptions, &sOptionsForBinary);
3477 1 : if (!osProcessingMode.empty())
3478 : {
3479 0 : const auto subParser = argParser->get_subparser(osProcessingMode);
3480 0 : if (subParser)
3481 : {
3482 0 : return subParser->usage();
3483 : }
3484 : else
3485 : {
3486 0 : CPLError(CE_Failure, CPLE_AppDefined,
3487 : "Invalid processing mode: %s",
3488 : osProcessingMode.c_str());
3489 : }
3490 : }
3491 1 : return argParser->usage();
3492 : }
3493 0 : catch (const std::exception &err)
3494 : {
3495 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
3496 0 : err.what());
3497 0 : return std::string();
3498 : }
3499 : }
3500 :
3501 : /************************************************************************/
3502 : /* GDALDEMProcessing() */
3503 : /************************************************************************/
3504 :
3505 : /**
3506 : * Apply a DEM processing.
3507 : *
3508 : * This is the equivalent of the <a href="/programs/gdaldem.html">gdaldem</a>
3509 : * utility.
3510 : *
3511 : * GDALDEMProcessingOptions* must be allocated and freed with
3512 : * GDALDEMProcessingOptionsNew() and GDALDEMProcessingOptionsFree()
3513 : * respectively.
3514 : *
3515 : * @param pszDest the destination dataset path.
3516 : * @param hSrcDataset the source dataset handle.
3517 : * @param pszProcessing the processing to apply (one of "hillshade", "slope",
3518 : * "aspect", "color-relief", "TRI", "TPI", "Roughness")
3519 : * @param pszColorFilename color file (mandatory for "color-relief" processing,
3520 : * should be NULL otherwise)
3521 : * @param psOptionsIn the options struct returned by
3522 : * GDALDEMProcessingOptionsNew() or NULL.
3523 : * @param pbUsageError pointer to a integer output variable to store if any
3524 : * usage error has occurred or NULL.
3525 : * @return the output dataset (new dataset that must be closed using
3526 : * GDALClose()) or NULL in case of error.
3527 : *
3528 : * @since GDAL 2.1
3529 : */
3530 :
3531 161 : GDALDatasetH GDALDEMProcessing(const char *pszDest, GDALDatasetH hSrcDataset,
3532 : const char *pszProcessing,
3533 : const char *pszColorFilename,
3534 : const GDALDEMProcessingOptions *psOptionsIn,
3535 : int *pbUsageError)
3536 : {
3537 161 : if (hSrcDataset == nullptr)
3538 : {
3539 0 : CPLError(CE_Failure, CPLE_AppDefined, "No source dataset specified.");
3540 :
3541 0 : if (pbUsageError)
3542 0 : *pbUsageError = TRUE;
3543 0 : return nullptr;
3544 : }
3545 161 : if (pszDest == nullptr)
3546 : {
3547 0 : CPLError(CE_Failure, CPLE_AppDefined, "No target dataset specified.");
3548 :
3549 0 : if (pbUsageError)
3550 0 : *pbUsageError = TRUE;
3551 0 : return nullptr;
3552 : }
3553 161 : if (pszProcessing == nullptr)
3554 : {
3555 0 : CPLError(CE_Failure, CPLE_AppDefined, "No target dataset specified.");
3556 :
3557 0 : if (pbUsageError)
3558 0 : *pbUsageError = TRUE;
3559 0 : return nullptr;
3560 : }
3561 :
3562 161 : Algorithm eUtilityMode = GetAlgorithm(pszProcessing);
3563 161 : if (eUtilityMode == INVALID)
3564 : {
3565 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Invalid processing mode: %s",
3566 : pszProcessing);
3567 0 : if (pbUsageError)
3568 0 : *pbUsageError = TRUE;
3569 0 : return nullptr;
3570 : }
3571 :
3572 161 : if (eUtilityMode == COLOR_RELIEF &&
3573 41 : (pszColorFilename == nullptr || strlen(pszColorFilename) == 0))
3574 : {
3575 0 : CPLError(CE_Failure, CPLE_AppDefined, "pszColorFilename == NULL.");
3576 :
3577 0 : if (pbUsageError)
3578 0 : *pbUsageError = TRUE;
3579 0 : return nullptr;
3580 : }
3581 161 : else if (eUtilityMode != COLOR_RELIEF && pszColorFilename != nullptr &&
3582 9 : strlen(pszColorFilename) > 0)
3583 : {
3584 0 : CPLError(CE_Failure, CPLE_AppDefined, "pszColorFilename != NULL.");
3585 :
3586 0 : if (pbUsageError)
3587 0 : *pbUsageError = TRUE;
3588 0 : return nullptr;
3589 : }
3590 :
3591 161 : if (psOptionsIn && psOptionsIn->bCombined && eUtilityMode != HILL_SHADE)
3592 : {
3593 0 : CPLError(CE_Failure, CPLE_NotSupported,
3594 : "-combined can only be used with hillshade");
3595 :
3596 0 : if (pbUsageError)
3597 0 : *pbUsageError = TRUE;
3598 0 : return nullptr;
3599 : }
3600 :
3601 161 : if (psOptionsIn && psOptionsIn->bIgor && eUtilityMode != HILL_SHADE)
3602 : {
3603 0 : CPLError(CE_Failure, CPLE_NotSupported,
3604 : "-igor can only be used with hillshade");
3605 :
3606 0 : if (pbUsageError)
3607 0 : *pbUsageError = TRUE;
3608 0 : return nullptr;
3609 : }
3610 :
3611 161 : if (psOptionsIn && psOptionsIn->bMultiDirectional &&
3612 : eUtilityMode != HILL_SHADE)
3613 : {
3614 0 : CPLError(CE_Failure, CPLE_NotSupported,
3615 : "-multidirectional can only be used with hillshade");
3616 :
3617 0 : if (pbUsageError)
3618 0 : *pbUsageError = TRUE;
3619 0 : return nullptr;
3620 : }
3621 :
3622 161 : std::unique_ptr<GDALDEMProcessingOptions> psOptionsToFree;
3623 161 : if (psOptionsIn)
3624 : {
3625 : psOptionsToFree =
3626 161 : std::make_unique<GDALDEMProcessingOptions>(*psOptionsIn);
3627 : }
3628 : else
3629 : {
3630 0 : psOptionsToFree.reset(GDALDEMProcessingOptionsNew(nullptr, nullptr));
3631 : }
3632 :
3633 161 : GDALDEMProcessingOptions *psOptions = psOptionsToFree.get();
3634 :
3635 161 : const int nXSize = GDALGetRasterXSize(hSrcDataset);
3636 161 : const int nYSize = GDALGetRasterYSize(hSrcDataset);
3637 :
3638 322 : if (psOptions->nBand <= 0 ||
3639 161 : psOptions->nBand > GDALGetRasterCount(hSrcDataset))
3640 : {
3641 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Unable to fetch band #%d",
3642 : psOptions->nBand);
3643 :
3644 0 : return nullptr;
3645 : }
3646 :
3647 161 : if (std::isnan(psOptions->xscale))
3648 : {
3649 124 : psOptions->xscale = 1;
3650 124 : psOptions->yscale = 1;
3651 :
3652 124 : double zunit = 1;
3653 :
3654 124 : auto poSrcDS = GDALDataset::FromHandle(hSrcDataset);
3655 :
3656 : const char *pszUnits =
3657 124 : poSrcDS->GetRasterBand(psOptions->nBand)->GetUnitType();
3658 124 : if (EQUAL(pszUnits, "m") || EQUAL(pszUnits, "metre") ||
3659 51 : EQUAL(pszUnits, "meter"))
3660 : {
3661 : }
3662 51 : else if (EQUAL(pszUnits, "ft") || EQUAL(pszUnits, "foot") ||
3663 49 : EQUAL(pszUnits, "foot (International)") ||
3664 49 : EQUAL(pszUnits, "feet"))
3665 : {
3666 2 : zunit = CPLAtof(SRS_UL_FOOT_CONV);
3667 : }
3668 49 : else if (EQUAL(pszUnits, "us-ft") || EQUAL(pszUnits, "Foot_US") ||
3669 47 : EQUAL(pszUnits, "US survey foot"))
3670 : {
3671 2 : zunit = CPLAtof(SRS_UL_US_FOOT_CONV);
3672 : }
3673 47 : else if (!EQUAL(pszUnits, ""))
3674 : {
3675 2 : CPLError(CE_Warning, CPLE_AppDefined,
3676 : "Unknown band unit '%s'. Assuming metre", pszUnits);
3677 : }
3678 :
3679 124 : const auto poSrcSRS = poSrcDS->GetSpatialRef();
3680 124 : if (poSrcSRS && poSrcSRS->IsGeographic())
3681 : {
3682 : double adfGT[6];
3683 87 : if (poSrcDS->GetGeoTransform(adfGT) == CE_None)
3684 : {
3685 87 : const double dfAngUnits = poSrcSRS->GetAngularUnits();
3686 : // Rough conversion of angular units to linear units.
3687 87 : psOptions->yscale =
3688 87 : dfAngUnits * poSrcSRS->GetSemiMajor() / zunit;
3689 : // Take the center latitude to compute the xscale.
3690 87 : const double dfMeanLat =
3691 87 : (adfGT[3] + nYSize * adfGT[5] / 2) * dfAngUnits;
3692 87 : if (std::fabs(dfMeanLat) / M_PI * 180 > 80)
3693 : {
3694 0 : CPLError(
3695 : CE_Warning, CPLE_AppDefined,
3696 : "Automatic computation of xscale at high latitudes may "
3697 : "lead to incorrect results. The source dataset should "
3698 : "likely be reprojected to a polar projection");
3699 : }
3700 87 : psOptions->xscale = psOptions->yscale * cos(dfMeanLat);
3701 : }
3702 : }
3703 37 : else if (poSrcSRS && poSrcSRS->IsProjected())
3704 : {
3705 7 : psOptions->xscale = poSrcSRS->GetLinearUnits() / zunit;
3706 7 : psOptions->yscale = psOptions->xscale;
3707 : }
3708 124 : CPLDebug("GDAL", "Using xscale=%f and yscale=%f", psOptions->xscale,
3709 : psOptions->yscale);
3710 : }
3711 :
3712 161 : if (psOptions->bGradientAlgSpecified &&
3713 26 : !(eUtilityMode == HILL_SHADE || eUtilityMode == SLOPE ||
3714 : eUtilityMode == ASPECT))
3715 : {
3716 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3717 : "This value of -alg is only valid for hillshade, slope or "
3718 : "aspect processing");
3719 :
3720 0 : return nullptr;
3721 : }
3722 :
3723 161 : if (psOptions->bTRIAlgSpecified && !(eUtilityMode == TRI))
3724 : {
3725 0 : CPLError(CE_Failure, CPLE_IllegalArg,
3726 : "This value of -alg is only valid for TRI processing");
3727 :
3728 0 : return nullptr;
3729 : }
3730 :
3731 161 : double adfGeoTransform[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3732 :
3733 161 : GDALRasterBandH hSrcBand = GDALGetRasterBand(hSrcDataset, psOptions->nBand);
3734 :
3735 161 : GDALGetGeoTransform(hSrcDataset, adfGeoTransform);
3736 :
3737 322 : CPLString osFormat;
3738 161 : if (psOptions->osFormat.empty())
3739 : {
3740 14 : osFormat = GetOutputDriverForRaster(pszDest);
3741 14 : if (osFormat.empty())
3742 : {
3743 0 : return nullptr;
3744 : }
3745 : }
3746 : else
3747 : {
3748 147 : osFormat = psOptions->osFormat;
3749 : }
3750 :
3751 161 : GDALDriverH hDriver = nullptr;
3752 161 : if (!EQUAL(osFormat.c_str(), "stream"))
3753 : {
3754 120 : hDriver = GDALGetDriverByName(osFormat);
3755 240 : if (hDriver == nullptr ||
3756 120 : (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) ==
3757 5 : nullptr &&
3758 5 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, nullptr) ==
3759 : nullptr))
3760 : {
3761 0 : CPLError(CE_Failure, CPLE_AppDefined,
3762 : "Output driver `%s' does not support writing.",
3763 : osFormat.c_str());
3764 0 : fprintf(stderr, "The following format drivers are enabled\n"
3765 : "and support writing:\n");
3766 :
3767 0 : for (int iDr = 0; iDr < GDALGetDriverCount(); iDr++)
3768 : {
3769 0 : hDriver = GDALGetDriver(iDr);
3770 :
3771 0 : if (GDALGetMetadataItem(hDriver, GDAL_DCAP_RASTER, nullptr) !=
3772 0 : nullptr &&
3773 0 : (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) !=
3774 0 : nullptr ||
3775 0 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY,
3776 : nullptr) != nullptr))
3777 : {
3778 0 : fprintf(stderr, " %s: %s\n",
3779 : GDALGetDriverShortName(hDriver),
3780 : GDALGetDriverLongName(hDriver));
3781 : }
3782 : }
3783 :
3784 0 : return nullptr;
3785 : }
3786 : }
3787 :
3788 161 : double dfDstNoDataValue = 0.0;
3789 161 : bool bDstHasNoData = false;
3790 161 : std::unique_ptr<AlgorithmParameters> pData;
3791 161 : GDALGeneric3x3ProcessingAlg<float>::type pfnAlgFloat = nullptr;
3792 161 : GDALGeneric3x3ProcessingAlg<GInt32>::type pfnAlgInt32 = nullptr;
3793 : GDALGeneric3x3ProcessingAlg_multisample<float>::type
3794 161 : pfnAlgFloat_multisample = nullptr;
3795 : GDALGeneric3x3ProcessingAlg_multisample<GInt32>::type
3796 161 : pfnAlgInt32_multisample = nullptr;
3797 :
3798 161 : if (eUtilityMode == HILL_SHADE && psOptions->bMultiDirectional)
3799 : {
3800 3 : dfDstNoDataValue = 0;
3801 3 : bDstHasNoData = true;
3802 6 : pData = GDALCreateHillshadeMultiDirectionalData(
3803 : adfGeoTransform, psOptions->z, psOptions->xscale, psOptions->yscale,
3804 3 : psOptions->alt, psOptions->eGradientAlg);
3805 3 : if (psOptions->eGradientAlg == GradientAlg::ZEVENBERGEN_THORNE)
3806 : {
3807 1 : pfnAlgFloat = GDALHillshadeMultiDirectionalAlg<
3808 : float, GradientAlg::ZEVENBERGEN_THORNE>;
3809 1 : pfnAlgInt32 = GDALHillshadeMultiDirectionalAlg<
3810 : GInt32, GradientAlg::ZEVENBERGEN_THORNE>;
3811 : }
3812 : else
3813 : {
3814 2 : pfnAlgFloat =
3815 : GDALHillshadeMultiDirectionalAlg<float, GradientAlg::HORN>;
3816 2 : pfnAlgInt32 =
3817 : GDALHillshadeMultiDirectionalAlg<GInt32, GradientAlg::HORN>;
3818 : }
3819 : }
3820 158 : else if (eUtilityMode == HILL_SHADE)
3821 : {
3822 69 : dfDstNoDataValue = 0;
3823 69 : bDstHasNoData = true;
3824 138 : pData = GDALCreateHillshadeData(
3825 : adfGeoTransform, psOptions->z, psOptions->xscale, psOptions->yscale,
3826 69 : psOptions->alt, psOptions->az, psOptions->eGradientAlg);
3827 69 : if (psOptions->eGradientAlg == GradientAlg::ZEVENBERGEN_THORNE)
3828 : {
3829 10 : if (psOptions->bCombined)
3830 : {
3831 4 : pfnAlgFloat =
3832 : GDALHillshadeCombinedAlg<float,
3833 : GradientAlg::ZEVENBERGEN_THORNE>;
3834 4 : pfnAlgInt32 =
3835 : GDALHillshadeCombinedAlg<GInt32,
3836 : GradientAlg::ZEVENBERGEN_THORNE>;
3837 : }
3838 6 : else if (psOptions->bIgor)
3839 : {
3840 1 : pfnAlgFloat =
3841 : GDALHillshadeIgorAlg<float,
3842 : GradientAlg::ZEVENBERGEN_THORNE>;
3843 1 : pfnAlgInt32 =
3844 : GDALHillshadeIgorAlg<GInt32,
3845 : GradientAlg::ZEVENBERGEN_THORNE>;
3846 : }
3847 : else
3848 : {
3849 5 : pfnAlgFloat =
3850 : GDALHillshadeAlg<float, GradientAlg::ZEVENBERGEN_THORNE>;
3851 5 : pfnAlgInt32 =
3852 : GDALHillshadeAlg<GInt32, GradientAlg::ZEVENBERGEN_THORNE>;
3853 : }
3854 : }
3855 : else
3856 : {
3857 59 : if (psOptions->bCombined)
3858 : {
3859 6 : pfnAlgFloat =
3860 : GDALHillshadeCombinedAlg<float, GradientAlg::HORN>;
3861 6 : pfnAlgInt32 =
3862 : GDALHillshadeCombinedAlg<GInt32, GradientAlg::HORN>;
3863 : }
3864 53 : else if (psOptions->bIgor)
3865 : {
3866 2 : pfnAlgFloat = GDALHillshadeIgorAlg<float, GradientAlg::HORN>;
3867 2 : pfnAlgInt32 = GDALHillshadeIgorAlg<GInt32, GradientAlg::HORN>;
3868 : }
3869 : else
3870 : {
3871 51 : if (adfGeoTransform[1] == -adfGeoTransform[5] &&
3872 51 : psOptions->xscale == psOptions->yscale)
3873 : {
3874 34 : pfnAlgFloat = GDALHillshadeAlg_same_res<float>;
3875 34 : pfnAlgInt32 = GDALHillshadeAlg_same_res<GInt32>;
3876 : #ifdef HAVE_16_SSE_REG
3877 34 : pfnAlgFloat_multisample =
3878 : GDALHillshadeAlg_same_res_multisample<float,
3879 : XMMReg4Float>;
3880 34 : pfnAlgInt32_multisample =
3881 : GDALHillshadeAlg_same_res_multisample<GInt32,
3882 : XMMReg4Int>;
3883 : #endif
3884 : }
3885 : else
3886 : {
3887 17 : pfnAlgFloat = GDALHillshadeAlg<float, GradientAlg::HORN>;
3888 17 : pfnAlgInt32 = GDALHillshadeAlg<GInt32, GradientAlg::HORN>;
3889 : }
3890 : }
3891 : }
3892 : }
3893 89 : else if (eUtilityMode == SLOPE)
3894 : {
3895 16 : dfDstNoDataValue = -9999;
3896 16 : bDstHasNoData = true;
3897 :
3898 16 : pData = GDALCreateSlopeData(adfGeoTransform, psOptions->xscale,
3899 : psOptions->yscale,
3900 16 : psOptions->bSlopeFormatUseDegrees);
3901 16 : if (psOptions->eGradientAlg == GradientAlg::ZEVENBERGEN_THORNE)
3902 : {
3903 6 : pfnAlgFloat = GDALSlopeZevenbergenThorneAlg<float>;
3904 6 : pfnAlgInt32 = GDALSlopeZevenbergenThorneAlg<GInt32>;
3905 : }
3906 : else
3907 : {
3908 10 : pfnAlgFloat = GDALSlopeHornAlg<float>;
3909 10 : pfnAlgInt32 = GDALSlopeHornAlg<GInt32>;
3910 : }
3911 : }
3912 :
3913 73 : else if (eUtilityMode == ASPECT)
3914 : {
3915 12 : if (!psOptions->bZeroForFlat)
3916 : {
3917 11 : dfDstNoDataValue = -9999;
3918 11 : bDstHasNoData = true;
3919 : }
3920 :
3921 12 : pData = GDALCreateAspectData(psOptions->bAngleAsAzimuth);
3922 12 : if (psOptions->eGradientAlg == GradientAlg::ZEVENBERGEN_THORNE)
3923 : {
3924 3 : pfnAlgFloat = GDALAspectZevenbergenThorneAlg<float>;
3925 3 : pfnAlgInt32 = GDALAspectZevenbergenThorneAlg<GInt32>;
3926 : }
3927 : else
3928 : {
3929 9 : pfnAlgFloat = GDALAspectAlg<float>;
3930 9 : pfnAlgInt32 = GDALAspectAlg<GInt32>;
3931 : }
3932 : }
3933 61 : else if (eUtilityMode == TRI)
3934 : {
3935 8 : dfDstNoDataValue = -9999;
3936 8 : bDstHasNoData = true;
3937 8 : if (psOptions->eTRIAlg == TRIAlg::WILSON)
3938 : {
3939 2 : pfnAlgFloat = GDALTRIAlgWilson<float>;
3940 2 : pfnAlgInt32 = GDALTRIAlgWilson<GInt32>;
3941 : }
3942 : else
3943 : {
3944 6 : pfnAlgFloat = GDALTRIAlgRiley<float>;
3945 6 : pfnAlgInt32 = GDALTRIAlgRiley<GInt32>;
3946 : }
3947 : }
3948 53 : else if (eUtilityMode == TPI)
3949 : {
3950 5 : dfDstNoDataValue = -9999;
3951 5 : bDstHasNoData = true;
3952 5 : pfnAlgFloat = GDALTPIAlg<float>;
3953 5 : pfnAlgInt32 = GDALTPIAlg<GInt32>;
3954 : }
3955 48 : else if (eUtilityMode == ROUGHNESS)
3956 : {
3957 7 : dfDstNoDataValue = -9999;
3958 7 : bDstHasNoData = true;
3959 7 : pfnAlgFloat = GDALRoughnessAlg<float>;
3960 7 : pfnAlgInt32 = GDALRoughnessAlg<GInt32>;
3961 : }
3962 :
3963 161 : const GDALDataType eDstDataType =
3964 89 : (eUtilityMode == HILL_SHADE || eUtilityMode == COLOR_RELIEF)
3965 250 : ? GDT_Byte
3966 : : GDT_Float32;
3967 :
3968 161 : if (EQUAL(osFormat, "VRT"))
3969 : {
3970 14 : if (eUtilityMode == COLOR_RELIEF)
3971 : {
3972 14 : const bool bTmpFile = pszDest[0] == 0;
3973 : const std::string osTmpFile =
3974 14 : VSIMemGenerateHiddenFilename("tmp.vrt");
3975 14 : if (bTmpFile)
3976 5 : pszDest = osTmpFile.c_str();
3977 14 : GDALDatasetH hDS = nullptr;
3978 14 : if (GDALGenerateVRTColorRelief(
3979 : pszDest, hSrcDataset, hSrcBand, pszColorFilename,
3980 14 : psOptions->eColorSelectionMode, psOptions->bAddAlpha))
3981 : {
3982 14 : if (bTmpFile)
3983 : {
3984 : const GByte *pabyData =
3985 5 : VSIGetMemFileBuffer(pszDest, nullptr, false);
3986 5 : hDS = GDALOpen(reinterpret_cast<const char *>(pabyData),
3987 : GA_Update);
3988 : }
3989 : else
3990 : {
3991 9 : hDS = GDALOpen(pszDest, GA_Update);
3992 : }
3993 : }
3994 14 : if (bTmpFile)
3995 5 : VSIUnlink(pszDest);
3996 14 : return hDS;
3997 : }
3998 : else
3999 : {
4000 0 : CPLError(CE_Failure, CPLE_NotSupported,
4001 : "VRT driver can only be used with color-relief utility.");
4002 :
4003 0 : return nullptr;
4004 : }
4005 : }
4006 :
4007 : // We might actually want to always go through the intermediate dataset
4008 147 : bool bForceUseIntermediateDataset = false;
4009 :
4010 147 : GDALProgressFunc pfnProgress = psOptions->pfnProgress;
4011 147 : void *pProgressData = psOptions->pProgressData;
4012 :
4013 147 : if (EQUAL(osFormat, "GTiff"))
4014 : {
4015 14 : if (!EQUAL(psOptions->aosCreationOptions.FetchNameValueDef("COMPRESS",
4016 : "NONE"),
4017 16 : "NONE") &&
4018 2 : CPLTestBool(
4019 : psOptions->aosCreationOptions.FetchNameValueDef("TILED", "NO")))
4020 : {
4021 1 : bForceUseIntermediateDataset = true;
4022 : }
4023 13 : else if (strcmp(pszDest, "/vsistdout/") == 0)
4024 : {
4025 0 : bForceUseIntermediateDataset = true;
4026 0 : pfnProgress = GDALDummyProgress;
4027 0 : pProgressData = nullptr;
4028 : }
4029 : #ifdef S_ISFIFO
4030 : else
4031 : {
4032 : VSIStatBufL sStat;
4033 13 : if (VSIStatExL(pszDest, &sStat,
4034 13 : VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG) == 0 &&
4035 0 : S_ISFIFO(sStat.st_mode))
4036 : {
4037 0 : bForceUseIntermediateDataset = true;
4038 : }
4039 : }
4040 : #endif
4041 : }
4042 :
4043 147 : const GDALDataType eSrcDT = GDALGetRasterDataType(hSrcBand);
4044 :
4045 253 : if (hDriver == nullptr ||
4046 106 : (GDALGetMetadataItem(hDriver, GDAL_DCAP_RASTER, nullptr) != nullptr &&
4047 105 : ((bForceUseIntermediateDataset ||
4048 105 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) ==
4049 6 : nullptr) &&
4050 6 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, nullptr) !=
4051 : nullptr)))
4052 : {
4053 47 : GDALDatasetH hIntermediateDataset = nullptr;
4054 :
4055 47 : if (eUtilityMode == COLOR_RELIEF)
4056 : {
4057 : auto poDS = std::make_unique<GDALColorReliefDataset>(
4058 : hSrcDataset, hSrcBand, pszColorFilename,
4059 2 : psOptions->eColorSelectionMode, psOptions->bAddAlpha);
4060 2 : if (!(poDS->InitOK()))
4061 : {
4062 0 : return nullptr;
4063 : }
4064 2 : hIntermediateDataset = GDALDataset::ToHandle(poDS.release());
4065 : }
4066 : else
4067 : {
4068 45 : if (eSrcDT == GDT_Byte || eSrcDT == GDT_Int16 ||
4069 : eSrcDT == GDT_UInt16)
4070 : {
4071 : auto poDS = std::make_unique<GDALGeneric3x3Dataset<GInt32>>(
4072 : hSrcDataset, hSrcBand, eDstDataType, bDstHasNoData,
4073 : dfDstNoDataValue, pfnAlgInt32, pfnAlgInt32_multisample,
4074 44 : std::move(pData), psOptions->bComputeAtEdges, true);
4075 :
4076 44 : if (!(poDS->InitOK()))
4077 : {
4078 0 : return nullptr;
4079 : }
4080 88 : hIntermediateDataset = GDALDataset::ToHandle(poDS.release());
4081 : }
4082 : else
4083 : {
4084 : auto poDS = std::make_unique<GDALGeneric3x3Dataset<float>>(
4085 : hSrcDataset, hSrcBand, eDstDataType, bDstHasNoData,
4086 : dfDstNoDataValue, pfnAlgFloat, pfnAlgFloat_multisample,
4087 1 : std::move(pData), psOptions->bComputeAtEdges, true);
4088 :
4089 1 : if (!(poDS->InitOK()))
4090 : {
4091 0 : return nullptr;
4092 : }
4093 1 : hIntermediateDataset = GDALDataset::ToHandle(poDS.release());
4094 : }
4095 : }
4096 :
4097 47 : if (!hDriver)
4098 : {
4099 41 : return hIntermediateDataset;
4100 : }
4101 :
4102 6 : GDALDatasetH hOutDS = GDALCreateCopy(
4103 : hDriver, pszDest, hIntermediateDataset, TRUE,
4104 6 : psOptions->aosCreationOptions.List(), pfnProgress, pProgressData);
4105 :
4106 6 : GDALClose(hIntermediateDataset);
4107 :
4108 6 : return hOutDS;
4109 : }
4110 :
4111 100 : const int nDstBands =
4112 100 : eUtilityMode == COLOR_RELIEF ? ((psOptions->bAddAlpha) ? 4 : 3) : 1;
4113 :
4114 : GDALDatasetH hDstDataset =
4115 100 : GDALCreate(hDriver, pszDest, nXSize, nYSize, nDstBands, eDstDataType,
4116 100 : psOptions->aosCreationOptions.List());
4117 :
4118 100 : if (hDstDataset == nullptr)
4119 : {
4120 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unable to create dataset %s",
4121 : pszDest);
4122 0 : return nullptr;
4123 : }
4124 :
4125 100 : GDALRasterBandH hDstBand = GDALGetRasterBand(hDstDataset, 1);
4126 :
4127 100 : GDALSetGeoTransform(hDstDataset, adfGeoTransform);
4128 100 : GDALSetProjection(hDstDataset, GDALGetProjectionRef(hSrcDataset));
4129 :
4130 100 : if (eUtilityMode == COLOR_RELIEF)
4131 : {
4132 25 : GDALColorRelief(hSrcBand, GDALGetRasterBand(hDstDataset, 1),
4133 : GDALGetRasterBand(hDstDataset, 2),
4134 : GDALGetRasterBand(hDstDataset, 3),
4135 25 : psOptions->bAddAlpha ? GDALGetRasterBand(hDstDataset, 4)
4136 : : nullptr,
4137 : pszColorFilename, psOptions->eColorSelectionMode,
4138 : pfnProgress, pProgressData);
4139 : }
4140 : else
4141 : {
4142 75 : if (bDstHasNoData)
4143 75 : GDALSetRasterNoDataValue(hDstBand, dfDstNoDataValue);
4144 :
4145 75 : if (eSrcDT == GDT_Byte || eSrcDT == GDT_Int16 || eSrcDT == GDT_UInt16)
4146 : {
4147 63 : GDALGeneric3x3Processing<GInt32>(
4148 : hSrcBand, hDstBand, pfnAlgInt32, pfnAlgInt32_multisample,
4149 63 : std::move(pData), psOptions->bComputeAtEdges, pfnProgress,
4150 : pProgressData);
4151 : }
4152 : else
4153 : {
4154 12 : GDALGeneric3x3Processing<float>(
4155 : hSrcBand, hDstBand, pfnAlgFloat, pfnAlgFloat_multisample,
4156 12 : std::move(pData), psOptions->bComputeAtEdges, pfnProgress,
4157 : pProgressData);
4158 : }
4159 : }
4160 :
4161 100 : return hDstDataset;
4162 : }
4163 :
4164 : /************************************************************************/
4165 : /* GDALDEMProcessingOptionsNew() */
4166 : /************************************************************************/
4167 :
4168 : /**
4169 : * Allocates a GDALDEMProcessingOptions struct.
4170 : *
4171 : * @param papszArgv NULL terminated list of options (potentially including
4172 : * filename and open options too), or NULL. The accepted options are the ones of
4173 : * the <a href="/programs/gdaldem.html">gdaldem</a> utility.
4174 : * @param psOptionsForBinary (output) may be NULL (and should generally be
4175 : * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
4176 : * GDALDEMProcessingOptionsForBinaryNew() prior to
4177 : * this function. Will be filled with potentially present filename, open
4178 : * options,...
4179 : * @return pointer to the allocated GDALDEMProcessingOptions struct. Must be
4180 : * freed with GDALDEMProcessingOptionsFree().
4181 : *
4182 : * @since GDAL 2.1
4183 : */
4184 :
4185 165 : GDALDEMProcessingOptions *GDALDEMProcessingOptionsNew(
4186 : char **papszArgv, GDALDEMProcessingOptionsForBinary *psOptionsForBinary)
4187 : {
4188 :
4189 330 : auto psOptions = std::make_unique<GDALDEMProcessingOptions>();
4190 : /* -------------------------------------------------------------------- */
4191 : /* Handle command line arguments. */
4192 : /* -------------------------------------------------------------------- */
4193 330 : CPLStringList aosArgv;
4194 :
4195 165 : if (papszArgv)
4196 : {
4197 165 : const int nArgc = CSLCount(papszArgv);
4198 1353 : for (int i = 0; i < nArgc; i++)
4199 1188 : aosArgv.AddString(papszArgv[i]);
4200 : }
4201 :
4202 : // Ugly hack: papszArgv may not contain the processing mode if coming from SWIG
4203 165 : const bool bNoProcessingMode(aosArgv.size() > 0 && aosArgv[0][0] == '-');
4204 :
4205 : auto argParser =
4206 330 : GDALDEMAppOptionsGetParser(psOptions.get(), psOptionsForBinary);
4207 :
4208 393 : auto tryHandleArgv = [&](const CPLStringList &args)
4209 : {
4210 393 : argParser->parse_args_without_binary_name(args);
4211 : // Validate the parsed options
4212 :
4213 165 : if (psOptions->nBand <= 0)
4214 : {
4215 0 : throw std::invalid_argument("Invalid value for -b");
4216 : }
4217 :
4218 165 : if (psOptions->z <= 0)
4219 : {
4220 0 : throw std::invalid_argument("Invalid value for -z");
4221 : }
4222 :
4223 165 : if (psOptions->globalScale <= 0)
4224 : {
4225 0 : throw std::invalid_argument("Invalid value for -s");
4226 : }
4227 :
4228 165 : if (psOptions->xscale <= 0)
4229 : {
4230 0 : throw std::invalid_argument("Invalid value for -xscale");
4231 : }
4232 :
4233 165 : if (psOptions->yscale <= 0)
4234 : {
4235 0 : throw std::invalid_argument("Invalid value for -yscale");
4236 : }
4237 :
4238 165 : if (psOptions->alt <= 0)
4239 : {
4240 0 : throw std::invalid_argument("Invalid value for -alt");
4241 : }
4242 :
4243 165 : if (psOptions->bMultiDirectional && argParser->is_used_globally("-az"))
4244 : {
4245 : throw std::invalid_argument(
4246 0 : "-multidirectional and -az cannot be used together");
4247 : }
4248 :
4249 165 : if (psOptions->bIgor && argParser->is_used_globally("-alt"))
4250 : {
4251 : throw std::invalid_argument(
4252 0 : "-igor and -alt cannot be used together");
4253 : }
4254 :
4255 165 : if (psOptionsForBinary && aosArgv.size() > 1)
4256 : {
4257 21 : psOptionsForBinary->osProcessing = aosArgv[0];
4258 : }
4259 165 : };
4260 :
4261 : try
4262 : {
4263 :
4264 : // Ugly hack: papszArgv may not contain the processing mode if coming from
4265 : // SWIG we have not other option than to check them all
4266 : const static std::list<std::string> processingModes{
4267 : {"hillshade", "slope", "aspect", "color-relief", "TRI", "TPI",
4268 319 : "roughness"}};
4269 :
4270 165 : if (bNoProcessingMode)
4271 : {
4272 : try
4273 : {
4274 144 : tryHandleArgv(aosArgv);
4275 : }
4276 288 : catch (std::exception &)
4277 : {
4278 144 : bool bSuccess = false;
4279 228 : for (const auto &processingMode : processingModes)
4280 : {
4281 228 : CPLStringList aosArgvTmp{aosArgv};
4282 228 : CPL_IGNORE_RET_VAL(aosArgv);
4283 228 : aosArgvTmp.InsertString(0, processingMode.c_str());
4284 : try
4285 : {
4286 228 : tryHandleArgv(aosArgvTmp);
4287 144 : bSuccess = true;
4288 144 : break;
4289 : }
4290 126 : catch (std::runtime_error &)
4291 : {
4292 63 : continue;
4293 : }
4294 42 : catch (std::invalid_argument &)
4295 : {
4296 21 : continue;
4297 : }
4298 0 : catch (std::logic_error &)
4299 : {
4300 0 : continue;
4301 : }
4302 : }
4303 :
4304 144 : if (!bSuccess)
4305 : {
4306 : throw std::invalid_argument(
4307 0 : "Argument(s) are not valid with any processing mode.");
4308 : }
4309 : }
4310 : }
4311 : else
4312 : {
4313 21 : tryHandleArgv(aosArgv);
4314 : }
4315 : }
4316 0 : catch (const std::exception &err)
4317 : {
4318 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
4319 0 : err.what());
4320 0 : return nullptr;
4321 : }
4322 :
4323 165 : if (!std::isnan(psOptions->globalScale))
4324 : {
4325 25 : if (!std::isnan(psOptions->xscale) || !std::isnan(psOptions->yscale))
4326 : {
4327 2 : CPLError(CE_Failure, CPLE_AppDefined,
4328 : "-scale and -xscale/-yscale are mutually exclusive.");
4329 2 : return nullptr;
4330 : }
4331 23 : psOptions->xscale = psOptions->globalScale;
4332 23 : psOptions->yscale = psOptions->globalScale;
4333 : }
4334 140 : else if ((!std::isnan(psOptions->xscale)) ^
4335 140 : (!std::isnan(psOptions->yscale)))
4336 : {
4337 2 : CPLError(CE_Failure, CPLE_AppDefined,
4338 : "When one of -xscale or -yscale is specified, both must be "
4339 : "specified.");
4340 2 : return nullptr;
4341 : }
4342 :
4343 161 : return psOptions.release();
4344 : }
4345 :
4346 : /************************************************************************/
4347 : /* GDALDEMProcessingOptionsFree() */
4348 : /************************************************************************/
4349 :
4350 : /**
4351 : * Frees the GDALDEMProcessingOptions struct.
4352 : *
4353 : * @param psOptions the options struct for GDALDEMProcessing().
4354 : *
4355 : * @since GDAL 2.1
4356 : */
4357 :
4358 161 : void GDALDEMProcessingOptionsFree(GDALDEMProcessingOptions *psOptions)
4359 : {
4360 161 : delete psOptions;
4361 161 : }
4362 :
4363 : /************************************************************************/
4364 : /* GDALDEMProcessingOptionsSetProgress() */
4365 : /************************************************************************/
4366 :
4367 : /**
4368 : * Set a progress function.
4369 : *
4370 : * @param psOptions the options struct for GDALDEMProcessing().
4371 : * @param pfnProgress the progress callback.
4372 : * @param pProgressData the user data for the progress callback.
4373 : *
4374 : * @since GDAL 2.1
4375 : */
4376 :
4377 42 : void GDALDEMProcessingOptionsSetProgress(GDALDEMProcessingOptions *psOptions,
4378 : GDALProgressFunc pfnProgress,
4379 : void *pProgressData)
4380 : {
4381 42 : psOptions->pfnProgress = pfnProgress;
4382 42 : psOptions->pProgressData = pProgressData;
4383 42 : }
|
