Coverage Report

Created: 2025-06-13 06:29

/src/gdal/frmts/vrt/vrtsources.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  Virtual GDAL Datasets
4
 * Purpose:  Implementation of VRTSimpleSource, VRTFuncSource and
5
 *           VRTAveragedSource.
6
 * Author:   Frank Warmerdam <warmerdam@pobox.com>
7
 *
8
 ******************************************************************************
9
 * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.com>
10
 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
11
 *
12
 * SPDX-License-Identifier: MIT
13
 ****************************************************************************/
14
15
#include "gdal_vrt.h"
16
#include "vrtdataset.h"
17
18
#include <cassert>
19
#include <climits>
20
#include <cmath>
21
#include <cstddef>
22
#include <cstdio>
23
#include <cstdlib>
24
#include <cstring>
25
#include <algorithm>
26
#include <limits>
27
#include <string>
28
29
#include "cpl_conv.h"
30
#include "cpl_error.h"
31
#include "cpl_hash_set.h"
32
#include "cpl_minixml.h"
33
#include "cpl_progress.h"
34
#include "cpl_string.h"
35
#include "cpl_vsi.h"
36
#include "gdal.h"
37
#include "gdal_priv.h"
38
#include "gdal_proxy.h"
39
#include "gdal_priv_templates.hpp"
40
41
/*! @cond Doxygen_Suppress */
42
43
// #define DEBUG_VERBOSE 1
44
45
/************************************************************************/
46
/* ==================================================================== */
47
/*                             VRTSource                                */
48
/* ==================================================================== */
49
/************************************************************************/
50
51
VRTSource::~VRTSource()
52
0
{
53
0
}
54
55
/************************************************************************/
56
/*                             GetFileList()                            */
57
/************************************************************************/
58
59
void VRTSource::GetFileList(char *** /* ppapszFileList */, int * /* pnSize */,
60
                            int * /* pnMaxSize */, CPLHashSet * /* hSetFiles */)
61
0
{
62
0
}
63
64
/************************************************************************/
65
/* ==================================================================== */
66
/*                          VRTSimpleSource                             */
67
/* ==================================================================== */
68
/************************************************************************/
69
70
/************************************************************************/
71
/*                          VRTSimpleSource()                           */
72
/************************************************************************/
73
74
0
VRTSimpleSource::VRTSimpleSource() = default;
75
76
/************************************************************************/
77
/*                          VRTSimpleSource()                           */
78
/************************************************************************/
79
80
VRTSimpleSource::VRTSimpleSource(const VRTSimpleSource *poSrcSource,
81
                                 double dfXDstRatio, double dfYDstRatio)
82
0
    : m_poMapSharedSources(poSrcSource->m_poMapSharedSources),
83
0
      m_poRasterBand(poSrcSource->m_poRasterBand),
84
0
      m_poMaskBandMainBand(poSrcSource->m_poMaskBandMainBand),
85
0
      m_aosOpenOptionsOri(poSrcSource->m_aosOpenOptionsOri),
86
0
      m_aosOpenOptions(poSrcSource->m_aosOpenOptions),
87
0
      m_bSrcDSNameFromVRT(poSrcSource->m_bSrcDSNameFromVRT),
88
0
      m_nBand(poSrcSource->m_nBand),
89
0
      m_bGetMaskBand(poSrcSource->m_bGetMaskBand),
90
0
      m_dfSrcXOff(poSrcSource->m_dfSrcXOff),
91
0
      m_dfSrcYOff(poSrcSource->m_dfSrcYOff),
92
0
      m_dfSrcXSize(poSrcSource->m_dfSrcXSize),
93
0
      m_dfSrcYSize(poSrcSource->m_dfSrcYSize),
94
0
      m_nMaxValue(poSrcSource->m_nMaxValue), m_bRelativeToVRTOri(-1),
95
0
      m_nExplicitSharedStatus(poSrcSource->m_nExplicitSharedStatus),
96
0
      m_osSrcDSName(poSrcSource->m_osSrcDSName),
97
0
      m_bDropRefOnSrcBand(poSrcSource->m_bDropRefOnSrcBand)
98
0
{
99
0
    if (!poSrcSource->IsSrcWinSet() && !poSrcSource->IsDstWinSet() &&
100
0
        (dfXDstRatio != 1.0 || dfYDstRatio != 1.0))
101
0
    {
102
0
        auto l_band = GetRasterBand();
103
0
        if (l_band)
104
0
        {
105
0
            m_dfSrcXOff = 0;
106
0
            m_dfSrcYOff = 0;
107
0
            m_dfSrcXSize = l_band->GetXSize();
108
0
            m_dfSrcYSize = l_band->GetYSize();
109
0
            m_dfDstXOff = 0;
110
0
            m_dfDstYOff = 0;
111
0
            m_dfDstXSize = l_band->GetXSize() * dfXDstRatio;
112
0
            m_dfDstYSize = l_band->GetYSize() * dfYDstRatio;
113
0
        }
114
0
    }
115
0
    else if (poSrcSource->IsDstWinSet())
116
0
    {
117
0
        m_dfDstXOff = poSrcSource->m_dfDstXOff * dfXDstRatio;
118
0
        m_dfDstYOff = poSrcSource->m_dfDstYOff * dfYDstRatio;
119
0
        m_dfDstXSize = poSrcSource->m_dfDstXSize * dfXDstRatio;
120
0
        m_dfDstYSize = poSrcSource->m_dfDstYSize * dfYDstRatio;
121
0
    }
122
0
}
123
124
/************************************************************************/
125
/*                          ~VRTSimpleSource()                          */
126
/************************************************************************/
127
128
VRTSimpleSource::~VRTSimpleSource()
129
130
0
{
131
0
    if (!m_bDropRefOnSrcBand)
132
0
        return;
133
134
0
    if (m_poMaskBandMainBand != nullptr)
135
0
    {
136
0
        if (m_poMaskBandMainBand->GetDataset() != nullptr)
137
0
        {
138
0
            m_poMaskBandMainBand->GetDataset()->ReleaseRef();
139
0
        }
140
0
    }
141
0
    else if (m_poRasterBand != nullptr &&
142
0
             m_poRasterBand->GetDataset() != nullptr)
143
0
    {
144
0
        m_poRasterBand->GetDataset()->ReleaseRef();
145
0
    }
146
0
}
147
148
/************************************************************************/
149
/*                           GetTypeStatic()                            */
150
/************************************************************************/
151
152
const char *VRTSimpleSource::GetTypeStatic()
153
0
{
154
0
    static const char *TYPE = "SimpleSource";
155
0
    return TYPE;
156
0
}
157
158
/************************************************************************/
159
/*                            GetType()                                 */
160
/************************************************************************/
161
162
const char *VRTSimpleSource::GetType() const
163
0
{
164
0
    return GetTypeStatic();
165
0
}
166
167
/************************************************************************/
168
/*                           FlushCache()                               */
169
/************************************************************************/
170
171
CPLErr VRTSimpleSource::FlushCache(bool bAtClosing)
172
173
0
{
174
0
    if (m_poMaskBandMainBand != nullptr)
175
0
    {
176
0
        return m_poMaskBandMainBand->FlushCache(bAtClosing);
177
0
    }
178
0
    else if (m_poRasterBand != nullptr)
179
0
    {
180
0
        return m_poRasterBand->FlushCache(bAtClosing);
181
0
    }
182
0
    return CE_None;
183
0
}
184
185
/************************************************************************/
186
/*                    UnsetPreservedRelativeFilenames()                 */
187
/************************************************************************/
188
189
void VRTSimpleSource::UnsetPreservedRelativeFilenames()
190
0
{
191
0
    if (m_bRelativeToVRTOri &&
192
0
        !STARTS_WITH(m_osSourceFileNameOri.c_str(), "http://") &&
193
0
        !STARTS_WITH(m_osSourceFileNameOri.c_str(), "https://"))
194
0
    {
195
0
        m_bRelativeToVRTOri = -1;
196
0
        m_osSourceFileNameOri = "";
197
0
    }
198
0
}
199
200
/************************************************************************/
201
/*                             SetSrcBand()                             */
202
/************************************************************************/
203
204
void VRTSimpleSource::SetSrcBand(const char *pszFilename, int nBand)
205
206
0
{
207
0
    m_nBand = nBand;
208
0
    m_osSrcDSName = pszFilename;
209
0
}
210
211
/************************************************************************/
212
/*                             SetSrcBand()                             */
213
/************************************************************************/
214
215
void VRTSimpleSource::SetSrcBand(GDALRasterBand *poNewSrcBand)
216
217
0
{
218
0
    m_poRasterBand = poNewSrcBand;
219
0
    m_nBand = m_poRasterBand->GetBand();
220
0
    auto poDS = poNewSrcBand->GetDataset();
221
0
    if (poDS != nullptr)
222
0
    {
223
0
        m_osSrcDSName = poDS->GetDescription();
224
0
        m_aosOpenOptions = CSLDuplicate(poDS->GetOpenOptions());
225
0
        m_aosOpenOptionsOri = m_aosOpenOptions;
226
0
    }
227
0
}
228
229
/************************************************************************/
230
/*                      SetSourceDatasetName()                          */
231
/************************************************************************/
232
233
void VRTSimpleSource::SetSourceDatasetName(const char *pszFilename,
234
                                           bool bRelativeToVRT)
235
0
{
236
0
    CPLAssert(m_nBand >= 0);
237
0
    m_osSrcDSName = pszFilename;
238
0
    m_osSourceFileNameOri = pszFilename;
239
0
    m_bRelativeToVRTOri = bRelativeToVRT;
240
0
}
241
242
/************************************************************************/
243
/*                          SetSrcMaskBand()                            */
244
/************************************************************************/
245
246
// poSrcBand is not the mask band, but the band from which the mask band is
247
// taken.
248
void VRTSimpleSource::SetSrcMaskBand(GDALRasterBand *poNewSrcBand)
249
250
0
{
251
0
    m_poRasterBand = poNewSrcBand->GetMaskBand();
252
0
    m_poMaskBandMainBand = poNewSrcBand;
253
0
    m_nBand = poNewSrcBand->GetBand();
254
0
    auto poDS = poNewSrcBand->GetDataset();
255
0
    if (poDS != nullptr)
256
0
    {
257
0
        m_osSrcDSName = poDS->GetDescription();
258
0
        m_aosOpenOptions = CSLDuplicate(poDS->GetOpenOptions());
259
0
        m_aosOpenOptionsOri = m_aosOpenOptions;
260
0
    }
261
0
    m_bGetMaskBand = true;
262
0
}
263
264
/************************************************************************/
265
/*                         RoundIfCloseToInt()                          */
266
/************************************************************************/
267
268
static double RoundIfCloseToInt(double dfValue)
269
0
{
270
0
    double dfClosestInt = floor(dfValue + 0.5);
271
0
    return (fabs(dfValue - dfClosestInt) < 1e-3) ? dfClosestInt : dfValue;
272
0
}
273
274
/************************************************************************/
275
/*                            SetSrcWindow()                            */
276
/************************************************************************/
277
278
void VRTSimpleSource::SetSrcWindow(double dfNewXOff, double dfNewYOff,
279
                                   double dfNewXSize, double dfNewYSize)
280
281
0
{
282
0
    m_dfSrcXOff = RoundIfCloseToInt(dfNewXOff);
283
0
    m_dfSrcYOff = RoundIfCloseToInt(dfNewYOff);
284
0
    m_dfSrcXSize = RoundIfCloseToInt(dfNewXSize);
285
0
    m_dfSrcYSize = RoundIfCloseToInt(dfNewYSize);
286
0
}
287
288
/************************************************************************/
289
/*                            SetDstWindow()                            */
290
/************************************************************************/
291
292
void VRTSimpleSource::SetDstWindow(double dfNewXOff, double dfNewYOff,
293
                                   double dfNewXSize, double dfNewYSize)
294
295
0
{
296
0
    m_dfDstXOff = RoundIfCloseToInt(dfNewXOff);
297
0
    m_dfDstYOff = RoundIfCloseToInt(dfNewYOff);
298
0
    m_dfDstXSize = RoundIfCloseToInt(dfNewXSize);
299
0
    m_dfDstYSize = RoundIfCloseToInt(dfNewYSize);
300
0
}
301
302
/************************************************************************/
303
/*                            GetDstWindow()                            */
304
/************************************************************************/
305
306
void VRTSimpleSource::GetDstWindow(double &dfDstXOff, double &dfDstYOff,
307
                                   double &dfDstXSize, double &dfDstYSize) const
308
0
{
309
0
    dfDstXOff = m_dfDstXOff;
310
0
    dfDstYOff = m_dfDstYOff;
311
0
    dfDstXSize = m_dfDstXSize;
312
0
    dfDstYSize = m_dfDstYSize;
313
0
}
314
315
/************************************************************************/
316
/*                        DstWindowIntersects()                         */
317
/************************************************************************/
318
319
bool VRTSimpleSource::DstWindowIntersects(double dfXOff, double dfYOff,
320
                                          double dfXSize, double dfYSize) const
321
0
{
322
0
    return IsDstWinSet() && m_dfDstXOff + m_dfDstXSize > dfXOff &&
323
0
           m_dfDstYOff + m_dfDstYSize > dfYOff &&
324
0
           m_dfDstXOff < dfXOff + dfXSize && m_dfDstYOff < dfYOff + dfYSize;
325
0
}
326
327
/************************************************************************/
328
/*                            IsSlowSource()                            */
329
/************************************************************************/
330
331
static bool IsSlowSource(const char *pszSrcName)
332
0
{
333
0
    return strstr(pszSrcName, "/vsicurl/http") != nullptr ||
334
0
           strstr(pszSrcName, "/vsicurl/ftp") != nullptr ||
335
0
           (strstr(pszSrcName, "/vsicurl?") != nullptr &&
336
0
            strstr(pszSrcName, "&url=http") != nullptr);
337
0
}
338
339
/************************************************************************/
340
/*                         AddSourceFilenameNode()                      */
341
/************************************************************************/
342
343
void VRTSimpleSource::AddSourceFilenameNode(const char *pszVRTPath,
344
                                            CPLXMLNode *psSrc)
345
0
{
346
347
0
    VSIStatBufL sStat;
348
0
    int bRelativeToVRT = FALSE;  // TODO(schwehr): Make this a bool?
349
0
    std::string osSourceFilename;
350
351
0
    if (m_bRelativeToVRTOri >= 0)
352
0
    {
353
0
        osSourceFilename = m_osSourceFileNameOri;
354
0
        bRelativeToVRT = m_bRelativeToVRTOri;
355
0
    }
356
0
    else if (IsSlowSource(m_osSrcDSName))
357
0
    {
358
        // Testing the existence of remote resources can be excruciating
359
        // slow, so let's just suppose they exist.
360
0
        osSourceFilename = m_osSrcDSName;
361
0
        bRelativeToVRT = FALSE;
362
0
    }
363
    // If this isn't actually a file, don't even try to know if it is a
364
    // relative path. It can't be !, and unfortunately CPLIsFilenameRelative()
365
    // can only work with strings that are filenames To be clear
366
    // NITF_TOC_ENTRY:CADRG_JOG-A_250K_1_0:some_path isn't a relative file
367
    // path.
368
0
    else if (VSIStatExL(m_osSrcDSName, &sStat, VSI_STAT_EXISTS_FLAG) != 0)
369
0
    {
370
0
        osSourceFilename = m_osSrcDSName;
371
0
        bRelativeToVRT = FALSE;
372
373
        // Try subdatasetinfo API first
374
        // Note: this will become the only branch when subdatasetinfo will become
375
        //       available for NITF_IM, RASTERLITE and TILEDB
376
0
        const auto oSubDSInfo{GDALGetSubdatasetInfo(osSourceFilename.c_str())};
377
0
        if (oSubDSInfo && !oSubDSInfo->GetPathComponent().empty())
378
0
        {
379
0
            auto path{oSubDSInfo->GetPathComponent()};
380
0
            std::string relPath{CPLExtractRelativePath(pszVRTPath, path.c_str(),
381
0
                                                       &bRelativeToVRT)};
382
0
            osSourceFilename = oSubDSInfo->ModifyPathComponent(relPath);
383
0
            GDALDestroySubdatasetInfo(oSubDSInfo);
384
0
        }
385
0
        else
386
0
        {
387
0
            for (const char *pszSyntax :
388
0
                 GDALDataset::apszSpecialSubDatasetSyntax)
389
0
            {
390
0
                CPLString osPrefix(pszSyntax);
391
0
                osPrefix.resize(strchr(pszSyntax, ':') - pszSyntax + 1);
392
0
                if (pszSyntax[osPrefix.size()] == '"')
393
0
                    osPrefix += '"';
394
0
                if (EQUALN(osSourceFilename.c_str(), osPrefix, osPrefix.size()))
395
0
                {
396
0
                    if (STARTS_WITH_CI(pszSyntax + osPrefix.size(), "{ANY}"))
397
0
                    {
398
0
                        const char *pszLastPart =
399
0
                            strrchr(osSourceFilename.c_str(), ':') + 1;
400
                        // CSV:z:/foo.xyz
401
0
                        if ((pszLastPart[0] == '/' || pszLastPart[0] == '\\') &&
402
0
                            pszLastPart - osSourceFilename.c_str() >= 3 &&
403
0
                            pszLastPart[-3] == ':')
404
0
                            pszLastPart -= 2;
405
0
                        CPLString osPrefixFilename(osSourceFilename);
406
0
                        osPrefixFilename.resize(pszLastPart -
407
0
                                                osSourceFilename.c_str());
408
0
                        osSourceFilename = CPLExtractRelativePath(
409
0
                            pszVRTPath, pszLastPart, &bRelativeToVRT);
410
0
                        osSourceFilename = osPrefixFilename + osSourceFilename;
411
0
                    }
412
0
                    else if (STARTS_WITH_CI(pszSyntax + osPrefix.size(),
413
0
                                            "{FILENAME}"))
414
0
                    {
415
0
                        CPLString osFilename(osSourceFilename.c_str() +
416
0
                                             osPrefix.size());
417
0
                        size_t nPos = 0;
418
0
                        if (osFilename.size() >= 3 && osFilename[1] == ':' &&
419
0
                            (osFilename[2] == '\\' || osFilename[2] == '/'))
420
0
                            nPos = 2;
421
0
                        nPos = osFilename.find(
422
0
                            pszSyntax[osPrefix.size() + strlen("{FILENAME}")],
423
0
                            nPos);
424
0
                        if (nPos != std::string::npos)
425
0
                        {
426
0
                            const CPLString osSuffix = osFilename.substr(nPos);
427
0
                            osFilename.resize(nPos);
428
0
                            osSourceFilename = CPLExtractRelativePath(
429
0
                                pszVRTPath, osFilename, &bRelativeToVRT);
430
0
                            osSourceFilename =
431
0
                                osPrefix + osSourceFilename + osSuffix;
432
0
                        }
433
0
                    }
434
0
                    break;
435
0
                }
436
0
            }
437
0
        }
438
0
    }
439
0
    else
440
0
    {
441
0
        std::string osVRTFilename = pszVRTPath;
442
0
        std::string osSourceDataset = m_osSrcDSName;
443
0
        char *pszCurDir = CPLGetCurrentDir();
444
0
        if (CPLIsFilenameRelative(osSourceDataset.c_str()) &&
445
0
            !CPLIsFilenameRelative(osVRTFilename.c_str()) &&
446
0
            pszCurDir != nullptr)
447
0
        {
448
0
            osSourceDataset = CPLFormFilenameSafe(
449
0
                pszCurDir, osSourceDataset.c_str(), nullptr);
450
0
        }
451
0
        else if (!CPLIsFilenameRelative(osSourceDataset.c_str()) &&
452
0
                 CPLIsFilenameRelative(osVRTFilename.c_str()) &&
453
0
                 pszCurDir != nullptr)
454
0
        {
455
0
            osVRTFilename =
456
0
                CPLFormFilenameSafe(pszCurDir, osVRTFilename.c_str(), nullptr);
457
0
        }
458
0
        CPLFree(pszCurDir);
459
0
        osSourceFilename = CPLExtractRelativePath(
460
0
            osVRTFilename.c_str(), osSourceDataset.c_str(), &bRelativeToVRT);
461
0
    }
462
463
0
    CPLSetXMLValue(psSrc, "SourceFilename", osSourceFilename.c_str());
464
465
0
    CPLCreateXMLNode(CPLCreateXMLNode(CPLGetXMLNode(psSrc, "SourceFilename"),
466
0
                                      CXT_Attribute, "relativeToVRT"),
467
0
                     CXT_Text, bRelativeToVRT ? "1" : "0");
468
469
    // Determine if we must write the shared attribute. The config option
470
    // will override the m_nExplicitSharedStatus value
471
0
    const char *pszShared = CPLGetConfigOption("VRT_SHARED_SOURCE", nullptr);
472
0
    if ((pszShared == nullptr && m_nExplicitSharedStatus == 0) ||
473
0
        (pszShared != nullptr && !CPLTestBool(pszShared)))
474
0
    {
475
0
        CPLCreateXMLNode(
476
0
            CPLCreateXMLNode(CPLGetXMLNode(psSrc, "SourceFilename"),
477
0
                             CXT_Attribute, "shared"),
478
0
            CXT_Text, "0");
479
0
    }
480
0
}
481
482
/************************************************************************/
483
/*                           SerializeToXML()                           */
484
/************************************************************************/
485
486
CPLXMLNode *VRTSimpleSource::SerializeToXML(const char *pszVRTPath)
487
488
0
{
489
0
    CPLXMLNode *const psSrc =
490
0
        CPLCreateXMLNode(nullptr, CXT_Element, GetTypeStatic());
491
492
0
    if (!m_osResampling.empty())
493
0
    {
494
0
        CPLCreateXMLNode(CPLCreateXMLNode(psSrc, CXT_Attribute, "resampling"),
495
0
                         CXT_Text, m_osResampling.c_str());
496
0
    }
497
498
0
    if (!m_osName.empty())
499
0
    {
500
0
        CPLAddXMLAttributeAndValue(psSrc, "name", m_osName);
501
0
    }
502
503
0
    if (m_bSrcDSNameFromVRT)
504
0
    {
505
0
        CPLAddXMLChild(psSrc, CPLParseXMLString(m_osSrcDSName.c_str()));
506
0
    }
507
0
    else
508
0
    {
509
0
        AddSourceFilenameNode(pszVRTPath, psSrc);
510
0
    }
511
512
0
    GDALSerializeOpenOptionsToXML(psSrc, m_aosOpenOptionsOri.List());
513
514
0
    if (m_bGetMaskBand)
515
0
        CPLSetXMLValue(psSrc, "SourceBand", CPLSPrintf("mask,%d", m_nBand));
516
0
    else
517
0
        CPLSetXMLValue(psSrc, "SourceBand", CPLSPrintf("%d", m_nBand));
518
519
    // TODO: in a later version, no longer emit SourceProperties, which
520
    // is no longer used by GDAL 3.4
521
0
    if (m_poRasterBand)
522
0
    {
523
        /* Write a few additional useful properties of the dataset */
524
        /* so that we can use a proxy dataset when re-opening. See XMLInit() */
525
        /* below */
526
0
        CPLSetXMLValue(psSrc, "SourceProperties.#RasterXSize",
527
0
                       CPLSPrintf("%d", m_poRasterBand->GetXSize()));
528
0
        CPLSetXMLValue(psSrc, "SourceProperties.#RasterYSize",
529
0
                       CPLSPrintf("%d", m_poRasterBand->GetYSize()));
530
0
        CPLSetXMLValue(
531
0
            psSrc, "SourceProperties.#DataType",
532
0
            GDALGetDataTypeName(m_poRasterBand->GetRasterDataType()));
533
534
0
        int nBlockXSize = 0;
535
0
        int nBlockYSize = 0;
536
0
        m_poRasterBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
537
538
0
        CPLSetXMLValue(psSrc, "SourceProperties.#BlockXSize",
539
0
                       CPLSPrintf("%d", nBlockXSize));
540
0
        CPLSetXMLValue(psSrc, "SourceProperties.#BlockYSize",
541
0
                       CPLSPrintf("%d", nBlockYSize));
542
0
    }
543
544
0
    if (IsSrcWinSet())
545
0
    {
546
0
        CPLSetXMLValue(psSrc, "SrcRect.#xOff",
547
0
                       CPLSPrintf("%.15g", m_dfSrcXOff));
548
0
        CPLSetXMLValue(psSrc, "SrcRect.#yOff",
549
0
                       CPLSPrintf("%.15g", m_dfSrcYOff));
550
0
        CPLSetXMLValue(psSrc, "SrcRect.#xSize",
551
0
                       CPLSPrintf("%.15g", m_dfSrcXSize));
552
0
        CPLSetXMLValue(psSrc, "SrcRect.#ySize",
553
0
                       CPLSPrintf("%.15g", m_dfSrcYSize));
554
0
    }
555
556
0
    if (IsDstWinSet())
557
0
    {
558
0
        CPLSetXMLValue(psSrc, "DstRect.#xOff",
559
0
                       CPLSPrintf("%.15g", m_dfDstXOff));
560
0
        CPLSetXMLValue(psSrc, "DstRect.#yOff",
561
0
                       CPLSPrintf("%.15g", m_dfDstYOff));
562
0
        CPLSetXMLValue(psSrc, "DstRect.#xSize",
563
0
                       CPLSPrintf("%.15g", m_dfDstXSize));
564
0
        CPLSetXMLValue(psSrc, "DstRect.#ySize",
565
0
                       CPLSPrintf("%.15g", m_dfDstYSize));
566
0
    }
567
568
0
    return psSrc;
569
0
}
570
571
/************************************************************************/
572
/*                              XMLInit()                               */
573
/************************************************************************/
574
575
CPLErr VRTSimpleSource::XMLInit(const CPLXMLNode *psSrc, const char *pszVRTPath,
576
                                VRTMapSharedResources &oMapSharedSources)
