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