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