577
578
0
{
579
0
    m_poMapSharedSources = &oMapSharedSources;
580
581
0
    m_osResampling = CPLGetXMLValue(psSrc, "resampling", "");
582
0
    m_osName = CPLGetXMLValue(psSrc, "name", "");
583
584
    /* -------------------------------------------------------------------- */
585
    /*      Prepare filename.                                               */
586
    /* -------------------------------------------------------------------- */
587
0
    const CPLXMLNode *psSourceFileNameNode =
588
0
        CPLGetXMLNode(psSrc, "SourceFilename");
589
0
    const CPLXMLNode *psSourceVRTDataset = CPLGetXMLNode(psSrc, "VRTDataset");
590
0
    const char *pszFilename =
591
0
        psSourceFileNameNode ? CPLGetXMLValue(psSourceFileNameNode, nullptr, "")
592
0
                             : "";
593
594
0
    if (pszFilename[0] == '\0' && !psSourceVRTDataset)
595
0
    {
596
0
        CPLError(CE_Warning, CPLE_AppDefined,
597
0
                 "Missing <SourceFilename> or <VRTDataset> element in <%s>.",
598
0
                 psSrc->pszValue);
599
0
        return CE_Failure;
600
0
    }
601
602
    // Backup original filename and relativeToVRT so as to be able to
603
    // serialize them identically again (#5985)
604
0
    m_osSourceFileNameOri = pszFilename;
605
0
    if (pszFilename[0])
606
0
    {
607
0
        m_bRelativeToVRTOri =
608
0
            atoi(CPLGetXMLValue(psSourceFileNameNode, "relativetoVRT", "0"));
609
0
        const char *pszShared =
610
0
            CPLGetXMLValue(psSourceFileNameNode, "shared", nullptr);
611
0
        if (pszShared == nullptr)
612
0
        {
613
0
            pszShared = CPLGetConfigOption("VRT_SHARED_SOURCE", nullptr);
614
0
        }
615
0
        if (pszShared != nullptr)
616
0
        {
617
0
            m_nExplicitSharedStatus = CPLTestBool(pszShared);
618
0
        }
619
620
0
        m_osSrcDSName = GDALDataset::BuildFilename(
621
0
            pszFilename, pszVRTPath, CPL_TO_BOOL(m_bRelativeToVRTOri));
622
0
    }
623
0
    else if (psSourceVRTDataset)
624
0
    {
625
0
        CPLXMLNode sNode;
626
0
        sNode.eType = psSourceVRTDataset->eType;
627
0
        sNode.pszValue = psSourceVRTDataset->pszValue;
628
0
        sNode.psNext = nullptr;
629
0
        sNode.psChild = psSourceVRTDataset->psChild;
630
0
        char *pszXML = CPLSerializeXMLTree(&sNode);
631
0
        if (pszXML)
632
0
        {
633
0
            m_bSrcDSNameFromVRT = true;
634
0
            m_osSrcDSName = pszXML;
635
0
            CPLFree(pszXML);
636
0
        }
637
0
    }
638
639
0
    const char *pszSourceBand = CPLGetXMLValue(psSrc, "SourceBand", "1");
640
0
    m_bGetMaskBand = false;
641
0
    if (STARTS_WITH_CI(pszSourceBand, "mask"))
642
0
    {
643
0
        m_bGetMaskBand = true;
644
0
        if (pszSourceBand[4] == ',')
645
0
            m_nBand = atoi(pszSourceBand + 5);
646
0
        else
647
0
            m_nBand = 1;
648
0
    }
649
0
    else
650
0
    {
651
0
        m_nBand = atoi(pszSourceBand);
652
0
    }
653
0
    if (!GDALCheckBandCount(m_nBand, 0))
654
0
    {
655
0
        CPLError(CE_Warning, CPLE_AppDefined,
656
0
                 "Invalid <SourceBand> element in VRTRasterBand.");
657
0
        return CE_Failure;
658
0
    }
659
660
0
    m_aosOpenOptions = GDALDeserializeOpenOptionsFromXML(psSrc);
661
0
    m_aosOpenOptionsOri = m_aosOpenOptions;
662
0
    if (strstr(m_osSrcDSName.c_str(), "<VRTDataset") != nullptr)
663
0
        m_aosOpenOptions.SetNameValue("ROOT_PATH", pszVRTPath);
664
665
0
    return ParseSrcRectAndDstRect(psSrc);
666
0
}
667
668
/************************************************************************/
669
/*                        ParseSrcRectAndDstRect()                      */
670
/************************************************************************/
671
672
CPLErr VRTSimpleSource::ParseSrcRectAndDstRect(const CPLXMLNode *psSrc)
673
0
{
674
0
    const auto GetAttrValue = [](const CPLXMLNode *psNode,
675
0
                                 const char *pszAttrName, double dfDefaultVal)
676
0
    {
677
0
        if (const char *pszVal = CPLGetXMLValue(psNode, pszAttrName, nullptr))
678
0
            return CPLAtof(pszVal);
679
0
        else
680
0
            return dfDefaultVal;
681
0
    };
682
683
    /* -------------------------------------------------------------------- */
684
    /*      Set characteristics.                                            */
685
    /* -------------------------------------------------------------------- */
686
0
    const CPLXMLNode *const psSrcRect = CPLGetXMLNode(psSrc, "SrcRect");
687
0
    if (psSrcRect)
688
0
    {
689
0
        double xOff = GetAttrValue(psSrcRect, "xOff", UNINIT_WINDOW);
690
0
        double yOff = GetAttrValue(psSrcRect, "yOff", UNINIT_WINDOW);
691
0
        double xSize = GetAttrValue(psSrcRect, "xSize", UNINIT_WINDOW);
692
0
        double ySize = GetAttrValue(psSrcRect, "ySize", UNINIT_WINDOW);
693
        // Test written that way to catch NaN values
694
0
        if (!(xOff >= INT_MIN && xOff <= INT_MAX) ||
695
0
            !(yOff >= INT_MIN && yOff <= INT_MAX) ||
696
0
            !(xSize > 0 || xSize == UNINIT_WINDOW) || xSize > INT_MAX ||
697
0
            !(ySize > 0 || ySize == UNINIT_WINDOW) || ySize > INT_MAX)
698
0
        {
699
0
            CPLError(CE_Failure, CPLE_AppDefined, "Wrong values in SrcRect");
700
0
            return CE_Failure;
701
0
        }
702
0
        SetSrcWindow(xOff, yOff, xSize, ySize);
703
0
    }
704
0
    else
705
0
    {
706
0
        m_dfSrcXOff = UNINIT_WINDOW;
707
0
        m_dfSrcYOff = UNINIT_WINDOW;
708
0
        m_dfSrcXSize = UNINIT_WINDOW;
709
0
        m_dfSrcYSize = UNINIT_WINDOW;
710
0
    }
711
712
0
    const CPLXMLNode *const psDstRect = CPLGetXMLNode(psSrc, "DstRect");
713
0
    if (psDstRect)
714
0
    {
715
0
        double xOff = GetAttrValue(psDstRect, "xOff", UNINIT_WINDOW);
716
0
        ;
717
0
        double yOff = GetAttrValue(psDstRect, "yOff", UNINIT_WINDOW);
718
0
        double xSize = GetAttrValue(psDstRect, "xSize", UNINIT_WINDOW);
719
0
        ;
720
0
        double ySize = GetAttrValue(psDstRect, "ySize", UNINIT_WINDOW);
721
        // Test written that way to catch NaN values
722
0
        if (!(xOff >= INT_MIN && xOff <= INT_MAX) ||
723
0
            !(yOff >= INT_MIN && yOff <= INT_MAX) ||
724
0
            !(xSize > 0 || xSize == UNINIT_WINDOW) || xSize > INT_MAX ||
725
0
            !(ySize > 0 || ySize == UNINIT_WINDOW) || ySize > INT_MAX)
726
0
        {
727
0
            CPLError(CE_Failure, CPLE_AppDefined, "Wrong values in DstRect");
728
0
            return CE_Failure;
729
0
        }
730
0
        SetDstWindow(xOff, yOff, xSize, ySize);
731
0
    }
732
0
    else
733
0
    {
734
0
        m_dfDstXOff = UNINIT_WINDOW;
735
0
        m_dfDstYOff = UNINIT_WINDOW;
736
0
        m_dfDstXSize = UNINIT_WINDOW;
737
0
        m_dfDstYSize = UNINIT_WINDOW;
738
0
    }
739
740
0
    return CE_None;
741
0
}
742
743
/************************************************************************/
744
/*                             GetFileList()                            */
745
/************************************************************************/
746
747
void VRTSimpleSource::GetFileList(char ***ppapszFileList, int *pnSize,
748
                                  int *pnMaxSize, CPLHashSet *hSetFiles)
749
0
{
750
0
    if (!m_osSrcDSName.empty() && !m_bSrcDSNameFromVRT)
751
0
    {
752
0
        const char *pszFilename = m_osSrcDSName.c_str();
753
754
        /* --------------------------------------------------------------------
755
         */
756
        /*      Is it already in the list ? */
757
        /* --------------------------------------------------------------------
758
         */
759
0
        if (CPLHashSetLookup(hSetFiles, pszFilename) != nullptr)
760
0
            return;
761
762
        /* --------------------------------------------------------------------
763
         */
764
        /*      Grow array if necessary */
765
        /* --------------------------------------------------------------------
766
         */
767
0
        if (*pnSize + 1 >= *pnMaxSize)
768
0
        {
769
0
            *pnMaxSize = std::max(*pnSize + 2, 2 + 2 * (*pnMaxSize));
770
0
            *ppapszFileList = static_cast<char **>(
771
0
                CPLRealloc(*ppapszFileList, sizeof(char *) * (*pnMaxSize)));
772
0
        }
773
774
        /* --------------------------------------------------------------------
775
         */
776
        /*      Add the string to the list */
777
        /* --------------------------------------------------------------------
778
         */
779
0
        (*ppapszFileList)[*pnSize] = CPLStrdup(pszFilename);
780
0
        (*ppapszFileList)[(*pnSize + 1)] = nullptr;
781
0
        CPLHashSetInsert(hSetFiles, (*ppapszFileList)[*pnSize]);
782
783
0
        (*pnSize)++;
784
0
    }
785
0
}
786
787
/************************************************************************/
788
/*                           OpenSource()                               */
789
/************************************************************************/
790
791
void VRTSimpleSource::OpenSource() const
792
0
{
793
0
    CPLAssert(m_poRasterBand == nullptr);
794
795
    /* ----------------------------------------------------------------- */
796
    /*      Create a proxy dataset                                       */
797
    /* ----------------------------------------------------------------- */
798
0
    GDALProxyPoolDataset *proxyDS = nullptr;
799
0
    std::string osKeyMapSharedSources;
800
0
    if (m_poMapSharedSources)
801
0
    {
802
0
        osKeyMapSharedSources = m_osSrcDSName;
803
0
        for (int i = 0; i < m_aosOpenOptions.size(); ++i)
804
0
        {
805
0
            osKeyMapSharedSources += "||";
806
0
            osKeyMapSharedSources += m_aosOpenOptions[i];
807
0
        }
808
809
0
        proxyDS = cpl::down_cast<GDALProxyPoolDataset *>(
810
0
            m_poMapSharedSources->Get(osKeyMapSharedSources));
811
0
    }
812
813
0
    if (proxyDS == nullptr)
814
0
    {
815
0
        int bShared = true;
816
0
        if (m_nExplicitSharedStatus != -1)
817
0
            bShared = m_nExplicitSharedStatus;
818
819
0
        const CPLString osUniqueHandle(CPLSPrintf("%p", m_poMapSharedSources));
820
0
        proxyDS = GDALProxyPoolDataset::Create(
821
0
            m_osSrcDSName, m_aosOpenOptions.List(), GA_ReadOnly, bShared,
822
0
            osUniqueHandle.c_str());
823
0
        if (proxyDS == nullptr)
824
0
            return;
825
0
    }
826
0
    else
827
0
    {
828
0
        proxyDS->Reference();
829
0
    }
830
831
0
    if (m_bGetMaskBand)
832
0
    {
833
0
        GDALProxyPoolRasterBand *poMaskBand =
834
0
            cpl::down_cast<GDALProxyPoolRasterBand *>(
835
0
                proxyDS->GetRasterBand(m_nBand));
836
0
        poMaskBand->AddSrcMaskBandDescriptionFromUnderlying();
837
0
    }
838
839
    /* -------------------------------------------------------------------- */
840
    /*      Get the raster band.                                            */
841
    /* -------------------------------------------------------------------- */
842
843
0
    m_poRasterBand = proxyDS->GetRasterBand(m_nBand);
844
0
    if (m_poRasterBand == nullptr || !ValidateOpenedBand(m_poRasterBand))
845
0
    {
846
0
        proxyDS->ReleaseRef();
847
0
        return;
848
0
    }
849
850
0
    if (m_bGetMaskBand)
851
0
    {
852
0
        m_poRasterBand = m_poRasterBand->GetMaskBand();
853
0
        if (m_poRasterBand == nullptr)
854
0
        {
855
0
            proxyDS->ReleaseRef();
856
0
            return;
857
0
        }
858
0
        m_poMaskBandMainBand = m_poRasterBand;
859
0
    }
860
861
0
    if (m_poMapSharedSources)
862
0
    {
863
0
        m_poMapSharedSources->Insert(osKeyMapSharedSources, proxyDS);
864
0
    }
865
0
}
866
867
/************************************************************************/
868
/*                         GetRasterBand()                              */
869
/************************************************************************/
870
871
GDALRasterBand *VRTSimpleSource::GetRasterBand() const
872
0
{
873
0
    if (m_poRasterBand == nullptr)
874
0
        OpenSource();
875
0
    return m_poRasterBand;
876
0
}
877
878
/************************************************************************/
879
/*                        GetMaskBandMainBand()                         */
880
/************************************************************************/
881
882
GDALRasterBand *VRTSimpleSource::GetMaskBandMainBand()
883
0
{
884
0
    if (m_poRasterBand == nullptr)
885
0
        OpenSource();
886
0
    return m_poMaskBandMainBand;
887
0
}
888
889
/************************************************************************/
890
/*                       IsSameExceptBandNumber()                       */
891
/************************************************************************/
892
893
bool VRTSimpleSource::IsSameExceptBandNumber(
894
    const VRTSimpleSource *poOtherSource) const
895
0
{
896
0
    return m_dfSrcXOff == poOtherSource->m_dfSrcXOff &&
897
0
           m_dfSrcYOff == poOtherSource->m_dfSrcYOff &&
898
0
           m_dfSrcXSize == poOtherSource->m_dfSrcXSize &&
899
0
           m_dfSrcYSize == poOtherSource->m_dfSrcYSize &&
900
0
           m_dfDstXOff == poOtherSource->m_dfDstXOff &&
901
0
           m_dfDstYOff == poOtherSource->m_dfDstYOff &&
902
0
           m_dfDstXSize == poOtherSource->m_dfDstXSize &&
903
0
           m_dfDstYSize == poOtherSource->m_dfDstYSize &&
904
0
           m_osSrcDSName == poOtherSource->m_osSrcDSName;
905
0
}
906
907
/************************************************************************/
908
/*                              SrcToDst()                              */
909
/*                                                                      */
910
/*      Note: this is a no-op if the both src and dst windows are unset */
911
/************************************************************************/
912
913
void VRTSimpleSource::SrcToDst(double dfX, double dfY, double &dfXOut,
914
                               double &dfYOut) const
915
916
0
{
917
0
    dfXOut = ((dfX - m_dfSrcXOff) / m_dfSrcXSize) * m_dfDstXSize + m_dfDstXOff;
918
0
    dfYOut = ((dfY - m_dfSrcYOff) / m_dfSrcYSize) * m_dfDstYSize + m_dfDstYOff;
919
0
}
920
921
/************************************************************************/
922
/*                              DstToSrc()                              */
923
/*                                                                      */
924
/*      Note: this is a no-op if the both src and dst windows are unset */
925
/************************************************************************/
926
927
void VRTSimpleSource::DstToSrc(double dfX, double dfY, double &dfXOut,
928
                               double &dfYOut) const
929
930
0
{
931
0
    dfXOut = ((dfX - m_dfDstXOff) / m_dfDstXSize) * m_dfSrcXSize + m_dfSrcXOff;
932
0
    dfYOut = ((dfY - m_dfDstYOff) / m_dfDstYSize) * m_dfSrcYSize + m_dfSrcYOff;
933
0
}
934
935
/************************************************************************/
936
/*                          GetSrcDstWindow()                           */
937
/************************************************************************/
938
939
int VRTSimpleSource::GetSrcDstWindow(
940
    double dfXOff, double dfYOff, double dfXSize, double dfYSize, int nBufXSize,
941
    int nBufYSize, double *pdfReqXOff, double *pdfReqYOff, double *pdfReqXSize,
942
    double *pdfReqYSize, int *pnReqXOff, int *pnReqYOff, int *pnReqXSize,
943
    int *pnReqYSize, int *pnOutXOff, int *pnOutYOff, int *pnOutXSize,
944
    int *pnOutYSize, bool &bErrorOut)
