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