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