Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: High Performance Image Reprojector
4 : * Purpose: Implementation of high level convenience APIs for warper.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdalwarper.h"
16 :
17 : #include <stdlib.h>
18 : #include <string.h>
19 :
20 : #include <algorithm>
21 : #include <cmath>
22 : #include <limits>
23 :
24 : #include "cpl_conv.h"
25 : #include "cpl_error.h"
26 : #include "cpl_float.h"
27 : #include "cpl_mask.h"
28 : #include "cpl_minixml.h"
29 : #include "cpl_progress.h"
30 : #include "cpl_string.h"
31 : #include "cpl_vsi.h"
32 : #include "gdal.h"
33 : #include "gdal_priv.h"
34 : #include "ogr_api.h"
35 : #include "ogr_core.h"
36 : #include "vrtdataset.h" // for VRTSerializeNoData
37 :
38 : #if (defined(__x86_64) || defined(_M_X64))
39 : #include <emmintrin.h>
40 : #endif
41 :
42 : /************************************************************************/
43 : /* GDALReprojectImage() */
44 : /************************************************************************/
45 :
46 : /**
47 : * Reproject image.
48 : *
49 : * This is a convenience function utilizing the GDALWarpOperation class to
50 : * reproject an image from a source to a destination. In particular, this
51 : * function takes care of establishing the transformation function to
52 : * implement the reprojection, and will default a variety of other
53 : * warp options.
54 : *
55 : * Nodata values set on destination dataset are taken into account.
56 : *
57 : * No metadata, projection info, or color tables are transferred
58 : * to the output file. Source overviews are not considered.
59 : *
60 : * For more advanced warping capabilities, consider using GDALWarp().
61 : *
62 : * @param hSrcDS the source image file.
63 : * @param pszSrcWKT the source projection. If NULL the source projection
64 : * is read from from hSrcDS.
65 : * @param hDstDS the destination image file.
66 : * @param pszDstWKT the destination projection. If NULL the destination
67 : * projection will be read from hDstDS.
68 : * @param eResampleAlg the type of resampling to use.
69 : * @param dfWarpMemoryLimit the amount of memory (in bytes) that the warp
70 : * API is allowed to use for caching. This is in addition to the memory
71 : * already allocated to the GDAL caching (as per GDALSetCacheMax()). May be
72 : * 0.0 to use default memory settings.
73 : * @param dfMaxError maximum error measured in input pixels that is allowed
74 : * in approximating the transformation (0.0 for exact calculations).
75 : * @param pfnProgress a GDALProgressFunc() compatible callback function for
76 : * reporting progress or NULL.
77 : * @param pProgressArg argument to be passed to pfnProgress. May be NULL.
78 : * @param psOptions warp options, normally NULL.
79 : *
80 : * @return CE_None on success or CE_Failure if something goes wrong.
81 : * @see GDALWarp()
82 : */
83 :
84 90 : CPLErr CPL_STDCALL GDALReprojectImage(
85 : GDALDatasetH hSrcDS, const char *pszSrcWKT, GDALDatasetH hDstDS,
86 : const char *pszDstWKT, GDALResampleAlg eResampleAlg,
87 : CPL_UNUSED double dfWarpMemoryLimit, double dfMaxError,
88 : GDALProgressFunc pfnProgress, void *pProgressArg,
89 : GDALWarpOptions *psOptions)
90 :
91 : {
92 : /* -------------------------------------------------------------------- */
93 : /* Setup a reprojection based transformer. */
94 : /* -------------------------------------------------------------------- */
95 90 : void *hTransformArg = GDALCreateGenImgProjTransformer(
96 : hSrcDS, pszSrcWKT, hDstDS, pszDstWKT, TRUE, 1000.0, 0);
97 :
98 90 : if (hTransformArg == nullptr)
99 0 : return CE_Failure;
100 :
101 : /* -------------------------------------------------------------------- */
102 : /* Create a copy of the user provided options, or a defaulted */
103 : /* options structure. */
104 : /* -------------------------------------------------------------------- */
105 : GDALWarpOptions *psWOptions = psOptions == nullptr
106 90 : ? GDALCreateWarpOptions()
107 1 : : GDALCloneWarpOptions(psOptions);
108 :
109 90 : psWOptions->eResampleAlg = eResampleAlg;
110 :
111 : /* -------------------------------------------------------------------- */
112 : /* Set transform. */
113 : /* -------------------------------------------------------------------- */
114 90 : if (dfMaxError > 0.0)
115 : {
116 3 : psWOptions->pTransformerArg = GDALCreateApproxTransformer(
117 : GDALGenImgProjTransform, hTransformArg, dfMaxError);
118 :
119 3 : psWOptions->pfnTransformer = GDALApproxTransform;
120 : }
121 : else
122 : {
123 87 : psWOptions->pfnTransformer = GDALGenImgProjTransform;
124 87 : psWOptions->pTransformerArg = hTransformArg;
125 : }
126 :
127 : /* -------------------------------------------------------------------- */
128 : /* Set file and band mapping. */
129 : /* -------------------------------------------------------------------- */
130 90 : psWOptions->hSrcDS = hSrcDS;
131 90 : psWOptions->hDstDS = hDstDS;
132 :
133 90 : int nSrcBands = GDALGetRasterCount(hSrcDS);
134 : {
135 90 : GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, nSrcBands);
136 90 : if (hBand && GDALGetRasterColorInterpretation(hBand) == GCI_AlphaBand)
137 : {
138 12 : psWOptions->nSrcAlphaBand = nSrcBands;
139 12 : nSrcBands--;
140 : }
141 : }
142 :
143 90 : int nDstBands = GDALGetRasterCount(hDstDS);
144 : {
145 90 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, nDstBands);
146 90 : if (hBand && GDALGetRasterColorInterpretation(hBand) == GCI_AlphaBand)
147 : {
148 12 : psWOptions->nDstAlphaBand = nDstBands;
149 12 : nDstBands--;
150 : }
151 : }
152 :
153 90 : GDALWarpInitDefaultBandMapping(psWOptions, std::min(nSrcBands, nDstBands));
154 :
155 : /* -------------------------------------------------------------------- */
156 : /* Set source nodata values if the source dataset seems to have */
157 : /* any. Same for target nodata values */
158 : /* -------------------------------------------------------------------- */
159 204 : for (int iBand = 0; iBand < psWOptions->nBandCount; iBand++)
160 : {
161 114 : GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, iBand + 1);
162 :
163 114 : int bGotNoData = FALSE;
164 114 : double dfNoDataValue = GDALGetRasterNoDataValue(hBand, &bGotNoData);
165 114 : if (bGotNoData)
166 : {
167 7 : GDALWarpInitSrcNoDataReal(psWOptions, -1.1e20);
168 7 : psWOptions->padfSrcNoDataReal[iBand] = dfNoDataValue;
169 : }
170 :
171 : // Deal with target band.
172 114 : hBand = GDALGetRasterBand(hDstDS, iBand + 1);
173 :
174 114 : dfNoDataValue = GDALGetRasterNoDataValue(hBand, &bGotNoData);
175 114 : if (bGotNoData)
176 : {
177 2 : GDALWarpInitDstNoDataReal(psWOptions, -1.1e20);
178 2 : psWOptions->padfDstNoDataReal[iBand] = dfNoDataValue;
179 : }
180 : }
181 :
182 : /* -------------------------------------------------------------------- */
183 : /* Set the progress function. */
184 : /* -------------------------------------------------------------------- */
185 90 : if (pfnProgress != nullptr)
186 : {
187 3 : psWOptions->pfnProgress = pfnProgress;
188 3 : psWOptions->pProgressArg = pProgressArg;
189 : }
190 :
191 : /* -------------------------------------------------------------------- */
192 : /* Create a warp options based on the options. */
193 : /* -------------------------------------------------------------------- */
194 90 : GDALWarpOperation oWarper;
195 90 : CPLErr eErr = oWarper.Initialize(psWOptions);
196 :
197 90 : if (eErr == CE_None)
198 90 : eErr = oWarper.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS),
199 : GDALGetRasterYSize(hDstDS));
200 :
201 : /* -------------------------------------------------------------------- */
202 : /* Cleanup. */
203 : /* -------------------------------------------------------------------- */
204 90 : GDALDestroyGenImgProjTransformer(hTransformArg);
205 :
206 90 : if (dfMaxError > 0.0)
207 3 : GDALDestroyApproxTransformer(psWOptions->pTransformerArg);
208 :
209 90 : GDALDestroyWarpOptions(psWOptions);
210 :
211 90 : return eErr;
212 : }
213 :
214 : /************************************************************************/
215 : /* GDALCreateAndReprojectImage() */
216 : /* */
217 : /* This is a "quicky" reprojection API. */
218 : /************************************************************************/
219 :
220 : /** Reproject an image and create the target reprojected image */
221 0 : CPLErr CPL_STDCALL GDALCreateAndReprojectImage(
222 : GDALDatasetH hSrcDS, const char *pszSrcWKT, const char *pszDstFilename,
223 : const char *pszDstWKT, GDALDriverH hDstDriver, char **papszCreateOptions,
224 : GDALResampleAlg eResampleAlg, double dfWarpMemoryLimit, double dfMaxError,
225 : GDALProgressFunc pfnProgress, void *pProgressArg,
226 : GDALWarpOptions *psOptions)
227 :
228 : {
229 0 : VALIDATE_POINTER1(hSrcDS, "GDALCreateAndReprojectImage", CE_Failure);
230 :
231 : /* -------------------------------------------------------------------- */
232 : /* Default a few parameters. */
233 : /* -------------------------------------------------------------------- */
234 0 : if (hDstDriver == nullptr)
235 : {
236 0 : hDstDriver = GDALGetDriverByName("GTiff");
237 0 : if (hDstDriver == nullptr)
238 : {
239 0 : CPLError(CE_Failure, CPLE_AppDefined,
240 : "GDALCreateAndReprojectImage needs GTiff driver");
241 0 : return CE_Failure;
242 : }
243 : }
244 :
245 0 : if (pszSrcWKT == nullptr)
246 0 : pszSrcWKT = GDALGetProjectionRef(hSrcDS);
247 :
248 0 : if (pszDstWKT == nullptr)
249 0 : pszDstWKT = pszSrcWKT;
250 :
251 : /* -------------------------------------------------------------------- */
252 : /* Create a transformation object from the source to */
253 : /* destination coordinate system. */
254 : /* -------------------------------------------------------------------- */
255 0 : void *hTransformArg = GDALCreateGenImgProjTransformer(
256 : hSrcDS, pszSrcWKT, nullptr, pszDstWKT, TRUE, 1000.0, 0);
257 :
258 0 : if (hTransformArg == nullptr)
259 0 : return CE_Failure;
260 :
261 : /* -------------------------------------------------------------------- */
262 : /* Get approximate output definition. */
263 : /* -------------------------------------------------------------------- */
264 0 : double adfDstGeoTransform[6] = {};
265 0 : int nPixels = 0;
266 0 : int nLines = 0;
267 :
268 0 : if (GDALSuggestedWarpOutput(hSrcDS, GDALGenImgProjTransform, hTransformArg,
269 : adfDstGeoTransform, &nPixels,
270 0 : &nLines) != CE_None)
271 0 : return CE_Failure;
272 :
273 0 : GDALDestroyGenImgProjTransformer(hTransformArg);
274 :
275 : /* -------------------------------------------------------------------- */
276 : /* Create the output file. */
277 : /* -------------------------------------------------------------------- */
278 0 : GDALDatasetH hDstDS = GDALCreate(
279 : hDstDriver, pszDstFilename, nPixels, nLines, GDALGetRasterCount(hSrcDS),
280 : GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1)),
281 : papszCreateOptions);
282 :
283 0 : if (hDstDS == nullptr)
284 0 : return CE_Failure;
285 :
286 : /* -------------------------------------------------------------------- */
287 : /* Write out the projection definition. */
288 : /* -------------------------------------------------------------------- */
289 0 : GDALSetProjection(hDstDS, pszDstWKT);
290 0 : GDALSetGeoTransform(hDstDS, adfDstGeoTransform);
291 :
292 : /* -------------------------------------------------------------------- */
293 : /* Perform the reprojection. */
294 : /* -------------------------------------------------------------------- */
295 0 : CPLErr eErr = GDALReprojectImage(
296 : hSrcDS, pszSrcWKT, hDstDS, pszDstWKT, eResampleAlg, dfWarpMemoryLimit,
297 : dfMaxError, pfnProgress, pProgressArg, psOptions);
298 :
299 0 : GDALClose(hDstDS);
300 :
301 0 : return eErr;
302 : }
303 :
304 : /************************************************************************/
305 : /* GDALWarpNoDataMaskerT() */
306 : /************************************************************************/
307 :
308 : template <class T>
309 714 : static CPLErr GDALWarpNoDataMaskerT(const double *padfNoData, size_t nPixels,
310 : const T *pData, GUInt32 *panValidityMask,
311 : int *pbOutAllValid)
312 : {
313 : // Nothing to do if value is out of range.
314 714 : if (padfNoData[0] < cpl::NumericLimits<T>::min() ||
315 1426 : padfNoData[0] > cpl::NumericLimits<T>::max() + 0.000001 ||
316 712 : padfNoData[1] != 0.0)
317 : {
318 2 : *pbOutAllValid = TRUE;
319 2 : return CE_None;
320 : }
321 :
322 712 : const int nNoData = static_cast<int>(floor(padfNoData[0] + 0.000001));
323 712 : int bAllValid = TRUE;
324 92288679 : for (size_t iOffset = 0; iOffset < nPixels; ++iOffset)
325 : {
326 92287921 : if (pData[iOffset] == nNoData)
327 : {
328 87182425 : bAllValid = FALSE;
329 87182425 : CPLMaskClear(panValidityMask, iOffset);
330 : }
331 : }
332 712 : *pbOutAllValid = bAllValid;
333 :
334 712 : return CE_None;
335 : }
336 :
337 : /************************************************************************/
338 : /* GDALWarpNoDataMasker() */
339 : /* */
340 : /* GDALMaskFunc for establishing a validity mask for a source */
341 : /* band based on a provided NODATA value. */
342 : /************************************************************************/
343 :
344 839 : CPLErr GDALWarpNoDataMasker(void *pMaskFuncArg, int nBandCount,
345 : GDALDataType eType, int /* nXOff */,
346 : int /* nYOff */, int nXSize, int nYSize,
347 : GByte **ppImageData, int bMaskIsFloat,
348 : void *pValidityMask, int *pbOutAllValid)
349 :
350 : {
351 839 : const double *padfNoData = static_cast<double *>(pMaskFuncArg);
352 839 : GUInt32 *panValidityMask = static_cast<GUInt32 *>(pValidityMask);
353 839 : const size_t nPixels = static_cast<size_t>(nXSize) * nYSize;
354 :
355 839 : *pbOutAllValid = FALSE;
356 :
357 839 : if (nBandCount != 1 || bMaskIsFloat)
358 : {
359 0 : CPLError(
360 : CE_Failure, CPLE_AppDefined,
361 : "Invalid nBandCount or bMaskIsFloat argument in SourceNoDataMask");
362 0 : return CE_Failure;
363 : }
364 :
365 839 : switch (eType)
366 : {
367 654 : case GDT_Byte:
368 654 : return GDALWarpNoDataMaskerT(padfNoData, nPixels,
369 : *ppImageData, // Already a GByte *.
370 654 : panValidityMask, pbOutAllValid);
371 :
372 56 : case GDT_Int16:
373 56 : return GDALWarpNoDataMaskerT(
374 : padfNoData, nPixels, reinterpret_cast<GInt16 *>(*ppImageData),
375 56 : panValidityMask, pbOutAllValid);
376 :
377 4 : case GDT_UInt16:
378 4 : return GDALWarpNoDataMaskerT(
379 : padfNoData, nPixels, reinterpret_cast<GUInt16 *>(*ppImageData),
380 4 : panValidityMask, pbOutAllValid);
381 :
382 78 : case GDT_Float32:
383 : {
384 78 : const float fNoData = static_cast<float>(padfNoData[0]);
385 78 : const float *pafData = reinterpret_cast<float *>(*ppImageData);
386 78 : const bool bIsNoDataNan = CPL_TO_BOOL(std::isnan(fNoData));
387 :
388 : // Nothing to do if value is out of range.
389 78 : if (padfNoData[1] != 0.0)
390 : {
391 0 : *pbOutAllValid = TRUE;
392 0 : return CE_None;
393 : }
394 :
395 78 : int bAllValid = TRUE;
396 45276 : for (size_t iOffset = 0; iOffset < nPixels; ++iOffset)
397 : {
398 45198 : float fVal = pafData[iOffset];
399 90166 : if ((bIsNoDataNan && std::isnan(fVal)) ||
400 44968 : (!bIsNoDataNan && ARE_REAL_EQUAL(fVal, fNoData)))
401 : {
402 37394 : bAllValid = FALSE;
403 37394 : CPLMaskClear(panValidityMask, iOffset);
404 : }
405 : }
406 78 : *pbOutAllValid = bAllValid;
407 : }
408 78 : break;
409 :
410 19 : case GDT_Float64:
411 : {
412 19 : const double dfNoData = padfNoData[0];
413 19 : const double *padfData = reinterpret_cast<double *>(*ppImageData);
414 19 : const bool bIsNoDataNan = CPL_TO_BOOL(std::isnan(dfNoData));
415 :
416 : // Nothing to do if value is out of range.
417 19 : if (padfNoData[1] != 0.0)
418 : {
419 0 : *pbOutAllValid = TRUE;
420 0 : return CE_None;
421 : }
422 :
423 19 : int bAllValid = TRUE;
424 2905 : for (size_t iOffset = 0; iOffset < nPixels; ++iOffset)
425 : {
426 2886 : double dfVal = padfData[iOffset];
427 3401 : if ((bIsNoDataNan && std::isnan(dfVal)) ||
428 515 : (!bIsNoDataNan && ARE_REAL_EQUAL(dfVal, dfNoData)))
429 : {
430 2502 : bAllValid = FALSE;
431 2502 : CPLMaskClear(panValidityMask, iOffset);
432 : }
433 : }
434 19 : *pbOutAllValid = bAllValid;
435 : }
436 19 : break;
437 :
438 28 : default:
439 : {
440 28 : const int nWordSize = GDALGetDataTypeSizeBytes(eType);
441 :
442 : const bool bIsNoDataRealNan =
443 28 : CPL_TO_BOOL(std::isnan(padfNoData[0]));
444 :
445 : double *padfWrk =
446 28 : static_cast<double *>(CPLMalloc(nXSize * sizeof(double) * 2));
447 28 : int bAllValid = TRUE;
448 314 : for (int iLine = 0; iLine < nYSize; iLine++)
449 : {
450 286 : GDALCopyWords((*ppImageData) + nWordSize * iLine * nXSize,
451 : eType, nWordSize, padfWrk, GDT_CFloat64, 16,
452 : nXSize);
453 :
454 6522 : for (int iPixel = 0; iPixel < nXSize; ++iPixel)
455 : {
456 0 : if (((bIsNoDataRealNan &&
457 12472 : std::isnan(padfWrk[iPixel * 2])) ||
458 12472 : (!bIsNoDataRealNan &&
459 6236 : ARE_REAL_EQUAL(padfWrk[iPixel * 2], padfNoData[0]))))
460 : {
461 933 : size_t iOffset =
462 933 : iPixel + static_cast<size_t>(iLine) * nXSize;
463 :
464 933 : bAllValid = FALSE;
465 933 : CPLMaskClear(panValidityMask, iOffset);
466 : }
467 : }
468 : }
469 28 : *pbOutAllValid = bAllValid;
470 :
471 28 : CPLFree(padfWrk);
472 : }
473 28 : break;
474 : }
475 :
476 125 : return CE_None;
477 : }
478 :
479 : /************************************************************************/
480 : /* GDALWarpSrcAlphaMasker() */
481 : /* */
482 : /* GDALMaskFunc for reading source simple 8bit alpha mask */
483 : /* information and building a floating point density mask from */
484 : /* it. */
485 : /************************************************************************/
486 :
487 79 : CPLErr GDALWarpSrcAlphaMasker(void *pMaskFuncArg, int /* nBandCount */,
488 : GDALDataType /* eType */, int nXOff, int nYOff,
489 : int nXSize, int nYSize, GByte ** /*ppImageData */,
490 : int bMaskIsFloat, void *pValidityMask,
491 : int *pbOutAllOpaque)
492 :
493 : {
494 79 : GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg);
495 79 : float *pafMask = static_cast<float *>(pValidityMask);
496 79 : *pbOutAllOpaque = FALSE;
497 79 : const size_t nPixels = static_cast<size_t>(nXSize) * nYSize;
498 :
499 : /* -------------------------------------------------------------------- */
500 : /* Do some minimal checking. */
501 : /* -------------------------------------------------------------------- */
502 79 : if (!bMaskIsFloat)
503 : {
504 0 : CPLAssert(false);
505 : return CE_Failure;
506 : }
507 :
508 79 : if (psWO == nullptr || psWO->nSrcAlphaBand < 1)
509 : {
510 0 : CPLAssert(false);
511 : return CE_Failure;
512 : }
513 :
514 : /* -------------------------------------------------------------------- */
515 : /* Read the alpha band. */
516 : /* -------------------------------------------------------------------- */
517 : GDALRasterBandH hAlphaBand =
518 79 : GDALGetRasterBand(psWO->hSrcDS, psWO->nSrcAlphaBand);
519 79 : if (hAlphaBand == nullptr)
520 0 : return CE_Failure;
521 :
522 : // Rescale.
523 79 : const float inv_alpha_max = static_cast<float>(
524 79 : 1.0 / CPLAtof(CSLFetchNameValueDef(psWO->papszWarpOptions,
525 79 : "SRC_ALPHA_MAX", "255")));
526 79 : bool bOutAllOpaque = true;
527 :
528 79 : size_t iPixel = 0;
529 : CPLErr eErr;
530 :
531 : #if (defined(__x86_64) || defined(_M_X64))
532 79 : GDALDataType eDT = GDALGetRasterDataType(hAlphaBand);
533 : // Make sure that pafMask is at least 8-byte aligned, which should
534 : // normally be always the case if being a ptr returned by malloc().
535 79 : if ((eDT == GDT_Byte || eDT == GDT_UInt16) && CPL_IS_ALIGNED(pafMask, 8))
536 : {
537 : // Read data.
538 136 : eErr = GDALRasterIOEx(
539 : hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize, pafMask, nXSize,
540 : nYSize, eDT, static_cast<GSpacing>(sizeof(int)),
541 68 : static_cast<GSpacing>(sizeof(int)) * nXSize, nullptr);
542 :
543 68 : if (eErr != CE_None)
544 0 : return eErr;
545 :
546 : // Make sure we have the correct alignment before doing SSE
547 : // On Linux x86_64, the alignment should be always correct due
548 : // the alignment of malloc() being 16 byte.
549 68 : const GUInt32 mask = (eDT == GDT_Byte) ? 0xff : 0xffff;
550 68 : if (!CPL_IS_ALIGNED(pafMask, 16))
551 : {
552 0 : pafMask[iPixel] =
553 0 : (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) *
554 0 : inv_alpha_max;
555 0 : if (pafMask[iPixel] >= 1.0f)
556 0 : pafMask[iPixel] = 1.0f;
557 : else
558 0 : bOutAllOpaque = false;
559 0 : iPixel++;
560 : }
561 68 : CPLAssert(CPL_IS_ALIGNED(pafMask + iPixel, 16));
562 68 : const __m128 xmm_inverse_alpha_max = _mm_load1_ps(&inv_alpha_max);
563 68 : const float one_single = 1.0f;
564 68 : const __m128 xmm_one = _mm_load1_ps(&one_single);
565 136 : const __m128i xmm_i_mask = _mm_set1_epi32(mask);
566 68 : __m128 xmmMaskNonOpaque0 = _mm_setzero_ps();
567 68 : __m128 xmmMaskNonOpaque1 = _mm_setzero_ps();
568 68 : __m128 xmmMaskNonOpaque2 = _mm_setzero_ps();
569 629236 : for (; iPixel + 6 * 4 - 1 < nPixels; iPixel += 6 * 4)
570 : {
571 1258340 : __m128 xmm_mask0 = _mm_cvtepi32_ps(_mm_and_si128(
572 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
573 629168 : pafMask + iPixel + 4 * 0))));
574 1258340 : __m128 xmm_mask1 = _mm_cvtepi32_ps(_mm_and_si128(
575 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
576 629168 : pafMask + iPixel + 4 * 1))));
577 1258340 : __m128 xmm_mask2 = _mm_cvtepi32_ps(_mm_and_si128(
578 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
579 629168 : pafMask + iPixel + 4 * 2))));
580 1258340 : __m128 xmm_mask3 = _mm_cvtepi32_ps(_mm_and_si128(
581 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
582 629168 : pafMask + iPixel + 4 * 3))));
583 1258340 : __m128 xmm_mask4 = _mm_cvtepi32_ps(_mm_and_si128(
584 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
585 629168 : pafMask + iPixel + 4 * 4))));
586 1887500 : __m128 xmm_mask5 = _mm_cvtepi32_ps(_mm_and_si128(
587 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
588 629168 : pafMask + iPixel + 4 * 5))));
589 629168 : xmm_mask0 = _mm_mul_ps(xmm_mask0, xmm_inverse_alpha_max);
590 629168 : xmm_mask1 = _mm_mul_ps(xmm_mask1, xmm_inverse_alpha_max);
591 629168 : xmm_mask2 = _mm_mul_ps(xmm_mask2, xmm_inverse_alpha_max);
592 629168 : xmm_mask3 = _mm_mul_ps(xmm_mask3, xmm_inverse_alpha_max);
593 629168 : xmm_mask4 = _mm_mul_ps(xmm_mask4, xmm_inverse_alpha_max);
594 629168 : xmm_mask5 = _mm_mul_ps(xmm_mask5, xmm_inverse_alpha_max);
595 : xmmMaskNonOpaque0 =
596 1258340 : _mm_or_ps(xmmMaskNonOpaque0, _mm_cmplt_ps(xmm_mask0, xmm_one));
597 : xmmMaskNonOpaque1 =
598 1258340 : _mm_or_ps(xmmMaskNonOpaque1, _mm_cmplt_ps(xmm_mask1, xmm_one));
599 : xmmMaskNonOpaque2 =
600 1258340 : _mm_or_ps(xmmMaskNonOpaque2, _mm_cmplt_ps(xmm_mask2, xmm_one));
601 : xmmMaskNonOpaque0 =
602 1258340 : _mm_or_ps(xmmMaskNonOpaque0, _mm_cmplt_ps(xmm_mask3, xmm_one));
603 : xmmMaskNonOpaque1 =
604 1258340 : _mm_or_ps(xmmMaskNonOpaque1, _mm_cmplt_ps(xmm_mask4, xmm_one));
605 : xmmMaskNonOpaque2 =
606 1258340 : _mm_or_ps(xmmMaskNonOpaque2, _mm_cmplt_ps(xmm_mask5, xmm_one));
607 629168 : xmm_mask0 = _mm_min_ps(xmm_mask0, xmm_one);
608 629168 : xmm_mask1 = _mm_min_ps(xmm_mask1, xmm_one);
609 629168 : xmm_mask2 = _mm_min_ps(xmm_mask2, xmm_one);
610 629168 : xmm_mask3 = _mm_min_ps(xmm_mask3, xmm_one);
611 629168 : xmm_mask4 = _mm_min_ps(xmm_mask4, xmm_one);
612 629168 : xmm_mask5 = _mm_min_ps(xmm_mask5, xmm_one);
613 629168 : _mm_store_ps(pafMask + iPixel + 4 * 0, xmm_mask0);
614 629168 : _mm_store_ps(pafMask + iPixel + 4 * 1, xmm_mask1);
615 629168 : _mm_store_ps(pafMask + iPixel + 4 * 2, xmm_mask2);
616 629168 : _mm_store_ps(pafMask + iPixel + 4 * 3, xmm_mask3);
617 629168 : _mm_store_ps(pafMask + iPixel + 4 * 4, xmm_mask4);
618 629168 : _mm_store_ps(pafMask + iPixel + 4 * 5, xmm_mask5);
619 : }
620 136 : if (_mm_movemask_ps(
621 : _mm_or_ps(_mm_or_ps(xmmMaskNonOpaque0, xmmMaskNonOpaque1),
622 68 : xmmMaskNonOpaque2)))
623 : {
624 43 : bOutAllOpaque = false;
625 : }
626 980 : for (; iPixel < nPixels; iPixel++)
627 : {
628 912 : pafMask[iPixel] =
629 912 : (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) *
630 912 : inv_alpha_max;
631 912 : if (pafMask[iPixel] >= 1.0f)
632 442 : pafMask[iPixel] = 1.0f;
633 : else
634 470 : bOutAllOpaque = false;
635 68 : }
636 : }
637 : else
638 : #endif
639 : {
640 : // Read data.
641 11 : eErr = GDALRasterIO(hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize,
642 : pafMask, nXSize, nYSize, GDT_Float32, 0, 0);
643 :
644 11 : if (eErr != CE_None)
645 0 : return eErr;
646 :
647 : // TODO(rouault): Is loop unrolling by hand (r34564) actually helpful?
648 12969 : for (; iPixel + 3 < nPixels; iPixel += 4)
649 : {
650 12958 : pafMask[iPixel] = pafMask[iPixel] * inv_alpha_max;
651 12958 : if (pafMask[iPixel] >= 1.0f)
652 3028 : pafMask[iPixel] = 1.0f;
653 : else
654 9930 : bOutAllOpaque = false;
655 12958 : pafMask[iPixel + 1] = pafMask[iPixel + 1] * inv_alpha_max;
656 12958 : if (pafMask[iPixel + 1] >= 1.0f)
657 3018 : pafMask[iPixel + 1] = 1.0f;
658 : else
659 9940 : bOutAllOpaque = false;
660 12958 : pafMask[iPixel + 2] = pafMask[iPixel + 2] * inv_alpha_max;
661 12958 : if (pafMask[iPixel + 2] >= 1.0f)
662 3066 : pafMask[iPixel + 2] = 1.0f;
663 : else
664 9892 : bOutAllOpaque = false;
665 12958 : pafMask[iPixel + 3] = pafMask[iPixel + 3] * inv_alpha_max;
666 12958 : if (pafMask[iPixel + 3] >= 1.0f)
667 3016 : pafMask[iPixel + 3] = 1.0f;
668 : else
669 9942 : bOutAllOpaque = false;
670 : }
671 :
672 12 : for (; iPixel < nPixels; iPixel++)
673 : {
674 1 : pafMask[iPixel] = pafMask[iPixel] * inv_alpha_max;
675 1 : if (pafMask[iPixel] >= 1.0f)
676 0 : pafMask[iPixel] = 1.0f;
677 : else
678 1 : bOutAllOpaque = false;
679 : }
680 : }
681 :
682 79 : *pbOutAllOpaque = bOutAllOpaque;
683 :
684 79 : return CE_None;
685 : }
686 :
687 : /************************************************************************/
688 : /* GDALWarpSrcMaskMasker() */
689 : /* */
690 : /* GDALMaskFunc for reading source simple 8bit validity mask */
691 : /* information and building a one bit validity mask. */
692 : /************************************************************************/
693 :
694 2 : CPLErr GDALWarpSrcMaskMasker(void *pMaskFuncArg, int /* nBandCount */,
695 : GDALDataType /* eType */, int nXOff, int nYOff,
696 : int nXSize, int nYSize, GByte ** /*ppImageData */,
697 : int bMaskIsFloat, void *pValidityMask)
698 :
699 : {
700 2 : GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg);
701 2 : GUInt32 *panMask = static_cast<GUInt32 *>(pValidityMask);
702 :
703 : /* -------------------------------------------------------------------- */
704 : /* Do some minimal checking. */
705 : /* -------------------------------------------------------------------- */
706 2 : if (bMaskIsFloat)
707 : {
708 0 : CPLAssert(false);
709 : return CE_Failure;
710 : }
711 :
712 2 : if (psWO == nullptr)
713 : {
714 0 : CPLAssert(false);
715 : return CE_Failure;
716 : }
717 :
718 : /* -------------------------------------------------------------------- */
719 : /* Allocate a temporary buffer to read mask byte data into. */
720 : /* -------------------------------------------------------------------- */
721 : GByte *pabySrcMask =
722 2 : static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nXSize, nYSize));
723 2 : if (pabySrcMask == nullptr)
724 : {
725 0 : return CE_Failure;
726 : }
727 :
728 : /* -------------------------------------------------------------------- */
729 : /* Fetch our mask band. */
730 : /* -------------------------------------------------------------------- */
731 2 : GDALRasterBandH hMaskBand = nullptr;
732 : GDALRasterBandH hSrcBand =
733 2 : GDALGetRasterBand(psWO->hSrcDS, psWO->panSrcBands[0]);
734 2 : if (hSrcBand != nullptr)
735 2 : hMaskBand = GDALGetMaskBand(hSrcBand);
736 :
737 2 : if (hMaskBand == nullptr)
738 : {
739 0 : CPLAssert(false);
740 : return CE_Failure;
741 : }
742 :
743 : /* -------------------------------------------------------------------- */
744 : /* Read the mask band. */
745 : /* -------------------------------------------------------------------- */
746 2 : CPLErr eErr = GDALRasterIO(hMaskBand, GF_Read, nXOff, nYOff, nXSize, nYSize,
747 : pabySrcMask, nXSize, nYSize, GDT_Byte, 0, 0);
748 :
749 2 : if (eErr != CE_None)
750 : {
751 0 : CPLFree(pabySrcMask);
752 0 : return eErr;
753 : }
754 :
755 : /* -------------------------------------------------------------------- */
756 : /* Pack into 1 bit per pixel for validity. */
757 : /* -------------------------------------------------------------------- */
758 2 : const size_t nPixels = static_cast<size_t>(nXSize) * nYSize;
759 234579 : for (size_t iPixel = 0; iPixel < nPixels; iPixel++)
760 : {
761 234577 : if (pabySrcMask[iPixel] == 0)
762 31660 : CPLMaskClear(panMask, iPixel);
763 : }
764 :
765 2 : CPLFree(pabySrcMask);
766 :
767 2 : return CE_None;
768 : }
769 :
770 : /************************************************************************/
771 : /* GDALWarpDstAlphaMasker() */
772 : /* */
773 : /* GDALMaskFunc for reading or writing the destination simple */
774 : /* 8bit alpha mask information and building a floating point */
775 : /* density mask from it. Note, writing is distinguished */
776 : /* negative bandcount. */
777 : /************************************************************************/
778 :
779 490 : CPLErr GDALWarpDstAlphaMasker(void *pMaskFuncArg, int nBandCount,
780 : CPL_UNUSED GDALDataType /* eType */, int nXOff,
781 : int nYOff, int nXSize, int nYSize,
782 : GByte ** /*ppImageData */, int bMaskIsFloat,
783 : void *pValidityMask)
784 : {
785 : /* -------------------------------------------------------------------- */
786 : /* Do some minimal checking. */
787 : /* -------------------------------------------------------------------- */
788 490 : if (!bMaskIsFloat)
789 : {
790 0 : CPLAssert(false);
791 : return CE_Failure;
792 : }
793 :
794 490 : GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg);
795 490 : if (psWO == nullptr || psWO->nDstAlphaBand < 1)
796 : {
797 0 : CPLAssert(false);
798 : return CE_Failure;
799 : }
800 :
801 490 : float *pafMask = static_cast<float *>(pValidityMask);
802 490 : const size_t nPixels = static_cast<size_t>(nXSize) * nYSize;
803 :
804 : GDALRasterBandH hAlphaBand =
805 490 : GDALGetRasterBand(psWO->hDstDS, psWO->nDstAlphaBand);
806 490 : if (hAlphaBand == nullptr)
807 0 : return CE_Failure;
808 :
809 490 : size_t iPixel = 0;
810 :
811 : /* -------------------------------------------------------------------- */
812 : /* Read alpha case. */
813 : /* -------------------------------------------------------------------- */
814 490 : if (nBandCount >= 0)
815 : {
816 : const char *pszInitDest =
817 245 : CSLFetchNameValue(psWO->papszWarpOptions, "INIT_DEST");
818 :
819 : // Special logic for destinations being initialized on-the-fly.
820 245 : if (pszInitDest != nullptr)
821 : {
822 155 : memset(pafMask, 0, nPixels * sizeof(float));
823 155 : return CE_None;
824 : }
825 :
826 : // Rescale.
827 90 : const float inv_alpha_max = static_cast<float>(
828 90 : 1.0 / CPLAtof(CSLFetchNameValueDef(psWO->papszWarpOptions,
829 90 : "DST_ALPHA_MAX", "255")));
830 :
831 : #if (defined(__x86_64) || defined(_M_X64))
832 90 : const GDALDataType eDT = GDALGetRasterDataType(hAlphaBand);
833 : // Make sure that pafMask is at least 8-byte aligned, which should
834 : // normally be always the case if being a ptr returned by malloc().
835 90 : if ((eDT == GDT_Byte || eDT == GDT_UInt16) &&
836 81 : CPL_IS_ALIGNED(pafMask, 8))
837 : {
838 : // Read data.
839 162 : const CPLErr eErr = GDALRasterIOEx(
840 : hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize, pafMask,
841 : nXSize, nYSize, eDT, static_cast<GSpacing>(sizeof(int)),
842 81 : static_cast<GSpacing>(sizeof(int)) * nXSize, nullptr);
843 :
844 81 : if (eErr != CE_None)
845 0 : return eErr;
846 :
847 : // Make sure we have the correct alignment before doing SSE
848 : // On Linux x86_64, the alignment should be always correct due
849 : // the alignment of malloc() being 16 byte.
850 81 : const GUInt32 mask = (eDT == GDT_Byte) ? 0xff : 0xffff;
851 81 : if (!CPL_IS_ALIGNED(pafMask, 16))
852 : {
853 0 : pafMask[iPixel] =
854 0 : (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) *
855 0 : inv_alpha_max;
856 0 : pafMask[iPixel] = std::min(1.0f, pafMask[iPixel]);
857 0 : iPixel++;
858 : }
859 81 : CPLAssert(CPL_IS_ALIGNED(pafMask + iPixel, 16));
860 81 : const __m128 xmm_inverse_alpha_max = _mm_load1_ps(&inv_alpha_max);
861 81 : const float one_single = 1.0f;
862 81 : const __m128 xmm_one = _mm_load1_ps(&one_single);
863 162 : const __m128i xmm_i_mask = _mm_set1_epi32(mask);
864 126073 : for (; iPixel + 31 < nPixels; iPixel += 32)
865 : {
866 251984 : __m128 xmm_mask0 = _mm_cvtepi32_ps(_mm_and_si128(
867 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
868 125992 : pafMask + iPixel + 4 * 0))));
869 251984 : __m128 xmm_mask1 = _mm_cvtepi32_ps(_mm_and_si128(
870 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
871 125992 : pafMask + iPixel + 4 * 1))));
872 251984 : __m128 xmm_mask2 = _mm_cvtepi32_ps(_mm_and_si128(
873 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
874 125992 : pafMask + iPixel + 4 * 2))));
875 251984 : __m128 xmm_mask3 = _mm_cvtepi32_ps(_mm_and_si128(
876 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
877 125992 : pafMask + iPixel + 4 * 3))));
878 251984 : __m128 xmm_mask4 = _mm_cvtepi32_ps(_mm_and_si128(
879 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
880 125992 : pafMask + iPixel + 4 * 4))));
881 251984 : __m128 xmm_mask5 = _mm_cvtepi32_ps(_mm_and_si128(
882 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
883 125992 : pafMask + iPixel + 4 * 5))));
884 251984 : __m128 xmm_mask6 = _mm_cvtepi32_ps(_mm_and_si128(
885 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
886 125992 : pafMask + iPixel + 4 * 6))));
887 377976 : __m128 xmm_mask7 = _mm_cvtepi32_ps(_mm_and_si128(
888 : xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>(
889 125992 : pafMask + iPixel + 4 * 7))));
890 125992 : xmm_mask0 = _mm_mul_ps(xmm_mask0, xmm_inverse_alpha_max);
891 125992 : xmm_mask1 = _mm_mul_ps(xmm_mask1, xmm_inverse_alpha_max);
892 125992 : xmm_mask2 = _mm_mul_ps(xmm_mask2, xmm_inverse_alpha_max);
893 125992 : xmm_mask3 = _mm_mul_ps(xmm_mask3, xmm_inverse_alpha_max);
894 125992 : xmm_mask4 = _mm_mul_ps(xmm_mask4, xmm_inverse_alpha_max);
895 125992 : xmm_mask5 = _mm_mul_ps(xmm_mask5, xmm_inverse_alpha_max);
896 125992 : xmm_mask6 = _mm_mul_ps(xmm_mask6, xmm_inverse_alpha_max);
897 125992 : xmm_mask7 = _mm_mul_ps(xmm_mask7, xmm_inverse_alpha_max);
898 125992 : xmm_mask0 = _mm_min_ps(xmm_mask0, xmm_one);
899 125992 : xmm_mask1 = _mm_min_ps(xmm_mask1, xmm_one);
900 125992 : xmm_mask2 = _mm_min_ps(xmm_mask2, xmm_one);
901 125992 : xmm_mask3 = _mm_min_ps(xmm_mask3, xmm_one);
902 125992 : xmm_mask4 = _mm_min_ps(xmm_mask4, xmm_one);
903 125992 : xmm_mask5 = _mm_min_ps(xmm_mask5, xmm_one);
904 125992 : xmm_mask6 = _mm_min_ps(xmm_mask6, xmm_one);
905 125992 : xmm_mask7 = _mm_min_ps(xmm_mask7, xmm_one);
906 125992 : _mm_store_ps(pafMask + iPixel + 4 * 0, xmm_mask0);
907 125992 : _mm_store_ps(pafMask + iPixel + 4 * 1, xmm_mask1);
908 125992 : _mm_store_ps(pafMask + iPixel + 4 * 2, xmm_mask2);
909 125992 : _mm_store_ps(pafMask + iPixel + 4 * 3, xmm_mask3);
910 125992 : _mm_store_ps(pafMask + iPixel + 4 * 4, xmm_mask4);
911 125992 : _mm_store_ps(pafMask + iPixel + 4 * 5, xmm_mask5);
912 125992 : _mm_store_ps(pafMask + iPixel + 4 * 6, xmm_mask6);
913 125992 : _mm_store_ps(pafMask + iPixel + 4 * 7, xmm_mask7);
914 : }
915 825 : for (; iPixel < nPixels; iPixel++)
916 : {
917 744 : pafMask[iPixel] =
918 744 : (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) *
919 744 : inv_alpha_max;
920 744 : pafMask[iPixel] = std::min(1.0f, pafMask[iPixel]);
921 81 : }
922 : }
923 : else
924 : #endif
925 : {
926 : // Read data.
927 : const CPLErr eErr =
928 9 : GDALRasterIO(hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize,
929 : pafMask, nXSize, nYSize, GDT_Float32, 0, 0);
930 :
931 9 : if (eErr != CE_None)
932 0 : return eErr;
933 :
934 842 : for (; iPixel < nPixels; iPixel++)
935 : {
936 833 : pafMask[iPixel] = pafMask[iPixel] * inv_alpha_max;
937 833 : pafMask[iPixel] = std::min(1.0f, pafMask[iPixel]);
938 : }
939 : }
940 :
941 90 : return CE_None;
942 : }
943 :
944 : /* -------------------------------------------------------------------- */
945 : /* Write alpha case. */
946 : /* -------------------------------------------------------------------- */
947 : else
948 : {
949 245 : GDALDataType eDT = GDALGetRasterDataType(hAlphaBand);
950 : const float cst_alpha_max =
951 245 : static_cast<float>(CPLAtof(CSLFetchNameValueDef(
952 245 : psWO->papszWarpOptions, "DST_ALPHA_MAX", "255"))) +
953 29 : ((eDT == GDT_Byte || eDT == GDT_Int16 || eDT == GDT_UInt16 ||
954 1 : eDT == GDT_Int32 || eDT == GDT_UInt32)
955 274 : ? 0.1f
956 245 : : 0.0f);
957 :
958 245 : CPLErr eErr = CE_None;
959 :
960 : #if (defined(__x86_64) || defined(_M_X64))
961 : // Make sure that pafMask is at least 8-byte aligned, which should
962 : // normally be always the case if being a ptr returned by malloc()
963 245 : if ((eDT == GDT_Byte || eDT == GDT_Int16 || eDT == GDT_UInt16) &&
964 244 : CPL_IS_ALIGNED(pafMask, 8))
965 : {
966 : // Make sure we have the correct alignment before doing SSE
967 : // On Linux x86_64, the alignment should be always correct due
968 : // the alignment of malloc() being 16 byte
969 244 : if (!CPL_IS_ALIGNED(pafMask, 16))
970 : {
971 0 : reinterpret_cast<int *>(pafMask)[iPixel] =
972 0 : static_cast<int>(pafMask[iPixel] * cst_alpha_max);
973 0 : iPixel++;
974 : }
975 244 : CPLAssert(CPL_IS_ALIGNED(pafMask + iPixel, 16));
976 244 : const __m128 xmm_alpha_max = _mm_load1_ps(&cst_alpha_max);
977 279059 : for (; iPixel + 31 < nPixels; iPixel += 32)
978 : {
979 278815 : __m128 xmm_mask0 = _mm_load_ps(pafMask + iPixel + 4 * 0);
980 278815 : __m128 xmm_mask1 = _mm_load_ps(pafMask + iPixel + 4 * 1);
981 278815 : __m128 xmm_mask2 = _mm_load_ps(pafMask + iPixel + 4 * 2);
982 278815 : __m128 xmm_mask3 = _mm_load_ps(pafMask + iPixel + 4 * 3);
983 278815 : __m128 xmm_mask4 = _mm_load_ps(pafMask + iPixel + 4 * 4);
984 278815 : __m128 xmm_mask5 = _mm_load_ps(pafMask + iPixel + 4 * 5);
985 278815 : __m128 xmm_mask6 = _mm_load_ps(pafMask + iPixel + 4 * 6);
986 557630 : __m128 xmm_mask7 = _mm_load_ps(pafMask + iPixel + 4 * 7);
987 278815 : xmm_mask0 = _mm_mul_ps(xmm_mask0, xmm_alpha_max);
988 278815 : xmm_mask1 = _mm_mul_ps(xmm_mask1, xmm_alpha_max);
989 278815 : xmm_mask2 = _mm_mul_ps(xmm_mask2, xmm_alpha_max);
990 278815 : xmm_mask3 = _mm_mul_ps(xmm_mask3, xmm_alpha_max);
991 278815 : xmm_mask4 = _mm_mul_ps(xmm_mask4, xmm_alpha_max);
992 278815 : xmm_mask5 = _mm_mul_ps(xmm_mask5, xmm_alpha_max);
993 278815 : xmm_mask6 = _mm_mul_ps(xmm_mask6, xmm_alpha_max);
994 278815 : xmm_mask7 = _mm_mul_ps(xmm_mask7, xmm_alpha_max);
995 : // Truncate to int.
996 278815 : _mm_store_si128(
997 278815 : reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 0),
998 : _mm_cvttps_epi32(xmm_mask0));
999 278815 : _mm_store_si128(
1000 278815 : reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 1),
1001 : _mm_cvttps_epi32(xmm_mask1));
1002 278815 : _mm_store_si128(
1003 278815 : reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 2),
1004 : _mm_cvttps_epi32(xmm_mask2));
1005 278815 : _mm_store_si128(
1006 278815 : reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 3),
1007 : _mm_cvttps_epi32(xmm_mask3));
1008 278815 : _mm_store_si128(
1009 278815 : reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 4),
1010 : _mm_cvttps_epi32(xmm_mask4));
1011 278815 : _mm_store_si128(
1012 278815 : reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 5),
1013 : _mm_cvttps_epi32(xmm_mask5));
1014 278815 : _mm_store_si128(
1015 278815 : reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 6),
1016 : _mm_cvttps_epi32(xmm_mask6));
1017 278815 : _mm_store_si128(
1018 278815 : reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 7),
1019 : _mm_cvttps_epi32(xmm_mask7));
1020 : }
1021 1804 : for (; iPixel < nPixels; iPixel++)
1022 1560 : reinterpret_cast<int *>(pafMask)[iPixel] =
1023 1560 : static_cast<int>(pafMask[iPixel] * cst_alpha_max);
1024 :
1025 : // Write data.
1026 : // Assumes little endianness here.
1027 488 : eErr = GDALRasterIOEx(
1028 : hAlphaBand, GF_Write, nXOff, nYOff, nXSize, nYSize, pafMask,
1029 : nXSize, nYSize, eDT, static_cast<GSpacing>(sizeof(int)),
1030 244 : static_cast<GSpacing>(sizeof(int)) * nXSize, nullptr);
1031 : }
1032 : else
1033 : #endif
1034 : {
1035 9 : for (; iPixel + 3 < nPixels; iPixel += 4)
1036 : {
1037 8 : pafMask[iPixel + 0] = static_cast<float>(
1038 8 : static_cast<int>(pafMask[iPixel + 0] * cst_alpha_max));
1039 8 : pafMask[iPixel + 1] = static_cast<float>(
1040 8 : static_cast<int>(pafMask[iPixel + 1] * cst_alpha_max));
1041 8 : pafMask[iPixel + 2] = static_cast<float>(
1042 8 : static_cast<int>(pafMask[iPixel + 2] * cst_alpha_max));
1043 8 : pafMask[iPixel + 3] = static_cast<float>(
1044 8 : static_cast<int>(pafMask[iPixel + 3] * cst_alpha_max));
1045 : }
1046 2 : for (; iPixel < nPixels; iPixel++)
1047 1 : pafMask[iPixel] = static_cast<float>(
1048 1 : static_cast<int>(pafMask[iPixel] * cst_alpha_max));
1049 :
1050 : // Write data.
1051 :
1052 : eErr =
1053 1 : GDALRasterIO(hAlphaBand, GF_Write, nXOff, nYOff, nXSize, nYSize,
1054 : pafMask, nXSize, nYSize, GDT_Float32, 0, 0);
1055 : }
1056 245 : return eErr;
1057 : }
1058 : }
1059 :
1060 : /************************************************************************/
1061 : /* ==================================================================== */
1062 : /* GDALWarpOptions */
1063 : /* ==================================================================== */
1064 : /************************************************************************/
1065 :
1066 : /**
1067 : * \var char **GDALWarpOptions::papszWarpOptions;
1068 : *
1069 : * A string list of additional options controlling the warp operation in
1070 : * name=value format. A suitable string list can be prepared with
1071 : * CSLSetNameValue().
1072 : *
1073 : * The following values are currently supported:
1074 : * <ul>
1075 : * <li>INIT_DEST=[value] or INIT_DEST=NO_DATA: This option forces the
1076 : * destination image to be initialized to the indicated value (for all bands)
1077 : * or indicates that it should be initialized to the NO_DATA value in
1078 : * padfDstNoDataReal/padfDstNoDataImag. If this value isn't set the
1079 : * destination image will be read and overlaid.</li>
1080 : *
1081 : * <li>WRITE_FLUSH=YES/NO: This option forces a flush to disk of data after
1082 : * each chunk is processed. In some cases this helps ensure a serial
1083 : * writing of the output data otherwise a block of data may be written to disk
1084 : * each time a block of data is read for the input buffer resulting in a lot
1085 : * of extra seeking around the disk, and reduced IO throughput. The default
1086 : * at this time is NO.</li>
1087 : *
1088 : * <li>SKIP_NOSOURCE=YES/NO: Skip all processing for chunks for which there
1089 : * is no corresponding input data. This will disable initializing the
1090 : * destination (INIT_DEST) and all other processing, and so should be used
1091 : * carefully. Mostly useful to short circuit a lot of extra work in mosaicing
1092 : * situations. Starting with GDAL 2.4, gdalwarp will automatically enable this
1093 : * option when it is assumed to be safe to do so.</li>
1094 : *
1095 : * <li>UNIFIED_SRC_NODATA=YES/NO/PARTIAL: This setting determines
1096 : * how to take into account nodata values when there are several input bands.
1097 : * <ul>
1098 : * <li>When YES, all bands are considered as nodata if and only if, all bands
1099 : * match the corresponding nodata values.
1100 : * Note: UNIFIED_SRC_NODATA=YES is set by default, when called from gdalwarp
1101 : * / GDALWarp() with an explicit -srcnodata setting.
1102 : *
1103 : * Example with nodata values at (1, 2, 3) and target alpha band requested.
1104 : * <ul>
1105 : * <li>input pixel = (1, 2, 3) ==> output pixel = (0, 0, 0, 0)</li>
1106 : * <li>input pixel = (1, 2, 127) ==> output pixel = (1, 2, 127, 255)</li>
1107 : * </ul>
1108 : * </li>
1109 : * <li>When NO, nodata masking values is considered independently for each band.
1110 : * A potential target alpha band will always be valid if there are multiple
1111 : * bands.
1112 : *
1113 : * Example with nodata values at (1, 2, 3) and target alpha band requested.
1114 : * <ul>
1115 : * <li>input pixel = (1, 2, 3) ==> output pixel = (0, 0, 0, 255)</li>
1116 : * <li>input pixel = (1, 2, 127) ==> output pixel = (0, 0, 127, 255)</li>
1117 : * </ul>
1118 : *
1119 : * Note: NO was the default behavior before GDAL 3.3.2
1120 : * </li>
1121 : * <li>When PARTIAL, or not specified at all (default behavior),
1122 : * nodata masking values is considered independently for each band.
1123 : * But, and this is the difference with NO, if for a given pixel, it
1124 : * evaluates to the nodata value of each band, the target pixel is
1125 : * considered as globally invalid, which impacts the value of a potential
1126 : * target alpha band.
1127 : *
1128 : * Note: PARTIAL is new to GDAL 3.3.2 and should not be used with
1129 : * earlier versions. The default behavior of GDAL < 3.3.2 was NO.
1130 : *
1131 : * Example with nodata values at (1, 2, 3) and target alpha band requested.
1132 : * <ul>
1133 : * <li>input pixel = (1, 2, 3) ==> output pixel = (0, 0, 0, 0)</li>
1134 : * <li>input pixel = (1, 2, 127) ==> output pixel = (0, 0, 127, 255)</li>
1135 : * </ul>
1136 : * </li>
1137 : * </ul>
1138 : * </li>
1139 : *
1140 : * <li>CUTLINE: This may contain the WKT geometry for a cutline. It will
1141 : * be converted into a geometry by GDALWarpOperation::Initialize() and assigned
1142 : * to the GDALWarpOptions hCutline field. The coordinates must be expressed
1143 : * in source pixel/line coordinates. Note: this is different from the
1144 : * assumptions made for the -cutline option of the gdalwarp utility !</li>
1145 : *
1146 : * <li>CUTLINE_BLEND_DIST: This may be set with a distance in pixels which
1147 : * will be assigned to the dfCutlineBlendDist field in the GDALWarpOptions.</li>
1148 : *
1149 : * <li>CUTLINE_ALL_TOUCHED: This defaults to FALSE, but may be set to TRUE
1150 : * to enable ALL_TOUCHEd mode when rasterizing cutline polygons. This is
1151 : * useful to ensure that that all pixels overlapping the cutline polygon
1152 : * will be selected, not just those whose center point falls within the
1153 : * polygon.</li>
1154 : *
1155 : * <li>XSCALE: Ratio expressing the resampling factor (number of destination
1156 : * pixels per source pixel) along the target horizontal axis.
1157 : * The scale is used to determine the number of source pixels along the x-axis
1158 : * that are considered by the resampling algorithm.
1159 : * Equals to one for no resampling, below one for downsampling
1160 : * and above one for upsampling. This is automatically computed, for each
1161 : * processing chunk, and may thus vary among them, depending on the
1162 : * shape of output regions vs input regions. Such variations can be undesired
1163 : * in some situations. If the resampling factor can be considered as constant
1164 : * over the warped area, setting a constant value can lead to more reproducible
1165 : * pixel output.</li>
1166 : *
1167 : * <li>YSCALE: Same as XSCALE, but along the horizontal axis.</li>
1168 : *
1169 : * <li>OPTIMIZE_SIZE: This defaults to FALSE, but may be set to TRUE
1170 : * typically when writing to a compressed dataset (GeoTIFF with
1171 : * COMPRESS creation option set for example) for achieving a smaller
1172 : * file size. This is achieved by writing at once data aligned on full
1173 : * blocks of the target dataset, which avoids partial writes of
1174 : * compressed blocks and lost space when they are rewritten at the end
1175 : * of the file. However sticking to target block size may cause major
1176 : * processing slowdown for some particular reprojections. Starting
1177 : * with GDAL 3.8, OPTIMIZE_SIZE mode is automatically enabled when it is safe
1178 : * to do so.
1179 : * As this parameter influences the shape of warping chunk, and by default the
1180 : * XSCALE and YSCALE parameters are computed per warping chunk, this parameter may
1181 : * influence the pixel output.
1182 : * </li>
1183 : *
1184 : * <li>NUM_THREADS: (GDAL >= 1.10) Can be set to a numeric value or ALL_CPUS to
1185 : * set the number of threads to use to parallelize the computation part of the
1186 : * warping. If not set, computation will be done in a single thread.</li>
1187 : *
1188 : * <li>STREAMABLE_OUTPUT: (GDAL >= 2.0) This defaults to FALSE, but may
1189 : * be set to TRUE typically when writing to a streamed file. The
1190 : * gdalwarp utility automatically sets this option when writing to
1191 : * /vsistdout/ or a named pipe (on Unix). This option has performance
1192 : * impacts for some reprojections. Note: band interleaved output is
1193 : * not currently supported by the warping algorithm in a streamable
1194 : * compatible way.</li>
1195 : *
1196 : * <li>SRC_COORD_PRECISION: (GDAL >= 2.0). Advanced setting. This
1197 : * defaults to 0, to indicate that no rounding of computing source
1198 : * image coordinates corresponding to the target image must be
1199 : * done. If greater than 0 (and typically below 1), this value,
1200 : * expressed in pixel, will be used to round computed source image
1201 : * coordinates. The purpose of this option is to make the results of
1202 : * warping with the approximated transformer more reproducible and not
1203 : * sensitive to changes in warping memory size. To achieve that,
1204 : * SRC_COORD_PRECISION must be at least 10 times greater than the
1205 : * error threshold. The higher the SRC_COORD_PRECISION/error_threshold
1206 : * ratio, the higher the performance will be, since exact
1207 : * reprojections must statistically be done with a frequency of
1208 : * 4*error_threshold/SRC_COORD_PRECISION.</li>
1209 : *
1210 : * <li>SRC_ALPHA_MAX: (GDAL >= 2.2). Maximum value for the alpha band of the
1211 : * source dataset. If the value is not set and the alpha band has a NBITS
1212 : * metadata item, it is used to set SRC_ALPHA_MAX = 2^NBITS-1. Otherwise, if the
1213 : * value is not set and the alpha band is of type UInt16 (resp Int16), 65535
1214 : * (resp 32767) is used. Otherwise, 255 is used.</li>
1215 : *
1216 : * <li>DST_ALPHA_MAX: (GDAL >= 2.2). Maximum value for the alpha band of the
1217 : * destination dataset. If the value is not set and the alpha band has a NBITS
1218 : * metadata item, it is used to set DST_ALPHA_MAX = 2^NBITS-1. Otherwise, if the
1219 : * value is not set and the alpha band is of type UInt16 (resp Int16), 65535
1220 : * (resp 32767) is used. Otherwise, 255 is used.</li>
1221 : * </ul>
1222 : *
1223 : * Normally when computing the source raster data to
1224 : * load to generate a particular output area, the warper samples transforms
1225 : * 21 points along each edge of the destination region back onto the source
1226 : * file, and uses this to compute a bounding window on the source image that
1227 : * is sufficient. Depending on the transformation in effect, the source
1228 : * window may be a bit too small, or even missing large areas. Problem
1229 : * situations are those where the transformation is very non-linear or
1230 : * "inside out". Examples are transforming from WGS84 to Polar Stereographic
1231 : * for areas around the pole, or transformations where some of the image is
1232 : * untransformable. The following options provide some additional control
1233 : * to deal with errors in computing the source window:
1234 : * <ul>
1235 : *
1236 : * <li>SAMPLE_GRID=YES/NO: Setting this option to YES will force the sampling to
1237 : * include internal points as well as edge points which can be important if
1238 : * the transformation is esoteric inside out, or if large sections of the
1239 : * destination image are not transformable into the source coordinate
1240 : * system.</li>
1241 : *
1242 : * <li>SAMPLE_STEPS: Modifies the density of the sampling grid. The default
1243 : * number of steps is 21. Increasing this can increase the computational
1244 : * cost, but improves the accuracy with which the source region is
1245 : * computed.
1246 : * Starting with GDAL 3.7, this can be set to ALL to mean to sample
1247 : * along all edge points of the destination region (if SAMPLE_GRID=NO or not
1248 : * specified), or all points of the destination region if SAMPLE_GRID=YES.</li>
1249 : *
1250 : * <li>SOURCE_EXTRA: This is a number of extra pixels added around the source
1251 : * window for a given request, and by default it is 1 to take care of rounding
1252 : * error. Setting this larger will increase the amount of data that needs to
1253 : * be read, but can avoid missing source data.</li>
1254 : * <li>APPLY_VERTICAL_SHIFT=YES/NO: Force the use of vertical shift.
1255 : * This option is generally not necessary, except when using an explicit
1256 : * coordinate transformation (COORDINATE_OPERATION), and not specifying
1257 : * an explicit source and target SRS.</li>
1258 : * <li>MULT_FACTOR_VERTICAL_SHIFT: Multiplication factor for the vertical
1259 : * shift. Default 1.0</li>
1260 : *
1261 : * <li>EXCLUDED_VALUES: (GDAL >= 3.9) Comma-separated tuple of values
1262 : * (thus typically "R,G,B"), that are ignored as contributing source
1263 : * pixels during resampling. The number of values in the tuple must be the same
1264 : * as the number of bands, excluding the alpha band.
1265 : * Several tuples of excluded values may be specified using the
1266 : * "(R1,G1,B2),(R2,G2,B2)" syntax.
1267 : * Only taken into account by Average currently.
1268 : * This concept is a bit similar to nodata/alpha, but the main difference is
1269 : * that pixels matching one of the excluded value tuples are still considered
1270 : * as valid, when determining the target pixel validity/density.
1271 : * </li>
1272 : *
1273 : * <li>EXCLUDED_VALUES_PCT_THRESHOLD=[0-100]: (GDAL >= 3.9) Minimum percentage
1274 : * of source pixels that must be set at one of the EXCLUDED_VALUES to cause
1275 : * the excluded value, that is in majority among source pixels, to be used as the
1276 : * target pixel value. Default value is 50 (%).
1277 : * Only taken into account by Average currently.</li>
1278 : *
1279 : * <li>NODATA_VALUES_PCT_THRESHOLD=[0-100]: (GDAL >= 3.9) Minimum percentage
1280 : * of source pixels that must be at nodata (or alpha=0 or any other way to express
1281 : * transparent pixel) to cause the target pixel value to not be set. Default
1282 : * value is 100 (%), which means that a target pixel is not set only if all
1283 : * contributing source pixels are not set.
1284 : * Note that NODATA_VALUES_PCT_THRESHOLD is taken into account before
1285 : * EXCLUDED_VALUES_PCT_THRESHOLD.
1286 : * Only taken into account by Average currently.</li>
1287 : *
1288 : * <li>MODE_TIES=FIRST/MIN/MAX: (GDAL >= 3.11) Strategy to use when breaking
1289 : * ties with MODE resampling. By default, the first value encountered will be used.
1290 : * Alternatively, the minimum or maximum value can be selected.</li>
1291 : *
1292 : * </ul>
1293 : */
1294 :
1295 : /************************************************************************/
1296 : /* GDALCreateWarpOptions() */
1297 : /************************************************************************/
1298 :
1299 : /** Create a warp options structure.
1300 : *
1301 : * Must be deallocated with GDALDestroyWarpOptions()
1302 : */
1303 3011 : GDALWarpOptions *CPL_STDCALL GDALCreateWarpOptions()
1304 :
1305 : {
1306 : GDALWarpOptions *psOptions =
1307 3011 : static_cast<GDALWarpOptions *>(CPLCalloc(sizeof(GDALWarpOptions), 1));
1308 :
1309 3011 : psOptions->nBandCount = 0;
1310 3011 : psOptions->eResampleAlg = GRA_NearestNeighbour;
1311 3011 : psOptions->pfnProgress = GDALDummyProgress;
1312 3011 : psOptions->eWorkingDataType = GDT_Unknown;
1313 3011 : psOptions->eTieStrategy = GWKTS_First;
1314 :
1315 3011 : return psOptions;
1316 : }
1317 :
1318 : /************************************************************************/
1319 : /* GDALDestroyWarpOptions() */
1320 : /************************************************************************/
1321 :
1322 : /** Destroy a warp options structure. */
1323 3011 : void CPL_STDCALL GDALDestroyWarpOptions(GDALWarpOptions *psOptions)
1324 :
1325 : {
1326 3011 : if (psOptions == nullptr)
1327 0 : return;
1328 :
1329 3011 : CSLDestroy(psOptions->papszWarpOptions);
1330 3011 : CPLFree(psOptions->panSrcBands);
1331 3011 : CPLFree(psOptions->panDstBands);
1332 3011 : CPLFree(psOptions->padfSrcNoDataReal);
1333 3011 : CPLFree(psOptions->padfSrcNoDataImag);
1334 3011 : CPLFree(psOptions->padfDstNoDataReal);
1335 3011 : CPLFree(psOptions->padfDstNoDataImag);
1336 3011 : CPLFree(psOptions->papfnSrcPerBandValidityMaskFunc);
1337 3011 : CPLFree(psOptions->papSrcPerBandValidityMaskFuncArg);
1338 :
1339 3011 : if (psOptions->hCutline != nullptr)
1340 48 : delete static_cast<OGRGeometry *>(psOptions->hCutline);
1341 :
1342 3011 : CPLFree(psOptions);
1343 : }
1344 :
1345 : #define COPY_MEM(target, type, count) \
1346 : do \
1347 : { \
1348 : if ((psSrcOptions->target) != nullptr && (count) != 0) \
1349 : { \
1350 : (psDstOptions->target) = \
1351 : static_cast<type *>(CPLMalloc(sizeof(type) * (count))); \
1352 : memcpy((psDstOptions->target), (psSrcOptions->target), \
1353 : sizeof(type) * (count)); \
1354 : } \
1355 : else \
1356 : (psDstOptions->target) = nullptr; \
1357 : } while (false)
1358 :
1359 : /************************************************************************/
1360 : /* GDALCloneWarpOptions() */
1361 : /************************************************************************/
1362 :
1363 : /** Clone a warp options structure.
1364 : *
1365 : * Must be deallocated with GDALDestroyWarpOptions()
1366 : */
1367 : GDALWarpOptions *CPL_STDCALL
1368 1599 : GDALCloneWarpOptions(const GDALWarpOptions *psSrcOptions)
1369 :
1370 : {
1371 1599 : GDALWarpOptions *psDstOptions = GDALCreateWarpOptions();
1372 :
1373 1599 : memcpy(psDstOptions, psSrcOptions, sizeof(GDALWarpOptions));
1374 :
1375 1599 : if (psSrcOptions->papszWarpOptions != nullptr)
1376 1338 : psDstOptions->papszWarpOptions =
1377 1338 : CSLDuplicate(psSrcOptions->papszWarpOptions);
1378 :
1379 1599 : COPY_MEM(panSrcBands, int, psSrcOptions->nBandCount);
1380 1599 : COPY_MEM(panDstBands, int, psSrcOptions->nBandCount);
1381 1599 : COPY_MEM(padfSrcNoDataReal, double, psSrcOptions->nBandCount);
1382 1599 : COPY_MEM(padfSrcNoDataImag, double, psSrcOptions->nBandCount);
1383 1599 : COPY_MEM(padfDstNoDataReal, double, psSrcOptions->nBandCount);
1384 1599 : COPY_MEM(padfDstNoDataImag, double, psSrcOptions->nBandCount);
1385 : // cppcheck-suppress pointerSize
1386 1599 : COPY_MEM(papfnSrcPerBandValidityMaskFunc, GDALMaskFunc,
1387 : psSrcOptions->nBandCount);
1388 1599 : psDstOptions->papSrcPerBandValidityMaskFuncArg = nullptr;
1389 :
1390 1599 : if (psSrcOptions->hCutline != nullptr)
1391 10 : psDstOptions->hCutline =
1392 10 : OGR_G_Clone(static_cast<OGRGeometryH>(psSrcOptions->hCutline));
1393 1599 : psDstOptions->dfCutlineBlendDist = psSrcOptions->dfCutlineBlendDist;
1394 :
1395 1599 : return psDstOptions;
1396 : }
1397 :
1398 : namespace
1399 : {
1400 143 : void InitNoData(int nBandCount, double **ppdNoDataReal, double dDataReal)
1401 : {
1402 143 : if (nBandCount <= 0)
1403 : {
1404 0 : return;
1405 : }
1406 143 : if (*ppdNoDataReal != nullptr)
1407 : {
1408 32 : return;
1409 : }
1410 :
1411 111 : *ppdNoDataReal =
1412 111 : static_cast<double *>(CPLMalloc(sizeof(double) * nBandCount));
1413 :
1414 254 : for (int i = 0; i < nBandCount; ++i)
1415 : {
1416 143 : (*ppdNoDataReal)[i] = dDataReal;
1417 : }
1418 : }
1419 : } // namespace
1420 :
1421 : /************************************************************************/
1422 : /* GDALWarpInitDstNoDataReal() */
1423 : /************************************************************************/
1424 :
1425 : /**
1426 : * \brief Initialize padfDstNoDataReal with specified value.
1427 : *
1428 : * @param psOptionsIn options to initialize.
1429 : * @param dNoDataReal value to initialize to.
1430 : *
1431 : */
1432 38 : void CPL_STDCALL GDALWarpInitDstNoDataReal(GDALWarpOptions *psOptionsIn,
1433 : double dNoDataReal)
1434 : {
1435 38 : VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitDstNoDataReal");
1436 38 : InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfDstNoDataReal,
1437 : dNoDataReal);
1438 : }
1439 :
1440 : /************************************************************************/
1441 : /* GDALWarpInitSrcNoDataReal() */
1442 : /************************************************************************/
1443 :
1444 : /**
1445 : * \brief Initialize padfSrcNoDataReal with specified value.
1446 : *
1447 : * @param psOptionsIn options to initialize.
1448 : * @param dNoDataReal value to initialize to.
1449 : *
1450 : */
1451 39 : void CPL_STDCALL GDALWarpInitSrcNoDataReal(GDALWarpOptions *psOptionsIn,
1452 : double dNoDataReal)
1453 : {
1454 39 : VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitSrcNoDataReal");
1455 39 : InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfSrcNoDataReal,
1456 : dNoDataReal);
1457 : }
1458 :
1459 : /************************************************************************/
1460 : /* GDALWarpInitNoDataReal() */
1461 : /************************************************************************/
1462 :
1463 : /**
1464 : * \brief Initialize padfSrcNoDataReal and padfDstNoDataReal with specified
1465 : * value.
1466 : *
1467 : * @param psOptionsIn options to initialize.
1468 : * @param dNoDataReal value to initialize to.
1469 : *
1470 : */
1471 1 : void CPL_STDCALL GDALWarpInitNoDataReal(GDALWarpOptions *psOptionsIn,
1472 : double dNoDataReal)
1473 : {
1474 1 : GDALWarpInitDstNoDataReal(psOptionsIn, dNoDataReal);
1475 1 : GDALWarpInitSrcNoDataReal(psOptionsIn, dNoDataReal);
1476 1 : }
1477 :
1478 : /************************************************************************/
1479 : /* GDALWarpInitDstNoDataImag() */
1480 : /************************************************************************/
1481 :
1482 : /**
1483 : * \brief Initialize padfDstNoDataImag with specified value.
1484 : *
1485 : * @param psOptionsIn options to initialize.
1486 : * @param dNoDataImag value to initialize to.
1487 : *
1488 : */
1489 35 : void CPL_STDCALL GDALWarpInitDstNoDataImag(GDALWarpOptions *psOptionsIn,
1490 : double dNoDataImag)
1491 : {
1492 35 : VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitDstNoDataImag");
1493 35 : InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfDstNoDataImag,
1494 : dNoDataImag);
1495 : }
1496 :
1497 : /************************************************************************/
1498 : /* GDALWarpInitSrcNoDataImag() */
1499 : /************************************************************************/
1500 :
1501 : /**
1502 : * \brief Initialize padfSrcNoDataImag with specified value.
1503 : *
1504 : * @param psOptionsIn options to initialize.
1505 : * @param dNoDataImag value to initialize to.
1506 : *
1507 : */
1508 31 : void CPL_STDCALL GDALWarpInitSrcNoDataImag(GDALWarpOptions *psOptionsIn,
1509 : double dNoDataImag)
1510 : {
1511 31 : VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitSrcNoDataImag");
1512 31 : InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfSrcNoDataImag,
1513 : dNoDataImag);
1514 : }
1515 :
1516 : /************************************************************************/
1517 : /* GDALWarpResolveWorkingDataType() */
1518 : /************************************************************************/
1519 :
1520 : /**
1521 : * \brief If the working data type is unknown, this method will determine
1522 : * a valid working data type to support the data in the src and dest
1523 : * data sets and any noData values.
1524 : *
1525 : * @param psOptions options to initialize.
1526 : *
1527 : */
1528 1495 : void CPL_STDCALL GDALWarpResolveWorkingDataType(GDALWarpOptions *psOptions)
1529 : {
1530 1495 : if (psOptions == nullptr)
1531 : {
1532 0 : return;
1533 : }
1534 : /* -------------------------------------------------------------------- */
1535 : /* If no working data type was provided, set one now. */
1536 : /* */
1537 : /* Ensure that the working data type can encapsulate any value */
1538 : /* in the target, source, and the no data for either. */
1539 : /* -------------------------------------------------------------------- */
1540 1495 : if (psOptions->eWorkingDataType != GDT_Unknown)
1541 : {
1542 326 : return;
1543 : }
1544 :
1545 1169 : psOptions->eWorkingDataType = GDT_Byte;
1546 :
1547 : // If none of the provided input nodata values can be represented in the
1548 : // data type of the corresponding source band, ignore them.
1549 1169 : if (psOptions->hSrcDS && psOptions->padfSrcNoDataReal)
1550 : {
1551 87 : int nCountInvalidSrcNoDataReal = 0;
1552 220 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
1553 : {
1554 266 : GDALRasterBandH hSrcBand = GDALGetRasterBand(
1555 133 : psOptions->hSrcDS, psOptions->panSrcBands[iBand]);
1556 :
1557 266 : if (hSrcBand &&
1558 133 : !GDALIsValueExactAs(psOptions->padfSrcNoDataReal[iBand],
1559 : GDALGetRasterDataType(hSrcBand)))
1560 : {
1561 2 : nCountInvalidSrcNoDataReal++;
1562 : }
1563 : }
1564 87 : if (nCountInvalidSrcNoDataReal == psOptions->nBandCount)
1565 : {
1566 2 : CPLFree(psOptions->padfSrcNoDataReal);
1567 2 : psOptions->padfSrcNoDataReal = nullptr;
1568 2 : CPLFree(psOptions->padfSrcNoDataImag);
1569 2 : psOptions->padfSrcNoDataImag = nullptr;
1570 : }
1571 : }
1572 :
1573 2618 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
1574 : {
1575 1449 : if (psOptions->hDstDS != nullptr)
1576 : {
1577 2798 : GDALRasterBandH hDstBand = GDALGetRasterBand(
1578 1399 : psOptions->hDstDS, psOptions->panDstBands[iBand]);
1579 :
1580 1399 : if (hDstBand != nullptr)
1581 : {
1582 1399 : psOptions->eWorkingDataType =
1583 1399 : GDALDataTypeUnion(psOptions->eWorkingDataType,
1584 : GDALGetRasterDataType(hDstBand));
1585 : }
1586 : }
1587 :
1588 1449 : if (psOptions->hSrcDS != nullptr)
1589 : {
1590 2858 : GDALRasterBandH hSrcBand = GDALGetRasterBand(
1591 1429 : psOptions->hSrcDS, psOptions->panSrcBands[iBand]);
1592 :
1593 1429 : if (hSrcBand != nullptr)
1594 : {
1595 1429 : psOptions->eWorkingDataType =
1596 1429 : GDALDataTypeUnion(psOptions->eWorkingDataType,
1597 : GDALGetRasterDataType(hSrcBand));
1598 : }
1599 : }
1600 :
1601 1449 : if (psOptions->padfSrcNoDataReal != nullptr)
1602 : {
1603 141 : psOptions->eWorkingDataType = GDALDataTypeUnionWithValue(
1604 : psOptions->eWorkingDataType,
1605 141 : psOptions->padfSrcNoDataReal[iBand], false);
1606 : }
1607 :
1608 1449 : if (psOptions->padfSrcNoDataImag != nullptr &&
1609 4 : psOptions->padfSrcNoDataImag[iBand] != 0.0)
1610 : {
1611 3 : psOptions->eWorkingDataType = GDALDataTypeUnionWithValue(
1612 : psOptions->eWorkingDataType,
1613 3 : psOptions->padfSrcNoDataImag[iBand], true);
1614 : }
1615 :
1616 1449 : if (psOptions->padfDstNoDataReal != nullptr)
1617 : {
1618 132 : psOptions->eWorkingDataType = GDALDataTypeUnionWithValue(
1619 : psOptions->eWorkingDataType,
1620 132 : psOptions->padfDstNoDataReal[iBand], false);
1621 : }
1622 :
1623 1449 : if (psOptions->padfDstNoDataImag != nullptr &&
1624 52 : psOptions->padfDstNoDataImag[iBand] != 0.0)
1625 : {
1626 3 : psOptions->eWorkingDataType = GDALDataTypeUnionWithValue(
1627 : psOptions->eWorkingDataType,
1628 3 : psOptions->padfDstNoDataImag[iBand], true);
1629 : }
1630 : }
1631 :
1632 2338 : const bool bApplyVerticalShift = CPLFetchBool(
1633 1169 : psOptions->papszWarpOptions, "APPLY_VERTICAL_SHIFT", false);
1634 1188 : if (bApplyVerticalShift &&
1635 19 : GDALDataTypeIsInteger(psOptions->eWorkingDataType))
1636 : {
1637 16 : const double dfMultFactorVerticalShift = CPLAtof(CSLFetchNameValueDef(
1638 16 : psOptions->papszWarpOptions, "MULT_FACTOR_VERTICAL_SHIFT", "1.0"));
1639 16 : if (dfMultFactorVerticalShift != 1)
1640 : {
1641 0 : psOptions->eWorkingDataType =
1642 0 : GDALDataTypeUnion(psOptions->eWorkingDataType, GDT_Float32);
1643 : }
1644 : }
1645 : }
1646 :
1647 : /************************************************************************/
1648 : /* GDALWarpInitDefaultBandMapping() */
1649 : /************************************************************************/
1650 :
1651 : /**
1652 : * \brief Init src and dst band mappings such that Bands[i] = i+1
1653 : * for nBandCount
1654 : * Does nothing if psOptionsIn->nBandCount is non-zero.
1655 : *
1656 : * @param psOptionsIn options to initialize.
1657 : * @param nBandCount bands to initialize for.
1658 : *
1659 : */
1660 312 : void CPL_STDCALL GDALWarpInitDefaultBandMapping(GDALWarpOptions *psOptionsIn,
1661 : int nBandCount)
1662 : {
1663 312 : if (psOptionsIn->nBandCount != 0)
1664 : {
1665 0 : return;
1666 : }
1667 :
1668 312 : psOptionsIn->nBandCount = nBandCount;
1669 :
1670 312 : psOptionsIn->panSrcBands =
1671 312 : static_cast<int *>(CPLMalloc(sizeof(int) * psOptionsIn->nBandCount));
1672 312 : psOptionsIn->panDstBands =
1673 312 : static_cast<int *>(CPLMalloc(sizeof(int) * psOptionsIn->nBandCount));
1674 :
1675 820 : for (int i = 0; i < psOptionsIn->nBandCount; i++)
1676 : {
1677 508 : psOptionsIn->panSrcBands[i] = i + 1;
1678 508 : psOptionsIn->panDstBands[i] = i + 1;
1679 : }
1680 : }
1681 :
1682 : /************************************************************************/
1683 : /* GDALSerializeWarpOptions() */
1684 : /************************************************************************/
1685 :
1686 77 : CPLXMLNode *CPL_STDCALL GDALSerializeWarpOptions(const GDALWarpOptions *psWO)
1687 :
1688 : {
1689 : /* -------------------------------------------------------------------- */
1690 : /* Create root. */
1691 : /* -------------------------------------------------------------------- */
1692 : CPLXMLNode *psTree =
1693 77 : CPLCreateXMLNode(nullptr, CXT_Element, "GDALWarpOptions");
1694 :
1695 : /* -------------------------------------------------------------------- */
1696 : /* WarpMemoryLimit */
1697 : /* -------------------------------------------------------------------- */
1698 77 : CPLCreateXMLElementAndValue(
1699 : psTree, "WarpMemoryLimit",
1700 154 : CPLString().Printf("%g", psWO->dfWarpMemoryLimit));
1701 :
1702 : /* -------------------------------------------------------------------- */
1703 : /* ResampleAlg */
1704 : /* -------------------------------------------------------------------- */
1705 77 : const char *pszAlgName = nullptr;
1706 :
1707 77 : if (psWO->eResampleAlg == GRA_NearestNeighbour)
1708 77 : pszAlgName = "NearestNeighbour";
1709 0 : else if (psWO->eResampleAlg == GRA_Bilinear)
1710 0 : pszAlgName = "Bilinear";
1711 0 : else if (psWO->eResampleAlg == GRA_Cubic)
1712 0 : pszAlgName = "Cubic";
1713 0 : else if (psWO->eResampleAlg == GRA_CubicSpline)
1714 0 : pszAlgName = "CubicSpline";
1715 0 : else if (psWO->eResampleAlg == GRA_Lanczos)
1716 0 : pszAlgName = "Lanczos";
1717 0 : else if (psWO->eResampleAlg == GRA_Average)
1718 0 : pszAlgName = "Average";
1719 0 : else if (psWO->eResampleAlg == GRA_RMS)
1720 0 : pszAlgName = "RootMeanSquare";
1721 0 : else if (psWO->eResampleAlg == GRA_Mode)
1722 0 : pszAlgName = "Mode";
1723 0 : else if (psWO->eResampleAlg == GRA_Max)
1724 0 : pszAlgName = "Maximum";
1725 0 : else if (psWO->eResampleAlg == GRA_Min)
1726 0 : pszAlgName = "Minimum";
1727 0 : else if (psWO->eResampleAlg == GRA_Med)
1728 0 : pszAlgName = "Median";
1729 0 : else if (psWO->eResampleAlg == GRA_Q1)
1730 0 : pszAlgName = "Quartile1";
1731 0 : else if (psWO->eResampleAlg == GRA_Q3)
1732 0 : pszAlgName = "Quartile3";
1733 0 : else if (psWO->eResampleAlg == GRA_Sum)
1734 0 : pszAlgName = "Sum";
1735 : else
1736 0 : pszAlgName = "Unknown";
1737 :
1738 77 : CPLCreateXMLElementAndValue(psTree, "ResampleAlg", pszAlgName);
1739 :
1740 : /* -------------------------------------------------------------------- */
1741 : /* Working Data Type */
1742 : /* -------------------------------------------------------------------- */
1743 77 : CPLCreateXMLElementAndValue(psTree, "WorkingDataType",
1744 77 : GDALGetDataTypeName(psWO->eWorkingDataType));
1745 :
1746 : /* -------------------------------------------------------------------- */
1747 : /* Name/value warp options. */
1748 : /* -------------------------------------------------------------------- */
1749 314 : for (int iWO = 0; psWO->papszWarpOptions != nullptr &&
1750 314 : psWO->papszWarpOptions[iWO] != nullptr;
1751 : iWO++)
1752 : {
1753 237 : char *pszName = nullptr;
1754 : const char *pszValue =
1755 237 : CPLParseNameValue(psWO->papszWarpOptions[iWO], &pszName);
1756 :
1757 : // EXTRA_ELTS is an internal detail that we will recover
1758 : // no need to serialize it.
1759 : // And CUTLINE is also serialized in a special way
1760 237 : if (pszName != nullptr && !EQUAL(pszName, "EXTRA_ELTS") &&
1761 160 : !EQUAL(pszName, "CUTLINE"))
1762 : {
1763 : CPLXMLNode *psOption =
1764 160 : CPLCreateXMLElementAndValue(psTree, "Option", pszValue);
1765 :
1766 160 : CPLCreateXMLNode(CPLCreateXMLNode(psOption, CXT_Attribute, "name"),
1767 : CXT_Text, pszName);
1768 : }
1769 :
1770 237 : CPLFree(pszName);
1771 : }
1772 :
1773 : /* -------------------------------------------------------------------- */
1774 : /* Source and Destination Data Source */
1775 : /* -------------------------------------------------------------------- */
1776 77 : if (psWO->hSrcDS != nullptr)
1777 : {
1778 77 : CPLCreateXMLElementAndValue(psTree, "SourceDataset",
1779 77 : GDALGetDescription(psWO->hSrcDS));
1780 :
1781 : CSLConstList papszOpenOptions =
1782 77 : GDALDataset::FromHandle(psWO->hSrcDS)->GetOpenOptions();
1783 77 : GDALSerializeOpenOptionsToXML(psTree, papszOpenOptions);
1784 : }
1785 :
1786 154 : if (psWO->hDstDS != nullptr &&
1787 77 : strlen(GDALGetDescription(psWO->hDstDS)) != 0)
1788 : {
1789 0 : CPLCreateXMLElementAndValue(psTree, "DestinationDataset",
1790 0 : GDALGetDescription(psWO->hDstDS));
1791 : }
1792 :
1793 : /* -------------------------------------------------------------------- */
1794 : /* Serialize transformer. */
1795 : /* -------------------------------------------------------------------- */
1796 77 : if (psWO->pfnTransformer != nullptr)
1797 : {
1798 : CPLXMLNode *psTransformerContainer =
1799 77 : CPLCreateXMLNode(psTree, CXT_Element, "Transformer");
1800 :
1801 154 : CPLXMLNode *psTransformerTree = GDALSerializeTransformer(
1802 77 : psWO->pfnTransformer, psWO->pTransformerArg);
1803 :
1804 77 : if (psTransformerTree != nullptr)
1805 77 : CPLAddXMLChild(psTransformerContainer, psTransformerTree);
1806 : }
1807 :
1808 : /* -------------------------------------------------------------------- */
1809 : /* Band count and lists. */
1810 : /* -------------------------------------------------------------------- */
1811 77 : CPLXMLNode *psBandList = nullptr;
1812 :
1813 77 : if (psWO->nBandCount != 0)
1814 77 : psBandList = CPLCreateXMLNode(psTree, CXT_Element, "BandList");
1815 :
1816 240 : for (int i = 0; i < psWO->nBandCount; i++)
1817 : {
1818 : CPLXMLNode *psBand;
1819 :
1820 163 : psBand = CPLCreateXMLNode(psBandList, CXT_Element, "BandMapping");
1821 163 : if (psWO->panSrcBands != nullptr)
1822 163 : CPLCreateXMLNode(CPLCreateXMLNode(psBand, CXT_Attribute, "src"),
1823 : CXT_Text,
1824 326 : CPLString().Printf("%d", psWO->panSrcBands[i]));
1825 163 : if (psWO->panDstBands != nullptr)
1826 163 : CPLCreateXMLNode(CPLCreateXMLNode(psBand, CXT_Attribute, "dst"),
1827 : CXT_Text,
1828 326 : CPLString().Printf("%d", psWO->panDstBands[i]));
1829 :
1830 163 : if (psWO->padfSrcNoDataReal != nullptr)
1831 : {
1832 13 : CPLCreateXMLElementAndValue(
1833 : psBand, "SrcNoDataReal",
1834 13 : VRTSerializeNoData(psWO->padfSrcNoDataReal[i],
1835 13 : psWO->eWorkingDataType, 16)
1836 : .c_str());
1837 : }
1838 :
1839 163 : if (psWO->padfSrcNoDataImag != nullptr)
1840 : {
1841 2 : if (std::isnan(psWO->padfSrcNoDataImag[i]))
1842 0 : CPLCreateXMLElementAndValue(psBand, "SrcNoDataImag", "nan");
1843 : else
1844 2 : CPLCreateXMLElementAndValue(
1845 : psBand, "SrcNoDataImag",
1846 4 : CPLString().Printf("%.16g", psWO->padfSrcNoDataImag[i]));
1847 : }
1848 : // Compatibility with GDAL <= 2.2: if we serialize a SrcNoDataReal,
1849 : // it needs a SrcNoDataImag as well
1850 161 : else if (psWO->padfSrcNoDataReal != nullptr)
1851 : {
1852 11 : CPLCreateXMLElementAndValue(psBand, "SrcNoDataImag", "0");
1853 : }
1854 :
1855 163 : if (psWO->padfDstNoDataReal != nullptr)
1856 : {
1857 13 : CPLCreateXMLElementAndValue(
1858 : psBand, "DstNoDataReal",
1859 13 : VRTSerializeNoData(psWO->padfDstNoDataReal[i],
1860 13 : psWO->eWorkingDataType, 16)
1861 : .c_str());
1862 : }
1863 :
1864 163 : if (psWO->padfDstNoDataImag != nullptr)
1865 : {
1866 2 : if (std::isnan(psWO->padfDstNoDataImag[i]))
1867 0 : CPLCreateXMLElementAndValue(psBand, "DstNoDataImag", "nan");
1868 : else
1869 2 : CPLCreateXMLElementAndValue(
1870 : psBand, "DstNoDataImag",
1871 4 : CPLString().Printf("%.16g", psWO->padfDstNoDataImag[i]));
1872 : }
1873 : // Compatibility with GDAL <= 2.2: if we serialize a DstNoDataReal,
1874 : // it needs a SrcNoDataImag as well
1875 161 : else if (psWO->padfDstNoDataReal != nullptr)
1876 : {
1877 11 : CPLCreateXMLElementAndValue(psBand, "DstNoDataImag", "0");
1878 : }
1879 : }
1880 :
1881 : /* -------------------------------------------------------------------- */
1882 : /* Alpha bands. */
1883 : /* -------------------------------------------------------------------- */
1884 77 : if (psWO->nSrcAlphaBand > 0)
1885 0 : CPLCreateXMLElementAndValue(
1886 : psTree, "SrcAlphaBand",
1887 0 : CPLString().Printf("%d", psWO->nSrcAlphaBand));
1888 :
1889 77 : if (psWO->nDstAlphaBand > 0)
1890 23 : CPLCreateXMLElementAndValue(
1891 : psTree, "DstAlphaBand",
1892 46 : CPLString().Printf("%d", psWO->nDstAlphaBand));
1893 :
1894 : /* -------------------------------------------------------------------- */
1895 : /* Cutline. */
1896 : /* -------------------------------------------------------------------- */
1897 77 : if (psWO->hCutline != nullptr)
1898 : {
1899 0 : char *pszWKT = nullptr;
1900 0 : if (OGR_G_ExportToWkt(static_cast<OGRGeometryH>(psWO->hCutline),
1901 0 : &pszWKT) == OGRERR_NONE)
1902 : {
1903 0 : CPLCreateXMLElementAndValue(psTree, "Cutline", pszWKT);
1904 : }
1905 0 : CPLFree(pszWKT);
1906 : }
1907 :
1908 77 : if (psWO->dfCutlineBlendDist != 0.0)
1909 0 : CPLCreateXMLElementAndValue(
1910 : psTree, "CutlineBlendDist",
1911 0 : CPLString().Printf("%.5g", psWO->dfCutlineBlendDist));
1912 :
1913 77 : return psTree;
1914 : }
1915 :
1916 : /************************************************************************/
1917 : /* GDALDeserializeWarpOptions() */
1918 : /************************************************************************/
1919 :
1920 198 : GDALWarpOptions *CPL_STDCALL GDALDeserializeWarpOptions(CPLXMLNode *psTree)
1921 :
1922 : {
1923 198 : CPLErrorReset();
1924 :
1925 : /* -------------------------------------------------------------------- */
1926 : /* Verify this is the right kind of object. */
1927 : /* -------------------------------------------------------------------- */
1928 198 : if (psTree == nullptr || psTree->eType != CXT_Element ||
1929 198 : !EQUAL(psTree->pszValue, "GDALWarpOptions"))
1930 : {
1931 0 : CPLError(CE_Failure, CPLE_AppDefined,
1932 : "Wrong node, unable to deserialize GDALWarpOptions.");
1933 0 : return nullptr;
1934 : }
1935 :
1936 : /* -------------------------------------------------------------------- */
1937 : /* Create pre-initialized warp options. */
1938 : /* -------------------------------------------------------------------- */
1939 198 : GDALWarpOptions *psWO = GDALCreateWarpOptions();
1940 :
1941 : /* -------------------------------------------------------------------- */
1942 : /* Warp memory limit. */
1943 : /* -------------------------------------------------------------------- */
1944 198 : psWO->dfWarpMemoryLimit =
1945 198 : CPLAtof(CPLGetXMLValue(psTree, "WarpMemoryLimit", "0.0"));
1946 :
1947 : /* -------------------------------------------------------------------- */
1948 : /* resample algorithm */
1949 : /* -------------------------------------------------------------------- */
1950 198 : const char *pszValue = CPLGetXMLValue(psTree, "ResampleAlg", "Default");
1951 :
1952 198 : if (EQUAL(pszValue, "NearestNeighbour"))
1953 140 : psWO->eResampleAlg = GRA_NearestNeighbour;
1954 58 : else if (EQUAL(pszValue, "Bilinear"))
1955 8 : psWO->eResampleAlg = GRA_Bilinear;
1956 50 : else if (EQUAL(pszValue, "Cubic"))
1957 9 : psWO->eResampleAlg = GRA_Cubic;
1958 41 : else if (EQUAL(pszValue, "CubicSpline"))
1959 9 : psWO->eResampleAlg = GRA_CubicSpline;
1960 32 : else if (EQUAL(pszValue, "Lanczos"))
1961 4 : psWO->eResampleAlg = GRA_Lanczos;
1962 28 : else if (EQUAL(pszValue, "Average"))
1963 6 : psWO->eResampleAlg = GRA_Average;
1964 22 : else if (EQUAL(pszValue, "RootMeanSquare"))
1965 5 : psWO->eResampleAlg = GRA_RMS;
1966 17 : else if (EQUAL(pszValue, "Mode"))
1967 4 : psWO->eResampleAlg = GRA_Mode;
1968 13 : else if (EQUAL(pszValue, "Maximum"))
1969 3 : psWO->eResampleAlg = GRA_Max;
1970 10 : else if (EQUAL(pszValue, "Minimum"))
1971 2 : psWO->eResampleAlg = GRA_Min;
1972 8 : else if (EQUAL(pszValue, "Median"))
1973 3 : psWO->eResampleAlg = GRA_Med;
1974 5 : else if (EQUAL(pszValue, "Quartile1"))
1975 2 : psWO->eResampleAlg = GRA_Q1;
1976 3 : else if (EQUAL(pszValue, "Quartile3"))
1977 2 : psWO->eResampleAlg = GRA_Q3;
1978 1 : else if (EQUAL(pszValue, "Sum"))
1979 1 : psWO->eResampleAlg = GRA_Sum;
1980 0 : else if (EQUAL(pszValue, "Default"))
1981 : /* leave as is */;
1982 : else
1983 : {
1984 0 : CPLError(CE_Failure, CPLE_AppDefined,
1985 : "Unrecognised ResampleAlg value '%s'.", pszValue);
1986 : }
1987 :
1988 : /* -------------------------------------------------------------------- */
1989 : /* Working data type. */
1990 : /* -------------------------------------------------------------------- */
1991 198 : psWO->eWorkingDataType = GDALGetDataTypeByName(
1992 : CPLGetXMLValue(psTree, "WorkingDataType", "Unknown"));
1993 :
1994 : /* -------------------------------------------------------------------- */
1995 : /* Name/value warp options. */
1996 : /* -------------------------------------------------------------------- */
1997 1849 : for (CPLXMLNode *psItem = psTree->psChild; psItem != nullptr;
1998 1651 : psItem = psItem->psNext)
1999 : {
2000 1651 : if (psItem->eType == CXT_Element && EQUAL(psItem->pszValue, "Option"))
2001 : {
2002 346 : const char *pszName = CPLGetXMLValue(psItem, "Name", nullptr);
2003 346 : pszValue = CPLGetXMLValue(psItem, "", nullptr);
2004 :
2005 346 : if (pszName != nullptr && pszValue != nullptr)
2006 : {
2007 346 : psWO->papszWarpOptions =
2008 346 : CSLSetNameValue(psWO->papszWarpOptions, pszName, pszValue);
2009 : }
2010 : }
2011 : }
2012 :
2013 : /* -------------------------------------------------------------------- */
2014 : /* Source Dataset. */
2015 : /* -------------------------------------------------------------------- */
2016 198 : pszValue = CPLGetXMLValue(psTree, "SourceDataset", nullptr);
2017 :
2018 198 : if (pszValue != nullptr)
2019 : {
2020 : CPLXMLNode *psGeoLocNode =
2021 198 : CPLSearchXMLNode(psTree, "GeoLocTransformer");
2022 198 : if (psGeoLocNode)
2023 : {
2024 1 : CPLCreateXMLElementAndValue(psGeoLocNode, "SourceDataset",
2025 : pszValue);
2026 : }
2027 :
2028 396 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
2029 :
2030 198 : char **papszOpenOptions = GDALDeserializeOpenOptionsFromXML(psTree);
2031 198 : psWO->hSrcDS =
2032 198 : GDALOpenEx(pszValue, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
2033 : nullptr, papszOpenOptions, nullptr);
2034 198 : CSLDestroy(papszOpenOptions);
2035 : }
2036 :
2037 : /* -------------------------------------------------------------------- */
2038 : /* Destination Dataset. */
2039 : /* -------------------------------------------------------------------- */
2040 198 : pszValue = CPLGetXMLValue(psTree, "DestinationDataset", nullptr);
2041 :
2042 198 : if (pszValue != nullptr)
2043 : {
2044 0 : psWO->hDstDS = GDALOpenShared(pszValue, GA_Update);
2045 : }
2046 :
2047 : /* -------------------------------------------------------------------- */
2048 : /* First, count band mappings so we can establish the bandcount. */
2049 : /* -------------------------------------------------------------------- */
2050 198 : CPLXMLNode *psBandTree = CPLGetXMLNode(psTree, "BandList");
2051 :
2052 198 : int nBandCount = 0;
2053 198 : CPLXMLNode *psBand = psBandTree ? psBandTree->psChild : nullptr;
2054 563 : for (; psBand != nullptr; psBand = psBand->psNext)
2055 : {
2056 365 : if (psBand->eType != CXT_Element ||
2057 365 : !EQUAL(psBand->pszValue, "BandMapping"))
2058 0 : continue;
2059 :
2060 365 : nBandCount++;
2061 : }
2062 :
2063 198 : GDALWarpInitDefaultBandMapping(psWO, nBandCount);
2064 :
2065 : /* ==================================================================== */
2066 : /* Now actually process each bandmapping. */
2067 : /* ==================================================================== */
2068 198 : int iBand = 0;
2069 :
2070 198 : psBand = psBandTree ? psBandTree->psChild : nullptr;
2071 :
2072 563 : for (; psBand != nullptr; psBand = psBand->psNext)
2073 : {
2074 365 : if (psBand->eType != CXT_Element ||
2075 365 : !EQUAL(psBand->pszValue, "BandMapping"))
2076 0 : continue;
2077 :
2078 : /* --------------------------------------------------------------------
2079 : */
2080 : /* Source band */
2081 : /* --------------------------------------------------------------------
2082 : */
2083 365 : pszValue = CPLGetXMLValue(psBand, "src", nullptr);
2084 365 : if (pszValue != nullptr)
2085 365 : psWO->panSrcBands[iBand] = atoi(pszValue);
2086 :
2087 : /* --------------------------------------------------------------------
2088 : */
2089 : /* Destination band. */
2090 : /* --------------------------------------------------------------------
2091 : */
2092 365 : pszValue = CPLGetXMLValue(psBand, "dst", nullptr);
2093 365 : if (pszValue != nullptr)
2094 365 : psWO->panDstBands[iBand] = atoi(pszValue);
2095 :
2096 66 : const auto NormalizeValue = [](const char *pszValueIn,
2097 : GDALDataType eDataType) -> double
2098 : {
2099 78 : if (eDataType == GDT_Float32 &&
2100 78 : CPLString().Printf(
2101 12 : "%.16g", -std::numeric_limits<float>::max()) == pszValueIn)
2102 : {
2103 0 : return std::numeric_limits<float>::lowest();
2104 : }
2105 78 : else if (eDataType == GDT_Float32 &&
2106 78 : CPLString().Printf("%.16g",
2107 12 : std::numeric_limits<float>::max()) ==
2108 : pszValueIn)
2109 : {
2110 0 : return std::numeric_limits<float>::max();
2111 : }
2112 : else
2113 : {
2114 66 : return CPLAtof(pszValueIn);
2115 : }
2116 : };
2117 :
2118 : /* --------------------------------------------------------------------
2119 : */
2120 : /* Source nodata. */
2121 : /* --------------------------------------------------------------------
2122 : */
2123 365 : pszValue = CPLGetXMLValue(psBand, "SrcNoDataReal", nullptr);
2124 365 : if (pszValue != nullptr)
2125 : {
2126 31 : GDALWarpInitSrcNoDataReal(psWO, -1.1e20);
2127 62 : psWO->padfSrcNoDataReal[iBand] =
2128 31 : NormalizeValue(pszValue, psWO->eWorkingDataType);
2129 : }
2130 :
2131 365 : pszValue = CPLGetXMLValue(psBand, "SrcNoDataImag", nullptr);
2132 365 : if (pszValue != nullptr)
2133 : {
2134 31 : GDALWarpInitSrcNoDataImag(psWO, 0);
2135 31 : psWO->padfSrcNoDataImag[iBand] = CPLAtof(pszValue);
2136 : }
2137 :
2138 : /* --------------------------------------------------------------------
2139 : */
2140 : /* Destination nodata. */
2141 : /* --------------------------------------------------------------------
2142 : */
2143 365 : pszValue = CPLGetXMLValue(psBand, "DstNoDataReal", nullptr);
2144 365 : if (pszValue != nullptr)
2145 : {
2146 35 : GDALWarpInitDstNoDataReal(psWO, -1.1e20);
2147 70 : psWO->padfDstNoDataReal[iBand] =
2148 35 : NormalizeValue(pszValue, psWO->eWorkingDataType);
2149 : }
2150 :
2151 365 : pszValue = CPLGetXMLValue(psBand, "DstNoDataImag", nullptr);
2152 365 : if (pszValue != nullptr)
2153 : {
2154 35 : GDALWarpInitDstNoDataImag(psWO, 0);
2155 35 : psWO->padfDstNoDataImag[iBand] = CPLAtof(pszValue);
2156 : }
2157 :
2158 365 : iBand++;
2159 : }
2160 :
2161 : /* -------------------------------------------------------------------- */
2162 : /* Alpha bands. */
2163 : /* -------------------------------------------------------------------- */
2164 198 : psWO->nSrcAlphaBand = atoi(CPLGetXMLValue(psTree, "SrcAlphaBand", "0"));
2165 198 : psWO->nDstAlphaBand = atoi(CPLGetXMLValue(psTree, "DstAlphaBand", "0"));
2166 :
2167 : /* -------------------------------------------------------------------- */
2168 : /* Cutline. */
2169 : /* -------------------------------------------------------------------- */
2170 198 : const char *pszWKT = CPLGetXMLValue(psTree, "Cutline", nullptr);
2171 198 : if (pszWKT)
2172 : {
2173 7 : char *pszWKTTemp = const_cast<char *>(pszWKT);
2174 7 : OGRGeometryH hCutline = nullptr;
2175 7 : OGR_G_CreateFromWkt(&pszWKTTemp, nullptr, &hCutline);
2176 7 : psWO->hCutline = hCutline;
2177 : }
2178 :
2179 198 : psWO->dfCutlineBlendDist =
2180 198 : CPLAtof(CPLGetXMLValue(psTree, "CutlineBlendDist", "0"));
2181 :
2182 : /* -------------------------------------------------------------------- */
2183 : /* Transformation. */
2184 : /* -------------------------------------------------------------------- */
2185 198 : CPLXMLNode *psTransformer = CPLGetXMLNode(psTree, "Transformer");
2186 :
2187 198 : if (psTransformer != nullptr && psTransformer->psChild != nullptr)
2188 : {
2189 198 : GDALDeserializeTransformer(psTransformer->psChild,
2190 : &(psWO->pfnTransformer),
2191 : &(psWO->pTransformerArg));
2192 : }
2193 :
2194 : /* -------------------------------------------------------------------- */
2195 : /* If any error has occurred, cleanup else return success. */
2196 : /* -------------------------------------------------------------------- */
2197 198 : if (CPLGetLastErrorType() != CE_None)
2198 : {
2199 0 : if (psWO->pTransformerArg)
2200 : {
2201 0 : GDALDestroyTransformer(psWO->pTransformerArg);
2202 0 : psWO->pTransformerArg = nullptr;
2203 : }
2204 0 : if (psWO->hSrcDS != nullptr)
2205 : {
2206 0 : GDALClose(psWO->hSrcDS);
2207 0 : psWO->hSrcDS = nullptr;
2208 : }
2209 0 : if (psWO->hDstDS != nullptr)
2210 : {
2211 0 : GDALClose(psWO->hDstDS);
2212 0 : psWO->hDstDS = nullptr;
2213 : }
2214 0 : GDALDestroyWarpOptions(psWO);
2215 0 : return nullptr;
2216 : }
2217 :
2218 198 : return psWO;
2219 : }
|