945
946
0
{
947
0
    bErrorOut = false;
948
949
0
    if (m_dfSrcXSize == 0.0 || m_dfSrcYSize == 0.0 || m_dfDstXSize == 0.0 ||
950
0
        m_dfDstYSize == 0.0)
951
0
    {
952
0
        return FALSE;
953
0
    }
954
955
0
    const bool bDstWinSet = IsDstWinSet();
956
957
0
#ifdef DEBUG
958
0
    const bool bSrcWinSet = IsSrcWinSet();
959
960
0
    if (bSrcWinSet != bDstWinSet)
961
0
    {
962
0
        return FALSE;
963
0
    }
964
0
#endif
965
966
    /* -------------------------------------------------------------------- */
967
    /*      If the input window completely misses the portion of the        */
968
    /*      virtual dataset provided by this source we have nothing to do.  */
969
    /* -------------------------------------------------------------------- */
970
0
    if (bDstWinSet)
971
0
    {
972
0
        if (dfXOff >= m_dfDstXOff + m_dfDstXSize ||
973
0
            dfYOff >= m_dfDstYOff + m_dfDstYSize ||
974
0
            dfXOff + dfXSize <= m_dfDstXOff || dfYOff + dfYSize <= m_dfDstYOff)
975
0
            return FALSE;
976
0
    }
977
978
    /* -------------------------------------------------------------------- */
979
    /*      This request window corresponds to the whole output buffer.     */
980
    /* -------------------------------------------------------------------- */
981
0
    *pnOutXOff = 0;
982
0
    *pnOutYOff = 0;
983
0
    *pnOutXSize = nBufXSize;
984
0
    *pnOutYSize = nBufYSize;
985
986
    /* -------------------------------------------------------------------- */
987
    /*      If the input window extents outside the portion of the on       */
988
    /*      the virtual file that this source can set, then clip down       */
989
    /*      the requested window.                                           */
990
    /* -------------------------------------------------------------------- */
991
0
    bool bModifiedX = false;
992
0
    bool bModifiedY = false;
993
0
    double dfRXOff = dfXOff;
994
0
    double dfRYOff = dfYOff;
995
0
    double dfRXSize = dfXSize;
996
0
    double dfRYSize = dfYSize;
997
998
0
    if (bDstWinSet)
999
0
    {
1000
0
        if (dfRXOff < m_dfDstXOff)
1001
0
        {
1002
0
            dfRXSize = dfRXSize + dfRXOff - m_dfDstXOff;
1003
0
            dfRXOff = m_dfDstXOff;
1004
0
            bModifiedX = true;
1005
0
        }
1006
1007
0
        if (dfRYOff < m_dfDstYOff)
1008
0
        {
1009
0
            dfRYSize = dfRYSize + dfRYOff - m_dfDstYOff;
1010
0
            dfRYOff = m_dfDstYOff;
1011
0
            bModifiedY = true;
1012
0
        }
1013
1014
0
        if (dfRXOff + dfRXSize > m_dfDstXOff + m_dfDstXSize)
1015
0
        {
1016
0
            dfRXSize = m_dfDstXOff + m_dfDstXSize - dfRXOff;
1017
0
            bModifiedX = true;
1018
0
        }
1019
1020
0
        if (dfRYOff + dfRYSize > m_dfDstYOff + m_dfDstYSize)
1021
0
        {
1022
0
            dfRYSize = m_dfDstYOff + m_dfDstYSize - dfRYOff;
1023
0
            bModifiedY = true;
1024
0
        }
1025
0
    }
1026
1027
    /* -------------------------------------------------------------------- */
1028
    /*      Translate requested region in virtual file into the source      */
1029
    /*      band coordinates.                                               */
1030
    /* -------------------------------------------------------------------- */
1031
0
    const double dfScaleX = m_dfSrcXSize / m_dfDstXSize;
1032
0
    const double dfScaleY = m_dfSrcYSize / m_dfDstYSize;
1033
1034
0
    *pdfReqXOff = (dfRXOff - m_dfDstXOff) * dfScaleX + m_dfSrcXOff;
1035
0
    *pdfReqYOff = (dfRYOff - m_dfDstYOff) * dfScaleY + m_dfSrcYOff;
1036
0
    *pdfReqXSize = dfRXSize * dfScaleX;
1037
0
    *pdfReqYSize = dfRYSize * dfScaleY;
1038
1039
0
    if (!std::isfinite(*pdfReqXOff) || !std::isfinite(*pdfReqYOff) ||
1040
0
        !std::isfinite(*pdfReqXSize) || !std::isfinite(*pdfReqYSize) ||
1041
0
        *pdfReqXOff > INT_MAX || *pdfReqYOff > INT_MAX || *pdfReqXSize < 0 ||
1042
0
        *pdfReqYSize < 0)
1043
0
    {
1044
0
        return FALSE;
1045
0
    }
1046
1047
    /* -------------------------------------------------------------------- */
1048
    /*      Clamp within the bounds of the available source data.           */
1049
    /* -------------------------------------------------------------------- */
1050
0
    if (*pdfReqXOff < 0)
1051
0
    {
1052
0
        *pdfReqXSize += *pdfReqXOff;
1053
0
        *pdfReqXOff = 0;
1054
0
        bModifiedX = true;
1055
0
    }
1056
0
    if (*pdfReqYOff < 0)
1057
0
    {
1058
0
        *pdfReqYSize += *pdfReqYOff;
1059
0
        *pdfReqYOff = 0;
1060
0
        bModifiedY = true;
1061
0
    }
1062
1063
0
    *pnReqXOff = static_cast<int>(floor(*pdfReqXOff));
1064
0
    *pnReqYOff = static_cast<int>(floor(*pdfReqYOff));
1065
1066
0
    constexpr double EPS = 1e-3;
1067
0
    constexpr double ONE_MINUS_EPS = 1.0 - EPS;
1068
0
    if (*pdfReqXOff - *pnReqXOff > ONE_MINUS_EPS)
1069
0
    {
1070
0
        (*pnReqXOff)++;
1071
0
        *pdfReqXOff = *pnReqXOff;
1072
0
    }
1073
0
    if (*pdfReqYOff - *pnReqYOff > ONE_MINUS_EPS)
1074
0
    {
1075
0
        (*pnReqYOff)++;
1076
0
        *pdfReqYOff = *pnReqYOff;
1077
0
    }
1078
1079
0
    if (*pdfReqXSize > INT_MAX)
1080
0
        *pnReqXSize = INT_MAX;
1081
0
    else
1082
0
        *pnReqXSize = static_cast<int>(floor(*pdfReqXSize + 0.5));
1083
1084
0
    if (*pdfReqYSize > INT_MAX)
1085
0
        *pnReqYSize = INT_MAX;
1086
0
    else
1087
0
        *pnReqYSize = static_cast<int>(floor(*pdfReqYSize + 0.5));
1088
1089
    /* -------------------------------------------------------------------- */
1090
    /*      Clamp within the bounds of the available source data.           */
1091
    /* -------------------------------------------------------------------- */
1092
1093
0
    if (*pnReqXSize == 0)
1094
0
        *pnReqXSize = 1;
1095
0
    if (*pnReqYSize == 0)
1096
0
        *pnReqYSize = 1;
1097
1098
0
    auto l_band = GetRasterBand();
1099
0
    if (!l_band)
1100
0
    {
1101
0
        bErrorOut = true;
1102
0
        return FALSE;
1103
0
    }
1104
0
    if (*pnReqXSize > INT_MAX - *pnReqXOff ||
1105
0
        *pnReqXOff + *pnReqXSize > l_band->GetXSize())
1106
0
    {
1107
0
        *pnReqXSize = l_band->GetXSize() - *pnReqXOff;
1108
0
        bModifiedX = true;
1109
0
    }
1110
0
    if (*pdfReqXOff + *pdfReqXSize > l_band->GetXSize())
1111
0
    {
1112
0
        *pdfReqXSize = l_band->GetXSize() - *pdfReqXOff;
1113
0
        bModifiedX = true;
1114
0
    }
1115
1116
0
    if (*pnReqYSize > INT_MAX - *pnReqYOff ||
1117
0
        *pnReqYOff + *pnReqYSize > l_band->GetYSize())
1118
0
    {
1119
0
        *pnReqYSize = l_band->GetYSize() - *pnReqYOff;
1120
0
        bModifiedY = true;
1121
0
    }
1122
0
    if (*pdfReqYOff + *pdfReqYSize > l_band->GetYSize())
1123
0
    {
1124
0
        *pdfReqYSize = l_band->GetYSize() - *pdfReqYOff;
1125
0
        bModifiedY = true;
1126
0
    }
1127
1128
    /* -------------------------------------------------------------------- */
1129
    /*      Don't do anything if the requesting region is completely off    */
1130
    /*      the source image.                                               */
1131
    /* -------------------------------------------------------------------- */
1132
0
    if (*pnReqXOff >= l_band->GetXSize() || *pnReqYOff >= l_band->GetYSize() ||
1133
0
        *pnReqXSize <= 0 || *pnReqYSize <= 0)
1134
0
    {
1135
0
        return FALSE;
1136
0
    }
1137
1138
    /* -------------------------------------------------------------------- */
1139
    /*      If we haven't had to modify the source rectangle, then the      */
1140
    /*      destination rectangle must be the whole region.                 */
1141
    /* -------------------------------------------------------------------- */
1142
0
    if (bModifiedX || bModifiedY)
1143
0
    {
1144
        /* --------------------------------------------------------------------
1145
         */
1146
        /*      Now transform this possibly reduced request back into the */
1147
        /*      destination buffer coordinates in case the output region is */
1148
        /*      less than the whole buffer. */
1149
        /* --------------------------------------------------------------------
1150
         */
1151
0
        double dfDstULX = 0.0;
1152
0
        double dfDstULY = 0.0;
1153
0
        double dfDstLRX = 0.0;
1154
0
        double dfDstLRY = 0.0;
1155
1156
0
        SrcToDst(*pdfReqXOff, *pdfReqYOff, dfDstULX, dfDstULY);
1157
0
        SrcToDst(*pdfReqXOff + *pdfReqXSize, *pdfReqYOff + *pdfReqYSize,
1158
0
                 dfDstLRX, dfDstLRY);
1159
#if DEBUG_VERBOSE
1160
        CPLDebug("VRT", "dfDstULX=%g dfDstULY=%g dfDstLRX=%g dfDstLRY=%g",
1161
                 dfDstULX, dfDstULY, dfDstLRX, dfDstLRY);
1162
#endif
1163
1164
0
        if (bModifiedX)
1165
0
        {
1166
0
            const double dfScaleWinToBufX = nBufXSize / dfXSize;
1167
1168
0
            const double dfOutXOff = (dfDstULX - dfXOff) * dfScaleWinToBufX;
1169
0
            if (dfOutXOff <= 0)
1170
0
                *pnOutXOff = 0;
1171
0
            else if (dfOutXOff > INT_MAX)
1172
0
                *pnOutXOff = INT_MAX;
1173
0
            else
1174
0
                *pnOutXOff = static_cast<int>(dfOutXOff + EPS);
1175
1176
            // Apply correction on floating-point source window
1177
0
            {
1178
0
                double dfDstDeltaX =
1179
0
                    (dfOutXOff - *pnOutXOff) / dfScaleWinToBufX;
1180
0
                double dfSrcDeltaX = dfDstDeltaX / m_dfDstXSize * m_dfSrcXSize;
1181
0
                *pdfReqXOff -= dfSrcDeltaX;
1182
0
                *pdfReqXSize = std::min(*pdfReqXSize + dfSrcDeltaX,
1183
0
                                        static_cast<double>(INT_MAX));
1184
0
            }
1185
1186
0
            double dfOutRightXOff = (dfDstLRX - dfXOff) * dfScaleWinToBufX;
1187
0
            if (dfOutRightXOff < dfOutXOff)
1188
0
                return FALSE;
1189
0
            if (dfOutRightXOff > INT_MAX)
1190
0
                dfOutRightXOff = INT_MAX;
1191
0
            const int nOutRightXOff =
1192
0
                static_cast<int>(ceil(dfOutRightXOff - EPS));
1193
0
            *pnOutXSize = nOutRightXOff - *pnOutXOff;
1194
1195
0
            if (*pnOutXSize > INT_MAX - *pnOutXOff ||
1196
0
                *pnOutXOff + *pnOutXSize > nBufXSize)
1197
0
                *pnOutXSize = nBufXSize - *pnOutXOff;
1198
1199
            // Apply correction on floating-point source window
1200
0
            {
1201
0
                double dfDstDeltaX =
1202
0
                    (nOutRightXOff - dfOutRightXOff) / dfScaleWinToBufX;
1203
0
                double dfSrcDeltaX = dfDstDeltaX / m_dfDstXSize * m_dfSrcXSize;
1204
0
                *pdfReqXSize = std::min(*pdfReqXSize + dfSrcDeltaX,
1205
0
                                        static_cast<double>(INT_MAX));
1206
0
            }
1207
0
        }
1208
1209
0
        if (bModifiedY)
1210
0
        {
1211
0
            const double dfScaleWinToBufY = nBufYSize / dfYSize;
1212
1213
0
            const double dfOutYOff = (dfDstULY - dfYOff) * dfScaleWinToBufY;
1214
0
            if (dfOutYOff <= 0)
1215
0
                *pnOutYOff = 0;
1216
0
            else if (dfOutYOff > INT_MAX)
1217
0
                *pnOutYOff = INT_MAX;
1218
0
            else
1219
0
                *pnOutYOff = static_cast<int>(dfOutYOff + EPS);
1220
1221
            // Apply correction on floating-point source window
1222
0
            {
1223
0
                double dfDstDeltaY =
1224
0
                    (dfOutYOff - *pnOutYOff) / dfScaleWinToBufY;
1225
0
                double dfSrcDeltaY = dfDstDeltaY / m_dfDstYSize * m_dfSrcYSize;
1226
0
                *pdfReqYOff -= dfSrcDeltaY;
1227
0
                *pdfReqYSize = std::min(*pdfReqYSize + dfSrcDeltaY,
1228
0
                                        static_cast<double>(INT_MAX));
1229
0
            }
1230
1231
0
            double dfOutTopYOff = (dfDstLRY - dfYOff) * dfScaleWinToBufY;
1232
0
            if (dfOutTopYOff < dfOutYOff)
1233
0
                return FALSE;
1234
0
            if (dfOutTopYOff > INT_MAX)
1235
0
                dfOutTopYOff = INT_MAX;
1236
0
            const int nOutTopYOff = static_cast<int>(ceil(dfOutTopYOff - EPS));
1237
0
            *pnOutYSize = nOutTopYOff - *pnOutYOff;
1238
1239
0
            if (*pnOutYSize > INT_MAX - *pnOutYOff ||
1240
0
                *pnOutYOff + *pnOutYSize > nBufYSize)
1241
0
                *pnOutYSize = nBufYSize - *pnOutYOff;
1242
1243
            // Apply correction on floating-point source window
1244
0
            {
1245
0
                double dfDstDeltaY =
1246
0
                    (nOutTopYOff - dfOutTopYOff) / dfScaleWinToBufY;
1247
0
                double dfSrcDeltaY = dfDstDeltaY / m_dfDstYSize * m_dfSrcYSize;
1248
0
                *pdfReqYSize = std::min(*pdfReqYSize + dfSrcDeltaY,
1249
0
                                        static_cast<double>(INT_MAX));
1250
0
            }
1251
0
        }
1252
1253
0
        if (*pnOutXSize < 1 || *pnOutYSize < 1)
1254
0
            return FALSE;
1255
0
    }
1256
1257
0
    *pdfReqXOff = RoundIfCloseToInt(*pdfReqXOff);
1258
0
    *pdfReqYOff = RoundIfCloseToInt(*pdfReqYOff);
1259
0
    *pdfReqXSize = RoundIfCloseToInt(*pdfReqXSize);
1260
0
    *pdfReqYSize = RoundIfCloseToInt(*pdfReqYSize);
1261
1262
0
    return TRUE;
1263
0
}
1264
1265
/************************************************************************/
1266
/*                          NeedMaxValAdjustment()                      */
1267
/************************************************************************/
1268
1269
int VRTSimpleSource::NeedMaxValAdjustment() const
1270
0
{
1271
0
    if (!m_nMaxValue)
1272
0
        return FALSE;
1273
1274
0
    auto l_band = GetRasterBand();
1275
0
    if (!l_band)
1276
0
        return FALSE;
1277
0
    const char *pszNBITS = l_band->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1278
0
    const int nBits = (pszNBITS) ? atoi(pszNBITS) : 0;
1279
0
    if (nBits >= 1 && nBits <= 31)
1280
0
    {
1281
0
        const int nBandMaxValue = static_cast<int>((1U << nBits) - 1);
1282
0
        return nBandMaxValue > m_nMaxValue;
1283
0
    }
1284
0
    return TRUE;
1285
0
}
1286
1287
/************************************************************************/
1288
/*                              RasterIO()                              */
1289
/************************************************************************/
1290
1291
CPLErr VRTSimpleSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
1292
                                 int nYOff, int nXSize, int nYSize, void *pData,
1293
                                 int nBufXSize, int nBufYSize,
1294
                                 GDALDataType eBufType, GSpacing nPixelSpace,
1295
                                 GSpacing nLineSpace,
1296
                                 GDALRasterIOExtraArg *psExtraArgIn,
1297
                                 WorkingState & /*oWorkingState*/)
1298
1299
0
{
1300
0
    GDALRasterIOExtraArg sExtraArg;
1301
0
    INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1302
0
    GDALRasterIOExtraArg *psExtraArg = &sExtraArg;
1303
1304
0
    double dfXOff = nXOff;
1305
0
    double dfYOff = nYOff;
1306
0
    double dfXSize = nXSize;
1307
0
    double dfYSize = nYSize;
1308
0
    if (psExtraArgIn != nullptr && psExtraArgIn->bFloatingPointWindowValidity)
1309
0
    {
1310
0
        dfXOff = psExtraArgIn->dfXOff;
1311
0
        dfYOff = psExtraArgIn->dfYOff;
1312
0
        dfXSize = psExtraArgIn->dfXSize;
1313
0
        dfYSize = psExtraArgIn->dfYSize;
1314
0
    }
1315
1316
    // The window we will actually request from the source raster band.
1317
0
    double dfReqXOff = 0.0;
1318
0
    double dfReqYOff = 0.0;
1319
0
    double dfReqXSize = 0.0;
1320
0
    double dfReqYSize = 0.0;
1321
0
    int nReqXOff = 0;
1322
0
    int nReqYOff = 0;
1323
0
    int nReqXSize = 0;
1324
0
    int nReqYSize = 0;
1325
1326
    // The window we will actual set _within_ the pData buffer.
1327
0
    int nOutXOff = 0;
1328
0
    int nOutYOff = 0;
1329
0
    int nOutXSize = 0;
1330
0
    int nOutYSize = 0;
1331
1332
0
    bool bError = false;
1333
0
    if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
1334
0
                         &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
1335
0
                         &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
1336
0
                         &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
1337
0
    {
1338
0
        return bError ? CE_Failure : CE_None;
1339
0
    }
1340
#if DEBUG_VERBOSE
1341
    CPLDebug("VRT",
1342
             "nXOff=%d, nYOff=%d, nXSize=%d, nYSize=%d, nBufXSize=%d, "
1343
             "nBufYSize=%d,\n"
1344
             "dfReqXOff=%g, dfReqYOff=%g, dfReqXSize=%g, dfReqYSize=%g,\n"
1345
             "nReqXOff=%d, nReqYOff=%d, nReqXSize=%d, nReqYSize=%d,\n"
1346
             "nOutXOff=%d, nOutYOff=%d, nOutXSize=%d, nOutYSize=%d",
1347
             nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize, dfReqXOff,
1348
             dfReqYOff, dfReqXSize, dfReqYSize, nReqXOff, nReqYOff, nReqXSize,
1349
             nReqYSize, nOutXOff, nOutYOff, nOutXSize, nOutYSize);
1350
#endif
1351
1352
    /* -------------------------------------------------------------------- */
1353
    /*      Actually perform the IO request.                                */
1354
    /* -------------------------------------------------------------------- */
1355
0
    if (!m_osResampling.empty())
1356
0
    {
1357
0
        psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
1358
0
    }
1359
0
    else if (psExtraArgIn != nullptr)
1360
0
    {
1361
0
        psExtraArg->eResampleAlg = psExtraArgIn->eResampleAlg;
1362
0
    }
1363
0
    psExtraArg->bFloatingPointWindowValidity = TRUE;
1364
0
    psExtraArg->dfXOff = dfReqXOff;
1365
0
    psExtraArg->dfYOff = dfReqYOff;
1366
0
    psExtraArg->dfXSize = dfReqXSize;
1367
0
    psExtraArg->dfYSize = dfReqYSize;
1368
0
    if (psExtraArgIn)
1369
0
    {
1370
0
        psExtraArg->pfnProgress = psExtraArgIn->pfnProgress;
1371
0
        psExtraArg->pProgressData = psExtraArgIn->pProgressData;
1372
0
        if (psExtraArgIn->nVersion >= 2)
1373
0
        {
1374
0
            psExtraArg->bUseOnlyThisScale = psExtraArgIn->bUseOnlyThisScale;
1375
0
        }
1376
0
    }
1377
1378
0
    GByte *pabyOut = static_cast<unsigned char *>(pData) +
1379
0
                     nOutXOff * nPixelSpace +
1380
0
                     static_cast<GPtrDiff_t>(nOutYOff) * nLineSpace;
1381
1382
0
    auto l_band = GetRasterBand();
1383
0
    if (!l_band)
1384
0
        return CE_Failure;
1385
1386
0
    CPLErr eErr = CE_Failure;
1387
0
    if (GDALDataTypeIsConversionLossy(l_band->GetRasterDataType(),
1388
0
                                      eVRTBandDataType))
1389
0
    {
1390
0
        const int nBandDTSize = GDALGetDataTypeSizeBytes(eVRTBandDataType);
1391
0
        void *pTemp = VSI_MALLOC3_VERBOSE(nOutXSize, nOutYSize, nBandDTSize);
1392
0
        if (pTemp)
1393
0
        {
1394
0
            eErr = l_band->RasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize,
1395
0
                                    nReqYSize, pTemp, nOutXSize, nOutYSize,
1396
0
                                    eVRTBandDataType, 0, 0, psExtraArg);
1397
0
            if (eErr == CE_None)
1398
0
            {
1399
0
                GByte *pabyTemp = static_cast<GByte *>(pTemp);
1400
0
                for (int iY = 0; iY < nOutYSize; iY++)
1401
0
                {
1402
0
                    GDALCopyWords(
1403
0
                        pabyTemp +
1404
0
                            static_cast<size_t>(iY) * nBandDTSize * nOutXSize,
1405
0
                        eVRTBandDataType, nBandDTSize,
1406
0
                        pabyOut + static_cast<GPtrDiff_t>(iY * nLineSpace),
1407
0
                        eBufType, static_cast<int>(nPixelSpace), nOutXSize);
1408
0
                }
1409
0
            }
1410
0
            VSIFree(pTemp);
1411
0
        }
1412
0
    }
1413
0
    else
1414
0
    {
1415
0
        eErr = l_band->RasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize,
1416
0
                                nReqYSize, pabyOut, nOutXSize, nOutYSize,
1417
0
                                eBufType, nPixelSpace, nLineSpace, psExtraArg);
1418
0
    }
1419
1420
0
    if (NeedMaxValAdjustment())
1421
0
    {
1422
0
        for (int j = 0; j < nOutYSize; j++)
1423
0
        {
1424
0
            for (int i = 0; i < nOutXSize; i++)
1425
0
            {
1426
0
                int nVal = 0;
1427
0
                GDALCopyWords(pabyOut + j * nLineSpace + i * nPixelSpace,
1428
0
                              eBufType, 0, &nVal, GDT_Int32, 0, 1);
1429
0
                if (nVal > m_nMaxValue)
1430
0
                    nVal = m_nMaxValue;
1431
0
                GDALCopyWords(&nVal, GDT_Int32, 0,
1432
0
                              pabyOut + j * nLineSpace + i * nPixelSpace,
1433
0
                              eBufType, 0, 1);
1434
0
            }
1435
0
        }
1436
0
    }
1437
1438
0
    if (psExtraArg->pfnProgress)
1439
0
        psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
1440
1441
0
    return eErr;
