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