Coverage Report

Created: 2025-06-09 07:43

/src/gdal/frmts/vrt/vrtprocesseddataset.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  Virtual GDAL Datasets
4
 * Purpose:  Implementation of VRTProcessedDataset.
5
 * Author:   Even Rouault <even.rouault at spatialys.com>
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 2024, Even Rouault <even.rouault at spatialys.com>
9
 *
10
 * SPDX-License-Identifier: MIT
11
 ****************************************************************************/
12
13
#include "cpl_minixml.h"
14
#include "cpl_string.h"
15
#include "gdal_utils.h"
16
#include "vrtdataset.h"
17
18
#include <algorithm>
19
#include <limits>
20
#include <map>
21
#include <vector>
22
23
/************************************************************************/
24
/*                        VRTProcessedDatasetFunc                       */
25
/************************************************************************/
26
27
//! Structure holding information for a VRTProcessedDataset function.
28
struct VRTProcessedDatasetFunc
29
{
30
    //! Processing function name
31
    std::string osFuncName{};
32
33
    //! User data to provide to pfnInit, pfnFree, pfnProcess callbacks.
34
    void *pUserData = nullptr;
35
36
    //! Whether XML metadata has been specified
37
    bool bMetadataSpecified = false;
38
39
    //! Map of (constant argument name, constant value)
40
    std::map<std::string, std::string> oMapConstantArguments{};
41
42
    //! Set of builtin argument names (e.g "offset", "scale", "nodata")
43
    std::set<std::string> oSetBuiltinArguments{};
44
45
    //! Arguments defined in the VRT
46
    struct OtherArgument
47
    {
48
        std::string osType{};
49
        bool bRequired = false;
50
    };
51
52
    std::map<std::string, OtherArgument> oOtherArguments{};
53
54
    //! Requested input data type.
55
    GDALDataType eRequestedInputDT = GDT_Unknown;
56
57
    //! List of supported input datatypes. Empty if no restriction.
58
    std::vector<GDALDataType> aeSupportedInputDT{};
59
60
    //! List of supported input band counts. Empty if no restriction.
61
    std::vector<int> anSupportedInputBandCount{};
62
63
    //! Optional initialization function
64
    GDALVRTProcessedDatasetFuncInit pfnInit = nullptr;
65
66
    //! Optional free function
67
    GDALVRTProcessedDatasetFuncFree pfnFree = nullptr;
68
69
    //! Required processing function
70
    GDALVRTProcessedDatasetFuncProcess pfnProcess = nullptr;
71
};
72
73
/************************************************************************/
74
/*                      GetGlobalMapProcessedDatasetFunc()              */
75
/************************************************************************/
76
77
/** Return the registry of VRTProcessedDatasetFunc functions */
78
static std::map<std::string, VRTProcessedDatasetFunc> &
79
GetGlobalMapProcessedDatasetFunc()
80
10
{
81
10
    static std::map<std::string, VRTProcessedDatasetFunc> goMap;
82
10
    return goMap;
83
10
}
84
85
/************************************************************************/
86
/*                            Step::~Step()                             */
87
/************************************************************************/
88
89
/*! @cond Doxygen_Suppress */
90
91
/** Step destructor */
92
VRTProcessedDataset::Step::~Step()
93
0
{
94
0
    deinit();
95
0
}
96
97
/************************************************************************/
98
/*                           Step::deinit()                             */
99
/************************************************************************/
100
101
/** Free pWorkingData */
102
void VRTProcessedDataset::Step::deinit()
103
0
{
104
0
    if (pWorkingData)
105
0
    {
106
0
        const auto &oMapFunctions = GetGlobalMapProcessedDatasetFunc();
107
0
        const auto oIterFunc = oMapFunctions.find(osAlgorithm);
108
0
        if (oIterFunc != oMapFunctions.end())
109
0
        {
110
0
            if (oIterFunc->second.pfnFree)
111
0
            {
112
0
                oIterFunc->second.pfnFree(osAlgorithm.c_str(),
113
0
                                          oIterFunc->second.pUserData,
114
0
                                          pWorkingData);
115
0
            }
116
0
        }
117
0
        else
118
0
        {
119
0
            CPLAssert(false);
120
0
        }
121
0
        pWorkingData = nullptr;
122
0
    }
123
0
}
124
125
/************************************************************************/
126
/*                        Step::Step(Step&& other)                      */
127
/************************************************************************/
128
129
/** Move constructor */
130
VRTProcessedDataset::Step::Step(Step &&other)
131
0
    : osAlgorithm(std::move(other.osAlgorithm)),
132
0
      aosArguments(std::move(other.aosArguments)), eInDT(other.eInDT),
133
0
      eOutDT(other.eOutDT), nInBands(other.nInBands),
134
0
      nOutBands(other.nOutBands), adfInNoData(other.adfInNoData),
135
0
      adfOutNoData(other.adfOutNoData), pWorkingData(other.pWorkingData)
136
0
{
137
0
    other.pWorkingData = nullptr;
138
0
}
139
140
/************************************************************************/
141
/*                      Step operator=(Step&& other)                    */
142
/************************************************************************/
143
144
/** Move assignment operator */
145
VRTProcessedDataset::Step &VRTProcessedDataset::Step::operator=(Step &&other)
146
0
{
147
0
    if (&other != this)
148
0
    {
149
0
        deinit();
150
0
        osAlgorithm = std::move(other.osAlgorithm);
151
0
        aosArguments = std::move(other.aosArguments);
152
0
        eInDT = other.eInDT;
153
0
        eOutDT = other.eOutDT;
154
0
        nInBands = other.nInBands;
155
0
        nOutBands = other.nOutBands;
156
0
        adfInNoData = std::move(other.adfInNoData);
157
0
        adfOutNoData = std::move(other.adfOutNoData);
158
0
        std::swap(pWorkingData, other.pWorkingData);
159
0
    }
160
0
    return *this;
161
0
}
162
163
/************************************************************************/
164
/*                        VRTProcessedDataset()                         */
165
/************************************************************************/
166
167
/** Constructor */
168
VRTProcessedDataset::VRTProcessedDataset(int nXSize, int nYSize)
169
0
    : VRTDataset(nXSize, nYSize)
170
0
{
171
0
}
172
173
/************************************************************************/
174
/*                       ~VRTProcessedDataset()                         */
175
/************************************************************************/
176
177
VRTProcessedDataset::~VRTProcessedDataset()
178
179
0
{
180
0
    VRTProcessedDataset::FlushCache(true);
181
0
    VRTProcessedDataset::CloseDependentDatasets();
182
0
}
183
184
/************************************************************************/
185
/*                              XMLInit()                               */
186
/************************************************************************/
187
188
/** Instantiate object from XML tree */
189
CPLErr VRTProcessedDataset::XMLInit(const CPLXMLNode *psTree,
190
                                    const char *pszVRTPathIn)
191
192
0
{
193
0
    if (Init(psTree, pszVRTPathIn, nullptr, nullptr, -1) != CE_None)
194
0
        return CE_Failure;
195
196
0
    const auto poSrcFirstBand = m_poSrcDS->GetRasterBand(1);
197
0
    const int nOvrCount = poSrcFirstBand->GetOverviewCount();
198
0
    for (int i = 0; i < nOvrCount; ++i)
199
0
    {
200
0
        auto poOvrDS = std::make_unique<VRTProcessedDataset>(0, 0);
201
0
        if (poOvrDS->Init(psTree, pszVRTPathIn, this, m_poSrcDS.get(), i) !=
202
0
            CE_None)
203
0
            break;
204
0
        m_apoOverviewDatasets.emplace_back(std::move(poOvrDS));
205
0
    }
206
207
0
    return CE_None;
208
0
}
209
210
static bool HasScaleOffset(GDALDataset &oSrcDS)
211
0
{
212
0
    for (int i = 1; i <= oSrcDS.GetRasterCount(); i++)
213
0
    {
214
0
        int pbSuccess;
215
0
        GDALRasterBand &oBand = *oSrcDS.GetRasterBand(i);
216
0
        double scale = oBand.GetScale(&pbSuccess);
217
0
        if (pbSuccess && scale != 1)
218
0
        {
219
0
            return true;
220
0
        }
221
0
        double offset = oBand.GetOffset(&pbSuccess);
222
0
        if (pbSuccess && offset != 0)
223
0
        {
224
0
            return true;
225
0
        }
226
0
    }
227
228
0
    return false;
229
0
}
230
231
/** Instantiate object from XML tree */
232
CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree,
233
                                 const char *pszVRTPathIn,
234
                                 const VRTProcessedDataset *poParentDS,
235
                                 GDALDataset *poParentSrcDS, int iOvrLevel)