1442
0
}
1443
1444
/************************************************************************/
1445
/*                             GetMinimum()                             */
1446
/************************************************************************/
1447
1448
double VRTSimpleSource::GetMinimum(int nXSize, int nYSize, int *pbSuccess)
1449
0
{
1450
    // The window we will actually request from the source raster band.
1451
0
    double dfReqXOff = 0.0;
1452
0
    double dfReqYOff = 0.0;
1453
0
    double dfReqXSize = 0.0;
1454
0
    double dfReqYSize = 0.0;
1455
0
    int nReqXOff = 0;
1456
0
    int nReqYOff = 0;
1457
0
    int nReqXSize = 0;
1458
0
    int nReqYSize = 0;
1459
1460
    // The window we will actual set _within_ the pData buffer.
1461
0
    int nOutXOff = 0;
1462
0
    int nOutYOff = 0;
1463
0
    int nOutXSize = 0;
1464
0
    int nOutYSize = 0;
1465
1466
0
    bool bError = false;
1467
0
    auto l_band = GetRasterBand();
1468
0
    if (!l_band ||
1469
0
        !GetSrcDstWindow(0, 0, nXSize, nYSize, nXSize, nYSize, &dfReqXOff,
1470
0
                         &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
1471
0
                         &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff,
1472
0
                         &nOutYOff, &nOutXSize, &nOutYSize, bError) ||
1473
0
        nReqXOff != 0 || nReqYOff != 0 || nReqXSize != l_band->GetXSize() ||
1474
0
        nReqYSize != l_band->GetYSize())
1475
0
    {
1476
0
        *pbSuccess = FALSE;
1477
0
        return 0;
1478
0
    }
1479
1480
0
    const double dfVal = l_band->GetMinimum(pbSuccess);
1481
0
    if (NeedMaxValAdjustment() && dfVal > m_nMaxValue)
1482
0
        return m_nMaxValue;
1483
0
    return dfVal;
1484
0
}
1485
1486
/************************************************************************/
1487
/*                             GetMaximum()                             */
1488
/************************************************************************/
1489
1490
double VRTSimpleSource::GetMaximum(int nXSize, int nYSize, int *pbSuccess)
1491
0
{
1492
    // The window we will actually request from the source raster band.
1493
0
    double dfReqXOff = 0.0;
1494
0
    double dfReqYOff = 0.0;
1495
0
    double dfReqXSize = 0.0;
1496
0
    double dfReqYSize = 0.0;
1497
0
    int nReqXOff = 0;
1498
0
    int nReqYOff = 0;
1499
0
    int nReqXSize = 0;
1500
0
    int nReqYSize = 0;
1501
1502
    // The window we will actual set _within_ the pData buffer.
1503
0
    int nOutXOff = 0;
1504
0
    int nOutYOff = 0;
1505
0
    int nOutXSize = 0;
1506
0
    int nOutYSize = 0;
1507
1508
0
    bool bError = false;
1509
0
    auto l_band = GetRasterBand();
1510
0
    if (!l_band ||
1511
0
        !GetSrcDstWindow(0, 0, nXSize, nYSize, nXSize, nYSize, &dfReqXOff,
1512
0
                         &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
1513
0
                         &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff,
1514
0
                         &nOutYOff, &nOutXSize, &nOutYSize, bError) ||
1515
0
        nReqXOff != 0 || nReqYOff != 0 || nReqXSize != l_band->GetXSize() ||
1516
0
        nReqYSize != l_band->GetYSize())
1517
0
    {
1518
0
        *pbSuccess = FALSE;
1519
0
        return 0;
1520
0
    }
1521
1522
0
    const double dfVal = l_band->GetMaximum(pbSuccess);
1523
0
    if (NeedMaxValAdjustment() && dfVal > m_nMaxValue)
1524
0
        return m_nMaxValue;
1525
0
    return dfVal;
1526
0
}
1527
1528
/************************************************************************/
1529
/*                            GetHistogram()                            */
1530
/************************************************************************/
1531
1532
CPLErr VRTSimpleSource::GetHistogram(int nXSize, int nYSize, double dfMin,
1533
                                     double dfMax, int nBuckets,
1534
                                     GUIntBig *panHistogram,
1535
                                     int bIncludeOutOfRange, int bApproxOK,
1536
                                     GDALProgressFunc pfnProgress,
1537
                                     void *pProgressData)
1538
0
{
1539
    // The window we will actually request from the source raster band.
1540
0
    double dfReqXOff = 0.0;
1541
0
    double dfReqYOff = 0.0;
1542
0
    double dfReqXSize = 0.0;
1543
0
    double dfReqYSize = 0.0;
1544
0
    int nReqXOff = 0;
1545
0
    int nReqYOff = 0;
1546
0
    int nReqXSize = 0;
1547
0
    int nReqYSize = 0;
1548
1549
    // The window we will actual set _within_ the pData buffer.
1550
0
    int nOutXOff = 0;
1551
0
    int nOutYOff = 0;
1552
0
    int nOutXSize = 0;
1553
0
    int nOutYSize = 0;
1554
1555
0
    bool bError = false;
1556
0
    auto l_band = GetRasterBand();
1557
0
    if (!l_band || NeedMaxValAdjustment() ||
1558
0
        !GetSrcDstWindow(0, 0, nXSize, nYSize, nXSize, nYSize, &dfReqXOff,
1559
0
                         &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
1560
0
                         &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff,
1561
0
                         &nOutYOff, &nOutXSize, &nOutYSize, bError) ||
1562
0
        nReqXOff != 0 || nReqYOff != 0 || nReqXSize != l_band->GetXSize() ||
1563
0
        nReqYSize != l_band->GetYSize())
1564
0
    {
1565
0
        return CE_Failure;
1566
0
    }
1567
1568
0
    return l_band->GetHistogram(dfMin, dfMax, nBuckets, panHistogram,
1569
0
                                bIncludeOutOfRange, bApproxOK, pfnProgress,
1570
0
                                pProgressData);
1571
0
}
1572
1573
/************************************************************************/
1574
/*                          DatasetRasterIO()                           */
1575
/************************************************************************/
1576
1577
CPLErr VRTSimpleSource::DatasetRasterIO(
1578
    GDALDataType eVRTBandDataType, int nXOff, int nYOff, int nXSize, int nYSize,
1579
    void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
1580
    int nBandCount, const int *panBandMap, GSpacing nPixelSpace,
1581
    GSpacing nLineSpace, GSpacing nBandSpace,
1582
    GDALRasterIOExtraArg *psExtraArgIn)
1583
0
{
1584
0
    if (GetType() != VRTSimpleSource::GetTypeStatic())
1585
0
    {
1586
0
        CPLError(CE_Failure, CPLE_NotSupported,
1587
0
                 "DatasetRasterIO() not implemented for %s", GetType());
1588
0
        return CE_Failure;
1589
0
    }
1590
1591
0
    GDALRasterIOExtraArg sExtraArg;
1592
0
    INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1593
0
    GDALRasterIOExtraArg *psExtraArg = &sExtraArg;
1594
1595
0
    double dfXOff = nXOff;
1596
0
    double dfYOff = nYOff;
1597
0
    double dfXSize = nXSize;
1598
0
    double dfYSize = nYSize;
1599
0
    if (psExtraArgIn != nullptr && psExtraArgIn->bFloatingPointWindowValidity)
1600
0
    {
1601
0
        dfXOff = psExtraArgIn->dfXOff;
1602
0
        dfYOff = psExtraArgIn->dfYOff;
1603
0
        dfXSize = psExtraArgIn->dfXSize;
1604
0
        dfYSize = psExtraArgIn->dfYSize;
1605
0
    }
1606
1607
    // The window we will actually request from the source raster band.
1608
0
    double dfReqXOff = 0.0;
1609
0
    double dfReqYOff = 0.0;
1610
0
    double dfReqXSize = 0.0;
1611
0
    double dfReqYSize = 0.0;
1612
0
    int nReqXOff = 0;
1613
0
    int nReqYOff = 0;
1614
0
    int nReqXSize = 0;
1615
0
    int nReqYSize = 0;
1616
1617
    // The window we will actual set _within_ the pData buffer.
1618
0
    int nOutXOff = 0;
1619
0
    int nOutYOff = 0;
1620
0
    int nOutXSize = 0;
1621
0
    int nOutYSize = 0;
1622
1623
0
    bool bError = false;
1624
0
    if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
1625
0
                         &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
1626
0
                         &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
1627
0
                         &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
1628
0
    {
1629
0
        return bError ? CE_Failure : CE_None;
1630
0
    }
1631
1632
0
    auto l_band = GetRasterBand();
1633
0
    if (!l_band)
1634
0
        return CE_Failure;
1635
1636
0
    GDALDataset *poDS = l_band->GetDataset();
1637
0
    if (poDS == nullptr)
1638
0
        return CE_Failure;
1639
1640
0
    if (!m_osResampling.empty())
1641
0
    {
1642
0
        psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
1643
0
    }
1644
0
    else if (psExtraArgIn != nullptr)
1645
0
    {
1646
0
        psExtraArg->eResampleAlg = psExtraArgIn->eResampleAlg;
1647
0
    }
1648
0
    psExtraArg->bFloatingPointWindowValidity = TRUE;
1649
0
    psExtraArg->dfXOff = dfReqXOff;
1650
0
    psExtraArg->dfYOff = dfReqYOff;
1651
0
    psExtraArg->dfXSize = dfReqXSize;
1652
0
    psExtraArg->dfYSize = dfReqYSize;
1653
0
    if (psExtraArgIn)
1654
0
    {
1655
0
        psExtraArg->pfnProgress = psExtraArgIn->pfnProgress;
1656
0
        psExtraArg->pProgressData = psExtraArgIn->pProgressData;
1657
0
    }
1658
1659
0
    GByte *pabyOut = static_cast<unsigned char *>(pData) +
1660
0
                     nOutXOff * nPixelSpace +
1661
0
                     static_cast<GPtrDiff_t>(nOutYOff) * nLineSpace;
1662
1663
0
    CPLErr eErr = CE_Failure;
1664
1665
0
    if (GDALDataTypeIsConversionLossy(l_band->GetRasterDataType(),
1666
0
                                      eVRTBandDataType))
1667
0
    {
1668
0
        const int nBandDTSize = GDALGetDataTypeSizeBytes(eVRTBandDataType);
1669
0
        void *pTemp = VSI_MALLOC3_VERBOSE(
1670
0
            nOutXSize, nOutYSize, cpl::fits_on<int>(nBandDTSize * nBandCount));
1671
0
        if (pTemp)
1672
0
        {
1673
0
            eErr = poDS->RasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize,
1674
0
                                  nReqYSize, pTemp, nOutXSize, nOutYSize,
1675
0
                                  eVRTBandDataType, nBandCount, panBandMap, 0,
1676
0
                                  0, 0, psExtraArg);
1677
0
            if (eErr == CE_None)
1678
0
            {
1679
0
                GByte *pabyTemp = static_cast<GByte *>(pTemp);
1680
0
                const size_t nSrcBandSpace =
1681
0
                    static_cast<size_t>(nOutYSize) * nOutXSize * nBandDTSize;
1682
0
                for (int iBand = 0; iBand < nBandCount; iBand++)
1683
0
                {
1684
0
                    for (int iY = 0; iY < nOutYSize; iY++)
1685
0
                    {
1686
0
                        GDALCopyWords(
1687
0
                            pabyTemp + iBand * nSrcBandSpace +
1688
0
                                static_cast<size_t>(iY) * nBandDTSize *
1689
0
                                    nOutXSize,
1690
0
                            eVRTBandDataType, nBandDTSize,
1691
0
                            pabyOut + static_cast<GPtrDiff_t>(
1692
0
                                          iY * nLineSpace + iBand * nBandSpace),
1693
0
                            eBufType, static_cast<int>(nPixelSpace), nOutXSize);
1694
0
                    }
1695
0
                }
1696
0
            }
1697
0
            VSIFree(pTemp);
1698
0
        }
1699
0
    }
1700
0
    else
1701
0
    {
1702
0
        eErr = poDS->RasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
1703
0
                              pabyOut, nOutXSize, nOutYSize, eBufType,
1704
0
                              nBandCount, panBandMap, nPixelSpace, nLineSpace,
1705
0
                              nBandSpace, psExtraArg);
1706
0
    }
1707
1708
0
    if (NeedMaxValAdjustment())
1709
0
    {
1710
0
        for (int k = 0; k < nBandCount; k++)
1711
0
        {
1712
0
            for (int j = 0; j < nOutYSize; j++)
1713
0
            {
1714
0
                for (int i = 0; i < nOutXSize; i++)
1715
0
                {
1716
0
                    int nVal = 0;
1717
0
                    GDALCopyWords(pabyOut + k * nBandSpace + j * nLineSpace +
1718
0
                                      i * nPixelSpace,
1719
0
                                  eBufType, 0, &nVal, GDT_Int32, 0, 1);
1720
1721
0
                    if (nVal > m_nMaxValue)
1722
0
                        nVal = m_nMaxValue;
1723
1724
0
                    GDALCopyWords(&nVal, GDT_Int32, 0,
1725
0
                                  pabyOut + k * nBandSpace + j * nLineSpace +
1726
0
                                      i * nPixelSpace,
1727
0
                                  eBufType, 0, 1);
1728
0
                }
1729
0
            }
1730
0
        }
1731
0
    }
1732
1733
0
    if (psExtraArg->pfnProgress)
1734
0
        psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
1735
1736
0
    return eErr;
1737
0
}
1738
1739
/************************************************************************/
1740
/*                          SetResampling()                             */
1741
/************************************************************************/
1742
1743
void VRTSimpleSource::SetResampling(const char *pszResampling)
1744
0
{
1745
0
    m_osResampling = (pszResampling) ? pszResampling : "";
1746
0
}
1747
1748
/************************************************************************/
1749
/* ==================================================================== */
1750
/*                         VRTAveragedSource                            */
1751
/* ==================================================================== */
1752
/************************************************************************/
1753
1754
/************************************************************************/
1755
/*                         VRTAveragedSource()                          */
1756
/************************************************************************/
1757
1758
VRTAveragedSource::VRTAveragedSource()
1759
0
{
1760
0
}
1761
1762
/************************************************************************/
1763
/*                           GetTypeStatic()                            */
1764
/************************************************************************/
1765
1766
const char *VRTAveragedSource::GetTypeStatic()
1767
0
{
1768
0
    static const char *TYPE = "AveragedSource";
1769
0
    return TYPE;
1770
0
}
1771
1772
/************************************************************************/
1773
/*                            GetType()                                 */
1774
/************************************************************************/
1775
1776
const char *VRTAveragedSource::GetType() const
1777
0
{
1778
0
    return GetTypeStatic();
1779
0
}
1780
1781
/************************************************************************/
1782
/*                           SerializeToXML()                           */
1783
/************************************************************************/
1784
1785
CPLXMLNode *VRTAveragedSource::SerializeToXML(const char *pszVRTPath)
1786
1787
0
{
1788
0
    CPLXMLNode *const psSrc = VRTSimpleSource::SerializeToXML(pszVRTPath);
1789
1790
0
    if (psSrc == nullptr)
1791
0
        return nullptr;
1792
1793
0
    CPLFree(psSrc->pszValue);
1794
0
    psSrc->pszValue = CPLStrdup(GetTypeStatic());
1795
1796
0
    return psSrc;
1797
0
}
1798
1799
/************************************************************************/
1800
/*                           SetNoDataValue()                           */
1801
/************************************************************************/
1802
1803
void VRTAveragedSource::SetNoDataValue(double dfNewNoDataValue)
1804
1805
0
{
1806
0
    if (dfNewNoDataValue == VRT_NODATA_UNSET)
1807
0
    {
1808
0
        m_bNoDataSet = FALSE;
1809
0
        m_dfNoDataValue = VRT_NODATA_UNSET;
1810
0
        return;
1811
0
    }
1812
1813
0
    m_bNoDataSet = TRUE;
1814
0
    m_dfNoDataValue = dfNewNoDataValue;
1815
0
}
1816
1817
/************************************************************************/
1818
/*                              RasterIO()                              */
1819
/************************************************************************/
1820
1821
CPLErr VRTAveragedSource::RasterIO(GDALDataType /*eVRTBandDataType*/, int nXOff,
1822
                                   int nYOff, int nXSize, int nYSize,
1823
                                   void *pData, int nBufXSize, int nBufYSize,
1824
                                   GDALDataType eBufType, GSpacing nPixelSpace,
1825
                                   GSpacing nLineSpace,
1826
                                   GDALRasterIOExtraArg *psExtraArgIn,
1827
                                   WorkingState & /*oWorkingState*/)
1828
1829
0
{
1830
0
    GDALRasterIOExtraArg sExtraArg;
1831
0
    INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1832
0
    GDALRasterIOExtraArg *psExtraArg = &sExtraArg;
1833
1834
0
    double dfXOff = nXOff;
1835
0
    double dfYOff = nYOff;
1836
0
    double dfXSize = nXSize;
1837
0
    double dfYSize = nYSize;
1838
0
    if (psExtraArgIn != nullptr && psExtraArgIn->bFloatingPointWindowValidity)
1839
0
    {
1840
0
        dfXOff = psExtraArgIn->dfXOff;
1841
0
        dfYOff = psExtraArgIn->dfYOff;
1842
0
        dfXSize = psExtraArgIn->dfXSize;
1843
0
        dfYSize = psExtraArgIn->dfYSize;
1844
0
    }
1845
1846
    // The window we will actually request from the source raster band.
1847
0
    double dfReqXOff = 0.0;
1848
0
    double dfReqYOff = 0.0;
1849
0
    double dfReqXSize = 0.0;
1850
0
    double dfReqYSize = 0.0;
1851
0
    int nReqXOff = 0;
1852
0
    int nReqYOff = 0;
1853
0
    int nReqXSize = 0;
1854
0
    int nReqYSize = 0;
1855
1856
    // The window we will actual set _within_ the pData buffer.
1857
0
    int nOutXOff = 0;
1858
0
    int nOutYOff = 0;
1859
0
    int nOutXSize = 0;
1860
0
    int nOutYSize = 0;
1861
1862
0
    bool bError = false;
1863
0
    if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
1864
0
                         &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
1865
0
                         &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
1866
0
                         &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
1867
0
    {
1868
0
        return bError ? CE_Failure : CE_None;
1869
0
    }
1870
1871
0
    auto l_band = GetRasterBand();
1872
0
    if (!l_band)
1873
0
        return CE_Failure;
1874
1875
    /* -------------------------------------------------------------------- */
1876
    /*      Allocate a temporary buffer to whole the full resolution        */
1877
    /*      data from the area of interest.                                 */
1878
    /* -------------------------------------------------------------------- */
1879
0
    float *const pafSrc = static_cast<float *>(
1880
0
        VSI_MALLOC3_VERBOSE(sizeof(float), nReqXSize, nReqYSize));
1881
0
    if (pafSrc == nullptr)
1882
0
    {
1883
0
        return CE_Failure;
1884
0
    }
1885
1886
    /* -------------------------------------------------------------------- */
1887
    /*      Load it.                                                        */
1888
    /* -------------------------------------------------------------------- */
1889
0
    if (!m_osResampling.empty())
1890
0
    {
1891
0
        psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
1892
0
    }
1893
0
    else if (psExtraArgIn != nullptr)
1894
0
    {
1895
0
        psExtraArg->eResampleAlg = psExtraArgIn->eResampleAlg;
1896
0
    }
1897
1898
0
    psExtraArg->bFloatingPointWindowValidity = TRUE;
1899
0
    psExtraArg->dfXOff = dfReqXOff;
1900
0
    psExtraArg->dfYOff = dfReqYOff;
1901
0
    psExtraArg->dfXSize = dfReqXSize;
1902
0
    psExtraArg->dfYSize = dfReqYSize;
1903
0
    if (psExtraArgIn)
1904
0
    {
1905
0
        psExtraArg->pfnProgress = psExtraArgIn->pfnProgress;
1906
0
        psExtraArg->pProgressData = psExtraArgIn->pProgressData;
1907
0
    }
1908
1909
0
    const CPLErr eErr = l_band->RasterIO(
1910
0
        GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize, pafSrc, nReqXSize,
1911
0
        nReqYSize, GDT_Float32, 0, 0, psExtraArg);
1912
1913
0
    if (eErr != CE_None)
1914
0
    {
1915
0
        VSIFree(pafSrc);
1916
0
        return eErr;
1917
0
    }
1918
1919
    /* -------------------------------------------------------------------- */
1920
    /*      Do the averaging.                                               */
1921
    /* -------------------------------------------------------------------- */
1922
0
    for (int iBufLine = nOutYOff; iBufLine < nOutYOff + nOutYSize; iBufLine++)
1923
0
    {
1924
0
        const double dfYDst =
1925
0
            (iBufLine / static_cast<double>(nBufYSize)) * nYSize + nYOff;
1926
1927
0
        for (int iBufPixel = nOutXOff; iBufPixel < nOutXOff + nOutXSize;
1928
0
             iBufPixel++)
1929
0
        {
1930
0
            double dfXSrcStart, dfXSrcEnd, dfYSrcStart, dfYSrcEnd;
1931
0
            int iXSrcStart, iYSrcStart, iXSrcEnd, iYSrcEnd;
1932
1933
0
            const double dfXDst =
1934
0
                (iBufPixel / static_cast<double>(nBufXSize)) * nXSize + nXOff;
1935
1936
            // Compute the source image rectangle needed for this pixel.
1937
0
            DstToSrc(dfXDst, dfYDst, dfXSrcStart, dfYSrcStart);
1938
0
            DstToSrc(dfXDst + 1.0, dfYDst + 1.0, dfXSrcEnd, dfYSrcEnd);
1939
1940
            // Convert to integers, assuming that the center of the source
1941
            // pixel must be in our rect to get included.
1942
0
            if (dfXSrcEnd >= dfXSrcStart + 1)
1943
0
            {
1944
0
                iXSrcStart = static_cast<int>(floor(dfXSrcStart + 0.5));
1945
0
                iXSrcEnd = static_cast<int>(floor(dfXSrcEnd + 0.5));
1946
0
            }
1947
0
            else
1948
0
            {
1949
                /* If the resampling factor is less than 100%, the distance */
1950
                /* between the source pixel is < 1, so we stick to nearest */
1951
                /* neighbour */
1952
0
                iXSrcStart = static_cast<int>(floor(dfXSrcStart));
1953
0
                iXSrcEnd = iXSrcStart + 1;
1954
0
            }
1955
0
            if (dfYSrcEnd >= dfYSrcStart + 1)
1956
0
            {
1957
0
                iYSrcStart = static_cast<int>(floor(dfYSrcStart + 0.5));
1958
0
                iYSrcEnd = static_cast<int>(floor(dfYSrcEnd + 0.5));
1959
0
            }
1960
0
            else
1961
0
            {
1962
0
                iYSrcStart = static_cast<int>(floor(dfYSrcStart));
1963
0
                iYSrcEnd = iYSrcStart + 1;
1964
0
            }
1965
1966
            // Transform into the coordinate system of the source *buffer*
1967
0
            iXSrcStart -= nReqXOff;
1968
0
            iYSrcStart -= nReqYOff;
1969
0
            iXSrcEnd -= nReqXOff;
1970
0
            iYSrcEnd -= nReqYOff;
1971
1972
0
            double dfSum = 0.0;
1973
0
            int nPixelCount = 0;
1974
1975
0
            for (int iY = iYSrcStart; iY < iYSrcEnd; iY++)
1976
0
            {
1977
0
                if (iY < 0 || iY >= nReqYSize)
1978
0
                    continue;
1979
1980
0
                for (int iX = iXSrcStart; iX < iXSrcEnd; iX++)
1981
0
                {
1982
0
                    if (iX < 0 || iX >= nReqXSize)
1983
0
                        continue;
1984
1985
0
                    const float fSampledValue =
1986
0
                        pafSrc[iX + static_cast<size_t>(iY) * nReqXSize];
1987
0
                    if (std::isnan(fSampledValue))
1988
0
                        continue;
1989
1990
0
                    if (m_bNoDataSet &&
1991
0
                        GDALIsValueInRange<float>(m_dfNoDataValue) &&
1992
0
                        ARE_REAL_EQUAL(fSampledValue,
1993
0
                                       static_cast<float>(m_dfNoDataValue)))
1994
0
                        continue;
1995
1996
0
                    nPixelCount++;
1997
0
                    dfSum += pafSrc[iX + static_cast<size_t>(iY) * nReqXSize];
1998
0
                }
1999
0
            }
2000
2001
0
            if (nPixelCount == 0)
2002
0
                continue;
2003
2004
            // Compute output value.
2005
0
            const float dfOutputValue = static_cast<float>(dfSum / nPixelCount);
2006
2007
            // Put it in the output buffer.
2008
0
            GByte *pDstLocation =
2009
0
                static_cast<GByte *>(pData) + nPixelSpace * iBufPixel +
2010
0
                static_cast<GPtrDiff_t>(nLineSpace) * iBufLine;
2011
2012
0
            if (eBufType == GDT_Byte)
2013
0
                *pDstLocation = static_cast<GByte>(
2014
0
                    std::min(255.0, std::max(0.0, dfOutputValue + 0.5)));
2015
0
            else
2016
0
                GDALCopyWords(&dfOutputValue, GDT_Float32, 4, pDstLocation,
2017
0
                              eBufType, 8, 1);
2018
0
        }
2019
0
    }
2020
2021
0
    VSIFree(pafSrc);
2022
2023
0
    if (psExtraArg->pfnProgress)
2024
0
        psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
2025
2026
0
    return CE_None;
2027
0
}
2028
2029
/************************************************************************/
2030
/*                             GetMinimum()                             */
2031
/************************************************************************/
2032
2033
double VRTAveragedSource::GetMinimum(int /* nXSize */, int /* nYSize */,
2034
                                     int *pbSuccess)
