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