236
237
0
{
238
0
    const CPLXMLNode *psInput = CPLGetXMLNode(psTree, "Input");
239
0
    if (!psInput)
240
0
    {
241
0
        CPLError(CE_Failure, CPLE_AppDefined, "Input element missing");
242
0
        return CE_Failure;
243
0
    }
244
245
0
    if (pszVRTPathIn)
246
0
        m_osVRTPath = pszVRTPathIn;
247
248
0
    if (poParentSrcDS)
249
0
    {
250
0
        m_poSrcDS.reset(
251
0
            GDALCreateOverviewDataset(poParentSrcDS, iOvrLevel, true));
252
0
    }
253
0
    else if (const CPLXMLNode *psSourceFileNameNode =
254
0
                 CPLGetXMLNode(psInput, "SourceFilename"))
255
0
    {
256
0
        const bool bRelativeToVRT = CPL_TO_BOOL(
257
0
            atoi(CPLGetXMLValue(psSourceFileNameNode, "relativetoVRT", "0")));
258
0
        const std::string osFilename = GDALDataset::BuildFilename(
259
0
            CPLGetXMLValue(psInput, "SourceFilename", ""), pszVRTPathIn,
260
0
            bRelativeToVRT);
261
0
        m_poSrcDS.reset(GDALDataset::Open(
262
0
            osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr,
263
0
            nullptr, nullptr));
264
0
    }
265
0
    else if (const CPLXMLNode *psVRTDataset =
266
0
                 CPLGetXMLNode(psInput, "VRTDataset"))
267
0
    {
268
0
        CPLXMLNode sVRTDatasetTmp = *psVRTDataset;
269
0
        sVRTDatasetTmp.psNext = nullptr;
270
0
        char *pszXML = CPLSerializeXMLTree(&sVRTDatasetTmp);
271
0
        m_poSrcDS = VRTDataset::OpenXML(pszXML, pszVRTPathIn, GA_ReadOnly);
272
0
        CPLFree(pszXML);
273
0
    }
274
0
    else
275
0
    {
276
0
        CPLError(
277
0
            CE_Failure, CPLE_AppDefined,
278
0
            "Input element should have a SourceFilename or VRTDataset element");
279
0
        return CE_Failure;
280
0
    }
281
282
0
    if (!m_poSrcDS)
283
0
        return CE_Failure;
284
285
0
    const char *pszUnscale = CPLGetXMLValue(psInput, "unscale", "AUTO");
286
0
    bool bUnscale = false;
287
0
    if (EQUAL(pszUnscale, "AUTO"))
288
0
    {
289
0
        if (HasScaleOffset(*m_poSrcDS))
290
0
        {
291
0
            bUnscale = true;
292
0
        }
293
0
    }
294
0
    else if (EQUAL(pszUnscale, "YES") || EQUAL(pszUnscale, "ON") ||
295
0
             EQUAL(pszUnscale, "TRUE") || EQUAL(pszUnscale, "1"))
296
0
    {
297
0
        bUnscale = true;
298
0
    }
299
0
    else if (!(EQUAL(pszUnscale, "NO") || EQUAL(pszUnscale, "OFF") ||
300
0
               EQUAL(pszUnscale, "FALSE") || EQUAL(pszUnscale, "0")))
301
0
    {
302
0
        CPLError(CE_Failure, CPLE_AppDefined, "Invalid value of 'unscale'");
303
0
        return CE_Failure;
304
0
    }
305
306
0
    if (bUnscale)
307
0
    {
308
0
        CPLStringList oArgs;
309
0
        oArgs.AddString("-unscale");
310
0
        oArgs.AddString("-ot");
311
0
        oArgs.AddString("Float64");
312
0
        oArgs.AddString("-of");
313
0
        oArgs.AddString("VRT");
314
0
        oArgs.AddString("-a_nodata");
315
0
        oArgs.AddString("nan");
316
0
        auto *poArgs = GDALTranslateOptionsNew(oArgs.List(), nullptr);
317
0
        int pbUsageError;
318
0
        CPLAssert(poArgs);
319
0
        m_poVRTSrcDS.reset(m_poSrcDS.release());
320
        // https://trac.cppcheck.net/ticket/11325
321
        // cppcheck-suppress accessMoved
322
0
        m_poSrcDS.reset(GDALDataset::FromHandle(
323
0
            GDALTranslate("", m_poVRTSrcDS.get(), poArgs, &pbUsageError)));
324
0
        GDALTranslateOptionsFree(poArgs);
325
326
0
        if (pbUsageError || !m_poSrcDS)
327
0
        {
328
0
            return CE_Failure;
329
0
        }
330
0
    }
331
332
0
    if (nRasterXSize == 0 && nRasterYSize == 0)
333
0
    {
334
0
        nRasterXSize = m_poSrcDS->GetRasterXSize();
335
0
        nRasterYSize = m_poSrcDS->GetRasterYSize();
336
0
    }
337
0
    else if (nRasterXSize != m_poSrcDS->GetRasterXSize() ||
338
0
             nRasterYSize != m_poSrcDS->GetRasterYSize())
339
0
    {
340
0
        CPLError(CE_Failure, CPLE_AppDefined,
341
0
                 "Inconsistent declared VRT dimensions with input dataset");
342
0
        return CE_Failure;
343
0
    }
344
345
0
    if (m_poSrcDS->GetRasterCount() == 0)
346
0
        return CE_Failure;
347
348
    // Inherit SRS from source if not explicitly defined in VRT
349
0
    if (!CPLGetXMLNode(psTree, "SRS"))
350
0
    {
351
0
        const OGRSpatialReference *poSRS = m_poSrcDS->GetSpatialRef();
352
0
        if (poSRS)
353
0
        {
354
0
            m_poSRS.reset(poSRS->Clone());
355
0
        }
356
0
    }
357
358
    // Inherit GeoTransform from source if not explicitly defined in VRT
359
0
    if (iOvrLevel < 0 && !CPLGetXMLNode(psTree, "GeoTransform"))
360
0
    {
361
0
        if (m_poSrcDS->GetGeoTransform(m_adfGeoTransform) == CE_None)
362
0
            m_bGeoTransformSet = true;
363
0
    }
364
365
    /* -------------------------------------------------------------------- */
366
    /*      Initialize blocksize before calling sub-init so that the        */
367
    /*      band initializers can get it from the dataset object when       */
368
    /*      they are created.                                               */
369
    /* -------------------------------------------------------------------- */
370
371
0
    const auto poSrcFirstBand = m_poSrcDS->GetRasterBand(1);
372
0
    poSrcFirstBand->GetBlockSize(&m_nBlockXSize, &m_nBlockYSize);
373
0
    bool bUserBlockSize = false;
374
0
    if (const char *pszBlockXSize =
375
0
            CPLGetXMLValue(psTree, "BlockXSize", nullptr))
376
0
    {
377
0
        bUserBlockSize = true;
378
0
        m_nBlockXSize = atoi(pszBlockXSize);
379
0
        if (m_nBlockXSize <= 1)
380
0
        {
381
0
            CPLError(CE_Failure, CPLE_AppDefined, "Invalid BlockXSize");
382
0
            return CE_Failure;
383
0
        }
384
0
    }
385
0
    if (const char *pszBlockYSize =
386
0
            CPLGetXMLValue(psTree, "BlockYSize", nullptr))
387
0
    {
388
0
        bUserBlockSize = true;
389
0
        m_nBlockYSize = atoi(pszBlockYSize);
390
0
        if (m_nBlockYSize <= 1)
391
0
        {
392
0
            CPLError(CE_Failure, CPLE_AppDefined, "Invalid BlockYSize");
393
0
            return CE_Failure;
394
0
        }
395
0
    }
396
397
    // Initialize all the general VRT stuff.
398
0
    if (VRTDataset::XMLInit(psTree, pszVRTPathIn) != CE_None)
399
0
    {
400
0
        return CE_Failure;
401
0
    }
402
403
    // Use geotransform from parent for overviews
404
0
    if (iOvrLevel >= 0 && poParentDS->m_bGeoTransformSet)
405
0
    {
406
0
        m_bGeoTransformSet = true;
407
0
        m_adfGeoTransform[0] = poParentDS->m_adfGeoTransform[0];
408
0
        m_adfGeoTransform[1] = poParentDS->m_adfGeoTransform[1];
409
0
        m_adfGeoTransform[2] = poParentDS->m_adfGeoTransform[2];
410
0
        m_adfGeoTransform[3] = poParentDS->m_adfGeoTransform[3];
411
0
        m_adfGeoTransform[4] = poParentDS->m_adfGeoTransform[4];
412
0
        m_adfGeoTransform[5] = poParentDS->m_adfGeoTransform[5];
413
414
0
        m_adfGeoTransform[1] *=
415
0
            static_cast<double>(poParentDS->GetRasterXSize()) / nRasterXSize;
416
0
        m_adfGeoTransform[2] *=
417
0
            static_cast<double>(poParentDS->GetRasterYSize()) / nRasterYSize;
418
0
        m_adfGeoTransform[4] *=
419
0
            static_cast<double>(poParentDS->GetRasterXSize()) / nRasterXSize;
420
0
        m_adfGeoTransform[5] *=
421
0
            static_cast<double>(poParentDS->GetRasterYSize()) / nRasterYSize;
422
0
    }
423
424
0
    const CPLXMLNode *psOutputBands = CPLGetXMLNode(psTree, "OutputBands");
425
0
    if (psOutputBands)
426
0
    {
427
0
        if (const char *pszCount =
428
0
                CPLGetXMLValue(psOutputBands, "count", nullptr))
429
0
        {
430
0
            if (EQUAL(pszCount, "FROM_LAST_STEP"))
431
0
            {
432
0
                m_outputBandCountProvenance = ValueProvenance::FROM_LAST_STEP;
433
0
            }
434
0
            else if (!EQUAL(pszCount, "FROM_SOURCE"))
435
0
            {
436
0
                if (CPLGetValueType(pszCount) == CPL_VALUE_INTEGER)
437
0
                {
438
0
                    m_outputBandCountProvenance =
439
0
                        ValueProvenance::USER_PROVIDED;
440
0
                    m_outputBandCountValue = atoi(pszCount);
441
0
                    if (!GDALCheckBandCount(m_outputBandCountValue,
442
0
                                            /* bIsZeroAllowed = */ false))
443
0
                    {
444
                        // Error emitted by above function
445
0
                        return CE_Failure;
446
0
                    }
447
0
                }
448
0
                else
449
0
                {
450
0
                    CPLError(CE_Failure, CPLE_AppDefined,
451
0
                             "Invalid value for OutputBands.count");
452
0
                    return CE_Failure;
453
0
                }
454
0
            }
455
0
        }
456
457
0
        if (const char *pszType =
458
0
                CPLGetXMLValue(psOutputBands, "dataType", nullptr))
459
0
        {
460
0
            if (EQUAL(pszType, "FROM_LAST_STEP"))
461
0
            {
462
0
                m_outputBandDataTypeProvenance =
463
0
                    ValueProvenance::FROM_LAST_STEP;
464
0
            }
465
0
            else if (!EQUAL(pszType, "FROM_SOURCE"))
466
0
            {
467
0
                m_outputBandDataTypeProvenance = ValueProvenance::USER_PROVIDED;
468
0
                m_outputBandDataTypeValue = GDALGetDataTypeByName(pszType);
469
0
                if (m_outputBandDataTypeValue == GDT_Unknown)
470
0
                {
471
0
                    CPLError(CE_Failure, CPLE_AppDefined,
472
0
                             "Invalid value for OutputBands.dataType");
473
0
                    return CE_Failure;
474
0
                }
475
0
            }
476
0
        }
477
0
    }
478
0
    else if (CPLGetXMLNode(psTree, "VRTRasterBand"))
479
0
    {
480
0
        m_outputBandCountProvenance = ValueProvenance::FROM_VRTRASTERBAND;
481
0
        m_outputBandDataTypeProvenance = ValueProvenance::FROM_VRTRASTERBAND;
482
0
    }
483
484
0
    int nOutputBandCount = 0;
485
0
    switch (m_outputBandCountProvenance)
486
0
    {
487
0
        case ValueProvenance::USER_PROVIDED:
488
0
            nOutputBandCount = m_outputBandCountValue;
489
0
            break;
490
0
        case ValueProvenance::FROM_SOURCE:
491
0
            nOutputBandCount = m_poSrcDS->GetRasterCount();
492
0
            break;
493
0
        case ValueProvenance::FROM_VRTRASTERBAND:
494
0
            nOutputBandCount = nBands;
495
0
            break;
496
0
        case ValueProvenance::FROM_LAST_STEP:
497
0
            break;
498
0
    }
499
500
0
    const CPLXMLNode *psProcessingSteps =
501
0
        CPLGetXMLNode(psTree, "ProcessingSteps");
502
0
    if (!psProcessingSteps)
503
0
    {
504
0
        CPLError(CE_Failure, CPLE_AppDefined,
505
0
                 "ProcessingSteps element missing");
506
0
        return CE_Failure;
507
0
    }
508
509
0
    const auto eInDT = poSrcFirstBand->GetRasterDataType();
510
0
    for (int i = 1; i < m_poSrcDS->GetRasterCount(); ++i)
511
0
    {
512
0
        const auto eDT = m_poSrcDS->GetRasterBand(i + 1)->GetRasterDataType();
513
0
        if (eDT != eInDT)
514
0
        {
515
0
            CPLError(CE_Warning, CPLE_AppDefined,
516
0
                     "Not all bands of the input dataset have the same data "
517
0
                     "type. The data type of the first band will be used as "
518
0
                     "the reference one.");
519
0
            break;
520
0
        }
521
0
    }
522
523
0
    GDALDataType eCurrentDT = eInDT;
524
0
    int nCurrentBandCount = m_poSrcDS->GetRasterCount();
525
526
0
    std::vector<double> adfNoData;
527
0
    for (int i = 1; i <= nCurrentBandCount; ++i)
528
0
    {
529
0
        int bHasVal = FALSE;
530
0
        const double dfVal =
531
0
            m_poSrcDS->GetRasterBand(i)->GetNoDataValue(&bHasVal);
532
0
        adfNoData.emplace_back(
533
0
            bHasVal ? dfVal : std::numeric_limits<double>::quiet_NaN());
534
0
    }
535
536
0
    int nStepCount = 0;
537
0
    for (const CPLXMLNode *psStep = psProcessingSteps->psChild; psStep;
538
0
         psStep = psStep->psNext)
539
0
    {
540
0
        if (psStep->eType == CXT_Element &&
541
0
            strcmp(psStep->pszValue, "Step") == 0)
542
0
        {
543
0
            ++nStepCount;
544
0
        }
545
0
    }
546
547
0
    int iStep = 0;
548
0
    for (const CPLXMLNode *psStep = psProcessingSteps->psChild; psStep;
549
0
         psStep = psStep->psNext)
550
0
    {
551
0
        if (psStep->eType == CXT_Element &&
552
0
            strcmp(psStep->pszValue, "Step") == 0)
553
0
        {
554
0
            ++iStep;
555
0
            const bool bIsFinalStep = (iStep == nStepCount);
556
0
            std::vector<double> adfOutNoData;
557
0
            if (bIsFinalStep)
558
0
            {
559
                // Initialize adfOutNoData with nodata value of *output* bands
560
                // for final step
561
0
                if (m_outputBandCountProvenance ==
562
0
                    ValueProvenance::FROM_VRTRASTERBAND)
563
0
                {
564
0
                    for (int i = 1; i <= nBands; ++i)
565
0
                    {
566
0
                        int bHasVal = FALSE;
567
0
                        const double dfVal =
568
0
                            GetRasterBand(i)->GetNoDataValue(&bHasVal);
569
0
                        adfOutNoData.emplace_back(
570
0
                            bHasVal ? dfVal
571
0
                                    : std::numeric_limits<double>::quiet_NaN());
572
0
                    }
573
0
                }
574
0
                else if (m_outputBandCountProvenance ==
575
0
                         ValueProvenance::FROM_SOURCE)
576
0
                {
577
0
                    for (int i = 1; i <= m_poSrcDS->GetRasterCount(); ++i)
578
0
                    {
579
0
                        int bHasVal = FALSE;
580
0
                        const double dfVal =
581
0
                            m_poSrcDS->GetRasterBand(i)->GetNoDataValue(
582
0
                                &bHasVal);
583
0
                        adfOutNoData.emplace_back(
584
0
                            bHasVal ? dfVal
585
0
                                    : std::numeric_limits<double>::quiet_NaN());
586
0
                    }
587
0
                }
588
0
                else if (m_outputBandCountProvenance ==
589
0
                         ValueProvenance::USER_PROVIDED)
590
0
                {
591
0
                    adfOutNoData.resize(
592
0
                        m_outputBandCountValue,
593
0
                        std::numeric_limits<double>::quiet_NaN());
594
0
                }
595
0
            }
596
0
            if (!ParseStep(psStep, bIsFinalStep, eCurrentDT, nCurrentBandCount,
597
0
                           adfNoData, adfOutNoData))
598
0
                return CE_Failure;
599
0
            adfNoData = std::move(adfOutNoData);
600
0
        }
601
0
    }
602
603
0
    if (m_aoSteps.empty())
604
0
    {
605
0
        CPLError(CE_Failure, CPLE_AppDefined,
606
0
                 "At least one step should be defined");
607
0
        return CE_Failure;
608
0
    }
609
610
0
    int nLargestInDTSizeTimesBand = 1;
611
0
    int nLargestOutDTSizeTimesBand = 1;
612
0
    for (const auto &oStep : m_aoSteps)
613
0
    {
614
0
        const int nInDTSizeTimesBand =
615
0
            GDALGetDataTypeSizeBytes(oStep.eInDT) * oStep.nInBands;
616
0
        nLargestInDTSizeTimesBand =
617
0
            std::max(nLargestInDTSizeTimesBand, nInDTSizeTimesBand);
618
0
        const int nOutDTSizeTimesBand =
619
0
            GDALGetDataTypeSizeBytes(oStep.eOutDT) * oStep.nOutBands;
620
0
        nLargestOutDTSizeTimesBand =
621
0
            std::max(nLargestOutDTSizeTimesBand, nOutDTSizeTimesBand);
622
0
    }
623
0
    m_nWorkingBytesPerPixel =
624
0
        nLargestInDTSizeTimesBand + nLargestOutDTSizeTimesBand;
625
626
    // Use only up to 40% of RAM to acquire source bands and generate the output
627
    // buffer.
628
0
    m_nAllowedRAMUsage = CPLGetUsablePhysicalRAM() / 10 * 4;
629
    // Only for tests now
630
0
    const char *pszMAX_RAM = "VRT_PROCESSED_DATASET_ALLOWED_RAM_USAGE";
631
0
    if (const char *pszVal = CPLGetConfigOption(pszMAX_RAM, nullptr))
632
0
    {
633
0
        CPL_IGNORE_RET_VAL(
634
0
            CPLParseMemorySize(pszVal, &m_nAllowedRAMUsage, nullptr));
635
0
    }
636
637
0
    if (m_nAllowedRAMUsage > 0)
638
0
    {
639
0
        bool bBlockSizeModified = false;
640
0
        while ((m_nBlockXSize >= 2 || m_nBlockYSize >= 2) &&
641
0
               static_cast<GIntBig>(m_nBlockXSize) * m_nBlockYSize >
642
0
                   m_nAllowedRAMUsage / m_nWorkingBytesPerPixel)
643
0
        {
644
0
            if ((m_nBlockXSize == nRasterXSize ||
645
0
                 m_nBlockYSize >= m_nBlockXSize) &&
646
0
                m_nBlockYSize >= 2)
647
0
            {
648
0
                m_nBlockYSize /= 2;
649
0
            }
650
0
            else
651
0
            {
652
0
                m_nBlockXSize /= 2;
653
0
            }
654
0
            bBlockSizeModified = true;
655
0
        }
656
0
        if (bBlockSizeModified)
657
0
        {
658
0
            if (bUserBlockSize)
659
0
            {
660
0
                CPLError(
661
0
                    CE_Warning, CPLE_AppDefined,
662
0
                    "Reducing block size to %d x %d to avoid consuming too "
663
0
                    "much RAM",
664
0
                    m_nBlockXSize, m_nBlockYSize);
665
0
            }
666
0
            else
667
0
            {
668
0
                CPLDebug(
669
0
                    "VRT",
670
0
                    "Reducing block size to %d x %d to avoid consuming too "
671
0
                    "much RAM",
672
0
                    m_nBlockXSize, m_nBlockYSize);
673
0
            }
674
0
        }
675
0
    }
676
677
0
    if (m_outputBandCountProvenance == ValueProvenance::FROM_LAST_STEP)
678
0
    {
679
0
        nOutputBandCount = nCurrentBandCount;
680
0
    }
681
0
    else if (nOutputBandCount != nCurrentBandCount)
682
0
    {
683
        // Should not happen frequently as pixel init functions are expected
684
        // to validate that they can accept the number of output bands provided
685
        // to them
686
0
        CPLError(
687
0
            CE_Failure, CPLE_AppDefined,
688
0
            "Number of output bands of last step (%d) is not consistent with "
689
0
            "number of VRTProcessedRasterBand's (%d)",
690
0
            nCurrentBandCount, nBands);
691
0
        return CE_Failure;
692
0
    }
693
694
0
    if (m_outputBandDataTypeProvenance == ValueProvenance::FROM_LAST_STEP)
695
0
    {
696
0
        m_outputBandDataTypeValue = eCurrentDT;
697
0
    }
698
699
0
    if (nBands != 0 &&
700
0
        (nBands != nOutputBandCount ||
701
0
         (m_outputBandDataTypeProvenance == ValueProvenance::FROM_LAST_STEP &&
702
0
          m_outputBandDataTypeValue != papoBands[0]->GetRasterDataType())))
703
0
    {
704
0
        for (int i = 0; i < nBands; ++i)
705
0
            delete papoBands[i];
706
0
        CPLFree(papoBands);
707
0
        papoBands = nullptr;
708
0
        nBands = 0;
709
0
    }
710
711
0
    const auto GetOutputBandType = [this, eCurrentDT](GDALDataType eSourceDT)
712
0
    {
713
0
        if (m_outputBandDataTypeProvenance == ValueProvenance::FROM_LAST_STEP)
714
0
            return eCurrentDT;
715
0
        else if (m_outputBandDataTypeProvenance ==
716
0
                 ValueProvenance::USER_PROVIDED)
717
0
            return m_outputBandDataTypeValue;
718
0
        else
719
0
            return eSourceDT;
720
0
    };
721
722
0
    if (m_outputBandCountProvenance == ValueProvenance::FROM_SOURCE)
723
0
    {
724
0
        for (int i = 0; i < m_poSrcDS->GetRasterCount(); ++i)
725
0
        {
726
0
            const auto poSrcBand = m_poSrcDS->GetRasterBand(i + 1);
727
0
            const GDALDataType eOutputBandType =
728
0
                GetOutputBandType(poSrcBand->GetRasterDataType());
729
0
            auto poBand =
730
0
                new VRTProcessedRasterBand(this, i + 1, eOutputBandType);
731
0
            poBand->CopyCommonInfoFrom(poSrcBand);
732
0
            SetBand(i + 1, poBand);
733
0
        }
734
0
    }
735
0
    else if (m_outputBandCountProvenance != ValueProvenance::FROM_VRTRASTERBAND)
736
0
    {
737
0
        const GDALDataType eOutputBandType = GetOutputBandType(eInDT);
738
0
        for (int i = 0; i < nOutputBandCount; ++i)
739
0
        {
740
0
            auto poBand =
741
0
                new VRTProcessedRasterBand(this, i + 1, eOutputBandType);
742
0
            SetBand(i + 1, poBand);
743
0
        }
744
0
    }
745
746
0
    if (nBands > 1)
747
0
        SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
748
749
0
    m_oXMLTree.reset(CPLCloneXMLTree(psTree));
750
751
0
    return CE_None;
752
0
}
753
754
/************************************************************************/
755
/*                            ParseStep()                               */
756
/************************************************************************/
757
758
/** Parse the current Step node and create a corresponding entry in m_aoSteps.
759
 *
760
 * @param psStep Step node
761
 * @param bIsFinalStep Whether this is the final step.
762
 * @param[in,out] eCurrentDT Input data type for this step.
763
 *                           Updated to output data type at end of method.
764
 * @param[in,out] nCurrentBandCount Input band count for this step.
765
 *                                  Updated to output band cout at end of
766
 *                                  method.
767
 * @param adfInNoData Input nodata values
768
 * @param[in,out] adfOutNoData Output nodata values, to be filled by this
769
 *                             method. When bIsFinalStep, this is also an
770
 *                             input parameter.
771
 * @return true on success.
772
 */
