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