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