773
bool VRTProcessedDataset::ParseStep(const CPLXMLNode *psStep, bool bIsFinalStep,
774
                                    GDALDataType &eCurrentDT,
775
                                    int &nCurrentBandCount,
776
                                    std::vector<double> &adfInNoData,
777
                                    std::vector<double> &adfOutNoData)
778
0
{
779
0
    const char *pszStepName = CPLGetXMLValue(
780
0
        psStep, "name", CPLSPrintf("nr %d", 1 + int(m_aoSteps.size())));
781
0
    const char *pszAlgorithm = CPLGetXMLValue(psStep, "Algorithm", nullptr);
782
0
    if (!pszAlgorithm)
783
0
    {
784
0
        CPLError(CE_Failure, CPLE_AppDefined,
785
0
                 "Step '%s' lacks a Algorithm element", pszStepName);
786
0
        return false;
787
0
    }
788
789
0
    const auto &oMapFunctions = GetGlobalMapProcessedDatasetFunc();
790
0
    const auto oIterFunc = oMapFunctions.find(pszAlgorithm);
791
0
    if (oIterFunc == oMapFunctions.end())
792
0
    {
793
0
        CPLError(CE_Failure, CPLE_AppDefined,
794
0
                 "Step '%s' uses unregistered algorithm '%s'", pszStepName,
795
0
                 pszAlgorithm);
796
0
        return false;
797
0
    }
798
799
0
    const auto &oFunc = oIterFunc->second;
800
801
0
    if (!oFunc.aeSupportedInputDT.empty())
802
0
    {
803
0
        if (std::find(oFunc.aeSupportedInputDT.begin(),
804
0
                      oFunc.aeSupportedInputDT.end(),
805
0
                      eCurrentDT) == oFunc.aeSupportedInputDT.end())
806
0
        {
807
0
            CPLError(CE_Failure, CPLE_AppDefined,
808
0
                     "Step '%s' (using algorithm '%s') does not "
809
0
                     "support input data type = '%s'",
810
0
                     pszStepName, pszAlgorithm,
811
0
                     GDALGetDataTypeName(eCurrentDT));
812
0
            return false;
813
0
        }
814
0
    }
815
816
0
    if (!oFunc.anSupportedInputBandCount.empty())
817
0
    {
818
0
        if (std::find(oFunc.anSupportedInputBandCount.begin(),
819
0
                      oFunc.anSupportedInputBandCount.end(),
820
0
                      nCurrentBandCount) ==
821
0
            oFunc.anSupportedInputBandCount.end())
822
0
        {
823
0
            CPLError(CE_Failure, CPLE_AppDefined,
824
0
                     "Step '%s' (using algorithm '%s') does not "
825
0
                     "support input band count = %d",
826
0
                     pszStepName, pszAlgorithm, nCurrentBandCount);
827
0
            return false;
828
0
        }
829
0
    }
830
831
0
    Step oStep;
832
0
    oStep.osAlgorithm = pszAlgorithm;
833
0
    oStep.eInDT = oFunc.eRequestedInputDT != GDT_Unknown
834
0
                      ? oFunc.eRequestedInputDT
835
0
                      : eCurrentDT;
836
0
    oStep.nInBands = nCurrentBandCount;
837
838
    // Unless modified by pfnInit...
839
0
    oStep.eOutDT = oStep.eInDT;
840
841
0
    oStep.adfInNoData = adfInNoData;
842
0
    oStep.adfOutNoData = bIsFinalStep ? adfOutNoData : adfInNoData;
843
844
    // Deal with constant arguments
845
0
    for (const auto &nameValuePair : oFunc.oMapConstantArguments)
846
0
    {
847
0
        oStep.aosArguments.AddNameValue(nameValuePair.first.c_str(),
848
0
                                        nameValuePair.second.c_str());
849
0
    }
850
851
    // Deal with built-in arguments
852
0
    if (oFunc.oSetBuiltinArguments.find("nodata") !=
853
0
        oFunc.oSetBuiltinArguments.end())
854
0
    {
855
0
        int bHasVal = false;
856
0
        const auto poSrcFirstBand = m_poSrcDS->GetRasterBand(1);
857
0
        const double dfVal = poSrcFirstBand->GetNoDataValue(&bHasVal);
858
0
        if (bHasVal)
859
0
        {
860
0
            oStep.aosArguments.AddNameValue("nodata",
861
0
                                            CPLSPrintf("%.17g", dfVal));
862
0
        }
863
0
    }
864
865
0
    if (oFunc.oSetBuiltinArguments.find("offset_{band}") !=
866
0
        oFunc.oSetBuiltinArguments.end())
867
0
    {
868
0
        for (int i = 1; i <= m_poSrcDS->GetRasterCount(); ++i)
869
0
        {
870
0
            int bHasVal = false;
871
0
            const double dfVal = GetRasterBand(i)->GetOffset(&bHasVal);
872
0
            oStep.aosArguments.AddNameValue(
873
0
                CPLSPrintf("offset_%d", i),
874
0
                CPLSPrintf("%.17g", bHasVal ? dfVal : 0.0));
875
0
        }
876
0
    }
877
878
0
    if (oFunc.oSetBuiltinArguments.find("scale_{band}") !=
879
0
        oFunc.oSetBuiltinArguments.end())
880
0
    {
881
0
        for (int i = 1; i <= m_poSrcDS->GetRasterCount(); ++i)
882
0
        {
883
0
            int bHasVal = false;
884
0
            const double dfVal = GetRasterBand(i)->GetScale(&bHasVal);
885
0
            oStep.aosArguments.AddNameValue(
886
0
                CPLSPrintf("scale_%d", i),
887
0
                CPLSPrintf("%.17g", bHasVal ? dfVal : 1.0));
888
0
        }
889
0
    }
890
891
    // Parse arguments specified in VRT
892
0
    std::set<std::string> oFoundArguments;
893
894
0
    for (const CPLXMLNode *psStepChild = psStep->psChild; psStepChild;
895
0
         psStepChild = psStepChild->psNext)
896
0
    {
897
0
        if (psStepChild->eType == CXT_Element &&
898
0
            strcmp(psStepChild->pszValue, "Argument") == 0)
899
0
        {
900
0
            const char *pszParamName =
901
0
                CPLGetXMLValue(psStepChild, "name", nullptr);
902
0
            if (!pszParamName)
903
0
            {
904
0
                CPLError(CE_Failure, CPLE_AppDefined,
905
0
                         "Step '%s' has a Argument without a name attribute",
906
0
                         pszStepName);
907
0
                return false;
908
0
            }
909
0
            const char *pszValue = CPLGetXMLValue(psStepChild, nullptr, "");
910
0
            auto oOtherArgIter =
911
0
                oFunc.oOtherArguments.find(CPLString(pszParamName).tolower());
912
0
            if (!oFunc.oOtherArguments.empty() &&
913
0
                oOtherArgIter == oFunc.oOtherArguments.end())
914
0
            {
915
                // If we got a parameter name like 'coefficients_1',
916
                // try to fetch the generic 'coefficients_{band}'
917
0
                std::string osParamName(pszParamName);
918
0
                const auto nPos = osParamName.rfind('_');
919
0
                if (nPos != std::string::npos)
920
0
                {
921
0
                    osParamName.resize(nPos + 1);
922
0
                    osParamName += "{band}";
923
0
                    oOtherArgIter = oFunc.oOtherArguments.find(
924
0
                        CPLString(osParamName).tolower());
925
0
                }
926
0
            }
927
0
            if (oOtherArgIter != oFunc.oOtherArguments.end())
928
0
            {
929
0
                oFoundArguments.insert(oOtherArgIter->first);
930
931
0
                const std::string &osType = oOtherArgIter->second.osType;
932
0
                if (osType == "boolean")
933
0
                {
934
0
                    if (!EQUAL(pszValue, "true") && !EQUAL(pszValue, "false"))
935
0
                    {
936
0
                        CPLError(CE_Failure, CPLE_NotSupported,
937
0
                                 "Step '%s' has a Argument '%s' whose "
938
0
                                 "value '%s' is not a boolean",
939
0
                                 pszStepName, pszParamName, pszValue);
940
0
                        return false;
941
0
                    }
942
0
                }
943
0
                else if (osType == "integer")
944
0
                {
945
0
                    if (CPLGetValueType(pszValue) != CPL_VALUE_INTEGER)
946
0
                    {
947
0
                        CPLError(CE_Failure, CPLE_NotSupported,
948
0
                                 "Step '%s' has a Argument '%s' whose "
949
0
                                 "value '%s' is not a integer",
950
0
                                 pszStepName, pszParamName, pszValue);
951
0
                        return false;
952
0
                    }
953
0
                }
954
0
                else if (osType == "double")
955
0
                {
956
0
                    const auto eType = CPLGetValueType(pszValue);
957
0
                    if (eType != CPL_VALUE_INTEGER && eType != CPL_VALUE_REAL)
958
0
                    {
959
0
                        CPLError(CE_Failure, CPLE_NotSupported,
960
0
                                 "Step '%s' has a Argument '%s' whose "
961
0
                                 "value '%s' is not a double",
962
0
                                 pszStepName, pszParamName, pszValue);
963
0
                        return false;
964
0
                    }
965
0
                }
966
0
                else if (osType == "double_list")
967
0
                {
968
0
                    const CPLStringList aosTokens(
969
0
                        CSLTokenizeString2(pszValue, ",", 0));
970
0
                    for (int i = 0; i < aosTokens.size(); ++i)
971
0
                    {
972
0
                        const auto eType = CPLGetValueType(aosTokens[i]);
973
0
                        if (eType != CPL_VALUE_INTEGER &&
974
0
                            eType != CPL_VALUE_REAL)
975
0
                        {
976
0
                            CPLError(CE_Failure, CPLE_NotSupported,
977
0
                                     "Step '%s' has a Argument '%s' "
978
0
                                     "whose value '%s' is not a "
979
0
                                     "comma-separated list of doubles",
980
0
                                     pszStepName, pszParamName, pszValue);
981
0
                            return false;
982
0
                        }
983
0
                    }
984
0
                }
985
0
                else if (osType != "string")
986
0
                {
987
0
                    CPLDebug("VRT", "Unhandled argument type '%s'",
988
0
                             osType.c_str());
989
0
                    CPLAssert(0);
990
0
                }
991
0
            }
992
0
            else if (oFunc.bMetadataSpecified &&
993
0
                     oFunc.oSetBuiltinArguments.find(
994
0
                         CPLString(pszParamName).tolower()) ==
995
0
                         oFunc.oSetBuiltinArguments.end() &&
996
0
                     oFunc.oMapConstantArguments.find(
997
0
                         CPLString(pszParamName).tolower()) ==
998
0
                         oFunc.oMapConstantArguments.end())
999
0
            {
1000
0
                CPLError(CE_Warning, CPLE_NotSupported,
1001
0
                         "Step '%s' has a Argument '%s' which is not "
1002
0
                         "supported",
1003
0
                         pszStepName, pszParamName);
1004
0
            }
1005
1006
0
            oStep.aosArguments.AddNameValue(pszParamName, pszValue);
1007
0
        }
1008
0
    }
1009
1010
    // Check that required arguments have been specified
1011
0
    for (const auto &oIter : oFunc.oOtherArguments)
1012
0
    {
1013
0
        if (oIter.second.bRequired &&
1014
0
            oFoundArguments.find(oIter.first) == oFoundArguments.end())
1015
0
        {
1016
0
            CPLError(CE_Failure, CPLE_AppDefined,
1017
0
                     "Step '%s' lacks required Argument '%s'", pszStepName,
1018
0
                     oIter.first.c_str());
1019
0
            return false;
1020
0
        }
1021
0
    }
1022
1023
0
    if (oFunc.pfnInit)
1024
0
    {
1025
0
        double *padfOutNoData = nullptr;
1026
0
        if (bIsFinalStep && !adfOutNoData.empty())
1027
0
        {
1028
0
            oStep.nOutBands = static_cast<int>(adfOutNoData.size());
1029
0
            padfOutNoData = static_cast<double *>(
1030
0
                CPLMalloc(adfOutNoData.size() * sizeof(double)));
1031
0
            memcpy(padfOutNoData, adfOutNoData.data(),
1032
0
                   adfOutNoData.size() * sizeof(double));
1033
0
        }
1034
0
        else
1035
0
        {
1036
0
            oStep.nOutBands = 0;
1037
0
        }
1038
1039
0
        if (oFunc.pfnInit(pszAlgorithm, oFunc.pUserData,
1040
0
                          oStep.aosArguments.List(), oStep.nInBands,
1041
0
                          oStep.eInDT, adfInNoData.data(), &(oStep.nOutBands),
1042
0
                          &(oStep.eOutDT), &padfOutNoData, m_osVRTPath.c_str(),
1043
0
                          &(oStep.pWorkingData)) != CE_None)
1044
0
        {
1045
0
            CPLError(CE_Failure, CPLE_AppDefined,
1046
0
                     "Step '%s' (using algorithm '%s') init() function "
1047
0
                     "failed",
1048
0
                     pszStepName, pszAlgorithm);
1049
0
            CPLFree(padfOutNoData);
1050
0
            return false;
1051
0
        }
1052
1053
        // Input nodata values may have been modified by pfnInit()
1054
0
        oStep.adfInNoData = adfInNoData;
1055
1056
0
        if (padfOutNoData)
1057
0
        {
1058
0
            adfOutNoData = std::vector<double>(padfOutNoData,
1059
0
                                               padfOutNoData + oStep.nOutBands);
1060
0
        }
1061
0
        else
1062
0
        {
1063
0
            adfOutNoData = std::vector<double>(
1064
0
                oStep.nOutBands, std::numeric_limits<double>::quiet_NaN());
1065
0
        }
1066
0
        CPLFree(padfOutNoData);
1067
1068
0
        oStep.adfOutNoData = adfOutNoData;
1069
0
    }
1070
0
    else
1071
0
    {
1072
0
        oStep.nOutBands = oStep.nInBands;
1073
0
        adfOutNoData = oStep.adfOutNoData;
1074
0
    }
1075
1076
0
    eCurrentDT = oStep.eOutDT;
1077
0
    nCurrentBandCount = oStep.nOutBands;
1078
1079
0
    m_aoSteps.emplace_back(std::move(oStep));
1080
1081
0
    return true;
1082
0
}
1083
1084
/************************************************************************/
1085
/*                           SerializeToXML()                           */
1086
/************************************************************************/
1087
1088
CPLXMLNode *VRTProcessedDataset::SerializeToXML(const char *pszVRTPathIn)
1089
1090
0
{
1091
0
    CPLXMLNode *psTree = CPLCloneXMLTree(m_oXMLTree.get());
1092
0
    if (psTree == nullptr)
1093
0
        return psTree;
1094
1095
    /* -------------------------------------------------------------------- */
1096
    /*      Remove VRTRasterBand nodes from the original tree and find the  */
1097
    /*      last child.                                                     */
1098
    /* -------------------------------------------------------------------- */
1099
0
    CPLXMLNode *psLastChild = psTree->psChild;
1100
0
    CPLXMLNode *psPrevChild = nullptr;
1101
0
    while (psLastChild)
1102
0
    {
1103
0
        CPLXMLNode *psNextChild = psLastChild->psNext;
1104
0
        if (psLastChild->eType == CXT_Element &&
1105
0
            strcmp(psLastChild->pszValue, "VRTRasterBand") == 0)
1106
0
        {
1107
0
            if (psPrevChild)
1108
0
                psPrevChild->psNext = psNextChild;
1109
0
            else
1110
0
                psTree->psChild = psNextChild;
1111
0
            psLastChild->psNext = nullptr;
1112
0
            CPLDestroyXMLNode(psLastChild);
1113
0
            psLastChild = psPrevChild ? psPrevChild : psTree->psChild;
1114
0
        }
1115
0
        else if (!psNextChild)
1116
0
        {
1117
0
            break;
1118
0
        }
1119
0
        else
1120
0
        {
1121
0
            psPrevChild = psLastChild;
1122
0
            psLastChild = psNextChild;
1123
0
        }
1124
0
    }
1125
0
    CPLAssert(psLastChild);  // we have at least Input
1126
1127
    /* -------------------------------------------------------------------- */
1128
    /*      Serialize bands.                                                */
1129
    /* -------------------------------------------------------------------- */
1130
0
    bool bHasWarnedAboutRAMUsage = false;
1131
0
    size_t nAccRAMUsage = 0;
1132
0
    for (int iBand = 0; iBand < nBands; iBand++)
1133
0
    {
1134
0
        CPLXMLNode *psBandTree =
1135
0
            static_cast<VRTRasterBand *>(papoBands[iBand])
1136
0
                ->SerializeToXML(pszVRTPathIn, bHasWarnedAboutRAMUsage,
1137
0
                                 nAccRAMUsage);
1138
1139
0
        if (psBandTree != nullptr)
1140
0
        {
1141
0
            psLastChild->psNext = psBandTree;
1142
0
            psLastChild = psBandTree;
1143
0
        }
1144
0
    }
1145
1146
0
    return psTree;
1147
0
}
1148
1149
/************************************************************************/
1150
/*                           SerializeToXML()                           */
1151
/************************************************************************/
1152
1153
CPLXMLNode *
1154
VRTProcessedRasterBand::SerializeToXML(const char *pszVRTPathIn,
1155
                                       bool &bHasWarnedAboutRAMUsage,
1156
                                       size_t &nAccRAMUsage)
