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