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