1157
1158
0
{
1159
0
    CPLXMLNode *psTree = VRTRasterBand::SerializeToXML(
1160
0
        pszVRTPathIn, bHasWarnedAboutRAMUsage, nAccRAMUsage);
1161
1162
    /* -------------------------------------------------------------------- */
1163
    /*      Set subclass.                                                   */
1164
    /* -------------------------------------------------------------------- */
1165
0
    CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
1166
0
                     CXT_Text, "VRTProcessedRasterBand");
1167
1168
0
    return psTree;
1169
0
}
1170
1171
/************************************************************************/
1172
/*                            GetBlockSize()                            */
1173
/************************************************************************/
1174
1175
/** Return block size */
1176
void VRTProcessedDataset::GetBlockSize(int *pnBlockXSize,
1177
                                       int *pnBlockYSize) const
1178
1179
0
{
1180
0
    *pnBlockXSize = m_nBlockXSize;
1181
0
    *pnBlockYSize = m_nBlockYSize;
1182
0
}
1183
1184
/************************************************************************/
1185
/*                            ProcessRegion()                           */
1186
/************************************************************************/
1187
1188
/** Compute pixel values for the specified region.
1189
 *
1190
 * The output is stored in m_abyInput in a pixel-interleaved way.
1191
 */
1192
bool VRTProcessedDataset::ProcessRegion(int nXOff, int nYOff, int nBufXSize,
1193
                                        int nBufYSize,
1194
                                        GDALProgressFunc pfnProgress,
1195
                                        void *pProgressData)
