Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Virtual GDAL Datasets
4 : * Purpose: Implementation of VRTSourcedRasterBand
5 : * Author: Frank Warmerdam <warmerdam@pobox.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdal_vrt.h"
16 : #include "vrtdataset.h"
17 :
18 : #include <algorithm>
19 : #include <cassert>
20 : #include <cmath>
21 : #include <cstddef>
22 : #include <cstdio>
23 : #include <cstdlib>
24 : #include <cstring>
25 : #include <limits>
26 : #include <mutex>
27 : #include <set>
28 : #include <string>
29 :
30 : #include "cpl_conv.h"
31 : #include "cpl_error.h"
32 : #include "cpl_error_internal.h"
33 : #include "cpl_hash_set.h"
34 : #include "cpl_minixml.h"
35 : #include "cpl_progress.h"
36 : #include "cpl_quad_tree.h"
37 : #include "cpl_string.h"
38 : #include "cpl_vsi.h"
39 : #include "cpl_worker_thread_pool.h"
40 : #include "gdal.h"
41 : #include "gdalantirecursion.h"
42 : #include "gdal_priv.h"
43 : #include "gdal_thread_pool.h"
44 : #include "ogr_geometry.h"
45 :
46 : /*! @cond Doxygen_Suppress */
47 :
48 : /************************************************************************/
49 : /* ==================================================================== */
50 : /* VRTSourcedRasterBand */
51 : /* ==================================================================== */
52 : /************************************************************************/
53 :
54 : /************************************************************************/
55 : /* VRTSourcedRasterBand() */
56 : /************************************************************************/
57 :
58 2599 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn)
59 : {
60 2599 : VRTRasterBand::Initialize(poDSIn->GetRasterXSize(),
61 : poDSIn->GetRasterYSize());
62 :
63 2599 : poDS = poDSIn;
64 2599 : nBand = nBandIn;
65 2599 : }
66 :
67 : /************************************************************************/
68 : /* VRTSourcedRasterBand() */
69 : /************************************************************************/
70 :
71 0 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataType eType, int nXSize,
72 0 : int nYSize)
73 : {
74 0 : VRTRasterBand::Initialize(nXSize, nYSize);
75 :
76 0 : eDataType = eType;
77 0 : }
78 :
79 : /************************************************************************/
80 : /* VRTSourcedRasterBand() */
81 : /************************************************************************/
82 :
83 325 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn,
84 : GDALDataType eType, int nXSize,
85 325 : int nYSize)
86 325 : : VRTSourcedRasterBand(poDSIn, nBandIn, eType, nXSize, nYSize, 0, 0)
87 : {
88 325 : }
89 :
90 : /************************************************************************/
91 : /* VRTSourcedRasterBand() */
92 : /************************************************************************/
93 :
94 139157 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn,
95 : GDALDataType eType, int nXSize,
96 : int nYSize, int nBlockXSizeIn,
97 139157 : int nBlockYSizeIn)
98 : {
99 139157 : VRTRasterBand::Initialize(nXSize, nYSize);
100 :
101 139157 : poDS = poDSIn;
102 139157 : nBand = nBandIn;
103 :
104 139157 : eDataType = eType;
105 139157 : if (nBlockXSizeIn > 0)
106 138779 : nBlockXSize = nBlockXSizeIn;
107 139157 : if (nBlockYSizeIn > 0)
108 138797 : nBlockYSize = nBlockYSizeIn;
109 139157 : }
110 :
111 : /************************************************************************/
112 : /* ~VRTSourcedRasterBand() */
113 : /************************************************************************/
114 :
115 281667 : VRTSourcedRasterBand::~VRTSourcedRasterBand()
116 :
117 : {
118 141756 : VRTSourcedRasterBand::CloseDependentDatasets();
119 281667 : }
120 :
121 : /************************************************************************/
122 : /* CopyForCloneWithoutSources() */
123 : /************************************************************************/
124 :
125 102 : void VRTSourcedRasterBand::CopyForCloneWithoutSources(
126 : const VRTSourcedRasterBand *poSrcBand)
127 : {
128 102 : CopyCommonInfoFrom(poSrcBand);
129 102 : m_bNoDataValueSet = poSrcBand->m_bNoDataValueSet;
130 102 : m_dfNoDataValue = poSrcBand->m_dfNoDataValue;
131 102 : m_bHideNoDataValue = poSrcBand->m_bHideNoDataValue;
132 102 : }
133 :
134 : /************************************************************************/
135 : /* CloneWithoutSources() */
136 : /************************************************************************/
137 :
138 : std::unique_ptr<VRTSourcedRasterBand>
139 101 : VRTSourcedRasterBand::CloneWithoutSources(GDALDataset *poNewDS, int nNewXSize,
140 : int nNewYSize) const
141 : {
142 : auto poClone = std::make_unique<VRTSourcedRasterBand>(
143 101 : poNewDS, GetBand(), GetRasterDataType(), nNewXSize, nNewYSize);
144 101 : poClone->CopyForCloneWithoutSources(this);
145 101 : return poClone;
146 : }
147 :
148 : /************************************************************************/
149 : /* CanIRasterIOBeForwardedToEachSource() */
150 : /************************************************************************/
151 :
152 289889 : bool VRTSourcedRasterBand::CanIRasterIOBeForwardedToEachSource(
153 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
154 : int nBufXSize, int nBufYSize, const GDALRasterIOExtraArg *psExtraArg) const
155 : {
156 10335 : const auto IsNonNearestInvolved = [this, psExtraArg]
157 : {
158 5175 : if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
159 : {
160 15 : return true;
161 : }
162 14478 : for (const auto &poSource : m_papoSources)
163 : {
164 9326 : if (poSource->GetType() == VRTComplexSource::GetTypeStatic())
165 : {
166 : const auto *const poComplexSource =
167 3425 : static_cast<const VRTComplexSource *>(poSource.get());
168 : const auto &osSourceResampling =
169 3425 : poComplexSource->GetResampling();
170 3433 : if (!osSourceResampling.empty() &&
171 8 : osSourceResampling != "nearest")
172 8 : return true;
173 : }
174 : }
175 5152 : return false;
176 289889 : };
177 :
178 : // If resampling with non-nearest neighbour, we need to be careful
179 : // if the VRT band exposes a nodata value, but the sources do not have it.
180 : // To also avoid edge effects on sources when downsampling, use the
181 : // base implementation of IRasterIO() (that is acquiring sources at their
182 : // nominal resolution, and then downsampling), but only if none of the
183 : // contributing sources have overviews.
184 289889 : if (eRWFlag == GF_Read && (nXSize != nBufXSize || nYSize != nBufYSize) &&
185 579778 : !m_papoSources.empty() && IsNonNearestInvolved())
186 : {
187 23 : bool bSourceHasOverviews = false;
188 23 : const bool bIsDownsampling = (nBufXSize < nXSize && nBufYSize < nYSize);
189 23 : int nContributingSources = 0;
190 23 : bool bSourceFullySatisfiesRequest = true;
191 37 : for (auto &poSource : m_papoSources)
192 : {
193 26 : if (!poSource->IsSimpleSource())
194 : {
195 12 : return false;
196 : }
197 : else
198 : {
199 : VRTSimpleSource *const poSimpleSource =
200 26 : static_cast<VRTSimpleSource *>(poSource.get());
201 :
202 26 : if (poSimpleSource->GetType() ==
203 26 : VRTComplexSource::GetTypeStatic())
204 : {
205 16 : auto *const poComplexSource =
206 : static_cast<VRTComplexSource *>(poSimpleSource);
207 : const auto &osSourceResampling =
208 16 : poComplexSource->GetResampling();
209 24 : if (!osSourceResampling.empty() &&
210 8 : osSourceResampling != "nearest")
211 : {
212 : const int lMaskFlags =
213 : const_cast<VRTSourcedRasterBand *>(this)
214 8 : ->GetMaskFlags();
215 6 : if ((lMaskFlags != GMF_ALL_VALID &&
216 16 : lMaskFlags != GMF_NODATA) ||
217 2 : IsMaskBand())
218 : {
219 : // Unfortunately this will prevent using overviews
220 : // of the sources, but it is unpractical to use
221 : // them without serious implementation complications
222 12 : return false;
223 : }
224 : }
225 : }
226 :
227 18 : double dfXOff = nXOff;
228 18 : double dfYOff = nYOff;
229 18 : double dfXSize = nXSize;
230 18 : double dfYSize = nYSize;
231 18 : if (psExtraArg->bFloatingPointWindowValidity)
232 : {
233 1 : dfXOff = psExtraArg->dfXOff;
234 1 : dfYOff = psExtraArg->dfYOff;
235 1 : dfXSize = psExtraArg->dfXSize;
236 1 : dfYSize = psExtraArg->dfYSize;
237 : }
238 :
239 : // The window we will actually request from the source raster
240 : // band.
241 18 : double dfReqXOff = 0.0;
242 18 : double dfReqYOff = 0.0;
243 18 : double dfReqXSize = 0.0;
244 18 : double dfReqYSize = 0.0;
245 18 : int nReqXOff = 0;
246 18 : int nReqYOff = 0;
247 18 : int nReqXSize = 0;
248 18 : int nReqYSize = 0;
249 :
250 : // The window we will actual set _within_ the pData buffer.
251 18 : int nOutXOff = 0;
252 18 : int nOutYOff = 0;
253 18 : int nOutXSize = 0;
254 18 : int nOutYSize = 0;
255 :
256 18 : bool bError = false;
257 18 : if (!poSimpleSource->GetSrcDstWindow(
258 : dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
259 18 : psExtraArg->eResampleAlg, &dfReqXOff, &dfReqYOff,
260 : &dfReqXSize, &dfReqYSize, &nReqXOff, &nReqYOff,
261 : &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
262 : &nOutXSize, &nOutYSize, bError))
263 : {
264 0 : continue;
265 : }
266 18 : auto poBand = poSimpleSource->GetRasterBand();
267 18 : if (poBand == nullptr)
268 : {
269 0 : return false;
270 : }
271 18 : ++nContributingSources;
272 18 : if (!(nOutXOff == 0 && nOutYOff == 0 &&
273 17 : nOutXSize == nBufXSize && nOutYSize == nBufYSize))
274 4 : bSourceFullySatisfiesRequest = false;
275 18 : if (m_bNoDataValueSet)
276 : {
277 7 : int bSrcHasNoData = FALSE;
278 : const double dfSrcNoData =
279 7 : poBand->GetNoDataValue(&bSrcHasNoData);
280 7 : if (!bSrcHasNoData || dfSrcNoData != m_dfNoDataValue)
281 : {
282 4 : return false;
283 : }
284 : }
285 14 : if (bIsDownsampling)
286 : {
287 11 : if (poBand->GetOverviewCount() != 0)
288 : {
289 4 : bSourceHasOverviews = true;
290 : }
291 : }
292 : }
293 : }
294 11 : if (bIsDownsampling && !bSourceHasOverviews &&
295 3 : (nContributingSources > 1 || !bSourceFullySatisfiesRequest))
296 : {
297 2 : return false;
298 : }
299 : }
300 289875 : return true;
301 : }
302 :
303 : /************************************************************************/
304 : /* CanMultiThreadRasterIO() */
305 : /************************************************************************/
306 :
307 513 : bool VRTSourcedRasterBand::CanMultiThreadRasterIO(
308 : double dfXOff, double dfYOff, double dfXSize, double dfYSize,
309 : int &nContributingSources) const
310 : {
311 513 : int iLastSource = 0;
312 : CPLRectObj sSourceBounds;
313 513 : CPLQuadTree *hQuadTree = nullptr;
314 513 : bool bRet = true;
315 513 : std::set<std::string> oSetDSName;
316 :
317 513 : nContributingSources = 0;
318 1321 : for (int iSource = 0; iSource < static_cast<int>(m_papoSources.size());
319 : iSource++)
320 : {
321 812 : const auto &poSource = m_papoSources[iSource];
322 812 : if (!poSource->IsSimpleSource())
323 : {
324 0 : bRet = false;
325 0 : break;
326 : }
327 : const auto poSimpleSource =
328 812 : cpl::down_cast<VRTSimpleSource *>(poSource.get());
329 812 : if (poSimpleSource->DstWindowIntersects(dfXOff, dfYOff, dfXSize,
330 : dfYSize))
331 : {
332 : // Only build hQuadTree if there are 2 or more sources
333 547 : if (nContributingSources == 1)
334 : {
335 : std::string &oFirstSrcDSName =
336 28 : cpl::down_cast<VRTSimpleSource *>(
337 28 : m_papoSources[iLastSource].get())
338 28 : ->m_osSrcDSName;
339 28 : oSetDSName.insert(oFirstSrcDSName);
340 :
341 : CPLRectObj sGlobalBounds;
342 28 : sGlobalBounds.minx = dfXOff;
343 28 : sGlobalBounds.miny = dfYOff;
344 28 : sGlobalBounds.maxx = dfXOff + dfXSize;
345 28 : sGlobalBounds.maxy = dfYOff + dfYSize;
346 28 : hQuadTree = CPLQuadTreeCreate(&sGlobalBounds, nullptr);
347 :
348 28 : CPLQuadTreeInsertWithBounds(
349 : hQuadTree,
350 : reinterpret_cast<void *>(
351 28 : static_cast<uintptr_t>(iLastSource)),
352 : &sSourceBounds);
353 : }
354 :
355 : // Check there are not several sources with the same name, to avoid
356 : // the same GDALDataset* to be used from multiple threads. We may
357 : // be a bit too pessimistic, for example if working with unnamed
358 : // Memory datasets, but that would involve comparing
359 : // poSource->GetRasterBandNoOpen()->GetDataset()
360 547 : if (oSetDSName.find(poSimpleSource->m_osSrcDSName) !=
361 1094 : oSetDSName.end())
362 : {
363 0 : bRet = false;
364 4 : break;
365 : }
366 547 : oSetDSName.insert(poSimpleSource->m_osSrcDSName);
367 :
368 : double dfSourceXOff;
369 : double dfSourceYOff;
370 : double dfSourceXSize;
371 : double dfSourceYSize;
372 547 : poSimpleSource->GetDstWindow(dfSourceXOff, dfSourceYOff,
373 : dfSourceXSize, dfSourceYSize);
374 547 : constexpr double EPSILON = 1e-1;
375 : // We floor/ceil to detect potential overlaps of sources whose
376 : // destination offset is not integer. Otherwise, in case of
377 : // multithreading, this could cause one line to be written
378 : // concurrently by multiple threads.
379 547 : sSourceBounds.minx = std::floor(dfSourceXOff) + EPSILON;
380 547 : sSourceBounds.miny = std::floor(dfSourceYOff) + EPSILON;
381 547 : sSourceBounds.maxx =
382 547 : std::ceil(dfSourceXOff + dfSourceXSize) - EPSILON;
383 547 : sSourceBounds.maxy =
384 547 : std::ceil(dfSourceYOff + dfSourceYSize) - EPSILON;
385 547 : iLastSource = iSource;
386 :
387 547 : if (hQuadTree)
388 : {
389 : // Check that the new source doesn't overlap an existing one.
390 280 : if (CPLQuadTreeHasMatch(hQuadTree, &sSourceBounds))
391 : {
392 4 : bRet = false;
393 4 : break;
394 : }
395 :
396 276 : CPLQuadTreeInsertWithBounds(
397 : hQuadTree,
398 276 : reinterpret_cast<void *>(static_cast<uintptr_t>(iSource)),
399 : &sSourceBounds);
400 : }
401 :
402 543 : ++nContributingSources;
403 : }
404 : }
405 :
406 513 : if (hQuadTree)
407 28 : CPLQuadTreeDestroy(hQuadTree);
408 :
409 1026 : return bRet;
410 : }
411 :
412 : /************************************************************************/
413 : /* VRTSourcedRasterBandRasterIOJob */
414 : /************************************************************************/
415 :
416 : /** Structure used to declare a threaded job to satisfy IRasterIO()
417 : * on a given source.
418 : */
419 : struct VRTSourcedRasterBandRasterIOJob
420 : {
421 : std::atomic<int> *pnCompletedJobs = nullptr;
422 : std::atomic<bool> *pbSuccess = nullptr;
423 : VRTDataset::QueueWorkingStates *poQueueWorkingStates = nullptr;
424 : CPLErrorAccumulator *poErrorAccumulator = nullptr;
425 :
426 : GDALDataType eVRTBandDataType = GDT_Unknown;
427 : int nXOff = 0;
428 : int nYOff = 0;
429 : int nXSize = 0;
430 : int nYSize = 0;
431 : void *pData = nullptr;
432 : int nBufXSize = 0;
433 : int nBufYSize = 0;
434 : GDALDataType eBufType = GDT_Unknown;
435 : GSpacing nPixelSpace = 0;
436 : GSpacing nLineSpace = 0;
437 : GDALRasterIOExtraArg sExtraArg{};
438 : VRTSimpleSource *poSource = nullptr;
439 :
440 : static void Func(void *pData);
441 : };
442 :
443 : /************************************************************************/
444 : /* VRTSourcedRasterBandRasterIOJob::Func() */
445 : /************************************************************************/
446 :
447 99 : void VRTSourcedRasterBandRasterIOJob::Func(void *pData)
448 : {
449 : auto psJob = std::unique_ptr<VRTSourcedRasterBandRasterIOJob>(
450 198 : static_cast<VRTSourcedRasterBandRasterIOJob *>(pData));
451 99 : if (*psJob->pbSuccess)
452 : {
453 99 : psJob->sExtraArg.pfnProgress = nullptr;
454 99 : psJob->sExtraArg.pProgressData = nullptr;
455 :
456 99 : std::unique_ptr<VRTSource::WorkingState> poWorkingState;
457 : {
458 198 : std::lock_guard oLock(psJob->poQueueWorkingStates->oMutex);
459 : poWorkingState =
460 99 : std::move(psJob->poQueueWorkingStates->oStates.back());
461 99 : psJob->poQueueWorkingStates->oStates.pop_back();
462 99 : CPLAssert(poWorkingState.get());
463 : }
464 :
465 198 : auto oAccumulator = psJob->poErrorAccumulator->InstallForCurrentScope();
466 99 : CPL_IGNORE_RET_VAL(oAccumulator);
467 :
468 198 : if (psJob->poSource->RasterIO(
469 99 : psJob->eVRTBandDataType, psJob->nXOff, psJob->nYOff,
470 99 : psJob->nXSize, psJob->nYSize, psJob->pData, psJob->nBufXSize,
471 99 : psJob->nBufYSize, psJob->eBufType, psJob->nPixelSpace,
472 99 : psJob->nLineSpace, &psJob->sExtraArg,
473 198 : *(poWorkingState.get())) != CE_None)
474 : {
475 0 : *psJob->pbSuccess = false;
476 : }
477 :
478 : {
479 198 : std::lock_guard oLock(psJob->poQueueWorkingStates->oMutex);
480 198 : psJob->poQueueWorkingStates->oStates.push_back(
481 99 : std::move(poWorkingState));
482 : }
483 : }
484 :
485 99 : ++(*psJob->pnCompletedJobs);
486 99 : }
487 :
488 : /************************************************************************/
489 : /* MayMultiBlockReadingBeMultiThreaded() */
490 : /************************************************************************/
491 :
492 0 : bool VRTSourcedRasterBand::MayMultiBlockReadingBeMultiThreaded() const
493 : {
494 : GDALRasterIOExtraArg sExtraArg;
495 0 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
496 0 : sExtraArg.dfXOff = 0;
497 0 : sExtraArg.dfYOff = 0;
498 0 : sExtraArg.dfXSize = nRasterXSize;
499 0 : sExtraArg.dfYSize = nRasterYSize;
500 0 : int nContributingSources = 0;
501 0 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
502 0 : return l_poDS &&
503 0 : CanIRasterIOBeForwardedToEachSource(GF_Read, 0, 0, nRasterXSize,
504 0 : nRasterYSize, nRasterXSize,
505 0 : nRasterYSize, &sExtraArg) &&
506 0 : CanMultiThreadRasterIO(0, 0, nRasterXSize, nRasterYSize,
507 0 : nContributingSources) &&
508 0 : nContributingSources > 1 && VRTDataset::GetNumThreads(l_poDS) > 1;
509 : }
510 :
511 : /************************************************************************/
512 : /* VRTSourcedRasterBand::InitializeOutputBuffer() */
513 : /************************************************************************/
514 :
515 61729 : void VRTSourcedRasterBand::InitializeOutputBuffer(void *pData, int nBufXSize,
516 : int nBufYSize,
517 : GDALDataType eBufType,
518 : GSpacing nPixelSpace,
519 : GSpacing nLineSpace) const
520 : {
521 61729 : if (nPixelSpace == GDALGetDataTypeSizeBytes(eBufType) &&
522 29966 : !(m_bNoDataValueSet && m_dfNoDataValue != 0.0) &&
523 121397 : !(m_bNoDataSetAsInt64 && m_nNoDataValueInt64 != 0) &&
524 29702 : !(m_bNoDataSetAsUInt64 && m_nNoDataValueUInt64 != 0))
525 : {
526 29701 : if (nLineSpace == nBufXSize * nPixelSpace)
527 : {
528 28621 : memset(pData, 0, static_cast<size_t>(nBufYSize * nLineSpace));
529 : }
530 : else
531 : {
532 13861 : for (int iLine = 0; iLine < nBufYSize; iLine++)
533 : {
534 12781 : memset(static_cast<GByte *>(pData) +
535 12781 : static_cast<GIntBig>(iLine) * nLineSpace,
536 12781 : 0, static_cast<size_t>(nBufXSize * nPixelSpace));
537 : }
538 : }
539 : }
540 32028 : else if (m_bNoDataSetAsInt64)
541 : {
542 7 : for (int iLine = 0; iLine < nBufYSize; iLine++)
543 : {
544 6 : GDALCopyWords(&m_nNoDataValueInt64, GDT_Int64, 0,
545 : static_cast<GByte *>(pData) +
546 6 : static_cast<GIntBig>(nLineSpace) * iLine,
547 : eBufType, static_cast<int>(nPixelSpace), nBufXSize);
548 : }
549 : }
550 32027 : else if (m_bNoDataSetAsUInt64)
551 : {
552 7 : for (int iLine = 0; iLine < nBufYSize; iLine++)
553 : {
554 6 : GDALCopyWords(&m_nNoDataValueUInt64, GDT_UInt64, 0,
555 : static_cast<GByte *>(pData) +
556 6 : static_cast<GIntBig>(nLineSpace) * iLine,
557 : eBufType, static_cast<int>(nPixelSpace), nBufXSize);
558 : }
559 : }
560 : else
561 : {
562 32026 : double dfWriteValue = 0.0;
563 32026 : if (m_bNoDataValueSet)
564 290 : dfWriteValue = m_dfNoDataValue;
565 :
566 131674 : for (int iLine = 0; iLine < nBufYSize; iLine++)
567 : {
568 99648 : GDALCopyWords(&dfWriteValue, GDT_Float64, 0,
569 : static_cast<GByte *>(pData) +
570 99648 : static_cast<GIntBig>(nLineSpace) * iLine,
571 : eBufType, static_cast<int>(nPixelSpace), nBufXSize);
572 : }
573 : }
574 61729 : }
575 :
576 : /************************************************************************/
577 : /* IRasterIO() */
578 : /************************************************************************/
579 :
580 289889 : CPLErr VRTSourcedRasterBand::IRasterIO(
581 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
582 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
583 : GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
584 :
585 : {
586 289889 : if (eRWFlag == GF_Write)
587 : {
588 0 : CPLError(CE_Failure, CPLE_AppDefined,
589 : "Writing through VRTSourcedRasterBand is not supported.");
590 0 : return CE_Failure;
591 : }
592 :
593 579778 : const std::string osFctId("VRTSourcedRasterBand::IRasterIO");
594 579778 : GDALAntiRecursionGuard oGuard(osFctId);
595 289889 : if (oGuard.GetCallDepth() >= 32)
596 : {
597 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
598 0 : return CE_Failure;
599 : }
600 :
601 869667 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
602 : // Allow 2 recursion depths on the same dataset for non-nearest resampling
603 289889 : if (oGuard2.GetCallDepth() > 2)
604 : {
605 1 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
606 1 : return CE_Failure;
607 : }
608 :
609 : /* ==================================================================== */
610 : /* Do we have overviews that would be appropriate to satisfy */
611 : /* this request? */
612 : /* ==================================================================== */
613 289888 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
614 289888 : if (l_poDS &&
615 579417 : l_poDS->m_apoOverviews.empty() && // do not use virtual overviews
616 579776 : (nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() > 0)
617 : {
618 1 : if (OverviewRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
619 : nBufXSize, nBufYSize, eBufType, nPixelSpace,
620 1 : nLineSpace, psExtraArg) == CE_None)
621 0 : return CE_None;
622 : }
623 :
624 : // If resampling with non-nearest neighbour, we need to be careful
625 : // if the VRT band exposes a nodata value, but the sources do not have it.
626 : // To also avoid edge effects on sources when downsampling, use the
627 : // base implementation of IRasterIO() (that is acquiring sources at their
628 : // nominal resolution, and then downsampling), but only if none of the
629 : // contributing sources have overviews.
630 289888 : if (l_poDS && !CanIRasterIOBeForwardedToEachSource(
631 : eRWFlag, nXOff, nYOff, nXSize, nYSize, nBufXSize,
632 : nBufYSize, psExtraArg))
633 : {
634 14 : const bool bBackupEnabledOverviews = l_poDS->AreOverviewsEnabled();
635 14 : if (!l_poDS->m_apoOverviews.empty() && l_poDS->AreOverviewsEnabled())
636 : {
637 : // Disable use of implicit overviews to avoid infinite
638 : // recursion
639 1 : l_poDS->SetEnableOverviews(false);
640 : }
641 :
642 14 : const auto eResampleAlgBackup = psExtraArg->eResampleAlg;
643 14 : if (psExtraArg->eResampleAlg == GRIORA_NearestNeighbour)
644 : {
645 16 : std::string osResampling;
646 24 : for (auto &poSource : m_papoSources)
647 : {
648 16 : if (poSource->GetType() == VRTComplexSource::GetTypeStatic())
649 : {
650 : auto *const poComplexSource =
651 16 : static_cast<VRTComplexSource *>(poSource.get());
652 16 : if (!poComplexSource->GetResampling().empty())
653 : {
654 16 : if (osResampling.empty())
655 8 : osResampling = poComplexSource->GetResampling();
656 8 : else if (osResampling !=
657 8 : poComplexSource->GetResampling())
658 : {
659 0 : osResampling.clear();
660 0 : break;
661 : }
662 : }
663 : }
664 : }
665 8 : if (!osResampling.empty())
666 8 : psExtraArg->eResampleAlg =
667 8 : GDALRasterIOGetResampleAlg(osResampling.c_str());
668 : }
669 :
670 14 : const auto eErr = GDALRasterBand::IRasterIO(
671 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
672 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
673 :
674 14 : psExtraArg->eResampleAlg = eResampleAlgBackup;
675 14 : l_poDS->SetEnableOverviews(bBackupEnabledOverviews);
676 14 : return eErr;
677 : }
678 :
679 : /* -------------------------------------------------------------------- */
680 : /* Initialize the buffer to some background value. Use the */
681 : /* nodata value if available. */
682 : /* -------------------------------------------------------------------- */
683 289874 : if (!SkipBufferInitialization())
684 : {
685 59121 : InitializeOutputBuffer(pData, nBufXSize, nBufYSize, eBufType,
686 : nPixelSpace, nLineSpace);
687 : }
688 :
689 : /* -------------------------------------------------------------------- */
690 : /* Overlay each source in turn over top this. */
691 : /* -------------------------------------------------------------------- */
692 289874 : CPLErr eErr = CE_None;
693 :
694 289874 : double dfXOff = nXOff;
695 289874 : double dfYOff = nYOff;
696 289874 : double dfXSize = nXSize;
697 289874 : double dfYSize = nYSize;
698 289874 : if (psExtraArg->bFloatingPointWindowValidity)
699 : {
700 311 : dfXOff = psExtraArg->dfXOff;
701 311 : dfYOff = psExtraArg->dfYOff;
702 311 : dfXSize = psExtraArg->dfXSize;
703 311 : dfYSize = psExtraArg->dfYSize;
704 : }
705 :
706 289874 : if (l_poDS)
707 289874 : l_poDS->m_bMultiThreadedRasterIOLastUsed = false;
708 :
709 289874 : int nContributingSources = 0;
710 289874 : int nMaxThreads = 0;
711 289874 : constexpr int MINIMUM_PIXEL_COUNT_FOR_THREADED_IO = 1000 * 1000;
712 579748 : if (l_poDS &&
713 289874 : (static_cast<int64_t>(nBufXSize) * nBufYSize >=
714 289541 : MINIMUM_PIXEL_COUNT_FOR_THREADED_IO ||
715 289541 : static_cast<int64_t>(nXSize) * nYSize >=
716 340 : MINIMUM_PIXEL_COUNT_FOR_THREADED_IO) &&
717 340 : CanMultiThreadRasterIO(dfXOff, dfYOff, dfXSize, dfYSize,
718 337 : nContributingSources) &&
719 579768 : nContributingSources > 1 &&
720 20 : (nMaxThreads = VRTDataset::GetNumThreads(l_poDS)) > 1)
721 : {
722 18 : l_poDS->m_bMultiThreadedRasterIOLastUsed = true;
723 18 : l_poDS->m_oMapSharedSources.InitMutex();
724 :
725 36 : CPLErrorAccumulator errorAccumulator;
726 18 : std::atomic<bool> bSuccess = true;
727 18 : CPLWorkerThreadPool *psThreadPool = GDALGetGlobalThreadPool(
728 18 : std::min(nContributingSources, nMaxThreads));
729 : const int nThreads =
730 18 : std::min(nContributingSources, psThreadPool->GetThreadCount());
731 18 : CPLDebugOnly("VRT",
732 : "IRasterIO(): use optimized "
733 : "multi-threaded code path for mosaic. "
734 : "Using %d threads",
735 : nThreads);
736 :
737 : {
738 36 : std::lock_guard oLock(l_poDS->m_oQueueWorkingStates.oMutex);
739 18 : if (l_poDS->m_oQueueWorkingStates.oStates.size() <
740 18 : static_cast<size_t>(nThreads))
741 : {
742 6 : l_poDS->m_oQueueWorkingStates.oStates.resize(nThreads);
743 : }
744 117 : for (int i = 0; i < nThreads; ++i)
745 : {
746 99 : if (!l_poDS->m_oQueueWorkingStates.oStates[i])
747 75 : l_poDS->m_oQueueWorkingStates.oStates[i] =
748 150 : std::make_unique<VRTSource::WorkingState>();
749 : }
750 : }
751 :
752 18 : auto oQueue = psThreadPool->CreateJobQueue();
753 18 : std::atomic<int> nCompletedJobs = 0;
754 180 : for (auto &poSource : m_papoSources)
755 : {
756 162 : if (!poSource->IsSimpleSource())
757 0 : continue;
758 : auto poSimpleSource =
759 162 : cpl::down_cast<VRTSimpleSource *>(poSource.get());
760 162 : if (poSimpleSource->DstWindowIntersects(dfXOff, dfYOff, dfXSize,
761 : dfYSize))
762 : {
763 99 : auto psJob = new VRTSourcedRasterBandRasterIOJob();
764 99 : psJob->pbSuccess = &bSuccess;
765 99 : psJob->pnCompletedJobs = &nCompletedJobs;
766 99 : psJob->poQueueWorkingStates = &(l_poDS->m_oQueueWorkingStates);
767 99 : psJob->poErrorAccumulator = &errorAccumulator;
768 99 : psJob->eVRTBandDataType = eDataType;
769 99 : psJob->nXOff = nXOff;
770 99 : psJob->nYOff = nYOff;
771 99 : psJob->nXSize = nXSize;
772 99 : psJob->nYSize = nYSize;
773 99 : psJob->pData = pData;
774 99 : psJob->nBufXSize = nBufXSize;
775 99 : psJob->nBufYSize = nBufYSize;
776 99 : psJob->eBufType = eBufType;
777 99 : psJob->nPixelSpace = nPixelSpace;
778 99 : psJob->nLineSpace = nLineSpace;
779 99 : psJob->sExtraArg = *psExtraArg;
780 99 : psJob->poSource = poSimpleSource;
781 :
782 99 : if (!oQueue->SubmitJob(VRTSourcedRasterBandRasterIOJob::Func,
783 : psJob))
784 : {
785 0 : delete psJob;
786 0 : bSuccess = false;
787 0 : break;
788 : }
789 : }
790 : }
791 :
792 73 : while (oQueue->WaitEvent())
793 : {
794 : // Quite rough progress callback. We could do better by counting
795 : // the number of contributing pixels.
796 55 : if (psExtraArg->pfnProgress)
797 : {
798 88 : psExtraArg->pfnProgress(double(nCompletedJobs.load()) /
799 : nContributingSources,
800 : "", psExtraArg->pProgressData);
801 : }
802 : }
803 :
804 18 : errorAccumulator.ReplayErrors();
805 18 : eErr = bSuccess ? CE_None : CE_Failure;
806 : }
807 : else
808 : {
809 289856 : GDALProgressFunc const pfnProgressGlobal = psExtraArg->pfnProgress;
810 289856 : void *const pProgressDataGlobal = psExtraArg->pProgressData;
811 :
812 289856 : VRTSource::WorkingState oWorkingState;
813 289856 : const int nSources = static_cast<int>(m_papoSources.size());
814 564584 : for (int iSource = 0; eErr == CE_None && iSource < nSources; iSource++)
815 : {
816 274728 : psExtraArg->pfnProgress = GDALScaledProgress;
817 549456 : psExtraArg->pProgressData = GDALCreateScaledProgress(
818 274728 : 1.0 * iSource / nSources, 1.0 * (iSource + 1) / nSources,
819 : pfnProgressGlobal, pProgressDataGlobal);
820 274728 : if (psExtraArg->pProgressData == nullptr)
821 274169 : psExtraArg->pfnProgress = nullptr;
822 :
823 549456 : eErr = m_papoSources[iSource]->RasterIO(
824 : eDataType, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
825 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg,
826 274728 : l_poDS ? l_poDS->m_oWorkingState : oWorkingState);
827 :
828 274728 : GDALDestroyScaledProgress(psExtraArg->pProgressData);
829 : }
830 :
831 289856 : psExtraArg->pfnProgress = pfnProgressGlobal;
832 289856 : psExtraArg->pProgressData = pProgressDataGlobal;
833 : }
834 :
835 289874 : if (eErr == CE_None && psExtraArg->pfnProgress)
836 : {
837 426 : psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
838 : }
839 :
840 289874 : return eErr;
841 : }
842 :
843 : /************************************************************************/
844 : /* IGetDataCoverageStatus() */
845 : /************************************************************************/
846 :
847 3040 : int VRTSourcedRasterBand::IGetDataCoverageStatus(int nXOff, int nYOff,
848 : int nXSize, int nYSize,
849 : int nMaskFlagStop,
850 : double *pdfDataPct)
851 : {
852 3040 : if (pdfDataPct)
853 8 : *pdfDataPct = -1.0;
854 :
855 : // Particular case for a single simple source covering the whole dataset
856 6065 : if (m_papoSources.size() == 1 && m_papoSources[0]->IsSimpleSource() &&
857 3025 : m_papoSources[0]->GetType() == VRTSimpleSource::GetTypeStatic())
858 : {
859 : VRTSimpleSource *poSource =
860 2985 : static_cast<VRTSimpleSource *>(m_papoSources[0].get());
861 :
862 2985 : GDALRasterBand *poBand = poSource->GetRasterBand();
863 2985 : if (!poBand)
864 0 : poBand = poSource->GetMaskBandMainBand();
865 2985 : if (!poBand)
866 : {
867 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
868 353 : GDAL_DATA_COVERAGE_STATUS_DATA;
869 : }
870 :
871 : /* Check that it uses the full source dataset */
872 2985 : double dfReqXOff = 0.0;
873 2985 : double dfReqYOff = 0.0;
874 2985 : double dfReqXSize = 0.0;
875 2985 : double dfReqYSize = 0.0;
876 2985 : int nReqXOff = 0;
877 2985 : int nReqYOff = 0;
878 2985 : int nReqXSize = 0;
879 2985 : int nReqYSize = 0;
880 2985 : int nOutXOff = 0;
881 2985 : int nOutYOff = 0;
882 2985 : int nOutXSize = 0;
883 2985 : int nOutYSize = 0;
884 2985 : bool bError = false;
885 2985 : if (poSource->GetSrcDstWindow(
886 2985 : 0, 0, GetXSize(), GetYSize(), GetXSize(), GetYSize(),
887 : GRIORA_NearestNeighbour, &dfReqXOff, &dfReqYOff, &dfReqXSize,
888 : &dfReqYSize, &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
889 2985 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError) &&
890 2985 : nReqXOff == 0 && nReqYOff == 0 && nReqXSize == GetXSize() &&
891 468 : nReqXSize == poBand->GetXSize() && nReqYSize == GetYSize() &&
892 362 : nReqYSize == poBand->GetYSize() && nOutXOff == 0 && nOutYOff == 0 &&
893 5970 : nOutXSize == GetXSize() && nOutYSize == GetYSize())
894 : {
895 353 : return poBand->GetDataCoverageStatus(nXOff, nYOff, nXSize, nYSize,
896 353 : nMaskFlagStop, pdfDataPct);
897 : }
898 : }
899 :
900 : #ifndef HAVE_GEOS
901 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
902 : GDAL_DATA_COVERAGE_STATUS_DATA;
903 : #else
904 2687 : int nStatus = 0;
905 :
906 5374 : auto poPolyNonCoveredBySources = std::make_unique<OGRPolygon>();
907 : {
908 5374 : auto poLR = std::make_unique<OGRLinearRing>();
909 2687 : poLR->addPoint(nXOff, nYOff);
910 2687 : poLR->addPoint(nXOff, nYOff + nYSize);
911 2687 : poLR->addPoint(nXOff + nXSize, nYOff + nYSize);
912 2687 : poLR->addPoint(nXOff + nXSize, nYOff);
913 2687 : poLR->addPoint(nXOff, nYOff);
914 2687 : poPolyNonCoveredBySources->addRingDirectly(poLR.release());
915 : }
916 :
917 2697 : for (auto &poSource : m_papoSources)
918 : {
919 2692 : if (!poSource->IsSimpleSource())
920 : {
921 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
922 2682 : GDAL_DATA_COVERAGE_STATUS_DATA;
923 : }
924 2688 : VRTSimpleSource *poSS = static_cast<VRTSimpleSource *>(poSource.get());
925 : // Check if the AOI is fully inside the source
926 2688 : double dfDstXOff = std::max(0.0, poSS->m_dfDstXOff);
927 2688 : double dfDstYOff = std::max(0.0, poSS->m_dfDstYOff);
928 2688 : double dfDstXSize = poSS->m_dfDstXSize;
929 2688 : double dfDstYSize = poSS->m_dfDstYSize;
930 2688 : auto l_poBand = poSS->GetRasterBand();
931 2688 : if (!l_poBand)
932 0 : continue;
933 2688 : if (dfDstXSize == -1)
934 8 : dfDstXSize = l_poBand->GetXSize() - dfDstXOff;
935 2688 : if (dfDstYSize == -1)
936 8 : dfDstYSize = l_poBand->GetYSize() - dfDstYOff;
937 :
938 2688 : if (nXOff >= dfDstXOff && nYOff >= dfDstYOff &&
939 2679 : nXOff + nXSize <= dfDstXOff + dfDstXSize &&
940 2655 : nYOff + nYSize <= dfDstYOff + dfDstYSize)
941 : {
942 2655 : if (pdfDataPct)
943 1 : *pdfDataPct = 100.0;
944 2655 : return GDAL_DATA_COVERAGE_STATUS_DATA;
945 : }
946 : // Check intersection of bounding boxes.
947 33 : if (dfDstXOff + dfDstXSize > nXOff && dfDstYOff + dfDstYSize > nYOff &&
948 27 : dfDstXOff < nXOff + nXSize && dfDstYOff < nYOff + nYSize)
949 : {
950 25 : nStatus |= GDAL_DATA_COVERAGE_STATUS_DATA;
951 25 : if (poPolyNonCoveredBySources)
952 : {
953 25 : OGRPolygon oPolySource;
954 25 : auto poLR = std::make_unique<OGRLinearRing>();
955 25 : poLR->addPoint(dfDstXOff, dfDstYOff);
956 25 : poLR->addPoint(dfDstXOff, dfDstYOff + dfDstYSize);
957 25 : poLR->addPoint(dfDstXOff + dfDstXSize, dfDstYOff + dfDstYSize);
958 25 : poLR->addPoint(dfDstXOff + dfDstXSize, dfDstYOff);
959 25 : poLR->addPoint(dfDstXOff, dfDstYOff);
960 25 : oPolySource.addRingDirectly(poLR.release());
961 : auto poRes = std::unique_ptr<OGRGeometry>(
962 25 : poPolyNonCoveredBySources->Difference(&oPolySource));
963 25 : if (poRes && poRes->IsEmpty())
964 : {
965 1 : if (pdfDataPct)
966 1 : *pdfDataPct = 100.0;
967 1 : return GDAL_DATA_COVERAGE_STATUS_DATA;
968 : }
969 24 : else if (poRes && poRes->getGeometryType() == wkbPolygon)
970 : {
971 24 : poPolyNonCoveredBySources.reset(
972 : poRes.release()->toPolygon());
973 : }
974 : else
975 : {
976 0 : poPolyNonCoveredBySources.reset();
977 : }
978 : }
979 : }
980 32 : if (nMaskFlagStop != 0 && (nStatus & nMaskFlagStop) != 0)
981 : {
982 22 : return nStatus;
983 : }
984 : }
985 5 : if (poPolyNonCoveredBySources)
986 : {
987 5 : if (!poPolyNonCoveredBySources->IsEmpty())
988 5 : nStatus |= GDAL_DATA_COVERAGE_STATUS_EMPTY;
989 5 : if (pdfDataPct)
990 2 : *pdfDataPct = 100.0 * (1.0 - poPolyNonCoveredBySources->get_Area() /
991 2 : nXSize / nYSize);
992 : }
993 5 : return nStatus;
994 : #endif // HAVE_GEOS
995 : }
996 :
997 : /************************************************************************/
998 : /* IReadBlock() */
999 : /************************************************************************/
1000 :
1001 938 : CPLErr VRTSourcedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
1002 : void *pImage)
1003 :
1004 : {
1005 938 : const int nPixelSize = GDALGetDataTypeSizeBytes(eDataType);
1006 :
1007 938 : int nReadXSize = 0;
1008 938 : if ((nBlockXOff + 1) * nBlockXSize > GetXSize())
1009 16 : nReadXSize = GetXSize() - nBlockXOff * nBlockXSize;
1010 : else
1011 922 : nReadXSize = nBlockXSize;
1012 :
1013 938 : int nReadYSize = 0;
1014 938 : if ((nBlockYOff + 1) * nBlockYSize > GetYSize())
1015 0 : nReadYSize = GetYSize() - nBlockYOff * nBlockYSize;
1016 : else
1017 938 : nReadYSize = nBlockYSize;
1018 :
1019 : GDALRasterIOExtraArg sExtraArg;
1020 938 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1021 :
1022 938 : return IRasterIO(
1023 938 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize, nReadXSize,
1024 : nReadYSize, pImage, nReadXSize, nReadYSize, eDataType, nPixelSize,
1025 1876 : static_cast<GSpacing>(nPixelSize) * nBlockXSize, &sExtraArg);
1026 : }
1027 :
1028 : /************************************************************************/
1029 : /* CPLGettimeofday() */
1030 : /************************************************************************/
1031 :
1032 : #if defined(_WIN32) && !defined(__CYGWIN__)
1033 : #include <sys/timeb.h>
1034 :
1035 : namespace
1036 : {
1037 : struct CPLTimeVal
1038 : {
1039 : time_t tv_sec; /* seconds */
1040 : long tv_usec; /* and microseconds */
1041 : };
1042 : } // namespace
1043 :
1044 : static int CPLGettimeofday(struct CPLTimeVal *tp, void * /* timezonep*/)
1045 : {
1046 : struct _timeb theTime;
1047 :
1048 : _ftime(&theTime);
1049 : tp->tv_sec = static_cast<time_t>(theTime.time);
1050 : tp->tv_usec = theTime.millitm * 1000;
1051 : return 0;
1052 : }
1053 : #else
1054 : #include <sys/time.h> /* for gettimeofday() */
1055 : #define CPLTimeVal timeval
1056 : #define CPLGettimeofday(t, u) gettimeofday(t, u)
1057 : #endif
1058 :
1059 : /************************************************************************/
1060 : /* CanUseSourcesMinMaxImplementations() */
1061 : /************************************************************************/
1062 :
1063 109 : bool VRTSourcedRasterBand::CanUseSourcesMinMaxImplementations()
1064 : {
1065 : const char *pszUseSources =
1066 109 : CPLGetConfigOption("VRT_MIN_MAX_FROM_SOURCES", nullptr);
1067 109 : if (pszUseSources)
1068 0 : return CPLTestBool(pszUseSources);
1069 :
1070 : // Use heuristics to determine if we are going to use the source
1071 : // GetMinimum() or GetMaximum() implementation: all the sources must be
1072 : // "simple" sources with a dataset description that match a "regular" file
1073 : // on the filesystem, whose open time and GetMinimum()/GetMaximum()
1074 : // implementations we hope to be fast enough.
1075 : // In case of doubt return FALSE.
1076 : struct CPLTimeVal tvStart;
1077 109 : memset(&tvStart, 0, sizeof(CPLTimeVal));
1078 109 : if (m_papoSources.size() > 1)
1079 10 : CPLGettimeofday(&tvStart, nullptr);
1080 237 : for (auto &poSource : m_papoSources)
1081 : {
1082 131 : if (!(poSource->IsSimpleSource()))
1083 3 : return false;
1084 : VRTSimpleSource *const poSimpleSource =
1085 129 : static_cast<VRTSimpleSource *>(poSource.get());
1086 129 : const char *pszFilename = poSimpleSource->m_osSrcDSName.c_str();
1087 : // /vsimem/ should be fast.
1088 129 : if (STARTS_WITH(pszFilename, "/vsimem/"))
1089 10 : continue;
1090 : // but not other /vsi filesystems
1091 119 : if (STARTS_WITH(pszFilename, "/vsi"))
1092 0 : return false;
1093 119 : char ch = '\0';
1094 : // We will assume that filenames that are only with ascii characters
1095 : // are real filenames and so we will not try to 'stat' them.
1096 1789 : for (int i = 0; (ch = pszFilename[i]) != '\0'; i++)
1097 : {
1098 1702 : if (!((ch >= 'a' && ch <= 'z') || (ch >= 'A' && ch <= 'Z') ||
1099 399 : (ch >= '0' && ch <= '9') || ch == ':' || ch == '/' ||
1100 137 : ch == '\\' || ch == ' ' || ch == '.' || ch == '_'))
1101 13 : break;
1102 : }
1103 119 : if (ch != '\0')
1104 : {
1105 : // Otherwise do a real filesystem check.
1106 : VSIStatBuf sStat;
1107 13 : if (VSIStat(pszFilename, &sStat) != 0)
1108 1 : return false;
1109 12 : if (m_papoSources.size() > 1)
1110 : {
1111 : struct CPLTimeVal tvCur;
1112 10 : CPLGettimeofday(&tvCur, nullptr);
1113 10 : if (tvCur.tv_sec - tvStart.tv_sec +
1114 10 : (tvCur.tv_usec - tvStart.tv_usec) * 1e-6 >
1115 : 1)
1116 0 : return false;
1117 : }
1118 : }
1119 : }
1120 106 : return true;
1121 : }
1122 :
1123 : /************************************************************************/
1124 : /* GetMinimum() */
1125 : /************************************************************************/
1126 :
1127 60 : double VRTSourcedRasterBand::GetMinimum(int *pbSuccess)
1128 : {
1129 60 : const char *const pszValue = GetMetadataItem("STATISTICS_MINIMUM");
1130 60 : if (pszValue != nullptr)
1131 : {
1132 5 : if (pbSuccess != nullptr)
1133 5 : *pbSuccess = TRUE;
1134 :
1135 5 : return CPLAtofM(pszValue);
1136 : }
1137 :
1138 55 : if (!CanUseSourcesMinMaxImplementations())
1139 2 : return GDALRasterBand::GetMinimum(pbSuccess);
1140 :
1141 106 : const std::string osFctId("VRTSourcedRasterBand::GetMinimum");
1142 106 : GDALAntiRecursionGuard oGuard(osFctId);
1143 53 : if (oGuard.GetCallDepth() >= 32)
1144 : {
1145 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1146 0 : if (pbSuccess != nullptr)
1147 0 : *pbSuccess = FALSE;
1148 0 : return 0;
1149 : }
1150 :
1151 159 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1152 53 : if (oGuard2.GetCallDepth() >= 2)
1153 : {
1154 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1155 0 : if (pbSuccess != nullptr)
1156 0 : *pbSuccess = FALSE;
1157 0 : return 0;
1158 : }
1159 :
1160 : struct CPLTimeVal tvStart;
1161 53 : memset(&tvStart, 0, sizeof(CPLTimeVal));
1162 53 : if (m_papoSources.size() > 1)
1163 5 : CPLGettimeofday(&tvStart, nullptr);
1164 53 : double dfMin = 0;
1165 58 : for (size_t iSource = 0; iSource < m_papoSources.size(); iSource++)
1166 : {
1167 53 : int bSuccess = FALSE;
1168 106 : double dfSourceMin = m_papoSources[iSource]->GetMinimum(
1169 53 : GetXSize(), GetYSize(), &bSuccess);
1170 53 : if (!bSuccess)
1171 : {
1172 48 : dfMin = GDALRasterBand::GetMinimum(pbSuccess);
1173 48 : return dfMin;
1174 : }
1175 :
1176 5 : if (iSource == 0 || dfSourceMin < dfMin)
1177 : {
1178 4 : dfMin = dfSourceMin;
1179 4 : if (dfMin == 0 && eDataType == GDT_UInt8)
1180 0 : break;
1181 : }
1182 5 : if (m_papoSources.size() > 1)
1183 : {
1184 : struct CPLTimeVal tvCur;
1185 2 : CPLGettimeofday(&tvCur, nullptr);
1186 2 : if (tvCur.tv_sec - tvStart.tv_sec +
1187 2 : (tvCur.tv_usec - tvStart.tv_usec) * 1e-6 >
1188 : 1)
1189 : {
1190 0 : return GDALRasterBand::GetMinimum(pbSuccess);
1191 : }
1192 : }
1193 : }
1194 :
1195 5 : if (pbSuccess != nullptr)
1196 5 : *pbSuccess = TRUE;
1197 :
1198 5 : return dfMin;
1199 : }
1200 :
1201 : /************************************************************************/
1202 : /* GetMaximum() */
1203 : /************************************************************************/
1204 :
1205 58 : double VRTSourcedRasterBand::GetMaximum(int *pbSuccess)
1206 : {
1207 58 : const char *const pszValue = GetMetadataItem("STATISTICS_MAXIMUM");
1208 58 : if (pszValue != nullptr)
1209 : {
1210 4 : if (pbSuccess != nullptr)
1211 4 : *pbSuccess = TRUE;
1212 :
1213 4 : return CPLAtofM(pszValue);
1214 : }
1215 :
1216 54 : if (!CanUseSourcesMinMaxImplementations())
1217 1 : return GDALRasterBand::GetMaximum(pbSuccess);
1218 :
1219 106 : const std::string osFctId("VRTSourcedRasterBand::GetMaximum");
1220 106 : GDALAntiRecursionGuard oGuard(osFctId);
1221 53 : if (oGuard.GetCallDepth() >= 32)
1222 : {
1223 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1224 0 : if (pbSuccess != nullptr)
1225 0 : *pbSuccess = FALSE;
1226 0 : return 0;
1227 : }
1228 :
1229 159 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1230 53 : if (oGuard2.GetCallDepth() >= 2)
1231 : {
1232 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1233 0 : if (pbSuccess != nullptr)
1234 0 : *pbSuccess = FALSE;
1235 0 : return 0;
1236 : }
1237 :
1238 : struct CPLTimeVal tvStart;
1239 53 : memset(&tvStart, 0, sizeof(CPLTimeVal));
1240 53 : if (m_papoSources.size() > 1)
1241 5 : CPLGettimeofday(&tvStart, nullptr);
1242 53 : double dfMax = 0;
1243 55 : for (size_t iSource = 0; iSource < m_papoSources.size(); iSource++)
1244 : {
1245 52 : int bSuccess = FALSE;
1246 104 : const double dfSourceMax = m_papoSources[iSource]->GetMaximum(
1247 52 : GetXSize(), GetYSize(), &bSuccess);
1248 52 : if (!bSuccess)
1249 : {
1250 48 : dfMax = GDALRasterBand::GetMaximum(pbSuccess);
1251 48 : return dfMax;
1252 : }
1253 :
1254 4 : if (iSource == 0 || dfSourceMax > dfMax)
1255 : {
1256 4 : dfMax = dfSourceMax;
1257 4 : if (dfMax == 255.0 && eDataType == GDT_UInt8)
1258 2 : break;
1259 : }
1260 2 : if (m_papoSources.size() > 1)
1261 : {
1262 : struct CPLTimeVal tvCur;
1263 0 : CPLGettimeofday(&tvCur, nullptr);
1264 0 : if (tvCur.tv_sec - tvStart.tv_sec +
1265 0 : (tvCur.tv_usec - tvStart.tv_usec) * 1e-6 >
1266 : 1)
1267 : {
1268 0 : return GDALRasterBand::GetMaximum(pbSuccess);
1269 : }
1270 : }
1271 : }
1272 :
1273 5 : if (pbSuccess != nullptr)
1274 5 : *pbSuccess = TRUE;
1275 :
1276 5 : return dfMax;
1277 : }
1278 :
1279 : /************************************************************************/
1280 : /*IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange() */
1281 : /************************************************************************/
1282 :
1283 : /* Returns true if the VRT raster band consists of non-overlapping simple
1284 : * sources or complex sources that don't change values, and use the full extent
1285 : * of the source band.
1286 : */
1287 56 : bool VRTSourcedRasterBand::
1288 : IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1289 : bool bAllowMaxValAdjustment) const
1290 : {
1291 56 : bool bRet = true;
1292 : CPLRectObj sGlobalBounds;
1293 56 : sGlobalBounds.minx = 0;
1294 56 : sGlobalBounds.miny = 0;
1295 56 : sGlobalBounds.maxx = nRasterXSize;
1296 56 : sGlobalBounds.maxy = nRasterYSize;
1297 56 : CPLQuadTree *hQuadTree = CPLQuadTreeCreate(&sGlobalBounds, nullptr);
1298 130 : for (int i = 0; i < static_cast<int>(m_papoSources.size()); ++i)
1299 : {
1300 83 : if (!m_papoSources[i]->IsSimpleSource())
1301 : {
1302 1 : bRet = false;
1303 9 : break;
1304 : }
1305 :
1306 : auto poSimpleSource =
1307 82 : cpl::down_cast<VRTSimpleSource *>(m_papoSources[i].get());
1308 82 : const char *pszType = poSimpleSource->GetType();
1309 82 : if (pszType == VRTSimpleSource::GetTypeStatic())
1310 : {
1311 : // ok
1312 : }
1313 5 : else if (pszType == VRTComplexSource::GetTypeStatic())
1314 : {
1315 : auto poComplexSource =
1316 5 : cpl::down_cast<VRTComplexSource *>(poSimpleSource);
1317 5 : if (!poComplexSource->AreValuesUnchanged())
1318 : {
1319 2 : bRet = false;
1320 2 : break;
1321 : }
1322 : }
1323 : else
1324 : {
1325 0 : bRet = false;
1326 0 : break;
1327 : }
1328 :
1329 80 : if (!bAllowMaxValAdjustment && poSimpleSource->NeedMaxValAdjustment())
1330 : {
1331 2 : bRet = false;
1332 2 : break;
1333 : }
1334 :
1335 78 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1336 156 : if (poSimpleSourceBand == nullptr ||
1337 78 : poSimpleSourceBand->GetRasterDataType() != eDataType)
1338 : {
1339 0 : bRet = false;
1340 0 : break;
1341 : }
1342 :
1343 78 : double dfReqXOff = 0.0;
1344 78 : double dfReqYOff = 0.0;
1345 78 : double dfReqXSize = 0.0;
1346 78 : double dfReqYSize = 0.0;
1347 78 : int nReqXOff = 0;
1348 78 : int nReqYOff = 0;
1349 78 : int nReqXSize = 0;
1350 78 : int nReqYSize = 0;
1351 78 : int nOutXOff = 0;
1352 78 : int nOutYOff = 0;
1353 78 : int nOutXSize = 0;
1354 78 : int nOutYSize = 0;
1355 :
1356 78 : bool bError = false;
1357 234 : if (!poSimpleSource->GetSrcDstWindow(
1358 78 : 0, 0, nRasterXSize, nRasterYSize, nRasterXSize, nRasterYSize,
1359 : GRIORA_NearestNeighbour, &dfReqXOff, &dfReqYOff, &dfReqXSize,
1360 : &dfReqYSize, &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
1361 78 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError) ||
1362 78 : nReqXOff != 0 || nReqYOff != 0 ||
1363 77 : nReqXSize != poSimpleSourceBand->GetXSize() ||
1364 77 : nReqYSize != poSimpleSourceBand->GetYSize() ||
1365 156 : nOutXSize != nReqXSize || nOutYSize != nReqYSize)
1366 : {
1367 1 : bRet = false;
1368 1 : break;
1369 : }
1370 :
1371 : CPLRectObj sBounds;
1372 77 : constexpr double EPSILON = 1e-1;
1373 77 : sBounds.minx = nOutXOff + EPSILON;
1374 77 : sBounds.miny = nOutYOff + EPSILON;
1375 77 : sBounds.maxx = nOutXOff + nOutXSize - EPSILON;
1376 77 : sBounds.maxy = nOutYOff + nOutYSize - EPSILON;
1377 :
1378 : // Check that the new source doesn't overlap an existing one.
1379 77 : if (CPLQuadTreeHasMatch(hQuadTree, &sBounds))
1380 : {
1381 3 : bRet = false;
1382 3 : break;
1383 : }
1384 :
1385 74 : CPLQuadTreeInsertWithBounds(
1386 74 : hQuadTree, reinterpret_cast<void *>(static_cast<uintptr_t>(i)),
1387 : &sBounds);
1388 : }
1389 56 : CPLQuadTreeDestroy(hQuadTree);
1390 :
1391 56 : return bRet;
1392 : }
1393 :
1394 : /************************************************************************/
1395 : /* ComputeRasterMinMax() */
1396 : /************************************************************************/
1397 :
1398 23 : CPLErr VRTSourcedRasterBand::ComputeRasterMinMax(int bApproxOK,
1399 : double *adfMinMax)
1400 : {
1401 : /* -------------------------------------------------------------------- */
1402 : /* Does the driver already know the min/max? */
1403 : /* -------------------------------------------------------------------- */
1404 23 : if (bApproxOK)
1405 : {
1406 3 : int bSuccessMin = FALSE;
1407 3 : int bSuccessMax = FALSE;
1408 :
1409 3 : const double dfMin = GetMinimum(&bSuccessMin);
1410 3 : const double dfMax = GetMaximum(&bSuccessMax);
1411 :
1412 3 : if (bSuccessMin && bSuccessMax)
1413 : {
1414 1 : adfMinMax[0] = dfMin;
1415 1 : adfMinMax[1] = dfMax;
1416 1 : return CE_None;
1417 : }
1418 : }
1419 :
1420 44 : const std::string osFctId("VRTSourcedRasterBand::ComputeRasterMinMax");
1421 44 : GDALAntiRecursionGuard oGuard(osFctId);
1422 22 : if (oGuard.GetCallDepth() >= 32)
1423 : {
1424 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1425 0 : return CE_Failure;
1426 : }
1427 :
1428 66 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1429 22 : if (oGuard2.GetCallDepth() >= 2)
1430 : {
1431 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1432 0 : return CE_Failure;
1433 : }
1434 :
1435 : /* -------------------------------------------------------------------- */
1436 : /* If we have overview bands, use them for min/max. */
1437 : /* -------------------------------------------------------------------- */
1438 22 : if (bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews())
1439 : {
1440 : GDALRasterBand *const poBand =
1441 0 : GetRasterSampleOverview(GDALSTAT_APPROX_NUMSAMPLES);
1442 :
1443 0 : if (poBand != nullptr && poBand != this)
1444 : {
1445 0 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
1446 0 : if (l_poDS && !l_poDS->m_apoOverviews.empty() &&
1447 0 : dynamic_cast<VRTSourcedRasterBand *>(poBand) != nullptr)
1448 : {
1449 0 : auto apoTmpOverviews = std::move(l_poDS->m_apoOverviews);
1450 0 : l_poDS->m_apoOverviews.clear();
1451 0 : auto eErr = poBand->GDALRasterBand::ComputeRasterMinMax(
1452 : TRUE, adfMinMax);
1453 0 : l_poDS->m_apoOverviews = std::move(apoTmpOverviews);
1454 0 : return eErr;
1455 : }
1456 : else
1457 : {
1458 0 : return poBand->ComputeRasterMinMax(TRUE, adfMinMax);
1459 : }
1460 : }
1461 : }
1462 :
1463 22 : if (IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1464 : /*bAllowMaxValAdjustment = */ true))
1465 : {
1466 19 : CPLDebugOnly(
1467 : "VRT", "ComputeRasterMinMax(): use optimized code path for mosaic");
1468 :
1469 19 : uint64_t nCoveredArea = 0;
1470 :
1471 : // If source bands have nodata value, we can't use source band's
1472 : // ComputeRasterMinMax() as we don't know if there are pixels actually
1473 : // at the nodata value, so use ComputeStatistics() instead that takes
1474 : // into account that aspect.
1475 19 : bool bUseComputeStatistics = false;
1476 38 : for (auto &poSource : m_papoSources)
1477 : {
1478 : auto poSimpleSource =
1479 24 : cpl::down_cast<VRTSimpleSource *>(poSource.get());
1480 24 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1481 : // already checked by IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange
1482 24 : assert(poSimpleSourceBand);
1483 24 : int bHasNoData = FALSE;
1484 24 : CPL_IGNORE_RET_VAL(poSimpleSourceBand->GetNoDataValue(&bHasNoData));
1485 24 : if (bHasNoData)
1486 : {
1487 5 : bUseComputeStatistics = true;
1488 5 : break;
1489 : }
1490 19 : nCoveredArea +=
1491 19 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
1492 19 : poSimpleSourceBand->GetYSize();
1493 : }
1494 :
1495 19 : if (bUseComputeStatistics)
1496 : {
1497 : CPLErr eErr;
1498 5 : std::string osLastErrorMsg;
1499 : {
1500 10 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1501 5 : CPLErrorReset();
1502 10 : eErr = ComputeStatistics(bApproxOK, &adfMinMax[0],
1503 : &adfMinMax[1], nullptr, nullptr,
1504 5 : nullptr, nullptr, nullptr);
1505 5 : if (eErr == CE_Failure)
1506 : {
1507 2 : osLastErrorMsg = CPLGetLastErrorMsg();
1508 : }
1509 : }
1510 5 : if (eErr == CE_Failure)
1511 : {
1512 2 : if (strstr(osLastErrorMsg.c_str(), "no valid pixels found") !=
1513 : nullptr)
1514 : {
1515 2 : ReportError(CE_Failure, CPLE_AppDefined,
1516 : "Failed to compute min/max, no valid pixels "
1517 : "found in sampling.");
1518 : }
1519 : else
1520 : {
1521 0 : ReportError(CE_Failure, CPLE_AppDefined, "%s",
1522 : osLastErrorMsg.c_str());
1523 : }
1524 : }
1525 5 : return eErr;
1526 : }
1527 :
1528 14 : bool bSignedByte = false;
1529 14 : if (eDataType == GDT_UInt8)
1530 : {
1531 8 : EnablePixelTypeSignedByteWarning(false);
1532 : const char *pszPixelType =
1533 8 : GetMetadataItem("PIXELTYPE", GDAL_MDD_IMAGE_STRUCTURE);
1534 8 : EnablePixelTypeSignedByteWarning(true);
1535 8 : bSignedByte =
1536 8 : pszPixelType != nullptr && EQUAL(pszPixelType, "SIGNEDBYTE");
1537 : }
1538 :
1539 14 : double dfGlobalMin = std::numeric_limits<double>::max();
1540 14 : double dfGlobalMax = -std::numeric_limits<double>::max();
1541 :
1542 : // If the mosaic doesn't cover the whole VRT raster, take into account
1543 : // VRT nodata value.
1544 14 : if (nCoveredArea < static_cast<uint64_t>(nRasterXSize) * nRasterYSize)
1545 : {
1546 2 : if (m_bNoDataValueSet && m_bHideNoDataValue)
1547 : {
1548 1 : if (IsNoDataValueInDataTypeRange())
1549 : {
1550 1 : dfGlobalMin = std::min(dfGlobalMin, m_dfNoDataValue);
1551 1 : dfGlobalMax = std::max(dfGlobalMax, m_dfNoDataValue);
1552 : }
1553 : }
1554 1 : else if (!m_bNoDataValueSet)
1555 : {
1556 0 : dfGlobalMin = std::min(dfGlobalMin, 0.0);
1557 0 : dfGlobalMax = std::max(dfGlobalMax, 0.0);
1558 : }
1559 : }
1560 :
1561 31 : for (auto &poSource : m_papoSources)
1562 : {
1563 : auto poSimpleSource =
1564 18 : cpl::down_cast<VRTSimpleSource *>(poSource.get());
1565 18 : double adfMinMaxSource[2] = {0};
1566 :
1567 18 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1568 : // already checked by IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange
1569 18 : assert(poSimpleSourceBand);
1570 36 : CPLErr eErr = poSimpleSourceBand->ComputeRasterMinMax(
1571 18 : bApproxOK, adfMinMaxSource);
1572 18 : if (eErr == CE_Failure)
1573 : {
1574 1 : return CE_Failure;
1575 : }
1576 : else
1577 : {
1578 18 : if (poSimpleSource->NeedMaxValAdjustment())
1579 : {
1580 2 : const double dfMaxValue =
1581 2 : static_cast<double>(poSimpleSource->m_nMaxValue);
1582 2 : adfMinMaxSource[0] =
1583 2 : std::min(adfMinMaxSource[0], dfMaxValue);
1584 2 : adfMinMaxSource[1] =
1585 2 : std::min(adfMinMaxSource[1], dfMaxValue);
1586 : }
1587 :
1588 18 : if (m_bNoDataValueSet && !m_bHideNoDataValue &&
1589 3 : m_dfNoDataValue >= adfMinMaxSource[0] &&
1590 1 : m_dfNoDataValue <= adfMinMaxSource[1])
1591 : {
1592 1 : return GDALRasterBand::ComputeRasterMinMax(bApproxOK,
1593 1 : adfMinMax);
1594 : }
1595 :
1596 17 : dfGlobalMin = std::min(dfGlobalMin, adfMinMaxSource[0]);
1597 17 : dfGlobalMax = std::max(dfGlobalMax, adfMinMaxSource[1]);
1598 : }
1599 :
1600 : // Early exit if we know we reached theoretical bounds
1601 17 : if (eDataType == GDT_UInt8 && !bSignedByte && dfGlobalMin == 0.0 &&
1602 0 : dfGlobalMax == 255.0)
1603 : {
1604 0 : break;
1605 : }
1606 : }
1607 :
1608 13 : if (dfGlobalMin > dfGlobalMax)
1609 : {
1610 0 : adfMinMax[0] = 0.0;
1611 0 : adfMinMax[1] = 0.0;
1612 0 : ReportError(CE_Failure, CPLE_AppDefined,
1613 : "Failed to compute min/max, no valid pixels found in "
1614 : "sampling.");
1615 0 : return CE_Failure;
1616 : }
1617 :
1618 13 : adfMinMax[0] = dfGlobalMin;
1619 13 : adfMinMax[1] = dfGlobalMax;
1620 13 : return CE_None;
1621 : }
1622 : else
1623 : {
1624 3 : return GDALRasterBand::ComputeRasterMinMax(bApproxOK, adfMinMax);
1625 : }
1626 : }
1627 :
1628 : // #define naive_update_not_used
1629 :
1630 : namespace
1631 : {
1632 : struct Context
1633 : {
1634 : CPL_DISALLOW_COPY_ASSIGN(Context)
1635 : Context() = default;
1636 :
1637 : // Protected by mutex
1638 : std::mutex oMutex{};
1639 : uint64_t nTotalIteratedPixels = 0;
1640 : uint64_t nLastReportedPixels = 0;
1641 : bool bFailure = false;
1642 : bool bFallbackToBase = false;
1643 : // End of protected by mutex
1644 :
1645 : int nSources = 0;
1646 :
1647 : bool bApproxOK = false;
1648 : GDALProgressFunc pfnProgress = nullptr;
1649 : void *pProgressData = nullptr;
1650 :
1651 : // VRTSourcedRasterBand parameters
1652 : double dfNoDataValue = 0;
1653 : bool bNoDataValueSet = false;
1654 : bool bHideNoDataValue = false;
1655 :
1656 : double dfGlobalMin = std::numeric_limits<double>::max();
1657 : double dfGlobalMax = -std::numeric_limits<double>::max();
1658 : #ifdef naive_update_not_used
1659 : // This native method uses the fact that stddev = sqrt(sum_of_squares/N -
1660 : // mean^2) and that thus sum_of_squares = N * (stddev^2 + mean^2)
1661 : double dfGlobalSum = 0;
1662 : double dfGlobalSumSquare = 0;
1663 : #else
1664 : // This method uses
1665 : // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Parallel_algorithm
1666 : // which is more numerically robust
1667 : double dfGlobalMean = 0;
1668 : double dfGlobalM2 = 0;
1669 : #endif
1670 : uint64_t nGlobalValidPixels = 0;
1671 : uint64_t nTotalPixelsOfSources = 0;
1672 :
1673 : // Keep original values from single source to avoid slight changes
1674 : // due to recomputation. Cf https://github.com/OSGeo/gdal/issues/12650
1675 : bool bUpdateStatsWithConstantValueRun = false;
1676 : double dfSingleSourceMin = 0;
1677 : double dfSingleSourceMax = 0;
1678 : double dfSingleSourceMean = 0;
1679 : double dfSingleSourceStdDev = 0;
1680 : };
1681 : } // namespace
1682 :
1683 3 : static void UpdateStatsWithConstantValue(Context &sContext, double dfVal,
1684 : uint64_t nPixelCount)
1685 : {
1686 3 : sContext.bUpdateStatsWithConstantValueRun = true;
1687 3 : sContext.dfGlobalMin = std::min(sContext.dfGlobalMin, dfVal);
1688 3 : sContext.dfGlobalMax = std::max(sContext.dfGlobalMax, dfVal);
1689 : #ifdef naive_update_not_used
1690 : sContext.dfGlobalSum += dfVal * nPixelCount;
1691 : sContext.dfGlobalSumSquare += dfVal * dfVal * nPixelCount;
1692 : #else
1693 3 : const auto nNewGlobalValidPixels =
1694 3 : sContext.nGlobalValidPixels + nPixelCount;
1695 3 : const double dfDelta = dfVal - sContext.dfGlobalMean;
1696 3 : sContext.dfGlobalMean += nPixelCount * dfDelta / nNewGlobalValidPixels;
1697 3 : sContext.dfGlobalM2 += dfDelta * dfDelta * nPixelCount *
1698 3 : sContext.nGlobalValidPixels / nNewGlobalValidPixels;
1699 : #endif
1700 3 : sContext.nGlobalValidPixels += nPixelCount;
1701 3 : }
1702 :
1703 : /************************************************************************/
1704 : /* ComputeStatistics() */
1705 : /************************************************************************/
1706 :
1707 36 : CPLErr VRTSourcedRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin,
1708 : double *pdfMax, double *pdfMean,
1709 : double *pdfStdDev,
1710 : GDALProgressFunc pfnProgress,
1711 : void *pProgressData,
1712 : CSLConstList papszOptions)
1713 :
1714 : {
1715 72 : const std::string osFctId("VRTSourcedRasterBand::ComputeStatistics");
1716 72 : GDALAntiRecursionGuard oGuard(osFctId);
1717 36 : if (oGuard.GetCallDepth() >= 32)
1718 : {
1719 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1720 0 : return CE_Failure;
1721 : }
1722 :
1723 108 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1724 36 : if (oGuard2.GetCallDepth() >= 2)
1725 : {
1726 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1727 0 : return CE_Failure;
1728 : }
1729 :
1730 : const bool bSetStatistics =
1731 36 : CPLFetchBool(papszOptions, "SET_STATISTICS", true);
1732 :
1733 : /* -------------------------------------------------------------------- */
1734 : /* If we have overview bands, use them for statistics. */
1735 : /* -------------------------------------------------------------------- */
1736 36 : if (bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews())
1737 : {
1738 : GDALRasterBand *const poBand =
1739 2 : GetRasterSampleOverview(GDALSTAT_APPROX_NUMSAMPLES);
1740 :
1741 2 : if (poBand != nullptr && poBand != this)
1742 : {
1743 2 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
1744 : CPLErr eErr;
1745 3 : if (l_poDS && !l_poDS->m_apoOverviews.empty() &&
1746 1 : dynamic_cast<VRTSourcedRasterBand *>(poBand) != nullptr)
1747 : {
1748 2 : auto apoTmpOverviews = std::move(l_poDS->m_apoOverviews);
1749 1 : l_poDS->m_apoOverviews.clear();
1750 1 : eErr = poBand->GDALRasterBand::ComputeStatistics(
1751 : TRUE, pdfMin, pdfMax, pdfMean, pdfStdDev, pfnProgress,
1752 : pProgressData, papszOptions);
1753 1 : l_poDS->m_apoOverviews = std::move(apoTmpOverviews);
1754 : }
1755 : else
1756 : {
1757 1 : eErr = poBand->ComputeStatistics(TRUE, pdfMin, pdfMax, pdfMean,
1758 : pdfStdDev, pfnProgress,
1759 1 : pProgressData, papszOptions);
1760 : }
1761 2 : if (eErr == CE_None && pdfMin && pdfMax && pdfMean && pdfStdDev &&
1762 : bSetStatistics)
1763 : {
1764 2 : SetMetadataItem("STATISTICS_APPROXIMATE", "YES");
1765 2 : SetMetadataItem(
1766 : "STATISTICS_VALID_PERCENT",
1767 2 : poBand->GetMetadataItem("STATISTICS_VALID_PERCENT"));
1768 2 : SetStatistics(*pdfMin, *pdfMax, *pdfMean, *pdfStdDev);
1769 : }
1770 2 : return eErr;
1771 : }
1772 : }
1773 :
1774 34 : if (IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1775 : /*bAllowMaxValAdjustment = */ false))
1776 : {
1777 28 : Context sContext;
1778 28 : sContext.bApproxOK = CPL_TO_BOOL(bApproxOK);
1779 28 : sContext.dfNoDataValue = m_dfNoDataValue;
1780 28 : sContext.bNoDataValueSet = m_bNoDataValueSet;
1781 28 : sContext.bHideNoDataValue = CPL_TO_BOOL(m_bHideNoDataValue);
1782 28 : sContext.pfnProgress = pfnProgress;
1783 28 : sContext.pProgressData = pProgressData;
1784 28 : sContext.nSources = static_cast<int>(m_papoSources.size());
1785 :
1786 : struct Job
1787 : {
1788 : Context *psContext = nullptr;
1789 : GDALRasterBand *poRasterBand = nullptr;
1790 : uint64_t nPixelCount = 0;
1791 : uint64_t nLastAddedPixels = 0;
1792 : uint64_t nValidPixels = 0;
1793 : double dfMin = 0;
1794 : double dfMax = 0;
1795 : double dfMean = 0;
1796 : double dfStdDev = 0;
1797 : CSLConstList papszOptions = nullptr;
1798 :
1799 58 : static int CPL_STDCALL ProgressFunc(double dfComplete,
1800 : const char *pszMessage,
1801 : void *pProgressArg)
1802 : {
1803 58 : auto psJob = static_cast<Job *>(pProgressArg);
1804 58 : auto psContext = psJob->psContext;
1805 58 : const uint64_t nNewAddedPixels =
1806 : dfComplete == 1.0
1807 58 : ? psJob->nPixelCount
1808 : : static_cast<uint64_t>(
1809 54 : dfComplete * psJob->nPixelCount + 0.5);
1810 : const auto nUpdateThreshold =
1811 116 : std::min(psContext->nTotalPixelsOfSources / 1000,
1812 58 : static_cast<uint64_t>(1000 * 1000));
1813 116 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1814 58 : psContext->nTotalIteratedPixels +=
1815 58 : (nNewAddedPixels - psJob->nLastAddedPixels);
1816 58 : psJob->nLastAddedPixels = nNewAddedPixels;
1817 58 : if (psContext->nTotalIteratedPixels ==
1818 58 : psContext->nTotalPixelsOfSources)
1819 : {
1820 2 : psContext->nLastReportedPixels =
1821 2 : psContext->nTotalIteratedPixels;
1822 2 : return psContext->pfnProgress(1.0, pszMessage,
1823 2 : psContext->pProgressData);
1824 : }
1825 56 : else if (psContext->nTotalIteratedPixels -
1826 56 : psContext->nLastReportedPixels >
1827 : nUpdateThreshold)
1828 : {
1829 48 : psContext->nLastReportedPixels =
1830 48 : psContext->nTotalIteratedPixels;
1831 96 : return psContext->pfnProgress(
1832 48 : static_cast<double>(psContext->nTotalIteratedPixels) /
1833 48 : psContext->nTotalPixelsOfSources,
1834 48 : pszMessage, psContext->pProgressData);
1835 : }
1836 8 : return 1;
1837 : }
1838 :
1839 39 : static void UpdateStats(const Job *psJob)
1840 : {
1841 39 : const auto nValidPixels = psJob->nValidPixels;
1842 39 : auto psContext = psJob->psContext;
1843 39 : if (nValidPixels > 0)
1844 : {
1845 29 : psContext->dfGlobalMin =
1846 29 : std::min(psContext->dfGlobalMin, psJob->dfMin);
1847 29 : psContext->dfGlobalMax =
1848 29 : std::max(psContext->dfGlobalMax, psJob->dfMax);
1849 : #ifdef naive_update_not_used
1850 : psContext->dfGlobalSum += nValidPixels * psJob->dfMean;
1851 : psContext->dfGlobalSumSquare +=
1852 : nValidPixels * (psJob->dfStdDev * psJob->dfStdDev +
1853 : psJob->dfMean * psJob->dfMean);
1854 : psContext->nGlobalValidPixels += nValidPixels;
1855 : #else
1856 29 : const auto nNewGlobalValidPixels =
1857 29 : psContext->nGlobalValidPixels + nValidPixels;
1858 29 : const double dfDelta =
1859 29 : psJob->dfMean - psContext->dfGlobalMean;
1860 29 : psContext->dfGlobalMean +=
1861 29 : nValidPixels * dfDelta / nNewGlobalValidPixels;
1862 29 : psContext->dfGlobalM2 +=
1863 29 : nValidPixels * psJob->dfStdDev * psJob->dfStdDev +
1864 29 : dfDelta * dfDelta * nValidPixels *
1865 29 : psContext->nGlobalValidPixels /
1866 : nNewGlobalValidPixels;
1867 29 : psContext->nGlobalValidPixels = nNewGlobalValidPixels;
1868 : #endif
1869 : }
1870 39 : int bHasNoData = FALSE;
1871 : const double dfNoDataValue =
1872 39 : psJob->poRasterBand->GetNoDataValue(&bHasNoData);
1873 11 : if (nValidPixels < psJob->nPixelCount && bHasNoData &&
1874 61 : !std::isnan(dfNoDataValue) &&
1875 11 : (!psContext->bNoDataValueSet ||
1876 9 : dfNoDataValue != psContext->dfNoDataValue))
1877 : {
1878 : const auto eBandDT =
1879 2 : psJob->poRasterBand->GetRasterDataType();
1880 : // Check that the band nodata value is in the range of the
1881 : // original raster type
1882 : GByte abyTempBuffer[2 * sizeof(double)];
1883 2 : CPLAssert(GDALGetDataTypeSizeBytes(eBandDT) <=
1884 : static_cast<int>(sizeof(abyTempBuffer)));
1885 2 : GDALCopyWords(&dfNoDataValue, GDT_Float64, 0,
1886 : &abyTempBuffer[0], eBandDT, 0, 1);
1887 2 : double dfNoDataValueAfter = dfNoDataValue;
1888 2 : GDALCopyWords(&abyTempBuffer[0], eBandDT, 0,
1889 : &dfNoDataValueAfter, GDT_Float64, 0, 1);
1890 4 : if (!std::isfinite(dfNoDataValue) ||
1891 2 : std::fabs(dfNoDataValueAfter - dfNoDataValue) < 1.0)
1892 : {
1893 2 : UpdateStatsWithConstantValue(
1894 : *psContext, dfNoDataValueAfter,
1895 2 : psJob->nPixelCount - nValidPixels);
1896 : }
1897 : }
1898 :
1899 39 : if (psContext->nSources == 1)
1900 : {
1901 13 : psContext->dfSingleSourceMin = psJob->dfMin;
1902 13 : psContext->dfSingleSourceMax = psJob->dfMax;
1903 13 : psContext->dfSingleSourceMean = psJob->dfMean;
1904 13 : psContext->dfSingleSourceStdDev = psJob->dfStdDev;
1905 : }
1906 39 : }
1907 : };
1908 :
1909 42 : const auto JobRunner = [](void *pData)
1910 : {
1911 42 : auto psJob = static_cast<Job *>(pData);
1912 42 : auto psContext = psJob->psContext;
1913 : {
1914 42 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1915 42 : if (psContext->bFallbackToBase || psContext->bFailure)
1916 1 : return;
1917 : }
1918 :
1919 41 : auto poSimpleSourceBand = psJob->poRasterBand;
1920 41 : psJob->nPixelCount =
1921 41 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
1922 41 : poSimpleSourceBand->GetYSize();
1923 :
1924 41 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1925 41 : CPLErr eErr = poSimpleSourceBand->ComputeStatistics(
1926 41 : psContext->bApproxOK, &psJob->dfMin, &psJob->dfMax,
1927 : &psJob->dfMean, &psJob->dfStdDev,
1928 47 : psContext->pfnProgress == nullptr ||
1929 6 : psContext->pfnProgress == GDALDummyProgress
1930 : ? GDALDummyProgress
1931 : : Job::ProgressFunc,
1932 41 : psJob, psJob->papszOptions);
1933 : const char *pszValidPercent =
1934 41 : poSimpleSourceBand->GetMetadataItem("STATISTICS_VALID_PERCENT");
1935 41 : psJob->nValidPixels =
1936 : pszValidPercent
1937 41 : ? static_cast<uint64_t>(CPLAtof(pszValidPercent) *
1938 41 : psJob->nPixelCount / 100.0)
1939 : : psJob->nPixelCount;
1940 41 : if (eErr == CE_Failure)
1941 : {
1942 20 : if (pszValidPercent != nullptr &&
1943 10 : CPLAtof(pszValidPercent) == 0.0)
1944 : {
1945 : // ok: no valid sample
1946 : }
1947 : else
1948 : {
1949 0 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1950 0 : psContext->bFailure = true;
1951 : }
1952 : }
1953 : else
1954 : {
1955 31 : int bHasNoData = FALSE;
1956 31 : CPL_IGNORE_RET_VAL(
1957 31 : psJob->poRasterBand->GetNoDataValue(&bHasNoData));
1958 31 : if (!bHasNoData && psContext->bNoDataValueSet &&
1959 8 : !psContext->bHideNoDataValue &&
1960 6 : psContext->dfNoDataValue >= psJob->dfMin &&
1961 2 : psContext->dfNoDataValue <= psJob->dfMax)
1962 : {
1963 2 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1964 2 : psJob->psContext->bFallbackToBase = true;
1965 2 : return;
1966 : }
1967 : }
1968 : };
1969 :
1970 28 : CPLWorkerThreadPool *poThreadPool = nullptr;
1971 : int nThreads =
1972 28 : m_papoSources.size() > 1
1973 28 : ? VRTDataset::GetNumThreads(dynamic_cast<VRTDataset *>(poDS))
1974 28 : : 0;
1975 28 : if (nThreads > 1024)
1976 0 : nThreads = 1024; // to please Coverity
1977 28 : if (nThreads > 1)
1978 : {
1979 : // Check that all sources refer to different datasets
1980 : // before allowing multithreaded access
1981 : // If the datasets belong to the MEM driver, check GDALDataset*
1982 : // pointer values. Otherwise use dataset name.
1983 26 : std::set<std::string> oSetDatasetNames;
1984 26 : std::set<GDALDataset *> oSetDatasetPointers;
1985 39 : for (auto &poSource : m_papoSources)
1986 : {
1987 : auto poSimpleSource =
1988 26 : static_cast<VRTSimpleSource *>(poSource.get());
1989 26 : assert(poSimpleSource);
1990 26 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1991 26 : assert(poSimpleSourceBand);
1992 26 : auto poSourceDataset = poSimpleSourceBand->GetDataset();
1993 26 : if (poSourceDataset == nullptr)
1994 : {
1995 0 : nThreads = 0;
1996 0 : break;
1997 : }
1998 26 : auto poDriver = poSourceDataset->GetDriver();
1999 26 : if (poDriver && EQUAL(poDriver->GetDescription(), "MEM"))
2000 : {
2001 26 : if (oSetDatasetPointers.find(poSourceDataset) !=
2002 52 : oSetDatasetPointers.end())
2003 : {
2004 0 : nThreads = 0;
2005 0 : break;
2006 : }
2007 26 : oSetDatasetPointers.insert(poSourceDataset);
2008 : }
2009 : else
2010 : {
2011 0 : if (oSetDatasetNames.find(
2012 0 : poSourceDataset->GetDescription()) !=
2013 0 : oSetDatasetNames.end())
2014 : {
2015 0 : nThreads = 0;
2016 0 : break;
2017 : }
2018 0 : oSetDatasetNames.insert(poSourceDataset->GetDescription());
2019 : }
2020 : }
2021 13 : if (nThreads > 1)
2022 : {
2023 13 : poThreadPool = GDALGetGlobalThreadPool(nThreads);
2024 : }
2025 : }
2026 :
2027 : // Compute total number of pixels of sources
2028 70 : for (auto &poSource : m_papoSources)
2029 : {
2030 : auto poSimpleSource =
2031 42 : static_cast<VRTSimpleSource *>(poSource.get());
2032 42 : assert(poSimpleSource);
2033 42 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
2034 42 : assert(poSimpleSourceBand);
2035 42 : sContext.nTotalPixelsOfSources +=
2036 42 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
2037 42 : poSimpleSourceBand->GetYSize();
2038 : }
2039 :
2040 28 : if (poThreadPool)
2041 : {
2042 13 : CPLDebugOnly("VRT", "ComputeStatistics(): use optimized "
2043 : "multi-threaded code path for mosaic");
2044 26 : std::vector<Job> asJobs(m_papoSources.size());
2045 26 : auto poQueue = poThreadPool->CreateJobQueue();
2046 39 : for (size_t i = 0; i < m_papoSources.size(); ++i)
2047 : {
2048 : auto poSimpleSource =
2049 26 : static_cast<VRTSimpleSource *>(m_papoSources[i].get());
2050 26 : assert(poSimpleSource);
2051 26 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
2052 26 : assert(poSimpleSourceBand);
2053 26 : asJobs[i].psContext = &sContext;
2054 26 : asJobs[i].poRasterBand = poSimpleSourceBand;
2055 26 : asJobs[i].papszOptions = papszOptions;
2056 26 : if (!poQueue->SubmitJob(JobRunner, &asJobs[i]))
2057 : {
2058 0 : sContext.bFailure = true;
2059 0 : break;
2060 : }
2061 : }
2062 13 : poQueue->WaitCompletion();
2063 13 : if (!(sContext.bFailure || sContext.bFallbackToBase))
2064 : {
2065 36 : for (size_t i = 0; i < m_papoSources.size(); ++i)
2066 : {
2067 24 : Job::UpdateStats(&asJobs[i]);
2068 : }
2069 : }
2070 : }
2071 : else
2072 : {
2073 15 : CPLDebugOnly(
2074 : "VRT",
2075 : "ComputeStatistics(): use optimized code path for mosaic");
2076 30 : for (auto &poSource : m_papoSources)
2077 : {
2078 : auto poSimpleSource =
2079 16 : static_cast<VRTSimpleSource *>(poSource.get());
2080 16 : assert(poSimpleSource);
2081 16 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
2082 16 : assert(poSimpleSourceBand);
2083 16 : Job sJob;
2084 16 : sJob.psContext = &sContext;
2085 16 : sJob.poRasterBand = poSimpleSourceBand;
2086 16 : sJob.papszOptions = papszOptions;
2087 16 : JobRunner(&sJob);
2088 16 : if (sContext.bFailure || sContext.bFallbackToBase)
2089 : break;
2090 15 : Job::UpdateStats(&sJob);
2091 : }
2092 : }
2093 :
2094 28 : if (sContext.bFailure)
2095 0 : return CE_Failure;
2096 28 : if (sContext.bFallbackToBase)
2097 : {
2098 : // If the VRT band nodata value is in the [min, max] range
2099 : // of the source and that the source has no nodata value set,
2100 : // then we can't use the optimization.
2101 2 : CPLDebugOnly("VRT", "ComputeStatistics(): revert back to "
2102 : "generic case because of nodata value in range "
2103 : "of source raster");
2104 2 : return GDALRasterBand::ComputeStatistics(
2105 : bApproxOK, pdfMin, pdfMax, pdfMean, pdfStdDev, pfnProgress,
2106 2 : pProgressData, papszOptions);
2107 : }
2108 :
2109 26 : const uint64_t nTotalPixels =
2110 26 : static_cast<uint64_t>(nRasterXSize) * nRasterYSize;
2111 8 : if (m_bNoDataValueSet && m_bHideNoDataValue &&
2112 34 : !std::isnan(m_dfNoDataValue) && IsNoDataValueInDataTypeRange())
2113 : {
2114 1 : UpdateStatsWithConstantValue(sContext, m_dfNoDataValue,
2115 : nTotalPixels -
2116 1 : sContext.nGlobalValidPixels);
2117 : }
2118 25 : else if (!m_bNoDataValueSet)
2119 : {
2120 18 : sContext.nGlobalValidPixels = nTotalPixels;
2121 : }
2122 :
2123 : #ifdef naive_update_not_used
2124 : double dfGlobalMean =
2125 : sContext.nGlobalValidPixels > 0
2126 : ? sContext.dfGlobalSum / sContext.nGlobalValidPixels
2127 : : 0;
2128 : double dfGlobalStdDev = sContext.nGlobalValidPixels > 0
2129 : ? sqrt(sContext.dfGlobalSumSquare /
2130 : sContext.nGlobalValidPixels -
2131 : dfGlobalMean * dfGlobalMean)
2132 : : 0;
2133 : #else
2134 26 : double dfGlobalMean = sContext.dfGlobalMean;
2135 : double dfGlobalStdDev =
2136 26 : sContext.nGlobalValidPixels > 0
2137 26 : ? sqrt(sContext.dfGlobalM2 / sContext.nGlobalValidPixels)
2138 26 : : 0;
2139 : #endif
2140 39 : if (m_papoSources.size() == 1 &&
2141 13 : !sContext.bUpdateStatsWithConstantValueRun)
2142 : {
2143 13 : sContext.dfGlobalMin = sContext.dfSingleSourceMin;
2144 13 : sContext.dfGlobalMax = sContext.dfSingleSourceMax;
2145 13 : dfGlobalMean = sContext.dfSingleSourceMean;
2146 13 : dfGlobalStdDev = sContext.dfSingleSourceStdDev;
2147 : }
2148 :
2149 26 : if (sContext.nGlobalValidPixels > 0)
2150 : {
2151 23 : if (bSetStatistics)
2152 : {
2153 23 : if (bApproxOK)
2154 : {
2155 1 : SetMetadataItem("STATISTICS_APPROXIMATE", "YES");
2156 : }
2157 22 : else if (GetMetadataItem("STATISTICS_APPROXIMATE"))
2158 : {
2159 0 : SetMetadataItem("STATISTICS_APPROXIMATE", nullptr);
2160 : }
2161 23 : SetStatistics(sContext.dfGlobalMin, sContext.dfGlobalMax,
2162 23 : dfGlobalMean, dfGlobalStdDev);
2163 : }
2164 : }
2165 : else
2166 : {
2167 3 : sContext.dfGlobalMin = 0.0;
2168 3 : sContext.dfGlobalMax = 0.0;
2169 : }
2170 :
2171 26 : if (bSetStatistics)
2172 26 : SetValidPercent(nTotalPixels, sContext.nGlobalValidPixels);
2173 :
2174 26 : if (pdfMin)
2175 18 : *pdfMin = sContext.dfGlobalMin;
2176 26 : if (pdfMax)
2177 18 : *pdfMax = sContext.dfGlobalMax;
2178 26 : if (pdfMean)
2179 13 : *pdfMean = dfGlobalMean;
2180 26 : if (pdfStdDev)
2181 13 : *pdfStdDev = dfGlobalStdDev;
2182 :
2183 26 : if (sContext.nGlobalValidPixels == 0)
2184 : {
2185 3 : ReportError(CE_Failure, CPLE_AppDefined,
2186 : "Failed to compute statistics, no valid pixels found "
2187 : "in sampling.");
2188 : }
2189 :
2190 26 : return sContext.nGlobalValidPixels > 0 ? CE_None : CE_Failure;
2191 : }
2192 : else
2193 : {
2194 6 : return GDALRasterBand::ComputeStatistics(
2195 : bApproxOK, pdfMin, pdfMax, pdfMean, pdfStdDev, pfnProgress,
2196 6 : pProgressData, papszOptions);
2197 : }
2198 : }
2199 :
2200 : /************************************************************************/
2201 : /* GetHistogram() */
2202 : /************************************************************************/
2203 :
2204 7 : CPLErr VRTSourcedRasterBand::GetHistogram(double dfMin, double dfMax,
2205 : int nBuckets, GUIntBig *panHistogram,
2206 : int bIncludeOutOfRange, int bApproxOK,
2207 : GDALProgressFunc pfnProgress,
2208 : void *pProgressData)
2209 :
2210 : {
2211 : /* -------------------------------------------------------------------- */
2212 : /* If we have overviews, use them for the histogram. */
2213 : /* -------------------------------------------------------------------- */
2214 7 : if (bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews())
2215 : {
2216 : // FIXME: Should we use the most reduced overview here or use some
2217 : // minimum number of samples like GDALRasterBand::ComputeStatistics()
2218 : // does?
2219 1 : GDALRasterBand *poBand = GetRasterSampleOverview(0);
2220 :
2221 1 : if (poBand != nullptr && poBand != this)
2222 : {
2223 1 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
2224 2 : if (l_poDS && !l_poDS->m_apoOverviews.empty() &&
2225 1 : dynamic_cast<VRTSourcedRasterBand *>(poBand) != nullptr)
2226 : {
2227 1 : auto apoTmpOverviews = std::move(l_poDS->m_apoOverviews);
2228 1 : l_poDS->m_apoOverviews.clear();
2229 1 : auto eErr = poBand->GDALRasterBand::GetHistogram(
2230 : dfMin, dfMax, nBuckets, panHistogram, bIncludeOutOfRange,
2231 : bApproxOK, pfnProgress, pProgressData);
2232 1 : l_poDS->m_apoOverviews = std::move(apoTmpOverviews);
2233 1 : return eErr;
2234 : }
2235 : else
2236 : {
2237 0 : return poBand->GetHistogram(
2238 : dfMin, dfMax, nBuckets, panHistogram, bIncludeOutOfRange,
2239 0 : bApproxOK, pfnProgress, pProgressData);
2240 : }
2241 : }
2242 : }
2243 :
2244 6 : if (m_papoSources.size() != 1)
2245 2 : return VRTRasterBand::GetHistogram(dfMin, dfMax, nBuckets, panHistogram,
2246 : bIncludeOutOfRange, bApproxOK,
2247 2 : pfnProgress, pProgressData);
2248 :
2249 4 : if (pfnProgress == nullptr)
2250 4 : pfnProgress = GDALDummyProgress;
2251 :
2252 8 : const std::string osFctId("VRTSourcedRasterBand::GetHistogram");
2253 8 : GDALAntiRecursionGuard oGuard(osFctId);
2254 4 : if (oGuard.GetCallDepth() >= 32)
2255 : {
2256 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
2257 0 : return CE_Failure;
2258 : }
2259 :
2260 12 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
2261 4 : if (oGuard2.GetCallDepth() >= 2)
2262 : {
2263 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
2264 0 : return CE_Failure;
2265 : }
2266 :
2267 : /* -------------------------------------------------------------------- */
2268 : /* Try with source bands. */
2269 : /* -------------------------------------------------------------------- */
2270 8 : const CPLErr eErr = m_papoSources[0]->GetHistogram(
2271 : GetXSize(), GetYSize(), dfMin, dfMax, nBuckets, panHistogram,
2272 4 : bIncludeOutOfRange, bApproxOK, pfnProgress, pProgressData);
2273 4 : if (eErr != CE_None)
2274 : {
2275 0 : const CPLErr eErr2 = GDALRasterBand::GetHistogram(
2276 : dfMin, dfMax, nBuckets, panHistogram, bIncludeOutOfRange, bApproxOK,
2277 : pfnProgress, pProgressData);
2278 0 : return eErr2;
2279 : }
2280 :
2281 4 : SetDefaultHistogram(dfMin, dfMax, nBuckets, panHistogram);
2282 :
2283 4 : return CE_None;
2284 : }
2285 :
2286 : /************************************************************************/
2287 : /* AddSource() */
2288 : /************************************************************************/
2289 :
2290 241950 : CPLErr VRTSourcedRasterBand::AddSource(VRTSource *poNewSource)
2291 :
2292 : {
2293 241950 : return AddSource(std::unique_ptr<VRTSource>(poNewSource));
2294 : }
2295 :
2296 : /************************************************************************/
2297 : /* AddSource() */
2298 : /************************************************************************/
2299 :
2300 247110 : CPLErr VRTSourcedRasterBand::AddSource(std::unique_ptr<VRTSource> poNewSource)
2301 :
2302 : {
2303 247110 : auto l_poDS = static_cast<VRTDataset *>(poDS);
2304 247110 : l_poDS->SetNeedsFlush();
2305 247110 : l_poDS->SourceAdded();
2306 :
2307 247110 : if (poNewSource->IsSimpleSource())
2308 : {
2309 : VRTSimpleSource *poSS =
2310 247075 : static_cast<VRTSimpleSource *>(poNewSource.get());
2311 247075 : if (GetMetadataItem(GDALMD_NBITS, GDAL_MDD_IMAGE_STRUCTURE) != nullptr)
2312 : {
2313 : int nBits =
2314 47 : atoi(GetMetadataItem(GDALMD_NBITS, GDAL_MDD_IMAGE_STRUCTURE));
2315 47 : if (nBits >= 1 && nBits <= 31)
2316 : {
2317 47 : poSS->SetMaxValue(static_cast<int>((1U << nBits) - 1));
2318 : }
2319 : }
2320 : }
2321 :
2322 247110 : m_papoSources.push_back(std::move(poNewSource));
2323 :
2324 247110 : return CE_None;
2325 : }
2326 :
2327 : /*! @endcond */
2328 :
2329 : /************************************************************************/
2330 : /* VRTAddSource() */
2331 : /************************************************************************/
2332 :
2333 : /**
2334 : * @see VRTSourcedRasterBand::AddSource().
2335 : */
2336 :
2337 0 : CPLErr CPL_STDCALL VRTAddSource(VRTSourcedRasterBandH hVRTBand,
2338 : VRTSourceH hNewSource)
2339 : {
2340 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddSource", CE_Failure);
2341 :
2342 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddSource(
2343 0 : reinterpret_cast<VRTSource *>(hNewSource));
2344 : }
2345 :
2346 : /*! @cond Doxygen_Suppress */
2347 :
2348 : /************************************************************************/
2349 : /* XMLInit() */
2350 : /************************************************************************/
2351 :
2352 2546 : CPLErr VRTSourcedRasterBand::XMLInit(const CPLXMLNode *psTree,
2353 : const char *pszVRTPath,
2354 : VRTMapSharedResources &oMapSharedSources)
2355 :
2356 : {
2357 : {
2358 : const CPLErr eErr =
2359 2546 : VRTRasterBand::XMLInit(psTree, pszVRTPath, oMapSharedSources);
2360 2546 : if (eErr != CE_None)
2361 3 : return eErr;
2362 : }
2363 :
2364 : /* -------------------------------------------------------------------- */
2365 : /* Process sources. */
2366 : /* -------------------------------------------------------------------- */
2367 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
2368 2543 : GetGDALDriverManager()->GetDriverByName("VRT"));
2369 :
2370 2543 : for (const CPLXMLNode *psChild = psTree->psChild;
2371 118410 : psChild != nullptr && poDriver != nullptr; psChild = psChild->psNext)
2372 : {
2373 115889 : if (psChild->eType != CXT_Element)
2374 6667 : continue;
2375 :
2376 109222 : CPLErrorReset();
2377 : VRTSource *const poSource =
2378 109222 : poDriver->ParseSource(psChild, pszVRTPath, oMapSharedSources);
2379 109222 : if (poSource != nullptr)
2380 104160 : AddSource(poSource);
2381 5062 : else if (CPLGetLastErrorType() != CE_None)
2382 22 : return CE_Failure;
2383 : }
2384 :
2385 : /* -------------------------------------------------------------------- */
2386 : /* Done. */
2387 : /* -------------------------------------------------------------------- */
2388 : const char *pszSubclass =
2389 2521 : CPLGetXMLValue(psTree, "subclass", "VRTSourcedRasterBand");
2390 2521 : if (m_papoSources.empty() && !EQUAL(pszSubclass, "VRTDerivedRasterBand"))
2391 132 : CPLDebug("VRT", "No valid sources found for band in VRT file %s",
2392 132 : GetDataset() ? GetDataset()->GetDescription() : "");
2393 :
2394 2521 : return CE_None;
2395 : }
2396 :
2397 : /************************************************************************/
2398 : /* SerializeToXML() */
2399 : /************************************************************************/
2400 :
2401 657 : CPLXMLNode *VRTSourcedRasterBand::SerializeToXML(const char *pszVRTPath,
2402 : bool &bHasWarnedAboutRAMUsage,
2403 : size_t &nAccRAMUsage)
2404 :
2405 : {
2406 657 : CPLXMLNode *psTree = VRTRasterBand::SerializeToXML(
2407 : pszVRTPath, bHasWarnedAboutRAMUsage, nAccRAMUsage);
2408 657 : CPLXMLNode *psLastChild = psTree->psChild;
2409 2393 : while (psLastChild != nullptr && psLastChild->psNext != nullptr)
2410 1736 : psLastChild = psLastChild->psNext;
2411 :
2412 : /* -------------------------------------------------------------------- */
2413 : /* Process Sources. */
2414 : /* -------------------------------------------------------------------- */
2415 :
2416 657 : GIntBig nUsableRAM = -1;
2417 :
2418 3274 : for (const auto &poSource : m_papoSources)
2419 : {
2420 2617 : CPLXMLNode *const psXMLSrc = poSource->SerializeToXML(pszVRTPath);
2421 :
2422 2617 : if (psXMLSrc == nullptr)
2423 0 : break;
2424 :
2425 : // Creating the CPLXMLNode tree representation of a VRT can easily
2426 : // take several times RAM usage than its string serialization, or its
2427 : // internal representation in the driver.
2428 : // We multiply the estimate by a factor of 2, experimentally found to
2429 : // be more realistic than the conservative raw estimate.
2430 2617 : nAccRAMUsage += 2 * CPLXMLNodeGetRAMUsageEstimate(psXMLSrc);
2431 2617 : if (!bHasWarnedAboutRAMUsage && nAccRAMUsage > 512 * 1024 * 1024)
2432 : {
2433 0 : if (nUsableRAM < 0)
2434 0 : nUsableRAM = CPLGetUsablePhysicalRAM();
2435 0 : if (nUsableRAM > 0 &&
2436 0 : nAccRAMUsage > static_cast<uint64_t>(nUsableRAM) / 10 * 8)
2437 : {
2438 0 : bHasWarnedAboutRAMUsage = true;
2439 0 : CPLError(CE_Warning, CPLE_AppDefined,
2440 : "Serialization of this VRT file has already consumed "
2441 : "at least %.02f GB of RAM over a total of %.02f. This "
2442 : "process may abort",
2443 0 : double(nAccRAMUsage) / (1024 * 1024 * 1024),
2444 0 : double(nUsableRAM) / (1024 * 1024 * 1024));
2445 : }
2446 : }
2447 :
2448 2617 : if (psLastChild == nullptr)
2449 0 : psTree->psChild = psXMLSrc;
2450 : else
2451 2617 : psLastChild->psNext = psXMLSrc;
2452 2617 : psLastChild = psXMLSrc;
2453 : }
2454 :
2455 657 : return psTree;
2456 : }
2457 :
2458 : /************************************************************************/
2459 : /* SkipBufferInitialization() */
2460 : /************************************************************************/
2461 :
2462 289874 : bool VRTSourcedRasterBand::SkipBufferInitialization()
2463 : {
2464 289874 : if (m_nSkipBufferInitialization >= 0)
2465 284598 : return m_nSkipBufferInitialization != 0;
2466 : /* -------------------------------------------------------------------- */
2467 : /* Check if we can avoid buffer initialization. */
2468 : /* -------------------------------------------------------------------- */
2469 :
2470 : // Note: if one day we do alpha compositing, we will need to check that.
2471 5276 : m_nSkipBufferInitialization = FALSE;
2472 5276 : if (m_papoSources.size() != 1 || !m_papoSources[0]->IsSimpleSource())
2473 : {
2474 3558 : return false;
2475 : }
2476 : VRTSimpleSource *poSS =
2477 1718 : static_cast<VRTSimpleSource *>(m_papoSources[0].get());
2478 1718 : if (poSS->GetType() == VRTSimpleSource::GetTypeStatic())
2479 : {
2480 1300 : auto l_poBand = poSS->GetRasterBand();
2481 1289 : if (l_poBand != nullptr && poSS->m_dfSrcXOff >= 0.0 &&
2482 1251 : poSS->m_dfSrcYOff >= 0.0 &&
2483 1249 : poSS->m_dfSrcXOff + poSS->m_dfSrcXSize <= l_poBand->GetXSize() &&
2484 1247 : poSS->m_dfSrcYOff + poSS->m_dfSrcYSize <= l_poBand->GetYSize() &&
2485 1247 : poSS->m_dfDstXOff <= 0.0 && poSS->m_dfDstYOff <= 0.0 &&
2486 3757 : poSS->m_dfDstXOff + poSS->m_dfDstXSize >= nRasterXSize &&
2487 1168 : poSS->m_dfDstYOff + poSS->m_dfDstYSize >= nRasterYSize)
2488 : {
2489 1166 : m_nSkipBufferInitialization = TRUE;
2490 : }
2491 : }
2492 1718 : return m_nSkipBufferInitialization != 0;
2493 : }
2494 :
2495 : /************************************************************************/
2496 : /* ConfigureSource() */
2497 : /************************************************************************/
2498 :
2499 142338 : void VRTSourcedRasterBand::ConfigureSource(VRTSimpleSource *poSimpleSource,
2500 : GDALRasterBand *poSrcBand,
2501 : int bAddAsMaskBand, double dfSrcXOff,
2502 : double dfSrcYOff, double dfSrcXSize,
2503 : double dfSrcYSize, double dfDstXOff,
2504 : double dfDstYOff, double dfDstXSize,
2505 : double dfDstYSize)
2506 : {
2507 : /* -------------------------------------------------------------------- */
2508 : /* Default source and dest rectangles. */
2509 : /* -------------------------------------------------------------------- */
2510 142338 : if (dfSrcYSize == -1)
2511 : {
2512 66381 : dfSrcXOff = 0;
2513 66381 : dfSrcYOff = 0;
2514 66381 : dfSrcXSize = poSrcBand->GetXSize();
2515 66381 : dfSrcYSize = poSrcBand->GetYSize();
2516 : }
2517 :
2518 142338 : if (dfDstYSize == -1)
2519 : {
2520 66381 : dfDstXOff = 0;
2521 66381 : dfDstYOff = 0;
2522 66381 : dfDstXSize = nRasterXSize;
2523 66381 : dfDstYSize = nRasterYSize;
2524 : }
2525 :
2526 142338 : if (bAddAsMaskBand)
2527 53 : poSimpleSource->SetSrcMaskBand(poSrcBand);
2528 : else
2529 142285 : poSimpleSource->SetSrcBand(poSrcBand);
2530 :
2531 142338 : poSimpleSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2532 142338 : poSimpleSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2533 :
2534 : /* -------------------------------------------------------------------- */
2535 : /* If we can get the associated GDALDataset, add a reference to it.*/
2536 : /* -------------------------------------------------------------------- */
2537 142338 : GDALDataset *poSrcBandDataset = poSrcBand->GetDataset();
2538 142338 : if (poSrcBandDataset != nullptr)
2539 : {
2540 : VRTDataset *poVRTSrcBandDataset =
2541 142338 : dynamic_cast<VRTDataset *>(poSrcBandDataset);
2542 142338 : if (poVRTSrcBandDataset && !poVRTSrcBandDataset->m_bCanTakeRef)
2543 : {
2544 : // Situation triggered by VRTDataset::AddVirtualOverview()
2545 : // We create an overview dataset that is a VRT of a reduction of
2546 : // ourselves. But we don't want to take a reference on ourselves,
2547 : // otherwise this will prevent us to be closed in number of
2548 : // circumstances
2549 65 : poSimpleSource->m_bDropRefOnSrcBand = false;
2550 : }
2551 : else
2552 : {
2553 142273 : poSrcBandDataset->Reference();
2554 : }
2555 : }
2556 142338 : }
2557 :
2558 : /************************************************************************/
2559 : /* AddSimpleSource() */
2560 : /************************************************************************/
2561 :
2562 377 : CPLErr VRTSourcedRasterBand::AddSimpleSource(
2563 : const char *pszFilename, int nBandIn, double dfSrcXOff, double dfSrcYOff,
2564 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2565 : double dfDstXSize, double dfDstYSize, const char *pszResampling,
2566 : double dfNoDataValueIn)
2567 :
2568 : {
2569 : /* -------------------------------------------------------------------- */
2570 : /* Create source. */
2571 : /* -------------------------------------------------------------------- */
2572 377 : VRTSimpleSource *poSimpleSource = nullptr;
2573 :
2574 377 : if (pszResampling != nullptr && STARTS_WITH_CI(pszResampling, "aver"))
2575 : {
2576 0 : auto poAveragedSource = new VRTAveragedSource();
2577 0 : poSimpleSource = poAveragedSource;
2578 0 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2579 0 : poAveragedSource->SetNoDataValue(dfNoDataValueIn);
2580 : }
2581 : else
2582 : {
2583 377 : poSimpleSource = new VRTSimpleSource();
2584 377 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2585 0 : CPLError(
2586 : CE_Warning, CPLE_AppDefined,
2587 : "NODATA setting not currently supported for nearest "
2588 : "neighbour sampled simple sources on Virtual Datasources.");
2589 : }
2590 :
2591 377 : poSimpleSource->SetSrcBand(pszFilename, nBandIn);
2592 :
2593 377 : poSimpleSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2594 377 : poSimpleSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2595 :
2596 : /* -------------------------------------------------------------------- */
2597 : /* add to list. */
2598 : /* -------------------------------------------------------------------- */
2599 377 : return AddSource(poSimpleSource);
2600 : }
2601 :
2602 : /************************************************************************/
2603 : /* AddSimpleSource() */
2604 : /************************************************************************/
2605 :
2606 67846 : CPLErr VRTSourcedRasterBand::AddSimpleSource(
2607 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2608 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2609 : double dfDstXSize, double dfDstYSize, const char *pszResampling,
2610 : double dfNoDataValueIn)
2611 :
2612 : {
2613 : /* -------------------------------------------------------------------- */
2614 : /* Create source. */
2615 : /* -------------------------------------------------------------------- */
2616 67846 : VRTSimpleSource *poSimpleSource = nullptr;
2617 :
2618 67846 : if (pszResampling != nullptr && STARTS_WITH_CI(pszResampling, "aver"))
2619 : {
2620 0 : auto poAveragedSource = new VRTAveragedSource();
2621 0 : poSimpleSource = poAveragedSource;
2622 0 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2623 0 : poAveragedSource->SetNoDataValue(dfNoDataValueIn);
2624 : }
2625 : else
2626 : {
2627 67846 : poSimpleSource = new VRTSimpleSource();
2628 67846 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2629 0 : CPLError(
2630 : CE_Warning, CPLE_AppDefined,
2631 : "NODATA setting not currently supported for "
2632 : "neighbour sampled simple sources on Virtual Datasources.");
2633 : }
2634 :
2635 67846 : ConfigureSource(poSimpleSource, poSrcBand, FALSE, dfSrcXOff, dfSrcYOff,
2636 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2637 : dfDstYSize);
2638 :
2639 : /* -------------------------------------------------------------------- */
2640 : /* add to list. */
2641 : /* -------------------------------------------------------------------- */
2642 67846 : return AddSource(poSimpleSource);
2643 : }
2644 :
2645 : /************************************************************************/
2646 : /* AddMaskBandSource() */
2647 : /************************************************************************/
2648 :
2649 : // poSrcBand is not the mask band, but the band from which the mask band is
2650 : // taken.
2651 37 : CPLErr VRTSourcedRasterBand::AddMaskBandSource(
2652 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2653 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2654 : double dfDstXSize, double dfDstYSize)
2655 : {
2656 : /* -------------------------------------------------------------------- */
2657 : /* Create source. */
2658 : /* -------------------------------------------------------------------- */
2659 37 : VRTSimpleSource *poSimpleSource = new VRTSimpleSource();
2660 :
2661 37 : ConfigureSource(poSimpleSource, poSrcBand, TRUE, dfSrcXOff, dfSrcYOff,
2662 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2663 : dfDstYSize);
2664 :
2665 : /* -------------------------------------------------------------------- */
2666 : /* add to list. */
2667 : /* -------------------------------------------------------------------- */
2668 37 : return AddSource(poSimpleSource);
2669 : }
2670 :
2671 : /*! @endcond */
2672 :
2673 : /************************************************************************/
2674 : /* VRTAddSimpleSource() */
2675 : /************************************************************************/
2676 :
2677 : /**
2678 : * @see VRTSourcedRasterBand::AddSimpleSource().
2679 : */
2680 :
2681 968 : CPLErr CPL_STDCALL VRTAddSimpleSource(VRTSourcedRasterBandH hVRTBand,
2682 : GDALRasterBandH hSrcBand, int nSrcXOff,
2683 : int nSrcYOff, int nSrcXSize,
2684 : int nSrcYSize, int nDstXOff, int nDstYOff,
2685 : int nDstXSize, int nDstYSize,
2686 : const char *pszResampling,
2687 : double dfNoDataValue)
2688 : {
2689 968 : VALIDATE_POINTER1(hVRTBand, "VRTAddSimpleSource", CE_Failure);
2690 :
2691 968 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddSimpleSource(
2692 : reinterpret_cast<GDALRasterBand *>(hSrcBand), nSrcXOff, nSrcYOff,
2693 : nSrcXSize, nSrcYSize, nDstXOff, nDstYOff, nDstXSize, nDstYSize,
2694 968 : pszResampling, dfNoDataValue);
2695 : }
2696 :
2697 : /*! @cond Doxygen_Suppress */
2698 :
2699 : /************************************************************************/
2700 : /* AddComplexSource() */
2701 : /************************************************************************/
2702 :
2703 29 : CPLErr VRTSourcedRasterBand::AddComplexSource(
2704 : const char *pszFilename, int nBandIn, double dfSrcXOff, double dfSrcYOff,
2705 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2706 : double dfDstXSize, double dfDstYSize, double dfScaleOff,
2707 : double dfScaleRatio, double dfNoDataValueIn, int nColorTableComponent)
2708 :
2709 : {
2710 : /* -------------------------------------------------------------------- */
2711 : /* Create source. */
2712 : /* -------------------------------------------------------------------- */
2713 29 : VRTComplexSource *const poSource = new VRTComplexSource();
2714 :
2715 29 : poSource->SetSrcBand(pszFilename, nBandIn);
2716 :
2717 29 : poSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2718 29 : poSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2719 :
2720 : /* -------------------------------------------------------------------- */
2721 : /* Set complex parameters. */
2722 : /* -------------------------------------------------------------------- */
2723 29 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2724 8 : poSource->SetNoDataValue(dfNoDataValueIn);
2725 :
2726 29 : if (dfScaleOff != 0.0 || dfScaleRatio != 1.0)
2727 4 : poSource->SetLinearScaling(dfScaleOff, dfScaleRatio);
2728 :
2729 29 : poSource->SetColorTableComponent(nColorTableComponent);
2730 :
2731 : /* -------------------------------------------------------------------- */
2732 : /* add to list. */
2733 : /* -------------------------------------------------------------------- */
2734 29 : return AddSource(poSource);
2735 : }
2736 :
2737 : /************************************************************************/
2738 : /* AddComplexSource() */
2739 : /************************************************************************/
2740 :
2741 30 : CPLErr VRTSourcedRasterBand::AddComplexSource(
2742 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2743 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2744 : double dfDstXSize, double dfDstYSize, double dfScaleOff,
2745 : double dfScaleRatio, double dfNoDataValueIn, int nColorTableComponent)
2746 :
2747 : {
2748 : /* -------------------------------------------------------------------- */
2749 : /* Create source. */
2750 : /* -------------------------------------------------------------------- */
2751 30 : VRTComplexSource *const poSource = new VRTComplexSource();
2752 :
2753 30 : ConfigureSource(poSource, poSrcBand, FALSE, dfSrcXOff, dfSrcYOff,
2754 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2755 : dfDstYSize);
2756 :
2757 : /* -------------------------------------------------------------------- */
2758 : /* Set complex parameters. */
2759 : /* -------------------------------------------------------------------- */
2760 30 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2761 19 : poSource->SetNoDataValue(dfNoDataValueIn);
2762 :
2763 30 : if (dfScaleOff != 0.0 || dfScaleRatio != 1.0)
2764 21 : poSource->SetLinearScaling(dfScaleOff, dfScaleRatio);
2765 :
2766 30 : poSource->SetColorTableComponent(nColorTableComponent);
2767 :
2768 : /* -------------------------------------------------------------------- */
2769 : /* add to list. */
2770 : /* -------------------------------------------------------------------- */
2771 30 : return AddSource(poSource);
2772 : }
2773 :
2774 : /*! @endcond */
2775 :
2776 : /************************************************************************/
2777 : /* VRTAddComplexSource() */
2778 : /************************************************************************/
2779 :
2780 : /**
2781 : * @see VRTSourcedRasterBand::AddComplexSource().
2782 : */
2783 :
2784 0 : CPLErr CPL_STDCALL VRTAddComplexSource(
2785 : VRTSourcedRasterBandH hVRTBand, GDALRasterBandH hSrcBand, int nSrcXOff,
2786 : int nSrcYOff, int nSrcXSize, int nSrcYSize, int nDstXOff, int nDstYOff,
2787 : int nDstXSize, int nDstYSize, double dfScaleOff, double dfScaleRatio,
2788 : double dfNoDataValue)
2789 : {
2790 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddComplexSource", CE_Failure);
2791 :
2792 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddComplexSource(
2793 : reinterpret_cast<GDALRasterBand *>(hSrcBand), nSrcXOff, nSrcYOff,
2794 : nSrcXSize, nSrcYSize, nDstXOff, nDstYOff, nDstXSize, nDstYSize,
2795 0 : dfScaleOff, dfScaleRatio, dfNoDataValue);
2796 : }
2797 :
2798 : /*! @cond Doxygen_Suppress */
2799 :
2800 : /************************************************************************/
2801 : /* AddFuncSource() */
2802 : /************************************************************************/
2803 :
2804 12 : CPLErr VRTSourcedRasterBand::AddFuncSource(VRTImageReadFunc pfnReadFunc,
2805 : void *pCBData,
2806 : double dfNoDataValueIn)
2807 :
2808 : {
2809 : /* -------------------------------------------------------------------- */
2810 : /* Create source. */
2811 : /* -------------------------------------------------------------------- */
2812 12 : VRTFuncSource *const poFuncSource = new VRTFuncSource;
2813 :
2814 12 : poFuncSource->fNoDataValue = static_cast<float>(dfNoDataValueIn);
2815 12 : poFuncSource->pfnReadFunc = pfnReadFunc;
2816 12 : poFuncSource->pCBData = pCBData;
2817 12 : poFuncSource->eType = GetRasterDataType();
2818 :
2819 : /* -------------------------------------------------------------------- */
2820 : /* add to list. */
2821 : /* -------------------------------------------------------------------- */
2822 12 : return AddSource(poFuncSource);
2823 : }
2824 :
2825 : /*! @endcond */
2826 :
2827 : /************************************************************************/
2828 : /* VRTAddFuncSource() */
2829 : /************************************************************************/
2830 :
2831 : /**
2832 : * @see VRTSourcedRasterBand::AddFuncSource().
2833 : */
2834 :
2835 0 : CPLErr CPL_STDCALL VRTAddFuncSource(VRTSourcedRasterBandH hVRTBand,
2836 : VRTImageReadFunc pfnReadFunc, void *pCBData,
2837 : double dfNoDataValue)
2838 : {
2839 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddFuncSource", CE_Failure);
2840 :
2841 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddFuncSource(
2842 0 : pfnReadFunc, pCBData, dfNoDataValue);
2843 : }
2844 :
2845 : /*! @cond Doxygen_Suppress */
2846 :
2847 : /************************************************************************/
2848 : /* GetMetadataDomainList() */
2849 : /************************************************************************/
2850 :
2851 0 : char **VRTSourcedRasterBand::GetMetadataDomainList()
2852 : {
2853 0 : return CSLAddString(GDALRasterBand::GetMetadataDomainList(),
2854 0 : "LocationInfo");
2855 : }
2856 :
2857 : /************************************************************************/
2858 : /* GetMetadataItem() */
2859 : /************************************************************************/
2860 :
2861 255845 : const char *VRTSourcedRasterBand::GetMetadataItem(const char *pszName,
2862 : const char *pszDomain)
2863 :
2864 : {
2865 : /* ==================================================================== */
2866 : /* LocationInfo handling. */
2867 : /* ==================================================================== */
2868 255845 : if (pszDomain != nullptr && EQUAL(pszDomain, "LocationInfo") &&
2869 3 : (STARTS_WITH_CI(pszName, "Pixel_") ||
2870 0 : STARTS_WITH_CI(pszName, "GeoPixel_")))
2871 : {
2872 : /* --------------------------------------------------------------------
2873 : */
2874 : /* What pixel are we aiming at? */
2875 : /* --------------------------------------------------------------------
2876 : */
2877 3 : int iPixel = 0;
2878 3 : int iLine = 0;
2879 :
2880 3 : if (STARTS_WITH_CI(pszName, "Pixel_"))
2881 : {
2882 : // TODO(schwehr): Replace sscanf.
2883 3 : if (sscanf(pszName + 6, "%d_%d", &iPixel, &iLine) != 2)
2884 0 : return nullptr;
2885 : }
2886 0 : else if (STARTS_WITH_CI(pszName, "GeoPixel_"))
2887 : {
2888 0 : const double dfGeoX = CPLAtof(pszName + 9);
2889 0 : const char *const pszUnderscore = strchr(pszName + 9, '_');
2890 0 : if (!pszUnderscore)
2891 0 : return nullptr;
2892 0 : const double dfGeoY = CPLAtof(pszUnderscore + 1);
2893 :
2894 0 : if (GetDataset() == nullptr)
2895 0 : return nullptr;
2896 :
2897 0 : GDALGeoTransform gt, invGT;
2898 0 : if (GetDataset()->GetGeoTransform(gt) != CE_None ||
2899 0 : !gt.GetInverse(invGT))
2900 : {
2901 0 : return nullptr;
2902 : }
2903 :
2904 0 : iPixel = static_cast<int>(
2905 0 : floor(invGT[0] + invGT[1] * dfGeoX + invGT[2] * dfGeoY));
2906 0 : iLine = static_cast<int>(
2907 0 : floor(invGT[3] + invGT[4] * dfGeoX + invGT[5] * dfGeoY));
2908 : }
2909 : else
2910 : {
2911 0 : return nullptr;
2912 : }
2913 :
2914 6 : if (iPixel < 0 || iLine < 0 || iPixel >= GetXSize() ||
2915 3 : iLine >= GetYSize())
2916 0 : return nullptr;
2917 :
2918 : /* --------------------------------------------------------------------
2919 : */
2920 : /* Find the file(s) at this location. */
2921 : /* --------------------------------------------------------------------
2922 : */
2923 3 : char **papszFileList = nullptr;
2924 3 : int nListSize = 0; // keep it in this scope
2925 3 : int nListMaxSize = 0; // keep it in this scope
2926 : CPLHashSet *const hSetFiles =
2927 3 : CPLHashSetNew(CPLHashSetHashStr, CPLHashSetEqualStr, nullptr);
2928 :
2929 6 : for (auto &poSource : m_papoSources)
2930 : {
2931 3 : if (!poSource->IsSimpleSource())
2932 0 : continue;
2933 :
2934 : VRTSimpleSource *const poSrc =
2935 3 : static_cast<VRTSimpleSource *>(poSource.get());
2936 :
2937 3 : double dfReqXOff = 0.0;
2938 3 : double dfReqYOff = 0.0;
2939 3 : double dfReqXSize = 0.0;
2940 3 : double dfReqYSize = 0.0;
2941 3 : int nReqXOff = 0;
2942 3 : int nReqYOff = 0;
2943 3 : int nReqXSize = 0;
2944 3 : int nReqYSize = 0;
2945 3 : int nOutXOff = 0;
2946 3 : int nOutYOff = 0;
2947 3 : int nOutXSize = 0;
2948 3 : int nOutYSize = 0;
2949 :
2950 3 : bool bError = false;
2951 3 : if (!poSrc->GetSrcDstWindow(
2952 : iPixel, iLine, 1, 1, 1, 1, GRIORA_NearestNeighbour,
2953 : &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
2954 : &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
2955 : &nOutXSize, &nOutYSize, bError))
2956 : {
2957 0 : if (bError)
2958 : {
2959 0 : CSLDestroy(papszFileList);
2960 0 : CPLHashSetDestroy(hSetFiles);
2961 0 : return nullptr;
2962 : }
2963 0 : continue;
2964 : }
2965 :
2966 3 : poSrc->GetFileList(&papszFileList, &nListSize, &nListMaxSize,
2967 3 : hSetFiles);
2968 : }
2969 :
2970 : /* --------------------------------------------------------------------
2971 : */
2972 : /* Format into XML. */
2973 : /* --------------------------------------------------------------------
2974 : */
2975 3 : m_osLastLocationInfo = "<LocationInfo>";
2976 6 : for (int i = 0; i < nListSize && papszFileList[i] != nullptr; i++)
2977 : {
2978 3 : m_osLastLocationInfo += "<File>";
2979 : char *const pszXMLEscaped =
2980 3 : CPLEscapeString(papszFileList[i], -1, CPLES_XML);
2981 3 : m_osLastLocationInfo += pszXMLEscaped;
2982 3 : CPLFree(pszXMLEscaped);
2983 3 : m_osLastLocationInfo += "</File>";
2984 : }
2985 3 : m_osLastLocationInfo += "</LocationInfo>";
2986 :
2987 3 : CSLDestroy(papszFileList);
2988 3 : CPLHashSetDestroy(hSetFiles);
2989 :
2990 3 : return m_osLastLocationInfo.c_str();
2991 : }
2992 :
2993 : /* ==================================================================== */
2994 : /* Other domains. */
2995 : /* ==================================================================== */
2996 :
2997 255842 : return GDALRasterBand::GetMetadataItem(pszName, pszDomain);
2998 : }
2999 :
3000 : /************************************************************************/
3001 : /* GetMetadata() */
3002 : /************************************************************************/
3003 :
3004 14637 : CSLConstList VRTSourcedRasterBand::GetMetadata(const char *pszDomain)
3005 :
3006 : {
3007 : /* ==================================================================== */
3008 : /* vrt_sources domain handling. */
3009 : /* ==================================================================== */
3010 14637 : if (pszDomain != nullptr && EQUAL(pszDomain, "vrt_sources"))
3011 : {
3012 1 : if (static_cast<size_t>(m_aosSourceList.size()) != m_papoSources.size())
3013 : {
3014 1 : m_aosSourceList.clear();
3015 :
3016 : // Process SimpleSources
3017 2 : for (int iSource = 0;
3018 2 : iSource < static_cast<int>(m_papoSources.size()); iSource++)
3019 : {
3020 : CPLXMLNode *const psXMLSrc =
3021 1 : m_papoSources[iSource]->SerializeToXML(nullptr);
3022 1 : if (psXMLSrc == nullptr)
3023 0 : continue;
3024 :
3025 1 : char *const pszXML = CPLSerializeXMLTree(psXMLSrc);
3026 :
3027 : m_aosSourceList.AddString(
3028 1 : CPLSPrintf("source_%d=%s", iSource, pszXML));
3029 1 : CPLFree(pszXML);
3030 1 : CPLDestroyXMLNode(psXMLSrc);
3031 : }
3032 : }
3033 1 : return m_aosSourceList.List();
3034 : }
3035 :
3036 : /* ==================================================================== */
3037 : /* Other domains. */
3038 : /* ==================================================================== */
3039 :
3040 14636 : return GDALRasterBand::GetMetadata(pszDomain);
3041 : }
3042 :
3043 : /************************************************************************/
3044 : /* SetMetadataItem() */
3045 : /************************************************************************/
3046 :
3047 133616 : CPLErr VRTSourcedRasterBand::SetMetadataItem(const char *pszName,
3048 : const char *pszValue,
3049 : const char *pszDomain)
3050 :
3051 : {
3052 : #if DEBUG_VERBOSE
3053 : CPLDebug("VRT", "VRTSourcedRasterBand::SetMetadataItem(%s,%s,%s)\n",
3054 : pszName, pszValue ? pszValue : "(null)",
3055 : pszDomain ? pszDomain : "(null)");
3056 : #endif
3057 :
3058 133616 : if (pszDomain != nullptr && EQUAL(pszDomain, "new_vrt_sources"))
3059 : {
3060 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
3061 1 : GetGDALDriverManager()->GetDriverByName("VRT"));
3062 1 : if (poDriver)
3063 : {
3064 1 : CPLXMLNode *const psTree = CPLParseXMLString(pszValue);
3065 1 : if (psTree == nullptr)
3066 0 : return CE_Failure;
3067 :
3068 1 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
3069 1 : if (l_poDS == nullptr)
3070 : {
3071 0 : CPLDestroyXMLNode(psTree);
3072 0 : return CE_Failure;
3073 : }
3074 2 : VRTSource *const poSource = poDriver->ParseSource(
3075 1 : psTree, nullptr, l_poDS->m_oMapSharedSources);
3076 1 : CPLDestroyXMLNode(psTree);
3077 :
3078 1 : if (poSource != nullptr)
3079 1 : return AddSource(poSource);
3080 : }
3081 :
3082 0 : return CE_Failure;
3083 : }
3084 133615 : else if (pszDomain != nullptr && EQUAL(pszDomain, "vrt_sources"))
3085 : {
3086 1 : int iSource = 0;
3087 : // TODO(schwehr): Replace sscanf.
3088 2 : if (sscanf(pszName, "source_%d", &iSource) != 1 || iSource < 0 ||
3089 1 : iSource >= static_cast<int>(m_papoSources.size()))
3090 : {
3091 0 : CPLError(CE_Failure, CPLE_AppDefined,
3092 : "%s metadata item name is not recognized. "
3093 : "Should be between source_0 and source_%d",
3094 0 : pszName, static_cast<int>(m_papoSources.size()) - 1);
3095 0 : return CE_Failure;
3096 : }
3097 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
3098 1 : GetGDALDriverManager()->GetDriverByName("VRT"));
3099 1 : if (poDriver)
3100 : {
3101 1 : CPLXMLNode *const psTree = CPLParseXMLString(pszValue);
3102 1 : if (psTree == nullptr)
3103 1 : return CE_Failure;
3104 :
3105 1 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
3106 1 : if (l_poDS == nullptr)
3107 : {
3108 0 : CPLDestroyXMLNode(psTree);
3109 0 : return CE_Failure;
3110 : }
3111 : auto poSource = std::unique_ptr<VRTSource>(poDriver->ParseSource(
3112 1 : psTree, nullptr, l_poDS->m_oMapSharedSources));
3113 1 : CPLDestroyXMLNode(psTree);
3114 :
3115 1 : if (poSource != nullptr)
3116 : {
3117 1 : m_papoSources[iSource] = std::move(poSource);
3118 1 : static_cast<VRTDataset *>(poDS)->SetNeedsFlush();
3119 1 : return CE_None;
3120 : }
3121 : }
3122 0 : return CE_Failure;
3123 : }
3124 :
3125 133614 : return VRTRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
3126 : }
3127 :
3128 : /************************************************************************/
3129 : /* SetMetadata() */
3130 : /************************************************************************/
3131 :
3132 70879 : CPLErr VRTSourcedRasterBand::SetMetadata(CSLConstList papszNewMD,
3133 : const char *pszDomain)
3134 :
3135 : {
3136 70879 : if (pszDomain != nullptr && (EQUAL(pszDomain, "new_vrt_sources") ||
3137 70879 : EQUAL(pszDomain, "vrt_sources")))
3138 : {
3139 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
3140 8 : GetGDALDriverManager()->GetDriverByName("VRT"));
3141 8 : if (!poDriver)
3142 0 : return CE_Failure;
3143 :
3144 8 : if (EQUAL(pszDomain, "vrt_sources"))
3145 : {
3146 8 : m_papoSources.clear();
3147 : }
3148 :
3149 6 : for (const char *const pszMDItem :
3150 20 : cpl::Iterate(CSLConstList(papszNewMD)))
3151 : {
3152 8 : const char *const pszXML = CPLParseNameValue(pszMDItem, nullptr);
3153 8 : CPLXMLTreeCloser psTree(CPLParseXMLString(pszXML));
3154 8 : if (!psTree)
3155 0 : return CE_Failure;
3156 :
3157 8 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
3158 8 : if (l_poDS == nullptr)
3159 : {
3160 0 : return CE_Failure;
3161 : }
3162 8 : VRTSource *const poSource = poDriver->ParseSource(
3163 8 : psTree.get(), nullptr, l_poDS->m_oMapSharedSources);
3164 8 : if (poSource == nullptr)
3165 2 : return CE_Failure;
3166 :
3167 6 : const CPLErr eErr = AddSource(poSource);
3168 : // cppcheck-suppress knownConditionTrueFalse
3169 6 : if (eErr != CE_None)
3170 0 : return eErr;
3171 : }
3172 :
3173 6 : return CE_None;
3174 : }
3175 :
3176 70871 : return VRTRasterBand::SetMetadata(papszNewMD, pszDomain);
3177 : }
3178 :
3179 : /************************************************************************/
3180 : /* GetFileList() */
3181 : /************************************************************************/
3182 :
3183 68 : void VRTSourcedRasterBand::GetFileList(char ***ppapszFileList, int *pnSize,
3184 : int *pnMaxSize, CPLHashSet *hSetFiles)
3185 : {
3186 179 : for (auto &poSource : m_papoSources)
3187 : {
3188 111 : poSource->GetFileList(ppapszFileList, pnSize, pnMaxSize, hSetFiles);
3189 : }
3190 :
3191 68 : VRTRasterBand::GetFileList(ppapszFileList, pnSize, pnMaxSize, hSetFiles);
3192 68 : }
3193 :
3194 : /************************************************************************/
3195 : /* CloseDependentDatasets() */
3196 : /************************************************************************/
3197 :
3198 141813 : int VRTSourcedRasterBand::CloseDependentDatasets()
3199 : {
3200 141813 : int ret = VRTRasterBand::CloseDependentDatasets();
3201 :
3202 141813 : if (m_papoSources.empty())
3203 384 : return ret;
3204 :
3205 141429 : m_papoSources.clear();
3206 :
3207 141429 : return TRUE;
3208 : }
3209 :
3210 : /************************************************************************/
3211 : /* FlushCache() */
3212 : /************************************************************************/
3213 :
3214 207966 : CPLErr VRTSourcedRasterBand::FlushCache(bool bAtClosing)
3215 : {
3216 207966 : CPLErr eErr = VRTRasterBand::FlushCache(bAtClosing);
3217 521253 : for (size_t i = 0; i < m_papoSources.size() && eErr == CE_None; i++)
3218 : {
3219 313287 : eErr = m_papoSources[i]->FlushCache(bAtClosing);
3220 : }
3221 207966 : return eErr;
3222 : }
3223 :
3224 : /************************************************************************/
3225 : /* RemoveCoveredSources() */
3226 : /************************************************************************/
3227 :
3228 : /** Remove sources that are covered by other sources.
3229 : *
3230 : * This method removes sources that are covered entirely by (one or several)
3231 : * sources of higher priority (even if they declare a nodata setting).
3232 : * This optimizes the size of the VRT and the rendering time.
3233 : */
3234 14 : void VRTSourcedRasterBand::RemoveCoveredSources(CSLConstList papszOptions)
3235 : {
3236 : #ifndef HAVE_GEOS
3237 : if (CPLTestBool(CSLFetchNameValueDef(
3238 : papszOptions, "EMIT_ERROR_IF_GEOS_NOT_AVAILABLE", "TRUE")))
3239 : {
3240 : CPLError(CE_Failure, CPLE_NotSupported,
3241 : "RemoveCoveredSources() not implemented in builds "
3242 : "without GEOS support");
3243 : }
3244 : #else
3245 : (void)papszOptions;
3246 :
3247 : CPLRectObj globalBounds;
3248 14 : globalBounds.minx = 0;
3249 14 : globalBounds.miny = 0;
3250 14 : globalBounds.maxx = nRasterXSize;
3251 14 : globalBounds.maxy = nRasterYSize;
3252 :
3253 : // Create an index with the bbox of all sources
3254 14 : CPLQuadTree *hTree = CPLQuadTreeCreate(&globalBounds, nullptr);
3255 38 : for (int i = 0; i < static_cast<int>(m_papoSources.size()); i++)
3256 : {
3257 24 : if (m_papoSources[i]->IsSimpleSource())
3258 : {
3259 : VRTSimpleSource *poSS =
3260 24 : cpl::down_cast<VRTSimpleSource *>(m_papoSources[i].get());
3261 24 : void *hFeature =
3262 24 : reinterpret_cast<void *>(static_cast<uintptr_t>(i));
3263 : CPLRectObj rect;
3264 24 : rect.minx = std::max(0.0, poSS->m_dfDstXOff);
3265 24 : rect.miny = std::max(0.0, poSS->m_dfDstYOff);
3266 48 : rect.maxx = std::min(double(nRasterXSize),
3267 24 : poSS->m_dfDstXOff + poSS->m_dfDstXSize);
3268 48 : rect.maxy = std::min(double(nRasterYSize),
3269 24 : poSS->m_dfDstYOff + poSS->m_dfDstYSize);
3270 24 : CPLQuadTreeInsertWithBounds(hTree, hFeature, &rect);
3271 : }
3272 : }
3273 :
3274 38 : for (int i = 0; i < static_cast<int>(m_papoSources.size()); i++)
3275 : {
3276 24 : if (m_papoSources[i]->IsSimpleSource())
3277 : {
3278 : VRTSimpleSource *poSS =
3279 24 : cpl::down_cast<VRTSimpleSource *>(m_papoSources[i].get());
3280 : CPLRectObj rect;
3281 24 : rect.minx = std::max(0.0, poSS->m_dfDstXOff);
3282 24 : rect.miny = std::max(0.0, poSS->m_dfDstYOff);
3283 48 : rect.maxx = std::min(double(nRasterXSize),
3284 24 : poSS->m_dfDstXOff + poSS->m_dfDstXSize);
3285 48 : rect.maxy = std::min(double(nRasterYSize),
3286 24 : poSS->m_dfDstYOff + poSS->m_dfDstYSize);
3287 :
3288 : // Find sources whose extent intersect with the current one
3289 24 : int nFeatureCount = 0;
3290 : void **pahFeatures =
3291 24 : CPLQuadTreeSearch(hTree, &rect, &nFeatureCount);
3292 :
3293 : // Compute the bounding box of those sources, only if they are
3294 : // on top of the current one
3295 : CPLRectObj rectIntersecting;
3296 24 : rectIntersecting.minx = std::numeric_limits<double>::max();
3297 24 : rectIntersecting.miny = std::numeric_limits<double>::max();
3298 24 : rectIntersecting.maxx = -std::numeric_limits<double>::max();
3299 24 : rectIntersecting.maxy = -std::numeric_limits<double>::max();
3300 61 : for (int j = 0; j < nFeatureCount; j++)
3301 : {
3302 37 : const int curFeature = static_cast<int>(
3303 37 : reinterpret_cast<uintptr_t>(pahFeatures[j]));
3304 37 : if (curFeature > i)
3305 : {
3306 : VRTSimpleSource *poOtherSS =
3307 13 : cpl::down_cast<VRTSimpleSource *>(
3308 13 : m_papoSources[curFeature].get());
3309 13 : rectIntersecting.minx =
3310 13 : std::min(rectIntersecting.minx, poOtherSS->m_dfDstXOff);
3311 13 : rectIntersecting.miny =
3312 13 : std::min(rectIntersecting.miny, poOtherSS->m_dfDstYOff);
3313 13 : rectIntersecting.maxx = std::max(
3314 : rectIntersecting.maxx,
3315 13 : poOtherSS->m_dfDstXOff + poOtherSS->m_dfDstXSize);
3316 13 : rectIntersecting.maxy = std::max(
3317 : rectIntersecting.maxy,
3318 13 : poOtherSS->m_dfDstYOff + poOtherSS->m_dfDstXSize);
3319 : }
3320 : }
3321 :
3322 : // If the boundinx box of those sources overlap the current one,
3323 : // then compute their union, and check if it contains the current
3324 : // source
3325 24 : if (rectIntersecting.minx <= rect.minx &&
3326 7 : rectIntersecting.miny <= rect.miny &&
3327 7 : rectIntersecting.maxx >= rect.maxx &&
3328 7 : rectIntersecting.maxy >= rect.maxy)
3329 : {
3330 14 : OGRPolygon oPoly;
3331 : {
3332 7 : auto poLR = new OGRLinearRing();
3333 7 : poLR->addPoint(rect.minx, rect.miny);
3334 7 : poLR->addPoint(rect.minx, rect.maxy);
3335 7 : poLR->addPoint(rect.maxx, rect.maxy);
3336 7 : poLR->addPoint(rect.maxx, rect.miny);
3337 7 : poLR->addPoint(rect.minx, rect.miny);
3338 7 : oPoly.addRingDirectly(poLR);
3339 : }
3340 :
3341 7 : std::unique_ptr<OGRGeometry> poUnion;
3342 24 : for (int j = 0; j < nFeatureCount; j++)
3343 : {
3344 17 : const int curFeature = static_cast<int>(
3345 17 : reinterpret_cast<uintptr_t>(pahFeatures[j]));
3346 17 : if (curFeature > i)
3347 : {
3348 : VRTSimpleSource *poOtherSS =
3349 10 : cpl::down_cast<VRTSimpleSource *>(
3350 10 : m_papoSources[curFeature].get());
3351 : CPLRectObj otherRect;
3352 10 : otherRect.minx = std::max(0.0, poOtherSS->m_dfDstXOff);
3353 10 : otherRect.miny = std::max(0.0, poOtherSS->m_dfDstYOff);
3354 20 : otherRect.maxx = std::min(double(nRasterXSize),
3355 20 : poOtherSS->m_dfDstXOff +
3356 10 : poOtherSS->m_dfDstXSize);
3357 20 : otherRect.maxy = std::min(double(nRasterYSize),
3358 20 : poOtherSS->m_dfDstYOff +
3359 10 : poOtherSS->m_dfDstYSize);
3360 20 : OGRPolygon oOtherPoly;
3361 : {
3362 10 : auto poLR = new OGRLinearRing();
3363 10 : poLR->addPoint(otherRect.minx, otherRect.miny);
3364 10 : poLR->addPoint(otherRect.minx, otherRect.maxy);
3365 10 : poLR->addPoint(otherRect.maxx, otherRect.maxy);
3366 10 : poLR->addPoint(otherRect.maxx, otherRect.miny);
3367 10 : poLR->addPoint(otherRect.minx, otherRect.miny);
3368 10 : oOtherPoly.addRingDirectly(poLR);
3369 : }
3370 10 : if (poUnion == nullptr)
3371 7 : poUnion.reset(oOtherPoly.clone());
3372 : else
3373 3 : poUnion.reset(oOtherPoly.Union(poUnion.get()));
3374 : }
3375 : }
3376 :
3377 7 : if (poUnion != nullptr && poUnion->Contains(&oPoly))
3378 : {
3379 : // We can remove the current source
3380 7 : m_papoSources[i].reset();
3381 : }
3382 : }
3383 24 : CPLFree(pahFeatures);
3384 :
3385 24 : void *hFeature =
3386 24 : reinterpret_cast<void *>(static_cast<uintptr_t>(i));
3387 24 : CPLQuadTreeRemove(hTree, hFeature, &rect);
3388 : }
3389 : }
3390 :
3391 : // Compact the papoSources array
3392 0 : m_papoSources.erase(std::remove_if(m_papoSources.begin(),
3393 : m_papoSources.end(),
3394 24 : [](const std::unique_ptr<VRTSource> &src)
3395 38 : { return src.get() == nullptr; }),
3396 28 : m_papoSources.end());
3397 :
3398 14 : CPLQuadTreeDestroy(hTree);
3399 : #endif
3400 14 : }
3401 :
3402 : /*! @endcond */
|