2035
0
{
2036
0
    *pbSuccess = FALSE;
2037
0
    return 0.0;
2038
0
}
2039
2040
/************************************************************************/
2041
/*                             GetMaximum()                             */
2042
/************************************************************************/
2043
2044
double VRTAveragedSource::GetMaximum(int /* nXSize */, int /* nYSize */,
2045
                                     int *pbSuccess)
2046
0
{
2047
0
    *pbSuccess = FALSE;
2048
0
    return 0.0;
2049
0
}
2050
2051
/************************************************************************/
2052
/*                            GetHistogram()                            */
2053
/************************************************************************/
2054
2055
CPLErr VRTAveragedSource::GetHistogram(
2056
    int /* nXSize */, int /* nYSize */, double /* dfMin */, double /* dfMax */,
2057
    int /* nBuckets */, GUIntBig * /* panHistogram */,
2058
    int /* bIncludeOutOfRange */, int /* bApproxOK */,
2059
    GDALProgressFunc /* pfnProgress */, void * /* pProgressData */)
2060
0
{
2061
0
    return CE_Failure;
2062
0
}
2063
2064
/************************************************************************/
2065
/* ==================================================================== */
2066
/*                     VRTNoDataFromMaskSource                          */
2067
/* ==================================================================== */
2068
/************************************************************************/
2069
2070
/************************************************************************/
2071
/*                     VRTNoDataFromMaskSource()                        */
2072
/************************************************************************/
2073
2074
VRTNoDataFromMaskSource::VRTNoDataFromMaskSource()
2075
0
{
2076
0
}
2077
2078
/************************************************************************/
2079
/*                              XMLInit()                               */
2080
/************************************************************************/
2081
2082
CPLErr
2083
VRTNoDataFromMaskSource::XMLInit(const CPLXMLNode *psSrc,
2084
                                 const char *pszVRTPath,
2085
                                 VRTMapSharedResources &oMapSharedSources)
2086
2087
0
{
2088
    /* -------------------------------------------------------------------- */
2089
    /*      Do base initialization.                                         */
2090
    /* -------------------------------------------------------------------- */
2091
0
    {
2092
0
        const CPLErr eErr =
2093
0
            VRTSimpleSource::XMLInit(psSrc, pszVRTPath, oMapSharedSources);
2094
0
        if (eErr != CE_None)
2095
0
            return eErr;
2096
0
    }
2097
2098
0
    if (const char *pszNODATA = CPLGetXMLValue(psSrc, "NODATA", nullptr))
2099
0
    {
2100
0
        m_bNoDataSet = true;
2101
0
        m_dfNoDataValue = CPLAtofM(pszNODATA);
2102
0
    }
2103
2104
0
    m_dfMaskValueThreshold =
2105
0
        CPLAtofM(CPLGetXMLValue(psSrc, "MaskValueThreshold", "0"));
2106
2107
0
    if (const char *pszRemappedValue =
2108
0
            CPLGetXMLValue(psSrc, "RemappedValue", nullptr))
2109
0
    {
2110
0
        m_bHasRemappedValue = true;
2111
0
        m_dfRemappedValue = CPLAtofM(pszRemappedValue);
2112
0
    }
2113
2114
0
    return CE_None;
2115
0
}
2116
2117
/************************************************************************/
2118
/*                           GetTypeStatic()                            */
2119
/************************************************************************/
2120
2121
const char *VRTNoDataFromMaskSource::GetTypeStatic()
2122
0
{
2123
0
    static const char *TYPE = "NoDataFromMaskSource";
2124
0
    return TYPE;
2125
0
}
2126
2127
/************************************************************************/
2128
/*                            GetType()                                 */
2129
/************************************************************************/
2130
2131
const char *VRTNoDataFromMaskSource::GetType() const
2132
0
{
2133
0
    return GetTypeStatic();
2134
0
}
2135
2136
/************************************************************************/
2137
/*                           SerializeToXML()                           */
2138
/************************************************************************/
2139
2140
CPLXMLNode *VRTNoDataFromMaskSource::SerializeToXML(const char *pszVRTPath)
2141
2142
0
{
2143
0
    CPLXMLNode *const psSrc = VRTSimpleSource::SerializeToXML(pszVRTPath);
2144
2145
0
    if (psSrc == nullptr)
2146
0
        return nullptr;
2147
2148
0
    CPLFree(psSrc->pszValue);
2149
0
    psSrc->pszValue = CPLStrdup(GetTypeStatic());
2150
2151
0
    if (m_bNoDataSet)
2152
0
    {
2153
0
        CPLSetXMLValue(psSrc, "MaskValueThreshold",
2154
0
                       CPLSPrintf("%.17g", m_dfMaskValueThreshold));
2155
2156
0
        GDALDataType eBandDT = GDT_Unknown;
2157
0
        double dfNoDataValue = m_dfNoDataValue;
2158
0
        const auto kMaxFloat = std::numeric_limits<float>::max();
2159
0
        if (std::fabs(std::fabs(m_dfNoDataValue) - kMaxFloat) <
2160
0
            1e-10 * kMaxFloat)
2161
0
        {
2162
0
            auto l_band = GetRasterBand();
2163
0
            if (l_band)
2164
0
            {
2165
0
                eBandDT = l_band->GetRasterDataType();
2166
0
                if (eBandDT == GDT_Float32)
2167
0
                {
2168
0
                    dfNoDataValue =
2169
0
                        GDALAdjustNoDataCloseToFloatMax(m_dfNoDataValue);
2170
0
                }
2171
0
            }
2172
0
        }
2173
0
        CPLSetXMLValue(psSrc, "NODATA",
2174
0
                       VRTSerializeNoData(dfNoDataValue, eBandDT, 18).c_str());
2175
0
    }
2176
2177
0
    if (m_bHasRemappedValue)
2178
0
    {
2179
0
        CPLSetXMLValue(psSrc, "RemappedValue",
2180
0
                       CPLSPrintf("%.17g", m_dfRemappedValue));
2181
0
    }
2182
2183
0
    return psSrc;
2184
0
}
2185
2186
/************************************************************************/
2187
/*                           SetParameters()                            */
2188
/************************************************************************/
2189
2190
void VRTNoDataFromMaskSource::SetParameters(double dfNoDataValue,
2191
                                            double dfMaskValueThreshold)
2192
0
{
2193
0
    m_bNoDataSet = true;
2194
0
    m_dfNoDataValue = dfNoDataValue;
2195
0
    m_dfMaskValueThreshold = dfMaskValueThreshold;
2196
0
    if (!m_bHasRemappedValue)
2197
0
        m_dfRemappedValue = m_dfNoDataValue;
2198
0
}
2199
2200
/************************************************************************/
2201
/*                           SetParameters()                            */
2202
/************************************************************************/
2203
2204
void VRTNoDataFromMaskSource::SetParameters(double dfNoDataValue,
2205
                                            double dfMaskValueThreshold,
2206
                                            double dfRemappedValue)
2207
0
{
2208
0
    SetParameters(dfNoDataValue, dfMaskValueThreshold);
2209
0
    m_bHasRemappedValue = true;
2210
0
    m_dfRemappedValue = dfRemappedValue;
2211
0
}
2212
2213
/************************************************************************/
2214
/*                              RasterIO()                              */
2215
/************************************************************************/
2216
2217
CPLErr VRTNoDataFromMaskSource::RasterIO(
2218
    GDALDataType eVRTBandDataType, int nXOff, int nYOff, int nXSize, int nYSize,
2219
    void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
2220
    GSpacing nPixelSpace, GSpacing nLineSpace,
2221
    GDALRasterIOExtraArg *psExtraArgIn, WorkingState &oWorkingState)
2222
2223
0
{
2224
0
    if (!m_bNoDataSet)
2225
0
    {
2226
0
        return VRTSimpleSource::RasterIO(eVRTBandDataType, nXOff, nYOff, nXSize,
2227
0
                                         nYSize, pData, nBufXSize, nBufYSize,
2228
0
                                         eBufType, nPixelSpace, nLineSpace,
2229
0
                                         psExtraArgIn, oWorkingState);
2230
0
    }
2231
2232
0
    GDALRasterIOExtraArg sExtraArg;
2233
0
    INIT_RASTERIO_EXTRA_ARG(sExtraArg);
2234
0
    GDALRasterIOExtraArg *psExtraArg = &sExtraArg;
2235
2236
0
    double dfXOff = nXOff;
2237
0
    double dfYOff = nYOff;
2238
0
    double dfXSize = nXSize;
2239
0
    double dfYSize = nYSize;
2240
0
    if (psExtraArgIn != nullptr && psExtraArgIn->bFloatingPointWindowValidity)
2241
0
    {
2242
0
        dfXOff = psExtraArgIn->dfXOff;
2243
0
        dfYOff = psExtraArgIn->dfYOff;
2244
0
        dfXSize = psExtraArgIn->dfXSize;
2245
0
        dfYSize = psExtraArgIn->dfYSize;
2246
0
    }
2247
2248
    // The window we will actually request from the source raster band.
2249
0
    double dfReqXOff = 0.0;
2250
0
    double dfReqYOff = 0.0;
2251
0
    double dfReqXSize = 0.0;
2252
0
    double dfReqYSize = 0.0;
2253
0
    int nReqXOff = 0;
2254
0
    int nReqYOff = 0;
2255
0
    int nReqXSize = 0;
2256
0
    int nReqYSize = 0;
2257
2258
    // The window we will actual set _within_ the pData buffer.
2259
0
    int nOutXOff = 0;
2260
0
    int nOutYOff = 0;
2261
0
    int nOutXSize = 0;
2262
0
    int nOutYSize = 0;
2263
2264
0
    bool bError = false;
2265
0
    if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
2266
0
                         &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
2267
0
                         &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
2268
0
                         &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
2269
0
    {
2270
0
        return bError ? CE_Failure : CE_None;
2271
0
    }
2272
2273
0
    auto l_band = GetRasterBand();
2274
0
    if (!l_band)
2275
0
        return CE_Failure;
2276
2277
    /* -------------------------------------------------------------------- */
2278
    /*      Allocate temporary buffer(s).                                   */
2279
    /* -------------------------------------------------------------------- */
2280
0
    const auto eSrcBandDT = l_band->GetRasterDataType();
2281
0
    const int nSrcBandDTSize = GDALGetDataTypeSizeBytes(eSrcBandDT);
2282
0
    const auto eSrcMaskBandDT = l_band->GetMaskBand()->GetRasterDataType();
2283
0
    const int nSrcMaskBandDTSize = GDALGetDataTypeSizeBytes(eSrcMaskBandDT);
2284
0
    double dfRemappedValue = m_dfRemappedValue;
2285
0
    if (!m_bHasRemappedValue)
2286
0
    {
2287
0
        if (eSrcBandDT == GDT_Byte &&
2288
0
            m_dfNoDataValue >= std::numeric_limits<GByte>::min() &&
2289
0
            m_dfNoDataValue <= std::numeric_limits<GByte>::max() &&
2290
0
            static_cast<int>(m_dfNoDataValue) == m_dfNoDataValue)
2291
0
        {
2292
0
            if (m_dfNoDataValue == std::numeric_limits<GByte>::max())
2293
0
                dfRemappedValue = m_dfNoDataValue - 1;
2294
0
            else
2295
0
                dfRemappedValue = m_dfNoDataValue + 1;
2296
0
        }
2297
0
        else if (eSrcBandDT == GDT_UInt16 &&
2298
0
                 m_dfNoDataValue >= std::numeric_limits<uint16_t>::min() &&
2299
0
                 m_dfNoDataValue <= std::numeric_limits<uint16_t>::max() &&
2300
0
                 static_cast<int>(m_dfNoDataValue) == m_dfNoDataValue)
2301
0
        {
2302
0
            if (m_dfNoDataValue == std::numeric_limits<uint16_t>::max())
2303
0
                dfRemappedValue = m_dfNoDataValue - 1;
2304
0
            else
2305
0
                dfRemappedValue = m_dfNoDataValue + 1;
2306
0
        }
2307
0
        else if (eSrcBandDT == GDT_Int16 &&
2308
0
                 m_dfNoDataValue >= std::numeric_limits<int16_t>::min() &&
2309
0
                 m_dfNoDataValue <= std::numeric_limits<int16_t>::max() &&
2310
0
                 static_cast<int>(m_dfNoDataValue) == m_dfNoDataValue)
2311
0
        {
2312
0
            if (m_dfNoDataValue == std::numeric_limits<int16_t>::max())
2313
0
                dfRemappedValue = m_dfNoDataValue - 1;
2314
0
            else
2315
0
                dfRemappedValue = m_dfNoDataValue + 1;
2316
0
        }
2317
0
        else
2318
0
        {
2319
0
            constexpr double EPS = 1e-3;
2320
0
            if (m_dfNoDataValue == 0)
2321
0
                dfRemappedValue = EPS;
2322
0
            else
2323
0
                dfRemappedValue = m_dfNoDataValue * (1 + EPS);
2324
0
        }
2325
0
    }
2326
0
    const bool bByteOptim =
2327
0
        (eSrcBandDT == GDT_Byte && eBufType == GDT_Byte &&
2328
0
         eSrcMaskBandDT == GDT_Byte && m_dfMaskValueThreshold >= 0 &&
2329
0
         m_dfMaskValueThreshold <= 255 &&
2330
0
         static_cast<int>(m_dfMaskValueThreshold) == m_dfMaskValueThreshold &&
2331
0
         m_dfNoDataValue >= 0 && m_dfNoDataValue <= 255 &&
2332
0
         static_cast<int>(m_dfNoDataValue) == m_dfNoDataValue &&
2333
0
         dfRemappedValue >= 0 && dfRemappedValue <= 255 &&
2334
0
         static_cast<int>(dfRemappedValue) == dfRemappedValue);
2335
0
    GByte *pabyWrkBuffer;
2336
0
    try
2337
0
    {
2338
0
        if (bByteOptim && nOutXOff == 0 && nOutYOff == 0 &&
2339
0
            nOutXSize == nBufXSize && nOutYSize == nBufYSize &&
2340
0
            eSrcBandDT == eBufType && nPixelSpace == nSrcBandDTSize &&
2341
0
            nLineSpace == nPixelSpace * nBufXSize)
2342
0
        {
2343
0
            pabyWrkBuffer = static_cast<GByte *>(pData);
2344
0
        }
2345
0
        else
2346
0
        {
2347
0
            oWorkingState.m_abyWrkBuffer.resize(static_cast<size_t>(nOutXSize) *
2348
0
                                                nOutYSize * nSrcBandDTSize);
2349
0
            pabyWrkBuffer =
2350
0
                reinterpret_cast<GByte *>(oWorkingState.m_abyWrkBuffer.data());
2351
0
        }
2352
0
        oWorkingState.m_abyWrkBufferMask.resize(static_cast<size_t>(nOutXSize) *
2353
0
                                                nOutYSize * nSrcMaskBandDTSize);
2354
0
    }
2355
0
    catch (const std::exception &)
2356
0
    {
2357
0
        CPLError(CE_Failure, CPLE_OutOfMemory,
2358
0
                 "Out of memory when allocating buffers");
2359
0
        return CE_Failure;
2360
0
    }
2361
2362
    /* -------------------------------------------------------------------- */
2363
    /*      Load data.                                                      */
2364
    /* -------------------------------------------------------------------- */
2365
0
    if (!m_osResampling.empty())
2366
0
    {
2367
0
        psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
2368
0
    }
2369
0
    else if (psExtraArgIn != nullptr)
2370
0
    {
2371
0
        psExtraArg->eResampleAlg = psExtraArgIn->eResampleAlg;
2372
0
    }
2373
2374
0
    psExtraArg->bFloatingPointWindowValidity = TRUE;
2375
0
    psExtraArg->dfXOff = dfReqXOff;
2376
0
    psExtraArg->dfYOff = dfReqYOff;
2377
0
    psExtraArg->dfXSize = dfReqXSize;
2378
0
    psExtraArg->dfYSize = dfReqYSize;
2379
0
    if (psExtraArgIn)
2380
0
    {
2381
0
        psExtraArg->pfnProgress = psExtraArgIn->pfnProgress;
2382
0
        psExtraArg->pProgressData = psExtraArgIn->pProgressData;
2383
0
    }
2384
2385
0
    if (l_band->RasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
2386
0
                         pabyWrkBuffer, nOutXSize, nOutYSize, eSrcBandDT, 0, 0,
2387
0
                         psExtraArg) != CE_None)
2388
0
    {
2389
0
        return CE_Failure;
2390
0
    }
2391
2392
0
    if (l_band->GetMaskBand()->RasterIO(
2393
0
            GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
2394
0
            oWorkingState.m_abyWrkBufferMask.data(), nOutXSize, nOutYSize,
2395
0
            eSrcMaskBandDT, 0, 0, psExtraArg) != CE_None)
2396
0
    {
2397
0
        return CE_Failure;
2398
0
    }
2399
2400
    /* -------------------------------------------------------------------- */
2401
    /*      Do the processing.                                              */
2402
    /* -------------------------------------------------------------------- */
2403
2404
0
    GByte *const pabyOut = static_cast<GByte *>(pData) +
2405
0
                           nPixelSpace * nOutXOff +
2406
0
                           static_cast<GPtrDiff_t>(nLineSpace) * nOutYOff;
2407
0
    if (bByteOptim)
2408
0
    {
2409
        // Special case when everything fits on Byte
2410
0
        const GByte nMaskValueThreshold =
2411
0
            static_cast<GByte>(m_dfMaskValueThreshold);
2412
0
        const GByte nNoDataValue = static_cast<GByte>(m_dfNoDataValue);
2413
0
        const GByte nRemappedValue = static_cast<GByte>(dfRemappedValue);
2414
0
        size_t nSrcIdx = 0;
2415
0
        for (int iY = 0; iY < nOutYSize; iY++)
2416
0
        {
2417
0
            GSpacing nDstOffset = iY * nLineSpace;
2418
0
            for (int iX = 0; iX < nOutXSize; iX++)
2419
0
            {
2420
0
                const GByte nMaskVal =
2421
0
                    oWorkingState.m_abyWrkBufferMask[nSrcIdx];
2422
0
                if (nMaskVal <= nMaskValueThreshold)
2423
0
                {
2424
0
                    pabyOut[static_cast<GPtrDiff_t>(nDstOffset)] = nNoDataValue;
2425
0
                }
2426
0
                else
2427
0
                {
2428
0
                    if (pabyWrkBuffer[nSrcIdx] == nNoDataValue)
2429
0
                    {
2430
0
                        pabyOut[static_cast<GPtrDiff_t>(nDstOffset)] =
2431
0
                            nRemappedValue;
2432
0
                    }
2433
0
                    else
2434
0
                    {
2435
0
                        pabyOut[static_cast<GPtrDiff_t>(nDstOffset)] =
2436
0
                            pabyWrkBuffer[nSrcIdx];
2437
0
                    }
2438
0
                }
2439
0
                nDstOffset += nPixelSpace;
2440
0
                nSrcIdx++;
2441
0
            }
2442
0
        }
2443
0
    }
2444
0
    else
2445
0
    {
2446
0
        size_t nSrcIdx = 0;
2447
0
        double dfMaskVal = 0;
2448
0
        const int nBufDTSize = GDALGetDataTypeSizeBytes(eBufType);
2449
0
        std::vector<GByte> abyDstNoData(nBufDTSize);
2450
0
        GDALCopyWords(&m_dfNoDataValue, GDT_Float64, 0, abyDstNoData.data(),
2451
0
                      eBufType, 0, 1);
2452
0
        std::vector<GByte> abyRemappedValue(nBufDTSize);
2453
0
        GDALCopyWords(&dfRemappedValue, GDT_Float64, 0, abyRemappedValue.data(),
2454
0
                      eBufType, 0, 1);
2455
0
        for (int iY = 0; iY < nOutYSize; iY++)
2456
0
        {
2457
0
            GSpacing nDstOffset = iY * nLineSpace;
2458
0
            for (int iX = 0; iX < nOutXSize; iX++)
2459
0
            {
2460
0
                if (eSrcMaskBandDT == GDT_Byte)
2461
0
                {
2462
0
                    dfMaskVal = oWorkingState.m_abyWrkBufferMask[nSrcIdx];
2463
0
                }
2464
0
                else
2465
0
                {
2466
0
                    GDALCopyWords(oWorkingState.m_abyWrkBufferMask.data() +
2467
0
                                      nSrcIdx * nSrcMaskBandDTSize,
2468
0
                                  eSrcMaskBandDT, 0, &dfMaskVal, GDT_Float64, 0,
2469
0
                                  1);
2470
0
                }
2471
0
                void *const pDst =
2472
0
                    pabyOut + static_cast<GPtrDiff_t>(nDstOffset);
2473
0
                if (!(dfMaskVal > m_dfMaskValueThreshold))
2474
0
                {
2475
0
                    memcpy(pDst, abyDstNoData.data(), nBufDTSize);
2476
0
                }
2477
0
                else
2478
0
                {
2479
0
                    const void *const pSrc =
2480
0
                        pabyWrkBuffer + nSrcIdx * nSrcBandDTSize;
2481
0
                    if (eSrcBandDT == eBufType)
2482
0
                    {
2483
                        // coverity[overrun-buffer-arg]
2484
0
                        memcpy(pDst, pSrc, nBufDTSize);
2485
0
                    }
2486
0
                    else
2487
0
                    {
2488
0
                        GDALCopyWords(pSrc, eSrcBandDT, 0, pDst, eBufType, 0,
2489
0
                                      1);
2490
0
                    }
2491
0
                    if (memcmp(pDst, abyDstNoData.data(), nBufDTSize) == 0)
2492
0
                        memcpy(pDst, abyRemappedValue.data(), nBufDTSize);
2493
0
                }
2494
0
                nDstOffset += nPixelSpace;
2495
0
                nSrcIdx++;
2496
0
            }
2497
0
        }
2498
0
    }
2499
2500
0
    if (psExtraArg->pfnProgress)
2501
0
        psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
2502
2503
0
    return CE_None;