1196
0
{
1197
1198
0
    CPLAssert(!m_aoSteps.empty());
1199
1200
0
    const size_t nPixelCount = static_cast<size_t>(nBufXSize) * nBufYSize;
1201
1202
0
    const int nFirstBandCount = m_aoSteps.front().nInBands;
1203
0
    CPLAssert(nFirstBandCount == m_poSrcDS->GetRasterCount());
1204
0
    const GDALDataType eFirstDT = m_aoSteps.front().eInDT;
1205
0
    const int nFirstDTSize = GDALGetDataTypeSizeBytes(eFirstDT);
1206
0
    auto &abyInput = m_abyInput;
1207
0
    auto &abyOutput = m_abyOutput;
1208
1209
0
    const char *pszInterleave =
1210
0
        m_poSrcDS->GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE");
1211
0
    if (nFirstBandCount > 1 && (!pszInterleave || EQUAL(pszInterleave, "BAND")))
1212
0
    {
1213
        // If there are several bands and the source dataset organization
1214
        // is apparently band interleaved, then first acquire data in
1215
        // a BSQ organization in the abyInput array use in the native
1216
        // data type.
1217
        // And then transpose it and convert it to the expected data type
1218
        // of the first step.
1219
0
        const auto eSrcDT = m_poSrcDS->GetRasterBand(1)->GetRasterDataType();
1220
0
        try
1221
0
        {
1222
0
            abyInput.resize(nPixelCount * nFirstBandCount *
1223
0
                            GDALGetDataTypeSizeBytes(eSrcDT));
1224
0
        }
1225
0
        catch (const std::bad_alloc &)
1226
0
        {
1227
0
            CPLError(CE_Failure, CPLE_OutOfMemory,
1228
0
                     "Out of memory allocating working buffer");
1229
0
            return false;
1230
0
        }
1231
1232
0
        try
1233
0
        {
1234
0
            abyOutput.resize(nPixelCount * nFirstBandCount * nFirstDTSize);
1235
0
        }
1236
0
        catch (const std::bad_alloc &)
1237
0
        {
1238
0
            CPLError(CE_Failure, CPLE_OutOfMemory,
1239
0
                     "Out of memory allocating working buffer");
1240
0
            return false;
1241
0
        }
1242
1243
0
        GDALRasterIOExtraArg sArg;
1244
0
        INIT_RASTERIO_EXTRA_ARG(sArg);
1245
0
        sArg.pfnProgress = GDALScaledProgress;
1246
0
        sArg.pProgressData =
1247
0
            GDALCreateScaledProgress(0, 0.5, pfnProgress, pProgressData);
1248
0
        if (sArg.pProgressData == nullptr)
1249
0
            sArg.pfnProgress = nullptr;
1250
1251
0
        CPLDebugOnly("VRT", "ProcessRegion(): start RasterIO()");
1252
0
        const bool bOK =
1253
0
            m_poSrcDS->RasterIO(GF_Read, nXOff, nYOff, nBufXSize, nBufYSize,
1254
0
                                abyInput.data(), nBufXSize, nBufYSize, eSrcDT,
1255
0
                                nFirstBandCount, nullptr, 0, 0, 0,
1256
0
                                &sArg) == CE_None;
1257
0
        CPLDebugOnly("VRT", "ProcessRegion(): end RasterIO()");
1258
0
        GDALDestroyScaledProgress(sArg.pProgressData);
1259
0
        if (!bOK)
1260
0
            return false;
1261
1262
0
        CPLDebugOnly("VRT", "ProcessRegion(): start GDALTranspose2D()");
1263
0
        GDALTranspose2D(abyInput.data(), eSrcDT, abyOutput.data(), eFirstDT,
1264
0
                        static_cast<size_t>(nBufXSize) * nBufYSize,
1265
0
                        nFirstBandCount);
1266
0
        CPLDebugOnly("VRT", "ProcessRegion(): end GDALTranspose2D()");
1267
1268
        // Swap arrays
1269
0
        std::swap(abyInput, abyOutput);
1270
0
    }
1271
0
    else
1272
0
    {
1273
0
        try
1274
0
        {
1275
0
            abyInput.resize(nPixelCount * nFirstBandCount * nFirstDTSize);
1276
0
        }
1277
0
        catch (const std::bad_alloc &)
1278
0
        {
1279
0
            CPLError(CE_Failure, CPLE_OutOfMemory,
1280
0
                     "Out of memory allocating working buffer");
1281
0
            return false;
1282
0
        }
1283
1284
0
        GDALRasterIOExtraArg sArg;
1285
0
        INIT_RASTERIO_EXTRA_ARG(sArg);
1286
0
        sArg.pfnProgress = GDALScaledProgress;
1287
0
        sArg.pProgressData =
1288
0
            GDALCreateScaledProgress(0, 0.5, pfnProgress, pProgressData);
1289
0
        if (sArg.pProgressData == nullptr)
1290
0
            sArg.pfnProgress = nullptr;
1291
1292
0
        const bool bOK =
1293
0
            m_poSrcDS->RasterIO(
1294
0
                GF_Read, nXOff, nYOff, nBufXSize, nBufYSize, abyInput.data(),
1295
0
                nBufXSize, nBufYSize, eFirstDT, nFirstBandCount, nullptr,
1296
0
                static_cast<GSpacing>(nFirstDTSize) * nFirstBandCount,
1297
0
                static_cast<GSpacing>(nFirstDTSize) * nFirstBandCount *
1298
0
                    nBufXSize,
1299
0
                nFirstDTSize, &sArg) == CE_None;
1300
1301
0
        GDALDestroyScaledProgress(sArg.pProgressData);
1302
0
        if (!bOK)
1303
0
            return false;
1304
0
    }
1305
1306
0
    const double dfSrcXOff = nXOff;
1307
0
    const double dfSrcYOff = nYOff;
1308
0
    const double dfSrcXSize = nBufXSize;
1309
0
    const double dfSrcYSize = nBufYSize;
1310
1311
0
    double adfSrcGT[6];
1312
0
    if (m_poSrcDS->GetGeoTransform(adfSrcGT) != CE_None)
1313
0
    {
1314
0
        adfSrcGT[0] = 0;
1315
0
        adfSrcGT[1] = 1;
1316
0
        adfSrcGT[2] = 0;
1317
0
        adfSrcGT[3] = 0;
1318
0
        adfSrcGT[4] = 0;
1319
0
        adfSrcGT[5] = 1;
1320
0
    }
1321
1322
0
    GDALDataType eLastDT = eFirstDT;
1323
0
    const auto &oMapFunctions = GetGlobalMapProcessedDatasetFunc();
1324
1325
0
    int iStep = 0;
1326
0
    for (const auto &oStep : m_aoSteps)
1327
0
    {
1328
0
        const auto oIterFunc = oMapFunctions.find(oStep.osAlgorithm);
1329
0
        CPLAssert(oIterFunc != oMapFunctions.end());
1330
1331
        // Data type adaptation
1332
0
        if (eLastDT != oStep.eInDT)
1333
0
        {
1334
0
            try
1335
0
            {
1336
0
                abyOutput.resize(nPixelCount * oStep.nInBands *
1337
0
                                 GDALGetDataTypeSizeBytes(oStep.eInDT));
1338
0
            }
1339
0
            catch (const std::bad_alloc &)
1340
0
            {
1341
0
                CPLError(CE_Failure, CPLE_OutOfMemory,
1342
0
                         "Out of memory allocating working buffer");
1343
0
                return false;
1344
0
            }
1345
1346
0
            GDALCopyWords64(abyInput.data(), eLastDT,
1347
0
                            GDALGetDataTypeSizeBytes(eLastDT), abyOutput.data(),
1348
0
                            oStep.eInDT, GDALGetDataTypeSizeBytes(oStep.eInDT),
1349
0
                            nPixelCount * oStep.nInBands);
1350
1351
0
            std::swap(abyInput, abyOutput);
1352
0
        }
1353
1354
0
        try
1355
0
        {
1356
0
            abyOutput.resize(nPixelCount * oStep.nOutBands *
1357
0
                             GDALGetDataTypeSizeBytes(oStep.eOutDT));
1358
0
        }
1359
0
        catch (const std::bad_alloc &)
1360
0
        {
1361
0
            CPLError(CE_Failure, CPLE_OutOfMemory,
1362
0
                     "Out of memory allocating working buffer");
1363
0
            return false;
1364
0
        }
1365
1366
0
        const auto &oFunc = oIterFunc->second;
1367
0
        if (oFunc.pfnProcess(
1368
0
                oStep.osAlgorithm.c_str(), oFunc.pUserData, oStep.pWorkingData,
1369
0
                oStep.aosArguments.List(), nBufXSize, nBufYSize,
1370
0
                abyInput.data(), abyInput.size(), oStep.eInDT, oStep.nInBands,
1371
0
                oStep.adfInNoData.data(), abyOutput.data(), abyOutput.size(),
1372
0
                oStep.eOutDT, oStep.nOutBands, oStep.adfOutNoData.data(),
1373
0
                dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize, adfSrcGT,
1374
0
                m_osVRTPath.c_str(),
1375
0
                /*papszExtra=*/nullptr) != CE_None)
1376
0
        {
1377
0
            return false;
1378
0
        }
1379
1380
0
        std::swap(abyInput, abyOutput);
1381
0
        eLastDT = oStep.eOutDT;
1382
1383
0
        ++iStep;
1384
0
        if (pfnProgress &&
1385
0
            !pfnProgress(0.5 + 0.5 * iStep / static_cast<int>(m_aoSteps.size()),
1386
0
                         "", pProgressData))
1387
0
            return false;
1388
0
    }
1389
1390
0
    return true;
1391
0
}
1392
1393
/************************************************************************/
1394
/*                        VRTProcessedRasterBand()                      */
1395
/************************************************************************/
1396
1397
/** Constructor */
1398
VRTProcessedRasterBand::VRTProcessedRasterBand(VRTProcessedDataset *poDSIn,
1399
                                               int nBandIn,
1400
                                               GDALDataType eDataTypeIn)
