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 2584 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn)
59 : {
60 2584 : VRTRasterBand::Initialize(poDSIn->GetRasterXSize(),
61 : poDSIn->GetRasterYSize());
62 :
63 2584 : poDS = poDSIn;
64 2584 : nBand = nBandIn;
65 2584 : }
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 139528 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn,
95 : GDALDataType eType, int nXSize,
96 : int nYSize, int nBlockXSizeIn,
97 139528 : int nBlockYSizeIn)
98 : {
99 139528 : VRTRasterBand::Initialize(nXSize, nYSize);
100 :
101 139528 : poDS = poDSIn;
102 139528 : nBand = nBandIn;
103 :
104 139528 : eDataType = eType;
105 139528 : if (nBlockXSizeIn > 0)
106 139150 : nBlockXSize = nBlockXSizeIn;
107 139528 : if (nBlockYSizeIn > 0)
108 139168 : nBlockYSize = nBlockYSizeIn;
109 139528 : }
110 :
111 : /************************************************************************/
112 : /* ~VRTSourcedRasterBand() */
113 : /************************************************************************/
114 :
115 281954 : VRTSourcedRasterBand::~VRTSourcedRasterBand()
116 :
117 : {
118 142112 : VRTSourcedRasterBand::CloseDependentDatasets();
119 281954 : }
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 279489 : 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 279489 : };
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 279489 : if (eRWFlag == GF_Read && (nXSize != nBufXSize || nYSize != nBufYSize) &&
185 558978 : !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 279475 : 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 51367 : void VRTSourcedRasterBand::InitializeOutputBuffer(void *pData, int nBufXSize,
516 : int nBufYSize,
517 : GDALDataType eBufType,
518 : GSpacing nPixelSpace,
519 : GSpacing nLineSpace) const
520 : {
521 51367 : if (nPixelSpace == GDALGetDataTypeSizeBytes(eBufType) &&
522 29844 : !(m_bNoDataValueSet && m_dfNoDataValue != 0.0) &&
523 110794 : !(m_bNoDataSetAsInt64 && m_nNoDataValueInt64 != 0) &&
524 29583 : !(m_bNoDataSetAsUInt64 && m_nNoDataValueUInt64 != 0))
525 : {
526 29582 : if (nLineSpace == nBufXSize * nPixelSpace)
527 : {
528 28502 : 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 21785 : 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 21784 : 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 21783 : double dfWriteValue = 0.0;
563 21783 : if (m_bNoDataValueSet)
564 287 : dfWriteValue = m_dfNoDataValue;
565 :
566 111188 : for (int iLine = 0; iLine < nBufYSize; iLine++)
567 : {
568 89405 : GDALCopyWords(&dfWriteValue, GDT_Float64, 0,
569 : static_cast<GByte *>(pData) +
570 89405 : static_cast<GIntBig>(nLineSpace) * iLine,
571 : eBufType, static_cast<int>(nPixelSpace), nBufXSize);
572 : }
573 : }
574 51367 : }
575 :
576 : /************************************************************************/
577 : /* IRasterIO() */
578 : /************************************************************************/
579 :
580 279490 : 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 279490 : 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 558980 : const std::string osFctId("VRTSourcedRasterBand::IRasterIO");
594 558980 : GDALAntiRecursionGuard oGuard(osFctId);
595 279490 : if (oGuard.GetCallDepth() >= 32)
596 : {
597 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
598 0 : return CE_Failure;
599 : }
600 :
601 838470 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
602 : // Allow 2 recursion depths on the same dataset for non-nearest resampling
603 279490 : 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 279489 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
614 279489 : if (l_poDS &&
615 558627 : l_poDS->m_apoOverviews.empty() && // do not use virtual overviews
616 558978 : (nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() > 0)
617 : {
618 2 : if (OverviewRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
619 : nBufXSize, nBufYSize, eBufType, nPixelSpace,
620 2 : nLineSpace, psExtraArg) == CE_None)
621 1 : 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 279488 : 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 279474 : if (!SkipBufferInitialization())
684 : {
685 48762 : InitializeOutputBuffer(pData, nBufXSize, nBufYSize, eBufType,
686 : nPixelSpace, nLineSpace);
687 : }
688 :
689 : /* -------------------------------------------------------------------- */
690 : /* Overlay each source in turn over top this. */
691 : /* -------------------------------------------------------------------- */
692 279474 : CPLErr eErr = CE_None;
693 :
694 279474 : double dfXOff = nXOff;
695 279474 : double dfYOff = nYOff;
696 279474 : double dfXSize = nXSize;
697 279474 : double dfYSize = nYSize;
698 279474 : if (psExtraArg->bFloatingPointWindowValidity)
699 : {
700 297 : dfXOff = psExtraArg->dfXOff;
701 297 : dfYOff = psExtraArg->dfYOff;
702 297 : dfXSize = psExtraArg->dfXSize;
703 297 : dfYSize = psExtraArg->dfYSize;
704 : }
705 :
706 279474 : if (l_poDS)
707 279474 : l_poDS->m_bMultiThreadedRasterIOLastUsed = false;
708 :
709 279474 : int nContributingSources = 0;
710 279474 : int nMaxThreads = 0;
711 279474 : constexpr int MINIMUM_PIXEL_COUNT_FOR_THREADED_IO = 1000 * 1000;
712 558948 : if (l_poDS &&
713 279474 : (static_cast<int64_t>(nBufXSize) * nBufYSize >=
714 279141 : MINIMUM_PIXEL_COUNT_FOR_THREADED_IO ||
715 279141 : 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 558968 : 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 77 : while (oQueue->WaitEvent())
793 : {
794 : // Quite rough progress callback. We could do better by counting
795 : // the number of contributing pixels.
796 59 : if (psExtraArg->pfnProgress)
797 : {
798 100 : 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 279456 : GDALProgressFunc const pfnProgressGlobal = psExtraArg->pfnProgress;
810 279456 : void *const pProgressDataGlobal = psExtraArg->pProgressData;
811 :
812 279456 : VRTSource::WorkingState oWorkingState;
813 279456 : const int nSources = static_cast<int>(m_papoSources.size());
814 554028 : for (int iSource = 0; eErr == CE_None && iSource < nSources; iSource++)
815 : {
816 274572 : psExtraArg->pfnProgress = GDALScaledProgress;
817 549144 : psExtraArg->pProgressData = GDALCreateScaledProgress(
818 274572 : 1.0 * iSource / nSources, 1.0 * (iSource + 1) / nSources,
819 : pfnProgressGlobal, pProgressDataGlobal);
820 274572 : if (psExtraArg->pProgressData == nullptr)
821 274022 : psExtraArg->pfnProgress = nullptr;
822 :
823 549144 : eErr = m_papoSources[iSource]->RasterIO(
824 : eDataType, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
825 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg,
826 274572 : l_poDS ? l_poDS->m_oWorkingState : oWorkingState);
827 :
828 274572 : GDALDestroyScaledProgress(psExtraArg->pProgressData);
829 : }
830 :
831 279456 : psExtraArg->pfnProgress = pfnProgressGlobal;
832 279456 : psExtraArg->pProgressData = pProgressDataGlobal;
833 : }
834 :
835 279474 : if (eErr == CE_None && psExtraArg->pfnProgress)
836 : {
837 417 : psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
838 : }
839 :
840 279474 : return eErr;
841 : }
842 :
843 : /************************************************************************/
844 : /* IGetDataCoverageStatus() */
845 : /************************************************************************/
846 :
847 3031 : int VRTSourcedRasterBand::IGetDataCoverageStatus(int nXOff, int nYOff,
848 : int nXSize, int nYSize,
849 : int nMaskFlagStop,
850 : double *pdfDataPct)
851 : {
852 3031 : if (pdfDataPct)
853 8 : *pdfDataPct = -1.0;
854 :
855 : // Particular case for a single simple source covering the whole dataset
856 6047 : if (m_papoSources.size() == 1 && m_papoSources[0]->IsSimpleSource() &&
857 3016 : m_papoSources[0]->GetType() == VRTSimpleSource::GetTypeStatic())
858 : {
859 : VRTSimpleSource *poSource =
860 2977 : static_cast<VRTSimpleSource *>(m_papoSources[0].get());
861 :
862 2977 : GDALRasterBand *poBand = poSource->GetRasterBand();
863 2977 : if (!poBand)
864 0 : poBand = poSource->GetMaskBandMainBand();
865 2977 : if (!poBand)
866 : {
867 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
868 352 : GDAL_DATA_COVERAGE_STATUS_DATA;
869 : }
870 :
871 : /* Check that it uses the full source dataset */
872 2977 : double dfReqXOff = 0.0;
873 2977 : double dfReqYOff = 0.0;
874 2977 : double dfReqXSize = 0.0;
875 2977 : double dfReqYSize = 0.0;
876 2977 : int nReqXOff = 0;
877 2977 : int nReqYOff = 0;
878 2977 : int nReqXSize = 0;
879 2977 : int nReqYSize = 0;
880 2977 : int nOutXOff = 0;
881 2977 : int nOutYOff = 0;
882 2977 : int nOutXSize = 0;
883 2977 : int nOutYSize = 0;
884 2977 : bool bError = false;
885 2977 : if (poSource->GetSrcDstWindow(
886 2977 : 0, 0, GetXSize(), GetYSize(), GetXSize(), GetYSize(),
887 : GRIORA_NearestNeighbour, &dfReqXOff, &dfReqYOff, &dfReqXSize,
888 : &dfReqYSize, &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
889 2977 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError) &&
890 2977 : nReqXOff == 0 && nReqYOff == 0 && nReqXSize == GetXSize() &&
891 467 : nReqXSize == poBand->GetXSize() && nReqYSize == GetYSize() &&
892 361 : nReqYSize == poBand->GetYSize() && nOutXOff == 0 && nOutYOff == 0 &&
893 5954 : nOutXSize == GetXSize() && nOutYSize == GetYSize())
894 : {
895 352 : return poBand->GetDataCoverageStatus(nXOff, nYOff, nXSize, nYSize,
896 352 : 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 2679 : int nStatus = 0;
905 :
906 5358 : auto poPolyNonCoveredBySources = std::make_unique<OGRPolygon>();
907 : {
908 5358 : auto poLR = std::make_unique<OGRLinearRing>();
909 2679 : poLR->addPoint(nXOff, nYOff);
910 2679 : poLR->addPoint(nXOff, nYOff + nYSize);
911 2679 : poLR->addPoint(nXOff + nXSize, nYOff + nYSize);
912 2679 : poLR->addPoint(nXOff + nXSize, nYOff);
913 2679 : poLR->addPoint(nXOff, nYOff);
914 2679 : poPolyNonCoveredBySources->addRingDirectly(poLR.release());
915 : }
916 :
917 2689 : for (auto &poSource : m_papoSources)
918 : {
919 2684 : if (!poSource->IsSimpleSource())
920 : {
921 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
922 2674 : GDAL_DATA_COVERAGE_STATUS_DATA;
923 : }
924 2680 : VRTSimpleSource *poSS = static_cast<VRTSimpleSource *>(poSource.get());
925 : // Check if the AOI is fully inside the source
926 2680 : double dfDstXOff = std::max(0.0, poSS->m_dfDstXOff);
927 2680 : double dfDstYOff = std::max(0.0, poSS->m_dfDstYOff);
928 2680 : double dfDstXSize = poSS->m_dfDstXSize;
929 2680 : double dfDstYSize = poSS->m_dfDstYSize;
930 2680 : auto l_poBand = poSS->GetRasterBand();
931 2680 : if (!l_poBand)
932 0 : continue;
933 2680 : if (dfDstXSize == -1)
934 8 : dfDstXSize = l_poBand->GetXSize() - dfDstXOff;
935 2680 : if (dfDstYSize == -1)
936 8 : dfDstYSize = l_poBand->GetYSize() - dfDstYOff;
937 :
938 2680 : if (nXOff >= dfDstXOff && nYOff >= dfDstYOff &&
939 2671 : nXOff + nXSize <= dfDstXOff + dfDstXSize &&
940 2647 : nYOff + nYSize <= dfDstYOff + dfDstYSize)
941 : {
942 2647 : if (pdfDataPct)
943 1 : *pdfDataPct = 100.0;
944 2647 : 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 54 : bool VRTSourcedRasterBand::
1288 : IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1289 : bool bAllowMaxValAdjustment) const
1290 : {
1291 54 : bool bRet = true;
1292 : CPLRectObj sGlobalBounds;
1293 54 : sGlobalBounds.minx = 0;
1294 54 : sGlobalBounds.miny = 0;
1295 54 : sGlobalBounds.maxx = nRasterXSize;
1296 54 : sGlobalBounds.maxy = nRasterYSize;
1297 54 : CPLQuadTree *hQuadTree = CPLQuadTreeCreate(&sGlobalBounds, nullptr);
1298 124 : for (int i = 0; i < static_cast<int>(m_papoSources.size()); ++i)
1299 : {
1300 79 : if (!m_papoSources[i]->IsSimpleSource())
1301 : {
1302 1 : bRet = false;
1303 9 : break;
1304 : }
1305 :
1306 : auto poSimpleSource =
1307 78 : cpl::down_cast<VRTSimpleSource *>(m_papoSources[i].get());
1308 78 : const char *pszType = poSimpleSource->GetType();
1309 78 : 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 76 : if (!bAllowMaxValAdjustment && poSimpleSource->NeedMaxValAdjustment())
1330 : {
1331 2 : bRet = false;
1332 2 : break;
1333 : }
1334 :
1335 74 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1336 148 : if (poSimpleSourceBand == nullptr ||
1337 74 : poSimpleSourceBand->GetRasterDataType() != eDataType)
1338 : {
1339 0 : bRet = false;
1340 0 : break;
1341 : }
1342 :
1343 74 : double dfReqXOff = 0.0;
1344 74 : double dfReqYOff = 0.0;
1345 74 : double dfReqXSize = 0.0;
1346 74 : double dfReqYSize = 0.0;
1347 74 : int nReqXOff = 0;
1348 74 : int nReqYOff = 0;
1349 74 : int nReqXSize = 0;
1350 74 : int nReqYSize = 0;
1351 74 : int nOutXOff = 0;
1352 74 : int nOutYOff = 0;
1353 74 : int nOutXSize = 0;
1354 74 : int nOutYSize = 0;
1355 :
1356 74 : bool bError = false;
1357 222 : if (!poSimpleSource->GetSrcDstWindow(
1358 74 : 0, 0, nRasterXSize, nRasterYSize, nRasterXSize, nRasterYSize,
1359 : GRIORA_NearestNeighbour, &dfReqXOff, &dfReqYOff, &dfReqXSize,
1360 : &dfReqYSize, &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
1361 74 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError) ||
1362 74 : nReqXOff != 0 || nReqYOff != 0 ||
1363 73 : nReqXSize != poSimpleSourceBand->GetXSize() ||
1364 73 : nReqYSize != poSimpleSourceBand->GetYSize() ||
1365 148 : nOutXSize != nReqXSize || nOutYSize != nReqYSize)
1366 : {
1367 1 : bRet = false;
1368 1 : break;
1369 : }
1370 :
1371 : CPLRectObj sBounds;
1372 73 : constexpr double EPSILON = 1e-1;
1373 73 : sBounds.minx = nOutXOff + EPSILON;
1374 73 : sBounds.miny = nOutYOff + EPSILON;
1375 73 : sBounds.maxx = nOutXOff + nOutXSize - EPSILON;
1376 73 : sBounds.maxy = nOutYOff + nOutYSize - EPSILON;
1377 :
1378 : // Check that the new source doesn't overlap an existing one.
1379 73 : if (CPLQuadTreeHasMatch(hQuadTree, &sBounds))
1380 : {
1381 3 : bRet = false;
1382 3 : break;
1383 : }
1384 :
1385 70 : CPLQuadTreeInsertWithBounds(
1386 70 : hQuadTree, reinterpret_cast<void *>(static_cast<uintptr_t>(i)),
1387 : &sBounds);
1388 : }
1389 54 : CPLQuadTreeDestroy(hQuadTree);
1390 :
1391 54 : return bRet;
1392 : }
1393 :
1394 : /************************************************************************/
1395 : /* ComputeRasterMinMax() */
1396 : /************************************************************************/
1397 :
1398 22 : CPLErr VRTSourcedRasterBand::ComputeRasterMinMax(int bApproxOK,
1399 : double *adfMinMax)
1400 : {
1401 : /* -------------------------------------------------------------------- */
1402 : /* Does the driver already know the min/max? */
1403 : /* -------------------------------------------------------------------- */
1404 22 : 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 42 : const std::string osFctId("VRTSourcedRasterBand::ComputeRasterMinMax");
1421 42 : GDALAntiRecursionGuard oGuard(osFctId);
1422 21 : if (oGuard.GetCallDepth() >= 32)
1423 : {
1424 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1425 0 : return CE_Failure;
1426 : }
1427 :
1428 63 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1429 21 : 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 21 : 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 21 : if (IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1464 : /*bAllowMaxValAdjustment = */ true))
1465 : {
1466 18 : CPLDebugOnly(
1467 : "VRT", "ComputeRasterMinMax(): use optimized code path for mosaic");
1468 :
1469 18 : 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 18 : bool bUseComputeStatistics = false;
1476 37 : for (auto &poSource : m_papoSources)
1477 : {
1478 : auto poSimpleSource =
1479 23 : cpl::down_cast<VRTSimpleSource *>(poSource.get());
1480 23 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1481 : // already checked by IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange
1482 23 : assert(poSimpleSourceBand);
1483 23 : int bHasNoData = FALSE;
1484 23 : CPL_IGNORE_RET_VAL(poSimpleSourceBand->GetNoDataValue(&bHasNoData));
1485 23 : if (bHasNoData)
1486 : {
1487 4 : bUseComputeStatistics = true;
1488 4 : break;
1489 : }
1490 19 : nCoveredArea +=
1491 19 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
1492 19 : poSimpleSourceBand->GetYSize();
1493 : }
1494 :
1495 18 : if (bUseComputeStatistics)
1496 : {
1497 : CPLErr eErr;
1498 4 : std::string osLastErrorMsg;
1499 : {
1500 8 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1501 4 : CPLErrorReset();
1502 : eErr =
1503 8 : ComputeStatistics(bApproxOK, &adfMinMax[0], &adfMinMax[1],
1504 4 : nullptr, nullptr, nullptr, nullptr);
1505 4 : if (eErr == CE_Failure)
1506 : {
1507 1 : osLastErrorMsg = CPLGetLastErrorMsg();
1508 : }
1509 : }
1510 4 : if (eErr == CE_Failure)
1511 : {
1512 1 : if (strstr(osLastErrorMsg.c_str(), "no valid pixels found") !=
1513 : nullptr)
1514 : {
1515 1 : 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 4 : 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", "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 35 : CPLErr VRTSourcedRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin,
1708 : double *pdfMax, double *pdfMean,
1709 : double *pdfStdDev,
1710 : GDALProgressFunc pfnProgress,
1711 : void *pProgressData)
1712 :
1713 : {
1714 70 : const std::string osFctId("VRTSourcedRasterBand::ComputeStatistics");
1715 70 : GDALAntiRecursionGuard oGuard(osFctId);
1716 35 : if (oGuard.GetCallDepth() >= 32)
1717 : {
1718 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1719 0 : return CE_Failure;
1720 : }
1721 :
1722 105 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1723 35 : 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 35 : 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 33 : if (IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1770 : /*bAllowMaxValAdjustment = */ false))
1771 : {
1772 27 : Context sContext;
1773 27 : sContext.bApproxOK = bApproxOK;
1774 27 : sContext.dfNoDataValue = m_dfNoDataValue;
1775 27 : sContext.bNoDataValueSet = m_bNoDataValueSet;
1776 27 : sContext.bHideNoDataValue = m_bHideNoDataValue;
1777 27 : sContext.pfnProgress = pfnProgress;
1778 27 : sContext.pProgressData = pProgressData;
1779 27 : 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 37 : static void UpdateStats(const Job *psJob)
1834 : {
1835 37 : const auto nValidPixels = psJob->nValidPixels;
1836 37 : auto psContext = psJob->psContext;
1837 37 : 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 37 : int bHasNoData = FALSE;
1865 : const double dfNoDataValue =
1866 37 : psJob->poRasterBand->GetNoDataValue(&bHasNoData);
1867 9 : if (nValidPixels < psJob->nPixelCount && bHasNoData &&
1868 55 : !std::isnan(dfNoDataValue) &&
1869 9 : (!psContext->bNoDataValueSet ||
1870 7 : 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 37 : 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 37 : }
1901 : };
1902 :
1903 40 : const auto JobRunner = [](void *pData)
1904 : {
1905 40 : auto psJob = static_cast<Job *>(pData);
1906 40 : auto psContext = psJob->psContext;
1907 : {
1908 40 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1909 40 : if (psContext->bFallbackToBase || psContext->bFailure)
1910 0 : return;
1911 : }
1912 :
1913 40 : auto poSimpleSourceBand = psJob->poRasterBand;
1914 40 : psJob->nPixelCount =
1915 40 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
1916 40 : poSimpleSourceBand->GetYSize();
1917 :
1918 40 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1919 40 : CPLErr eErr = poSimpleSourceBand->ComputeStatistics(
1920 40 : psContext->bApproxOK, &psJob->dfMin, &psJob->dfMax,
1921 : &psJob->dfMean, &psJob->dfStdDev,
1922 46 : psContext->pfnProgress == nullptr ||
1923 6 : psContext->pfnProgress == GDALDummyProgress
1924 : ? GDALDummyProgress
1925 : : Job::ProgressFunc,
1926 40 : psJob);
1927 : const char *pszValidPercent =
1928 40 : poSimpleSourceBand->GetMetadataItem("STATISTICS_VALID_PERCENT");
1929 40 : psJob->nValidPixels =
1930 : pszValidPercent
1931 40 : ? static_cast<uint64_t>(CPLAtof(pszValidPercent) *
1932 40 : psJob->nPixelCount / 100.0)
1933 : : psJob->nPixelCount;
1934 40 : if (eErr == CE_Failure)
1935 : {
1936 16 : if (pszValidPercent != nullptr &&
1937 8 : 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 27 : CPLWorkerThreadPool *poThreadPool = nullptr;
1965 : int nThreads =
1966 27 : m_papoSources.size() > 1
1967 27 : ? VRTDataset::GetNumThreads(dynamic_cast<VRTDataset *>(poDS))
1968 27 : : 0;
1969 27 : if (nThreads > 1024)
1970 0 : nThreads = 1024; // to please Coverity
1971 27 : 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 24 : std::set<std::string> oSetDatasetNames;
1978 24 : std::set<GDALDataset *> oSetDatasetPointers;
1979 36 : for (auto &poSource : m_papoSources)
1980 : {
1981 : auto poSimpleSource =
1982 24 : static_cast<VRTSimpleSource *>(poSource.get());
1983 24 : assert(poSimpleSource);
1984 24 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1985 24 : assert(poSimpleSourceBand);
1986 24 : auto poSourceDataset = poSimpleSourceBand->GetDataset();
1987 24 : if (poSourceDataset == nullptr)
1988 : {
1989 0 : nThreads = 0;
1990 0 : break;
1991 : }
1992 24 : auto poDriver = poSourceDataset->GetDriver();
1993 24 : if (poDriver && EQUAL(poDriver->GetDescription(), "MEM"))
1994 : {
1995 24 : if (oSetDatasetPointers.find(poSourceDataset) !=
1996 48 : oSetDatasetPointers.end())
1997 : {
1998 0 : nThreads = 0;
1999 0 : break;
2000 : }
2001 24 : 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 12 : if (nThreads > 1)
2016 : {
2017 12 : poThreadPool = GDALGetGlobalThreadPool(nThreads);
2018 : }
2019 : }
2020 :
2021 : // Compute total number of pixels of sources
2022 67 : for (auto &poSource : m_papoSources)
2023 : {
2024 : auto poSimpleSource =
2025 40 : static_cast<VRTSimpleSource *>(poSource.get());
2026 40 : assert(poSimpleSource);
2027 40 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
2028 40 : assert(poSimpleSourceBand);
2029 40 : sContext.nTotalPixelsOfSources +=
2030 40 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
2031 40 : poSimpleSourceBand->GetYSize();
2032 : }
2033 :
2034 27 : if (poThreadPool)
2035 : {
2036 12 : CPLDebugOnly("VRT", "ComputeStatistics(): use optimized "
2037 : "multi-threaded code path for mosaic");
2038 24 : std::vector<Job> asJobs(m_papoSources.size());
2039 24 : auto poQueue = poThreadPool->CreateJobQueue();
2040 36 : for (size_t i = 0; i < m_papoSources.size(); ++i)
2041 : {
2042 : auto poSimpleSource =
2043 24 : static_cast<VRTSimpleSource *>(m_papoSources[i].get());
2044 24 : assert(poSimpleSource);
2045 24 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
2046 24 : assert(poSimpleSourceBand);
2047 24 : asJobs[i].psContext = &sContext;
2048 24 : asJobs[i].poRasterBand = poSimpleSourceBand;
2049 24 : if (!poQueue->SubmitJob(JobRunner, &asJobs[i]))
2050 : {
2051 0 : sContext.bFailure = true;
2052 0 : break;
2053 : }
2054 : }
2055 12 : poQueue->WaitCompletion();
2056 12 : if (!(sContext.bFailure || sContext.bFallbackToBase))
2057 : {
2058 33 : for (size_t i = 0; i < m_papoSources.size(); ++i)
2059 : {
2060 22 : 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 27 : if (sContext.bFailure)
2087 0 : return CE_Failure;
2088 27 : 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 25 : const uint64_t nTotalPixels =
2102 25 : static_cast<uint64_t>(nRasterXSize) * nRasterYSize;
2103 7 : if (m_bNoDataValueSet && m_bHideNoDataValue &&
2104 32 : !std::isnan(m_dfNoDataValue) && IsNoDataValueInDataTypeRange())
2105 : {
2106 1 : UpdateStatsWithConstantValue(sContext, m_dfNoDataValue,
2107 : nTotalPixels -
2108 1 : sContext.nGlobalValidPixels);
2109 : }
2110 24 : 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 25 : double dfGlobalMean = sContext.dfGlobalMean;
2127 : double dfGlobalStdDev =
2128 25 : sContext.nGlobalValidPixels > 0
2129 25 : ? sqrt(sContext.dfGlobalM2 / sContext.nGlobalValidPixels)
2130 25 : : 0;
2131 : #endif
2132 38 : 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 25 : 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 2 : sContext.dfGlobalMin = 0.0;
2157 2 : sContext.dfGlobalMax = 0.0;
2158 : }
2159 :
2160 25 : SetValidPercent(nTotalPixels, sContext.nGlobalValidPixels);
2161 :
2162 25 : if (pdfMin)
2163 17 : *pdfMin = sContext.dfGlobalMin;
2164 25 : if (pdfMax)
2165 17 : *pdfMax = sContext.dfGlobalMax;
2166 25 : if (pdfMean)
2167 13 : *pdfMean = dfGlobalMean;
2168 25 : if (pdfStdDev)
2169 13 : *pdfStdDev = dfGlobalStdDev;
2170 :
2171 25 : if (sContext.nGlobalValidPixels == 0)
2172 : {
2173 2 : ReportError(CE_Failure, CPLE_AppDefined,
2174 : "Failed to compute statistics, no valid pixels found "
2175 : "in sampling.");
2176 : }
2177 :
2178 25 : 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 242512 : CPLErr VRTSourcedRasterBand::AddSource(VRTSource *poNewSource)
2279 :
2280 : {
2281 242512 : return AddSource(std::unique_ptr<VRTSource>(poNewSource));
2282 : }
2283 :
2284 : /************************************************************************/
2285 : /* AddSource() */
2286 : /************************************************************************/
2287 :
2288 247589 : CPLErr VRTSourcedRasterBand::AddSource(std::unique_ptr<VRTSource> poNewSource)
2289 :
2290 : {
2291 247589 : auto l_poDS = static_cast<VRTDataset *>(poDS);
2292 247589 : l_poDS->SetNeedsFlush();
2293 247589 : l_poDS->SourceAdded();
2294 :
2295 247589 : if (poNewSource->IsSimpleSource())
2296 : {
2297 : VRTSimpleSource *poSS =
2298 247554 : static_cast<VRTSimpleSource *>(poNewSource.get());
2299 247554 : if (GetMetadataItem("NBITS", "IMAGE_STRUCTURE") != nullptr)
2300 : {
2301 47 : int nBits = atoi(GetMetadataItem("NBITS", "IMAGE_STRUCTURE"));
2302 47 : if (nBits >= 1 && nBits <= 31)
2303 : {
2304 47 : poSS->SetMaxValue(static_cast<int>((1U << nBits) - 1));
2305 : }
2306 : }
2307 : }
2308 :
2309 247589 : m_papoSources.push_back(std::move(poNewSource));
2310 :
2311 247589 : return CE_None;
2312 : }
2313 :
2314 : /*! @endcond */
2315 :
2316 : /************************************************************************/
2317 : /* VRTAddSource() */
2318 : /************************************************************************/
2319 :
2320 : /**
2321 : * @see VRTSourcedRasterBand::AddSource().
2322 : */
2323 :
2324 0 : CPLErr CPL_STDCALL VRTAddSource(VRTSourcedRasterBandH hVRTBand,
2325 : VRTSourceH hNewSource)
2326 : {
2327 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddSource", CE_Failure);
2328 :
2329 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddSource(
2330 0 : reinterpret_cast<VRTSource *>(hNewSource));
2331 : }
2332 :
2333 : /*! @cond Doxygen_Suppress */
2334 :
2335 : /************************************************************************/
2336 : /* XMLInit() */
2337 : /************************************************************************/
2338 :
2339 2531 : CPLErr VRTSourcedRasterBand::XMLInit(const CPLXMLNode *psTree,
2340 : const char *pszVRTPath,
2341 : VRTMapSharedResources &oMapSharedSources)
2342 :
2343 : {
2344 : {
2345 : const CPLErr eErr =
2346 2531 : VRTRasterBand::XMLInit(psTree, pszVRTPath, oMapSharedSources);
2347 2531 : if (eErr != CE_None)
2348 0 : return eErr;
2349 : }
2350 :
2351 : /* -------------------------------------------------------------------- */
2352 : /* Process sources. */
2353 : /* -------------------------------------------------------------------- */
2354 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
2355 2531 : GetGDALDriverManager()->GetDriverByName("VRT"));
2356 :
2357 2531 : for (const CPLXMLNode *psChild = psTree->psChild;
2358 118348 : psChild != nullptr && poDriver != nullptr; psChild = psChild->psNext)
2359 : {
2360 115839 : if (psChild->eType != CXT_Element)
2361 6644 : continue;
2362 :
2363 109195 : CPLErrorReset();
2364 : VRTSource *const poSource =
2365 109195 : poDriver->ParseSource(psChild, pszVRTPath, oMapSharedSources);
2366 109195 : if (poSource != nullptr)
2367 104148 : AddSource(poSource);
2368 5047 : else if (CPLGetLastErrorType() != CE_None)
2369 22 : return CE_Failure;
2370 : }
2371 :
2372 : /* -------------------------------------------------------------------- */
2373 : /* Done. */
2374 : /* -------------------------------------------------------------------- */
2375 : const char *pszSubclass =
2376 2509 : CPLGetXMLValue(psTree, "subclass", "VRTSourcedRasterBand");
2377 2509 : if (m_papoSources.empty() && !EQUAL(pszSubclass, "VRTDerivedRasterBand"))
2378 132 : CPLDebug("VRT", "No valid sources found for band in VRT file %s",
2379 132 : GetDataset() ? GetDataset()->GetDescription() : "");
2380 :
2381 2509 : return CE_None;
2382 : }
2383 :
2384 : /************************************************************************/
2385 : /* SerializeToXML() */
2386 : /************************************************************************/
2387 :
2388 695 : CPLXMLNode *VRTSourcedRasterBand::SerializeToXML(const char *pszVRTPath,
2389 : bool &bHasWarnedAboutRAMUsage,
2390 : size_t &nAccRAMUsage)
2391 :
2392 : {
2393 695 : CPLXMLNode *psTree = VRTRasterBand::SerializeToXML(
2394 : pszVRTPath, bHasWarnedAboutRAMUsage, nAccRAMUsage);
2395 695 : CPLXMLNode *psLastChild = psTree->psChild;
2396 2502 : while (psLastChild != nullptr && psLastChild->psNext != nullptr)
2397 1807 : psLastChild = psLastChild->psNext;
2398 :
2399 : /* -------------------------------------------------------------------- */
2400 : /* Process Sources. */
2401 : /* -------------------------------------------------------------------- */
2402 :
2403 695 : GIntBig nUsableRAM = -1;
2404 :
2405 3364 : for (const auto &poSource : m_papoSources)
2406 : {
2407 2669 : CPLXMLNode *const psXMLSrc = poSource->SerializeToXML(pszVRTPath);
2408 :
2409 2669 : if (psXMLSrc == nullptr)
2410 0 : break;
2411 :
2412 : // Creating the CPLXMLNode tree representation of a VRT can easily
2413 : // take several times RAM usage than its string serialization, or its
2414 : // internal representation in the driver.
2415 : // We multiply the estimate by a factor of 2, experimentally found to
2416 : // be more realistic than the conservative raw estimate.
2417 2669 : nAccRAMUsage += 2 * CPLXMLNodeGetRAMUsageEstimate(psXMLSrc);
2418 2669 : if (!bHasWarnedAboutRAMUsage && nAccRAMUsage > 512 * 1024 * 1024)
2419 : {
2420 0 : if (nUsableRAM < 0)
2421 0 : nUsableRAM = CPLGetUsablePhysicalRAM();
2422 0 : if (nUsableRAM > 0 &&
2423 0 : nAccRAMUsage > static_cast<uint64_t>(nUsableRAM) / 10 * 8)
2424 : {
2425 0 : bHasWarnedAboutRAMUsage = true;
2426 0 : CPLError(CE_Warning, CPLE_AppDefined,
2427 : "Serialization of this VRT file has already consumed "
2428 : "at least %.02f GB of RAM over a total of %.02f. This "
2429 : "process may abort",
2430 0 : double(nAccRAMUsage) / (1024 * 1024 * 1024),
2431 0 : double(nUsableRAM) / (1024 * 1024 * 1024));
2432 : }
2433 : }
2434 :
2435 2669 : if (psLastChild == nullptr)
2436 0 : psTree->psChild = psXMLSrc;
2437 : else
2438 2669 : psLastChild->psNext = psXMLSrc;
2439 2669 : psLastChild = psXMLSrc;
2440 : }
2441 :
2442 695 : return psTree;
2443 : }
2444 :
2445 : /************************************************************************/
2446 : /* SkipBufferInitialization() */
2447 : /************************************************************************/
2448 :
2449 279474 : bool VRTSourcedRasterBand::SkipBufferInitialization()
2450 : {
2451 279474 : if (m_nSkipBufferInitialization >= 0)
2452 274276 : return m_nSkipBufferInitialization != 0;
2453 : /* -------------------------------------------------------------------- */
2454 : /* Check if we can avoid buffer initialization. */
2455 : /* -------------------------------------------------------------------- */
2456 :
2457 : // Note: if one day we do alpha compositing, we will need to check that.
2458 5198 : m_nSkipBufferInitialization = FALSE;
2459 5198 : if (m_papoSources.size() != 1 || !m_papoSources[0]->IsSimpleSource())
2460 : {
2461 3517 : return false;
2462 : }
2463 : VRTSimpleSource *poSS =
2464 1681 : static_cast<VRTSimpleSource *>(m_papoSources[0].get());
2465 1681 : if (poSS->GetType() == VRTSimpleSource::GetTypeStatic())
2466 : {
2467 1276 : auto l_poBand = poSS->GetRasterBand();
2468 1265 : if (l_poBand != nullptr && poSS->m_dfSrcXOff >= 0.0 &&
2469 1235 : poSS->m_dfSrcYOff >= 0.0 &&
2470 1233 : poSS->m_dfSrcXOff + poSS->m_dfSrcXSize <= l_poBand->GetXSize() &&
2471 1231 : poSS->m_dfSrcYOff + poSS->m_dfSrcYSize <= l_poBand->GetYSize() &&
2472 1231 : poSS->m_dfDstXOff <= 0.0 && poSS->m_dfDstYOff <= 0.0 &&
2473 3693 : poSS->m_dfDstXOff + poSS->m_dfDstXSize >= nRasterXSize &&
2474 1152 : poSS->m_dfDstYOff + poSS->m_dfDstYSize >= nRasterYSize)
2475 : {
2476 1150 : m_nSkipBufferInitialization = TRUE;
2477 : }
2478 : }
2479 1681 : return m_nSkipBufferInitialization != 0;
2480 : }
2481 :
2482 : /************************************************************************/
2483 : /* ConfigureSource() */
2484 : /************************************************************************/
2485 :
2486 142829 : void VRTSourcedRasterBand::ConfigureSource(VRTSimpleSource *poSimpleSource,
2487 : GDALRasterBand *poSrcBand,
2488 : int bAddAsMaskBand, double dfSrcXOff,
2489 : double dfSrcYOff, double dfSrcXSize,
2490 : double dfSrcYSize, double dfDstXOff,
2491 : double dfDstYOff, double dfDstXSize,
2492 : double dfDstYSize)
2493 : {
2494 : /* -------------------------------------------------------------------- */
2495 : /* Default source and dest rectangles. */
2496 : /* -------------------------------------------------------------------- */
2497 142829 : if (dfSrcYSize == -1)
2498 : {
2499 66962 : dfSrcXOff = 0;
2500 66962 : dfSrcYOff = 0;
2501 66962 : dfSrcXSize = poSrcBand->GetXSize();
2502 66962 : dfSrcYSize = poSrcBand->GetYSize();
2503 : }
2504 :
2505 142829 : if (dfDstYSize == -1)
2506 : {
2507 66962 : dfDstXOff = 0;
2508 66962 : dfDstYOff = 0;
2509 66962 : dfDstXSize = nRasterXSize;
2510 66962 : dfDstYSize = nRasterYSize;
2511 : }
2512 :
2513 142829 : if (bAddAsMaskBand)
2514 48 : poSimpleSource->SetSrcMaskBand(poSrcBand);
2515 : else
2516 142781 : poSimpleSource->SetSrcBand(poSrcBand);
2517 :
2518 142829 : poSimpleSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2519 142829 : poSimpleSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2520 :
2521 : /* -------------------------------------------------------------------- */
2522 : /* If we can get the associated GDALDataset, add a reference to it.*/
2523 : /* -------------------------------------------------------------------- */
2524 142829 : GDALDataset *poSrcBandDataset = poSrcBand->GetDataset();
2525 142829 : if (poSrcBandDataset != nullptr)
2526 : {
2527 : VRTDataset *poVRTSrcBandDataset =
2528 142829 : dynamic_cast<VRTDataset *>(poSrcBandDataset);
2529 142829 : if (poVRTSrcBandDataset && !poVRTSrcBandDataset->m_bCanTakeRef)
2530 : {
2531 : // Situation triggered by VRTDataset::AddVirtualOverview()
2532 : // We create an overview dataset that is a VRT of a reduction of
2533 : // ourselves. But we don't want to take a reference on ourselves,
2534 : // otherwise this will prevent us to be closed in number of
2535 : // circumstances
2536 65 : poSimpleSource->m_bDropRefOnSrcBand = false;
2537 : }
2538 : else
2539 : {
2540 142764 : poSrcBandDataset->Reference();
2541 : }
2542 : }
2543 142829 : }
2544 :
2545 : /************************************************************************/
2546 : /* AddSimpleSource() */
2547 : /************************************************************************/
2548 :
2549 377 : CPLErr VRTSourcedRasterBand::AddSimpleSource(
2550 : const char *pszFilename, int nBandIn, double dfSrcXOff, double dfSrcYOff,
2551 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2552 : double dfDstXSize, double dfDstYSize, const char *pszResampling,
2553 : double dfNoDataValueIn)
2554 :
2555 : {
2556 : /* -------------------------------------------------------------------- */
2557 : /* Create source. */
2558 : /* -------------------------------------------------------------------- */
2559 377 : VRTSimpleSource *poSimpleSource = nullptr;
2560 :
2561 377 : if (pszResampling != nullptr && STARTS_WITH_CI(pszResampling, "aver"))
2562 : {
2563 0 : auto poAveragedSource = new VRTAveragedSource();
2564 0 : poSimpleSource = poAveragedSource;
2565 0 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2566 0 : poAveragedSource->SetNoDataValue(dfNoDataValueIn);
2567 : }
2568 : else
2569 : {
2570 377 : poSimpleSource = new VRTSimpleSource();
2571 377 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2572 0 : CPLError(
2573 : CE_Warning, CPLE_AppDefined,
2574 : "NODATA setting not currently supported for nearest "
2575 : "neighbour sampled simple sources on Virtual Datasources.");
2576 : }
2577 :
2578 377 : poSimpleSource->SetSrcBand(pszFilename, nBandIn);
2579 :
2580 377 : poSimpleSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2581 377 : poSimpleSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2582 :
2583 : /* -------------------------------------------------------------------- */
2584 : /* add to list. */
2585 : /* -------------------------------------------------------------------- */
2586 377 : return AddSource(poSimpleSource);
2587 : }
2588 :
2589 : /************************************************************************/
2590 : /* AddSimpleSource() */
2591 : /************************************************************************/
2592 :
2593 68427 : CPLErr VRTSourcedRasterBand::AddSimpleSource(
2594 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2595 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2596 : double dfDstXSize, double dfDstYSize, const char *pszResampling,
2597 : double dfNoDataValueIn)
2598 :
2599 : {
2600 : /* -------------------------------------------------------------------- */
2601 : /* Create source. */
2602 : /* -------------------------------------------------------------------- */
2603 68427 : VRTSimpleSource *poSimpleSource = nullptr;
2604 :
2605 68427 : if (pszResampling != nullptr && STARTS_WITH_CI(pszResampling, "aver"))
2606 : {
2607 0 : auto poAveragedSource = new VRTAveragedSource();
2608 0 : poSimpleSource = poAveragedSource;
2609 0 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2610 0 : poAveragedSource->SetNoDataValue(dfNoDataValueIn);
2611 : }
2612 : else
2613 : {
2614 68427 : poSimpleSource = new VRTSimpleSource();
2615 68427 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2616 0 : CPLError(
2617 : CE_Warning, CPLE_AppDefined,
2618 : "NODATA setting not currently supported for "
2619 : "neighbour sampled simple sources on Virtual Datasources.");
2620 : }
2621 :
2622 68427 : ConfigureSource(poSimpleSource, poSrcBand, FALSE, dfSrcXOff, dfSrcYOff,
2623 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2624 : dfDstYSize);
2625 :
2626 : /* -------------------------------------------------------------------- */
2627 : /* add to list. */
2628 : /* -------------------------------------------------------------------- */
2629 68427 : return AddSource(poSimpleSource);
2630 : }
2631 :
2632 : /************************************************************************/
2633 : /* AddMaskBandSource() */
2634 : /************************************************************************/
2635 :
2636 : // poSrcBand is not the mask band, but the band from which the mask band is
2637 : // taken.
2638 32 : CPLErr VRTSourcedRasterBand::AddMaskBandSource(
2639 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2640 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2641 : double dfDstXSize, double dfDstYSize)
2642 : {
2643 : /* -------------------------------------------------------------------- */
2644 : /* Create source. */
2645 : /* -------------------------------------------------------------------- */
2646 32 : VRTSimpleSource *poSimpleSource = new VRTSimpleSource();
2647 :
2648 32 : ConfigureSource(poSimpleSource, poSrcBand, TRUE, dfSrcXOff, dfSrcYOff,
2649 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2650 : dfDstYSize);
2651 :
2652 : /* -------------------------------------------------------------------- */
2653 : /* add to list. */
2654 : /* -------------------------------------------------------------------- */
2655 32 : return AddSource(poSimpleSource);
2656 : }
2657 :
2658 : /*! @endcond */
2659 :
2660 : /************************************************************************/
2661 : /* VRTAddSimpleSource() */
2662 : /************************************************************************/
2663 :
2664 : /**
2665 : * @see VRTSourcedRasterBand::AddSimpleSource().
2666 : */
2667 :
2668 968 : CPLErr CPL_STDCALL VRTAddSimpleSource(VRTSourcedRasterBandH hVRTBand,
2669 : GDALRasterBandH hSrcBand, int nSrcXOff,
2670 : int nSrcYOff, int nSrcXSize,
2671 : int nSrcYSize, int nDstXOff, int nDstYOff,
2672 : int nDstXSize, int nDstYSize,
2673 : const char *pszResampling,
2674 : double dfNoDataValue)
2675 : {
2676 968 : VALIDATE_POINTER1(hVRTBand, "VRTAddSimpleSource", CE_Failure);
2677 :
2678 968 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddSimpleSource(
2679 : reinterpret_cast<GDALRasterBand *>(hSrcBand), nSrcXOff, nSrcYOff,
2680 : nSrcXSize, nSrcYSize, nDstXOff, nDstYOff, nDstXSize, nDstYSize,
2681 968 : pszResampling, dfNoDataValue);
2682 : }
2683 :
2684 : /*! @cond Doxygen_Suppress */
2685 :
2686 : /************************************************************************/
2687 : /* AddComplexSource() */
2688 : /************************************************************************/
2689 :
2690 29 : CPLErr VRTSourcedRasterBand::AddComplexSource(
2691 : const char *pszFilename, int nBandIn, double dfSrcXOff, double dfSrcYOff,
2692 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2693 : double dfDstXSize, double dfDstYSize, double dfScaleOff,
2694 : double dfScaleRatio, double dfNoDataValueIn, int nColorTableComponent)
2695 :
2696 : {
2697 : /* -------------------------------------------------------------------- */
2698 : /* Create source. */
2699 : /* -------------------------------------------------------------------- */
2700 29 : VRTComplexSource *const poSource = new VRTComplexSource();
2701 :
2702 29 : poSource->SetSrcBand(pszFilename, nBandIn);
2703 :
2704 29 : poSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2705 29 : poSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2706 :
2707 : /* -------------------------------------------------------------------- */
2708 : /* Set complex parameters. */
2709 : /* -------------------------------------------------------------------- */
2710 29 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2711 8 : poSource->SetNoDataValue(dfNoDataValueIn);
2712 :
2713 29 : if (dfScaleOff != 0.0 || dfScaleRatio != 1.0)
2714 4 : poSource->SetLinearScaling(dfScaleOff, dfScaleRatio);
2715 :
2716 29 : poSource->SetColorTableComponent(nColorTableComponent);
2717 :
2718 : /* -------------------------------------------------------------------- */
2719 : /* add to list. */
2720 : /* -------------------------------------------------------------------- */
2721 29 : return AddSource(poSource);
2722 : }
2723 :
2724 : /************************************************************************/
2725 : /* AddComplexSource() */
2726 : /************************************************************************/
2727 :
2728 30 : CPLErr VRTSourcedRasterBand::AddComplexSource(
2729 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2730 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2731 : double dfDstXSize, double dfDstYSize, double dfScaleOff,
2732 : double dfScaleRatio, double dfNoDataValueIn, int nColorTableComponent)
2733 :
2734 : {
2735 : /* -------------------------------------------------------------------- */
2736 : /* Create source. */
2737 : /* -------------------------------------------------------------------- */
2738 30 : VRTComplexSource *const poSource = new VRTComplexSource();
2739 :
2740 30 : ConfigureSource(poSource, poSrcBand, FALSE, dfSrcXOff, dfSrcYOff,
2741 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2742 : dfDstYSize);
2743 :
2744 : /* -------------------------------------------------------------------- */
2745 : /* Set complex parameters. */
2746 : /* -------------------------------------------------------------------- */
2747 30 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2748 19 : poSource->SetNoDataValue(dfNoDataValueIn);
2749 :
2750 30 : if (dfScaleOff != 0.0 || dfScaleRatio != 1.0)
2751 21 : poSource->SetLinearScaling(dfScaleOff, dfScaleRatio);
2752 :
2753 30 : poSource->SetColorTableComponent(nColorTableComponent);
2754 :
2755 : /* -------------------------------------------------------------------- */
2756 : /* add to list. */
2757 : /* -------------------------------------------------------------------- */
2758 30 : return AddSource(poSource);
2759 : }
2760 :
2761 : /*! @endcond */
2762 :
2763 : /************************************************************************/
2764 : /* VRTAddComplexSource() */
2765 : /************************************************************************/
2766 :
2767 : /**
2768 : * @see VRTSourcedRasterBand::AddComplexSource().
2769 : */
2770 :
2771 0 : CPLErr CPL_STDCALL VRTAddComplexSource(
2772 : VRTSourcedRasterBandH hVRTBand, GDALRasterBandH hSrcBand, int nSrcXOff,
2773 : int nSrcYOff, int nSrcXSize, int nSrcYSize, int nDstXOff, int nDstYOff,
2774 : int nDstXSize, int nDstYSize, double dfScaleOff, double dfScaleRatio,
2775 : double dfNoDataValue)
2776 : {
2777 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddComplexSource", CE_Failure);
2778 :
2779 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddComplexSource(
2780 : reinterpret_cast<GDALRasterBand *>(hSrcBand), nSrcXOff, nSrcYOff,
2781 : nSrcXSize, nSrcYSize, nDstXOff, nDstYOff, nDstXSize, nDstYSize,
2782 0 : dfScaleOff, dfScaleRatio, dfNoDataValue);
2783 : }
2784 :
2785 : /*! @cond Doxygen_Suppress */
2786 :
2787 : /************************************************************************/
2788 : /* AddFuncSource() */
2789 : /************************************************************************/
2790 :
2791 12 : CPLErr VRTSourcedRasterBand::AddFuncSource(VRTImageReadFunc pfnReadFunc,
2792 : void *pCBData,
2793 : double dfNoDataValueIn)
2794 :
2795 : {
2796 : /* -------------------------------------------------------------------- */
2797 : /* Create source. */
2798 : /* -------------------------------------------------------------------- */
2799 12 : VRTFuncSource *const poFuncSource = new VRTFuncSource;
2800 :
2801 12 : poFuncSource->fNoDataValue = static_cast<float>(dfNoDataValueIn);
2802 12 : poFuncSource->pfnReadFunc = pfnReadFunc;
2803 12 : poFuncSource->pCBData = pCBData;
2804 12 : poFuncSource->eType = GetRasterDataType();
2805 :
2806 : /* -------------------------------------------------------------------- */
2807 : /* add to list. */
2808 : /* -------------------------------------------------------------------- */
2809 12 : return AddSource(poFuncSource);
2810 : }
2811 :
2812 : /*! @endcond */
2813 :
2814 : /************************************************************************/
2815 : /* VRTAddFuncSource() */
2816 : /************************************************************************/
2817 :
2818 : /**
2819 : * @see VRTSourcedRasterBand::AddFuncSource().
2820 : */
2821 :
2822 0 : CPLErr CPL_STDCALL VRTAddFuncSource(VRTSourcedRasterBandH hVRTBand,
2823 : VRTImageReadFunc pfnReadFunc, void *pCBData,
2824 : double dfNoDataValue)
2825 : {
2826 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddFuncSource", CE_Failure);
2827 :
2828 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddFuncSource(
2829 0 : pfnReadFunc, pCBData, dfNoDataValue);
2830 : }
2831 :
2832 : /*! @cond Doxygen_Suppress */
2833 :
2834 : /************************************************************************/
2835 : /* GetMetadataDomainList() */
2836 : /************************************************************************/
2837 :
2838 0 : char **VRTSourcedRasterBand::GetMetadataDomainList()
2839 : {
2840 0 : return CSLAddString(GDALRasterBand::GetMetadataDomainList(),
2841 0 : "LocationInfo");
2842 : }
2843 :
2844 : /************************************************************************/
2845 : /* GetMetadataItem() */
2846 : /************************************************************************/
2847 :
2848 256279 : const char *VRTSourcedRasterBand::GetMetadataItem(const char *pszName,
2849 : const char *pszDomain)
2850 :
2851 : {
2852 : /* ==================================================================== */
2853 : /* LocationInfo handling. */
2854 : /* ==================================================================== */
2855 256279 : if (pszDomain != nullptr && EQUAL(pszDomain, "LocationInfo") &&
2856 3 : (STARTS_WITH_CI(pszName, "Pixel_") ||
2857 0 : STARTS_WITH_CI(pszName, "GeoPixel_")))
2858 : {
2859 : /* --------------------------------------------------------------------
2860 : */
2861 : /* What pixel are we aiming at? */
2862 : /* --------------------------------------------------------------------
2863 : */
2864 3 : int iPixel = 0;
2865 3 : int iLine = 0;
2866 :
2867 3 : if (STARTS_WITH_CI(pszName, "Pixel_"))
2868 : {
2869 : // TODO(schwehr): Replace sscanf.
2870 3 : if (sscanf(pszName + 6, "%d_%d", &iPixel, &iLine) != 2)
2871 0 : return nullptr;
2872 : }
2873 0 : else if (STARTS_WITH_CI(pszName, "GeoPixel_"))
2874 : {
2875 0 : const double dfGeoX = CPLAtof(pszName + 9);
2876 0 : const char *const pszUnderscore = strchr(pszName + 9, '_');
2877 0 : if (!pszUnderscore)
2878 0 : return nullptr;
2879 0 : const double dfGeoY = CPLAtof(pszUnderscore + 1);
2880 :
2881 0 : if (GetDataset() == nullptr)
2882 0 : return nullptr;
2883 :
2884 0 : GDALGeoTransform gt, invGT;
2885 0 : if (GetDataset()->GetGeoTransform(gt) != CE_None ||
2886 0 : !gt.GetInverse(invGT))
2887 : {
2888 0 : return nullptr;
2889 : }
2890 :
2891 0 : iPixel = static_cast<int>(
2892 0 : floor(invGT[0] + invGT[1] * dfGeoX + invGT[2] * dfGeoY));
2893 0 : iLine = static_cast<int>(
2894 0 : floor(invGT[3] + invGT[4] * dfGeoX + invGT[5] * dfGeoY));
2895 : }
2896 : else
2897 : {
2898 0 : return nullptr;
2899 : }
2900 :
2901 6 : if (iPixel < 0 || iLine < 0 || iPixel >= GetXSize() ||
2902 3 : iLine >= GetYSize())
2903 0 : return nullptr;
2904 :
2905 : /* --------------------------------------------------------------------
2906 : */
2907 : /* Find the file(s) at this location. */
2908 : /* --------------------------------------------------------------------
2909 : */
2910 3 : char **papszFileList = nullptr;
2911 3 : int nListSize = 0; // keep it in this scope
2912 3 : int nListMaxSize = 0; // keep it in this scope
2913 : CPLHashSet *const hSetFiles =
2914 3 : CPLHashSetNew(CPLHashSetHashStr, CPLHashSetEqualStr, nullptr);
2915 :
2916 6 : for (auto &poSource : m_papoSources)
2917 : {
2918 3 : if (!poSource->IsSimpleSource())
2919 0 : continue;
2920 :
2921 : VRTSimpleSource *const poSrc =
2922 3 : static_cast<VRTSimpleSource *>(poSource.get());
2923 :
2924 3 : double dfReqXOff = 0.0;
2925 3 : double dfReqYOff = 0.0;
2926 3 : double dfReqXSize = 0.0;
2927 3 : double dfReqYSize = 0.0;
2928 3 : int nReqXOff = 0;
2929 3 : int nReqYOff = 0;
2930 3 : int nReqXSize = 0;
2931 3 : int nReqYSize = 0;
2932 3 : int nOutXOff = 0;
2933 3 : int nOutYOff = 0;
2934 3 : int nOutXSize = 0;
2935 3 : int nOutYSize = 0;
2936 :
2937 3 : bool bError = false;
2938 3 : if (!poSrc->GetSrcDstWindow(
2939 : iPixel, iLine, 1, 1, 1, 1, GRIORA_NearestNeighbour,
2940 : &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
2941 : &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
2942 : &nOutXSize, &nOutYSize, bError))
2943 : {
2944 0 : if (bError)
2945 : {
2946 0 : CSLDestroy(papszFileList);
2947 0 : CPLHashSetDestroy(hSetFiles);
2948 0 : return nullptr;
2949 : }
2950 0 : continue;
2951 : }
2952 :
2953 3 : poSrc->GetFileList(&papszFileList, &nListSize, &nListMaxSize,
2954 3 : hSetFiles);
2955 : }
2956 :
2957 : /* --------------------------------------------------------------------
2958 : */
2959 : /* Format into XML. */
2960 : /* --------------------------------------------------------------------
2961 : */
2962 3 : m_osLastLocationInfo = "<LocationInfo>";
2963 6 : for (int i = 0; i < nListSize && papszFileList[i] != nullptr; i++)
2964 : {
2965 3 : m_osLastLocationInfo += "<File>";
2966 : char *const pszXMLEscaped =
2967 3 : CPLEscapeString(papszFileList[i], -1, CPLES_XML);
2968 3 : m_osLastLocationInfo += pszXMLEscaped;
2969 3 : CPLFree(pszXMLEscaped);
2970 3 : m_osLastLocationInfo += "</File>";
2971 : }
2972 3 : m_osLastLocationInfo += "</LocationInfo>";
2973 :
2974 3 : CSLDestroy(papszFileList);
2975 3 : CPLHashSetDestroy(hSetFiles);
2976 :
2977 3 : return m_osLastLocationInfo.c_str();
2978 : }
2979 :
2980 : /* ==================================================================== */
2981 : /* Other domains. */
2982 : /* ==================================================================== */
2983 :
2984 256276 : return GDALRasterBand::GetMetadataItem(pszName, pszDomain);
2985 : }
2986 :
2987 : /************************************************************************/
2988 : /* GetMetadata() */
2989 : /************************************************************************/
2990 :
2991 14548 : CSLConstList VRTSourcedRasterBand::GetMetadata(const char *pszDomain)
2992 :
2993 : {
2994 : /* ==================================================================== */
2995 : /* vrt_sources domain handling. */
2996 : /* ==================================================================== */
2997 14548 : if (pszDomain != nullptr && EQUAL(pszDomain, "vrt_sources"))
2998 : {
2999 1 : if (static_cast<size_t>(m_aosSourceList.size()) != m_papoSources.size())
3000 : {
3001 1 : m_aosSourceList.clear();
3002 :
3003 : // Process SimpleSources
3004 2 : for (int iSource = 0;
3005 2 : iSource < static_cast<int>(m_papoSources.size()); iSource++)
3006 : {
3007 : CPLXMLNode *const psXMLSrc =
3008 1 : m_papoSources[iSource]->SerializeToXML(nullptr);
3009 1 : if (psXMLSrc == nullptr)
3010 0 : continue;
3011 :
3012 1 : char *const pszXML = CPLSerializeXMLTree(psXMLSrc);
3013 :
3014 : m_aosSourceList.AddString(
3015 1 : CPLSPrintf("source_%d=%s", iSource, pszXML));
3016 1 : CPLFree(pszXML);
3017 1 : CPLDestroyXMLNode(psXMLSrc);
3018 : }
3019 : }
3020 1 : return m_aosSourceList.List();
3021 : }
3022 :
3023 : /* ==================================================================== */
3024 : /* Other domains. */
3025 : /* ==================================================================== */
3026 :
3027 14547 : return GDALRasterBand::GetMetadata(pszDomain);
3028 : }
3029 :
3030 : /************************************************************************/
3031 : /* SetMetadataItem() */
3032 : /************************************************************************/
3033 :
3034 133607 : CPLErr VRTSourcedRasterBand::SetMetadataItem(const char *pszName,
3035 : const char *pszValue,
3036 : const char *pszDomain)
3037 :
3038 : {
3039 : #if DEBUG_VERBOSE
3040 : CPLDebug("VRT", "VRTSourcedRasterBand::SetMetadataItem(%s,%s,%s)\n",
3041 : pszName, pszValue ? pszValue : "(null)",
3042 : pszDomain ? pszDomain : "(null)");
3043 : #endif
3044 :
3045 133607 : if (pszDomain != nullptr && EQUAL(pszDomain, "new_vrt_sources"))
3046 : {
3047 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
3048 1 : GetGDALDriverManager()->GetDriverByName("VRT"));
3049 1 : if (poDriver)
3050 : {
3051 1 : CPLXMLNode *const psTree = CPLParseXMLString(pszValue);
3052 1 : if (psTree == nullptr)
3053 0 : return CE_Failure;
3054 :
3055 1 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
3056 1 : if (l_poDS == nullptr)
3057 : {
3058 0 : CPLDestroyXMLNode(psTree);
3059 0 : return CE_Failure;
3060 : }
3061 2 : VRTSource *const poSource = poDriver->ParseSource(
3062 1 : psTree, nullptr, l_poDS->m_oMapSharedSources);
3063 1 : CPLDestroyXMLNode(psTree);
3064 :
3065 1 : if (poSource != nullptr)
3066 1 : return AddSource(poSource);
3067 : }
3068 :
3069 0 : return CE_Failure;
3070 : }
3071 133606 : else if (pszDomain != nullptr && EQUAL(pszDomain, "vrt_sources"))
3072 : {
3073 1 : int iSource = 0;
3074 : // TODO(schwehr): Replace sscanf.
3075 2 : if (sscanf(pszName, "source_%d", &iSource) != 1 || iSource < 0 ||
3076 1 : iSource >= static_cast<int>(m_papoSources.size()))
3077 : {
3078 0 : CPLError(CE_Failure, CPLE_AppDefined,
3079 : "%s metadata item name is not recognized. "
3080 : "Should be between source_0 and source_%d",
3081 0 : pszName, static_cast<int>(m_papoSources.size()) - 1);
3082 0 : return CE_Failure;
3083 : }
3084 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
3085 1 : GetGDALDriverManager()->GetDriverByName("VRT"));
3086 1 : if (poDriver)
3087 : {
3088 1 : CPLXMLNode *const psTree = CPLParseXMLString(pszValue);
3089 1 : if (psTree == nullptr)
3090 1 : return CE_Failure;
3091 :
3092 1 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
3093 1 : if (l_poDS == nullptr)
3094 : {
3095 0 : CPLDestroyXMLNode(psTree);
3096 0 : return CE_Failure;
3097 : }
3098 : auto poSource = std::unique_ptr<VRTSource>(poDriver->ParseSource(
3099 1 : psTree, nullptr, l_poDS->m_oMapSharedSources));
3100 1 : CPLDestroyXMLNode(psTree);
3101 :
3102 1 : if (poSource != nullptr)
3103 : {
3104 1 : m_papoSources[iSource] = std::move(poSource);
3105 1 : static_cast<VRTDataset *>(poDS)->SetNeedsFlush();
3106 1 : return CE_None;
3107 : }
3108 : }
3109 0 : return CE_Failure;
3110 : }
3111 :
3112 133605 : return VRTRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
3113 : }
3114 :
3115 : /************************************************************************/
3116 : /* SetMetadata() */
3117 : /************************************************************************/
3118 :
3119 70795 : CPLErr VRTSourcedRasterBand::SetMetadata(CSLConstList papszNewMD,
3120 : const char *pszDomain)
3121 :
3122 : {
3123 70795 : if (pszDomain != nullptr && (EQUAL(pszDomain, "new_vrt_sources") ||
3124 70795 : EQUAL(pszDomain, "vrt_sources")))
3125 : {
3126 0 : VRTDriver *const poDriver = dynamic_cast<VRTDriver *>(
3127 8 : GetGDALDriverManager()->GetDriverByName("VRT"));
3128 8 : if (!poDriver)
3129 0 : return CE_Failure;
3130 :
3131 8 : if (EQUAL(pszDomain, "vrt_sources"))
3132 : {
3133 8 : m_papoSources.clear();
3134 : }
3135 :
3136 6 : for (const char *const pszMDItem :
3137 20 : cpl::Iterate(CSLConstList(papszNewMD)))
3138 : {
3139 8 : const char *const pszXML = CPLParseNameValue(pszMDItem, nullptr);
3140 8 : CPLXMLTreeCloser psTree(CPLParseXMLString(pszXML));
3141 8 : if (!psTree)
3142 0 : return CE_Failure;
3143 :
3144 8 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
3145 8 : if (l_poDS == nullptr)
3146 : {
3147 0 : return CE_Failure;
3148 : }
3149 8 : VRTSource *const poSource = poDriver->ParseSource(
3150 8 : psTree.get(), nullptr, l_poDS->m_oMapSharedSources);
3151 8 : if (poSource == nullptr)
3152 2 : return CE_Failure;
3153 :
3154 6 : const CPLErr eErr = AddSource(poSource);
3155 : // cppcheck-suppress knownConditionTrueFalse
3156 6 : if (eErr != CE_None)
3157 0 : return eErr;
3158 : }
3159 :
3160 6 : return CE_None;
3161 : }
3162 :
3163 70787 : return VRTRasterBand::SetMetadata(papszNewMD, pszDomain);
3164 : }
3165 :
3166 : /************************************************************************/
3167 : /* GetFileList() */
3168 : /************************************************************************/
3169 :
3170 68 : void VRTSourcedRasterBand::GetFileList(char ***ppapszFileList, int *pnSize,
3171 : int *pnMaxSize, CPLHashSet *hSetFiles)
3172 : {
3173 179 : for (auto &poSource : m_papoSources)
3174 : {
3175 111 : poSource->GetFileList(ppapszFileList, pnSize, pnMaxSize, hSetFiles);
3176 : }
3177 :
3178 68 : VRTRasterBand::GetFileList(ppapszFileList, pnSize, pnMaxSize, hSetFiles);
3179 68 : }
3180 :
3181 : /************************************************************************/
3182 : /* CloseDependentDatasets() */
3183 : /************************************************************************/
3184 :
3185 142177 : int VRTSourcedRasterBand::CloseDependentDatasets()
3186 : {
3187 142177 : int ret = VRTRasterBand::CloseDependentDatasets();
3188 :
3189 142177 : if (m_papoSources.empty())
3190 389 : return ret;
3191 :
3192 141788 : m_papoSources.clear();
3193 :
3194 141788 : return TRUE;
3195 : }
3196 :
3197 : /************************************************************************/
3198 : /* FlushCache() */
3199 : /************************************************************************/
3200 :
3201 208343 : CPLErr VRTSourcedRasterBand::FlushCache(bool bAtClosing)
3202 : {
3203 208343 : CPLErr eErr = VRTRasterBand::FlushCache(bAtClosing);
3204 522123 : for (size_t i = 0; i < m_papoSources.size() && eErr == CE_None; i++)
3205 : {
3206 313780 : eErr = m_papoSources[i]->FlushCache(bAtClosing);
3207 : }
3208 208343 : return eErr;
3209 : }
3210 :
3211 : /************************************************************************/
3212 : /* RemoveCoveredSources() */
3213 : /************************************************************************/
3214 :
3215 : /** Remove sources that are covered by other sources.
3216 : *
3217 : * This method removes sources that are covered entirely by (one or several)
3218 : * sources of higher priority (even if they declare a nodata setting).
3219 : * This optimizes the size of the VRT and the rendering time.
3220 : */
3221 14 : void VRTSourcedRasterBand::RemoveCoveredSources(CSLConstList papszOptions)
3222 : {
3223 : #ifndef HAVE_GEOS
3224 : if (CPLTestBool(CSLFetchNameValueDef(
3225 : papszOptions, "EMIT_ERROR_IF_GEOS_NOT_AVAILABLE", "TRUE")))
3226 : {
3227 : CPLError(CE_Failure, CPLE_NotSupported,
3228 : "RemoveCoveredSources() not implemented in builds "
3229 : "without GEOS support");
3230 : }
3231 : #else
3232 : (void)papszOptions;
3233 :
3234 : CPLRectObj globalBounds;
3235 14 : globalBounds.minx = 0;
3236 14 : globalBounds.miny = 0;
3237 14 : globalBounds.maxx = nRasterXSize;
3238 14 : globalBounds.maxy = nRasterYSize;
3239 :
3240 : // Create an index with the bbox of all sources
3241 14 : CPLQuadTree *hTree = CPLQuadTreeCreate(&globalBounds, nullptr);
3242 38 : for (int i = 0; i < static_cast<int>(m_papoSources.size()); i++)
3243 : {
3244 24 : if (m_papoSources[i]->IsSimpleSource())
3245 : {
3246 : VRTSimpleSource *poSS =
3247 24 : cpl::down_cast<VRTSimpleSource *>(m_papoSources[i].get());
3248 24 : void *hFeature =
3249 24 : reinterpret_cast<void *>(static_cast<uintptr_t>(i));
3250 : CPLRectObj rect;
3251 24 : rect.minx = std::max(0.0, poSS->m_dfDstXOff);
3252 24 : rect.miny = std::max(0.0, poSS->m_dfDstYOff);
3253 48 : rect.maxx = std::min(double(nRasterXSize),
3254 24 : poSS->m_dfDstXOff + poSS->m_dfDstXSize);
3255 48 : rect.maxy = std::min(double(nRasterYSize),
3256 24 : poSS->m_dfDstYOff + poSS->m_dfDstYSize);
3257 24 : CPLQuadTreeInsertWithBounds(hTree, hFeature, &rect);
3258 : }
3259 : }
3260 :
3261 38 : for (int i = 0; i < static_cast<int>(m_papoSources.size()); i++)
3262 : {
3263 24 : if (m_papoSources[i]->IsSimpleSource())
3264 : {
3265 : VRTSimpleSource *poSS =
3266 24 : cpl::down_cast<VRTSimpleSource *>(m_papoSources[i].get());
3267 : CPLRectObj rect;
3268 24 : rect.minx = std::max(0.0, poSS->m_dfDstXOff);
3269 24 : rect.miny = std::max(0.0, poSS->m_dfDstYOff);
3270 48 : rect.maxx = std::min(double(nRasterXSize),
3271 24 : poSS->m_dfDstXOff + poSS->m_dfDstXSize);
3272 48 : rect.maxy = std::min(double(nRasterYSize),
3273 24 : poSS->m_dfDstYOff + poSS->m_dfDstYSize);
3274 :
3275 : // Find sources whose extent intersect with the current one
3276 24 : int nFeatureCount = 0;
3277 : void **pahFeatures =
3278 24 : CPLQuadTreeSearch(hTree, &rect, &nFeatureCount);
3279 :
3280 : // Compute the bounding box of those sources, only if they are
3281 : // on top of the current one
3282 : CPLRectObj rectIntersecting;
3283 24 : rectIntersecting.minx = std::numeric_limits<double>::max();
3284 24 : rectIntersecting.miny = std::numeric_limits<double>::max();
3285 24 : rectIntersecting.maxx = -std::numeric_limits<double>::max();
3286 24 : rectIntersecting.maxy = -std::numeric_limits<double>::max();
3287 61 : for (int j = 0; j < nFeatureCount; j++)
3288 : {
3289 37 : const int curFeature = static_cast<int>(
3290 37 : reinterpret_cast<uintptr_t>(pahFeatures[j]));
3291 37 : if (curFeature > i)
3292 : {
3293 : VRTSimpleSource *poOtherSS =
3294 13 : cpl::down_cast<VRTSimpleSource *>(
3295 13 : m_papoSources[curFeature].get());
3296 13 : rectIntersecting.minx =
3297 13 : std::min(rectIntersecting.minx, poOtherSS->m_dfDstXOff);
3298 13 : rectIntersecting.miny =
3299 13 : std::min(rectIntersecting.miny, poOtherSS->m_dfDstYOff);
3300 13 : rectIntersecting.maxx = std::max(
3301 : rectIntersecting.maxx,
3302 13 : poOtherSS->m_dfDstXOff + poOtherSS->m_dfDstXSize);
3303 13 : rectIntersecting.maxy = std::max(
3304 : rectIntersecting.maxy,
3305 13 : poOtherSS->m_dfDstYOff + poOtherSS->m_dfDstXSize);
3306 : }
3307 : }
3308 :
3309 : // If the boundinx box of those sources overlap the current one,
3310 : // then compute their union, and check if it contains the current
3311 : // source
3312 24 : if (rectIntersecting.minx <= rect.minx &&
3313 7 : rectIntersecting.miny <= rect.miny &&
3314 7 : rectIntersecting.maxx >= rect.maxx &&
3315 7 : rectIntersecting.maxy >= rect.maxy)
3316 : {
3317 14 : OGRPolygon oPoly;
3318 : {
3319 7 : auto poLR = new OGRLinearRing();
3320 7 : poLR->addPoint(rect.minx, rect.miny);
3321 7 : poLR->addPoint(rect.minx, rect.maxy);
3322 7 : poLR->addPoint(rect.maxx, rect.maxy);
3323 7 : poLR->addPoint(rect.maxx, rect.miny);
3324 7 : poLR->addPoint(rect.minx, rect.miny);
3325 7 : oPoly.addRingDirectly(poLR);
3326 : }
3327 :
3328 7 : std::unique_ptr<OGRGeometry> poUnion;
3329 24 : for (int j = 0; j < nFeatureCount; j++)
3330 : {
3331 17 : const int curFeature = static_cast<int>(
3332 17 : reinterpret_cast<uintptr_t>(pahFeatures[j]));
3333 17 : if (curFeature > i)
3334 : {
3335 : VRTSimpleSource *poOtherSS =
3336 10 : cpl::down_cast<VRTSimpleSource *>(
3337 10 : m_papoSources[curFeature].get());
3338 : CPLRectObj otherRect;
3339 10 : otherRect.minx = std::max(0.0, poOtherSS->m_dfDstXOff);
3340 10 : otherRect.miny = std::max(0.0, poOtherSS->m_dfDstYOff);
3341 20 : otherRect.maxx = std::min(double(nRasterXSize),
3342 20 : poOtherSS->m_dfDstXOff +
3343 10 : poOtherSS->m_dfDstXSize);
3344 20 : otherRect.maxy = std::min(double(nRasterYSize),
3345 20 : poOtherSS->m_dfDstYOff +
3346 10 : poOtherSS->m_dfDstYSize);
3347 20 : OGRPolygon oOtherPoly;
3348 : {
3349 10 : auto poLR = new OGRLinearRing();
3350 10 : poLR->addPoint(otherRect.minx, otherRect.miny);
3351 10 : poLR->addPoint(otherRect.minx, otherRect.maxy);
3352 10 : poLR->addPoint(otherRect.maxx, otherRect.maxy);
3353 10 : poLR->addPoint(otherRect.maxx, otherRect.miny);
3354 10 : poLR->addPoint(otherRect.minx, otherRect.miny);
3355 10 : oOtherPoly.addRingDirectly(poLR);
3356 : }
3357 10 : if (poUnion == nullptr)
3358 7 : poUnion.reset(oOtherPoly.clone());
3359 : else
3360 3 : poUnion.reset(oOtherPoly.Union(poUnion.get()));
3361 : }
3362 : }
3363 :
3364 7 : if (poUnion != nullptr && poUnion->Contains(&oPoly))
3365 : {
3366 : // We can remove the current source
3367 7 : m_papoSources[i].reset();
3368 : }
3369 : }
3370 24 : CPLFree(pahFeatures);
3371 :
3372 24 : void *hFeature =
3373 24 : reinterpret_cast<void *>(static_cast<uintptr_t>(i));
3374 24 : CPLQuadTreeRemove(hTree, hFeature, &rect);
3375 : }
3376 : }
3377 :
3378 : // Compact the papoSources array
3379 0 : m_papoSources.erase(std::remove_if(m_papoSources.begin(),
3380 : m_papoSources.end(),
3381 24 : [](const std::unique_ptr<VRTSource> &src)
3382 38 : { return src.get() == nullptr; }),
3383 28 : m_papoSources.end());
3384 :
3385 14 : CPLQuadTreeDestroy(hTree);
3386 : #endif
3387 14 : }
3388 :
3389 : /*! @endcond */
|
