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