1401
0
{
1402
0
    Initialize(poDSIn->GetRasterXSize(), poDSIn->GetRasterYSize());
1403
1404
0
    poDS = poDSIn;
1405
0
    nBand = nBandIn;
1406
0
    eAccess = GA_Update;
1407
0
    eDataType = eDataTypeIn;
1408
1409
0
    poDSIn->GetBlockSize(&nBlockXSize, &nBlockYSize);
1410
0
}
1411
1412
/************************************************************************/
1413
/*                           GetOverviewCount()                         */
1414
/************************************************************************/
1415
1416
int VRTProcessedRasterBand::GetOverviewCount()
1417
0
{
1418
0
    auto poVRTDS = cpl::down_cast<VRTProcessedDataset *>(poDS);
1419
0
    return static_cast<int>(poVRTDS->m_apoOverviewDatasets.size());
1420
0
}
1421
1422
/************************************************************************/
1423
/*                              GetOverview()                           */
1424
/************************************************************************/
1425
1426
GDALRasterBand *VRTProcessedRasterBand::GetOverview(int iOvr)
1427
0
{
1428
0
    auto poVRTDS = cpl::down_cast<VRTProcessedDataset *>(poDS);
1429
0
    if (iOvr < 0 ||
1430
0
        iOvr >= static_cast<int>(poVRTDS->m_apoOverviewDatasets.size()))
1431
0
        return nullptr;
1432
0
    return poVRTDS->m_apoOverviewDatasets[iOvr]->GetRasterBand(nBand);
1433
0
}
1434
1435
/************************************************************************/
1436
/*                             IReadBlock()                             */
1437
/************************************************************************/
1438
1439
CPLErr VRTProcessedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
1440
                                          void *pImage)
