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 2598 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn)
59 : {
60 2598 : VRTRasterBand::Initialize(poDSIn->GetRasterXSize(),
61 : poDSIn->GetRasterYSize());
62 :
63 2598 : poDS = poDSIn;
64 2598 : nBand = nBandIn;
65 2598 : }
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 139615 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn,
95 : GDALDataType eType, int nXSize,
96 : int nYSize, int nBlockXSizeIn,
97 139615 : int nBlockYSizeIn)
98 : {
99 139615 : VRTRasterBand::Initialize(nXSize, nYSize);
100 :
101 139615 : poDS = poDSIn;
102 139615 : nBand = nBandIn;
103 :
104 139615 : eDataType = eType;
105 139615 : if (nBlockXSizeIn > 0)
106 139237 : nBlockXSize = nBlockXSizeIn;
107 139615 : if (nBlockYSizeIn > 0)
108 139255 : nBlockYSize = nBlockYSizeIn;
109 139615 : }
110 :
111 : /************************************************************************/
112 : /* ~VRTSourcedRasterBand() */
113 : /************************************************************************/
114 :
115 282153 : VRTSourcedRasterBand::~VRTSourcedRasterBand()
116 :
117 : {
118 142213 : VRTSourcedRasterBand::CloseDependentDatasets();
119 282153 : }
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 289790 : 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 289790 : };
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 289790 : if (eRWFlag == GF_Read && (nXSize != nBufXSize || nYSize != nBufYSize) &&
185 579580 : !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 289776 : 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 61628 : void VRTSourcedRasterBand::InitializeOutputBuffer(void *pData, int nBufXSize,
516 : int nBufYSize,
517 : GDALDataType eBufType,
518 : GSpacing nPixelSpace,
519 : GSpacing nLineSpace) const
520 : {
521 61628 : if (nPixelSpace == GDALGetDataTypeSizeBytes(eBufType) &&
522 29865 : !(m_bNoDataValueSet && m_dfNoDataValue != 0.0) &&
523 121094 : !(m_bNoDataSetAsInt64 && m_nNoDataValueInt64 != 0) &&
524 29601 : !(m_bNoDataSetAsUInt64 && m_nNoDataValueUInt64 != 0))
525 : {
526 29600 : if (nLineSpace == nBufXSize * nPixelSpace)
527 : {
528 28520 : 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 61628 : }
575 :
576 : /************************************************************************/
577 : /* IRasterIO() */
578 : /************************************************************************/
579 :
580 289790 : 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 289790 : 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 579580 : const std::string osFctId("VRTSourcedRasterBand::IRasterIO");
594 579580 : GDALAntiRecursionGuard oGuard(osFctId);
595 289790 : if (oGuard.GetCallDepth() >= 32)
596 : {
597 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
598 0 : return CE_Failure;
599 : }
600 :
601 869370 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
602 : // Allow 2 recursion depths on the same dataset for non-nearest resampling
603 289790 : 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 289789 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
614 289789 : if (l_poDS &&
615 579219 : l_poDS->m_apoOverviews.empty() && // do not use virtual overviews
616 579578 : (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 289789 : 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 289775 : if (!SkipBufferInitialization())
684 : {
685 59020 : InitializeOutputBuffer(pData, nBufXSize, nBufYSize, eBufType,
686 : nPixelSpace, nLineSpace);
687 : }
688 :
689 : /* -------------------------------------------------------------------- */
690 : /* Overlay each source in turn over top this. */
691 : /* -------------------------------------------------------------------- */
692 289775 : CPLErr eErr = CE_None;
693 :
694 289775 : double dfXOff = nXOff;
695 289775 : double dfYOff = nYOff;
696 289775 : double dfXSize = nXSize;
697 289775 : double dfYSize = nYSize;
698 289775 : if (psExtraArg->bFloatingPointWindowValidity)
699 : {
700 309 : dfXOff = psExtraArg->dfXOff;
701 309 : dfYOff = psExtraArg->dfYOff;
702 309 : dfXSize = psExtraArg->dfXSize;
703 309 : dfYSize = psExtraArg->dfYSize;
704 : }
705 :
706 289775 : if (l_poDS)
707 289775 : l_poDS->m_bMultiThreadedRasterIOLastUsed = false;
708 :
709 289775 : int nContributingSources = 0;
710 289775 : int nMaxThreads = 0;
711 289775 : constexpr int MINIMUM_PIXEL_COUNT_FOR_THREADED_IO = 1000 * 1000;
712 579550 : if (l_poDS &&
713 289775 : (static_cast<int64_t>(nBufXSize) * nBufYSize >=
714 289442 : MINIMUM_PIXEL_COUNT_FOR_THREADED_IO ||
715 289442 : 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 579570 : 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 86 : while (oQueue->WaitEvent())
793 : {
794 : // Quite rough progress callback. We could do better by counting
795 : // the number of contributing pixels.
796 68 : if (psExtraArg->pfnProgress)
797 : {
798 110 : 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 289757 : GDALProgressFunc const pfnProgressGlobal = psExtraArg->pfnProgress;
810 289757 : void *const pProgressDataGlobal = psExtraArg->pProgressData;
811 :
812 289757 : VRTSource::WorkingState oWorkingState;
813 289757 : const int nSources = static_cast<int>(m_papoSources.size());
814 564386 : for (int iSource = 0; eErr == CE_None && iSource < nSources; iSource++)
815 : {
816 274629 : psExtraArg->pfnProgress = GDALScaledProgress;
817 549258 : psExtraArg->pProgressData = GDALCreateScaledProgress(
818 274629 : 1.0 * iSource / nSources, 1.0 * (iSource + 1) / nSources,
819 : pfnProgressGlobal, pProgressDataGlobal);
820 274629 : if (psExtraArg->pProgressData == nullptr)
821 274071 : psExtraArg->pfnProgress = nullptr;
822 :
823 549258 : eErr = m_papoSources[iSource]->RasterIO(
824 : eDataType, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
825 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg,
826 274629 : l_poDS ? l_poDS->m_oWorkingState : oWorkingState);
827 :
828 274629 : GDALDestroyScaledProgress(psExtraArg->pProgressData);
829 : }
830 :
831 289757 : psExtraArg->pfnProgress = pfnProgressGlobal;
832 289757 : psExtraArg->pProgressData = pProgressDataGlobal;
833 : }
834 :
835 289775 : if (eErr == CE_None && psExtraArg->pfnProgress)
836 : {
837 425 : psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
838 : }
839 :
840 289775 : return eErr;
841 : }
842 :
843 : /************************************************************************/
844 : /* IGetDataCoverageStatus() */
845 : /************************************************************************/
846 :
847 3039 : int VRTSourcedRasterBand::IGetDataCoverageStatus(int nXOff, int nYOff,
848 : int nXSize, int nYSize,
849 : int nMaskFlagStop,
850 : double *pdfDataPct)
851 : {
852 3039 : if (pdfDataPct)
853 8 : *pdfDataPct = -1.0;
854 :
855 : // Particular case for a single simple source covering the whole dataset
856 6063 : if (m_papoSources.size() == 1 && m_papoSources[0]->IsSimpleSource() &&
857 3024 : 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 2686 : int nStatus = 0;
905 :
906 5372 : auto poPolyNonCoveredBySources = std::make_unique<OGRPolygon>();
907 : {
908 5372 : auto poLR = std::make_unique<OGRLinearRing>();
909 2686 : poLR->addPoint(nXOff, nYOff);
910 2686 : poLR->addPoint(nXOff, nYOff + nYSize);
911 2686 : poLR->addPoint(nXOff + nXSize, nYOff + nYSize);
912 2686 : poLR->addPoint(nXOff + nXSize, nYOff);
913 2686 : poLR->addPoint(nXOff, nYOff);
914 2686 : poPolyNonCoveredBySources->addRingDirectly(poLR.release());
915 : }
916 :
917 2696 : for (auto &poSource : m_papoSources)
918 : {
919 2691 : if (!poSource->IsSimpleSource())
920 : {
921 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
922 2681 : GDAL_DATA_COVERAGE_STATUS_DATA;
923 : }
924 2687 : VRTSimpleSource *poSS = static_cast<VRTSimpleSource *>(poSource.get());
925 : // Check if the AOI is fully inside the source
926 2687 : double dfDstXOff = std::max(0.0, poSS->m_dfDstXOff);
927 2687 : double dfDstYOff = std::max(0.0, poSS->m_dfDstYOff);
928 2687 : double dfDstXSize = poSS->m_dfDstXSize;
929 2687 : double dfDstYSize = poSS->m_dfDstYSize;
930 2687 : auto l_poBand = poSS->GetRasterBand();
931 2687 : if (!l_poBand)
932 0 : continue;
933 2687 : if (dfDstXSize == -1)
934 8 : dfDstXSize = l_poBand->GetXSize() - dfDstXOff;
935 2687 : if (dfDstYSize == -1)
936 8 : dfDstYSize = l_poBand->GetYSize() - dfDstYOff;
937 :
938 2687 : if (nXOff >= dfDstXOff && nYOff >= dfDstYOff &&
939 2678 : nXOff + nXSize <= dfDstXOff + dfDstXSize &&
940 2654 : nYOff + nYSize <= dfDstYOff + dfDstYSize)
941 : {
942 2654 : if (pdfDataPct)
943 1 : *pdfDataPct = 100.0;
944 2654 : 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 930 : CPLErr VRTSourcedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
1002 : void *pImage)
1003 :
1004 : {
1005 930 : const int nPixelSize = GDALGetDataTypeSizeBytes(eDataType);
1006 :
1007 930 : int nReadXSize = 0;
1008 930 : if ((nBlockXOff + 1) * nBlockXSize > GetXSize())
1009 16 : nReadXSize = GetXSize() - nBlockXOff * nBlockXSize;
1010 : else
1011 914 : nReadXSize = nBlockXSize;
1012 :
1013 930 : int nReadYSize = 0;
1014 930 : if ((nBlockYOff + 1) * nBlockYSize > GetYSize())
1015 0 : nReadYSize = GetYSize() - nBlockYOff * nBlockYSize;
1016 : else
1017 930 : nReadYSize = nBlockYSize;
1018 :
1019 : GDALRasterIOExtraArg sExtraArg;
1020 930 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1021 :
1022 930 : return IRasterIO(
1023 930 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize, nReadXSize,
1024 : nReadYSize, pImage, nReadXSize, nReadYSize, eDataType, nPixelSize,
1025 1860 : 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 : eErr =
1503 10 : ComputeStatistics(bApproxOK, &adfMinMax[0], &adfMinMax[1],
1504 5 : nullptr, 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 :
1713 : {
1714 72 : const std::string osFctId("VRTSourcedRasterBand::ComputeStatistics");
1715 72 : GDALAntiRecursionGuard oGuard(osFctId);
1716 36 : if (oGuard.GetCallDepth() >= 32)
1717 : {
1718 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1719 0 : return CE_Failure;
1720 : }
1721 :
1722 108 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1723 36 : if (oGuard2.GetCallDepth() >= 2)
1724 : {
1725 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1726 0 : return CE_Failure;
1727 : }
1728 :
1729 : /* -------------------------------------------------------------------- */
1730 : /* If we have overview bands, use them for statistics. */
1731 : /* -------------------------------------------------------------------- */
1732 36 : if (bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews())
1733 : {
1734 : GDALRasterBand *const poBand =
1735 2 : GetRasterSampleOverview(GDALSTAT_APPROX_NUMSAMPLES);
1736 :
1737 2 : if (poBand != nullptr && poBand != this)
1738 : {
1739 2 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
1740 : CPLErr eErr;
1741 3 : if (l_poDS && !l_poDS->m_apoOverviews.empty() &&
1742 1 : dynamic_cast<VRTSourcedRasterBand *>(poBand) != nullptr)
1743 : {
1744 2 : auto apoTmpOverviews = std::move(l_poDS->m_apoOverviews);
1745 1 : l_poDS->m_apoOverviews.clear();
1746 1 : eErr = poBand->GDALRasterBand::ComputeStatistics(
1747 : TRUE, pdfMin, pdfMax, pdfMean, pdfStdDev, pfnProgress,
1748 : pProgressData);
1749 1 : l_poDS->m_apoOverviews = std::move(apoTmpOverviews);
1750 : }
1751 : else
1752 : {
1753 1 : eErr = poBand->ComputeStatistics(TRUE, pdfMin, pdfMax, pdfMean,
1754 : pdfStdDev, pfnProgress,
1755 1 : pProgressData);
1756 : }
1757 2 : if (eErr == CE_None && pdfMin && pdfMax && pdfMean && pdfStdDev)
1758 : {
1759 2 : SetMetadataItem("STATISTICS_APPROXIMATE", "YES");
1760 2 : SetMetadataItem(
1761 : "STATISTICS_VALID_PERCENT",
1762 2 : poBand->GetMetadataItem("STATISTICS_VALID_PERCENT"));
1763 2 : SetStatistics(*pdfMin, *pdfMax, *pdfMean, *pdfStdDev);
1764 : }
1765 2 : return eErr;
1766 : }
1767 : }
1768 :
1769 34 : if (IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1770 : /*bAllowMaxValAdjustment = */ false))
1771 : {
1772 28 : Context sContext;
1773 28 : sContext.bApproxOK = CPL_TO_BOOL(bApproxOK);
1774 28 : sContext.dfNoDataValue = m_dfNoDataValue;
1775 28 : sContext.bNoDataValueSet = m_bNoDataValueSet;
1776 28 : sContext.bHideNoDataValue = CPL_TO_BOOL(m_bHideNoDataValue);
1777 28 : sContext.pfnProgress = pfnProgress;
1778 28 : sContext.pProgressData = pProgressData;
1779 28 : sContext.nSources = static_cast<int>(m_papoSources.size());
1780 :
1781 : struct Job
1782 : {
1783 : Context *psContext = nullptr;
1784 : GDALRasterBand *poRasterBand = nullptr;
1785 : uint64_t nPixelCount = 0;
1786 : uint64_t nLastAddedPixels = 0;
1787 : uint64_t nValidPixels = 0;
1788 : double dfMin = 0;
1789 : double dfMax = 0;
1790 : double dfMean = 0;
1791 : double dfStdDev = 0;
1792 :
1793 58 : static int CPL_STDCALL ProgressFunc(double dfComplete,
1794 : const char *pszMessage,
1795 : void *pProgressArg)
1796 : {
1797 58 : auto psJob = static_cast<Job *>(pProgressArg);
1798 58 : auto psContext = psJob->psContext;
1799 58 : const uint64_t nNewAddedPixels =
1800 : dfComplete == 1.0
1801 58 : ? psJob->nPixelCount
1802 : : static_cast<uint64_t>(
1803 54 : dfComplete * psJob->nPixelCount + 0.5);
1804 : const auto nUpdateThreshold =
1805 116 : std::min(psContext->nTotalPixelsOfSources / 1000,
1806 58 : static_cast<uint64_t>(1000 * 1000));
1807 116 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1808 58 : psContext->nTotalIteratedPixels +=
1809 58 : (nNewAddedPixels - psJob->nLastAddedPixels);
1810 58 : psJob->nLastAddedPixels = nNewAddedPixels;
1811 58 : if (psContext->nTotalIteratedPixels ==
1812 58 : psContext->nTotalPixelsOfSources)
1813 : {
1814 2 : psContext->nLastReportedPixels =
1815 2 : psContext->nTotalIteratedPixels;
1816 2 : return psContext->pfnProgress(1.0, pszMessage,
1817 2 : psContext->pProgressData);
1818 : }
1819 56 : else if (psContext->nTotalIteratedPixels -
1820 56 : psContext->nLastReportedPixels >
1821 : nUpdateThreshold)
1822 : {
1823 48 : psContext->nLastReportedPixels =
1824 48 : psContext->nTotalIteratedPixels;
1825 96 : return psContext->pfnProgress(
1826 48 : static_cast<double>(psContext->nTotalIteratedPixels) /
1827 48 : psContext->nTotalPixelsOfSources,
1828 48 : pszMessage, psContext->pProgressData);
1829 : }
1830 8 : return 1;
1831 : }
1832 :
1833 39 : static void UpdateStats(const Job *psJob)
1834 : {
1835 39 : const auto nValidPixels = psJob->nValidPixels;
1836 39 : auto psContext = psJob->psContext;
1837 39 : if (nValidPixels > 0)
1838 : {
1839 29 : psContext->dfGlobalMin =
1840 29 : std::min(psContext->dfGlobalMin, psJob->dfMin);
1841 29 : psContext->dfGlobalMax =
1842 29 : std::max(psContext->dfGlobalMax, psJob->dfMax);
1843 : #ifdef naive_update_not_used
1844 : psContext->dfGlobalSum += nValidPixels * psJob->dfMean;
1845 : psContext->dfGlobalSumSquare +=
1846 : nValidPixels * (psJob->dfStdDev * psJob->dfStdDev +
1847 : psJob->dfMean * psJob->dfMean);
1848 : psContext->nGlobalValidPixels += nValidPixels;
1849 : #else
1850 29 : const auto nNewGlobalValidPixels =
1851 29 : psContext->nGlobalValidPixels + nValidPixels;
1852 29 : const double dfDelta =
1853 29 : psJob->dfMean - psContext->dfGlobalMean;
1854 29 : psContext->dfGlobalMean +=
1855 29 : nValidPixels * dfDelta / nNewGlobalValidPixels;
1856 29 : psContext->dfGlobalM2 +=
1857 29 : nValidPixels * psJob->dfStdDev * psJob->dfStdDev +
1858 29 : dfDelta * dfDelta * nValidPixels *
1859 29 : psContext->nGlobalValidPixels /
1860 : nNewGlobalValidPixels;
1861 29 : psContext->nGlobalValidPixels = nNewGlobalValidPixels;
1862 : #endif
1863 : }
1864 39 : int bHasNoData = FALSE;
1865 : const double dfNoDataValue =
1866 39 : psJob->poRasterBand->GetNoDataValue(&bHasNoData);
1867 11 : if (nValidPixels < psJob->nPixelCount && bHasNoData &&
1868 61 : !std::isnan(dfNoDataValue) &&
1869 11 : (!psContext->bNoDataValueSet ||
1870 9 : dfNoDataValue != psContext->dfNoDataValue))
1871 : {
1872 : const auto eBandDT =
1873 2 : psJob->poRasterBand->GetRasterDataType();
1874 : // Check that the band nodata value is in the range of the
1875 : // original raster type
1876 : GByte abyTempBuffer[2 * sizeof(double)];
1877 2 : CPLAssert(GDALGetDataTypeSizeBytes(eBandDT) <=
1878 : static_cast<int>(sizeof(abyTempBuffer)));
1879 2 : GDALCopyWords(&dfNoDataValue, GDT_Float64, 0,
1880 : &abyTempBuffer[0], eBandDT, 0, 1);
1881 2 : double dfNoDataValueAfter = dfNoDataValue;
1882 2 : GDALCopyWords(&abyTempBuffer[0], eBandDT, 0,
1883 : &dfNoDataValueAfter, GDT_Float64, 0, 1);
1884 4 : if (!std::isfinite(dfNoDataValue) ||
1885 2 : std::fabs(dfNoDataValueAfter - dfNoDataValue) < 1.0)
1886 : {
1887 2 : UpdateStatsWithConstantValue(
1888 : *psContext, dfNoDataValueAfter,
1889 2 : psJob->nPixelCount - nValidPixels);
1890 : }
1891 : }
1892 :
1893 39 : if (psContext->nSources == 1)
1894 : {
1895 13 : psContext->dfSingleSourceMin = psJob->dfMin;
1896 13 : psContext->dfSingleSourceMax = psJob->dfMax;
1897 13 : psContext->dfSingleSourceMean = psJob->dfMean;
1898 13 : psContext->dfSingleSourceStdDev = psJob->dfStdDev;
1899 : }
1900 39 : }
1901 : };
1902 :
1903 42 : const auto JobRunner = [](void *pData)
1904 : {
1905 42 : auto psJob = static_cast<Job *>(pData);
1906 42 : auto psContext = psJob->psContext;
1907 : {
1908 42 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1909 42 : if (psContext->bFallbackToBase || psContext->bFailure)
1910 0 : return;
1911 : }
1912 :
1913 42 : auto poSimpleSourceBand = psJob->poRasterBand;
1914 42 : psJob->nPixelCount =
1915 42 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
1916 42 : poSimpleSourceBand->GetYSize();
1917 :
1918 42 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1919 42 : CPLErr eErr = poSimpleSourceBand->ComputeStatistics(
1920 42 : psContext->bApproxOK, &psJob->dfMin, &psJob->dfMax,
1921 : &psJob->dfMean, &psJob->dfStdDev,
1922 48 : psContext->pfnProgress == nullptr ||
1923 6 : psContext->pfnProgress == GDALDummyProgress
1924 : ? GDALDummyProgress
1925 : : Job::ProgressFunc,
1926 42 : psJob);
1927 : const char *pszValidPercent =
1928 42 : poSimpleSourceBand->GetMetadataItem("STATISTICS_VALID_PERCENT");
1929 42 : psJob->nValidPixels =
1930 : pszValidPercent
1931 42 : ? static_cast<uint64_t>(CPLAtof(pszValidPercent) *
1932 42 : psJob->nPixelCount / 100.0)
1933 : : psJob->nPixelCount;
1934 42 : if (eErr == CE_Failure)
1935 : {
1936 20 : if (pszValidPercent != nullptr &&
1937 10 : CPLAtof(pszValidPercent) == 0.0)
1938 : {
1939 : // ok: no valid sample
1940 : }
1941 : else
1942 : {
1943 0 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1944 0 : psContext->bFailure = true;
1945 : }
1946 : }
1947 : else
1948 : {
1949 32 : int bHasNoData = FALSE;
1950 32 : CPL_IGNORE_RET_VAL(
1951 32 : psJob->poRasterBand->GetNoDataValue(&bHasNoData));
1952 32 : if (!bHasNoData && psContext->bNoDataValueSet &&
1953 9 : !psContext->bHideNoDataValue &&
1954 7 : psContext->dfNoDataValue >= psJob->dfMin &&
1955 2 : psContext->dfNoDataValue <= psJob->dfMax)
1956 : {
1957 2 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1958 2 : psJob->psContext->bFallbackToBase = true;
1959 2 : return;
1960 : }
1961 : }
1962 : };
1963 :
1964 28 : CPLWorkerThreadPool *poThreadPool = nullptr;
1965 : int nThreads =
1966 28 : m_papoSources.size() > 1
1967 28 : ? VRTDataset::GetNumThreads(dynamic_cast<VRTDataset *>(poDS))
1968 28 : : 0;
1969 28 : if (nThreads > 1024)
1970 0 : nThreads = 1024; // to please Coverity
1971 28 : if (nThreads > 1)
1972 : {
1973 : // Check that all sources refer to different datasets
1974 : // before allowing multithreaded access
1975 : // If the datasets belong to the MEM driver, check GDALDataset*
1976 : // pointer values. Otherwise use dataset name.
1977 26 : std::set<std::string> oSetDatasetNames;
1978 26 : std::set<GDALDataset *> oSetDatasetPointers;
1979 39 : for (auto &poSource : m_papoSources)
1980 : {
1981 : auto poSimpleSource =
1982 26 : static_cast<VRTSimpleSource *>(poSource.get());
1983 26 : assert(poSimpleSource);
1984 26 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1985 26 : assert(poSimpleSourceBand);
1986 26 : auto poSourceDataset = poSimpleSourceBand->GetDataset();
1987 26 : if (poSourceDataset == nullptr)
1988 : {
1989 0 : nThreads = 0;
1990 0 : break;
1991 : }
1992 26 : auto poDriver = poSourceDataset->GetDriver();
1993 26 : if (poDriver && EQUAL(poDriver->GetDescription(), "MEM"))
1994 : {
1995 26 : if (oSetDatasetPointers.find(poSourceDataset) !=
1996 52 : oSetDatasetPointers.end())
1997 : {
1998 0 : nThreads = 0;
1999 0 : break;
2000 : }
2001 26 : oSetDatasetPointers.insert(poSourceDataset);
2002 : }
2003 : else
2004 : {
2005 0 : if (oSetDatasetNames.find(
2006 0 : poSourceDataset->GetDescription()) !=
2007 0 : oSetDatasetNames.end())
2008 : {
2009 0 : nThreads = 0;
2010 0 : break;
2011 : }
2012 0 : oSetDatasetNames.insert(poSourceDataset->GetDescription());
2013 : }
2014 : }
2015 13 : if (nThreads > 1)
2016 : {
2017 13 : poThreadPool = GDALGetGlobalThreadPool(nThreads);
2018 : }
2019 : }
2020 :
2021 : // Compute total number of pixels of sources
2022 70 : for (auto &poSource : m_papoSources)
2023 : {
2024 : auto poSimpleSource =
2025 42 : static_cast<VRTSimpleSource *>(poSource.get());
2026 42 : assert(poSimpleSource);
2027 42 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
2028 42 : assert(poSimpleSourceBand);
2029 42 : sContext.nTotalPixelsOfSources +=
2030 42 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
2031 42 : poSimpleSourceBand->GetYSize();
2032 : }
2033 :
2034 28 : if (poThreadPool)
2035 : {
2036 13 : CPLDebugOnly("VRT", "ComputeStatistics(): use optimized "
2037 : "multi-threaded code path for mosaic");
2038 26 : std::vector<Job> asJobs(m_papoSources.size());
2039 26 : auto poQueue = poThreadPool->CreateJobQueue();
2040 39 : for (size_t i = 0; i < m_papoSources.size(); ++i)
2041 : {
2042 : auto poSimpleSource =
2043 26 : static_cast<VRTSimpleSource *>(m_papoSources[i].get());
2044 26 : assert(poSimpleSource);
2045 26 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
2046 26 : assert(poSimpleSourceBand);
2047 26 : asJobs[i].psContext = &sContext;
2048 26 : asJobs[i].poRasterBand = poSimpleSourceBand;
2049 26 : if (!poQueue->SubmitJob(JobRunner, &asJobs[i]))
2050 : {
2051 0 : sContext.bFailure = true;
2052 0 : break;
2053 : }
2054 : }
2055 13 : poQueue->WaitCompletion();
2056 13 : if (!(sContext.bFailure || sContext.bFallbackToBase))
2057 : {
2058 36 : for (size_t i = 0; i < m_papoSources.size(); ++i)
2059 : {
2060 24 : Job::UpdateStats(&asJobs[i]);
2061 : }
2062 : }
2063 : }
2064 : else
2065 : {
2066 15 : CPLDebugOnly(
2067 : "VRT",
2068 : "ComputeStatistics(): use optimized code path for mosaic");
2069 30 : for (auto &poSource : m_papoSources)
2070 : {
2071 : auto poSimpleSource =
2072 16 : static_cast<VRTSimpleSource *>(poSource.get());
2073 16 : assert(poSimpleSource);
2074 16 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
2075 16 : assert(poSimpleSourceBand);
2076 16 : Job sJob;
2077 16 : sJob.psContext = &sContext;
2078 16 : sJob.poRasterBand = poSimpleSourceBand;
2079 16 : JobRunner(&sJob);
2080 16 : if (sContext.bFailure || sContext.bFallbackToBase)
2081 : break;
2082 15 : Job::UpdateStats(&sJob);
2083 : }
2084 : }
2085 :
2086 28 : if (sContext.bFailure)
2087 0 : return CE_Failure;
2088 28 : if (sContext.bFallbackToBase)
2089 : {
2090 : // If the VRT band nodata value is in the [min, max] range
2091 : // of the source and that the source has no nodata value set,
2092 : // then we can't use the optimization.
2093 2 : CPLDebugOnly("VRT", "ComputeStatistics(): revert back to "
2094 : "generic case because of nodata value in range "
2095 : "of source raster");
2096 2 : return GDALRasterBand::ComputeStatistics(
2097 : bApproxOK, pdfMin, pdfMax, pdfMean, pdfStdDev, pfnProgress,
2098 2 : pProgressData);
2099 : }
2100 :
2101 26 : const uint64_t nTotalPixels =
2102 26 : static_cast<uint64_t>(nRasterXSize) * nRasterYSize;
2103 8 : if (m_bNoDataValueSet && m_bHideNoDataValue &&
2104 34 : !std::isnan(m_dfNoDataValue) && IsNoDataValueInDataTypeRange())
2105 : {
2106 1 : UpdateStatsWithConstantValue(sContext, m_dfNoDataValue,
2107 : nTotalPixels -
2108 1 : sContext.nGlobalValidPixels);
2109 : }
2110 25 : else if (!m_bNoDataValueSet)
2111 : {
2112 18 : sContext.nGlobalValidPixels = nTotalPixels;
2113 : }
2114 :
2115 : #ifdef naive_update_not_used
2116 : double dfGlobalMean =
2117 : sContext.nGlobalValidPixels > 0
2118 : ? sContext.dfGlobalSum / sContext.nGlobalValidPixels
2119 : : 0;
2120 : double dfGlobalStdDev = sContext.nGlobalValidPixels > 0
2121 : ? sqrt(sContext.dfGlobalSumSquare /
2122 : sContext.nGlobalValidPixels -
2123 : dfGlobalMean * dfGlobalMean)
2124 : : 0;
2125 : #else
2126 26 : double dfGlobalMean = sContext.dfGlobalMean;
2127 : double dfGlobalStdDev =
2128 26 : sContext.nGlobalValidPixels > 0
2129 26 : ? sqrt(sContext.dfGlobalM2 / sContext.nGlobalValidPixels)
2130 26 : : 0;
2131 : #endif
2132 39 : if (m_papoSources.size() == 1 &&
2133 13 : !sContext.bUpdateStatsWithConstantValueRun)
2134 : {
2135 13 : sContext.dfGlobalMin = sContext.dfSingleSourceMin;
2136 13 : sContext.dfGlobalMax = sContext.dfSingleSourceMax;
2137 13 : dfGlobalMean = sContext.dfSingleSourceMean;
2138 13 : dfGlobalStdDev = sContext.dfSingleSourceStdDev;
2139 : }
2140 :
2141 26 : if (sContext.nGlobalValidPixels > 0)
2142 : {
2143 23 : if (bApproxOK)
2144 : {
2145 1 : SetMetadataItem("STATISTICS_APPROXIMATE", "YES");
2146 : }
2147 22 : else if (GetMetadataItem("STATISTICS_APPROXIMATE"))
2148 : {
2149 0 : SetMetadataItem("STATISTICS_APPROXIMATE", nullptr);
2150 : }
2151 23 : SetStatistics(sContext.dfGlobalMin, sContext.dfGlobalMax,
2152 23 : dfGlobalMean, dfGlobalStdDev);
2153 : }
2154 : else
2155 : {
2156 3 : sContext.dfGlobalMin = 0.0;
2157 3 : sContext.dfGlobalMax = 0.0;
2158 : }
2159 :
2160 26 : SetValidPercent(nTotalPixels, sContext.nGlobalValidPixels);
2161 :
2162 26 : if (pdfMin)
2163 18 : *pdfMin = sContext.dfGlobalMin;
2164 26 : if (pdfMax)
2165 18 : *pdfMax = sContext.dfGlobalMax;
2166 26 : if (pdfMean)
2167 13 : *pdfMean = dfGlobalMean;
2168 26 : if (pdfStdDev)
2169 13 : *pdfStdDev = dfGlobalStdDev;
2170 :
2171 26 : if (sContext.nGlobalValidPixels == 0)
2172 : {
2173 3 : ReportError(CE_Failure, CPLE_AppDefined,
2174 : "Failed to compute statistics, no valid pixels found "
2175 : "in sampling.");
2176 : }
2177 :
2178 26 : return sContext.nGlobalValidPixels > 0 ? CE_None : CE_Failure;
2179 : }
2180 : else
2181 : {
2182 6 : return GDALRasterBand::ComputeStatistics(bApproxOK, pdfMin, pdfMax,
2183 : pdfMean, pdfStdDev,
2184 6 : pfnProgress, pProgressData);
2185 : }
2186 : }
2187 :
2188 : /************************************************************************/
2189 : /* GetHistogram() */
2190 : /************************************************************************/
2191 :
2192 7 : CPLErr VRTSourcedRasterBand::GetHistogram(double dfMin, double dfMax,
2193 : int nBuckets, GUIntBig *panHistogram,
2194 : int bIncludeOutOfRange, int bApproxOK,
2195 : GDALProgressFunc pfnProgress,
2196 : void *pProgressData)
2197 :
2198 : {
2199 : /* -------------------------------------------------------------------- */
2200 : /* If we have overviews, use them for the histogram. */
2201 : /* -------------------------------------------------------------------- */
2202 7 : if (bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews())
2203 : {
2204 : // FIXME: Should we use the most reduced overview here or use some
2205 : // minimum number of samples like GDALRasterBand::ComputeStatistics()
2206 : // does?
2207 1 : GDALRasterBand *poBand = GetRasterSampleOverview(0);
2208 :
2209 1 : if (poBand != nullptr && poBand != this)
2210 : {
2211 1 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
2212 2 : if (l_poDS && !l_poDS->m_apoOverviews.empty() &&
2213 1 : dynamic_cast<VRTSourcedRasterBand *>(poBand) != nullptr)
2214 : {
2215 1 : auto apoTmpOverviews = std::move(l_poDS->m_apoOverviews);
2216 1 : l_poDS->m_apoOverviews.clear();
2217 1 : auto eErr = poBand->GDALRasterBand::GetHistogram(
2218 : dfMin, dfMax, nBuckets, panHistogram, bIncludeOutOfRange,
2219 : bApproxOK, pfnProgress, pProgressData);
2220 1 : l_poDS->m_apoOverviews = std::move(apoTmpOverviews);
2221 1 : return eErr;
2222 : }
2223 : else
2224 : {
2225 0 : return poBand->GetHistogram(
2226 : dfMin, dfMax, nBuckets, panHistogram, bIncludeOutOfRange,
2227 0 : bApproxOK, pfnProgress, pProgressData);
2228 : }
2229 : }
2230 : }
2231 :
2232 6 : if (m_papoSources.size() != 1)
2233 2 : return VRTRasterBand::GetHistogram(dfMin, dfMax, nBuckets, panHistogram,
2234 : bIncludeOutOfRange, bApproxOK,
2235 2 : pfnProgress, pProgressData);
2236 :
2237 4 : if (pfnProgress == nullptr)
2238 4 : pfnProgress = GDALDummyProgress;
2239 :
2240 8 : const std::string osFctId("VRTSourcedRasterBand::GetHistogram");
2241 8 : GDALAntiRecursionGuard oGuard(osFctId);
2242 4 : if (oGuard.GetCallDepth() >= 32)
2243 : {
2244 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
2245 0 : return CE_Failure;
2246 : }
2247 :
2248 12 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
2249 4 : if (oGuard2.GetCallDepth() >= 2)
2250 : {
2251 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
2252 0 : return CE_Failure;
2253 : }
2254 :
2255 : /* -------------------------------------------------------------------- */
2256 : /* Try with source bands. */
2257 : /* -------------------------------------------------------------------- */
2258 8 : const CPLErr eErr = m_papoSources[0]->GetHistogram(
2259 : GetXSize(), GetYSize(), dfMin, dfMax, nBuckets, panHistogram,
2260 4 : bIncludeOutOfRange, bApproxOK, pfnProgress, pProgressData);
2261 4 : if (eErr != CE_None)
2262 : {
2263 0 : const CPLErr eErr2 = GDALRasterBand::GetHistogram(
2264 : dfMin, dfMax, nBuckets, panHistogram, bIncludeOutOfRange, bApproxOK,
2265 : pfnProgress, pProgressData);
2266 0 : return eErr2;
2267 : }
2268 :
2269 4 : SetDefaultHistogram(dfMin, dfMax, nBuckets, panHistogram);
2270 :
2271 4 : return CE_None;
2272 : }
2273 :
2274 : /************************************************************************/
2275 : /* AddSource() */
2276 : /************************************************************************/
2277 :
2278 242537 : CPLErr VRTSourcedRasterBand::AddSource(VRTSource *poNewSource)
2279 :
2280 : {
2281 242537 : return AddSource(std::unique_ptr<VRTSource>(poNewSource));
2282 : }
2283 :
2284 : /************************************************************************/
2285 : /* AddSource() */
2286 : /************************************************************************/
2287 :
2288 247690 : CPLErr VRTSourcedRasterBand::AddSource(std::unique_ptr<VRTSource> poNewSource)
2289 :
2290 : {
2291 247690 : auto l_poDS = static_cast<VRTDataset *>(poDS);
2292 247690 : l_poDS->SetNeedsFlush();
2293 247690 : l_poDS->SourceAdded();
2294 :
2295 247690 : if (poNewSource->IsSimpleSource())
2296 : {
2297 : VRTSimpleSource *poSS =
2298 247655 : static_cast<VRTSimpleSource *>(poNewSource.get());
2299 247655 : if (GetMetadataItem(GDALMD_NBITS, GDAL_MDD_IMAGE_STRUCTURE) != nullptr)
2300 : {
2301 : int nBits =
2302 47 : atoi(GetMetadataItem(GDALMD_NBITS, GDAL_MDD_IMAGE_STRUCTURE));
2303 47 : if (nBits >= 1 && nBits <= 31)
2304 : {
2305 47 : poSS->SetMaxValue(static_cast<int>((1U << nBits) - 1));
2306 : }
2307 : }
2308 : }
2309 :
2310 247690 : m_papoSources.push_back(std::move(poNewSource));
2311 :
2312 247690 : return CE_None;
2313 : }
2314 :
2315 : /*! @endcond */
2316 :
2317 : /************************************************************************/
2318 : /* VRTAddSource() */
2319 : /************************************************************************/
2320 :
2321 : /**
2322 : * @see VRTSourcedRasterBand::AddSource().
2323 : */
2324 :
2325 0 : CPLErr CPL_STDCALL VRTAddSource(VRTSourcedRasterBandH hVRTBand,
2326 : VRTSourceH hNewSource)
2327 : {
2328 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddSource", CE_Failure);
2329 :
2330 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddSource(
2331 0 : reinterpret_cast<VRTSource *>(hNewSource));
2332 : }
2333 :
2334 : /*! @cond Doxygen_Suppress */
2335 :
2336 : /************************************************************************/
2337 : /* XMLInit() */
2338 : /************************************************************************/
2339 :
2340 2545 : CPLErr VRTSourcedRasterBand::XMLInit(const CPLXMLNode *psTree,
2341 : const char *pszVRTPath,
2342 : VRTMapSharedResources &oMapSharedSources)
2343 :
2344 : {
2345 : {
2346 : const CPLErr eErr =
2347 2545 : VRTRasterBand::XMLInit(psTree, pszVRTPath, oMapSharedSources);
2348 2545 : if (eErr != CE_None)
2349 0 : return eErr;
2350 : }
2351 :
2352 : /* -------------------------------------------------------------------- */
2353 : /* Process sources. */
2354 : /* -------------------------------------------------------------------- */
2355 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
2356 2545 : GetGDALDriverManager()->GetDriverByName("VRT"));
2357 :
2358 2545 : for (const CPLXMLNode *psChild = psTree->psChild;
2359 118417 : psChild != nullptr && poDriver != nullptr; psChild = psChild->psNext)
2360 : {
2361 115894 : if (psChild->eType != CXT_Element)
2362 6673 : continue;
2363 :
2364 109221 : CPLErrorReset();
2365 : VRTSource *const poSource =
2366 109221 : poDriver->ParseSource(psChild, pszVRTPath, oMapSharedSources);
2367 109221 : if (poSource != nullptr)
2368 104162 : AddSource(poSource);
2369 5059 : else if (CPLGetLastErrorType() != CE_None)
2370 22 : return CE_Failure;
2371 : }
2372 :
2373 : /* -------------------------------------------------------------------- */
2374 : /* Done. */
2375 : /* -------------------------------------------------------------------- */
2376 : const char *pszSubclass =
2377 2523 : CPLGetXMLValue(psTree, "subclass", "VRTSourcedRasterBand");
2378 2523 : if (m_papoSources.empty() && !EQUAL(pszSubclass, "VRTDerivedRasterBand"))
2379 132 : CPLDebug("VRT", "No valid sources found for band in VRT file %s",
2380 132 : GetDataset() ? GetDataset()->GetDescription() : "");
2381 :
2382 2523 : return CE_None;
2383 : }
2384 :
2385 : /************************************************************************/
2386 : /* SerializeToXML() */
2387 : /************************************************************************/
2388 :
2389 695 : CPLXMLNode *VRTSourcedRasterBand::SerializeToXML(const char *pszVRTPath,
2390 : bool &bHasWarnedAboutRAMUsage,
2391 : size_t &nAccRAMUsage)
2392 :
2393 : {
2394 695 : CPLXMLNode *psTree = VRTRasterBand::SerializeToXML(
2395 : pszVRTPath, bHasWarnedAboutRAMUsage, nAccRAMUsage);
2396 695 : CPLXMLNode *psLastChild = psTree->psChild;
2397 2504 : while (psLastChild != nullptr && psLastChild->psNext != nullptr)
2398 1809 : psLastChild = psLastChild->psNext;
2399 :
2400 : /* -------------------------------------------------------------------- */
2401 : /* Process Sources. */
2402 : /* -------------------------------------------------------------------- */
2403 :
2404 695 : GIntBig nUsableRAM = -1;
2405 :
2406 3364 : for (const auto &poSource : m_papoSources)
2407 : {
2408 2669 : CPLXMLNode *const psXMLSrc = poSource->SerializeToXML(pszVRTPath);
2409 :
2410 2669 : if (psXMLSrc == nullptr)
2411 0 : break;
2412 :
2413 : // Creating the CPLXMLNode tree representation of a VRT can easily
2414 : // take several times RAM usage than its string serialization, or its
2415 : // internal representation in the driver.
2416 : // We multiply the estimate by a factor of 2, experimentally found to
2417 : // be more realistic than the conservative raw estimate.
2418 2669 : nAccRAMUsage += 2 * CPLXMLNodeGetRAMUsageEstimate(psXMLSrc);
2419 2669 : if (!bHasWarnedAboutRAMUsage && nAccRAMUsage > 512 * 1024 * 1024)
2420 : {
2421 0 : if (nUsableRAM < 0)
2422 0 : nUsableRAM = CPLGetUsablePhysicalRAM();
2423 0 : if (nUsableRAM > 0 &&
2424 0 : nAccRAMUsage > static_cast<uint64_t>(nUsableRAM) / 10 * 8)
2425 : {
2426 0 : bHasWarnedAboutRAMUsage = true;
2427 0 : CPLError(CE_Warning, CPLE_AppDefined,
2428 : "Serialization of this VRT file has already consumed "
2429 : "at least %.02f GB of RAM over a total of %.02f. This "
2430 : "process may abort",
2431 0 : double(nAccRAMUsage) / (1024 * 1024 * 1024),
2432 0 : double(nUsableRAM) / (1024 * 1024 * 1024));
2433 : }
2434 : }
2435 :
2436 2669 : if (psLastChild == nullptr)
2437 0 : psTree->psChild = psXMLSrc;
2438 : else
2439 2669 : psLastChild->psNext = psXMLSrc;
2440 2669 : psLastChild = psXMLSrc;
2441 : }
2442 :
2443 695 : return psTree;
2444 : }
2445 :
2446 : /************************************************************************/
2447 : /* SkipBufferInitialization() */
2448 : /************************************************************************/
2449 :
2450 289775 : bool VRTSourcedRasterBand::SkipBufferInitialization()
2451 : {
2452 289775 : if (m_nSkipBufferInitialization >= 0)
2453 284490 : return m_nSkipBufferInitialization != 0;
2454 : /* -------------------------------------------------------------------- */
2455 : /* Check if we can avoid buffer initialization. */
2456 : /* -------------------------------------------------------------------- */
2457 :
2458 : // Note: if one day we do alpha compositing, we will need to check that.
2459 5285 : m_nSkipBufferInitialization = FALSE;
2460 5285 : if (m_papoSources.size() != 1 || !m_papoSources[0]->IsSimpleSource())
2461 : {
2462 3558 : return false;
2463 : }
2464 : VRTSimpleSource *poSS =
2465 1727 : static_cast<VRTSimpleSource *>(m_papoSources[0].get());
2466 1727 : if (poSS->GetType() == VRTSimpleSource::GetTypeStatic())
2467 : {
2468 1315 : auto l_poBand = poSS->GetRasterBand();
2469 1304 : if (l_poBand != nullptr && poSS->m_dfSrcXOff >= 0.0 &&
2470 1266 : poSS->m_dfSrcYOff >= 0.0 &&
2471 1264 : poSS->m_dfSrcXOff + poSS->m_dfSrcXSize <= l_poBand->GetXSize() &&
2472 1262 : poSS->m_dfSrcYOff + poSS->m_dfSrcYSize <= l_poBand->GetYSize() &&
2473 1262 : poSS->m_dfDstXOff <= 0.0 && poSS->m_dfDstYOff <= 0.0 &&
2474 3802 : poSS->m_dfDstXOff + poSS->m_dfDstXSize >= nRasterXSize &&
2475 1183 : poSS->m_dfDstYOff + poSS->m_dfDstYSize >= nRasterYSize)
2476 : {
2477 1181 : m_nSkipBufferInitialization = TRUE;
2478 : }
2479 : }
2480 1727 : return m_nSkipBufferInitialization != 0;
2481 : }
2482 :
2483 : /************************************************************************/
2484 : /* ConfigureSource() */
2485 : /************************************************************************/
2486 :
2487 142916 : void VRTSourcedRasterBand::ConfigureSource(VRTSimpleSource *poSimpleSource,
2488 : GDALRasterBand *poSrcBand,
2489 : int bAddAsMaskBand, double dfSrcXOff,
2490 : double dfSrcYOff, double dfSrcXSize,
2491 : double dfSrcYSize, double dfDstXOff,
2492 : double dfDstYOff, double dfDstXSize,
2493 : double dfDstYSize)
2494 : {
2495 : /* -------------------------------------------------------------------- */
2496 : /* Default source and dest rectangles. */
2497 : /* -------------------------------------------------------------------- */
2498 142916 : if (dfSrcYSize == -1)
2499 : {
2500 66966 : dfSrcXOff = 0;
2501 66966 : dfSrcYOff = 0;
2502 66966 : dfSrcXSize = poSrcBand->GetXSize();
2503 66966 : dfSrcYSize = poSrcBand->GetYSize();
2504 : }
2505 :
2506 142916 : if (dfDstYSize == -1)
2507 : {
2508 66966 : dfDstXOff = 0;
2509 66966 : dfDstYOff = 0;
2510 66966 : dfDstXSize = nRasterXSize;
2511 66966 : dfDstYSize = nRasterYSize;
2512 : }
2513 :
2514 142916 : if (bAddAsMaskBand)
2515 53 : poSimpleSource->SetSrcMaskBand(poSrcBand);
2516 : else
2517 142863 : poSimpleSource->SetSrcBand(poSrcBand);
2518 :
2519 142916 : poSimpleSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2520 142916 : poSimpleSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2521 :
2522 : /* -------------------------------------------------------------------- */
2523 : /* If we can get the associated GDALDataset, add a reference to it.*/
2524 : /* -------------------------------------------------------------------- */
2525 142916 : GDALDataset *poSrcBandDataset = poSrcBand->GetDataset();
2526 142916 : if (poSrcBandDataset != nullptr)
2527 : {
2528 : VRTDataset *poVRTSrcBandDataset =
2529 142916 : dynamic_cast<VRTDataset *>(poSrcBandDataset);
2530 142916 : if (poVRTSrcBandDataset && !poVRTSrcBandDataset->m_bCanTakeRef)
2531 : {
2532 : // Situation triggered by VRTDataset::AddVirtualOverview()
2533 : // We create an overview dataset that is a VRT of a reduction of
2534 : // ourselves. But we don't want to take a reference on ourselves,
2535 : // otherwise this will prevent us to be closed in number of
2536 : // circumstances
2537 65 : poSimpleSource->m_bDropRefOnSrcBand = false;
2538 : }
2539 : else
2540 : {
2541 142851 : poSrcBandDataset->Reference();
2542 : }
2543 : }
2544 142916 : }
2545 :
2546 : /************************************************************************/
2547 : /* AddSimpleSource() */
2548 : /************************************************************************/
2549 :
2550 377 : CPLErr VRTSourcedRasterBand::AddSimpleSource(
2551 : const char *pszFilename, int nBandIn, double dfSrcXOff, double dfSrcYOff,
2552 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2553 : double dfDstXSize, double dfDstYSize, const char *pszResampling,
2554 : double dfNoDataValueIn)
2555 :
2556 : {
2557 : /* -------------------------------------------------------------------- */
2558 : /* Create source. */
2559 : /* -------------------------------------------------------------------- */
2560 377 : VRTSimpleSource *poSimpleSource = nullptr;
2561 :
2562 377 : if (pszResampling != nullptr && STARTS_WITH_CI(pszResampling, "aver"))
2563 : {
2564 0 : auto poAveragedSource = new VRTAveragedSource();
2565 0 : poSimpleSource = poAveragedSource;
2566 0 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2567 0 : poAveragedSource->SetNoDataValue(dfNoDataValueIn);
2568 : }
2569 : else
2570 : {
2571 377 : poSimpleSource = new VRTSimpleSource();
2572 377 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2573 0 : CPLError(
2574 : CE_Warning, CPLE_AppDefined,
2575 : "NODATA setting not currently supported for nearest "
2576 : "neighbour sampled simple sources on Virtual Datasources.");
2577 : }
2578 :
2579 377 : poSimpleSource->SetSrcBand(pszFilename, nBandIn);
2580 :
2581 377 : poSimpleSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2582 377 : poSimpleSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2583 :
2584 : /* -------------------------------------------------------------------- */
2585 : /* add to list. */
2586 : /* -------------------------------------------------------------------- */
2587 377 : return AddSource(poSimpleSource);
2588 : }
2589 :
2590 : /************************************************************************/
2591 : /* AddSimpleSource() */
2592 : /************************************************************************/
2593 :
2594 68431 : CPLErr VRTSourcedRasterBand::AddSimpleSource(
2595 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2596 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2597 : double dfDstXSize, double dfDstYSize, const char *pszResampling,
2598 : double dfNoDataValueIn)
2599 :
2600 : {
2601 : /* -------------------------------------------------------------------- */
2602 : /* Create source. */
2603 : /* -------------------------------------------------------------------- */
2604 68431 : VRTSimpleSource *poSimpleSource = nullptr;
2605 :
2606 68431 : if (pszResampling != nullptr && STARTS_WITH_CI(pszResampling, "aver"))
2607 : {
2608 0 : auto poAveragedSource = new VRTAveragedSource();
2609 0 : poSimpleSource = poAveragedSource;
2610 0 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2611 0 : poAveragedSource->SetNoDataValue(dfNoDataValueIn);
2612 : }
2613 : else
2614 : {
2615 68431 : poSimpleSource = new VRTSimpleSource();
2616 68431 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2617 0 : CPLError(
2618 : CE_Warning, CPLE_AppDefined,
2619 : "NODATA setting not currently supported for "
2620 : "neighbour sampled simple sources on Virtual Datasources.");
2621 : }
2622 :
2623 68431 : ConfigureSource(poSimpleSource, poSrcBand, FALSE, dfSrcXOff, dfSrcYOff,
2624 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2625 : dfDstYSize);
2626 :
2627 : /* -------------------------------------------------------------------- */
2628 : /* add to list. */
2629 : /* -------------------------------------------------------------------- */
2630 68431 : return AddSource(poSimpleSource);
2631 : }
2632 :
2633 : /************************************************************************/
2634 : /* AddMaskBandSource() */
2635 : /************************************************************************/
2636 :
2637 : // poSrcBand is not the mask band, but the band from which the mask band is
2638 : // taken.
2639 37 : CPLErr VRTSourcedRasterBand::AddMaskBandSource(
2640 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2641 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2642 : double dfDstXSize, double dfDstYSize)
2643 : {
2644 : /* -------------------------------------------------------------------- */
2645 : /* Create source. */
2646 : /* -------------------------------------------------------------------- */
2647 37 : VRTSimpleSource *poSimpleSource = new VRTSimpleSource();
2648 :
2649 37 : ConfigureSource(poSimpleSource, poSrcBand, TRUE, dfSrcXOff, dfSrcYOff,
2650 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2651 : dfDstYSize);
2652 :
2653 : /* -------------------------------------------------------------------- */
2654 : /* add to list. */
2655 : /* -------------------------------------------------------------------- */
2656 37 : return AddSource(poSimpleSource);
2657 : }
2658 :
2659 : /*! @endcond */
2660 :
2661 : /************************************************************************/
2662 : /* VRTAddSimpleSource() */
2663 : /************************************************************************/
2664 :
2665 : /**
2666 : * @see VRTSourcedRasterBand::AddSimpleSource().
2667 : */
2668 :
2669 968 : CPLErr CPL_STDCALL VRTAddSimpleSource(VRTSourcedRasterBandH hVRTBand,
2670 : GDALRasterBandH hSrcBand, int nSrcXOff,
2671 : int nSrcYOff, int nSrcXSize,
2672 : int nSrcYSize, int nDstXOff, int nDstYOff,
2673 : int nDstXSize, int nDstYSize,
2674 : const char *pszResampling,
2675 : double dfNoDataValue)
2676 : {
2677 968 : VALIDATE_POINTER1(hVRTBand, "VRTAddSimpleSource", CE_Failure);
2678 :
2679 968 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddSimpleSource(
2680 : reinterpret_cast<GDALRasterBand *>(hSrcBand), nSrcXOff, nSrcYOff,
2681 : nSrcXSize, nSrcYSize, nDstXOff, nDstYOff, nDstXSize, nDstYSize,
2682 968 : pszResampling, dfNoDataValue);
2683 : }
2684 :
2685 : /*! @cond Doxygen_Suppress */
2686 :
2687 : /************************************************************************/
2688 : /* AddComplexSource() */
2689 : /************************************************************************/
2690 :
2691 29 : CPLErr VRTSourcedRasterBand::AddComplexSource(
2692 : const char *pszFilename, int nBandIn, double dfSrcXOff, double dfSrcYOff,
2693 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2694 : double dfDstXSize, double dfDstYSize, double dfScaleOff,
2695 : double dfScaleRatio, double dfNoDataValueIn, int nColorTableComponent)
2696 :
2697 : {
2698 : /* -------------------------------------------------------------------- */
2699 : /* Create source. */
2700 : /* -------------------------------------------------------------------- */
2701 29 : VRTComplexSource *const poSource = new VRTComplexSource();
2702 :
2703 29 : poSource->SetSrcBand(pszFilename, nBandIn);
2704 :
2705 29 : poSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2706 29 : poSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2707 :
2708 : /* -------------------------------------------------------------------- */
2709 : /* Set complex parameters. */
2710 : /* -------------------------------------------------------------------- */
2711 29 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2712 8 : poSource->SetNoDataValue(dfNoDataValueIn);
2713 :
2714 29 : if (dfScaleOff != 0.0 || dfScaleRatio != 1.0)
2715 4 : poSource->SetLinearScaling(dfScaleOff, dfScaleRatio);
2716 :
2717 29 : poSource->SetColorTableComponent(nColorTableComponent);
2718 :
2719 : /* -------------------------------------------------------------------- */
2720 : /* add to list. */
2721 : /* -------------------------------------------------------------------- */
2722 29 : return AddSource(poSource);
2723 : }
2724 :
2725 : /************************************************************************/
2726 : /* AddComplexSource() */
2727 : /************************************************************************/
2728 :
2729 30 : CPLErr VRTSourcedRasterBand::AddComplexSource(
2730 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2731 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2732 : double dfDstXSize, double dfDstYSize, double dfScaleOff,
2733 : double dfScaleRatio, double dfNoDataValueIn, int nColorTableComponent)
2734 :
2735 : {
2736 : /* -------------------------------------------------------------------- */
2737 : /* Create source. */
2738 : /* -------------------------------------------------------------------- */
2739 30 : VRTComplexSource *const poSource = new VRTComplexSource();
2740 :
2741 30 : ConfigureSource(poSource, poSrcBand, FALSE, dfSrcXOff, dfSrcYOff,
2742 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2743 : dfDstYSize);
2744 :
2745 : /* -------------------------------------------------------------------- */
2746 : /* Set complex parameters. */
2747 : /* -------------------------------------------------------------------- */
2748 30 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2749 19 : poSource->SetNoDataValue(dfNoDataValueIn);
2750 :
2751 30 : if (dfScaleOff != 0.0 || dfScaleRatio != 1.0)
2752 21 : poSource->SetLinearScaling(dfScaleOff, dfScaleRatio);
2753 :
2754 30 : poSource->SetColorTableComponent(nColorTableComponent);
2755 :
2756 : /* -------------------------------------------------------------------- */
2757 : /* add to list. */
2758 : /* -------------------------------------------------------------------- */
2759 30 : return AddSource(poSource);
2760 : }
2761 :
2762 : /*! @endcond */
2763 :
2764 : /************************************************************************/
2765 : /* VRTAddComplexSource() */
2766 : /************************************************************************/
2767 :
2768 : /**
2769 : * @see VRTSourcedRasterBand::AddComplexSource().
2770 : */
2771 :
2772 0 : CPLErr CPL_STDCALL VRTAddComplexSource(
2773 : VRTSourcedRasterBandH hVRTBand, GDALRasterBandH hSrcBand, int nSrcXOff,
2774 : int nSrcYOff, int nSrcXSize, int nSrcYSize, int nDstXOff, int nDstYOff,
2775 : int nDstXSize, int nDstYSize, double dfScaleOff, double dfScaleRatio,
2776 : double dfNoDataValue)
2777 : {
2778 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddComplexSource", CE_Failure);
2779 :
2780 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddComplexSource(
2781 : reinterpret_cast<GDALRasterBand *>(hSrcBand), nSrcXOff, nSrcYOff,
2782 : nSrcXSize, nSrcYSize, nDstXOff, nDstYOff, nDstXSize, nDstYSize,
2783 0 : dfScaleOff, dfScaleRatio, dfNoDataValue);
2784 : }
2785 :
2786 : /*! @cond Doxygen_Suppress */
2787 :
2788 : /************************************************************************/
2789 : /* AddFuncSource() */
2790 : /************************************************************************/
2791 :
2792 12 : CPLErr VRTSourcedRasterBand::AddFuncSource(VRTImageReadFunc pfnReadFunc,
2793 : void *pCBData,
2794 : double dfNoDataValueIn)
2795 :
2796 : {
2797 : /* -------------------------------------------------------------------- */
2798 : /* Create source. */
2799 : /* -------------------------------------------------------------------- */
2800 12 : VRTFuncSource *const poFuncSource = new VRTFuncSource;
2801 :
2802 12 : poFuncSource->fNoDataValue = static_cast<float>(dfNoDataValueIn);
2803 12 : poFuncSource->pfnReadFunc = pfnReadFunc;
2804 12 : poFuncSource->pCBData = pCBData;
2805 12 : poFuncSource->eType = GetRasterDataType();
2806 :
2807 : /* -------------------------------------------------------------------- */
2808 : /* add to list. */
2809 : /* -------------------------------------------------------------------- */
2810 12 : return AddSource(poFuncSource);
2811 : }
2812 :
2813 : /*! @endcond */
2814 :
2815 : /************************************************************************/
2816 : /* VRTAddFuncSource() */
2817 : /************************************************************************/
2818 :
2819 : /**
2820 : * @see VRTSourcedRasterBand::AddFuncSource().
2821 : */
2822 :
2823 0 : CPLErr CPL_STDCALL VRTAddFuncSource(VRTSourcedRasterBandH hVRTBand,
2824 : VRTImageReadFunc pfnReadFunc, void *pCBData,
2825 : double dfNoDataValue)
2826 : {
2827 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddFuncSource", CE_Failure);
2828 :
2829 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddFuncSource(
2830 0 : pfnReadFunc, pCBData, dfNoDataValue);
2831 : }
2832 :
2833 : /*! @cond Doxygen_Suppress */
2834 :
2835 : /************************************************************************/
2836 : /* GetMetadataDomainList() */
2837 : /************************************************************************/
2838 :
2839 0 : char **VRTSourcedRasterBand::GetMetadataDomainList()
2840 : {
2841 0 : return CSLAddString(GDALRasterBand::GetMetadataDomainList(),
2842 0 : "LocationInfo");
2843 : }
2844 :
2845 : /************************************************************************/
2846 : /* GetMetadataItem() */
2847 : /************************************************************************/
2848 :
2849 256415 : const char *VRTSourcedRasterBand::GetMetadataItem(const char *pszName,
2850 : const char *pszDomain)
2851 :
2852 : {
2853 : /* ==================================================================== */
2854 : /* LocationInfo handling. */
2855 : /* ==================================================================== */
2856 256415 : if (pszDomain != nullptr && EQUAL(pszDomain, "LocationInfo") &&
2857 3 : (STARTS_WITH_CI(pszName, "Pixel_") ||
2858 0 : STARTS_WITH_CI(pszName, "GeoPixel_")))
2859 : {
2860 : /* --------------------------------------------------------------------
2861 : */
2862 : /* What pixel are we aiming at? */
2863 : /* --------------------------------------------------------------------
2864 : */
2865 3 : int iPixel = 0;
2866 3 : int iLine = 0;
2867 :
2868 3 : if (STARTS_WITH_CI(pszName, "Pixel_"))
2869 : {
2870 : // TODO(schwehr): Replace sscanf.
2871 3 : if (sscanf(pszName + 6, "%d_%d", &iPixel, &iLine) != 2)
2872 0 : return nullptr;
2873 : }
2874 0 : else if (STARTS_WITH_CI(pszName, "GeoPixel_"))
2875 : {
2876 0 : const double dfGeoX = CPLAtof(pszName + 9);
2877 0 : const char *const pszUnderscore = strchr(pszName + 9, '_');
2878 0 : if (!pszUnderscore)
2879 0 : return nullptr;
2880 0 : const double dfGeoY = CPLAtof(pszUnderscore + 1);
2881 :
2882 0 : if (GetDataset() == nullptr)
2883 0 : return nullptr;
2884 :
2885 0 : GDALGeoTransform gt, invGT;
2886 0 : if (GetDataset()->GetGeoTransform(gt) != CE_None ||
2887 0 : !gt.GetInverse(invGT))
2888 : {
2889 0 : return nullptr;
2890 : }
2891 :
2892 0 : iPixel = static_cast<int>(
2893 0 : floor(invGT[0] + invGT[1] * dfGeoX + invGT[2] * dfGeoY));
2894 0 : iLine = static_cast<int>(
2895 0 : floor(invGT[3] + invGT[4] * dfGeoX + invGT[5] * dfGeoY));
2896 : }
2897 : else
2898 : {
2899 0 : return nullptr;
2900 : }
2901 :
2902 6 : if (iPixel < 0 || iLine < 0 || iPixel >= GetXSize() ||
2903 3 : iLine >= GetYSize())
2904 0 : return nullptr;
2905 :
2906 : /* --------------------------------------------------------------------
2907 : */
2908 : /* Find the file(s) at this location. */
2909 : /* --------------------------------------------------------------------
2910 : */
2911 3 : char **papszFileList = nullptr;
2912 3 : int nListSize = 0; // keep it in this scope
2913 3 : int nListMaxSize = 0; // keep it in this scope
2914 : CPLHashSet *const hSetFiles =
2915 3 : CPLHashSetNew(CPLHashSetHashStr, CPLHashSetEqualStr, nullptr);
2916 :
2917 6 : for (auto &poSource : m_papoSources)
2918 : {
2919 3 : if (!poSource->IsSimpleSource())
2920 0 : continue;
2921 :
2922 : VRTSimpleSource *const poSrc =
2923 3 : static_cast<VRTSimpleSource *>(poSource.get());
2924 :
2925 3 : double dfReqXOff = 0.0;
2926 3 : double dfReqYOff = 0.0;
2927 3 : double dfReqXSize = 0.0;
2928 3 : double dfReqYSize = 0.0;
2929 3 : int nReqXOff = 0;
2930 3 : int nReqYOff = 0;
2931 3 : int nReqXSize = 0;
2932 3 : int nReqYSize = 0;
2933 3 : int nOutXOff = 0;
2934 3 : int nOutYOff = 0;
2935 3 : int nOutXSize = 0;
2936 3 : int nOutYSize = 0;
2937 :
2938 3 : bool bError = false;
2939 3 : if (!poSrc->GetSrcDstWindow(
2940 : iPixel, iLine, 1, 1, 1, 1, GRIORA_NearestNeighbour,
2941 : &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
2942 : &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
2943 : &nOutXSize, &nOutYSize, bError))
2944 : {
2945 0 : if (bError)
2946 : {
2947 0 : CSLDestroy(papszFileList);
2948 0 : CPLHashSetDestroy(hSetFiles);
2949 0 : return nullptr;
2950 : }
2951 0 : continue;
2952 : }
2953 :
2954 3 : poSrc->GetFileList(&papszFileList, &nListSize, &nListMaxSize,
2955 3 : hSetFiles);
2956 : }
2957 :
2958 : /* --------------------------------------------------------------------
2959 : */
2960 : /* Format into XML. */
2961 : /* --------------------------------------------------------------------
2962 : */
2963 3 : m_osLastLocationInfo = "<LocationInfo>";
2964 6 : for (int i = 0; i < nListSize && papszFileList[i] != nullptr; i++)
2965 : {
2966 3 : m_osLastLocationInfo += "<File>";
2967 : char *const pszXMLEscaped =
2968 3 : CPLEscapeString(papszFileList[i], -1, CPLES_XML);
2969 3 : m_osLastLocationInfo += pszXMLEscaped;
2970 3 : CPLFree(pszXMLEscaped);
2971 3 : m_osLastLocationInfo += "</File>";
2972 : }
2973 3 : m_osLastLocationInfo += "</LocationInfo>";
2974 :
2975 3 : CSLDestroy(papszFileList);
2976 3 : CPLHashSetDestroy(hSetFiles);
2977 :
2978 3 : return m_osLastLocationInfo.c_str();
2979 : }
2980 :
2981 : /* ==================================================================== */
2982 : /* Other domains. */
2983 : /* ==================================================================== */
2984 :
2985 256412 : return GDALRasterBand::GetMetadataItem(pszName, pszDomain);
2986 : }
2987 :
2988 : /************************************************************************/
2989 : /* GetMetadata() */
2990 : /************************************************************************/
2991 :
2992 14627 : CSLConstList VRTSourcedRasterBand::GetMetadata(const char *pszDomain)
2993 :
2994 : {
2995 : /* ==================================================================== */
2996 : /* vrt_sources domain handling. */
2997 : /* ==================================================================== */
2998 14627 : if (pszDomain != nullptr && EQUAL(pszDomain, "vrt_sources"))
2999 : {
3000 1 : if (static_cast<size_t>(m_aosSourceList.size()) != m_papoSources.size())
3001 : {
3002 1 : m_aosSourceList.clear();
3003 :
3004 : // Process SimpleSources
3005 2 : for (int iSource = 0;
3006 2 : iSource < static_cast<int>(m_papoSources.size()); iSource++)
3007 : {
3008 : CPLXMLNode *const psXMLSrc =
3009 1 : m_papoSources[iSource]->SerializeToXML(nullptr);
3010 1 : if (psXMLSrc == nullptr)
3011 0 : continue;
3012 :
3013 1 : char *const pszXML = CPLSerializeXMLTree(psXMLSrc);
3014 :
3015 : m_aosSourceList.AddString(
3016 1 : CPLSPrintf("source_%d=%s", iSource, pszXML));
3017 1 : CPLFree(pszXML);
3018 1 : CPLDestroyXMLNode(psXMLSrc);
3019 : }
3020 : }
3021 1 : return m_aosSourceList.List();
3022 : }
3023 :
3024 : /* ==================================================================== */
3025 : /* Other domains. */
3026 : /* ==================================================================== */
3027 :
3028 14626 : return GDALRasterBand::GetMetadata(pszDomain);
3029 : }
3030 :
3031 : /************************************************************************/
3032 : /* SetMetadataItem() */
3033 : /************************************************************************/
3034 :
3035 133616 : CPLErr VRTSourcedRasterBand::SetMetadataItem(const char *pszName,
3036 : const char *pszValue,
3037 : const char *pszDomain)
3038 :
3039 : {
3040 : #if DEBUG_VERBOSE
3041 : CPLDebug("VRT", "VRTSourcedRasterBand::SetMetadataItem(%s,%s,%s)\n",
3042 : pszName, pszValue ? pszValue : "(null)",
3043 : pszDomain ? pszDomain : "(null)");
3044 : #endif
3045 :
3046 133616 : if (pszDomain != nullptr && EQUAL(pszDomain, "new_vrt_sources"))
3047 : {
3048 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
3049 1 : GetGDALDriverManager()->GetDriverByName("VRT"));
3050 1 : if (poDriver)
3051 : {
3052 1 : CPLXMLNode *const psTree = CPLParseXMLString(pszValue);
3053 1 : if (psTree == nullptr)
3054 0 : return CE_Failure;
3055 :
3056 1 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
3057 1 : if (l_poDS == nullptr)
3058 : {
3059 0 : CPLDestroyXMLNode(psTree);
3060 0 : return CE_Failure;
3061 : }
3062 2 : VRTSource *const poSource = poDriver->ParseSource(
3063 1 : psTree, nullptr, l_poDS->m_oMapSharedSources);
3064 1 : CPLDestroyXMLNode(psTree);
3065 :
3066 1 : if (poSource != nullptr)
3067 1 : return AddSource(poSource);
3068 : }
3069 :
3070 0 : return CE_Failure;
3071 : }
3072 133615 : else if (pszDomain != nullptr && EQUAL(pszDomain, "vrt_sources"))
3073 : {
3074 1 : int iSource = 0;
3075 : // TODO(schwehr): Replace sscanf.
3076 2 : if (sscanf(pszName, "source_%d", &iSource) != 1 || iSource < 0 ||
3077 1 : iSource >= static_cast<int>(m_papoSources.size()))
3078 : {
3079 0 : CPLError(CE_Failure, CPLE_AppDefined,
3080 : "%s metadata item name is not recognized. "
3081 : "Should be between source_0 and source_%d",
3082 0 : pszName, static_cast<int>(m_papoSources.size()) - 1);
3083 0 : return CE_Failure;
3084 : }
3085 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
3086 1 : GetGDALDriverManager()->GetDriverByName("VRT"));
3087 1 : if (poDriver)
3088 : {
3089 1 : CPLXMLNode *const psTree = CPLParseXMLString(pszValue);
3090 1 : if (psTree == nullptr)
3091 1 : return CE_Failure;
3092 :
3093 1 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
3094 1 : if (l_poDS == nullptr)
3095 : {
3096 0 : CPLDestroyXMLNode(psTree);
3097 0 : return CE_Failure;
3098 : }
3099 : auto poSource = std::unique_ptr<VRTSource>(poDriver->ParseSource(
3100 1 : psTree, nullptr, l_poDS->m_oMapSharedSources));
3101 1 : CPLDestroyXMLNode(psTree);
3102 :
3103 1 : if (poSource != nullptr)
3104 : {
3105 1 : m_papoSources[iSource] = std::move(poSource);
3106 1 : static_cast<VRTDataset *>(poDS)->SetNeedsFlush();
3107 1 : return CE_None;
3108 : }
3109 : }
3110 0 : return CE_Failure;
3111 : }
3112 :
3113 133614 : return VRTRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
3114 : }
3115 :
3116 : /************************************************************************/
3117 : /* SetMetadata() */
3118 : /************************************************************************/
3119 :
3120 70872 : CPLErr VRTSourcedRasterBand::SetMetadata(CSLConstList papszNewMD,
3121 : const char *pszDomain)
3122 :
3123 : {
3124 70872 : if (pszDomain != nullptr && (EQUAL(pszDomain, "new_vrt_sources") ||
3125 70872 : EQUAL(pszDomain, "vrt_sources")))
3126 : {
3127 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
3128 8 : GetGDALDriverManager()->GetDriverByName("VRT"));
3129 8 : if (!poDriver)
3130 0 : return CE_Failure;
3131 :
3132 8 : if (EQUAL(pszDomain, "vrt_sources"))
3133 : {
3134 8 : m_papoSources.clear();
3135 : }
3136 :
3137 6 : for (const char *const pszMDItem :
3138 20 : cpl::Iterate(CSLConstList(papszNewMD)))
3139 : {
3140 8 : const char *const pszXML = CPLParseNameValue(pszMDItem, nullptr);
3141 8 : CPLXMLTreeCloser psTree(CPLParseXMLString(pszXML));
3142 8 : if (!psTree)
3143 0 : return CE_Failure;
3144 :
3145 8 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
3146 8 : if (l_poDS == nullptr)
3147 : {
3148 0 : return CE_Failure;
3149 : }
3150 8 : VRTSource *const poSource = poDriver->ParseSource(
3151 8 : psTree.get(), nullptr, l_poDS->m_oMapSharedSources);
3152 8 : if (poSource == nullptr)
3153 2 : return CE_Failure;
3154 :
3155 6 : const CPLErr eErr = AddSource(poSource);
3156 : // cppcheck-suppress knownConditionTrueFalse
3157 6 : if (eErr != CE_None)
3158 0 : return eErr;
3159 : }
3160 :
3161 6 : return CE_None;
3162 : }
3163 :
3164 70864 : return VRTRasterBand::SetMetadata(papszNewMD, pszDomain);
3165 : }
3166 :
3167 : /************************************************************************/
3168 : /* GetFileList() */
3169 : /************************************************************************/
3170 :
3171 68 : void VRTSourcedRasterBand::GetFileList(char ***ppapszFileList, int *pnSize,
3172 : int *pnMaxSize, CPLHashSet *hSetFiles)
3173 : {
3174 179 : for (auto &poSource : m_papoSources)
3175 : {
3176 111 : poSource->GetFileList(ppapszFileList, pnSize, pnMaxSize, hSetFiles);
3177 : }
3178 :
3179 68 : VRTRasterBand::GetFileList(ppapszFileList, pnSize, pnMaxSize, hSetFiles);
3180 68 : }
3181 :
3182 : /************************************************************************/
3183 : /* CloseDependentDatasets() */
3184 : /************************************************************************/
3185 :
3186 142284 : int VRTSourcedRasterBand::CloseDependentDatasets()
3187 : {
3188 142284 : int ret = VRTRasterBand::CloseDependentDatasets();
3189 :
3190 142284 : if (m_papoSources.empty())
3191 395 : return ret;
3192 :
3193 141889 : m_papoSources.clear();
3194 :
3195 141889 : return TRUE;
3196 : }
3197 :
3198 : /************************************************************************/
3199 : /* FlushCache() */
3200 : /************************************************************************/
3201 :
3202 208450 : CPLErr VRTSourcedRasterBand::FlushCache(bool bAtClosing)
3203 : {
3204 208450 : CPLErr eErr = VRTRasterBand::FlushCache(bAtClosing);
3205 522330 : for (size_t i = 0; i < m_papoSources.size() && eErr == CE_None; i++)
3206 : {
3207 313880 : eErr = m_papoSources[i]->FlushCache(bAtClosing);
3208 : }
3209 208450 : return eErr;
3210 : }
3211 :
3212 : /************************************************************************/
3213 : /* RemoveCoveredSources() */
3214 : /************************************************************************/
3215 :
3216 : /** Remove sources that are covered by other sources.
3217 : *
3218 : * This method removes sources that are covered entirely by (one or several)
3219 : * sources of higher priority (even if they declare a nodata setting).
3220 : * This optimizes the size of the VRT and the rendering time.
3221 : */
3222 14 : void VRTSourcedRasterBand::RemoveCoveredSources(CSLConstList papszOptions)
3223 : {
3224 : #ifndef HAVE_GEOS
3225 : if (CPLTestBool(CSLFetchNameValueDef(
3226 : papszOptions, "EMIT_ERROR_IF_GEOS_NOT_AVAILABLE", "TRUE")))
3227 : {
3228 : CPLError(CE_Failure, CPLE_NotSupported,
3229 : "RemoveCoveredSources() not implemented in builds "
3230 : "without GEOS support");
3231 : }
3232 : #else
3233 : (void)papszOptions;
3234 :
3235 : CPLRectObj globalBounds;
3236 14 : globalBounds.minx = 0;
3237 14 : globalBounds.miny = 0;
3238 14 : globalBounds.maxx = nRasterXSize;
3239 14 : globalBounds.maxy = nRasterYSize;
3240 :
3241 : // Create an index with the bbox of all sources
3242 14 : CPLQuadTree *hTree = CPLQuadTreeCreate(&globalBounds, nullptr);
3243 38 : for (int i = 0; i < static_cast<int>(m_papoSources.size()); i++)
3244 : {
3245 24 : if (m_papoSources[i]->IsSimpleSource())
3246 : {
3247 : VRTSimpleSource *poSS =
3248 24 : cpl::down_cast<VRTSimpleSource *>(m_papoSources[i].get());
3249 24 : void *hFeature =
3250 24 : reinterpret_cast<void *>(static_cast<uintptr_t>(i));
3251 : CPLRectObj rect;
3252 24 : rect.minx = std::max(0.0, poSS->m_dfDstXOff);
3253 24 : rect.miny = std::max(0.0, poSS->m_dfDstYOff);
3254 48 : rect.maxx = std::min(double(nRasterXSize),
3255 24 : poSS->m_dfDstXOff + poSS->m_dfDstXSize);
3256 48 : rect.maxy = std::min(double(nRasterYSize),
3257 24 : poSS->m_dfDstYOff + poSS->m_dfDstYSize);
3258 24 : CPLQuadTreeInsertWithBounds(hTree, hFeature, &rect);
3259 : }
3260 : }
3261 :
3262 38 : for (int i = 0; i < static_cast<int>(m_papoSources.size()); i++)
3263 : {
3264 24 : if (m_papoSources[i]->IsSimpleSource())
3265 : {
3266 : VRTSimpleSource *poSS =
3267 24 : cpl::down_cast<VRTSimpleSource *>(m_papoSources[i].get());
3268 : CPLRectObj rect;
3269 24 : rect.minx = std::max(0.0, poSS->m_dfDstXOff);
3270 24 : rect.miny = std::max(0.0, poSS->m_dfDstYOff);
3271 48 : rect.maxx = std::min(double(nRasterXSize),
3272 24 : poSS->m_dfDstXOff + poSS->m_dfDstXSize);
3273 48 : rect.maxy = std::min(double(nRasterYSize),
3274 24 : poSS->m_dfDstYOff + poSS->m_dfDstYSize);
3275 :
3276 : // Find sources whose extent intersect with the current one
3277 24 : int nFeatureCount = 0;
3278 : void **pahFeatures =
3279 24 : CPLQuadTreeSearch(hTree, &rect, &nFeatureCount);
3280 :
3281 : // Compute the bounding box of those sources, only if they are
3282 : // on top of the current one
3283 : CPLRectObj rectIntersecting;
3284 24 : rectIntersecting.minx = std::numeric_limits<double>::max();
3285 24 : rectIntersecting.miny = std::numeric_limits<double>::max();
3286 24 : rectIntersecting.maxx = -std::numeric_limits<double>::max();
3287 24 : rectIntersecting.maxy = -std::numeric_limits<double>::max();
3288 61 : for (int j = 0; j < nFeatureCount; j++)
3289 : {
3290 37 : const int curFeature = static_cast<int>(
3291 37 : reinterpret_cast<uintptr_t>(pahFeatures[j]));
3292 37 : if (curFeature > i)
3293 : {
3294 : VRTSimpleSource *poOtherSS =
3295 13 : cpl::down_cast<VRTSimpleSource *>(
3296 13 : m_papoSources[curFeature].get());
3297 13 : rectIntersecting.minx =
3298 13 : std::min(rectIntersecting.minx, poOtherSS->m_dfDstXOff);
3299 13 : rectIntersecting.miny =
3300 13 : std::min(rectIntersecting.miny, poOtherSS->m_dfDstYOff);
3301 13 : rectIntersecting.maxx = std::max(
3302 : rectIntersecting.maxx,
3303 13 : poOtherSS->m_dfDstXOff + poOtherSS->m_dfDstXSize);
3304 13 : rectIntersecting.maxy = std::max(
3305 : rectIntersecting.maxy,
3306 13 : poOtherSS->m_dfDstYOff + poOtherSS->m_dfDstXSize);
3307 : }
3308 : }
3309 :
3310 : // If the boundinx box of those sources overlap the current one,
3311 : // then compute their union, and check if it contains the current
3312 : // source
3313 24 : if (rectIntersecting.minx <= rect.minx &&
3314 7 : rectIntersecting.miny <= rect.miny &&
3315 7 : rectIntersecting.maxx >= rect.maxx &&
3316 7 : rectIntersecting.maxy >= rect.maxy)
3317 : {
3318 14 : OGRPolygon oPoly;
3319 : {
3320 7 : auto poLR = new OGRLinearRing();
3321 7 : poLR->addPoint(rect.minx, rect.miny);
3322 7 : poLR->addPoint(rect.minx, rect.maxy);
3323 7 : poLR->addPoint(rect.maxx, rect.maxy);
3324 7 : poLR->addPoint(rect.maxx, rect.miny);
3325 7 : poLR->addPoint(rect.minx, rect.miny);
3326 7 : oPoly.addRingDirectly(poLR);
3327 : }
3328 :
3329 7 : std::unique_ptr<OGRGeometry> poUnion;
3330 24 : for (int j = 0; j < nFeatureCount; j++)
3331 : {
3332 17 : const int curFeature = static_cast<int>(
3333 17 : reinterpret_cast<uintptr_t>(pahFeatures[j]));
3334 17 : if (curFeature > i)
3335 : {
3336 : VRTSimpleSource *poOtherSS =
3337 10 : cpl::down_cast<VRTSimpleSource *>(
3338 10 : m_papoSources[curFeature].get());
3339 : CPLRectObj otherRect;
3340 10 : otherRect.minx = std::max(0.0, poOtherSS->m_dfDstXOff);
3341 10 : otherRect.miny = std::max(0.0, poOtherSS->m_dfDstYOff);
3342 20 : otherRect.maxx = std::min(double(nRasterXSize),
3343 20 : poOtherSS->m_dfDstXOff +
3344 10 : poOtherSS->m_dfDstXSize);
3345 20 : otherRect.maxy = std::min(double(nRasterYSize),
3346 20 : poOtherSS->m_dfDstYOff +
3347 10 : poOtherSS->m_dfDstYSize);
3348 20 : OGRPolygon oOtherPoly;
3349 : {
3350 10 : auto poLR = new OGRLinearRing();
3351 10 : poLR->addPoint(otherRect.minx, otherRect.miny);
3352 10 : poLR->addPoint(otherRect.minx, otherRect.maxy);
3353 10 : poLR->addPoint(otherRect.maxx, otherRect.maxy);
3354 10 : poLR->addPoint(otherRect.maxx, otherRect.miny);
3355 10 : poLR->addPoint(otherRect.minx, otherRect.miny);
3356 10 : oOtherPoly.addRingDirectly(poLR);
3357 : }
3358 10 : if (poUnion == nullptr)
3359 7 : poUnion.reset(oOtherPoly.clone());
3360 : else
3361 3 : poUnion.reset(oOtherPoly.Union(poUnion.get()));
3362 : }
3363 : }
3364 :
3365 7 : if (poUnion != nullptr && poUnion->Contains(&oPoly))
3366 : {
3367 : // We can remove the current source
3368 7 : m_papoSources[i].reset();
3369 : }
3370 : }
3371 24 : CPLFree(pahFeatures);
3372 :
3373 24 : void *hFeature =
3374 24 : reinterpret_cast<void *>(static_cast<uintptr_t>(i));
3375 24 : CPLQuadTreeRemove(hTree, hFeature, &rect);
3376 : }
3377 : }
3378 :
3379 : // Compact the papoSources array
3380 0 : m_papoSources.erase(std::remove_if(m_papoSources.begin(),
3381 : m_papoSources.end(),
3382 24 : [](const std::unique_ptr<VRTSource> &src)
3383 38 : { return src.get() == nullptr; }),
3384 28 : m_papoSources.end());
3385 :
3386 14 : CPLQuadTreeDestroy(hTree);
3387 : #endif
3388 14 : }
3389 :
3390 : /*! @endcond */
|
