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