2504
0
}
2505
2506
/************************************************************************/
2507
/*                             GetMinimum()                             */
2508
/************************************************************************/
2509
2510
double VRTNoDataFromMaskSource::GetMinimum(int /* nXSize */, int /* nYSize */,
2511
                                           int *pbSuccess)
2512
0
{
2513
0
    *pbSuccess = FALSE;
2514
0
    return 0.0;
2515
0
}
2516
2517
/************************************************************************/
2518
/*                             GetMaximum()                             */
2519
/************************************************************************/
2520
2521
double VRTNoDataFromMaskSource::GetMaximum(int /* nXSize */, int /* nYSize */,
2522
                                           int *pbSuccess)
2523
0
{
2524
0
    *pbSuccess = FALSE;
2525
0
    return 0.0;
2526
0
}
2527
2528
/************************************************************************/
2529
/*                            GetHistogram()                            */
2530
/************************************************************************/
2531
2532
CPLErr VRTNoDataFromMaskSource::GetHistogram(
2533
    int /* nXSize */, int /* nYSize */, double /* dfMin */, double /* dfMax */,
2534
    int /* nBuckets */, GUIntBig * /* panHistogram */,
2535
    int /* bIncludeOutOfRange */, int /* bApproxOK */,
2536
    GDALProgressFunc /* pfnProgress */, void * /* pProgressData */)
2537
0
{
2538
0
    return CE_Failure;
2539
0
}
2540
2541
/************************************************************************/
2542
/* ==================================================================== */
2543
/*                          VRTComplexSource                            */
2544
/* ==================================================================== */
2545
/************************************************************************/
2546
2547
/************************************************************************/
2548
/*                          VRTComplexSource()                          */
2549
/************************************************************************/
2550
2551
VRTComplexSource::VRTComplexSource(const VRTComplexSource *poSrcSource,
2552
                                   double dfXDstRatio, double dfYDstRatio)
2553
0
    : VRTSimpleSource(poSrcSource, dfXDstRatio, dfYDstRatio),
2554
0
      m_nProcessingFlags(poSrcSource->m_nProcessingFlags),
2555
0
      m_dfNoDataValue(poSrcSource->m_dfNoDataValue),
2556
0
      m_osNoDataValueOri(poSrcSource->m_osNoDataValueOri),
2557
0
      m_dfScaleOff(poSrcSource->m_dfScaleOff),
2558
0
      m_dfScaleRatio(poSrcSource->m_dfScaleRatio),
2559
0
      m_bSrcMinMaxDefined(poSrcSource->m_bSrcMinMaxDefined),
2560
0
      m_dfSrcMin(poSrcSource->m_dfSrcMin), m_dfSrcMax(poSrcSource->m_dfSrcMax),
2561
0
      m_dfDstMin(poSrcSource->m_dfDstMin), m_dfDstMax(poSrcSource->m_dfDstMax),
2562
0
      m_dfExponent(poSrcSource->m_dfExponent), m_bClip(poSrcSource->m_bClip),
2563
0
      m_nColorTableComponent(poSrcSource->m_nColorTableComponent),
2564
0
      m_adfLUTInputs(poSrcSource->m_adfLUTInputs),
2565
0
      m_adfLUTOutputs(poSrcSource->m_adfLUTOutputs)
2566
0
{
2567
0
}
2568
2569
/************************************************************************/
2570
/*                           GetTypeStatic()                            */
2571
/************************************************************************/
2572
2573
const char *VRTComplexSource::GetTypeStatic()
2574
0
{
2575
0
    static const char *TYPE = "ComplexSource";
2576
0
    return TYPE;
2577
0
}
2578
2579
/************************************************************************/
2580
/*                            GetType()                                 */
2581
/************************************************************************/
2582
2583
const char *VRTComplexSource::GetType() const
2584
0
{
2585
0
    return GetTypeStatic();
2586
0
}
2587
2588
/************************************************************************/
2589
/*                           SetNoDataValue()                           */
2590
/************************************************************************/
2591
2592
void VRTComplexSource::SetNoDataValue(double dfNewNoDataValue)
2593
2594
0
{
2595
0
    if (dfNewNoDataValue == VRT_NODATA_UNSET)
2596
0
    {
2597
0
        m_nProcessingFlags &= ~PROCESSING_FLAG_NODATA;
2598
0
        m_dfNoDataValue = VRT_NODATA_UNSET;
2599
0
        return;
2600
0
    }
2601
2602
0
    m_nProcessingFlags |= PROCESSING_FLAG_NODATA;
2603
0
    m_dfNoDataValue = dfNewNoDataValue;
2604
0
}
2605
2606
/************************************************************************/
2607
/*                      GetAdjustedNoDataValue()                        */
2608
/************************************************************************/
2609
2610
double VRTComplexSource::GetAdjustedNoDataValue() const
2611
0
{
2612
0
    if ((m_nProcessingFlags & PROCESSING_FLAG_NODATA) != 0)
2613
0
    {
2614
0
        auto l_band = GetRasterBand();
2615
0
        if (l_band && l_band->GetRasterDataType() == GDT_Float32)
2616
0
        {
2617
0
            return GDALAdjustNoDataCloseToFloatMax(m_dfNoDataValue);
2618
0
        }
2619
0
    }
2620
0
    return m_dfNoDataValue;
2621
0
}
2622
2623
/************************************************************************/
2624
/*                           SerializeToXML()                           */
2625
/************************************************************************/
2626
2627
CPLXMLNode *VRTComplexSource::SerializeToXML(const char *pszVRTPath)
2628
2629
0
{
2630
0
    CPLXMLNode *psSrc = VRTSimpleSource::SerializeToXML(pszVRTPath);
2631
2632
0
    if (psSrc == nullptr)
2633
0
        return nullptr;
2634
2635
0
    CPLFree(psSrc->pszValue);
2636
0
    psSrc->pszValue = CPLStrdup(GetTypeStatic());
2637
2638
0
    if ((m_nProcessingFlags & PROCESSING_FLAG_USE_MASK_BAND) != 0)
2639
0
    {
2640
0
        CPLSetXMLValue(psSrc, "UseMaskBand", "true");
2641
0
    }
2642
2643
0
    if ((m_nProcessingFlags & PROCESSING_FLAG_NODATA) != 0)
2644
0
    {
2645
0
        if (!m_osNoDataValueOri.empty() && GetRasterBandNoOpen() == nullptr)
2646
0
        {
2647
0
            CPLSetXMLValue(psSrc, "NODATA", m_osNoDataValueOri.c_str());
2648
0
        }
2649
0
        else
2650
0
        {
2651
0
            GDALDataType eBandDT = GDT_Unknown;
2652
0
            double dfNoDataValue = m_dfNoDataValue;
2653
0
            const auto kMaxFloat = std::numeric_limits<float>::max();
2654
0
            if (std::fabs(std::fabs(m_dfNoDataValue) - kMaxFloat) <
2655
0
                1e-10 * kMaxFloat)
2656
0
            {
2657
0
                auto l_band = GetRasterBand();
2658
0
                if (l_band)
2659
0
                {
2660
0
                    dfNoDataValue = GetAdjustedNoDataValue();
2661
0
                    eBandDT = l_band->GetRasterDataType();
2662
0
                }
2663
0
            }
2664
0
            CPLSetXMLValue(
2665
0
                psSrc, "NODATA",
2666
0
                VRTSerializeNoData(dfNoDataValue, eBandDT, 18).c_str());
2667
0
        }
2668
0
    }
2669
2670
0
    if ((m_nProcessingFlags & PROCESSING_FLAG_SCALING_LINEAR) != 0)
2671
0
    {
2672
0
        CPLSetXMLValue(psSrc, "ScaleOffset", CPLSPrintf("%g", m_dfScaleOff));
2673
0
        CPLSetXMLValue(psSrc, "ScaleRatio", CPLSPrintf("%g", m_dfScaleRatio));
2674
0
    }
2675
0
    else if ((m_nProcessingFlags & PROCESSING_FLAG_SCALING_EXPONENTIAL) != 0)
2676
0
    {
2677
0
        CPLSetXMLValue(psSrc, "Exponent", CPLSPrintf("%g", m_dfExponent));
2678
0
        if (m_bSrcMinMaxDefined)
2679
0
        {
2680
0
            CPLSetXMLValue(psSrc, "SrcMin", CPLSPrintf("%g", m_dfSrcMin));
2681
0
            CPLSetXMLValue(psSrc, "SrcMax", CPLSPrintf("%g", m_dfSrcMax));
2682
0
        }
2683
0
        CPLSetXMLValue(psSrc, "DstMin", CPLSPrintf("%g", m_dfDstMin));
2684
0
        CPLSetXMLValue(psSrc, "DstMax", CPLSPrintf("%g", m_dfDstMax));
2685
0
        CPLSetXMLValue(psSrc, "Clip", m_bClip ? "true" : "false");
2686
0
    }
2687
2688
0
    if (!m_adfLUTInputs.empty())
2689
0
    {
2690
        // Make sure we print with sufficient precision to address really close
2691
        // entries (#6422).
2692
0
        CPLString osLUT;
2693
0
        if (m_adfLUTInputs.size() >= 2 &&
2694
0
            CPLString().Printf("%g", m_adfLUTInputs[0]) ==
2695
0
                CPLString().Printf("%g", m_adfLUTInputs[1]))
2696
0
        {
2697
0
            osLUT = CPLString().Printf("%.17g:%g", m_adfLUTInputs[0],
2698
0
                                       m_adfLUTOutputs[0]);
2699
0
        }
2700
0
        else
2701
0
        {
2702
0
            osLUT = CPLString().Printf("%g:%g", m_adfLUTInputs[0],
2703
0
                                       m_adfLUTOutputs[0]);
2704
0
        }
2705
0
        for (size_t i = 1; i < m_adfLUTInputs.size(); i++)
2706
0
        {
2707
0
            if (CPLString().Printf("%g", m_adfLUTInputs[i]) ==
2708
0
                    CPLString().Printf("%g", m_adfLUTInputs[i - 1]) ||
2709
0
                (i + 1 < m_adfLUTInputs.size() &&
2710
0
                 CPLString().Printf("%g", m_adfLUTInputs[i]) ==
2711
0
                     CPLString().Printf("%g", m_adfLUTInputs[i + 1])))
2712
0
            {
2713
                // TODO(schwehr): An explanation of the 18 would be helpful.
2714
                // Can someone distill the issue down to a quick comment?
2715
                // https://trac.osgeo.org/gdal/ticket/6422
2716
0
                osLUT += CPLString().Printf(",%.17g:%g", m_adfLUTInputs[i],
2717
0
                                            m_adfLUTOutputs[i]);
2718
0
            }
2719
0
            else
2720
0
            {
2721
0
                osLUT += CPLString().Printf(",%g:%g", m_adfLUTInputs[i],
2722
0
                                            m_adfLUTOutputs[i]);
2723
0
            }
2724
0
        }
2725
0
        CPLSetXMLValue(psSrc, "LUT", osLUT);
2726
0
    }
2727
2728
0
    if (m_nColorTableComponent)
2729
0
    {
2730
0
        CPLSetXMLValue(psSrc, "ColorTableComponent",
2731
0
                       CPLSPrintf("%d", m_nColorTableComponent));
2732
0
    }
2733
2734
0
    return psSrc;
2735
0
}
2736
2737
/************************************************************************/
2738
/*                              XMLInit()                               */
2739
/************************************************************************/
2740
2741
CPLErr VRTComplexSource::XMLInit(const CPLXMLNode *psSrc,
2742
                                 const char *pszVRTPath,
2743
                                 VRTMapSharedResources &oMapSharedSources)
2744
2745
0
{
2746
    /* -------------------------------------------------------------------- */
2747
    /*      Do base initialization.                                         */
2748
    /* -------------------------------------------------------------------- */
2749
0
    {
2750
0
        const CPLErr eErr =
2751
0
            VRTSimpleSource::XMLInit(psSrc, pszVRTPath, oMapSharedSources);
2752
0
        if (eErr != CE_None)
2753
0
            return eErr;
2754
0
    }
2755
2756
    /* -------------------------------------------------------------------- */
2757
    /*      Complex parameters.                                             */
2758
    /* -------------------------------------------------------------------- */
2759
0
    const char *pszScaleOffset = CPLGetXMLValue(psSrc, "ScaleOffset", nullptr);
2760
0
    const char *pszScaleRatio = CPLGetXMLValue(psSrc, "ScaleRatio", nullptr);
2761
0
    if (pszScaleOffset || pszScaleRatio)
2762
0
    {
2763
0
        m_nProcessingFlags |= PROCESSING_FLAG_SCALING_LINEAR;
2764
0
        if (pszScaleOffset)
2765
0
            m_dfScaleOff = CPLAtof(pszScaleOffset);
2766
0
        if (pszScaleRatio)
2767
0
            m_dfScaleRatio = CPLAtof(pszScaleRatio);
2768
0
    }
2769
0
    else if (CPLGetXMLValue(psSrc, "Exponent", nullptr) != nullptr &&
2770
0
             CPLGetXMLValue(psSrc, "DstMin", nullptr) != nullptr &&
2771
0
             CPLGetXMLValue(psSrc, "DstMax", nullptr) != nullptr)
2772
0
    {
2773
0
        m_nProcessingFlags |= PROCESSING_FLAG_SCALING_EXPONENTIAL;
2774
0
        m_dfExponent = CPLAtof(CPLGetXMLValue(psSrc, "Exponent", "1.0"));
2775
2776
0
        const char *pszSrcMin = CPLGetXMLValue(psSrc, "SrcMin", nullptr);
2777
0
        const char *pszSrcMax = CPLGetXMLValue(psSrc, "SrcMax", nullptr);
2778
0
        if (pszSrcMin && pszSrcMax)
2779
0
        {
2780
0
            m_dfSrcMin = CPLAtof(pszSrcMin);
2781
0
            m_dfSrcMax = CPLAtof(pszSrcMax);
2782
0
            m_bSrcMinMaxDefined = true;
2783
0
        }
2784
2785
0
        m_dfDstMin = CPLAtof(CPLGetXMLValue(psSrc, "DstMin", "0.0"));
2786
0
        m_dfDstMax = CPLAtof(CPLGetXMLValue(psSrc, "DstMax", "0.0"));
2787
0
        m_bClip = CPLTestBool(CPLGetXMLValue(psSrc, "Clip", "true"));
2788
0
    }
2789
2790
0
    if (const char *pszNODATA = CPLGetXMLValue(psSrc, "NODATA", nullptr))
2791
0
    {
2792
0
        m_nProcessingFlags |= PROCESSING_FLAG_NODATA;
2793
0
        m_osNoDataValueOri = pszNODATA;
2794
0
        m_dfNoDataValue = CPLAtofM(m_osNoDataValueOri.c_str());
2795
0
    }
2796
2797
0
    const char *pszUseMaskBand = CPLGetXMLValue(psSrc, "UseMaskBand", nullptr);
2798
0
    if (pszUseMaskBand && CPLTestBool(pszUseMaskBand))
2799
0
    {
2800
0
        m_nProcessingFlags |= PROCESSING_FLAG_USE_MASK_BAND;
2801
0
    }
2802
2803
0
    const char *pszLUT = CPLGetXMLValue(psSrc, "LUT", nullptr);
2804
0
    if (pszLUT)
2805
0
    {
2806
0
        const CPLStringList aosValues(
2807
0
            CSLTokenizeString2(pszLUT, ",:", CSLT_ALLOWEMPTYTOKENS));
2808
2809
0
        const int nLUTItemCount = aosValues.size() / 2;
2810
0
        try
2811
0
        {
2812
0
            m_adfLUTInputs.resize(nLUTItemCount);
2813
0
            m_adfLUTOutputs.resize(nLUTItemCount);
2814
0
        }
2815
0
        catch (const std::bad_alloc &e)
2816
0
        {
2817
0
            CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
2818
0
            m_adfLUTInputs.clear();
2819
0
            m_adfLUTOutputs.clear();
2820
0
            return CE_Failure;
2821
0
        }
2822
2823
0
        for (int nIndex = 0; nIndex < nLUTItemCount; nIndex++)
2824
0
        {
2825
0
            m_adfLUTInputs[nIndex] = CPLAtof(aosValues[nIndex * 2]);
2826
0
            m_adfLUTOutputs[nIndex] = CPLAtof(aosValues[nIndex * 2 + 1]);
2827
2828
            // Enforce the requirement that the LUT input array is
2829
            // monotonically non-decreasing.
2830
0
            if (std::isnan(m_adfLUTInputs[nIndex]) && nIndex != 0)
2831
0
            {
2832
0
                CPLError(CE_Failure, CPLE_AppDefined,
2833
0
                         "A Not-A-Number (NaN) source value should be the "
2834
0
                         "first one of the LUT.");
2835
0
                m_adfLUTInputs.clear();
2836
0
                m_adfLUTOutputs.clear();
2837
0
                return CE_Failure;
2838
0
            }
2839
0
            else if (nIndex > 0 &&
2840
0
                     m_adfLUTInputs[nIndex] < m_adfLUTInputs[nIndex - 1])
2841
0
            {
2842
0
                CPLError(CE_Failure, CPLE_AppDefined,
2843
0
                         "Source values of the LUT are not listed in a "
2844
0
                         "monotonically non-decreasing order");
2845
0
                m_adfLUTInputs.clear();
2846
0
                m_adfLUTOutputs.clear();
2847
0
                return CE_Failure;
2848
0
            }
2849
0
        }
2850
2851
0
        m_nProcessingFlags |= PROCESSING_FLAG_LUT;
2852
0
    }
2853
2854
0
    const char *pszColorTableComponent =
2855
0
        CPLGetXMLValue(psSrc, "ColorTableComponent", nullptr);
2856
0
    if (pszColorTableComponent)
2857
0
    {
2858
0
        m_nColorTableComponent = atoi(pszColorTableComponent);
2859
0
        m_nProcessingFlags |= PROCESSING_FLAG_COLOR_TABLE_EXPANSION;
2860
0
    }
2861
2862
0
    return CE_None;
2863
0
}
2864
2865
/************************************************************************/
2866
/*                              LookupValue()                           */
2867
/************************************************************************/
2868
2869
double VRTComplexSource::LookupValue(double dfInput)
2870
0
{
2871
0
    auto beginIter = m_adfLUTInputs.begin();
2872
0
    auto endIter = m_adfLUTInputs.end();
2873
0
    size_t offset = 0;
2874
0
    if (std::isnan(m_adfLUTInputs[0]))
2875
0
    {
2876
0
        if (std::isnan(dfInput) || m_adfLUTInputs.size() == 1)
2877
0
            return m_adfLUTOutputs[0];
2878
0
        ++beginIter;
2879
0
        offset = 1;
2880
0
    }
2881
2882
    // Find the index of the first element in the LUT input array that
2883
    // is not smaller than the input value.
2884
0
    const size_t i =
2885
0
        offset +
2886
0
        std::distance(beginIter, std::lower_bound(beginIter, endIter, dfInput));
2887
2888
0
    if (i == offset)
2889
0
        return m_adfLUTOutputs[offset];
2890
2891
    // If the index is beyond the end of the LUT input array, the input
2892
    // value is larger than all the values in the array.
2893
0
    if (i == m_adfLUTInputs.size())
2894
0
        return m_adfLUTOutputs.back();
2895
2896
0
    if (m_adfLUTInputs[i] == dfInput)
2897
0
        return m_adfLUTOutputs[i];
2898
2899
    // Otherwise, interpolate.
2900
0
    return m_adfLUTOutputs[i - 1] +
2901
0
           (dfInput - m_adfLUTInputs[i - 1]) *
2902
0
               ((m_adfLUTOutputs[i] - m_adfLUTOutputs[i - 1]) /
2903
0
                (m_adfLUTInputs[i] - m_adfLUTInputs[i - 1]));
2904
0
}
2905
2906
/************************************************************************/
2907
/*                         SetLinearScaling()                           */
2908
/************************************************************************/
2909
2910
void VRTComplexSource::SetLinearScaling(double dfOffset, double dfScale)
2911
0
{
2912
0
    m_nProcessingFlags &= ~PROCESSING_FLAG_SCALING_EXPONENTIAL;
2913
0
    m_nProcessingFlags |= PROCESSING_FLAG_SCALING_LINEAR;
2914
0
    m_dfScaleOff = dfOffset;
2915
0
    m_dfScaleRatio = dfScale;
2916
0
}
2917
2918
/************************************************************************/
2919
/*                         SetPowerScaling()                           */
2920
/************************************************************************/
2921
2922
void VRTComplexSource::SetPowerScaling(double dfExponentIn, double dfSrcMinIn,
2923
                                       double dfSrcMaxIn, double dfDstMinIn,
2924
                                       double dfDstMaxIn, bool bClip)
2925
0
{
2926
0
    m_nProcessingFlags &= ~PROCESSING_FLAG_SCALING_LINEAR;
2927
0
    m_nProcessingFlags |= PROCESSING_FLAG_SCALING_EXPONENTIAL;
2928
0
    m_dfExponent = dfExponentIn;
2929
0
    m_dfSrcMin = dfSrcMinIn;
2930
0
    m_dfSrcMax = dfSrcMaxIn;
2931
0
    m_dfDstMin = dfDstMinIn;
2932
0
    m_dfDstMax = dfDstMaxIn;
2933
0
    m_bSrcMinMaxDefined = true;
2934
0
    m_bClip = bClip;
2935
0
}
2936
2937
/************************************************************************/
2938
/*                    SetColorTableComponent()                          */
2939
/************************************************************************/
2940
2941
void VRTComplexSource::SetColorTableComponent(int nComponent)
2942
0
{
2943
0
    m_nProcessingFlags |= PROCESSING_FLAG_COLOR_TABLE_EXPANSION;
2944
0
    m_nColorTableComponent = nComponent;
2945
0
}
2946
2947
/************************************************************************/
2948
/*                              RasterIO()                              */
2949
/************************************************************************/
2950
2951
CPLErr VRTComplexSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
2952
                                  int nYOff, int nXSize, int nYSize,
2953
                                  void *pData, int nBufXSize, int nBufYSize,
2954
                                  GDALDataType eBufType, GSpacing nPixelSpace,
2955
                                  GSpacing nLineSpace,
2956
                                  GDALRasterIOExtraArg *psExtraArgIn,
2957
                                  WorkingState &oWorkingState)
