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