Coverage Report

Created: 2025-06-13 06:29

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