1441
1442
0
{
1443
0
    auto poVRTDS = cpl::down_cast<VRTProcessedDataset *>(poDS);
1444
1445
0
    int nBufXSize = 0;
1446
0
    int nBufYSize = 0;
1447
0
    GetActualBlockSize(nBlockXOff, nBlockYOff, &nBufXSize, &nBufYSize);
1448
1449
0
    const int nXPixelOff = nBlockXOff * nBlockXSize;
1450
0
    const int nYPixelOff = nBlockYOff * nBlockYSize;
1451
0
    if (!poVRTDS->ProcessRegion(nXPixelOff, nYPixelOff, nBufXSize, nBufYSize,
1452
0
                                nullptr, nullptr))
1453
0
    {
1454
0
        return CE_Failure;
1455
0
    }
1456
1457
0
    const int nOutBands = poVRTDS->m_aoSteps.back().nOutBands;
1458
0
    CPLAssert(nOutBands == poVRTDS->GetRasterCount());
1459
0
    const auto eLastDT = poVRTDS->m_aoSteps.back().eOutDT;
1460
0
    const int nLastDTSize = GDALGetDataTypeSizeBytes(eLastDT);
1461
0
    const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
1462
1463
    // Dispatch final output buffer to cached blocks of output bands
1464
0
    for (int iDstBand = 0; iDstBand < nOutBands; ++iDstBand)
1465
0
    {
1466
0
        GDALRasterBlock *poBlock = nullptr;
1467
0
        GByte *pDst;
1468
0
        if (iDstBand + 1 == nBand)
1469
0
        {
1470
0
            pDst = static_cast<GByte *>(pImage);
1471
0
        }
1472
0
        else
1473
0
        {
1474
0
            auto poOtherBand = poVRTDS->papoBands[iDstBand];
1475
0
            poBlock = poOtherBand->TryGetLockedBlockRef(nBlockXOff, nBlockYOff);
1476
0
            if (poBlock)
1477
0
            {
1478
0
                poBlock->DropLock();
1479
0
                continue;
1480
0
            }
1481
0
            poBlock = poOtherBand->GetLockedBlockRef(
1482
0
                nBlockXOff, nBlockYOff, /* bJustInitialized = */ true);
1483
0
            if (!poBlock)
1484
0
                continue;
1485
0
            pDst = static_cast<GByte *>(poBlock->GetDataRef());
1486
0
        }
1487
0
        for (int iY = 0; iY < nBufYSize; ++iY)
1488
0
        {
1489
0
            GDALCopyWords64(poVRTDS->m_abyInput.data() +
1490
0
                                (iDstBand + static_cast<size_t>(iY) *
1491
0
                                                nBufXSize * nOutBands) *
1492
0
                                    nLastDTSize,
1493
0
                            eLastDT, nLastDTSize * nOutBands,
1494
0
                            pDst +
1495
0
                                static_cast<size_t>(iY) * nBlockXSize * nDTSize,
1496
0
                            eDataType, nDTSize, nBufXSize);
1497
0
        }
1498
0
        if (poBlock)
1499
0
            poBlock->DropLock();
1500
0
    }
1501
1502
0
    return CE_None;
1503
0
}
1504
1505
/************************************************************************/
1506
/*                VRTProcessedDataset::IRasterIO()                      */
1507
/************************************************************************/
1508
1509
CPLErr VRTProcessedDataset::IRasterIO(
1510
    GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
1511
    void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
1512
    int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
1513
    GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
1514
0
{
1515
    // Try to pass the request to the most appropriate overview dataset.
1516
0
    if (nBufXSize < nXSize && nBufYSize < nYSize)
1517
0
    {
1518
0
        int bTried = FALSE;
1519
0
        const CPLErr eErr = TryOverviewRasterIO(
1520
0
            eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1521
0
            eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
1522
0
            nBandSpace, psExtraArg, &bTried);
1523
0
        if (bTried)
1524
0
            return eErr;
1525
0
    }
1526
1527
    // Optimize reading of all bands at nominal resolution for BIP-like or
1528
    // BSQ-like buffer spacing.
1529
0
    if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
1530
0
        nBandCount == nBands)
1531
0
    {
1532
0
        const auto IsSequentialBandMap = [panBandMap, nBandCount]()
1533
0
        {
1534
0
            for (int i = 0; i < nBandCount; ++i)
1535
0
            {
1536
0
                if (panBandMap[i] != i + 1)
1537
0
                {
1538
0
                    return false;
1539
0
                }
1540
0
            }
1541
0
            return true;
1542
0
        };
1543
1544
0
        const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
1545
0
        const bool bIsBIPLike =
1546
0
            nBandSpace == nBufTypeSize && nPixelSpace == nBandSpace * nBands &&
1547
0
            nLineSpace >= nPixelSpace * nBufXSize && IsSequentialBandMap();
1548
0
        const bool bIsBSQLike = nPixelSpace == nBufTypeSize &&
1549
0
                                nLineSpace >= nPixelSpace * nBufXSize &&
1550
0
                                nBandSpace >= nLineSpace * nBufYSize &&
1551
0
                                IsSequentialBandMap();
1552
0
        if (bIsBIPLike || bIsBSQLike)
1553
0
        {
1554
0
            GByte *pabyData = static_cast<GByte *>(pData);
1555
            // If acquiring the region of interest in a single time is going
1556
            // to consume too much RAM, split in halves.
1557
0
            if (m_nAllowedRAMUsage > 0 &&
1558
0
                static_cast<GIntBig>(nBufXSize) * nBufYSize >
1559
0
                    m_nAllowedRAMUsage / m_nWorkingBytesPerPixel)
1560
0
            {
1561
0
                if ((nBufXSize == nRasterXSize || nBufYSize >= nBufXSize) &&
1562
0
                    nBufYSize >= 2)
1563
0
                {
1564
0
                    GDALRasterIOExtraArg sArg;
1565
0
                    INIT_RASTERIO_EXTRA_ARG(sArg);
1566
0
                    const int nHalfHeight = nBufYSize / 2;
1567
1568
0
                    sArg.pfnProgress = GDALScaledProgress;
1569
0
                    sArg.pProgressData = GDALCreateScaledProgress(
1570
0
                        0, 0.5, psExtraArg->pfnProgress,
1571
0
                        psExtraArg->pProgressData);
1572
0
                    if (sArg.pProgressData == nullptr)
1573
0
                        sArg.pfnProgress = nullptr;
1574
0
                    bool bOK =
1575
0
                        IRasterIO(eRWFlag, nXOff, nYOff, nBufXSize, nHalfHeight,
1576
0
                                  pabyData, nBufXSize, nHalfHeight, eBufType,
1577
0
                                  nBandCount, panBandMap, nPixelSpace,
1578
0
                                  nLineSpace, nBandSpace, &sArg) == CE_None;
1579
0
                    GDALDestroyScaledProgress(sArg.pProgressData);
1580
1581
0
                    if (bOK)
1582
0
                    {
1583
0
                        sArg.pfnProgress = GDALScaledProgress;
1584
0
                        sArg.pProgressData = GDALCreateScaledProgress(
1585
0
                            0.5, 1, psExtraArg->pfnProgress,
1586
0
                            psExtraArg->pProgressData);
1587
0
                        if (sArg.pProgressData == nullptr)
1588
0
                            sArg.pfnProgress = nullptr;
1589
0
                        bOK = IRasterIO(eRWFlag, nXOff, nYOff + nHalfHeight,
1590
0
                                        nBufXSize, nBufYSize - nHalfHeight,
1591
0
                                        pabyData + nHalfHeight * nLineSpace,
1592
0
                                        nBufXSize, nBufYSize - nHalfHeight,
1593
0
                                        eBufType, nBandCount, panBandMap,
1594
0
                                        nPixelSpace, nLineSpace, nBandSpace,
1595
0
                                        &sArg) == CE_None;
1596
0
                        GDALDestroyScaledProgress(sArg.pProgressData);
1597
0
                    }
1598
0
                    return bOK ? CE_None : CE_Failure;
1599
0
                }
1600
0
                else if (nBufXSize >= 2)
1601
0
                {
1602
0
                    GDALRasterIOExtraArg sArg;
1603
0
                    INIT_RASTERIO_EXTRA_ARG(sArg);
1604
0
                    const int nHalfWidth = nBufXSize / 2;
1605
1606
0
                    sArg.pfnProgress = GDALScaledProgress;
1607
0
                    sArg.pProgressData = GDALCreateScaledProgress(
1608
0
                        0, 0.5, psExtraArg->pfnProgress,
1609
0
                        psExtraArg->pProgressData);
1610
0
                    if (sArg.pProgressData == nullptr)
1611
0
                        sArg.pfnProgress = nullptr;
1612
0
                    bool bOK =
1613
0
                        IRasterIO(eRWFlag, nXOff, nYOff, nHalfWidth, nBufYSize,
1614
0
                                  pabyData, nHalfWidth, nBufYSize, eBufType,
1615
0
                                  nBandCount, panBandMap, nPixelSpace,
1616
0
                                  nLineSpace, nBandSpace, &sArg) == CE_None;
1617
0
                    GDALDestroyScaledProgress(sArg.pProgressData);
1618
1619
0
                    if (bOK)
1620
0
                    {
1621
0
                        sArg.pfnProgress = GDALScaledProgress;
1622
0
                        sArg.pProgressData = GDALCreateScaledProgress(
1623
0
                            0.5, 1, psExtraArg->pfnProgress,
1624
0
                            psExtraArg->pProgressData);
1625
0
                        if (sArg.pProgressData == nullptr)
1626
0
                            sArg.pfnProgress = nullptr;
1627
0
                        bOK = IRasterIO(eRWFlag, nXOff + nHalfWidth, nYOff,
1628
0
                                        nBufXSize - nHalfWidth, nBufYSize,
1629
0
                                        pabyData + nHalfWidth * nPixelSpace,
1630
0
                                        nBufXSize - nHalfWidth, nBufYSize,
1631
0
                                        eBufType, nBandCount, panBandMap,
1632
0
                                        nPixelSpace, nLineSpace, nBandSpace,
1633
0
                                        &sArg) == CE_None;
1634
0
                        GDALDestroyScaledProgress(sArg.pProgressData);
1635
0
                    }
1636
0
                    return bOK ? CE_None : CE_Failure;
1637
0
                }
1638
0
            }
1639
1640
0
            if (!ProcessRegion(nXOff, nYOff, nBufXSize, nBufYSize,
1641
0
                               psExtraArg->pfnProgress,
1642
0
                               psExtraArg->pProgressData))
1643
0
            {
1644
0
                return CE_Failure;
1645
0
            }
1646
0
            const auto eLastDT = m_aoSteps.back().eOutDT;
1647
0
            const int nLastDTSize = GDALGetDataTypeSizeBytes(eLastDT);
1648
0
            if (bIsBIPLike)
1649
0
            {
1650
0
                for (int iY = 0; iY < nBufYSize; ++iY)
1651
0
                {
1652
0
                    GDALCopyWords64(
1653
0
                        m_abyInput.data() + static_cast<size_t>(iY) * nBands *
1654
0
                                                nBufXSize * nLastDTSize,
1655
0
                        eLastDT, nLastDTSize, pabyData + iY * nLineSpace,
1656
0
                        eBufType, GDALGetDataTypeSizeBytes(eBufType),
1657
0
                        static_cast<size_t>(nBufXSize) * nBands);
1658
0
                }
1659
0
            }
1660
0
            else
1661
0
            {
1662
0
                CPLAssert(bIsBSQLike);
1663
0
                for (int iBand = 0; iBand < nBands; ++iBand)
1664
0
                {
1665
0
                    for (int iY = 0; iY < nBufYSize; ++iY)
1666
0
                    {
1667
0
                        GDALCopyWords64(
1668
0
                            m_abyInput.data() +
1669
0
                                (static_cast<size_t>(iY) * nBands * nBufXSize +
1670
0
                                 iBand) *
1671
0
                                    nLastDTSize,
1672
0
                            eLastDT, nLastDTSize * nBands,
1673
0
                            pabyData + iBand * nBandSpace + iY * nLineSpace,
1674
0
                            eBufType, nBufTypeSize, nBufXSize);
1675
0
                    }
1676
0
                }
1677
0
            }
1678
0
            return CE_None;
1679
0
        }
1680
0
    }
1681
1682
0
    return VRTDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
1683
0
                                 nBufXSize, nBufYSize, eBufType, nBandCount,
1684
0
                                 panBandMap, nPixelSpace, nLineSpace,
1685
0
                                 nBandSpace, psExtraArg);
