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