2958
2959
0
{
2960
0
    GDALRasterIOExtraArg sExtraArg;
2961
0
    INIT_RASTERIO_EXTRA_ARG(sExtraArg);
2962
0
    GDALRasterIOExtraArg *psExtraArg = &sExtraArg;
2963
2964
0
    double dfXOff = nXOff;
2965
0
    double dfYOff = nYOff;
2966
0
    double dfXSize = nXSize;
2967
0
    double dfYSize = nYSize;
2968
0
    if (psExtraArgIn != nullptr && psExtraArgIn->bFloatingPointWindowValidity)
2969
0
    {
2970
0
        dfXOff = psExtraArgIn->dfXOff;
2971
0
        dfYOff = psExtraArgIn->dfYOff;
2972
0
        dfXSize = psExtraArgIn->dfXSize;
2973
0
        dfYSize = psExtraArgIn->dfYSize;
2974
0
    }
2975
2976
    // The window we will actually request from the source raster band.
2977
0
    double dfReqXOff = 0.0;
2978
0
    double dfReqYOff = 0.0;
2979
0
    double dfReqXSize = 0.0;
2980
0
    double dfReqYSize = 0.0;
2981
0
    int nReqXOff = 0;
2982
0
    int nReqYOff = 0;
2983
0
    int nReqXSize = 0;
2984
0
    int nReqYSize = 0;
2985
2986
    // The window we will actual set _within_ the pData buffer.
2987
0
    int nOutXOff = 0;
2988
0
    int nOutYOff = 0;
2989
0
    int nOutXSize = 0;
2990
0
    int nOutYSize = 0;
2991
2992
0
    bool bError = false;
2993
0
    if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
2994
0
                         &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
2995
0
                         &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
2996
0
                         &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
2997
0
    {
2998
0
        return bError ? CE_Failure : CE_None;
2999
0
    }
3000
#if DEBUG_VERBOSE
3001
    CPLDebug("VRT",
3002
             "nXOff=%d, nYOff=%d, nXSize=%d, nYSize=%d, nBufXSize=%d, "
3003
             "nBufYSize=%d,\n"
3004
             "dfReqXOff=%g, dfReqYOff=%g, dfReqXSize=%g, dfReqYSize=%g,\n"
3005
             "nReqXOff=%d, nReqYOff=%d, nReqXSize=%d, nReqYSize=%d,\n"
3006
             "nOutXOff=%d, nOutYOff=%d, nOutXSize=%d, nOutYSize=%d",
3007
             nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize, dfReqXOff,
3008
             dfReqYOff, dfReqXSize, dfReqYSize, nReqXOff, nReqYOff, nReqXSize,
3009
             nReqYSize, nOutXOff, nOutYOff, nOutXSize, nOutYSize);
3010
#endif
3011
3012
0
    auto poSourceBand = GetRasterBand();
3013
0
    if (!poSourceBand)
3014
0
        return CE_Failure;
3015
3016
0
    if (!m_osResampling.empty())
3017
0
    {
3018
0
        psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
3019
0
    }
3020
0
    else if (psExtraArgIn != nullptr)
3021
0
    {
3022
0
        psExtraArg->eResampleAlg = psExtraArgIn->eResampleAlg;
3023
0
    }
3024
0
    psExtraArg->bFloatingPointWindowValidity = TRUE;
3025
0
    psExtraArg->dfXOff = dfReqXOff;
3026
0
    psExtraArg->dfYOff = dfReqYOff;
3027
0
    psExtraArg->dfXSize = dfReqXSize;
3028
0
    psExtraArg->dfYSize = dfReqYSize;
3029
0
    if (psExtraArgIn)
3030
0
    {
3031
0
        psExtraArg->pfnProgress = psExtraArgIn->pfnProgress;
3032
0
        psExtraArg->pProgressData = psExtraArgIn->pProgressData;
3033
0
    }
3034
3035
0
    GByte *const pabyOut = static_cast<GByte *>(pData) +
3036
0
                           nPixelSpace * nOutXOff +
3037
0
                           static_cast<GPtrDiff_t>(nLineSpace) * nOutYOff;
3038
0
    const auto eSourceType = poSourceBand->GetRasterDataType();
3039
0
    if (m_nProcessingFlags == PROCESSING_FLAG_NODATA)
3040
0
    {
3041
        // Optimization if doing only nodata processing
3042
0
        if (eSourceType == GDT_Byte)
3043
0
        {
3044
0
            if (!GDALIsValueInRange<GByte>(m_dfNoDataValue))
3045
0
            {
3046
0
                return VRTSimpleSource::RasterIO(
3047
0
                    eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData,
3048
0
                    nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace,
3049
0
                    psExtraArgIn, oWorkingState);
3050
0
            }
3051
3052
0
            return RasterIOProcessNoData<GByte, GDT_Byte>(
3053
0
                poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3054
0
                nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3055
0
                nLineSpace, psExtraArg, oWorkingState);
3056
0
        }
3057
0
        else if (eSourceType == GDT_Int16)
3058
0
        {
3059
0
            if (!GDALIsValueInRange<int16_t>(m_dfNoDataValue))
3060
0
            {
3061
0
                return VRTSimpleSource::RasterIO(
3062
0
                    eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData,
3063
0
                    nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace,
3064
0
                    psExtraArgIn, oWorkingState);
3065
0
            }
3066
3067
0
            return RasterIOProcessNoData<int16_t, GDT_Int16>(
3068
0
                poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3069
0
                nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3070
0
                nLineSpace, psExtraArg, oWorkingState);
3071
0
        }
3072
0
        else if (eSourceType == GDT_UInt16)
3073
0
        {
3074
0
            if (!GDALIsValueInRange<uint16_t>(m_dfNoDataValue))
3075
0
            {
3076
0
                return VRTSimpleSource::RasterIO(
3077
0
                    eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData,
3078
0
                    nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace,
3079
0
                    psExtraArgIn, oWorkingState);
3080
0
            }
3081
3082
0
            return RasterIOProcessNoData<uint16_t, GDT_UInt16>(
3083
0
                poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3084
0
                nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3085
0
                nLineSpace, psExtraArg, oWorkingState);
3086
0
        }
3087
0
    }
3088
3089
0
    const bool bIsComplex = CPL_TO_BOOL(GDALDataTypeIsComplex(eBufType));
3090
0
    CPLErr eErr;
3091
    // For Int32, float32 isn't sufficiently precise as working data type
3092
0
    if (eVRTBandDataType == GDT_CInt32 || eVRTBandDataType == GDT_CFloat64 ||
3093
0
        eVRTBandDataType == GDT_Int32 || eVRTBandDataType == GDT_UInt32 ||
3094
0
        eVRTBandDataType == GDT_Int64 || eVRTBandDataType == GDT_UInt64 ||
3095
0
        eVRTBandDataType == GDT_Float64 || eSourceType == GDT_Int32 ||
3096
0
        eSourceType == GDT_UInt32 || eSourceType == GDT_Int64 ||
3097
0
        eSourceType == GDT_UInt64)
3098
0
    {
3099
0
        eErr = RasterIOInternal<double>(
3100
0
            poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3101
0
            nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3102
0
            nLineSpace, psExtraArg, bIsComplex ? GDT_CFloat64 : GDT_Float64,
3103
0
            oWorkingState);
3104
0
    }
3105
0
    else
3106
0
    {
3107
0
        eErr = RasterIOInternal<float>(
3108
0
            poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3109
0
            nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3110
0
            nLineSpace, psExtraArg, bIsComplex ? GDT_CFloat32 : GDT_Float32,
3111
0
            oWorkingState);
3112
0
    }
3113
3114
0
    if (psExtraArg->pfnProgress)
3115
0
        psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
3116
3117
0
    return eErr;
3118
0
}
3119
3120
/************************************************************************/
3121
/*                              hasZeroByte()                           */
3122
/************************************************************************/
3123
3124
CPL_NOSANITIZE_UNSIGNED_INT_OVERFLOW
3125
static inline bool hasZeroByte(uint32_t v)
3126
0
{
3127
    // Cf https://graphics.stanford.edu/~seander/bithacks.html#ZeroInWord
3128
0
    return (((v)-0x01010101U) & ~(v)&0x80808080U) != 0;
3129
0
}
3130
3131
/************************************************************************/
3132
/*                                CopyWord()                            */
3133
/************************************************************************/
3134
3135
template <class SrcType>
3136
static void CopyWord(const SrcType *pSrcVal, GDALDataType eSrcType,
3137
                     void *pDstVal, GDALDataType eDstType)
3138
0
{
3139
0
    switch (eDstType)
3140
0
    {
3141
0
        case GDT_Byte:
3142
0
            GDALCopyWord(*pSrcVal, *static_cast<uint8_t *>(pDstVal));
3143
0
            break;
3144
0
        case GDT_Int8:
3145
0
            GDALCopyWord(*pSrcVal, *static_cast<int8_t *>(pDstVal));
3146
0
            break;
3147
0
        case GDT_UInt16:
3148
0
            GDALCopyWord(*pSrcVal, *static_cast<uint16_t *>(pDstVal));
3149
0
            break;
3150
0
        case GDT_Int16:
3151
0
            GDALCopyWord(*pSrcVal, *static_cast<int16_t *>(pDstVal));
3152
0
            break;
3153
0
        case GDT_UInt32:
3154
0
            GDALCopyWord(*pSrcVal, *static_cast<uint32_t *>(pDstVal));
3155
0
            break;
3156
0
        case GDT_Int32:
3157
0
            GDALCopyWord(*pSrcVal, *static_cast<int32_t *>(pDstVal));
3158
0
            break;
3159
0
        case GDT_UInt64:
3160
0
            GDALCopyWord(*pSrcVal, *static_cast<uint64_t *>(pDstVal));
3161
0
            break;
3162
0
        case GDT_Int64:
3163
0
            GDALCopyWord(*pSrcVal, *static_cast<int64_t *>(pDstVal));
3164
0
            break;
3165
0
        case GDT_Float32:
3166
0
            GDALCopyWord(*pSrcVal, *static_cast<float *>(pDstVal));
3167
0
            break;
3168
0
        case GDT_Float64:
3169
0
            GDALCopyWord(*pSrcVal, *static_cast<double *>(pDstVal));
3170
0
            break;
3171
0
        default:
3172
0
            GDALCopyWords(pSrcVal, eSrcType, 0, pDstVal, eDstType, 0, 1);
3173
0
            break;
3174
0
    }
3175
0
}
Unexecuted instantiation: vrtsources.cpp:void CopyWord<unsigned char>(unsigned char const*, GDALDataType, void*, GDALDataType)
Unexecuted instantiation: vrtsources.cpp:void CopyWord<short>(short const*, GDALDataType, void*, GDALDataType)
Unexecuted instantiation: vrtsources.cpp:void CopyWord<unsigned short>(unsigned short const*, GDALDataType, void*, GDALDataType)
Unexecuted instantiation: vrtsources.cpp:void CopyWord<double>(double const*, GDALDataType, void*, GDALDataType)
Unexecuted instantiation: vrtsources.cpp:void CopyWord<float>(float const*, GDALDataType, void*, GDALDataType)
3176
3177
/************************************************************************/
3178
/*                       RasterIOProcessNoData()                        */
3179
/************************************************************************/
3180
3181
// This method is an optimization of the generic RasterIOInternal()
3182
// that deals with a VRTComplexSource with only a NODATA value in it, and
3183
// no other processing flags.
3184
3185
// nReqXOff, nReqYOff, nReqXSize, nReqYSize are expressed in source band
3186
// referential.
3187
template <class SourceDT, GDALDataType eSourceType>
3188
CPLErr VRTComplexSource::RasterIOProcessNoData(
3189
    GDALRasterBand *poSourceBand, GDALDataType eVRTBandDataType, int nReqXOff,
3190
    int nReqYOff, int nReqXSize, int nReqYSize, void *pData, int nOutXSize,
3191
    int nOutYSize, GDALDataType eBufType, GSpacing nPixelSpace,
3192
    GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg,
3193
    WorkingState &oWorkingState)
3194
0
{
3195
0
    CPLAssert(m_nProcessingFlags == PROCESSING_FLAG_NODATA);
3196
0
    CPLAssert(GDALIsValueInRange<SourceDT>(m_dfNoDataValue));
3197
3198
    /* -------------------------------------------------------------------- */
3199
    /*      Read into a temporary buffer.                                   */
3200
    /* -------------------------------------------------------------------- */
3201
0
    try
3202
0
    {
3203
        // Cannot overflow since pData should at least have that number of
3204
        // elements
3205
0
        const size_t nPixelCount = static_cast<size_t>(nOutXSize) * nOutYSize;
3206
0
        if (nPixelCount >
3207
0
            static_cast<size_t>(std::numeric_limits<ptrdiff_t>::max()) /
3208
0
                sizeof(SourceDT))
3209
0
        {
3210
0
            CPLError(CE_Failure, CPLE_OutOfMemory,
3211
0
                     "Too large temporary buffer");
3212
0
            return CE_Failure;
3213
0
        }
3214
0
        oWorkingState.m_abyWrkBuffer.resize(sizeof(SourceDT) * nPixelCount);
3215
0
    }
3216
0
    catch (const std::bad_alloc &e)
3217
0
    {
3218
0
        CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
3219
0
        return CE_Failure;
3220
0
    }
3221
0
    const auto paSrcData =
3222
0
        reinterpret_cast<const SourceDT *>(oWorkingState.m_abyWrkBuffer.data());
3223
3224
0
    const GDALRIOResampleAlg eResampleAlgBack = psExtraArg->eResampleAlg;
3225
0
    if (!m_osResampling.empty())
3226
0
    {
3227
0
        psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
3228
0
    }
3229
3230
0
    const CPLErr eErr = poSourceBand->RasterIO(
3231
0
        GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
3232
0
        oWorkingState.m_abyWrkBuffer.data(), nOutXSize, nOutYSize, eSourceType,
3233
0
        sizeof(SourceDT), sizeof(SourceDT) * static_cast<GSpacing>(nOutXSize),
3234
0
        psExtraArg);
3235
0
    if (!m_osResampling.empty())
3236
0
        psExtraArg->eResampleAlg = eResampleAlgBack;
3237
3238
0
    if (eErr != CE_None)
3239
0
    {
3240
0
        return eErr;
3241
0
    }
3242
3243
0
    const auto nNoDataValue = static_cast<SourceDT>(m_dfNoDataValue);
3244
0
    size_t idxBuffer = 0;
3245
0
    if (eSourceType == eBufType &&
3246
0
        !GDALDataTypeIsConversionLossy(eSourceType, eVRTBandDataType))
3247
0
    {
3248
        // Most optimized case: the output type is the same as the source type,
3249
        // and conversion from the source type to the VRT band data type is
3250
        // not lossy
3251
0
        for (int iY = 0; iY < nOutYSize; iY++)
3252
0
        {
3253
0
            GByte *pDstLocation = static_cast<GByte *>(pData) +
3254
0
                                  static_cast<GPtrDiff_t>(nLineSpace) * iY;
3255
3256
0
            int iX = 0;
3257
0
            if (sizeof(SourceDT) == 1 && nPixelSpace == 1)
3258
0
            {
3259
                // Optimization to detect more quickly if source pixels are
3260
                // at nodata.
3261
0
                const GByte byNoDataValue = static_cast<GByte>(nNoDataValue);
3262
0
                const uint32_t wordNoData =
3263
0
                    (static_cast<uint32_t>(byNoDataValue) << 24) |
3264
0
                    (byNoDataValue << 16) | (byNoDataValue << 8) |
3265
0
                    byNoDataValue;
3266
3267
                // Warning: hasZeroByte() assumes WORD_SIZE = 4
3268
0
                constexpr int WORD_SIZE = 4;
3269
0
                for (; iX < nOutXSize - (WORD_SIZE - 1); iX += WORD_SIZE)
3270
0
                {
3271
0
                    uint32_t v;
3272
0
                    static_assert(sizeof(v) == WORD_SIZE,
3273
0
                                  "sizeof(v) == WORD_SIZE");
3274
0
                    memcpy(&v, paSrcData + idxBuffer, sizeof(v));
3275
                    // Cf https://graphics.stanford.edu/~seander/bithacks.html#ValueInWord
3276
0
                    if (!hasZeroByte(v ^ wordNoData))
3277
0
                    {
3278
                        // No bytes are at nodata
3279
0
                        memcpy(pDstLocation, &v, WORD_SIZE);
3280
0
                        idxBuffer += WORD_SIZE;
3281
0
                        pDstLocation += WORD_SIZE;
3282
0
                    }
3283
0
                    else if (v == wordNoData)
3284
0
                    {
3285
                        // All bytes are at nodata
3286
0
                        idxBuffer += WORD_SIZE;
3287
0
                        pDstLocation += WORD_SIZE;
3288
0
                    }
3289
0
                    else
3290
0
                    {
3291
                        // There are both bytes at nodata and valid bytes
3292
0
                        for (int k = 0; k < WORD_SIZE; ++k)
3293
0
                        {
3294
0
                            if (paSrcData[idxBuffer] != nNoDataValue)
3295
0
                            {
3296
0
                                memcpy(pDstLocation, &paSrcData[idxBuffer],
3297
0
                                       sizeof(SourceDT));
3298
0
                            }
3299
0
                            idxBuffer++;
3300
0
                            pDstLocation += nPixelSpace;
3301
0
                        }
3302
0
                    }
3303
0
                }
3304
0
            }
3305
3306
0
            for (; iX < nOutXSize;
3307
0
                 iX++, pDstLocation += nPixelSpace, idxBuffer++)
3308
0
            {
3309
0
                if (paSrcData[idxBuffer] != nNoDataValue)
3310
0
                {
3311
0
                    memcpy(pDstLocation, &paSrcData[idxBuffer],
3312
0
                           sizeof(SourceDT));
3313
0
                }
3314
0
            }
3315
0
        }
3316
0
    }
3317
0
    else if (!GDALDataTypeIsConversionLossy(eSourceType, eVRTBandDataType))
3318
0
    {
3319
        // Conversion from the source type to the VRT band data type is
3320
        // not lossy, so we can directly convert from the source type to
3321
        // the the output type
3322
0
        for (int iY = 0; iY < nOutYSize; iY++)
3323
0
        {
3324
0
            GByte *pDstLocation = static_cast<GByte *>(pData) +
3325
0
                                  static_cast<GPtrDiff_t>(nLineSpace) * iY;
3326
3327
0
            for (int iX = 0; iX < nOutXSize;
3328
0
                 iX++, pDstLocation += nPixelSpace, idxBuffer++)
3329
0
            {
3330
0
                if (paSrcData[idxBuffer] != nNoDataValue)
3331
0
                {
3332
0
                    CopyWord(&paSrcData[idxBuffer], eSourceType, pDstLocation,
3333
0
                             eBufType);
3334
0
                }
3335
0
            }
3336
0
        }
3337
0
    }
3338
0
    else
3339
0
    {
3340
0
        GByte abyTemp[2 * sizeof(double)];
3341
0
        for (int iY = 0; iY < nOutYSize; iY++)
3342
0
        {
3343
0
            GByte *pDstLocation = static_cast<GByte *>(pData) +
3344
0
                                  static_cast<GPtrDiff_t>(nLineSpace) * iY;
3345
3346
0
            for (int iX = 0; iX < nOutXSize;
3347
0
                 iX++, pDstLocation += nPixelSpace, idxBuffer++)
3348
0
            {
3349
0
                if (paSrcData[idxBuffer] != nNoDataValue)
3350
0
                {
3351
                    // Convert first to the VRTRasterBand data type
3352
                    // to get its clamping, before outputting to buffer data type
3353
0
                    CopyWord(&paSrcData[idxBuffer], eSourceType, abyTemp,
3354
0
                             eVRTBandDataType);
3355
0
                    GDALCopyWords(abyTemp, eVRTBandDataType, 0, pDstLocation,
3356
0
                                  eBufType, 0, 1);
3357
0
                }
3358
0
            }
3359
0
        }
3360
0
    }
3361
3362
0
    return CE_None;
3363
0
}
Unexecuted instantiation: CPLErr VRTComplexSource::RasterIOProcessNoData<unsigned char, (GDALDataType)1>(GDALRasterBand*, GDALDataType, int, int, int, int, void*, int, int, GDALDataType, long long, long long, GDALRasterIOExtraArg*, VRTSource::WorkingState&)
Unexecuted instantiation: CPLErr VRTComplexSource::RasterIOProcessNoData<short, (GDALDataType)3>(GDALRasterBand*, GDALDataType, int, int, int, int, void*, int, int, GDALDataType, long long, long long, GDALRasterIOExtraArg*, VRTSource::WorkingState&)
Unexecuted instantiation: CPLErr VRTComplexSource::RasterIOProcessNoData<unsigned short, (GDALDataType)2>(GDALRasterBand*, GDALDataType, int, int, int, int, void*, int, int, GDALDataType, long long, long long, GDALRasterIOExtraArg*, VRTSource::WorkingState&)
3364
3365
/************************************************************************/
3366
/*                          RasterIOInternal()                          */
3367
/************************************************************************/
3368
3369
// nReqXOff, nReqYOff, nReqXSize, nReqYSize are expressed in source band
3370
// referential.
3371
template <class WorkingDT>
3372
CPLErr VRTComplexSource::RasterIOInternal(
3373
    GDALRasterBand *poSourceBand, GDALDataType eVRTBandDataType, int nReqXOff,
3374
    int nReqYOff, int nReqXSize, int nReqYSize, void *pData, int nOutXSize,
3375
    int nOutYSize, GDALDataType eBufType, GSpacing nPixelSpace,
3376
    GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg,
3377
    GDALDataType eWrkDataType, WorkingState &oWorkingState)
