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