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