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