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