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