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