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