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