3378
0
{
3379
0
    const GDALColorTable *poColorTable = nullptr;
3380
0
    const bool bIsComplex = CPL_TO_BOOL(GDALDataTypeIsComplex(eBufType));
3381
0
    const int nWordSize = GDALGetDataTypeSizeBytes(eWrkDataType);
3382
0
    assert(nWordSize != 0);
3383
3384
    // If no explicit <NODATA> is set, but UseMaskBand is set, and the band
3385
    // has a nodata value, then use it as if it was set as <NODATA>
3386
0
    int bNoDataSet = (m_nProcessingFlags & PROCESSING_FLAG_NODATA) != 0;
3387
0
    double dfNoDataValue = GetAdjustedNoDataValue();
3388
3389
0
    if ((m_nProcessingFlags & PROCESSING_FLAG_USE_MASK_BAND) != 0 &&
3390
0
        poSourceBand->GetMaskFlags() == GMF_NODATA)
3391
0
    {
3392
0
        dfNoDataValue = poSourceBand->GetNoDataValue(&bNoDataSet);
3393
0
    }
3394
3395
0
    const bool bNoDataSetIsNan = bNoDataSet && std::isnan(dfNoDataValue);
3396
0
    const bool bNoDataSetAndNotNan =
3397
0
        bNoDataSet && !std::isnan(dfNoDataValue) &&
3398
0
        GDALIsValueInRange<WorkingDT>(dfNoDataValue);
3399
0
    const auto fWorkingDataTypeNoData = static_cast<WorkingDT>(dfNoDataValue);
3400
3401
0
    const GByte *pabyMask = nullptr;
3402
0
    const WorkingDT *pafData = nullptr;
3403
0
    if ((m_nProcessingFlags & PROCESSING_FLAG_SCALING_LINEAR) != 0 &&
3404
0
        m_dfScaleRatio == 0 && bNoDataSet == FALSE &&
3405
0
        (m_nProcessingFlags & PROCESSING_FLAG_USE_MASK_BAND) == 0)
3406
0
    {
3407
        /* ------------------------------------------------------------------ */
3408
        /*      Optimization when writing a constant value */
3409
        /*      (used by the -addalpha option of gdalbuildvrt) */
3410
        /* ------------------------------------------------------------------ */
3411
        // Already set to NULL when defined.
3412
        // pafData = NULL;
3413
0
    }
3414
0
    else
3415
0
    {
3416
        /* ---------------------------------------------------------------- */
3417
        /*      Read into a temporary buffer.                               */
3418
        /* ---------------------------------------------------------------- */
3419
0
        const size_t nPixelCount = static_cast<size_t>(nOutXSize) * nOutYSize;
3420
0
        try
3421
0
        {
3422
            // Cannot overflow since pData should at least have that number of
3423
            // elements
3424
0
            if (nPixelCount >
3425
0
                static_cast<size_t>(std::numeric_limits<ptrdiff_t>::max()) /
3426
0
                    static_cast<size_t>(nWordSize))
3427
0
            {
3428
0
                CPLError(CE_Failure, CPLE_OutOfMemory,
3429
0
                         "Too large temporary buffer");
3430
0
                return CE_Failure;
3431
0
            }
3432
0
            oWorkingState.m_abyWrkBuffer.resize(nWordSize * nPixelCount);
3433
0
        }
3434
0
        catch (const std::bad_alloc &e)
3435
0
        {
3436
0
            CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
3437
0
            return CE_Failure;
3438
0
        }
3439
0
        pafData = reinterpret_cast<const WorkingDT *>(
3440
0
            oWorkingState.m_abyWrkBuffer.data());
3441
3442
0
        const GDALRIOResampleAlg eResampleAlgBack = psExtraArg->eResampleAlg;
3443
0
        if (!m_osResampling.empty())
3444
0
        {
3445
0
            psExtraArg->eResampleAlg =
3446
0
                GDALRasterIOGetResampleAlg(m_osResampling);
3447
0
        }
3448
3449
0
        const CPLErr eErr = poSourceBand->RasterIO(
3450
0
            GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
3451
0
            oWorkingState.m_abyWrkBuffer.data(), nOutXSize, nOutYSize,
3452
0
            eWrkDataType, nWordSize,
3453
0
            nWordSize * static_cast<GSpacing>(nOutXSize), psExtraArg);
3454
0
        if (!m_osResampling.empty())
3455
0
            psExtraArg->eResampleAlg = eResampleAlgBack;
3456
3457
0
        if (eErr != CE_None)
3458
0
        {
3459
0
            return eErr;
3460
0
        }
3461
3462
        // Allocate and read mask band if needed
3463
0
        if (!bNoDataSet &&
3464
0
            (m_nProcessingFlags & PROCESSING_FLAG_USE_MASK_BAND) != 0 &&
3465
0
            (poSourceBand->GetMaskFlags() != GMF_ALL_VALID ||
3466
0
             poSourceBand->GetColorInterpretation() == GCI_AlphaBand ||
3467
0
             GetMaskBandMainBand() != nullptr))
3468
0
        {
3469
0
            try
3470
0
            {
3471
0
                oWorkingState.m_abyWrkBufferMask.resize(nPixelCount);
3472
0
            }
3473
0
            catch (const std::exception &)
3474
0
            {
3475
0
                CPLError(CE_Failure, CPLE_OutOfMemory,
3476
0
                         "Out of memory when allocating mask buffer");
3477
0
                return CE_Failure;
3478
0
            }
3479
0
            pabyMask = reinterpret_cast<const GByte *>(
3480
0
                oWorkingState.m_abyWrkBufferMask.data());
3481
0
            auto poMaskBand =
3482
0
                (poSourceBand->GetColorInterpretation() == GCI_AlphaBand ||
3483
0
                 GetMaskBandMainBand() != nullptr)
3484
0
                    ? poSourceBand
3485
0
                    : poSourceBand->GetMaskBand();
3486
0
            if (poMaskBand->RasterIO(
3487
0
                    GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
3488
0
                    oWorkingState.m_abyWrkBufferMask.data(), nOutXSize,
3489
0
                    nOutYSize, GDT_Byte, 1, static_cast<GSpacing>(nOutXSize),
3490
0
                    psExtraArg) != CE_None)
3491
0
            {
3492
0
                return CE_Failure;
3493
0
            }
3494
0
        }
3495
3496
0
        if (m_nColorTableComponent != 0)
3497
0
        {
3498
0
            poColorTable = poSourceBand->GetColorTable();
3499
0
            if (poColorTable == nullptr)
3500
0
            {
3501
0
                CPLError(CE_Failure, CPLE_AppDefined,
3502
0
                         "Source band has no color table.");
3503
0
                return CE_Failure;
3504
0
            }
3505
0
        }
3506
0
    }
3507
3508
    /* -------------------------------------------------------------------- */
3509
    /*      Selectively copy into output buffer with nodata masking,        */
3510
    /*      and/or scaling.                                                 */
3511
    /* -------------------------------------------------------------------- */
3512
3513
0
    const bool bTwoStepDataTypeConversion =
3514
0
        CPL_TO_BOOL(
3515
0
            GDALDataTypeIsConversionLossy(eWrkDataType, eVRTBandDataType)) &&
3516
0
        !CPL_TO_BOOL(GDALDataTypeIsConversionLossy(eVRTBandDataType, eBufType));
3517
3518
0
    size_t idxBuffer = 0;
3519
0
    for (int iY = 0; iY < nOutYSize; iY++)
3520
0
    {
3521
0
        GByte *pDstLocation = static_cast<GByte *>(pData) +
3522
0
                              static_cast<GPtrDiff_t>(nLineSpace) * iY;
3523
3524
0
        for (int iX = 0; iX < nOutXSize;
3525
0
             iX++, pDstLocation += nPixelSpace, idxBuffer++)
3526
0
        {
3527
0
            WorkingDT afResult[2];
3528
0
            if (pafData && !bIsComplex)
3529
0
            {
3530
0
                WorkingDT fResult = pafData[idxBuffer];
3531
0
                if (bNoDataSetIsNan && std::isnan(fResult))
3532
0
                    continue;
3533
0
                if (bNoDataSetAndNotNan &&
3534
0
                    ARE_REAL_EQUAL(fResult, fWorkingDataTypeNoData))
3535
0
                    continue;
3536
0
                if (pabyMask && pabyMask[idxBuffer] == 0)
3537
0
                    continue;
3538
3539
0
                if (poColorTable)
3540
0
                {
3541
0
                    const GDALColorEntry *poEntry =
3542
0
                        poColorTable->GetColorEntry(static_cast<int>(fResult));
3543
0
                    if (poEntry)
3544
0
                    {
3545
0
                        if (m_nColorTableComponent == 1)
3546
0
                            fResult = poEntry->c1;
3547
0
                        else if (m_nColorTableComponent == 2)
3548
0
                            fResult = poEntry->c2;
3549
0
                        else if (m_nColorTableComponent == 3)
3550
0
                            fResult = poEntry->c3;
3551
0
                        else if (m_nColorTableComponent == 4)
3552
0
                            fResult = poEntry->c4;
3553
0
                    }
3554
0
                    else
3555
0
                    {
3556
0
                        CPLErrorOnce(CE_Failure, CPLE_AppDefined,
3557
0
                                     "No entry %d.", static_cast<int>(fResult));
3558
0
                        continue;
3559
0
                    }
3560
0
                }
3561
3562
0
                if ((m_nProcessingFlags & PROCESSING_FLAG_SCALING_LINEAR) != 0)
3563
0
                {
3564
0
                    fResult = static_cast<WorkingDT>(fResult * m_dfScaleRatio +
3565
0
                                                     m_dfScaleOff);
3566
0
                }
3567
0
                else if ((m_nProcessingFlags &
3568
0
                          PROCESSING_FLAG_SCALING_EXPONENTIAL) != 0)
3569
0
                {
3570
0
                    if (!m_bSrcMinMaxDefined)
3571
0
                    {
3572
0
                        int bSuccessMin = FALSE;
3573
0
                        int bSuccessMax = FALSE;
3574
0
                        double adfMinMax[2] = {
3575
0
                            poSourceBand->GetMinimum(&bSuccessMin),
3576
0
                            poSourceBand->GetMaximum(&bSuccessMax)};
3577
0
                        if ((bSuccessMin && bSuccessMax) ||
3578
0
                            poSourceBand->ComputeRasterMinMax(
3579
0
                                TRUE, adfMinMax) == CE_None)
3580
0
                        {
3581
0
                            m_dfSrcMin = adfMinMax[0];
3582
0
                            m_dfSrcMax = adfMinMax[1];
3583
0
                            m_bSrcMinMaxDefined = true;
3584
0
                        }
3585
0
                        else
3586
0
                        {
3587
0
                            CPLError(CE_Failure, CPLE_AppDefined,
3588
0
                                     "Cannot determine source min/max value");
3589
0
                            return CE_Failure;
3590
0
                        }
3591
0
                    }
3592
3593
0
                    double dfPowVal = (m_dfSrcMin == m_dfSrcMax)
3594
0
                                          ? 0
3595
0
                                          : (fResult - m_dfSrcMin) /
3596
0
                                                (m_dfSrcMax - m_dfSrcMin);
3597
0
                    if (m_bClip)
3598
0
                    {
3599
0
                        if (dfPowVal < 0.0)
3600
0
                            dfPowVal = 0.0;
3601
0
                        else if (dfPowVal > 1.0)
3602
0
                            dfPowVal = 1.0;
3603
0
                    }
3604
0
                    fResult =
3605
0
                        static_cast<WorkingDT>((m_dfDstMax - m_dfDstMin) *
3606
0
                                                   pow(dfPowVal, m_dfExponent) +
3607
0
                                               m_dfDstMin);
3608
0
                }
3609
3610
0
                if (!m_adfLUTInputs.empty())
3611
0
                    fResult = static_cast<WorkingDT>(LookupValue(fResult));
3612
3613
0
                if (m_nMaxValue != 0 && fResult > m_nMaxValue)
3614
0
                    fResult = static_cast<WorkingDT>(m_nMaxValue);
3615
3616
0
                afResult[0] = fResult;
3617
0
                afResult[1] = 0;
3618
0
            }
3619
0
            else if (pafData && bIsComplex)
3620
0
            {
3621
0
                afResult[0] = pafData[2 * idxBuffer];
3622
0
                afResult[1] = pafData[2 * idxBuffer + 1];
3623
3624
                // Do not use color table.
3625
0
                if ((m_nProcessingFlags & PROCESSING_FLAG_SCALING_LINEAR) != 0)
3626
0
                {
3627
0
                    afResult[0] = static_cast<WorkingDT>(
3628
0
                        afResult[0] * m_dfScaleRatio + m_dfScaleOff);
3629
0
                    afResult[1] = static_cast<WorkingDT>(
3630
0
                        afResult[1] * m_dfScaleRatio + m_dfScaleOff);
3631
0
                }
3632
3633
                /* Do not use LUT */
3634
0
            }
3635
0
            else
3636
0
            {
3637
0
                afResult[0] = static_cast<WorkingDT>(m_dfScaleOff);
3638
0
                afResult[1] = 0;
3639
3640
0
                if (!m_adfLUTInputs.empty())
3641
0
                    afResult[0] =
3642
0
                        static_cast<WorkingDT>(LookupValue(afResult[0]));
3643
3644
0
                if (m_nMaxValue != 0 && afResult[0] > m_nMaxValue)
3645
0
                    afResult[0] = static_cast<WorkingDT>(m_nMaxValue);
3646
0
            }
3647
3648
0
            if (eBufType == GDT_Byte && eVRTBandDataType == GDT_Byte)
3649
0
            {
3650
0
                *pDstLocation = static_cast<GByte>(std::min(
3651
0
                    255.0f,
3652
0
                    std::max(0.0f, static_cast<float>(afResult[0]) + 0.5f)));
3653
0
            }
3654
0
            else if (!bTwoStepDataTypeConversion)
3655
0
            {
3656
0
                CopyWord<WorkingDT>(afResult, eWrkDataType, pDstLocation,
3657
0
                                    eBufType);
3658
0
            }
3659
0
            else if (eVRTBandDataType == GDT_Float32 && eBufType == GDT_Float64)
3660
0
            {
3661
                // Particular case of the below 2-step conversion.
3662
                // Helps a bit for some geolocation based warping with Sentinel3
3663
                // data where the longitude/latitude arrays are Int32 bands,
3664
                // rescaled in VRT as Float32 and requested as Float64
3665
0
                float fVal;
3666
0
                GDALCopyWord(afResult[0], fVal);
3667
0
                *reinterpret_cast<double *>(pDstLocation) = fVal;
3668
0
            }
3669
0
            else
3670
0
            {
3671
0
                GByte abyTemp[2 * sizeof(double)];
3672
                // Convert first to the VRTRasterBand data type
3673
                // to get its clamping, before outputting to buffer data type
3674
0
                CopyWord<WorkingDT>(afResult, eWrkDataType, abyTemp,
3675
0
                                    eVRTBandDataType);
3676
0
                GDALCopyWords(abyTemp, eVRTBandDataType, 0, pDstLocation,
3677
0
                              eBufType, 0, 1);
3678
0
            }
3679
0
        }
3680
0
    }
3681
3682
0
    return CE_None;
3683
0
}
Unexecuted instantiation: CPLErr VRTComplexSource::RasterIOInternal<float>(GDALRasterBand*, GDALDataType, int, int, int, int, void*, int, int, GDALDataType, long long, long long, GDALRasterIOExtraArg*, GDALDataType, VRTSource::WorkingState&)
Unexecuted instantiation: CPLErr VRTComplexSource::RasterIOInternal<double>(GDALRasterBand*, GDALDataType, int, int, int, int, void*, int, int, GDALDataType, long long, long long, GDALRasterIOExtraArg*, GDALDataType, VRTSource::WorkingState&)
3684
3685
// Explicitly instantiate template method, as it is used in another file.
3686
template CPLErr VRTComplexSource::RasterIOInternal<float>(
3687
    GDALRasterBand *poSourceBand, GDALDataType eVRTBandDataType, int nReqXOff,
3688
    int nReqYOff, int nReqXSize, int nReqYSize, void *pData, int nOutXSize,
3689
    int nOutYSize, GDALDataType eBufType, GSpacing nPixelSpace,
3690
    GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg,
3691
    GDALDataType eWrkDataType, WorkingState &oWorkingState);
3692
3693
/************************************************************************/
3694
/*                        AreValuesUnchanged()                          */
3695
/************************************************************************/
3696
3697
bool VRTComplexSource::AreValuesUnchanged() const
3698
0
{
3699
0
    return m_dfScaleOff == 0.0 && m_dfScaleRatio == 1.0 &&
3700
0
           m_adfLUTInputs.empty() && m_nColorTableComponent == 0 &&
3701
0
           (m_nProcessingFlags & PROCESSING_FLAG_SCALING_EXPONENTIAL) == 0;
3702
0
}
3703
3704
/************************************************************************/
3705
/*                             GetMinimum()                             */
3706
/************************************************************************/
3707
3708
double VRTComplexSource::GetMinimum(int nXSize, int nYSize, int *pbSuccess)
3709
0
{
3710
0
    if (AreValuesUnchanged())
3711
0
    {
3712
0
        return VRTSimpleSource::GetMinimum(nXSize, nYSize, pbSuccess);
3713
0
    }
3714
3715
0
    *pbSuccess = FALSE;
3716
0
    return 0;
3717
0
}
3718
3719
/************************************************************************/
3720
/*                             GetMaximum()                             */
3721
/************************************************************************/
3722
3723
double VRTComplexSource::GetMaximum(int nXSize, int nYSize, int *pbSuccess)
3724
0
{
3725
0
    if (AreValuesUnchanged())
3726
0
    {
3727
0
        return VRTSimpleSource::GetMaximum(nXSize, nYSize, pbSuccess);
3728
0
    }
3729
3730
0
    *pbSuccess = FALSE;
3731
0
    return 0;
3732
0
}
3733
3734
/************************************************************************/
3735
/*                            GetHistogram()                            */
3736
/************************************************************************/
3737
3738
CPLErr VRTComplexSource::GetHistogram(int nXSize, int nYSize, double dfMin,
3739
                                      double dfMax, int nBuckets,
3740
                                      GUIntBig *panHistogram,
3741
                                      int bIncludeOutOfRange, int bApproxOK,
3742
                                      GDALProgressFunc pfnProgress,
3743
                                      void *pProgressData)
3744
0
{
3745
0
    if (AreValuesUnchanged())
3746
0
    {
3747
0
        return VRTSimpleSource::GetHistogram(
3748
0
            nXSize, nYSize, dfMin, dfMax, nBuckets, panHistogram,
3749
0
            bIncludeOutOfRange, bApproxOK, pfnProgress, pProgressData);
3750
0
    }
3751
3752
0
    return CE_Failure;
3753
0
}
3754
3755
/************************************************************************/
3756
/* ==================================================================== */
3757
/*                          VRTFuncSource                               */
3758
/* ==================================================================== */
3759
/************************************************************************/
3760
3761
/************************************************************************/
3762
/*                           VRTFuncSource()                            */
3763
/************************************************************************/
3764
3765
VRTFuncSource::VRTFuncSource()
3766
0
    : pfnReadFunc(nullptr), pCBData(nullptr), eType(GDT_Byte),
3767
0
      fNoDataValue(static_cast<float>(VRT_NODATA_UNSET))
3768
0
{
3769
0
}
3770
3771
/************************************************************************/
3772
/*                           ~VRTFuncSource()                           */
3773
/************************************************************************/
3774
3775
VRTFuncSource::~VRTFuncSource()
3776
{
3777
}
3778
3779
/************************************************************************/
3780
/*                            GetType()                                 */
3781
/************************************************************************/
3782
3783
const char *VRTFuncSource::GetType() const
3784
0
{
3785
0
    static const char *TYPE = "FuncSource";
3786
0
    return TYPE;
3787
0
}
3788
3789
/************************************************************************/
3790
/*                           SerializeToXML()                           */
3791
/************************************************************************/
3792
3793
CPLXMLNode *VRTFuncSource::SerializeToXML(CPL_UNUSED const char *pszVRTPath)
3794
0
{
3795
0
    return nullptr;
3796
0
}
3797
3798
/************************************************************************/
3799
/*                              RasterIO()                              */
3800
/************************************************************************/
3801
3802
CPLErr VRTFuncSource::RasterIO(GDALDataType /*eVRTBandDataType*/, int nXOff,
3803
                               int nYOff, int nXSize, int nYSize, void *pData,
3804
                               int nBufXSize, int nBufYSize,
3805
                               GDALDataType eBufType, GSpacing nPixelSpace,
3806
                               GSpacing nLineSpace,
3807
                               GDALRasterIOExtraArg * /* psExtraArg */,
3808
                               WorkingState & /* oWorkingState */)
3809
0
{
3810
0
    if (nPixelSpace * 8 == GDALGetDataTypeSize(eBufType) &&
3811
0
        nLineSpace == nPixelSpace * nXSize && nBufXSize == nXSize &&
3812
0
        nBufYSize == nYSize && eBufType == eType)
3813
0
    {
3814
0
        return pfnReadFunc(pCBData, nXOff, nYOff, nXSize, nYSize, pData);
3815
0
    }
3816
0
    else
3817
0
    {
3818
0
        CPLError(CE_Failure, CPLE_AppDefined,
3819
0
                 "VRTFuncSource::RasterIO() - Irregular request.");
3820
0
        CPLDebug("VRT", "Irregular request: %d,%d  %d,%d, %d,%d %d,%d %d,%d",
3821
0
                 static_cast<int>(nPixelSpace) * 8,
3822
0
                 GDALGetDataTypeSize(eBufType), static_cast<int>(nLineSpace),
3823
0
                 static_cast<int>(nPixelSpace) * nXSize, nBufXSize, nXSize,
3824
0
                 nBufYSize, nYSize, static_cast<int>(eBufType),
3825
0
                 static_cast<int>(eType));
3826
3827
0
        return CE_Failure;
3828
0
    }
3829
0
}
3830
3831
/************************************************************************/
3832
/*                             GetMinimum()                             */
3833
/************************************************************************/
3834
3835
double VRTFuncSource::GetMinimum(int /* nXSize */, int /* nYSize */,
3836
                                 int *pbSuccess)
3837
0
{
3838
0
    *pbSuccess = FALSE;
3839
0
    return 0;
3840
0
}
3841
3842
/************************************************************************/
3843
/*                             GetMaximum()                             */
3844
/************************************************************************/
3845
3846
double VRTFuncSource::GetMaximum(int /* nXSize */, int /* nYSize */,
3847
                                 int *pbSuccess)
3848
0
{
3849
0
    *pbSuccess = FALSE;
3850
0
    return 0;
3851
0
}
3852
3853
/************************************************************************/
3854
/*                            GetHistogram()                            */
3855
/************************************************************************/
3856
3857
CPLErr VRTFuncSource::GetHistogram(
3858
    int /* nXSize */, int /* nYSize */, double /* dfMin */, double /* dfMax */,
3859
    int /* nBuckets */, GUIntBig * /* panHistogram */,
3860
    int /* bIncludeOutOfRange */, int /* bApproxOK */,
3861
    GDALProgressFunc /* pfnProgress */, void * /* pProgressData */)
3862
0
{
3863
0
    return CE_Failure;
3864
0
}
3865
3866
/************************************************************************/
3867
/*                        VRTParseCoreSources()                         */
3868
/************************************************************************/
3869
3870
VRTSource *VRTParseCoreSources(const CPLXMLNode *psChild,
3871
                               const char *pszVRTPath,
3872
                               VRTMapSharedResources &oMapSharedSources)
3873
3874
0
{
3875
0
    VRTSource *poSource = nullptr;
3876
3877
0
    if (EQUAL(psChild->pszValue, VRTAveragedSource::GetTypeStatic()) ||
3878
0
        (EQUAL(psChild->pszValue, VRTSimpleSource::GetTypeStatic()) &&
3879
0
         STARTS_WITH_CI(CPLGetXMLValue(psChild, "Resampling", "Nearest"),
3880
0
                        "Aver")))
3881
0
    {
3882
0
        poSource = new VRTAveragedSource();
3883
0
    }
3884
0
    else if (EQUAL(psChild->pszValue, VRTSimpleSource::GetTypeStatic()))
3885
0
    {
3886
0
        poSource = new VRTSimpleSource();
3887
0
    }
3888
0
    else if (EQUAL(psChild->pszValue, VRTComplexSource::GetTypeStatic()))
3889
0
    {
3890
0
        poSource = new VRTComplexSource();
3891
0
    }
3892
0
    else if (EQUAL(psChild->pszValue, VRTNoDataFromMaskSource::GetTypeStatic()))
3893
0
    {
3894
0
        poSource = new VRTNoDataFromMaskSource();
3895
0
    }
3896
0
    else
3897
0
    {
3898
0
        CPLError(CE_Failure, CPLE_AppDefined,
3899
0
                 "VRTParseCoreSources() - Unknown source : %s",
3900
0
                 psChild->pszValue);
3901
0
        return nullptr;
3902
0
    }
3903
3904
0
    if (poSource->XMLInit(psChild, pszVRTPath, oMapSharedSources) == CE_None)
3905
0
        return poSource;
3906
3907
0
    delete poSource;
3908
0
    return nullptr;
3909
0
}
3910
3911
/*! @endcond */