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