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