1686
0
}
1687
1688
/*! @endcond */
1689
1690
/************************************************************************/
1691
/*                GDALVRTRegisterProcessedDatasetFunc()                 */
1692
/************************************************************************/
1693
1694
/** Register a function to be used by VRTProcessedDataset.
1695
1696
 An example of content for pszXMLMetadata is:
1697
 \verbatim
1698
  <ProcessedDatasetFunctionArgumentsList>
1699
     <Argument name='src_nodata' type='double' description='Override input nodata value'/>
1700
     <Argument name='dst_nodata' type='double' description='Override output nodata value'/>
1701
     <Argument name='replacement_nodata' description='value to substitute to a valid computed value that would be nodata' type='double'/>
1702
     <Argument name='dst_intended_datatype' type='string' description='Intented datatype of output (which might be different than the working data type)'/>
1703
     <Argument name='coefficients_{band}' description='Comma-separated coefficients for combining bands. First one is constant term' type='double_list' required='true'/>
1704
  </ProcessedDatasetFunctionArgumentsList>
1705
 \endverbatim
1706
1707
 @param pszFuncName Function name. Must be unique and not null.
1708
 @param pUserData User data. May be nullptr. Must remain valid during the
1709
                  lifetime of GDAL.
1710
 @param pszXMLMetadata XML metadata describing the function arguments. May be
1711
                       nullptr if there are no arguments.
1712
 @param eRequestedInputDT If the pfnProcess callback only supports a single
1713
                          data type, it should be specified in this parameter.
1714
                          Otherwise set it to GDT_Unknown.
1715
 @param paeSupportedInputDT List of supported input data types. May be nullptr
1716
                            if all are supported or if eRequestedInputDT is
1717
                            set to a non GDT_Unknown value.
1718
 @param nSupportedInputDTSize Size of paeSupportedInputDT
1719
 @param panSupportedInputBandCount List of supported band count. May be nullptr
1720
                                   if any source band count is supported.
1721
 @param nSupportedInputBandCountSize Size of panSupportedInputBandCount
1722
 @param pfnInit Initialization function called when a VRTProcessedDataset
1723
                step uses the register function. This initialization function
1724
                will return the output data type, output band count and
1725
                potentially initialize a working structure, typically parsing
1726
                arguments. May be nullptr.
1727
                If not specified, it will be assumed that the input and output
1728
                data types are the same, and that the input number of bands
1729
                and output number of bands are the same.
1730
 @param pfnFree Free function that will free the working structure allocated
1731
                by pfnInit. May be nullptr.
1732
 @param pfnProcess Processing function called to compute pixel values. Must
1733
                   not be nullptr.
1734
 @param papszOptions Unused currently. Must be nullptr.
1735
 @return CE_None in case of success, error otherwise.
1736
 @since 3.9
1737
 */
1738
CPLErr GDALVRTRegisterProcessedDatasetFunc(
1739
    const char *pszFuncName, void *pUserData, const char *pszXMLMetadata,
1740
    GDALDataType eRequestedInputDT, const GDALDataType *paeSupportedInputDT,
1741
    size_t nSupportedInputDTSize, const int *panSupportedInputBandCount,
1742
    size_t nSupportedInputBandCountSize,
1743
    GDALVRTProcessedDatasetFuncInit pfnInit,
1744
    GDALVRTProcessedDatasetFuncFree pfnFree,
1745
    GDALVRTProcessedDatasetFuncProcess pfnProcess,
1746
    CPL_UNUSED CSLConstList papszOptions)
1747
10
{
1748
10
    if (pszFuncName == nullptr || pszFuncName[0] == '\0')
1749
0
    {
1750
0
        CPLError(CE_Failure, CPLE_AppDefined,
1751
0
                 "pszFuncName should be non-empty");
1752
0
        return CE_Failure;
1753
0
    }
1754
1755
10
    auto &oMap = GetGlobalMapProcessedDatasetFunc();
1756
10
    if (oMap.find(pszFuncName) != oMap.end())
1757
0
    {
1758
0
        CPLError(CE_Failure, CPLE_AppDefined, "%s already registered",
1759
0
                 pszFuncName);
1760
0
        return CE_Failure;
1761
0
    }
1762
1763
10
    if (!pfnProcess)
1764
0
    {
1765
0
        CPLError(CE_Failure, CPLE_AppDefined, "pfnProcess should not be null");
1766
0
        return CE_Failure;
1767
0
    }
1768
1769
10
    VRTProcessedDatasetFunc oFunc;
1770
10
    oFunc.osFuncName = pszFuncName;
1771
10
    oFunc.pUserData = pUserData;
1772
10
    if (pszXMLMetadata)
1773
10
    {
1774
10
        oFunc.bMetadataSpecified = true;
1775
10
        auto psTree = CPLXMLTreeCloser(CPLParseXMLString(pszXMLMetadata));
1776
10
        if (!psTree)
1777
0
        {
1778
0
            CPLError(CE_Failure, CPLE_AppDefined,
1779
0
                     "Cannot parse pszXMLMetadata=%s for %s", pszXMLMetadata,
1780
0
                     pszFuncName);
1781
0
            return CE_Failure;
1782
0
        }
1783
10
        const CPLXMLNode *psRoot = CPLGetXMLNode(
1784
10
            psTree.get(), "=ProcessedDatasetFunctionArgumentsList");
1785
10
        if (!psRoot)
1786
0
        {
1787
0
            CPLError(CE_Failure, CPLE_AppDefined,
1788
0
                     "No root ProcessedDatasetFunctionArgumentsList element in "
1789
0
                     "pszXMLMetadata=%s for %s",
1790
0
                     pszXMLMetadata, pszFuncName);
1791
0
            return CE_Failure;
1792
0
        }
1793
76
        for (const CPLXMLNode *psIter = psRoot->psChild; psIter;
1794
66
             psIter = psIter->psNext)
1795
66
        {
1796
66
            if (psIter->eType == CXT_Element &&
1797
66
                strcmp(psIter->pszValue, "Argument") == 0)
1798
66
            {
1799
66
                const char *pszName = CPLGetXMLValue(psIter, "name", nullptr);
1800
66
                if (!pszName)
1801
0
                {
1802
0
                    CPLError(CE_Failure, CPLE_AppDefined,
1803
0
                             "Missing Argument.name attribute in "
1804
0
                             "pszXMLMetadata=%s for %s",
1805
0
                             pszXMLMetadata, pszFuncName);
1806
0
                    return CE_Failure;
1807
0
                }
1808
66
                const char *pszType = CPLGetXMLValue(psIter, "type", nullptr);
1809
66
                if (!pszType)
1810
0
                {
1811
0
                    CPLError(CE_Failure, CPLE_AppDefined,
1812
0
                             "Missing Argument.type attribute in "
1813
0
                             "pszXMLMetadata=%s for %s",
1814
0
                             pszXMLMetadata, pszFuncName);
1815
0
                    return CE_Failure;
1816
0
                }
1817
66
                if (strcmp(pszType, "constant") == 0)
1818
0
                {
1819
0
                    const char *pszValue =
1820
0
                        CPLGetXMLValue(psIter, "value", nullptr);
1821
0
                    if (!pszValue)
1822
0
                    {
1823
0
                        CPLError(CE_Failure, CPLE_AppDefined,
1824
0
                                 "Missing Argument.value attribute in "
1825
0
                                 "pszXMLMetadata=%s for %s",
1826
0
                                 pszXMLMetadata, pszFuncName);
1827
0
                        return CE_Failure;
1828
0
                    }
1829
0
                    oFunc.oMapConstantArguments[CPLString(pszName).tolower()] =
1830
0
                        pszValue;
1831
0
                }
1832
66
                else if (strcmp(pszType, "builtin") == 0)
1833
0
                {
1834
0
                    if (EQUAL(pszName, "nodata") ||
1835
0
                        EQUAL(pszName, "offset_{band}") ||
1836
0
                        EQUAL(pszName, "scale_{band}"))
1837
0
                    {
1838
0
                        oFunc.oSetBuiltinArguments.insert(
1839
0
                            CPLString(pszName).tolower());
1840
0
                    }
1841
0
                    else
1842
0
                    {
1843
0
                        CPLError(CE_Failure, CPLE_NotSupported,
1844
0
                                 "Unsupported builtin parameter name %s in "
1845
0
                                 "pszXMLMetadata=%s for %s. Only nodata, "
1846
0
                                 "offset_{band} and scale_{band} are supported",
1847
0
                                 pszName, pszXMLMetadata, pszFuncName);
1848
0
                        return CE_Failure;
1849
0
                    }
1850
0
                }
1851
66
                else if (strcmp(pszType, "boolean") == 0 ||
1852
66
                         strcmp(pszType, "string") == 0 ||
1853
66
                         strcmp(pszType, "integer") == 0 ||
1854
66
                         strcmp(pszType, "double") == 0 ||
1855
66
                         strcmp(pszType, "double_list") == 0)
1856
66
                {
1857
66
                    VRTProcessedDatasetFunc::OtherArgument otherArgument;
1858
66
                    otherArgument.bRequired = CPLTestBool(
1859
66
                        CPLGetXMLValue(psIter, "required", "false"));
1860
66
                    otherArgument.osType = pszType;
1861
66
                    oFunc.oOtherArguments[CPLString(pszName).tolower()] =
1862
66
                        std::move(otherArgument);
1863
66
                }
1864
0
                else
1865
0
                {
1866
0
                    CPLError(CE_Failure, CPLE_NotSupported,
1867
0
                             "Unsupported type for parameter %s in "
1868
0
                             "pszXMLMetadata=%s for %s. Only boolean, string, "
1869
0
                             "integer, double and double_list are supported",
1870
0
                             pszName, pszXMLMetadata, pszFuncName);
1871
0
                    return CE_Failure;
1872
0
                }
1873
66
            }
1874
66
        }
1875
10
    }
1876
10
    oFunc.eRequestedInputDT = eRequestedInputDT;
1877
10
    if (nSupportedInputDTSize)
1878
0
    {
1879
0
        oFunc.aeSupportedInputDT.insert(
1880
0
            oFunc.aeSupportedInputDT.end(), paeSupportedInputDT,
1881
0
            paeSupportedInputDT + nSupportedInputDTSize);
1882
0
    }
1883
10
    if (nSupportedInputBandCountSize)
1884
0
    {
1885
0
        oFunc.anSupportedInputBandCount.insert(
1886
0
            oFunc.anSupportedInputBandCount.end(), panSupportedInputBandCount,
1887
0
            panSupportedInputBandCount + nSupportedInputBandCountSize);
1888
0
    }
1889
10
    oFunc.pfnInit = pfnInit;
1890
10
    oFunc.pfnFree = pfnFree;
1891
10
    oFunc.pfnProcess = pfnProcess;
1892
1893
10
    oMap[pszFuncName] = std::move(oFunc);
1894
1895
10
    return CE_None;
1896
10
}