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 247479 : VRTSource::~VRTSource()
53 : {
54 247479 : }
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 494079 : VRTSimpleSource::~VRTSimpleSource()
137 :
138 : {
139 247427 : if (m_bDropRefOnSrcBand)
140 : {
141 247013 : GDALDataset *poDS = GetSourceDataset();
142 247013 : if (poDS)
143 244394 : poDS->ReleaseRef();
144 : }
145 494079 : }
146 :
147 : /************************************************************************/
148 : /* GetSourceDataset() */
149 : /************************************************************************/
150 :
151 247118 : GDALDataset *VRTSimpleSource::GetSourceDataset() const
152 : {
153 247118 : GDALDataset *poDS = nullptr;
154 247118 : if (m_poMaskBandMainBand)
155 72 : poDS = m_poMaskBandMainBand->GetDataset();
156 247046 : else if (m_poRasterBand)
157 244423 : poDS = m_poRasterBand->GetDataset();
158 247118 : return poDS;
159 : }
160 :
161 : /************************************************************************/
162 : /* GetTypeStatic() */
163 : /************************************************************************/
164 :
165 251562 : const char *VRTSimpleSource::GetTypeStatic()
166 : {
167 : static const char *TYPE = "SimpleSource";
168 251562 : return TYPE;
169 : }
170 :
171 : /************************************************************************/
172 : /* GetType() */
173 : /************************************************************************/
174 :
175 22967 : const char *VRTSimpleSource::GetType() const
176 : {
177 22967 : return GetTypeStatic();
178 : }
179 :
180 : /************************************************************************/
181 : /* FlushCache() */
182 : /************************************************************************/
183 :
184 313198 : CPLErr VRTSimpleSource::FlushCache(bool bAtClosing)
185 :
186 : {
187 313198 : if (m_poMaskBandMainBand != nullptr)
188 : {
189 18 : return m_poMaskBandMainBand->FlushCache(bAtClosing);
190 : }
191 313180 : else if (m_poRasterBand != nullptr)
192 : {
193 310524 : return m_poRasterBand->FlushCache(bAtClosing);
194 : }
195 2656 : return CE_None;
196 : }
197 :
198 : /************************************************************************/
199 : /* UnsetPreservedRelativeFilenames() */
200 : /************************************************************************/
201 :
202 106 : void VRTSimpleSource::UnsetPreservedRelativeFilenames()
203 : {
204 290 : if (m_bRelativeToVRTOri &&
205 184 : !STARTS_WITH(m_osSourceFileNameOri.c_str(), "http://") &&
206 78 : !STARTS_WITH(m_osSourceFileNameOri.c_str(), "https://"))
207 : {
208 78 : m_bRelativeToVRTOri = -1;
209 78 : m_osSourceFileNameOri = "";
210 : }
211 106 : }
212 :
213 : /************************************************************************/
214 : /* SetSrcBand() */
215 : /************************************************************************/
216 :
217 463 : void VRTSimpleSource::SetSrcBand(const char *pszFilename, int nBand)
218 :
219 : {
220 463 : m_nBand = nBand;
221 463 : m_osSrcDSName = pszFilename;
222 463 : }
223 :
224 : /************************************************************************/
225 : /* SetSrcBand() */
226 : /************************************************************************/
227 :
228 142313 : void VRTSimpleSource::SetSrcBand(GDALRasterBand *poNewSrcBand)
229 :
230 : {
231 142313 : m_poRasterBand = poNewSrcBand;
232 142313 : m_nBand = m_poRasterBand->GetBand();
233 142313 : auto poDS = poNewSrcBand->GetDataset();
234 142313 : if (poDS != nullptr)
235 : {
236 142313 : m_osSrcDSName = poDS->GetDescription();
237 142313 : m_aosOpenOptions = CSLDuplicate(poDS->GetOpenOptions());
238 142313 : m_aosOpenOptionsOri = m_aosOpenOptions;
239 : }
240 142313 : }
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 2607480 : static double RoundIfCloseToInt(double dfValue)
282 : {
283 2607480 : double dfClosestInt = floor(dfValue + 0.5);
284 2607480 : return (fabs(dfValue - dfClosestInt) < 1e-3) ? dfClosestInt : dfValue;
285 : }
286 :
287 : /************************************************************************/
288 : /* SetSrcWindow() */
289 : /************************************************************************/
290 :
291 145953 : void VRTSimpleSource::SetSrcWindow(double dfNewXOff, double dfNewYOff,
292 : double dfNewXSize, double dfNewYSize)
293 :
294 : {
295 145953 : m_dfSrcXOff = RoundIfCloseToInt(dfNewXOff);
296 145953 : m_dfSrcYOff = RoundIfCloseToInt(dfNewYOff);
297 145953 : m_dfSrcXSize = RoundIfCloseToInt(dfNewXSize);
298 145953 : m_dfSrcYSize = RoundIfCloseToInt(dfNewYSize);
299 145953 : }
300 :
301 : /************************************************************************/
302 : /* SetDstWindow() */
303 : /************************************************************************/
304 :
305 145952 : void VRTSimpleSource::SetDstWindow(double dfNewXOff, double dfNewYOff,
306 : double dfNewXSize, double dfNewYSize)
307 :
308 : {
309 145952 : m_dfDstXOff = RoundIfCloseToInt(dfNewXOff);
310 145952 : m_dfDstYOff = RoundIfCloseToInt(dfNewYOff);
311 145952 : m_dfDstXSize = RoundIfCloseToInt(dfNewXSize);
312 145952 : m_dfDstYSize = RoundIfCloseToInt(dfNewYSize);
313 145952 : }
314 :
315 : /************************************************************************/
316 : /* GetDstWindow() */
317 : /************************************************************************/
318 :
319 539 : void VRTSimpleSource::GetDstWindow(double &dfDstXOff, double &dfDstYOff,
320 : double &dfDstXSize, double &dfDstYSize) const
321 : {
322 539 : dfDstXOff = m_dfDstXOff;
323 539 : dfDstYOff = m_dfDstYOff;
324 539 : dfDstXSize = m_dfDstXSize;
325 539 : dfDstYSize = m_dfDstYSize;
326 539 : }
327 :
328 : /************************************************************************/
329 : /* DstWindowIntersects() */
330 : /************************************************************************/
331 :
332 1094 : bool VRTSimpleSource::DstWindowIntersects(double dfXOff, double dfYOff,
333 : double dfXSize, double dfYSize) const
334 : {
335 2179 : return IsDstWinSet() && m_dfDstXOff + m_dfDstXSize > dfXOff &&
336 1085 : m_dfDstYOff + m_dfDstYSize > dfYOff &&
337 2179 : m_dfDstXOff < dfXOff + dfXSize && m_dfDstYOff < dfYOff + dfYSize;
338 : }
339 :
340 : /************************************************************************/
341 : /* IsSlowSource() */
342 : /************************************************************************/
343 :
344 2548 : static bool IsSlowSource(const char *pszSrcName)
345 : {
346 5096 : return strstr(pszSrcName, "/vsicurl/http") != nullptr ||
347 5096 : strstr(pszSrcName, "/vsicurl/ftp") != nullptr ||
348 2548 : (strstr(pszSrcName, "/vsicurl?") != nullptr &&
349 2548 : strstr(pszSrcName, "&url=http") != nullptr);
350 : }
351 :
352 : /************************************************************************/
353 : /* AddSourceFilenameNode() */
354 : /************************************************************************/
355 :
356 2614 : void VRTSimpleSource::AddSourceFilenameNode(const char *pszVRTPath,
357 : CPLXMLNode *psSrc)
358 : {
359 :
360 : VSIStatBufL sStat;
361 2614 : int bRelativeToVRT = FALSE; // TODO(schwehr): Make this a bool?
362 5228 : std::string osSourceFilename;
363 :
364 2614 : if (m_bRelativeToVRTOri >= 0)
365 : {
366 66 : osSourceFilename = m_osSourceFileNameOri;
367 66 : bRelativeToVRT = m_bRelativeToVRTOri;
368 : }
369 2548 : 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 2548 : 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 4950 : std::string osVRTFilename = pszVRTPath;
455 4950 : std::string osSourceDataset = m_osSrcDSName;
456 2475 : char *pszCurDir = CPLGetCurrentDir();
457 2475 : if (CPLIsFilenameRelative(osSourceDataset.c_str()) &&
458 2475 : !CPLIsFilenameRelative(osVRTFilename.c_str()) &&
459 : pszCurDir != nullptr)
460 : {
461 184 : osSourceDataset = CPLFormFilenameSafe(
462 92 : pszCurDir, osSourceDataset.c_str(), nullptr);
463 : }
464 2383 : else if (!CPLIsFilenameRelative(osSourceDataset.c_str()) &&
465 2383 : CPLIsFilenameRelative(osVRTFilename.c_str()) &&
466 : pszCurDir != nullptr)
467 : {
468 : osVRTFilename =
469 84 : CPLFormFilenameSafe(pszCurDir, osVRTFilename.c_str(), nullptr);
470 : }
471 2475 : CPLFree(pszCurDir);
472 : osSourceFilename = CPLExtractRelativePath(
473 2475 : osVRTFilename.c_str(), osSourceDataset.c_str(), &bRelativeToVRT);
474 : }
475 :
476 2614 : CPLSetXMLValue(psSrc, "SourceFilename", osSourceFilename.c_str());
477 :
478 2614 : CPLCreateXMLNode(CPLCreateXMLNode(CPLGetXMLNode(psSrc, "SourceFilename"),
479 : CXT_Attribute, "relativeToVRT"),
480 2614 : 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 2614 : const char *pszShared = CPLGetConfigOption("VRT_SHARED_SOURCE", nullptr);
485 2614 : 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 2614 : }
494 :
495 : /************************************************************************/
496 : /* SerializeToXML() */
497 : /************************************************************************/
498 :
499 2655 : CPLXMLNode *VRTSimpleSource::SerializeToXML(const char *pszVRTPath)
500 :
501 : {
502 : CPLXMLNode *const psSrc =
503 2655 : CPLCreateXMLNode(nullptr, CXT_Element, GetTypeStatic());
504 :
505 2655 : if (!m_osResampling.empty())
506 : {
507 24 : CPLCreateXMLNode(CPLCreateXMLNode(psSrc, CXT_Attribute, "resampling"),
508 : CXT_Text, m_osResampling.c_str());
509 : }
510 :
511 2655 : if (!m_osName.empty())
512 : {
513 73 : CPLAddXMLAttributeAndValue(psSrc, "name", m_osName.c_str());
514 : }
515 :
516 2655 : if (m_bSrcDSNameFromVRT)
517 : {
518 1 : CPLAddXMLChild(psSrc, CPLParseXMLString(m_osSrcDSName.c_str()));
519 : }
520 : else
521 : {
522 2654 : bool bDone = false;
523 2654 : if (m_osSrcDSName.empty() && m_poRasterBand)
524 : {
525 100 : auto poSrcDS = m_poRasterBand->GetDataset();
526 100 : if (poSrcDS)
527 : {
528 100 : VRTDataset *poSrcVRTDS = nullptr;
529 : // For GDALComputedDataset
530 100 : void *pHandle = poSrcDS->GetInternalHandle("VRT_DATASET");
531 100 : if (pHandle && poSrcDS->GetInternalHandle(nullptr) == nullptr)
532 : {
533 40 : poSrcVRTDS = static_cast<VRTDataset *>(pHandle);
534 : }
535 : else
536 : {
537 60 : poSrcVRTDS = dynamic_cast<VRTDataset *>(poSrcDS);
538 : }
539 100 : if (poSrcVRTDS)
540 : {
541 40 : poSrcVRTDS->UnsetPreservedRelativeFilenames();
542 40 : CPLAddXMLChild(psSrc,
543 40 : poSrcVRTDS->SerializeToXML(pszVRTPath));
544 40 : bDone = true;
545 : }
546 : }
547 : }
548 2654 : if (!bDone)
549 : {
550 2614 : AddSourceFilenameNode(pszVRTPath, psSrc);
551 : }
552 : }
553 :
554 2655 : GDALSerializeOpenOptionsToXML(psSrc, m_aosOpenOptionsOri.List());
555 :
556 2655 : if (m_bGetMaskBand)
557 19 : CPLSetXMLValue(psSrc, "SourceBand", CPLSPrintf("mask,%d", m_nBand));
558 : else
559 2636 : 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 2655 : 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 2536 : CPLSetXMLValue(psSrc, "SourceProperties.#RasterXSize",
569 2536 : CPLSPrintf("%d", m_poRasterBand->GetXSize()));
570 2536 : CPLSetXMLValue(psSrc, "SourceProperties.#RasterYSize",
571 2536 : CPLSPrintf("%d", m_poRasterBand->GetYSize()));
572 2536 : CPLSetXMLValue(
573 : psSrc, "SourceProperties.#DataType",
574 2536 : GDALGetDataTypeName(m_poRasterBand->GetRasterDataType()));
575 :
576 2536 : int nBlockXSize = 0;
577 2536 : int nBlockYSize = 0;
578 2536 : m_poRasterBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
579 :
580 2536 : CPLSetXMLValue(psSrc, "SourceProperties.#BlockXSize",
581 : CPLSPrintf("%d", nBlockXSize));
582 2536 : CPLSetXMLValue(psSrc, "SourceProperties.#BlockYSize",
583 : CPLSPrintf("%d", nBlockYSize));
584 : }
585 :
586 2655 : if (IsSrcWinSet())
587 : {
588 2641 : CPLSetXMLValue(psSrc, "SrcRect.#xOff",
589 : CPLSPrintf("%.15g", m_dfSrcXOff));
590 2641 : CPLSetXMLValue(psSrc, "SrcRect.#yOff",
591 : CPLSPrintf("%.15g", m_dfSrcYOff));
592 2641 : CPLSetXMLValue(psSrc, "SrcRect.#xSize",
593 : CPLSPrintf("%.15g", m_dfSrcXSize));
594 2641 : CPLSetXMLValue(psSrc, "SrcRect.#ySize",
595 : CPLSPrintf("%.15g", m_dfSrcYSize));
596 : }
597 :
598 2655 : if (IsDstWinSet())
599 : {
600 2641 : CPLSetXMLValue(psSrc, "DstRect.#xOff",
601 : CPLSPrintf("%.15g", m_dfDstXOff));
602 2641 : CPLSetXMLValue(psSrc, "DstRect.#yOff",
603 : CPLSPrintf("%.15g", m_dfDstYOff));
604 2641 : CPLSetXMLValue(psSrc, "DstRect.#xSize",
605 : CPLSPrintf("%.15g", m_dfDstXSize));
606 2641 : CPLSetXMLValue(psSrc, "DstRect.#ySize",
607 : CPLSPrintf("%.15g", m_dfDstYSize));
608 : }
609 :
610 2655 : return psSrc;
611 : }
612 :
613 : /************************************************************************/
614 : /* XMLInit() */
615 : /************************************************************************/
616 :
617 104140 : CPLErr VRTSimpleSource::XMLInit(const CPLXMLNode *psSrc, const char *pszVRTPath,
618 : VRTMapSharedResources &oMapSharedSources)
619 :
620 : {
621 104140 : m_poMapSharedSources = &oMapSharedSources;
622 :
623 104140 : m_osResampling = CPLGetXMLValue(psSrc, "resampling", "");
624 104140 : m_osName = CPLGetXMLValue(psSrc, "name", "");
625 :
626 : /* -------------------------------------------------------------------- */
627 : /* Prepare filename. */
628 : /* -------------------------------------------------------------------- */
629 : const CPLXMLNode *psSourceFileNameNode =
630 104140 : CPLGetXMLNode(psSrc, "SourceFilename");
631 104140 : const CPLXMLNode *psSourceVRTDataset = CPLGetXMLNode(psSrc, "VRTDataset");
632 : const char *pszFilename =
633 104140 : psSourceFileNameNode ? CPLGetXMLValue(psSourceFileNameNode, nullptr, "")
634 104140 : : "";
635 :
636 104140 : 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 104139 : m_osSourceFileNameOri = pszFilename;
647 104139 : if (pszFilename[0])
648 : {
649 104130 : m_bRelativeToVRTOri =
650 104130 : atoi(CPLGetXMLValue(psSourceFileNameNode, "relativetoVRT", "0"));
651 : const char *pszShared =
652 104130 : CPLGetXMLValue(psSourceFileNameNode, "shared", nullptr);
653 104130 : if (pszShared == nullptr)
654 : {
655 104128 : pszShared = CPLGetConfigOption("VRT_SHARED_SOURCE", nullptr);
656 : }
657 104130 : if (pszShared != nullptr)
658 : {
659 8 : m_nExplicitSharedStatus = CPLTestBool(pszShared);
660 : }
661 :
662 208260 : m_osSrcDSName = GDALDataset::BuildFilename(
663 208260 : pszFilename, pszVRTPath, CPL_TO_BOOL(m_bRelativeToVRTOri));
664 : }
665 9 : else if (psSourceVRTDataset)
666 : {
667 : CPLXMLNode sNode;
668 9 : sNode.eType = psSourceVRTDataset->eType;
669 9 : sNode.pszValue = psSourceVRTDataset->pszValue;
670 9 : sNode.psNext = nullptr;
671 9 : sNode.psChild = psSourceVRTDataset->psChild;
672 9 : char *pszXML = CPLSerializeXMLTree(&sNode);
673 9 : if (pszXML)
674 : {
675 9 : m_bSrcDSNameFromVRT = true;
676 9 : m_osSrcDSName = pszXML;
677 9 : CPLFree(pszXML);
678 : }
679 : }
680 :
681 104139 : const char *pszSourceBand = CPLGetXMLValue(psSrc, "SourceBand", "1");
682 104139 : m_bGetMaskBand = false;
683 104139 : 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 104116 : m_nBand = atoi(pszSourceBand);
694 : }
695 104139 : 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 104139 : m_aosOpenOptions = GDALDeserializeOpenOptionsFromXML(psSrc);
703 104139 : m_aosOpenOptionsOri = m_aosOpenOptions;
704 104139 : if (strstr(m_osSrcDSName.c_str(), "<VRTDataset") != nullptr)
705 11 : m_aosOpenOptions.SetNameValue("ROOT_PATH", pszVRTPath);
706 :
707 104139 : return ParseSrcRectAndDstRect(psSrc);
708 : }
709 :
710 : /************************************************************************/
711 : /* ParseSrcRectAndDstRect() */
712 : /************************************************************************/
713 :
714 104161 : CPLErr VRTSimpleSource::ParseSrcRectAndDstRect(const CPLXMLNode *psSrc)
715 : {
716 25276 : const auto GetAttrValue = [](const CPLXMLNode *psNode,
717 : const char *pszAttrName, double dfDefaultVal)
718 : {
719 25276 : if (const char *pszVal = CPLGetXMLValue(psNode, pszAttrName, nullptr))
720 25276 : return CPLAtof(pszVal);
721 : else
722 0 : return dfDefaultVal;
723 : };
724 :
725 : /* -------------------------------------------------------------------- */
726 : /* Set characteristics. */
727 : /* -------------------------------------------------------------------- */
728 104161 : const CPLXMLNode *const psSrcRect = CPLGetXMLNode(psSrc, "SrcRect");
729 104161 : if (psSrcRect)
730 : {
731 3160 : double xOff = GetAttrValue(psSrcRect, "xOff", UNINIT_WINDOW);
732 3160 : double yOff = GetAttrValue(psSrcRect, "yOff", UNINIT_WINDOW);
733 3160 : double xSize = GetAttrValue(psSrcRect, "xSize", UNINIT_WINDOW);
734 3160 : double ySize = GetAttrValue(psSrcRect, "ySize", UNINIT_WINDOW);
735 : // Test written that way to catch NaN values
736 3160 : if (!(xOff >= INT_MIN && xOff <= INT_MAX) ||
737 3160 : !(yOff >= INT_MIN && yOff <= INT_MAX) ||
738 3160 : !(xSize > 0 || xSize == UNINIT_WINDOW) || xSize > INT_MAX ||
739 3159 : !(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 3159 : SetSrcWindow(xOff, yOff, xSize, ySize);
745 : }
746 : else
747 : {
748 101001 : m_dfSrcXOff = UNINIT_WINDOW;
749 101001 : m_dfSrcYOff = UNINIT_WINDOW;
750 101001 : m_dfSrcXSize = UNINIT_WINDOW;
751 101001 : m_dfSrcYSize = UNINIT_WINDOW;
752 : }
753 :
754 104160 : const CPLXMLNode *const psDstRect = CPLGetXMLNode(psSrc, "DstRect");
755 104160 : if (psDstRect)
756 : {
757 3159 : double xOff = GetAttrValue(psDstRect, "xOff", UNINIT_WINDOW);
758 : ;
759 3159 : double yOff = GetAttrValue(psDstRect, "yOff", UNINIT_WINDOW);
760 3159 : double xSize = GetAttrValue(psDstRect, "xSize", UNINIT_WINDOW);
761 : ;
762 3159 : double ySize = GetAttrValue(psDstRect, "ySize", UNINIT_WINDOW);
763 : // Test written that way to catch NaN values
764 3159 : if (!(xOff >= INT_MIN && xOff <= INT_MAX) ||
765 3159 : !(yOff >= INT_MIN && yOff <= INT_MAX) ||
766 3159 : !(xSize > 0 || xSize == UNINIT_WINDOW) || xSize > INT_MAX ||
767 3159 : !(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 3158 : SetDstWindow(xOff, yOff, xSize, ySize);
773 : }
774 : else
775 : {
776 101001 : m_dfDstXOff = UNINIT_WINDOW;
777 101001 : m_dfDstYOff = UNINIT_WINDOW;
778 101001 : m_dfDstXSize = UNINIT_WINDOW;
779 101001 : m_dfDstYSize = UNINIT_WINDOW;
780 : }
781 :
782 104159 : return CE_None;
783 : }
784 :
785 : /************************************************************************/
786 : /* GetFileList() */
787 : /************************************************************************/
788 :
789 113 : void VRTSimpleSource::GetFileList(char ***ppapszFileList, int *pnSize,
790 : int *pnMaxSize, CPLHashSet *hSetFiles)
791 : {
792 113 : if (!m_osSrcDSName.empty() && !m_bSrcDSNameFromVRT)
793 : {
794 113 : const char *pszFilename = m_osSrcDSName.c_str();
795 :
796 : /* --------------------------------------------------------------------
797 : */
798 : /* Is it already in the list ? */
799 : /* --------------------------------------------------------------------
800 : */
801 113 : if (CPLHashSetLookup(hSetFiles, pszFilename) != nullptr)
802 15 : return;
803 :
804 : /* --------------------------------------------------------------------
805 : */
806 : /* Grow array if necessary */
807 : /* --------------------------------------------------------------------
808 : */
809 98 : if (*pnSize + 1 >= *pnMaxSize)
810 : {
811 67 : *pnMaxSize = std::max(*pnSize + 2, 2 + 2 * (*pnMaxSize));
812 67 : *ppapszFileList = static_cast<char **>(
813 67 : CPLRealloc(*ppapszFileList, sizeof(char *) * (*pnMaxSize)));
814 : }
815 :
816 : /* --------------------------------------------------------------------
817 : */
818 : /* Add the string to the list */
819 : /* --------------------------------------------------------------------
820 : */
821 98 : (*ppapszFileList)[*pnSize] = CPLStrdup(pszFilename);
822 98 : (*ppapszFileList)[(*pnSize + 1)] = nullptr;
823 98 : CPLHashSetInsert(hSetFiles, (*ppapszFileList)[*pnSize]);
824 :
825 98 : (*pnSize)++;
826 : }
827 : }
828 :
829 : /************************************************************************/
830 : /* OpenSource() */
831 : /************************************************************************/
832 :
833 102251 : void VRTSimpleSource::OpenSource() const
834 : {
835 102251 : CPLAssert(m_poRasterBand == nullptr);
836 :
837 : /* ----------------------------------------------------------------- */
838 : /* Create a proxy dataset */
839 : /* ----------------------------------------------------------------- */
840 102251 : GDALProxyPoolDataset *proxyDS = nullptr;
841 102251 : std::string osKeyMapSharedSources;
842 102251 : if (m_poMapSharedSources)
843 : {
844 102203 : osKeyMapSharedSources = m_osSrcDSName;
845 102217 : for (int i = 0; i < m_aosOpenOptions.size(); ++i)
846 : {
847 14 : osKeyMapSharedSources += "||";
848 14 : osKeyMapSharedSources += m_aosOpenOptions[i];
849 : }
850 :
851 102203 : proxyDS = cpl::down_cast<GDALProxyPoolDataset *>(
852 102203 : m_poMapSharedSources->Get(osKeyMapSharedSources));
853 : }
854 :
855 102251 : if (proxyDS == nullptr)
856 : {
857 3275 : int bShared = true;
858 3275 : if (m_nExplicitSharedStatus != -1)
859 8 : bShared = m_nExplicitSharedStatus;
860 :
861 3275 : const CPLString osUniqueHandle(CPLSPrintf("%p", m_poMapSharedSources));
862 3274 : proxyDS = GDALProxyPoolDataset::Create(
863 : m_osSrcDSName, m_aosOpenOptions.List(), GA_ReadOnly, bShared,
864 : osUniqueHandle.c_str());
865 3275 : if (proxyDS == nullptr)
866 252 : return;
867 : }
868 : else
869 : {
870 98976 : proxyDS->Reference();
871 : }
872 :
873 101999 : 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 101999 : m_poRasterBand = proxyDS->GetRasterBand(m_nBand);
886 101999 : if (m_poRasterBand == nullptr || !ValidateOpenedBand(m_poRasterBand))
887 : {
888 2 : proxyDS->ReleaseRef();
889 2 : return;
890 : }
891 :
892 101997 : 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 101997 : if (m_poMapSharedSources)
904 : {
905 101949 : m_poMapSharedSources->Insert(osKeyMapSharedSources, proxyDS);
906 : }
907 : }
908 :
909 : /************************************************************************/
910 : /* GetRasterBand() */
911 : /************************************************************************/
912 :
913 936317 : GDALRasterBand *VRTSimpleSource::GetRasterBand() const
914 : {
915 936317 : if (m_poRasterBand == nullptr)
916 102102 : OpenSource();
917 936317 : return m_poRasterBand;
918 : }
919 :
920 : /************************************************************************/
921 : /* GetMaskBandMainBand() */
922 : /************************************************************************/
923 :
924 5638 : GDALRasterBand *VRTSimpleSource::GetMaskBandMainBand()
925 : {
926 5638 : if (m_poRasterBand == nullptr)
927 149 : OpenSource();
928 5638 : return m_poMaskBandMainBand;
929 : }
930 :
931 : /************************************************************************/
932 : /* IsSameExceptBandNumber() */
933 : /************************************************************************/
934 :
935 2889 : bool VRTSimpleSource::IsSameExceptBandNumber(
936 : const VRTSimpleSource *poOtherSource) const
937 : {
938 5778 : return m_dfSrcXOff == poOtherSource->m_dfSrcXOff &&
939 2889 : m_dfSrcYOff == poOtherSource->m_dfSrcYOff &&
940 2889 : m_dfSrcXSize == poOtherSource->m_dfSrcXSize &&
941 2889 : m_dfSrcYSize == poOtherSource->m_dfSrcYSize &&
942 2889 : m_dfDstXOff == poOtherSource->m_dfDstXOff &&
943 2889 : m_dfDstYOff == poOtherSource->m_dfDstYOff &&
944 2889 : m_dfDstXSize == poOtherSource->m_dfDstXSize &&
945 8667 : m_dfDstYSize == poOtherSource->m_dfDstYSize &&
946 5778 : m_osSrcDSName == poOtherSource->m_osSrcDSName;
947 : }
948 :
949 : /************************************************************************/
950 : /* SrcToDst() */
951 : /* */
952 : /* Note: this is a no-op if the both src and dst windows are unset */
953 : /************************************************************************/
954 :
955 10320 : void VRTSimpleSource::SrcToDst(double dfX, double dfY, double &dfXOut,
956 : double &dfYOut) const
957 :
958 : {
959 10320 : dfXOut = ((dfX - m_dfSrcXOff) / m_dfSrcXSize) * m_dfDstXSize + m_dfDstXOff;
960 10320 : dfYOut = ((dfY - m_dfSrcYOff) / m_dfSrcYSize) * m_dfDstYSize + m_dfDstYOff;
961 10320 : }
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 381280 : 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 381280 : bErrorOut = false;
1007 :
1008 381280 : if (m_dfSrcXSize == 0.0 || m_dfSrcYSize == 0.0 || m_dfDstXSize == 0.0 ||
1009 381280 : m_dfDstYSize == 0.0)
1010 : {
1011 0 : return FALSE;
1012 : }
1013 :
1014 381280 : const bool bDstWinSet = IsDstWinSet();
1015 :
1016 : #ifdef DEBUG
1017 381280 : const bool bSrcWinSet = IsSrcWinSet();
1018 :
1019 381280 : 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 381280 : if (bDstWinSet)
1030 : {
1031 278193 : if (dfXOff >= m_dfDstXOff + m_dfDstXSize ||
1032 267765 : dfYOff >= m_dfDstYOff + m_dfDstYSize ||
1033 266333 : 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 360073 : *pnOutXOff = 0;
1041 360073 : *pnOutYOff = 0;
1042 360073 : *pnOutXSize = nBufXSize;
1043 360073 : *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 360073 : bool bModifiedX = false;
1051 360073 : bool bModifiedY = false;
1052 360073 : double dfRXOff = dfXOff;
1053 360073 : double dfRYOff = dfYOff;
1054 360073 : double dfRXSize = dfXSize;
1055 360073 : double dfRYSize = dfYSize;
1056 :
1057 360073 : if (bDstWinSet)
1058 : {
1059 256986 : if (dfRXOff < m_dfDstXOff)
1060 : {
1061 1907 : dfRXSize = dfRXSize + dfRXOff - m_dfDstXOff;
1062 1907 : dfRXOff = m_dfDstXOff;
1063 1907 : bModifiedX = true;
1064 : }
1065 :
1066 256986 : if (dfRYOff < m_dfDstYOff)
1067 : {
1068 172 : dfRYSize = dfRYSize + dfRYOff - m_dfDstYOff;
1069 172 : dfRYOff = m_dfDstYOff;
1070 172 : bModifiedY = true;
1071 : }
1072 :
1073 256986 : if (dfRXOff + dfRXSize > m_dfDstXOff + m_dfDstXSize)
1074 : {
1075 3580 : dfRXSize = m_dfDstXOff + m_dfDstXSize - dfRXOff;
1076 3580 : bModifiedX = true;
1077 : }
1078 :
1079 256986 : 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 360073 : const double dfScaleX = m_dfSrcXSize / m_dfDstXSize;
1091 360073 : const double dfScaleY = m_dfSrcYSize / m_dfDstYSize;
1092 :
1093 360073 : *pdfReqXOff = (dfRXOff - m_dfDstXOff) * dfScaleX + m_dfSrcXOff;
1094 360073 : *pdfReqYOff = (dfRYOff - m_dfDstYOff) * dfScaleY + m_dfSrcYOff;
1095 360073 : *pdfReqXSize = dfRXSize * dfScaleX;
1096 360073 : *pdfReqYSize = dfRYSize * dfScaleY;
1097 :
1098 720146 : if (!std::isfinite(*pdfReqXOff) || !std::isfinite(*pdfReqYOff) ||
1099 360073 : !std::isfinite(*pdfReqXSize) || !std::isfinite(*pdfReqYSize) ||
1100 1080220 : *pdfReqXOff > INT_MAX || *pdfReqYOff > INT_MAX || *pdfReqXSize < 0 ||
1101 360073 : *pdfReqYSize < 0)
1102 : {
1103 0 : return FALSE;
1104 : }
1105 :
1106 : /* -------------------------------------------------------------------- */
1107 : /* Clamp within the bounds of the available source data. */
1108 : /* -------------------------------------------------------------------- */
1109 360073 : if (*pdfReqXOff < 0)
1110 : {
1111 9 : *pdfReqXSize += *pdfReqXOff;
1112 9 : *pdfReqXOff = 0;
1113 9 : bModifiedX = true;
1114 : }
1115 360073 : if (*pdfReqYOff < 0)
1116 : {
1117 9 : *pdfReqYSize += *pdfReqYOff;
1118 9 : *pdfReqYOff = 0;
1119 9 : bModifiedY = true;
1120 : }
1121 :
1122 360073 : constexpr double EPSILON = 1e-10;
1123 360073 : if (eResampleAlg == GRIORA_NearestNeighbour &&
1124 357493 : (std::fabs(m_dfSrcXOff - std::round(m_dfSrcXOff)) > EPSILON ||
1125 357429 : std::fabs(m_dfSrcYOff - std::round(m_dfSrcYOff)) > EPSILON ||
1126 357427 : std::fabs(m_dfDstXOff - std::round(m_dfDstXOff)) > EPSILON ||
1127 357421 : 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 359994 : *pnReqXOff = static_cast<int>(*pdfReqXOff);
1139 359994 : *pnReqYOff = static_cast<int>(*pdfReqYOff);
1140 : }
1141 :
1142 360073 : constexpr double EPS = 1e-3;
1143 360073 : constexpr double ONE_MINUS_EPS = 1.0 - EPS;
1144 360073 : if (*pdfReqXOff - *pnReqXOff > ONE_MINUS_EPS)
1145 : {
1146 0 : (*pnReqXOff)++;
1147 0 : *pdfReqXOff = *pnReqXOff;
1148 : }
1149 360073 : if (*pdfReqYOff - *pnReqYOff > ONE_MINUS_EPS)
1150 : {
1151 16 : (*pnReqYOff)++;
1152 16 : *pdfReqYOff = *pnReqYOff;
1153 : }
1154 :
1155 360073 : if (*pdfReqXSize > INT_MAX)
1156 0 : *pnReqXSize = INT_MAX;
1157 : else
1158 360073 : *pnReqXSize = static_cast<int>(floor(*pdfReqXSize + 0.5));
1159 :
1160 360073 : if (*pdfReqYSize > INT_MAX)
1161 0 : *pnReqYSize = INT_MAX;
1162 : else
1163 360073 : *pnReqYSize = static_cast<int>(floor(*pdfReqYSize + 0.5));
1164 :
1165 : /* -------------------------------------------------------------------- */
1166 : /* Clamp within the bounds of the available source data. */
1167 : /* -------------------------------------------------------------------- */
1168 :
1169 360073 : if (*pnReqXSize == 0)
1170 681 : *pnReqXSize = 1;
1171 360073 : if (*pnReqYSize == 0)
1172 22684 : *pnReqYSize = 1;
1173 :
1174 360073 : auto l_band = GetRasterBand();
1175 360073 : if (!l_band)
1176 : {
1177 89 : bErrorOut = true;
1178 89 : return FALSE;
1179 : }
1180 719968 : if (*pnReqXSize > INT_MAX - *pnReqXOff ||
1181 359984 : *pnReqXOff + *pnReqXSize > l_band->GetXSize())
1182 : {
1183 253 : *pnReqXSize = l_band->GetXSize() - *pnReqXOff;
1184 253 : bModifiedX = true;
1185 : }
1186 359984 : if (*pdfReqXOff + *pdfReqXSize > l_band->GetXSize())
1187 : {
1188 251 : *pdfReqXSize = l_band->GetXSize() - *pdfReqXOff;
1189 251 : bModifiedX = true;
1190 : }
1191 :
1192 719968 : if (*pnReqYSize > INT_MAX - *pnReqYOff ||
1193 359984 : *pnReqYOff + *pnReqYSize > l_band->GetYSize())
1194 : {
1195 48 : *pnReqYSize = l_band->GetYSize() - *pnReqYOff;
1196 48 : bModifiedY = true;
1197 : }
1198 359984 : 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 719961 : if (*pnReqXOff >= l_band->GetXSize() || *pnReqYOff >= l_band->GetYSize() ||
1209 719961 : *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 359969 : 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 5160 : double dfDstULX = 0.0;
1228 5160 : double dfDstULY = 0.0;
1229 5160 : double dfDstLRX = 0.0;
1230 5160 : double dfDstLRY = 0.0;
1231 :
1232 5160 : SrcToDst(*pdfReqXOff, *pdfReqYOff, dfDstULX, dfDstULY);
1233 5160 : 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 5160 : if (bModifiedX)
1241 : {
1242 5106 : const double dfScaleWinToBufX = nBufXSize / dfXSize;
1243 :
1244 5106 : const double dfOutXOff = (dfDstULX - dfXOff) * dfScaleWinToBufX;
1245 5106 : if (dfOutXOff <= 0)
1246 3201 : *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 5106 : double dfDstDeltaX =
1255 5106 : (dfOutXOff - *pnOutXOff) / dfScaleWinToBufX;
1256 5106 : double dfSrcDeltaX = dfDstDeltaX / m_dfDstXSize * m_dfSrcXSize;
1257 5106 : *pdfReqXOff -= dfSrcDeltaX;
1258 10212 : *pdfReqXSize = std::min(*pdfReqXSize + dfSrcDeltaX,
1259 5106 : static_cast<double>(INT_MAX));
1260 : }
1261 :
1262 5106 : double dfOutRightXOff = (dfDstLRX - dfXOff) * dfScaleWinToBufX;
1263 5106 : if (dfOutRightXOff < dfOutXOff)
1264 3 : return FALSE;
1265 5106 : if (dfOutRightXOff > INT_MAX)
1266 0 : dfOutRightXOff = INT_MAX;
1267 5106 : const int nOutRightXOff =
1268 5106 : static_cast<int>(ceil(dfOutRightXOff - EPS));
1269 5106 : *pnOutXSize = nOutRightXOff - *pnOutXOff;
1270 :
1271 5106 : if (*pnOutXSize > INT_MAX - *pnOutXOff ||
1272 5106 : *pnOutXOff + *pnOutXSize > nBufXSize)
1273 0 : *pnOutXSize = nBufXSize - *pnOutXOff;
1274 :
1275 : // Apply correction on floating-point source window
1276 : {
1277 5106 : double dfDstDeltaX =
1278 5106 : (nOutRightXOff - dfOutRightXOff) / dfScaleWinToBufX;
1279 5106 : double dfSrcDeltaX = dfDstDeltaX / m_dfDstXSize * m_dfSrcXSize;
1280 10212 : *pdfReqXSize = std::min(*pdfReqXSize + dfSrcDeltaX,
1281 5106 : static_cast<double>(INT_MAX));
1282 : }
1283 : }
1284 :
1285 5160 : 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 5160 : if (*pnOutXSize < 1 || *pnOutYSize < 1)
1330 3 : return FALSE;
1331 : }
1332 :
1333 359966 : *pdfReqXOff = RoundIfCloseToInt(*pdfReqXOff);
1334 359966 : *pdfReqYOff = RoundIfCloseToInt(*pdfReqYOff);
1335 359966 : *pdfReqXSize = RoundIfCloseToInt(*pdfReqXSize);
1336 359966 : *pdfReqYSize = RoundIfCloseToInt(*pdfReqYSize);
1337 :
1338 359966 : return TRUE;
1339 : }
1340 :
1341 : /************************************************************************/
1342 : /* NeedMaxValAdjustment() */
1343 : /************************************************************************/
1344 :
1345 347635 : int VRTSimpleSource::NeedMaxValAdjustment() const
1346 : {
1347 347635 : if (!m_nMaxValue)
1348 347613 : 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 352409 : 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 352409 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1378 352409 : GDALRasterIOExtraArg *psExtraArg = &sExtraArg;
1379 :
1380 352409 : double dfXOff = nXOff;
1381 352409 : double dfYOff = nYOff;
1382 352409 : double dfXSize = nXSize;
1383 352409 : double dfYSize = nYSize;
1384 352409 : if (psExtraArgIn != nullptr && psExtraArgIn->bFloatingPointWindowValidity)
1385 : {
1386 351 : dfXOff = psExtraArgIn->dfXOff;
1387 351 : dfYOff = psExtraArgIn->dfYOff;
1388 351 : dfXSize = psExtraArgIn->dfXSize;
1389 351 : dfYSize = psExtraArgIn->dfYSize;
1390 : }
1391 :
1392 : // The window we will actually request from the source raster band.
1393 352409 : double dfReqXOff = 0.0;
1394 352409 : double dfReqYOff = 0.0;
1395 352409 : double dfReqXSize = 0.0;
1396 352409 : double dfReqYSize = 0.0;
1397 352409 : int nReqXOff = 0;
1398 352409 : int nReqYOff = 0;
1399 352409 : int nReqXSize = 0;
1400 352409 : int nReqYSize = 0;
1401 :
1402 : // The window we will actual set _within_ the pData buffer.
1403 352409 : int nOutXOff = 0;
1404 352409 : int nOutYOff = 0;
1405 352409 : int nOutXSize = 0;
1406 352409 : int nOutYSize = 0;
1407 :
1408 352409 : if (!m_osResampling.empty())
1409 : {
1410 202747 : psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
1411 : }
1412 149662 : else if (psExtraArgIn != nullptr)
1413 : {
1414 149662 : psExtraArg->eResampleAlg = psExtraArgIn->eResampleAlg;
1415 : }
1416 :
1417 352409 : bool bError = false;
1418 352409 : 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 9446 : 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 342963 : psExtraArg->bFloatingPointWindowValidity = TRUE;
1442 342963 : psExtraArg->dfXOff = dfReqXOff;
1443 342963 : psExtraArg->dfYOff = dfReqYOff;
1444 342963 : psExtraArg->dfXSize = dfReqXSize;
1445 342963 : psExtraArg->dfYSize = dfReqYSize;
1446 342963 : if (psExtraArgIn)
1447 : {
1448 342963 : psExtraArg->pfnProgress = psExtraArgIn->pfnProgress;
1449 342963 : psExtraArg->pProgressData = psExtraArgIn->pProgressData;
1450 342963 : if (psExtraArgIn->nVersion >= 2)
1451 : {
1452 342963 : psExtraArg->bUseOnlyThisScale = psExtraArgIn->bUseOnlyThisScale;
1453 : }
1454 : }
1455 :
1456 342963 : GByte *pabyOut = static_cast<unsigned char *>(pData) +
1457 342963 : nOutXOff * nPixelSpace +
1458 342963 : static_cast<GPtrDiff_t>(nOutYOff) * nLineSpace;
1459 :
1460 342963 : auto l_band = GetRasterBand();
1461 342963 : if (!l_band)
1462 0 : return CE_Failure;
1463 :
1464 342963 : CPLErr eErr = CE_Failure;
1465 342963 : if (GDALDataTypeIsConversionLossy(l_band->GetRasterDataType(),
1466 342963 : 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 341733 : eErr = l_band->RasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize,
1494 : nReqYSize, pabyOut, nOutXSize, nOutYSize,
1495 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
1496 : }
1497 :
1498 342963 : 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 342963 : if (psExtraArg->pfnProgress)
1517 499 : psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
1518 :
1519 342963 : 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 5341 : 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 5341 : 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 5341 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1674 5341 : GDALRasterIOExtraArg *psExtraArg = &sExtraArg;
1675 :
1676 5341 : double dfXOff = nXOff;
1677 5341 : double dfYOff = nYOff;
1678 5341 : double dfXSize = nXSize;
1679 5341 : double dfYSize = nYSize;
1680 5341 : 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 5341 : double dfReqXOff = 0.0;
1690 5341 : double dfReqYOff = 0.0;
1691 5341 : double dfReqXSize = 0.0;
1692 5341 : double dfReqYSize = 0.0;
1693 5341 : int nReqXOff = 0;
1694 5341 : int nReqYOff = 0;
1695 5341 : int nReqXSize = 0;
1696 5341 : int nReqYSize = 0;
1697 :
1698 : // The window we will actual set _within_ the pData buffer.
1699 5341 : int nOutXOff = 0;
1700 5341 : int nOutYOff = 0;
1701 5341 : int nOutXSize = 0;
1702 5341 : int nOutYSize = 0;
1703 :
1704 5341 : if (!m_osResampling.empty())
1705 : {
1706 1556 : psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
1707 : }
1708 3785 : else if (psExtraArgIn != nullptr)
1709 : {
1710 3785 : psExtraArg->eResampleAlg = psExtraArgIn->eResampleAlg;
1711 : }
1712 :
1713 5341 : bool bError = false;
1714 5341 : 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 4502 : auto l_band = GetRasterBand();
1724 4502 : if (!l_band)
1725 0 : return CE_Failure;
1726 :
1727 4502 : GDALDataset *poDS = l_band->GetDataset();
1728 4502 : if (poDS == nullptr)
1729 0 : return CE_Failure;
1730 :
1731 4502 : psExtraArg->bFloatingPointWindowValidity = TRUE;
1732 4502 : psExtraArg->dfXOff = dfReqXOff;
1733 4502 : psExtraArg->dfYOff = dfReqYOff;
1734 4502 : psExtraArg->dfXSize = dfReqXSize;
1735 4502 : psExtraArg->dfYSize = dfReqYSize;
1736 4502 : if (psExtraArgIn)
1737 : {
1738 4502 : psExtraArg->pfnProgress = psExtraArgIn->pfnProgress;
1739 4502 : psExtraArg->pProgressData = psExtraArgIn->pProgressData;
1740 : }
1741 :
1742 4502 : GByte *pabyOut = static_cast<unsigned char *>(pData) +
1743 4502 : nOutXOff * nPixelSpace +
1744 4502 : static_cast<GPtrDiff_t>(nOutYOff) * nLineSpace;
1745 :
1746 4502 : CPLErr eErr = CE_Failure;
1747 :
1748 4502 : if (GDALDataTypeIsConversionLossy(l_band->GetRasterDataType(),
1749 4501 : 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 4499 : 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 4502 : 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 4502 : if (psExtraArg->pfnProgress)
1817 4024 : psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
1818 :
1819 4502 : return eErr;
1820 : }
1821 :
1822 : /************************************************************************/
1823 : /* SetResampling() */
1824 : /************************************************************************/
1825 :
1826 70241 : void VRTSimpleSource::SetResampling(const char *pszResampling)
1827 : {
1828 70241 : m_osResampling = (pszResampling) ? pszResampling : "";
1829 70241 : }
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 104137 : const char *VRTAveragedSource::GetTypeStatic()
1850 : {
1851 : static const char *TYPE = "AveragedSource";
1852 104137 : 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 13827 : const char *VRTComplexSource::GetTypeStatic()
2659 : {
2660 : static const char *TYPE = "ComplexSource";
2661 13827 : return TYPE;
2662 : }
2663 :
2664 : /************************************************************************/
2665 : /* GetType() */
2666 : /************************************************************************/
2667 :
2668 4021 : const char *VRTComplexSource::GetType() const
2669 : {
2670 4021 : 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 5057 : double VRTComplexSource::GetAdjustedNoDataValue() const
2696 : {
2697 5057 : 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 5034 : return m_dfNoDataValue;
2706 : }
2707 :
2708 : /************************************************************************/
2709 : /* SerializeToXML() */
2710 : /************************************************************************/
2711 :
2712 67 : CPLXMLNode *VRTComplexSource::SerializeToXML(const char *pszVRTPath)
2713 :
2714 : {
2715 67 : CPLXMLNode *psSrc = VRTSimpleSource::SerializeToXML(pszVRTPath);
2716 :
2717 67 : if (psSrc == nullptr)
2718 0 : return nullptr;
2719 :
2720 67 : CPLFree(psSrc->pszValue);
2721 67 : psSrc->pszValue = CPLStrdup(GetTypeStatic());
2722 :
2723 67 : if ((m_nProcessingFlags & PROCESSING_FLAG_USE_MASK_BAND) != 0)
2724 : {
2725 33 : CPLSetXMLValue(psSrc, "UseMaskBand", "true");
2726 : }
2727 :
2728 67 : 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 67 : 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 66 : 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 67 : if (!m_adfLUTInputs.empty())
2774 : {
2775 : // Make sure we print with sufficient precision to address really close
2776 : // entries (#6422).
2777 0 : CPLString osLUT;
2778 0 : if (m_adfLUTInputs.size() >= 2 &&
2779 0 : CPLString().Printf("%g", m_adfLUTInputs[0]) ==
2780 0 : 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 0 : osLUT = CPLString().Printf("%g:%g", m_adfLUTInputs[0],
2788 0 : m_adfLUTOutputs[0]);
2789 : }
2790 0 : for (size_t i = 1; i < m_adfLUTInputs.size(); i++)
2791 : {
2792 0 : if (CPLString().Printf("%g", m_adfLUTInputs[i]) ==
2793 0 : CPLString().Printf("%g", m_adfLUTInputs[i - 1]) ||
2794 0 : (i + 1 < m_adfLUTInputs.size() &&
2795 0 : CPLString().Printf("%g", m_adfLUTInputs[i]) ==
2796 0 : 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 0 : osLUT += CPLString().Printf(",%.17g:%g", m_adfLUTInputs[i],
2802 0 : m_adfLUTOutputs[i]);
2803 : }
2804 : else
2805 : {
2806 0 : osLUT += CPLString().Printf(",%g:%g", m_adfLUTInputs[i],
2807 0 : m_adfLUTOutputs[i]);
2808 : }
2809 : }
2810 0 : CPLSetXMLValue(psSrc, "LUT", osLUT);
2811 : }
2812 :
2813 67 : if (m_nColorTableComponent)
2814 : {
2815 7 : CPLSetXMLValue(psSrc, "ColorTableComponent",
2816 : CPLSPrintf("%d", m_nColorTableComponent));
2817 : }
2818 :
2819 67 : return psSrc;
2820 : }
2821 :
2822 : /************************************************************************/
2823 : /* XMLInit() */
2824 : /************************************************************************/
2825 :
2826 268 : 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 268 : VRTSimpleSource::XMLInit(psSrc, pszVRTPath, oMapSharedSources);
2837 268 : if (eErr != CE_None)
2838 0 : return eErr;
2839 : }
2840 :
2841 : /* -------------------------------------------------------------------- */
2842 : /* Complex parameters. */
2843 : /* -------------------------------------------------------------------- */
2844 268 : const char *pszScaleOffset = CPLGetXMLValue(psSrc, "ScaleOffset", nullptr);
2845 268 : const char *pszScaleRatio = CPLGetXMLValue(psSrc, "ScaleRatio", nullptr);
2846 268 : 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 240 : else if (CPLGetXMLValue(psSrc, "Exponent", nullptr) != nullptr &&
2855 241 : 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 268 : if (const char *pszNODATA = CPLGetXMLValue(psSrc, "NODATA", nullptr))
2876 : {
2877 66 : m_nProcessingFlags |= PROCESSING_FLAG_NODATA;
2878 66 : m_osNoDataValueOri = pszNODATA;
2879 66 : m_dfNoDataValue = CPLAtofM(m_osNoDataValueOri.c_str());
2880 : }
2881 :
2882 268 : const char *pszUseMaskBand = CPLGetXMLValue(psSrc, "UseMaskBand", nullptr);
2883 268 : if (pszUseMaskBand && CPLTestBool(pszUseMaskBand))
2884 : {
2885 41 : m_nProcessingFlags |= PROCESSING_FLAG_USE_MASK_BAND;
2886 : }
2887 :
2888 268 : const char *pszLUT = CPLGetXMLValue(psSrc, "LUT", nullptr);
2889 268 : if (pszLUT)
2890 : {
2891 : const CPLStringList aosValues(
2892 81 : CSLTokenizeString2(pszLUT, ",:", CSLT_ALLOWEMPTYTOKENS));
2893 :
2894 81 : const int nLUTItemCount = aosValues.size() / 2;
2895 : try
2896 : {
2897 81 : m_adfLUTInputs.resize(nLUTItemCount);
2898 81 : 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 780 : for (int nIndex = 0; nIndex < nLUTItemCount; nIndex++)
2909 : {
2910 699 : m_adfLUTInputs[nIndex] = CPLAtof(aosValues[nIndex * 2]);
2911 699 : m_adfLUTOutputs[nIndex] = CPLAtof(aosValues[nIndex * 2 + 1]);
2912 :
2913 : // Enforce the requirement that the LUT input array is
2914 : // monotonically non-decreasing.
2915 699 : 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 1317 : else if (nIndex > 0 &&
2925 618 : 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 81 : m_nProcessingFlags |= PROCESSING_FLAG_LUT;
2937 : }
2938 :
2939 : const char *pszColorTableComponent =
2940 268 : CPLGetXMLValue(psSrc, "ColorTableComponent", nullptr);
2941 268 : if (pszColorTableComponent)
2942 : {
2943 15 : m_nColorTableComponent = atoi(pszColorTableComponent);
2944 15 : m_nProcessingFlags |= PROCESSING_FLAG_COLOR_TABLE_EXPANSION;
2945 : }
2946 :
2947 268 : return CE_None;
2948 : }
2949 :
2950 : /************************************************************************/
2951 : /* LookupValue() */
2952 : /************************************************************************/
2953 :
2954 905438 : double VRTComplexSource::LookupValue(double dfInput)
2955 : {
2956 905438 : auto beginIter = m_adfLUTInputs.begin();
2957 905438 : auto endIter = m_adfLUTInputs.end();
2958 905438 : size_t offset = 0;
2959 905438 : if (std::isnan(m_adfLUTInputs[0]))
2960 : {
2961 6 : if (std::isnan(dfInput) || m_adfLUTInputs.size() == 1)
2962 1 : return m_adfLUTOutputs[0];
2963 5 : ++beginIter;
2964 5 : offset = 1;
2965 : }
2966 :
2967 : // Find the index of the first element in the LUT input array that
2968 : // is not smaller than the input value.
2969 : const size_t i =
2970 : offset +
2971 905437 : std::distance(beginIter, std::lower_bound(beginIter, endIter, dfInput));
2972 :
2973 905437 : if (i == offset)
2974 129039 : return m_adfLUTOutputs[offset];
2975 :
2976 : // If the index is beyond the end of the LUT input array, the input
2977 : // value is larger than all the values in the array.
2978 776398 : if (i == m_adfLUTInputs.size())
2979 8 : return m_adfLUTOutputs.back();
2980 :
2981 776390 : if (m_adfLUTInputs[i] == dfInput)
2982 268324 : return m_adfLUTOutputs[i];
2983 :
2984 : // Otherwise, interpolate.
2985 508066 : return m_adfLUTOutputs[i - 1] +
2986 508066 : (dfInput - m_adfLUTInputs[i - 1]) *
2987 508066 : ((m_adfLUTOutputs[i] - m_adfLUTOutputs[i - 1]) /
2988 508066 : (m_adfLUTInputs[i] - m_adfLUTInputs[i - 1]));
2989 : }
2990 :
2991 : /************************************************************************/
2992 : /* SetLinearScaling() */
2993 : /************************************************************************/
2994 :
2995 116 : void VRTComplexSource::SetLinearScaling(double dfOffset, double dfScale)
2996 : {
2997 116 : m_nProcessingFlags &= ~PROCESSING_FLAG_SCALING_EXPONENTIAL;
2998 116 : m_nProcessingFlags |= PROCESSING_FLAG_SCALING_LINEAR;
2999 116 : m_dfScaleOff = dfOffset;
3000 116 : m_dfScaleRatio = dfScale;
3001 116 : }
3002 :
3003 : /************************************************************************/
3004 : /* SetPowerScaling() */
3005 : /************************************************************************/
3006 :
3007 26 : void VRTComplexSource::SetPowerScaling(double dfExponentIn, double dfSrcMinIn,
3008 : double dfSrcMaxIn, double dfDstMinIn,
3009 : double dfDstMaxIn, bool bClip)
3010 : {
3011 26 : m_nProcessingFlags &= ~PROCESSING_FLAG_SCALING_LINEAR;
3012 26 : m_nProcessingFlags |= PROCESSING_FLAG_SCALING_EXPONENTIAL;
3013 26 : m_dfExponent = dfExponentIn;
3014 26 : m_dfSrcMin = dfSrcMinIn;
3015 26 : m_dfSrcMax = dfSrcMaxIn;
3016 26 : m_dfDstMin = dfDstMinIn;
3017 26 : m_dfDstMax = dfDstMaxIn;
3018 26 : m_bSrcMinMaxDefined = true;
3019 26 : m_bClip = bClip;
3020 26 : }
3021 :
3022 : /************************************************************************/
3023 : /* SetColorTableComponent() */
3024 : /************************************************************************/
3025 :
3026 255 : void VRTComplexSource::SetColorTableComponent(int nComponent)
3027 : {
3028 255 : m_nProcessingFlags |= PROCESSING_FLAG_COLOR_TABLE_EXPANSION;
3029 255 : m_nColorTableComponent = nComponent;
3030 255 : }
3031 :
3032 : /************************************************************************/
3033 : /* RasterIO() */
3034 : /************************************************************************/
3035 :
3036 16022 : CPLErr VRTComplexSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
3037 : int nYOff, int nXSize, int nYSize,
3038 : void *pData, int nBufXSize, int nBufYSize,
3039 : GDALDataType eBufType, GSpacing nPixelSpace,
3040 : GSpacing nLineSpace,
3041 : GDALRasterIOExtraArg *psExtraArgIn,
3042 : WorkingState &oWorkingState)
3043 :
3044 : {
3045 : GDALRasterIOExtraArg sExtraArg;
3046 16022 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
3047 16022 : GDALRasterIOExtraArg *psExtraArg = &sExtraArg;
3048 :
3049 16022 : double dfXOff = nXOff;
3050 16022 : double dfYOff = nYOff;
3051 16022 : double dfXSize = nXSize;
3052 16022 : double dfYSize = nYSize;
3053 16022 : if (psExtraArgIn != nullptr && psExtraArgIn->bFloatingPointWindowValidity)
3054 : {
3055 87 : dfXOff = psExtraArgIn->dfXOff;
3056 87 : dfYOff = psExtraArgIn->dfYOff;
3057 87 : dfXSize = psExtraArgIn->dfXSize;
3058 87 : dfYSize = psExtraArgIn->dfYSize;
3059 : }
3060 :
3061 : // The window we will actually request from the source raster band.
3062 16022 : double dfReqXOff = 0.0;
3063 16022 : double dfReqYOff = 0.0;
3064 16022 : double dfReqXSize = 0.0;
3065 16022 : double dfReqYSize = 0.0;
3066 16022 : int nReqXOff = 0;
3067 16022 : int nReqYOff = 0;
3068 16022 : int nReqXSize = 0;
3069 16022 : int nReqYSize = 0;
3070 :
3071 : // The window we will actual set _within_ the pData buffer.
3072 16022 : int nOutXOff = 0;
3073 16022 : int nOutYOff = 0;
3074 16022 : int nOutXSize = 0;
3075 16022 : int nOutYSize = 0;
3076 :
3077 16022 : if (!m_osResampling.empty())
3078 : {
3079 28 : psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
3080 : }
3081 15994 : else if (psExtraArgIn != nullptr)
3082 : {
3083 15994 : psExtraArg->eResampleAlg = psExtraArgIn->eResampleAlg;
3084 : }
3085 :
3086 16022 : bool bError = false;
3087 16022 : if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
3088 : psExtraArg->eResampleAlg, &dfReqXOff, &dfReqYOff,
3089 : &dfReqXSize, &dfReqYSize, &nReqXOff, &nReqYOff,
3090 : &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
3091 : &nOutXSize, &nOutYSize, bError))
3092 : {
3093 10949 : return bError ? CE_Failure : CE_None;
3094 : }
3095 : #if DEBUG_VERBOSE
3096 : CPLDebug("VRT",
3097 : "nXOff=%d, nYOff=%d, nXSize=%d, nYSize=%d, nBufXSize=%d, "
3098 : "nBufYSize=%d,\n"
3099 : "dfReqXOff=%g, dfReqYOff=%g, dfReqXSize=%g, dfReqYSize=%g,\n"
3100 : "nReqXOff=%d, nReqYOff=%d, nReqXSize=%d, nReqYSize=%d,\n"
3101 : "nOutXOff=%d, nOutYOff=%d, nOutXSize=%d, nOutYSize=%d",
3102 : nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize, dfReqXOff,
3103 : dfReqYOff, dfReqXSize, dfReqYSize, nReqXOff, nReqYOff, nReqXSize,
3104 : nReqYSize, nOutXOff, nOutYOff, nOutXSize, nOutYSize);
3105 : #endif
3106 :
3107 5073 : auto poSourceBand = GetRasterBand();
3108 5073 : if (!poSourceBand)
3109 0 : return CE_Failure;
3110 :
3111 5073 : psExtraArg->bFloatingPointWindowValidity = TRUE;
3112 5073 : psExtraArg->dfXOff = dfReqXOff;
3113 5073 : psExtraArg->dfYOff = dfReqYOff;
3114 5073 : psExtraArg->dfXSize = dfReqXSize;
3115 5073 : psExtraArg->dfYSize = dfReqYSize;
3116 5073 : if (psExtraArgIn)
3117 : {
3118 5073 : psExtraArg->pfnProgress = psExtraArgIn->pfnProgress;
3119 5073 : psExtraArg->pProgressData = psExtraArgIn->pProgressData;
3120 : }
3121 :
3122 5073 : GByte *const pabyOut = static_cast<GByte *>(pData) +
3123 5073 : nPixelSpace * nOutXOff +
3124 5073 : static_cast<GPtrDiff_t>(nLineSpace) * nOutYOff;
3125 5073 : const auto eSourceType = poSourceBand->GetRasterDataType();
3126 5073 : if (m_nProcessingFlags == PROCESSING_FLAG_NODATA)
3127 : {
3128 : // Optimization if doing only nodata processing
3129 87 : if (eSourceType == GDT_UInt8)
3130 : {
3131 45 : if (!GDALIsValueInRange<GByte>(m_dfNoDataValue))
3132 : {
3133 7 : return VRTSimpleSource::RasterIO(
3134 : eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData,
3135 : nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace,
3136 7 : psExtraArgIn, oWorkingState);
3137 : }
3138 :
3139 38 : return RasterIOProcessNoData<GByte, GDT_UInt8>(
3140 : poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3141 : nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3142 38 : nLineSpace, psExtraArg, oWorkingState);
3143 : }
3144 42 : else if (eSourceType == GDT_Int16)
3145 : {
3146 6 : if (!GDALIsValueInRange<int16_t>(m_dfNoDataValue))
3147 : {
3148 1 : return VRTSimpleSource::RasterIO(
3149 : eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData,
3150 : nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace,
3151 1 : psExtraArgIn, oWorkingState);
3152 : }
3153 :
3154 5 : return RasterIOProcessNoData<int16_t, GDT_Int16>(
3155 : poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3156 : nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3157 5 : nLineSpace, psExtraArg, oWorkingState);
3158 : }
3159 36 : else if (eSourceType == GDT_UInt16)
3160 : {
3161 11 : if (!GDALIsValueInRange<uint16_t>(m_dfNoDataValue))
3162 : {
3163 1 : return VRTSimpleSource::RasterIO(
3164 : eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData,
3165 : nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace,
3166 1 : psExtraArgIn, oWorkingState);
3167 : }
3168 :
3169 10 : return RasterIOProcessNoData<uint16_t, GDT_UInt16>(
3170 : poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3171 : nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3172 10 : nLineSpace, psExtraArg, oWorkingState);
3173 : }
3174 : }
3175 :
3176 5011 : const bool bIsComplex = CPL_TO_BOOL(GDALDataTypeIsComplex(eBufType));
3177 : CPLErr eErr;
3178 : // For Int32, float32 isn't sufficiently precise as working data type
3179 5011 : if (eVRTBandDataType == GDT_CInt32 || eVRTBandDataType == GDT_CFloat64 ||
3180 5007 : eVRTBandDataType == GDT_Int32 || eVRTBandDataType == GDT_UInt32 ||
3181 5003 : eVRTBandDataType == GDT_Int64 || eVRTBandDataType == GDT_UInt64 ||
3182 4978 : eVRTBandDataType == GDT_Float64 || eSourceType == GDT_Int32 ||
3183 4977 : eSourceType == GDT_UInt32 || eSourceType == GDT_Int64 ||
3184 : eSourceType == GDT_UInt64)
3185 : {
3186 35 : eErr = RasterIOInternal<double>(
3187 : poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3188 : nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3189 : nLineSpace, psExtraArg, bIsComplex ? GDT_CFloat64 : GDT_Float64,
3190 : oWorkingState);
3191 : }
3192 : else
3193 : {
3194 4976 : eErr = RasterIOInternal<float>(
3195 : poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
3196 : nReqYSize, pabyOut, nOutXSize, nOutYSize, eBufType, nPixelSpace,
3197 : nLineSpace, psExtraArg, bIsComplex ? GDT_CFloat32 : GDT_Float32,
3198 : oWorkingState);
3199 : }
3200 :
3201 5011 : if (psExtraArg->pfnProgress)
3202 56 : psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
3203 :
3204 5011 : return eErr;
3205 : }
3206 :
3207 : /************************************************************************/
3208 : /* hasZeroByte() */
3209 : /************************************************************************/
3210 :
3211 : CPL_NOSANITIZE_UNSIGNED_INT_OVERFLOW
3212 4392380 : static inline bool hasZeroByte(uint32_t v)
3213 : {
3214 : // Cf https://graphics.stanford.edu/~seander/bithacks.html#ZeroInWord
3215 4392380 : return (((v)-0x01010101U) & ~(v)&0x80808080U) != 0;
3216 : }
3217 :
3218 : /************************************************************************/
3219 : /* CopyWord() */
3220 : /************************************************************************/
3221 :
3222 : template <class SrcType>
3223 7644159 : static void CopyWord(const SrcType *pSrcVal, GDALDataType eSrcType,
3224 : void *pDstVal, GDALDataType eDstType)
3225 : {
3226 7644159 : switch (eDstType)
3227 : {
3228 2384950 : case GDT_UInt8:
3229 2384950 : GDALCopyWord(*pSrcVal, *static_cast<uint8_t *>(pDstVal));
3230 2384950 : break;
3231 82 : case GDT_Int8:
3232 82 : GDALCopyWord(*pSrcVal, *static_cast<int8_t *>(pDstVal));
3233 82 : break;
3234 2072520 : case GDT_UInt16:
3235 2072520 : GDALCopyWord(*pSrcVal, *static_cast<uint16_t *>(pDstVal));
3236 2072520 : break;
3237 2925180 : case GDT_Int16:
3238 2925180 : GDALCopyWord(*pSrcVal, *static_cast<int16_t *>(pDstVal));
3239 2925180 : break;
3240 22 : case GDT_UInt32:
3241 22 : GDALCopyWord(*pSrcVal, *static_cast<uint32_t *>(pDstVal));
3242 22 : break;
3243 641 : case GDT_Int32:
3244 641 : GDALCopyWord(*pSrcVal, *static_cast<int32_t *>(pDstVal));
3245 641 : break;
3246 22 : case GDT_UInt64:
3247 22 : GDALCopyWord(*pSrcVal, *static_cast<uint64_t *>(pDstVal));
3248 22 : break;
3249 42 : case GDT_Int64:
3250 42 : GDALCopyWord(*pSrcVal, *static_cast<int64_t *>(pDstVal));
3251 42 : break;
3252 252738 : case GDT_Float32:
3253 252738 : GDALCopyWord(*pSrcVal, *static_cast<float *>(pDstVal));
3254 252738 : break;
3255 6940 : case GDT_Float64:
3256 6940 : GDALCopyWord(*pSrcVal, *static_cast<double *>(pDstVal));
3257 6940 : break;
3258 1019 : default:
3259 1019 : GDALCopyWords(pSrcVal, eSrcType, 0, pDstVal, eDstType, 0, 1);
3260 1019 : break;
3261 : }
3262 7644159 : }
3263 :
3264 : /************************************************************************/
3265 : /* RasterIOProcessNoData() */
3266 : /************************************************************************/
3267 :
3268 : // This method is an optimization of the generic RasterIOInternal()
3269 : // that deals with a VRTComplexSource with only a NODATA value in it, and
3270 : // no other processing flags.
3271 :
3272 : // nReqXOff, nReqYOff, nReqXSize, nReqYSize are expressed in source band
3273 : // referential.
3274 : template <class SourceDT, GDALDataType eSourceType>
3275 53 : CPLErr VRTComplexSource::RasterIOProcessNoData(
3276 : GDALRasterBand *poSourceBand, GDALDataType eVRTBandDataType, int nReqXOff,
3277 : int nReqYOff, int nReqXSize, int nReqYSize, void *pData, int nOutXSize,
3278 : int nOutYSize, GDALDataType eBufType, GSpacing nPixelSpace,
3279 : GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg,
3280 : WorkingState &oWorkingState)
3281 : {
3282 53 : CPLAssert(m_nProcessingFlags == PROCESSING_FLAG_NODATA);
3283 53 : CPLAssert(GDALIsValueInRange<SourceDT>(m_dfNoDataValue));
3284 :
3285 : /* -------------------------------------------------------------------- */
3286 : /* Read into a temporary buffer. */
3287 : /* -------------------------------------------------------------------- */
3288 : try
3289 : {
3290 : // Cannot overflow since pData should at least have that number of
3291 : // elements
3292 53 : const size_t nPixelCount = static_cast<size_t>(nOutXSize) * nOutYSize;
3293 53 : if (nPixelCount >
3294 53 : static_cast<size_t>(std::numeric_limits<ptrdiff_t>::max()) /
3295 : sizeof(SourceDT))
3296 : {
3297 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
3298 : "Too large temporary buffer");
3299 0 : return CE_Failure;
3300 : }
3301 53 : oWorkingState.m_abyWrkBuffer.resize(sizeof(SourceDT) * nPixelCount);
3302 : }
3303 0 : catch (const std::bad_alloc &e)
3304 : {
3305 0 : CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
3306 0 : return CE_Failure;
3307 : }
3308 : const auto paSrcData =
3309 53 : reinterpret_cast<const SourceDT *>(oWorkingState.m_abyWrkBuffer.data());
3310 :
3311 53 : const GDALRIOResampleAlg eResampleAlgBack = psExtraArg->eResampleAlg;
3312 53 : if (!m_osResampling.empty())
3313 : {
3314 0 : psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(m_osResampling);
3315 : }
3316 :
3317 53 : const CPLErr eErr = poSourceBand->RasterIO(
3318 : GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
3319 53 : oWorkingState.m_abyWrkBuffer.data(), nOutXSize, nOutYSize, eSourceType,
3320 15 : sizeof(SourceDT), sizeof(SourceDT) * static_cast<GSpacing>(nOutXSize),
3321 : psExtraArg);
3322 53 : if (!m_osResampling.empty())
3323 0 : psExtraArg->eResampleAlg = eResampleAlgBack;
3324 :
3325 53 : if (eErr != CE_None)
3326 : {
3327 0 : return eErr;
3328 : }
3329 :
3330 53 : const auto nNoDataValue = static_cast<SourceDT>(m_dfNoDataValue);
3331 53 : size_t idxBuffer = 0;
3332 87 : if (eSourceType == eBufType &&
3333 34 : !GDALDataTypeIsConversionLossy(eSourceType, eVRTBandDataType))
3334 : {
3335 : // Most optimized case: the output type is the same as the source type,
3336 : // and conversion from the source type to the VRT band data type is
3337 : // not lossy
3338 346742 : for (int iY = 0; iY < nOutYSize; iY++)
3339 : {
3340 19624 : GByte *pDstLocation = static_cast<GByte *>(pData) +
3341 19624 : static_cast<GPtrDiff_t>(nLineSpace) * iY;
3342 :
3343 19624 : int iX = 0;
3344 19597 : if (sizeof(SourceDT) == 1 && nPixelSpace == 1)
3345 : {
3346 : // Optimization to detect more quickly if source pixels are
3347 : // at nodata.
3348 19600 : const GByte byNoDataValue = static_cast<GByte>(nNoDataValue);
3349 19600 : const uint32_t wordNoData =
3350 19600 : (static_cast<uint32_t>(byNoDataValue) << 24) |
3351 19600 : (byNoDataValue << 16) | (byNoDataValue << 8) |
3352 19600 : byNoDataValue;
3353 :
3354 : // Warning: hasZeroByte() assumes WORD_SIZE = 4
3355 19600 : constexpr int WORD_SIZE = 4;
3356 4420120 : for (; iX < nOutXSize - (WORD_SIZE - 1); iX += WORD_SIZE)
3357 : {
3358 : uint32_t v;
3359 : static_assert(sizeof(v) == WORD_SIZE,
3360 : "sizeof(v) == WORD_SIZE");
3361 4073440 : memcpy(&v, paSrcData + idxBuffer, sizeof(v));
3362 : // Cf https://graphics.stanford.edu/~seander/bithacks.html#ValueInWord
3363 4073440 : if (!hasZeroByte(v ^ wordNoData))
3364 : {
3365 : // No bytes are at nodata
3366 4291900 : memcpy(pDstLocation, &v, WORD_SIZE);
3367 4291900 : idxBuffer += WORD_SIZE;
3368 4291900 : pDstLocation += WORD_SIZE;
3369 : }
3370 108622 : else if (v == wordNoData)
3371 : {
3372 : // All bytes are at nodata
3373 131075 : idxBuffer += WORD_SIZE;
3374 131075 : pDstLocation += WORD_SIZE;
3375 : }
3376 : else
3377 : {
3378 : // There are both bytes at nodata and valid bytes
3379 0 : for (int k = 0; k < WORD_SIZE; ++k)
3380 : {
3381 48 : if (paSrcData[idxBuffer] != nNoDataValue)
3382 : {
3383 36 : memcpy(pDstLocation, &paSrcData[idxBuffer],
3384 : sizeof(SourceDT));
3385 : }
3386 48 : idxBuffer++;
3387 48 : pDstLocation += nPixelSpace;
3388 : }
3389 : }
3390 : }
3391 : }
3392 :
3393 350553 : for (; iX < nOutXSize;
3394 3844 : iX++, pDstLocation += nPixelSpace, idxBuffer++)
3395 : {
3396 3844 : if (paSrcData[idxBuffer] != nNoDataValue)
3397 : {
3398 3814 : memcpy(pDstLocation, &paSrcData[idxBuffer],
3399 : sizeof(SourceDT));
3400 : }
3401 : }
3402 : }
3403 : }
3404 20 : else if (!GDALDataTypeIsConversionLossy(eSourceType, eVRTBandDataType))
3405 : {
3406 : // Conversion from the source type to the VRT band data type is
3407 : // not lossy, so we can directly convert from the source type to
3408 : // the the output type
3409 124 : for (int iY = 0; iY < nOutYSize; iY++)
3410 : {
3411 105 : GByte *pDstLocation = static_cast<GByte *>(pData) +
3412 105 : static_cast<GPtrDiff_t>(nLineSpace) * iY;
3413 :
3414 1586 : for (int iX = 0; iX < nOutXSize;
3415 1481 : iX++, pDstLocation += nPixelSpace, idxBuffer++)
3416 : {
3417 1481 : if (paSrcData[idxBuffer] != nNoDataValue)
3418 : {
3419 672 : CopyWord(&paSrcData[idxBuffer], eSourceType, pDstLocation,
3420 : eBufType);
3421 : }
3422 : }
3423 : }
3424 : }
3425 : else
3426 : {
3427 : GByte abyTemp[2 * sizeof(double)];
3428 7 : for (int iY = 0; iY < nOutYSize; iY++)
3429 : {
3430 6 : GByte *pDstLocation = static_cast<GByte *>(pData) +
3431 6 : static_cast<GPtrDiff_t>(nLineSpace) * iY;
3432 :
3433 36 : for (int iX = 0; iX < nOutXSize;
3434 30 : iX++, pDstLocation += nPixelSpace, idxBuffer++)
3435 : {
3436 30 : if (paSrcData[idxBuffer] != nNoDataValue)
3437 : {
3438 : // Convert first to the VRTRasterBand data type
3439 : // to get its clamping, before outputting to buffer data type
3440 20 : CopyWord(&paSrcData[idxBuffer], eSourceType, abyTemp,
3441 : eVRTBandDataType);
3442 20 : GDALCopyWords(abyTemp, eVRTBandDataType, 0, pDstLocation,
3443 : eBufType, 0, 1);
3444 : }
3445 : }
3446 : }
3447 : }
3448 :
3449 327138 : return CE_None;
3450 : }
3451 :
3452 : /************************************************************************/
3453 : /* RasterIOInternal() */
3454 : /************************************************************************/
3455 :
3456 : // nReqXOff, nReqYOff, nReqXSize, nReqYSize are expressed in source band
3457 : // referential.
3458 : template <class WorkingDT>
3459 5056 : CPLErr VRTComplexSource::RasterIOInternal(
3460 : GDALRasterBand *poSourceBand, GDALDataType eVRTBandDataType, int nReqXOff,
3461 : int nReqYOff, int nReqXSize, int nReqYSize, void *pData, int nOutXSize,
3462 : int nOutYSize, GDALDataType eBufType, GSpacing nPixelSpace,
3463 : GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg,
3464 : GDALDataType eWrkDataType, WorkingState &oWorkingState)
3465 : {
3466 5056 : const GDALColorTable *poColorTable = nullptr;
3467 5056 : const bool bIsComplex = CPL_TO_BOOL(GDALDataTypeIsComplex(eBufType));
3468 5056 : const int nWordSize = GDALGetDataTypeSizeBytes(eWrkDataType);
3469 5056 : assert(nWordSize != 0);
3470 :
3471 : // If no explicit <NODATA> is set, but UseMaskBand is set, and the band
3472 : // has a nodata value, then use it as if it was set as <NODATA>
3473 5056 : int bNoDataSet = (m_nProcessingFlags & PROCESSING_FLAG_NODATA) != 0;
3474 5056 : double dfNoDataValue = GetAdjustedNoDataValue();
3475 :
3476 5163 : if ((m_nProcessingFlags & PROCESSING_FLAG_USE_MASK_BAND) != 0 &&
3477 107 : poSourceBand->GetMaskFlags() == GMF_NODATA)
3478 : {
3479 0 : dfNoDataValue = poSourceBand->GetNoDataValue(&bNoDataSet);
3480 : }
3481 :
3482 5056 : const bool bNoDataSetIsNan = bNoDataSet && std::isnan(dfNoDataValue);
3483 5056 : const bool bNoDataSetAndNotNan =
3484 5100 : bNoDataSet && !std::isnan(dfNoDataValue) &&
3485 44 : GDALIsValueInRange<WorkingDT>(dfNoDataValue);
3486 5056 : const auto fWorkingDataTypeNoData = static_cast<WorkingDT>(dfNoDataValue);
3487 :
3488 5056 : const GByte *pabyMask = nullptr;
3489 5056 : const WorkingDT *pafData = nullptr;
3490 5056 : if ((m_nProcessingFlags & PROCESSING_FLAG_SCALING_LINEAR) != 0 &&
3491 4462 : m_dfScaleRatio == 0 && bNoDataSet == FALSE &&
3492 4051 : (m_nProcessingFlags & PROCESSING_FLAG_USE_MASK_BAND) == 0)
3493 : {
3494 : /* ------------------------------------------------------------------ */
3495 : /* Optimization when writing a constant value */
3496 : /* (used by the -addalpha option of gdalbuildvrt) */
3497 : /* ------------------------------------------------------------------ */
3498 : // Already set to NULL when defined.
3499 : // pafData = NULL;
3500 : }
3501 : else
3502 : {
3503 : /* ---------------------------------------------------------------- */
3504 : /* Read into a temporary buffer. */
3505 : /* ---------------------------------------------------------------- */
3506 1005 : const size_t nPixelCount = static_cast<size_t>(nOutXSize) * nOutYSize;
3507 : try
3508 : {
3509 : // Cannot overflow since pData should at least have that number of
3510 : // elements
3511 1005 : if (nPixelCount >
3512 1005 : static_cast<size_t>(std::numeric_limits<ptrdiff_t>::max()) /
3513 1005 : static_cast<size_t>(nWordSize))
3514 : {
3515 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
3516 : "Too large temporary buffer");
3517 0 : return CE_Failure;
3518 : }
3519 1005 : oWorkingState.m_abyWrkBuffer.resize(nWordSize * nPixelCount);
3520 : }
3521 0 : catch (const std::bad_alloc &e)
3522 : {
3523 0 : CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
3524 0 : return CE_Failure;
3525 : }
3526 : pafData = reinterpret_cast<const WorkingDT *>(
3527 1005 : oWorkingState.m_abyWrkBuffer.data());
3528 :
3529 1005 : const GDALRIOResampleAlg eResampleAlgBack = psExtraArg->eResampleAlg;
3530 1005 : if (!m_osResampling.empty())
3531 : {
3532 28 : psExtraArg->eResampleAlg =
3533 28 : GDALRasterIOGetResampleAlg(m_osResampling);
3534 : }
3535 :
3536 1005 : const CPLErr eErr = poSourceBand->RasterIO(
3537 : GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
3538 1005 : oWorkingState.m_abyWrkBuffer.data(), nOutXSize, nOutYSize,
3539 : eWrkDataType, nWordSize,
3540 1005 : nWordSize * static_cast<GSpacing>(nOutXSize), psExtraArg);
3541 1005 : if (!m_osResampling.empty())
3542 28 : psExtraArg->eResampleAlg = eResampleAlgBack;
3543 :
3544 1005 : if (eErr != CE_None)
3545 : {
3546 0 : return eErr;
3547 : }
3548 :
3549 : // Allocate and read mask band if needed
3550 2962 : if (!bNoDataSet &&
3551 1112 : (m_nProcessingFlags & PROCESSING_FLAG_USE_MASK_BAND) != 0 &&
3552 107 : (poSourceBand->GetMaskFlags() != GMF_ALL_VALID ||
3553 52 : poSourceBand->GetColorInterpretation() == GCI_AlphaBand ||
3554 32 : GetMaskBandMainBand() != nullptr))
3555 : {
3556 : try
3557 : {
3558 107 : oWorkingState.m_abyWrkBufferMask.resize(nPixelCount);
3559 : }
3560 0 : catch (const std::exception &)
3561 : {
3562 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
3563 : "Out of memory when allocating mask buffer");
3564 0 : return CE_Failure;
3565 : }
3566 : pabyMask = reinterpret_cast<const GByte *>(
3567 107 : oWorkingState.m_abyWrkBufferMask.data());
3568 52 : auto poMaskBand =
3569 107 : (poSourceBand->GetColorInterpretation() == GCI_AlphaBand ||
3570 87 : GetMaskBandMainBand() != nullptr)
3571 : ? poSourceBand
3572 55 : : poSourceBand->GetMaskBand();
3573 107 : if (poMaskBand->RasterIO(
3574 : GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
3575 107 : oWorkingState.m_abyWrkBufferMask.data(), nOutXSize,
3576 : nOutYSize, GDT_UInt8, 1, static_cast<GSpacing>(nOutXSize),
3577 107 : psExtraArg) != CE_None)
3578 : {
3579 0 : return CE_Failure;
3580 : }
3581 : }
3582 :
3583 1005 : if (m_nColorTableComponent != 0)
3584 : {
3585 274 : poColorTable = poSourceBand->GetColorTable();
3586 274 : if (poColorTable == nullptr)
3587 : {
3588 0 : CPLError(CE_Failure, CPLE_AppDefined,
3589 : "Source band has no color table.");
3590 0 : return CE_Failure;
3591 : }
3592 : }
3593 : }
3594 :
3595 : /* -------------------------------------------------------------------- */
3596 : /* Selectively copy into output buffer with nodata masking, */
3597 : /* and/or scaling. */
3598 : /* -------------------------------------------------------------------- */
3599 :
3600 5056 : const bool bTwoStepDataTypeConversion =
3601 5056 : CPL_TO_BOOL(
3602 10008 : GDALDataTypeIsConversionLossy(eWrkDataType, eVRTBandDataType)) &&
3603 4952 : !CPL_TO_BOOL(GDALDataTypeIsConversionLossy(eVRTBandDataType, eBufType));
3604 :
3605 5056 : size_t idxBuffer = 0;
3606 81870 : for (int iY = 0; iY < nOutYSize; iY++)
3607 : {
3608 76814 : GByte *pDstLocation = static_cast<GByte *>(pData) +
3609 76814 : static_cast<GPtrDiff_t>(nLineSpace) * iY;
3610 :
3611 13305784 : for (int iX = 0; iX < nOutXSize;
3612 13228959 : iX++, pDstLocation += nPixelSpace, idxBuffer++)
3613 : {
3614 : WorkingDT afResult[2];
3615 13228959 : if (pafData && !bIsComplex)
3616 : {
3617 7594266 : WorkingDT fResult = pafData[idxBuffer];
3618 7594266 : if (bNoDataSetIsNan && std::isnan(fResult))
3619 1188642 : continue;
3620 7595996 : if (bNoDataSetAndNotNan &&
3621 1756 : ARE_REAL_EQUAL(fResult, fWorkingDataTypeNoData))
3622 273 : continue;
3623 7593964 : if (pabyMask && pabyMask[idxBuffer] == 0)
3624 1188350 : continue;
3625 :
3626 6405624 : if (poColorTable)
3627 : {
3628 : const GDALColorEntry *poEntry =
3629 1933900 : poColorTable->GetColorEntry(static_cast<int>(fResult));
3630 1933900 : if (poEntry)
3631 : {
3632 1933900 : if (m_nColorTableComponent == 1)
3633 685985 : fResult = poEntry->c1;
3634 1247920 : else if (m_nColorTableComponent == 2)
3635 531408 : fResult = poEntry->c2;
3636 716509 : else if (m_nColorTableComponent == 3)
3637 531408 : fResult = poEntry->c3;
3638 185101 : else if (m_nColorTableComponent == 4)
3639 185101 : fResult = poEntry->c4;
3640 : }
3641 : else
3642 : {
3643 0 : CPLErrorOnce(CE_Failure, CPLE_AppDefined,
3644 : "No entry %d.", static_cast<int>(fResult));
3645 0 : continue;
3646 : }
3647 : }
3648 :
3649 6405624 : if ((m_nProcessingFlags & PROCESSING_FLAG_SCALING_LINEAR) != 0)
3650 : {
3651 2368452 : fResult = static_cast<WorkingDT>(fResult * m_dfScaleRatio +
3652 2368452 : m_dfScaleOff);
3653 : }
3654 4037162 : else if ((m_nProcessingFlags &
3655 : PROCESSING_FLAG_SCALING_EXPONENTIAL) != 0)
3656 : {
3657 3579 : if (!m_bSrcMinMaxDefined)
3658 : {
3659 1 : int bSuccessMin = FALSE;
3660 1 : int bSuccessMax = FALSE;
3661 2 : double adfMinMax[2] = {
3662 1 : poSourceBand->GetMinimum(&bSuccessMin),
3663 1 : poSourceBand->GetMaximum(&bSuccessMax)};
3664 2 : if ((bSuccessMin && bSuccessMax) ||
3665 1 : poSourceBand->ComputeRasterMinMax(
3666 : TRUE, adfMinMax) == CE_None)
3667 : {
3668 1 : m_dfSrcMin = adfMinMax[0];
3669 1 : m_dfSrcMax = adfMinMax[1];
3670 1 : m_bSrcMinMaxDefined = true;
3671 : }
3672 : else
3673 : {
3674 0 : CPLError(CE_Failure, CPLE_AppDefined,
3675 : "Cannot determine source min/max value");
3676 0 : return CE_Failure;
3677 : }
3678 : }
3679 :
3680 7158 : double dfPowVal = (m_dfSrcMin == m_dfSrcMax)
3681 3579 : ? 0
3682 3579 : : (fResult - m_dfSrcMin) /
3683 3579 : (m_dfSrcMax - m_dfSrcMin);
3684 3579 : if (m_bClip)
3685 : {
3686 3574 : if (dfPowVal < 0.0)
3687 1 : dfPowVal = 0.0;
3688 3573 : else if (dfPowVal > 1.0)
3689 700 : dfPowVal = 1.0;
3690 : }
3691 3579 : fResult =
3692 3579 : static_cast<WorkingDT>((m_dfDstMax - m_dfDstMin) *
3693 3579 : pow(dfPowVal, m_dfExponent) +
3694 3579 : m_dfDstMin);
3695 : }
3696 :
3697 6405624 : if (!m_adfLUTInputs.empty())
3698 905438 : fResult = static_cast<WorkingDT>(LookupValue(fResult));
3699 :
3700 6405624 : if (m_nMaxValue != 0 && fResult > m_nMaxValue)
3701 800 : fResult = static_cast<WorkingDT>(m_nMaxValue);
3702 :
3703 6405624 : afResult[0] = fResult;
3704 6405624 : afResult[1] = 0;
3705 : }
3706 5634663 : else if (pafData && bIsComplex)
3707 : {
3708 1018 : afResult[0] = pafData[2 * idxBuffer];
3709 1018 : afResult[1] = pafData[2 * idxBuffer + 1];
3710 :
3711 : // Do not use color table.
3712 1018 : if ((m_nProcessingFlags & PROCESSING_FLAG_SCALING_LINEAR) != 0)
3713 : {
3714 4 : afResult[0] = static_cast<WorkingDT>(
3715 4 : afResult[0] * m_dfScaleRatio + m_dfScaleOff);
3716 4 : afResult[1] = static_cast<WorkingDT>(
3717 4 : afResult[1] * m_dfScaleRatio + m_dfScaleOff);
3718 : }
3719 :
3720 : /* Do not use LUT */
3721 : }
3722 : else
3723 : {
3724 5633642 : afResult[0] = static_cast<WorkingDT>(m_dfScaleOff);
3725 5633642 : afResult[1] = 0;
3726 :
3727 5633642 : if (!m_adfLUTInputs.empty())
3728 0 : afResult[0] =
3729 0 : static_cast<WorkingDT>(LookupValue(afResult[0]));
3730 :
3731 5633642 : if (m_nMaxValue != 0 && afResult[0] > m_nMaxValue)
3732 0 : afResult[0] = static_cast<WorkingDT>(m_nMaxValue);
3733 : }
3734 :
3735 12040257 : if (eBufType == GDT_UInt8 && eVRTBandDataType == GDT_UInt8)
3736 : {
3737 4396790 : *pDstLocation = static_cast<GByte>(std::min(
3738 8793580 : 255.0f,
3739 4396790 : std::max(0.0f, static_cast<float>(afResult[0]) + 0.5f)));
3740 : }
3741 7643487 : else if (!bTwoStepDataTypeConversion)
3742 : {
3743 260514 : CopyWord<WorkingDT>(afResult, eWrkDataType, pDstLocation,
3744 : eBufType);
3745 : }
3746 7382978 : else if (eVRTBandDataType == GDT_Float32 && eBufType == GDT_Float64)
3747 : {
3748 : // Particular case of the below 2-step conversion.
3749 : // Helps a bit for some geolocation based warping with Sentinel3
3750 : // data where the longitude/latitude arrays are Int32 bands,
3751 : // rescaled in VRT as Float32 and requested as Float64
3752 : float fVal;
3753 20 : GDALCopyWord(afResult[0], fVal);
3754 20 : *reinterpret_cast<double *>(pDstLocation) = fVal;
3755 : }
3756 : else
3757 : {
3758 : GByte abyTemp[2 * sizeof(double)];
3759 : // Convert first to the VRTRasterBand data type
3760 : // to get its clamping, before outputting to buffer data type
3761 7382958 : CopyWord<WorkingDT>(afResult, eWrkDataType, abyTemp,
3762 : eVRTBandDataType);
3763 7382958 : GDALCopyWords(abyTemp, eVRTBandDataType, 0, pDstLocation,
3764 : eBufType, 0, 1);
3765 : }
3766 : }
3767 : }
3768 :
3769 5056 : return CE_None;
3770 : }
3771 :
3772 : // Explicitly instantiate template method, as it is used in another file.
3773 : template CPLErr VRTComplexSource::RasterIOInternal<float>(
3774 : GDALRasterBand *poSourceBand, GDALDataType eVRTBandDataType, int nReqXOff,
3775 : int nReqYOff, int nReqXSize, int nReqYSize, void *pData, int nOutXSize,
3776 : int nOutYSize, GDALDataType eBufType, GSpacing nPixelSpace,
3777 : GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg,
3778 : GDALDataType eWrkDataType, WorkingState &oWorkingState);
3779 :
3780 : /************************************************************************/
3781 : /* AreValuesUnchanged() */
3782 : /************************************************************************/
3783 :
3784 9 : bool VRTComplexSource::AreValuesUnchanged() const
3785 : {
3786 10 : return m_dfScaleOff == 0.0 && m_dfScaleRatio == 1.0 &&
3787 24 : m_adfLUTInputs.empty() && m_nColorTableComponent == 0 &&
3788 14 : (m_nProcessingFlags & PROCESSING_FLAG_SCALING_EXPONENTIAL) == 0;
3789 : }
3790 :
3791 : /************************************************************************/
3792 : /* GetMinimum() */
3793 : /************************************************************************/
3794 :
3795 2 : double VRTComplexSource::GetMinimum(int nXSize, int nYSize, int *pbSuccess)
3796 : {
3797 2 : if (AreValuesUnchanged())
3798 : {
3799 1 : return VRTSimpleSource::GetMinimum(nXSize, nYSize, pbSuccess);
3800 : }
3801 :
3802 1 : *pbSuccess = FALSE;
3803 1 : return 0;
3804 : }
3805 :
3806 : /************************************************************************/
3807 : /* GetMaximum() */
3808 : /************************************************************************/
3809 :
3810 2 : double VRTComplexSource::GetMaximum(int nXSize, int nYSize, int *pbSuccess)
3811 : {
3812 2 : if (AreValuesUnchanged())
3813 : {
3814 1 : return VRTSimpleSource::GetMaximum(nXSize, nYSize, pbSuccess);
3815 : }
3816 :
3817 1 : *pbSuccess = FALSE;
3818 1 : return 0;
3819 : }
3820 :
3821 : /************************************************************************/
3822 : /* GetHistogram() */
3823 : /************************************************************************/
3824 :
3825 0 : CPLErr VRTComplexSource::GetHistogram(int nXSize, int nYSize, double dfMin,
3826 : double dfMax, int nBuckets,
3827 : GUIntBig *panHistogram,
3828 : int bIncludeOutOfRange, int bApproxOK,
3829 : GDALProgressFunc pfnProgress,
3830 : void *pProgressData)
3831 : {
3832 0 : if (AreValuesUnchanged())
3833 : {
3834 0 : return VRTSimpleSource::GetHistogram(
3835 : nXSize, nYSize, dfMin, dfMax, nBuckets, panHistogram,
3836 0 : bIncludeOutOfRange, bApproxOK, pfnProgress, pProgressData);
3837 : }
3838 :
3839 0 : return CE_Failure;
3840 : }
3841 :
3842 : /************************************************************************/
3843 : /* ==================================================================== */
3844 : /* VRTFuncSource */
3845 : /* ==================================================================== */
3846 : /************************************************************************/
3847 :
3848 : /************************************************************************/
3849 : /* VRTFuncSource() */
3850 : /************************************************************************/
3851 :
3852 12 : VRTFuncSource::VRTFuncSource()
3853 : : pfnReadFunc(nullptr), pCBData(nullptr), eType(GDT_UInt8),
3854 12 : fNoDataValue(static_cast<float>(VRT_NODATA_UNSET))
3855 : {
3856 12 : }
3857 :
3858 : /************************************************************************/
3859 : /* ~VRTFuncSource() */
3860 : /************************************************************************/
3861 :
3862 24 : VRTFuncSource::~VRTFuncSource()
3863 : {
3864 24 : }
3865 :
3866 : /************************************************************************/
3867 : /* GetType() */
3868 : /************************************************************************/
3869 :
3870 0 : const char *VRTFuncSource::GetType() const
3871 : {
3872 : static const char *TYPE = "FuncSource";
3873 0 : return TYPE;
3874 : }
3875 :
3876 : /************************************************************************/
3877 : /* SerializeToXML() */
3878 : /************************************************************************/
3879 :
3880 0 : CPLXMLNode *VRTFuncSource::SerializeToXML(CPL_UNUSED const char *pszVRTPath)
3881 : {
3882 0 : return nullptr;
3883 : }
3884 :
3885 : /************************************************************************/
3886 : /* RasterIO() */
3887 : /************************************************************************/
3888 :
3889 10 : CPLErr VRTFuncSource::RasterIO(GDALDataType /*eVRTBandDataType*/, int nXOff,
3890 : int nYOff, int nXSize, int nYSize, void *pData,
3891 : int nBufXSize, int nBufYSize,
3892 : GDALDataType eBufType, GSpacing nPixelSpace,
3893 : GSpacing nLineSpace,
3894 : GDALRasterIOExtraArg * /* psExtraArg */,
3895 : WorkingState & /* oWorkingState */)
3896 : {
3897 10 : if (nPixelSpace == GDALGetDataTypeSizeBytes(eBufType) &&
3898 10 : nLineSpace == nPixelSpace * nXSize && nBufXSize == nXSize &&
3899 20 : nBufYSize == nYSize && eBufType == eType)
3900 : {
3901 10 : return pfnReadFunc(pCBData, nXOff, nYOff, nXSize, nYSize, pData);
3902 : }
3903 : else
3904 : {
3905 0 : CPLError(CE_Failure, CPLE_AppDefined,
3906 : "VRTFuncSource::RasterIO() - Irregular request.");
3907 0 : CPLDebug("VRT", "Irregular request: %d,%d %d,%d, %d,%d %d,%d %d,%d",
3908 : static_cast<int>(nPixelSpace),
3909 : GDALGetDataTypeSizeBytes(eBufType),
3910 : static_cast<int>(nLineSpace),
3911 : static_cast<int>(nPixelSpace) * nXSize, nBufXSize, nXSize,
3912 : nBufYSize, nYSize, static_cast<int>(eBufType),
3913 0 : static_cast<int>(eType));
3914 :
3915 0 : return CE_Failure;
3916 : }
3917 : }
3918 :
3919 : /************************************************************************/
3920 : /* GetMinimum() */
3921 : /************************************************************************/
3922 :
3923 0 : double VRTFuncSource::GetMinimum(int /* nXSize */, int /* nYSize */,
3924 : int *pbSuccess)
3925 : {
3926 0 : *pbSuccess = FALSE;
3927 0 : return 0;
3928 : }
3929 :
3930 : /************************************************************************/
3931 : /* GetMaximum() */
3932 : /************************************************************************/
3933 :
3934 0 : double VRTFuncSource::GetMaximum(int /* nXSize */, int /* nYSize */,
3935 : int *pbSuccess)
3936 : {
3937 0 : *pbSuccess = FALSE;
3938 0 : return 0;
3939 : }
3940 :
3941 : /************************************************************************/
3942 : /* GetHistogram() */
3943 : /************************************************************************/
3944 :
3945 0 : CPLErr VRTFuncSource::GetHistogram(
3946 : int /* nXSize */, int /* nYSize */, double /* dfMin */, double /* dfMax */,
3947 : int /* nBuckets */, GUIntBig * /* panHistogram */,
3948 : int /* bIncludeOutOfRange */, int /* bApproxOK */,
3949 : GDALProgressFunc /* pfnProgress */, void * /* pProgressData */)
3950 : {
3951 0 : return CE_Failure;
3952 : }
3953 :
3954 : /************************************************************************/
3955 : /* VRTParseCoreSources() */
3956 : /************************************************************************/
3957 :
3958 104126 : VRTSource *VRTParseCoreSources(const CPLXMLNode *psChild,
3959 : const char *pszVRTPath,
3960 : VRTMapSharedResources &oMapSharedSources)
3961 :
3962 : {
3963 104126 : VRTSource *poSource = nullptr;
3964 :
3965 208245 : if (EQUAL(psChild->pszValue, VRTAveragedSource::GetTypeStatic()) ||
3966 104119 : (EQUAL(psChild->pszValue, VRTSimpleSource::GetTypeStatic()) &&
3967 103857 : STARTS_WITH_CI(CPLGetXMLValue(psChild, "Resampling", "Nearest"),
3968 : "Aver")))
3969 : {
3970 16 : poSource = new VRTAveragedSource();
3971 : }
3972 104110 : else if (EQUAL(psChild->pszValue, VRTSimpleSource::GetTypeStatic()))
3973 : {
3974 103848 : poSource = new VRTSimpleSource();
3975 : }
3976 262 : else if (EQUAL(psChild->pszValue, VRTComplexSource::GetTypeStatic()))
3977 : {
3978 254 : poSource = new VRTComplexSource();
3979 : }
3980 8 : else if (EQUAL(psChild->pszValue, VRTNoDataFromMaskSource::GetTypeStatic()))
3981 : {
3982 8 : poSource = new VRTNoDataFromMaskSource();
3983 : }
3984 : else
3985 : {
3986 0 : CPLError(CE_Failure, CPLE_AppDefined,
3987 : "VRTParseCoreSources() - Unknown source : %s",
3988 0 : psChild->pszValue);
3989 0 : return nullptr;
3990 : }
3991 :
3992 104126 : if (poSource->XMLInit(psChild, pszVRTPath, oMapSharedSources) == CE_None)
3993 104123 : return poSource;
3994 :
3995 3 : delete poSource;
3996 3 : return nullptr;
3997 : }
3998 :
3999 : /*! @endcond */
|
