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