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