Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Virtual GDAL Datasets
4 : * Purpose: Implementation of VRTPansharpenedRasterBand and
5 : *VRTPansharpenedDataset. Author: Even Rouault <even.rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2015, Even Rouault <even.rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_port.h"
14 : #include "gdal_vrt.h"
15 : #include "vrtdataset.h"
16 :
17 : #include <cassert>
18 : #include <cmath>
19 : #include <cstddef>
20 : #include <cstdlib>
21 : #include <cstring>
22 :
23 : #include <algorithm>
24 : #include <limits>
25 : #include <map>
26 : #include <set>
27 : #include <string>
28 : #include <utility>
29 : #include <vector>
30 :
31 : #include "cpl_conv.h"
32 : #include "cpl_error.h"
33 : #include "cpl_minixml.h"
34 : #include "cpl_string.h"
35 : #include "cpl_vsi.h"
36 : #include "gdal.h"
37 : #include "gdal_priv.h"
38 : #include "gdalpansharpen.h"
39 : #include "ogr_core.h"
40 : #include "ogr_spatialref.h"
41 :
42 : /************************************************************************/
43 : /* GDALCreatePansharpenedVRT() */
44 : /************************************************************************/
45 :
46 : /**
47 : * Create a virtual pansharpened dataset.
48 : *
49 : * This function will create a virtual pansharpened dataset.
50 : *
51 : * Starting with GDAL 3.12, a reference will be taken on the dataset owing
52 : * each input bands.
53 : * Before 3.12, no reference were taken on the passed bands. Consequently,
54 : * they or their dataset to which they belong had to be kept open until
55 : * this virtual pansharpened dataset is closed.
56 : *
57 : * The returned dataset will have no associated filename for itself. If you
58 : * want to write the virtual dataset description to a file, use the
59 : * GDALSetDescription() function (or SetDescription() method) on the dataset
60 : * to assign a filename before it is closed.
61 : *
62 : * @param pszXML Pansharpened VRT XML where <SpectralBand> elements have
63 : * no explicit SourceFilename and SourceBand. The spectral bands in the XML will
64 : * be assigned the successive values of the pahInputSpectralBands array. Must
65 : * not be NULL.
66 : *
67 : * @param hPanchroBand Panchromatic band. Must not be NULL.
68 : *
69 : * @param nInputSpectralBands Number of input spectral bands. Must be greater
70 : * than zero.
71 : *
72 : * @param pahInputSpectralBands Array of nInputSpectralBands spectral bands.
73 : *
74 : * @return NULL on failure, or a new virtual dataset handle on success to be
75 : * closed with GDALClose().
76 : *
77 : * @since GDAL 2.1
78 : */
79 :
80 16 : GDALDatasetH GDALCreatePansharpenedVRT(const char *pszXML,
81 : GDALRasterBandH hPanchroBand,
82 : int nInputSpectralBands,
83 : GDALRasterBandH *pahInputSpectralBands)
84 : {
85 16 : VALIDATE_POINTER1(pszXML, "GDALCreatePansharpenedVRT", nullptr);
86 16 : VALIDATE_POINTER1(hPanchroBand, "GDALCreatePansharpenedVRT", nullptr);
87 16 : VALIDATE_POINTER1(pahInputSpectralBands, "GDALCreatePansharpenedVRT",
88 : nullptr);
89 :
90 16 : CPLXMLNode *psTree = CPLParseXMLString(pszXML);
91 16 : if (psTree == nullptr)
92 1 : return nullptr;
93 15 : VRTPansharpenedDataset *poDS = new VRTPansharpenedDataset(0, 0);
94 15 : CPLErr eErr = poDS->XMLInit(psTree, nullptr, hPanchroBand,
95 : nInputSpectralBands, pahInputSpectralBands);
96 15 : CPLDestroyXMLNode(psTree);
97 15 : if (eErr != CE_None)
98 : {
99 3 : delete poDS;
100 3 : return nullptr;
101 : }
102 12 : return GDALDataset::ToHandle(poDS);
103 : }
104 :
105 : /*! @cond Doxygen_Suppress */
106 :
107 : /************************************************************************/
108 : /* ==================================================================== */
109 : /* VRTPansharpenedDataset */
110 : /* ==================================================================== */
111 : /************************************************************************/
112 :
113 : /************************************************************************/
114 : /* VRTPansharpenedDataset() */
115 : /************************************************************************/
116 :
117 92 : VRTPansharpenedDataset::VRTPansharpenedDataset(int nXSize, int nYSize,
118 92 : int nBlockXSize, int nBlockYSize)
119 : : VRTDataset(nXSize, nYSize,
120 92 : nBlockXSize > 0 ? nBlockXSize : std::min(nXSize, 512),
121 92 : nBlockYSize > 0 ? nBlockYSize : std::min(nYSize, 512)),
122 : m_poPansharpener(nullptr), m_poMainDataset(nullptr),
123 : m_bLoadingOtherBands(FALSE), m_pabyLastBufferBandRasterIO(nullptr),
124 : m_nLastBandRasterIOXOff(0), m_nLastBandRasterIOYOff(0),
125 : m_nLastBandRasterIOXSize(0), m_nLastBandRasterIOYSize(0),
126 : m_eLastBandRasterIODataType(GDT_Unknown), m_eGTAdjustment(GTAdjust_Union),
127 276 : m_bNoDataDisabled(FALSE)
128 : {
129 92 : eAccess = GA_Update;
130 92 : m_poMainDataset = this;
131 92 : }
132 :
133 : /************************************************************************/
134 : /* ~VRTPansharpenedDataset() */
135 : /************************************************************************/
136 :
137 184 : VRTPansharpenedDataset::~VRTPansharpenedDataset()
138 :
139 : {
140 92 : VRTPansharpenedDataset::FlushCache(true);
141 92 : VRTPansharpenedDataset::CloseDependentDatasets();
142 92 : CPLFree(m_pabyLastBufferBandRasterIO);
143 184 : }
144 :
145 : /************************************************************************/
146 : /* CloseDependentDatasets() */
147 : /************************************************************************/
148 :
149 95 : int VRTPansharpenedDataset::CloseDependentDatasets()
150 : {
151 95 : if (m_poMainDataset == nullptr)
152 3 : return FALSE;
153 :
154 92 : VRTPansharpenedDataset *poMainDatasetLocal = m_poMainDataset;
155 92 : m_poMainDataset = nullptr;
156 92 : int bHasDroppedRef = VRTDataset::CloseDependentDatasets();
157 :
158 : /* -------------------------------------------------------------------- */
159 : /* Destroy the raster bands if they exist. */
160 : /* -------------------------------------------------------------------- */
161 288 : for (int iBand = 0; iBand < nBands; iBand++)
162 : {
163 196 : delete papoBands[iBand];
164 : }
165 92 : nBands = 0;
166 :
167 : // Destroy the overviews before m_poPansharpener as they might reference
168 : // files that are in m_apoDatasetsToReleaseRef.
169 92 : bHasDroppedRef |= !m_apoOverviewDatasets.empty();
170 92 : m_apoOverviewDatasets.clear();
171 :
172 92 : if (m_poPansharpener != nullptr)
173 : {
174 : // Delete the pansharper object before closing the dataset
175 : // because it may have warped the bands into an intermediate VRT
176 69 : delete m_poPansharpener;
177 69 : m_poPansharpener = nullptr;
178 :
179 : // Close in reverse order (VRT firsts and real datasets after)
180 324 : for (int i = static_cast<int>(m_apoDatasetsToReleaseRef.size()) - 1;
181 324 : i >= 0; i--)
182 : {
183 255 : bHasDroppedRef = TRUE;
184 255 : m_apoDatasetsToReleaseRef[i].reset();
185 : }
186 69 : m_apoDatasetsToReleaseRef.clear();
187 : }
188 :
189 92 : if (poMainDatasetLocal != this)
190 : {
191 : // To avoid killing us
192 3 : for (size_t i = 0; i < poMainDatasetLocal->m_apoOverviewDatasets.size();
193 : i++)
194 : {
195 3 : if (poMainDatasetLocal->m_apoOverviewDatasets[i].get() == this)
196 : {
197 : // coverity[leaked_storage]
198 3 : CPL_IGNORE_RET_VAL(
199 3 : poMainDatasetLocal->m_apoOverviewDatasets[i].release());
200 3 : break;
201 : }
202 : }
203 3 : bHasDroppedRef |= poMainDatasetLocal->CloseDependentDatasets();
204 : }
205 :
206 92 : return bHasDroppedRef;
207 : }
208 :
209 : /************************************************************************/
210 : /* GetFileList() */
211 : /************************************************************************/
212 :
213 1 : char **VRTPansharpenedDataset::GetFileList()
214 : {
215 1 : char **papszFileList = GDALDataset::GetFileList();
216 :
217 1 : if (m_poPansharpener != nullptr)
218 : {
219 1 : GDALPansharpenOptions *psOptions = m_poPansharpener->GetOptions();
220 1 : if (psOptions != nullptr)
221 : {
222 2 : std::set<CPLString> oSetNames;
223 1 : if (psOptions->hPanchroBand != nullptr)
224 : {
225 1 : GDALDatasetH hDS = GDALGetBandDataset(psOptions->hPanchroBand);
226 1 : if (hDS != nullptr)
227 : {
228 : papszFileList =
229 1 : CSLAddString(papszFileList, GDALGetDescription(hDS));
230 1 : oSetNames.insert(GDALGetDescription(hDS));
231 : }
232 : }
233 4 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
234 : {
235 3 : if (psOptions->pahInputSpectralBands[i] != nullptr)
236 : {
237 : GDALDatasetH hDS =
238 3 : GDALGetBandDataset(psOptions->pahInputSpectralBands[i]);
239 6 : if (hDS != nullptr && oSetNames.find(GDALGetDescription(
240 6 : hDS)) == oSetNames.end())
241 : {
242 1 : papszFileList = CSLAddString(papszFileList,
243 : GDALGetDescription(hDS));
244 1 : oSetNames.insert(GDALGetDescription(hDS));
245 : }
246 : }
247 : }
248 : }
249 : }
250 :
251 1 : return papszFileList;
252 : }
253 :
254 : /************************************************************************/
255 : /* XMLInit() */
256 : /************************************************************************/
257 :
258 74 : CPLErr VRTPansharpenedDataset::XMLInit(const CPLXMLNode *psTree,
259 : const char *pszVRTPathIn)
260 :
261 : {
262 74 : return XMLInit(psTree, pszVRTPathIn, nullptr, 0, nullptr);
263 : }
264 :
265 89 : CPLErr VRTPansharpenedDataset::XMLInit(const CPLXMLNode *psTree,
266 : const char *pszVRTPathIn,
267 : GDALRasterBandH hPanchroBandIn,
268 : int nInputSpectralBandsIn,
269 : GDALRasterBandH *pahInputSpectralBandsIn)
270 : {
271 : CPLErr eErr;
272 : GDALPansharpenOptions *psPanOptions;
273 :
274 : /* -------------------------------------------------------------------- */
275 : /* Initialize blocksize before calling sub-init so that the */
276 : /* band initializers can get it from the dataset object when */
277 : /* they are created. */
278 : /* -------------------------------------------------------------------- */
279 89 : m_nBlockXSize = atoi(CPLGetXMLValue(psTree, "BlockXSize", "512"));
280 89 : m_nBlockYSize = atoi(CPLGetXMLValue(psTree, "BlockYSize", "512"));
281 :
282 : /* -------------------------------------------------------------------- */
283 : /* Parse PansharpeningOptions */
284 : /* -------------------------------------------------------------------- */
285 :
286 89 : const CPLXMLNode *psOptions = CPLGetXMLNode(psTree, "PansharpeningOptions");
287 89 : if (psOptions == nullptr)
288 : {
289 1 : CPLError(CE_Failure, CPLE_AppDefined, "Missing PansharpeningOptions");
290 1 : return CE_Failure;
291 : }
292 :
293 176 : if (CPLGetXMLValue(psOptions, "MSShiftX", nullptr) != nullptr ||
294 88 : CPLGetXMLValue(psOptions, "MSShiftY", nullptr) != nullptr)
295 : {
296 0 : CPLError(CE_Failure, CPLE_AppDefined,
297 : "Options MSShiftX and MSShiftY are no longer supported. "
298 : "Spatial adjustment between panchromatic and multispectral "
299 : "bands is done through their geotransform.");
300 0 : return CE_Failure;
301 : }
302 :
303 176 : CPLString osSourceFilename;
304 88 : GDALDataset *poPanDataset = nullptr;
305 88 : GDALRasterBand *poPanBand = nullptr;
306 176 : std::map<CPLString, GDALDataset *> oMapNamesToDataset;
307 : int nPanBand;
308 :
309 88 : if (hPanchroBandIn == nullptr)
310 : {
311 : const CPLXMLNode *psPanchroBand =
312 73 : CPLGetXMLNode(psOptions, "PanchroBand");
313 73 : if (psPanchroBand == nullptr)
314 : {
315 1 : CPLError(CE_Failure, CPLE_AppDefined, "PanchroBand missing");
316 4 : return CE_Failure;
317 : }
318 :
319 72 : osSourceFilename = CPLGetXMLValue(psPanchroBand, "SourceFilename", "");
320 72 : if (osSourceFilename.empty())
321 : {
322 1 : CPLError(CE_Failure, CPLE_AppDefined,
323 : "PanchroBand.SourceFilename missing");
324 1 : return CE_Failure;
325 : }
326 71 : const bool bRelativeToVRT = CPL_TO_BOOL(atoi(CPLGetXMLValue(
327 : psPanchroBand, "SourceFilename.relativetoVRT", "0")));
328 71 : if (bRelativeToVRT)
329 : {
330 : const std::string osAbs = CPLProjectRelativeFilenameSafe(
331 41 : pszVRTPathIn, osSourceFilename.c_str());
332 41 : m_oMapToRelativeFilenames[osAbs] = osSourceFilename;
333 41 : osSourceFilename = osAbs;
334 : }
335 :
336 : const CPLStringList aosOpenOptions(
337 71 : GDALDeserializeOpenOptionsFromXML(psPanchroBand));
338 :
339 71 : poPanDataset = GDALDataset::Open(osSourceFilename,
340 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
341 : nullptr, aosOpenOptions.List());
342 71 : if (poPanDataset == nullptr)
343 : {
344 1 : return CE_Failure;
345 : }
346 :
347 : const char *pszSourceBand =
348 70 : CPLGetXMLValue(psPanchroBand, "SourceBand", "1");
349 70 : nPanBand = atoi(pszSourceBand);
350 70 : if (poPanBand == nullptr)
351 70 : poPanBand = poPanDataset->GetRasterBand(nPanBand);
352 70 : if (poPanBand == nullptr)
353 : {
354 1 : CPLError(CE_Failure, CPLE_AppDefined, "%s invalid band of %s",
355 : pszSourceBand, osSourceFilename.c_str());
356 1 : GDALClose(poPanDataset);
357 1 : return CE_Failure;
358 : }
359 69 : oMapNamesToDataset[osSourceFilename] = poPanDataset;
360 : }
361 : else
362 : {
363 15 : poPanBand = GDALRasterBand::FromHandle(hPanchroBandIn);
364 15 : nPanBand = poPanBand->GetBand();
365 15 : poPanDataset = poPanBand->GetDataset();
366 15 : if (poPanDataset == nullptr)
367 : {
368 0 : CPLError(
369 : CE_Failure, CPLE_AppDefined,
370 : "Cannot retrieve dataset associated with panchromatic band");
371 0 : return CE_Failure;
372 : }
373 15 : oMapNamesToDataset[CPLSPrintf("%p", poPanDataset)] = poPanDataset;
374 15 : poPanDataset->Reference();
375 : }
376 84 : m_apoDatasetsToReleaseRef.push_back(
377 168 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>(
378 : poPanDataset));
379 84 : poPanDataset = m_apoDatasetsToReleaseRef.back().get();
380 : // Figure out which kind of adjustment we should do if the pan and spectral
381 : // bands do not share the same geotransform
382 : const char *pszGTAdjustment =
383 84 : CPLGetXMLValue(psOptions, "SpatialExtentAdjustment", "Union");
384 84 : if (EQUAL(pszGTAdjustment, "Union"))
385 78 : m_eGTAdjustment = GTAdjust_Union;
386 6 : else if (EQUAL(pszGTAdjustment, "Intersection"))
387 2 : m_eGTAdjustment = GTAdjust_Intersection;
388 4 : else if (EQUAL(pszGTAdjustment, "None"))
389 2 : m_eGTAdjustment = GTAdjust_None;
390 2 : else if (EQUAL(pszGTAdjustment, "NoneWithoutWarning"))
391 1 : m_eGTAdjustment = GTAdjust_NoneWithoutWarning;
392 : else
393 : {
394 1 : m_eGTAdjustment = GTAdjust_Union;
395 1 : CPLError(CE_Failure, CPLE_AppDefined,
396 : "Unsupported value for GeoTransformAdjustment. Defaulting to "
397 : "Union");
398 : }
399 :
400 : const char *pszNumThreads =
401 84 : CPLGetXMLValue(psOptions, "NumThreads", nullptr);
402 84 : int nThreads = 0;
403 84 : if (pszNumThreads != nullptr)
404 : {
405 37 : if (EQUAL(pszNumThreads, "ALL_CPUS"))
406 37 : nThreads = -1;
407 : else
408 0 : nThreads = atoi(pszNumThreads);
409 : }
410 :
411 : const char *pszAlgorithm =
412 84 : CPLGetXMLValue(psOptions, "Algorithm", "WeightedBrovey");
413 84 : if (!EQUAL(pszAlgorithm, "WeightedBrovey"))
414 : {
415 1 : CPLError(CE_Failure, CPLE_AppDefined, "Algorithm %s unsupported",
416 : pszAlgorithm);
417 1 : return CE_Failure;
418 : }
419 :
420 166 : std::vector<double> adfWeights;
421 83 : if (const CPLXMLNode *psAlgOptions =
422 83 : CPLGetXMLNode(psOptions, "AlgorithmOptions"))
423 : {
424 : const char *pszWeights =
425 35 : CPLGetXMLValue(psAlgOptions, "Weights", nullptr);
426 35 : if (pszWeights != nullptr)
427 : {
428 35 : char **papszTokens = CSLTokenizeString2(pszWeights, " ,", 0);
429 124 : for (int i = 0; papszTokens && papszTokens[i]; i++)
430 89 : adfWeights.push_back(CPLAtof(papszTokens[i]));
431 35 : CSLDestroy(papszTokens);
432 : }
433 : }
434 :
435 83 : GDALRIOResampleAlg eResampleAlg = GDALRasterIOGetResampleAlg(
436 : CPLGetXMLValue(psOptions, "Resampling", "Cubic"));
437 :
438 166 : std::vector<GDALRasterBand *> ahSpectralBands;
439 166 : std::map<int, int> aMapDstBandToSpectralBand;
440 83 : std::map<int, int>::iterator aMapDstBandToSpectralBandIter;
441 83 : int nBitDepth = 0;
442 83 : bool bFoundNonMatchingGT = false;
443 83 : double adfPanGT[6] = {0, 0, 0, 0, 0, 0};
444 : const bool bPanGeoTransformValid =
445 83 : (poPanDataset->GetGeoTransform(adfPanGT) == CE_None);
446 83 : if (!bPanGeoTransformValid)
447 : {
448 1 : CPLError(CE_Failure, CPLE_AppDefined,
449 : "Panchromatic band has no associated geotransform");
450 1 : return CE_Failure;
451 : }
452 82 : int nPanXSize = poPanBand->GetXSize();
453 82 : int nPanYSize = poPanBand->GetYSize();
454 82 : int bHasNoData = FALSE;
455 82 : double dfNoData = poPanBand->GetNoDataValue(&bHasNoData);
456 82 : double dfLRPanX =
457 82 : adfPanGT[0] + nPanXSize * adfPanGT[1] + nPanYSize * adfPanGT[2];
458 82 : double dfLRPanY =
459 82 : adfPanGT[3] + nPanXSize * adfPanGT[4] + nPanYSize * adfPanGT[5];
460 82 : bool bFoundRotatingTerms = (adfPanGT[2] != 0.0 || adfPanGT[4] != 0.0);
461 82 : double dfMinX = adfPanGT[0];
462 82 : double dfMaxX = dfLRPanX;
463 82 : double dfMaxY = adfPanGT[3];
464 82 : double dfMinY = dfLRPanY;
465 82 : if (dfMinY > dfMaxY)
466 11 : std::swap(dfMinY, dfMaxY);
467 82 : const double dfPanMinY = dfMinY;
468 82 : const double dfPanMaxY = dfMaxY;
469 :
470 82 : const auto poPanSRS = poPanDataset->GetSpatialRef();
471 :
472 : /* -------------------------------------------------------------------- */
473 : /* First pass on spectral datasets to check their georeferencing. */
474 : /* -------------------------------------------------------------------- */
475 82 : int iSpectralBand = 0;
476 528 : for (const CPLXMLNode *psIter = psOptions->psChild; psIter;
477 446 : psIter = psIter->psNext)
478 : {
479 : GDALDataset *poDataset;
480 450 : if (psIter->eType != CXT_Element ||
481 450 : !EQUAL(psIter->pszValue, "SpectralBand"))
482 246 : continue;
483 :
484 204 : if (nInputSpectralBandsIn && pahInputSpectralBandsIn != nullptr)
485 : {
486 40 : if (iSpectralBand == nInputSpectralBandsIn)
487 : {
488 1 : CPLError(CE_Failure, CPLE_AppDefined,
489 : "More SpectralBand elements than in source array");
490 4 : goto error;
491 : }
492 : poDataset = GDALRasterBand::FromHandle(
493 39 : pahInputSpectralBandsIn[iSpectralBand])
494 39 : ->GetDataset();
495 39 : if (poDataset == nullptr)
496 : {
497 0 : CPLError(CE_Failure, CPLE_AppDefined,
498 : "SpectralBand has no associated dataset");
499 0 : goto error;
500 : }
501 39 : osSourceFilename = poDataset->GetDescription();
502 :
503 39 : oMapNamesToDataset[CPLSPrintf("%p", poDataset)] = poDataset;
504 39 : poDataset->Reference();
505 : }
506 : else
507 : {
508 164 : osSourceFilename = CPLGetXMLValue(psIter, "SourceFilename", "");
509 164 : if (osSourceFilename.empty())
510 : {
511 1 : CPLError(CE_Failure, CPLE_AppDefined,
512 : "SpectralBand.SourceFilename missing");
513 1 : goto error;
514 : }
515 163 : int bRelativeToVRT = atoi(
516 : CPLGetXMLValue(psIter, "SourceFilename.relativetoVRT", "0"));
517 163 : if (bRelativeToVRT)
518 : {
519 : const std::string osAbs = CPLProjectRelativeFilenameSafe(
520 107 : pszVRTPathIn, osSourceFilename.c_str());
521 107 : m_oMapToRelativeFilenames[osAbs] = osSourceFilename;
522 107 : osSourceFilename = osAbs;
523 : }
524 163 : poDataset = oMapNamesToDataset[osSourceFilename];
525 163 : if (poDataset == nullptr)
526 : {
527 : const CPLStringList aosOpenOptions(
528 68 : GDALDeserializeOpenOptionsFromXML(psIter));
529 :
530 68 : poDataset = GDALDataset::Open(
531 : osSourceFilename, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
532 : nullptr, aosOpenOptions.List());
533 68 : if (poDataset == nullptr)
534 : {
535 1 : goto error;
536 : }
537 67 : oMapNamesToDataset[osSourceFilename] = poDataset;
538 : }
539 : else
540 : {
541 95 : poDataset->Reference();
542 : }
543 : }
544 201 : m_apoDatasetsToReleaseRef.push_back(
545 402 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>(
546 : poDataset));
547 201 : poDataset = m_apoDatasetsToReleaseRef.back().get();
548 :
549 : // Check that the spectral band has a georeferencing consistent
550 : // of the pan band. Allow an error of at most the size of one pixel
551 : // of the spectral band.
552 201 : const auto poSpectralSRS = poDataset->GetSpatialRef();
553 201 : if (poPanSRS)
554 : {
555 164 : if (poSpectralSRS)
556 : {
557 164 : if (!poPanSRS->IsSame(poSpectralSRS))
558 : {
559 1 : CPLError(CE_Warning, CPLE_AppDefined,
560 : "Pan dataset and %s do not seem to "
561 : "have same projection. Results might "
562 : "be incorrect",
563 : osSourceFilename.c_str());
564 : }
565 : }
566 : else
567 : {
568 0 : CPLError(CE_Warning, CPLE_AppDefined,
569 : "Pan dataset has a projection, whereas %s "
570 : "not. Results might be incorrect",
571 : osSourceFilename.c_str());
572 : }
573 : }
574 37 : else if (poSpectralSRS)
575 : {
576 0 : CPLError(CE_Warning, CPLE_AppDefined,
577 : "Pan dataset has no projection, whereas %s has "
578 : "one. Results might be incorrect",
579 : osSourceFilename.c_str());
580 : }
581 :
582 : double adfSpectralGeoTransform[6];
583 201 : if (poDataset->GetGeoTransform(adfSpectralGeoTransform) != CE_None)
584 : {
585 0 : CPLError(CE_Failure, CPLE_AppDefined,
586 : "Spectral band has no associated geotransform");
587 0 : goto error;
588 : }
589 201 : if (adfSpectralGeoTransform[3] * adfPanGT[3] < 0)
590 : {
591 1 : CPLError(CE_Failure, CPLE_NotSupported,
592 : "Spectral band vertical orientation is "
593 : "different from pan one");
594 1 : goto error;
595 : }
596 :
597 200 : int bIsThisOneNonMatching = FALSE;
598 : double dfPixelSize = std::max(adfSpectralGeoTransform[1],
599 200 : std::abs(adfSpectralGeoTransform[5]));
600 376 : if (std::abs(adfPanGT[0] - adfSpectralGeoTransform[0]) > dfPixelSize ||
601 176 : std::abs(adfPanGT[3] - adfSpectralGeoTransform[3]) > dfPixelSize)
602 : {
603 24 : bIsThisOneNonMatching = TRUE;
604 24 : if (m_eGTAdjustment == GTAdjust_None)
605 : {
606 2 : CPLError(CE_Warning, CPLE_AppDefined,
607 : "Georeferencing of top-left corner of pan "
608 : "dataset and %s do not match",
609 : osSourceFilename.c_str());
610 : }
611 : }
612 200 : bFoundRotatingTerms =
613 400 : bFoundRotatingTerms || (adfSpectralGeoTransform[2] != 0.0 ||
614 200 : adfSpectralGeoTransform[4] != 0.0);
615 : double dfLRSpectralX =
616 400 : adfSpectralGeoTransform[0] +
617 200 : poDataset->GetRasterXSize() * adfSpectralGeoTransform[1] +
618 200 : poDataset->GetRasterYSize() * adfSpectralGeoTransform[2];
619 : double dfLRSpectralY =
620 400 : adfSpectralGeoTransform[3] +
621 200 : poDataset->GetRasterXSize() * adfSpectralGeoTransform[4] +
622 200 : poDataset->GetRasterYSize() * adfSpectralGeoTransform[5];
623 376 : if (std::abs(dfLRPanX - dfLRSpectralX) > dfPixelSize ||
624 176 : std::abs(dfLRPanY - dfLRSpectralY) > dfPixelSize)
625 : {
626 24 : bIsThisOneNonMatching = TRUE;
627 24 : if (m_eGTAdjustment == GTAdjust_None)
628 : {
629 2 : CPLError(CE_Warning, CPLE_AppDefined,
630 : "Georeferencing of bottom-right corner of "
631 : "pan dataset and %s do not match",
632 : osSourceFilename.c_str());
633 : }
634 : }
635 :
636 200 : double dfThisMinY = dfLRSpectralY;
637 200 : double dfThisMaxY = adfSpectralGeoTransform[3];
638 200 : if (dfThisMinY > dfThisMaxY)
639 23 : std::swap(dfThisMinY, dfThisMaxY);
640 200 : if (bIsThisOneNonMatching && m_eGTAdjustment == GTAdjust_Union)
641 : {
642 20 : dfMinX = std::min(dfMinX, adfSpectralGeoTransform[0]);
643 20 : dfMinY = std::min(dfMinY, dfThisMinY);
644 20 : dfMaxX = std::max(dfMaxX, dfLRSpectralX);
645 20 : dfMaxY = std::max(dfMaxY, dfThisMaxY);
646 : }
647 180 : else if (bIsThisOneNonMatching &&
648 4 : m_eGTAdjustment == GTAdjust_Intersection)
649 : {
650 1 : dfMinX = std::max(dfMinX, adfSpectralGeoTransform[0]);
651 1 : dfMinY = std::max(dfMinY, dfThisMinY);
652 1 : dfMaxX = std::min(dfMaxX, dfLRSpectralX);
653 1 : dfMaxY = std::min(dfMaxY, dfThisMaxY);
654 : }
655 :
656 200 : bFoundNonMatchingGT = bFoundNonMatchingGT || bIsThisOneNonMatching;
657 :
658 200 : iSpectralBand++;
659 : }
660 :
661 78 : if (nInputSpectralBandsIn && iSpectralBand != nInputSpectralBandsIn)
662 : {
663 1 : CPLError(CE_Failure, CPLE_AppDefined,
664 : "Less SpectralBand elements than in source array");
665 1 : goto error;
666 : }
667 :
668 : /* -------------------------------------------------------------------- */
669 : /* On-the-fly spatial extent adjustment if needed and asked. */
670 : /* -------------------------------------------------------------------- */
671 :
672 77 : if (bFoundNonMatchingGT && (m_eGTAdjustment == GTAdjust_Union ||
673 4 : m_eGTAdjustment == GTAdjust_Intersection))
674 : {
675 8 : if (bFoundRotatingTerms)
676 : {
677 0 : CPLError(
678 : CE_Failure, CPLE_NotSupported,
679 : "One of the panchromatic or spectral datasets has rotating "
680 : "terms in their geotransform matrix. Adjustment not possible");
681 0 : goto error;
682 : }
683 8 : if (m_eGTAdjustment == GTAdjust_Intersection &&
684 1 : (dfMinX >= dfMaxX || dfMinY >= dfMaxY))
685 : {
686 0 : CPLError(
687 : CE_Failure, CPLE_NotSupported,
688 : "One of the panchromatic or spectral datasets has rotating "
689 : "terms in their geotransform matrix. Adjustment not possible");
690 0 : goto error;
691 : }
692 8 : if (m_eGTAdjustment == GTAdjust_Union)
693 7 : CPLDebug("VRT", "Do union of bounding box of panchromatic and "
694 : "spectral datasets");
695 : else
696 1 : CPLDebug("VRT", "Do intersection of bounding box of panchromatic "
697 : "and spectral datasets");
698 :
699 : // If the pandataset needs adjustments, make sure the coordinates of the
700 : // union/intersection properly align with the grid of the pandataset
701 : // to avoid annoying sub-pixel shifts on the panchro band.
702 8 : double dfPixelSize = std::max(adfPanGT[1], std::abs(adfPanGT[5]));
703 8 : if (std::abs(adfPanGT[0] - dfMinX) > dfPixelSize ||
704 3 : std::abs(dfPanMaxY - dfMaxY) > dfPixelSize ||
705 13 : std::abs(dfLRPanX - dfMaxX) > dfPixelSize ||
706 2 : std::abs(dfPanMinY - dfMinY) > dfPixelSize)
707 : {
708 6 : dfMinX = adfPanGT[0] +
709 6 : std::floor((dfMinX - adfPanGT[0]) / adfPanGT[1] + 0.5) *
710 6 : adfPanGT[1];
711 6 : dfMaxX =
712 6 : dfLRPanX + std::floor((dfMaxX - dfLRPanX) / adfPanGT[1] + 0.5) *
713 6 : adfPanGT[1];
714 6 : if (adfPanGT[5] > 0)
715 : {
716 0 : dfMinY = dfPanMinY +
717 0 : std::floor((dfMinY - dfPanMinY) / adfPanGT[5] + 0.5) *
718 0 : adfPanGT[5];
719 0 : dfMaxY = dfPanMinY +
720 0 : std::floor((dfMaxY - dfPanMinY) / adfPanGT[5] + 0.5) *
721 0 : adfPanGT[5];
722 : }
723 : else
724 : {
725 6 : dfMinY = dfPanMaxY +
726 6 : std::floor((dfMinY - dfPanMaxY) / adfPanGT[5] + 0.5) *
727 6 : adfPanGT[5];
728 6 : dfMaxY = dfPanMaxY +
729 6 : std::floor((dfMaxY - dfPanMaxY) / adfPanGT[5] + 0.5) *
730 6 : adfPanGT[5];
731 : }
732 : }
733 :
734 : std::map<CPLString, GDALDataset *>::iterator oIter =
735 8 : oMapNamesToDataset.begin();
736 24 : for (; oIter != oMapNamesToDataset.end(); ++oIter)
737 : {
738 16 : GDALDataset *poSrcDS = oIter->second;
739 : double adfGT[6];
740 16 : if (poSrcDS->GetGeoTransform(adfGT) != CE_None)
741 3 : continue;
742 :
743 : // Check if this dataset needs adjustments
744 16 : dfPixelSize = std::max(adfGT[1], std::abs(adfGT[5]));
745 16 : dfPixelSize = std::max(adfPanGT[1], dfPixelSize);
746 16 : dfPixelSize = std::max(std::abs(adfPanGT[5]), dfPixelSize);
747 16 : double dfThisMinX = adfGT[0];
748 16 : double dfThisMaxX = adfGT[0] + poSrcDS->GetRasterXSize() * adfGT[1];
749 16 : double dfThisMaxY = adfGT[3];
750 16 : double dfThisMinY = adfGT[3] + poSrcDS->GetRasterYSize() * adfGT[5];
751 16 : if (dfThisMinY > dfThisMaxY)
752 2 : std::swap(dfThisMinY, dfThisMaxY);
753 16 : if (std::abs(dfThisMinX - dfMinX) <= dfPixelSize &&
754 8 : std::abs(dfThisMaxY - dfMaxY) <= dfPixelSize &&
755 27 : std::abs(dfThisMaxX - dfMaxX) <= dfPixelSize &&
756 3 : std::abs(dfThisMinY - dfMinY) <= dfPixelSize)
757 : {
758 3 : continue;
759 : }
760 :
761 : double adfAdjustedGT[6];
762 13 : adfAdjustedGT[0] = dfMinX;
763 13 : adfAdjustedGT[1] = adfGT[1];
764 13 : adfAdjustedGT[2] = 0;
765 13 : adfAdjustedGT[3] = adfGT[5] > 0 ? dfMinY : dfMaxY;
766 13 : adfAdjustedGT[4] = 0;
767 13 : adfAdjustedGT[5] = adfGT[5];
768 13 : int nAdjustRasterXSize =
769 13 : static_cast<int>(0.5 + (dfMaxX - dfMinX) / adfAdjustedGT[1]);
770 : int nAdjustRasterYSize = static_cast<int>(
771 13 : 0.5 + (dfMaxY - dfMinY) / std::abs(adfAdjustedGT[5]));
772 :
773 : VRTDataset *poVDS =
774 13 : new VRTDataset(nAdjustRasterXSize, nAdjustRasterYSize);
775 13 : poVDS->SetWritable(FALSE);
776 26 : poVDS->SetDescription(std::string("Shifted ")
777 13 : .append(poSrcDS->GetDescription())
778 13 : .c_str());
779 13 : poVDS->SetGeoTransform(adfAdjustedGT);
780 13 : poVDS->SetProjection(poPanDataset->GetProjectionRef());
781 :
782 32 : for (int i = 0; i < poSrcDS->GetRasterCount(); i++)
783 : {
784 19 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(i + 1);
785 19 : poVDS->AddBand(poSrcBand->GetRasterDataType(), nullptr);
786 : VRTSourcedRasterBand *poVRTBand =
787 : static_cast<VRTSourcedRasterBand *>(
788 19 : poVDS->GetRasterBand(i + 1));
789 :
790 : const char *pszNBITS =
791 19 : poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
792 19 : if (pszNBITS)
793 0 : poVRTBand->SetMetadataItem("NBITS", pszNBITS,
794 0 : "IMAGE_STRUCTURE");
795 :
796 19 : VRTSimpleSource *poSimpleSource = new VRTSimpleSource();
797 76 : poVRTBand->ConfigureSource(
798 : poSimpleSource, poSrcBand, FALSE,
799 19 : static_cast<int>(
800 19 : std::floor((dfMinX - adfGT[0]) / adfGT[1] + 0.001)),
801 19 : adfGT[5] < 0
802 16 : ? static_cast<int>(std::floor(
803 16 : (dfMaxY - dfThisMaxY) / adfGT[5] + 0.001))
804 3 : : static_cast<int>(std::floor(
805 3 : (dfMinY - dfThisMinY) / adfGT[5] + 0.001)),
806 19 : static_cast<int>(0.5 + (dfMaxX - dfMinX) / adfGT[1]),
807 19 : static_cast<int>(0.5 +
808 19 : (dfMaxY - dfMinY) / std::abs(adfGT[5])),
809 : 0, 0, nAdjustRasterXSize, nAdjustRasterYSize);
810 :
811 19 : poVRTBand->AddSource(poSimpleSource);
812 : }
813 :
814 13 : oIter->second = poVDS;
815 13 : if (poSrcDS == poPanDataset)
816 : {
817 6 : memcpy(adfPanGT, adfAdjustedGT, 6 * sizeof(double));
818 6 : poPanDataset = poVDS;
819 6 : poPanBand = poPanDataset->GetRasterBand(nPanBand);
820 6 : nPanXSize = poPanDataset->GetRasterXSize();
821 6 : nPanYSize = poPanDataset->GetRasterYSize();
822 : }
823 :
824 13 : m_apoDatasetsToReleaseRef.push_back(
825 26 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>(
826 : poVDS));
827 : }
828 : }
829 :
830 77 : if (nRasterXSize == 0 && nRasterYSize == 0)
831 : {
832 57 : nRasterXSize = nPanXSize;
833 57 : nRasterYSize = nPanYSize;
834 : }
835 20 : else if (nRasterXSize != nPanXSize || nRasterYSize != nPanYSize)
836 : {
837 1 : CPLError(CE_Failure, CPLE_AppDefined,
838 : "Inconsistent declared VRT dimensions with panchro dataset");
839 1 : goto error;
840 : }
841 :
842 : /* -------------------------------------------------------------------- */
843 : /* Initialize all the general VRT stuff. This will even */
844 : /* create the VRTPansharpenedRasterBands and initialize them. */
845 : /* -------------------------------------------------------------------- */
846 76 : eErr = VRTDataset::XMLInit(psTree, pszVRTPathIn);
847 :
848 76 : if (eErr != CE_None)
849 : {
850 1 : goto error;
851 : }
852 :
853 : /* -------------------------------------------------------------------- */
854 : /* Inherit georeferencing info from panchro band if not defined */
855 : /* in VRT. */
856 : /* -------------------------------------------------------------------- */
857 :
858 : {
859 : double adfOutGT[6];
860 140 : if (GetGeoTransform(adfOutGT) != CE_None && GetGCPCount() == 0 &&
861 65 : GetSpatialRef() == nullptr)
862 : {
863 65 : if (bPanGeoTransformValid)
864 : {
865 65 : SetGeoTransform(adfPanGT);
866 : }
867 65 : if (poPanSRS)
868 : {
869 47 : SetSpatialRef(poPanSRS);
870 : }
871 : }
872 : }
873 :
874 : /* -------------------------------------------------------------------- */
875 : /* Parse rest of PansharpeningOptions */
876 : /* -------------------------------------------------------------------- */
877 75 : iSpectralBand = 0;
878 476 : for (const CPLXMLNode *psIter = psOptions->psChild; psIter;
879 401 : psIter = psIter->psNext)
880 : {
881 404 : if (psIter->eType != CXT_Element ||
882 404 : !EQUAL(psIter->pszValue, "SpectralBand"))
883 220 : continue;
884 :
885 : GDALDataset *poDataset;
886 : GDALRasterBand *poBand;
887 :
888 184 : const char *pszDstBand = CPLGetXMLValue(psIter, "dstBand", nullptr);
889 184 : int nDstBand = -1;
890 184 : if (pszDstBand != nullptr)
891 : {
892 178 : nDstBand = atoi(pszDstBand);
893 178 : if (nDstBand <= 0)
894 : {
895 1 : CPLError(CE_Failure, CPLE_AppDefined,
896 : "SpectralBand.dstBand = '%s' invalid", pszDstBand);
897 3 : goto error;
898 : }
899 : }
900 :
901 183 : if (nInputSpectralBandsIn && pahInputSpectralBandsIn != nullptr)
902 : {
903 68 : poBand = GDALRasterBand::FromHandle(
904 34 : pahInputSpectralBandsIn[iSpectralBand]);
905 34 : poDataset = poBand->GetDataset();
906 34 : if (poDataset)
907 : {
908 34 : poDataset = oMapNamesToDataset[CPLSPrintf("%p", poDataset)];
909 34 : CPLAssert(poDataset);
910 34 : if (poDataset)
911 34 : poBand = poDataset->GetRasterBand(poBand->GetBand());
912 : }
913 : }
914 : else
915 : {
916 149 : osSourceFilename = CPLGetXMLValue(psIter, "SourceFilename", "");
917 149 : const bool bRelativeToVRT = CPL_TO_BOOL(atoi(
918 : CPLGetXMLValue(psIter, "SourceFilename.relativetoVRT", "0")));
919 149 : if (bRelativeToVRT)
920 188 : osSourceFilename = CPLProjectRelativeFilenameSafe(
921 94 : pszVRTPathIn, osSourceFilename.c_str());
922 149 : poDataset = oMapNamesToDataset[osSourceFilename];
923 149 : CPLAssert(poDataset);
924 : const char *pszSourceBand =
925 149 : CPLGetXMLValue(psIter, "SourceBand", "1");
926 149 : const int nBand = atoi(pszSourceBand);
927 149 : poBand = poDataset->GetRasterBand(nBand);
928 149 : if (poBand == nullptr)
929 : {
930 1 : CPLError(CE_Failure, CPLE_AppDefined, "%s invalid band of %s",
931 : pszSourceBand, osSourceFilename.c_str());
932 1 : goto error;
933 : }
934 : }
935 :
936 182 : if (bHasNoData)
937 : {
938 3 : double dfSpectralNoData = poPanBand->GetNoDataValue(&bHasNoData);
939 3 : if (bHasNoData && dfSpectralNoData != dfNoData)
940 0 : bHasNoData = FALSE;
941 : }
942 :
943 182 : ahSpectralBands.push_back(poBand);
944 182 : if (nDstBand >= 1)
945 : {
946 176 : if (aMapDstBandToSpectralBand.find(nDstBand - 1) !=
947 352 : aMapDstBandToSpectralBand.end())
948 : {
949 1 : CPLError(
950 : CE_Failure, CPLE_AppDefined,
951 : "Another spectral band is already mapped to output band %d",
952 : nDstBand);
953 1 : goto error;
954 : }
955 175 : aMapDstBandToSpectralBand[nDstBand - 1] =
956 175 : static_cast<int>(ahSpectralBands.size() - 1);
957 : }
958 :
959 181 : iSpectralBand++;
960 : }
961 :
962 72 : if (ahSpectralBands.empty())
963 : {
964 1 : CPLError(CE_Failure, CPLE_AppDefined, "No spectral band defined");
965 1 : goto error;
966 : }
967 :
968 : {
969 71 : const char *pszNoData = CPLGetXMLValue(psOptions, "NoData", nullptr);
970 71 : if (pszNoData)
971 : {
972 7 : if (EQUAL(pszNoData, "NONE"))
973 : {
974 0 : m_bNoDataDisabled = TRUE;
975 0 : bHasNoData = FALSE;
976 : }
977 : else
978 : {
979 7 : bHasNoData = TRUE;
980 7 : dfNoData = CPLAtof(pszNoData);
981 : }
982 : }
983 : }
984 :
985 71 : if (GetRasterCount() == 0)
986 : {
987 156 : for (int i = 0; i < static_cast<int>(aMapDstBandToSpectralBand.size());
988 : i++)
989 : {
990 109 : if (aMapDstBandToSpectralBand.find(i) ==
991 218 : aMapDstBandToSpectralBand.end())
992 : {
993 1 : CPLError(CE_Failure, CPLE_AppDefined,
994 : "Hole in SpectralBand.dstBand numbering");
995 1 : goto error;
996 : }
997 108 : GDALRasterBand *poInputBand = GDALRasterBand::FromHandle(
998 108 : ahSpectralBands[aMapDstBandToSpectralBand[i]]);
999 : GDALRasterBand *poBand = new VRTPansharpenedRasterBand(
1000 108 : this, i + 1, poInputBand->GetRasterDataType());
1001 108 : poBand->SetColorInterpretation(
1002 108 : poInputBand->GetColorInterpretation());
1003 108 : if (bHasNoData)
1004 13 : poBand->SetNoDataValue(dfNoData);
1005 108 : SetBand(i + 1, poBand);
1006 : }
1007 : }
1008 : else
1009 : {
1010 23 : int nIdxAsPansharpenedBand = 0;
1011 91 : for (int i = 0; i < nBands; i++)
1012 : {
1013 69 : if (static_cast<VRTRasterBand *>(GetRasterBand(i + 1))
1014 69 : ->IsPansharpenRasterBand())
1015 : {
1016 65 : if (aMapDstBandToSpectralBand.find(i) ==
1017 130 : aMapDstBandToSpectralBand.end())
1018 : {
1019 1 : CPLError(CE_Failure, CPLE_AppDefined,
1020 : "Band %d of type VRTPansharpenedRasterBand, but "
1021 : "no corresponding SpectralBand",
1022 : i + 1);
1023 1 : goto error;
1024 : }
1025 : else
1026 : {
1027 : static_cast<VRTPansharpenedRasterBand *>(
1028 64 : GetRasterBand(i + 1))
1029 64 : ->SetIndexAsPansharpenedBand(nIdxAsPansharpenedBand);
1030 64 : nIdxAsPansharpenedBand++;
1031 : }
1032 : }
1033 : }
1034 : }
1035 :
1036 : // Figure out bit depth
1037 : {
1038 : const char *pszBitDepth =
1039 69 : CPLGetXMLValue(psOptions, "BitDepth", nullptr);
1040 69 : if (pszBitDepth == nullptr)
1041 56 : pszBitDepth = GDALRasterBand::FromHandle(ahSpectralBands[0])
1042 56 : ->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1043 69 : if (pszBitDepth)
1044 15 : nBitDepth = atoi(pszBitDepth);
1045 69 : if (nBitDepth)
1046 : {
1047 44 : for (int i = 0; i < nBands; i++)
1048 : {
1049 29 : if (!static_cast<VRTRasterBand *>(GetRasterBand(i + 1))
1050 29 : ->IsPansharpenRasterBand())
1051 2 : continue;
1052 27 : if (GetRasterBand(i + 1)->GetMetadataItem(
1053 27 : "NBITS", "IMAGE_STRUCTURE") == nullptr)
1054 : {
1055 27 : if (nBitDepth != 8 && nBitDepth != 16 && nBitDepth != 32)
1056 : {
1057 12 : GetRasterBand(i + 1)->SetMetadataItem(
1058 : "NBITS", CPLSPrintf("%d", nBitDepth),
1059 6 : "IMAGE_STRUCTURE");
1060 : }
1061 : }
1062 0 : else if (nBitDepth == 8 || nBitDepth == 16 || nBitDepth == 32)
1063 : {
1064 0 : GetRasterBand(i + 1)->SetMetadataItem("NBITS", nullptr,
1065 0 : "IMAGE_STRUCTURE");
1066 : }
1067 : }
1068 : }
1069 : }
1070 :
1071 69 : if (GDALGetRasterBandXSize(ahSpectralBands[0]) >
1072 138 : GDALGetRasterBandXSize(poPanBand) ||
1073 69 : GDALGetRasterBandYSize(ahSpectralBands[0]) >
1074 69 : GDALGetRasterBandYSize(poPanBand))
1075 : {
1076 0 : CPLError(CE_Warning, CPLE_AppDefined,
1077 : "Dimensions of spectral band larger than panchro band");
1078 : }
1079 :
1080 69 : aMapDstBandToSpectralBandIter = aMapDstBandToSpectralBand.begin();
1081 236 : for (; aMapDstBandToSpectralBandIter != aMapDstBandToSpectralBand.end();
1082 167 : ++aMapDstBandToSpectralBandIter)
1083 : {
1084 168 : const int nDstBand = 1 + aMapDstBandToSpectralBandIter->first;
1085 335 : if (nDstBand > nBands ||
1086 167 : !static_cast<VRTRasterBand *>(GetRasterBand(nDstBand))
1087 167 : ->IsPansharpenRasterBand())
1088 : {
1089 1 : CPLError(CE_Failure, CPLE_AppDefined,
1090 : "SpectralBand.dstBand = '%d' invalid", nDstBand);
1091 1 : goto error;
1092 : }
1093 : }
1094 :
1095 68 : if (adfWeights.empty())
1096 : {
1097 158 : for (int i = 0; i < static_cast<int>(ahSpectralBands.size()); i++)
1098 : {
1099 114 : adfWeights.push_back(1.0 / ahSpectralBands.size());
1100 : }
1101 : }
1102 24 : else if (adfWeights.size() != ahSpectralBands.size())
1103 : {
1104 1 : CPLError(CE_Failure, CPLE_AppDefined,
1105 : "%d weights defined, but %d input spectral bands",
1106 1 : static_cast<int>(adfWeights.size()),
1107 1 : static_cast<int>(ahSpectralBands.size()));
1108 1 : goto error;
1109 : }
1110 :
1111 67 : if (aMapDstBandToSpectralBand.empty())
1112 : {
1113 1 : CPLError(CE_Warning, CPLE_AppDefined,
1114 : "No spectral band is mapped to an output band");
1115 : }
1116 :
1117 : /* -------------------------------------------------------------------- */
1118 : /* Instantiate poPansharpener */
1119 : /* -------------------------------------------------------------------- */
1120 67 : psPanOptions = GDALCreatePansharpenOptions();
1121 67 : psPanOptions->ePansharpenAlg = GDAL_PSH_WEIGHTED_BROVEY;
1122 67 : psPanOptions->eResampleAlg = eResampleAlg;
1123 67 : psPanOptions->nBitDepth = nBitDepth;
1124 67 : psPanOptions->nWeightCount = static_cast<int>(adfWeights.size());
1125 67 : psPanOptions->padfWeights =
1126 67 : static_cast<double *>(CPLMalloc(sizeof(double) * adfWeights.size()));
1127 67 : memcpy(psPanOptions->padfWeights, &adfWeights[0],
1128 67 : sizeof(double) * adfWeights.size());
1129 67 : psPanOptions->hPanchroBand = poPanBand;
1130 67 : psPanOptions->nInputSpectralBands =
1131 67 : static_cast<int>(ahSpectralBands.size());
1132 67 : psPanOptions->pahInputSpectralBands = static_cast<GDALRasterBandH *>(
1133 67 : CPLMalloc(sizeof(GDALRasterBandH) * ahSpectralBands.size()));
1134 67 : memcpy(psPanOptions->pahInputSpectralBands, &ahSpectralBands[0],
1135 67 : sizeof(GDALRasterBandH) * ahSpectralBands.size());
1136 67 : psPanOptions->nOutPansharpenedBands =
1137 67 : static_cast<int>(aMapDstBandToSpectralBand.size());
1138 67 : psPanOptions->panOutPansharpenedBands = static_cast<int *>(
1139 67 : CPLMalloc(sizeof(int) * aMapDstBandToSpectralBand.size()));
1140 67 : aMapDstBandToSpectralBandIter = aMapDstBandToSpectralBand.begin();
1141 229 : for (int i = 0;
1142 229 : aMapDstBandToSpectralBandIter != aMapDstBandToSpectralBand.end();
1143 162 : ++aMapDstBandToSpectralBandIter, ++i)
1144 : {
1145 162 : psPanOptions->panOutPansharpenedBands[i] =
1146 162 : aMapDstBandToSpectralBandIter->second;
1147 : }
1148 67 : psPanOptions->bHasNoData = bHasNoData;
1149 67 : psPanOptions->dfNoData = dfNoData;
1150 67 : psPanOptions->nThreads = nThreads;
1151 :
1152 67 : if (nBands == psPanOptions->nOutPansharpenedBands)
1153 63 : SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
1154 :
1155 67 : m_poPansharpener = new GDALPansharpenOperation();
1156 67 : eErr = m_poPansharpener->Initialize(psPanOptions);
1157 67 : if (eErr != CE_None)
1158 : {
1159 : // Delete the pansharper object before closing the dataset
1160 : // because it may have warped the bands into an intermediate VRT
1161 1 : delete m_poPansharpener;
1162 1 : m_poPansharpener = nullptr;
1163 :
1164 : // Close in reverse order (VRT firsts and real datasets after)
1165 4 : for (int i = static_cast<int>(m_apoDatasetsToReleaseRef.size()) - 1;
1166 4 : i >= 0; i--)
1167 : {
1168 3 : m_apoDatasetsToReleaseRef[i].reset();
1169 : }
1170 1 : m_apoDatasetsToReleaseRef.clear();
1171 : }
1172 67 : GDALDestroyPansharpenOptions(psPanOptions);
1173 :
1174 67 : return eErr;
1175 :
1176 15 : error:
1177 : // Close in reverse order (VRT firsts and real datasets after)
1178 63 : for (int i = static_cast<int>(m_apoDatasetsToReleaseRef.size()) - 1; i >= 0;
1179 : i--)
1180 : {
1181 48 : m_apoDatasetsToReleaseRef[i].reset();
1182 : }
1183 15 : m_apoDatasetsToReleaseRef.clear();
1184 15 : return CE_Failure;
1185 : }
1186 :
1187 : /************************************************************************/
1188 : /* SerializeToXML() */
1189 : /************************************************************************/
1190 :
1191 7 : CPLXMLNode *VRTPansharpenedDataset::SerializeToXML(const char *pszVRTPathIn)
1192 :
1193 : {
1194 7 : CPLXMLNode *psTree = VRTDataset::SerializeToXML(pszVRTPathIn);
1195 :
1196 7 : if (psTree == nullptr)
1197 0 : return psTree;
1198 :
1199 : /* -------------------------------------------------------------------- */
1200 : /* Set subclass. */
1201 : /* -------------------------------------------------------------------- */
1202 7 : CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
1203 : CXT_Text, "VRTPansharpenedDataset");
1204 :
1205 : /* -------------------------------------------------------------------- */
1206 : /* Serialize the block size. */
1207 : /* -------------------------------------------------------------------- */
1208 7 : CPLCreateXMLElementAndValue(psTree, "BlockXSize",
1209 : CPLSPrintf("%d", m_nBlockXSize));
1210 7 : CPLCreateXMLElementAndValue(psTree, "BlockYSize",
1211 : CPLSPrintf("%d", m_nBlockYSize));
1212 :
1213 : /* -------------------------------------------------------------------- */
1214 : /* Serialize the options. */
1215 : /* -------------------------------------------------------------------- */
1216 7 : if (m_poPansharpener == nullptr)
1217 0 : return psTree;
1218 7 : GDALPansharpenOptions *psOptions = m_poPansharpener->GetOptions();
1219 7 : if (psOptions == nullptr)
1220 0 : return psTree;
1221 :
1222 : CPLXMLNode *psOptionsNode =
1223 7 : CPLCreateXMLNode(psTree, CXT_Element, "PansharpeningOptions");
1224 :
1225 7 : if (psOptions->ePansharpenAlg == GDAL_PSH_WEIGHTED_BROVEY)
1226 : {
1227 7 : CPLCreateXMLElementAndValue(psOptionsNode, "Algorithm",
1228 : "WeightedBrovey");
1229 : }
1230 : else
1231 : {
1232 0 : CPLAssert(false);
1233 : }
1234 7 : if (psOptions->nWeightCount)
1235 : {
1236 14 : CPLString osWeights;
1237 27 : for (int i = 0; i < psOptions->nWeightCount; i++)
1238 : {
1239 20 : if (i)
1240 13 : osWeights += ",";
1241 20 : osWeights += CPLSPrintf("%.16g", psOptions->padfWeights[i]);
1242 : }
1243 7 : CPLCreateXMLElementAndValue(
1244 : CPLCreateXMLNode(psOptionsNode, CXT_Element, "AlgorithmOptions"),
1245 : "Weights", osWeights.c_str());
1246 : }
1247 7 : CPLCreateXMLElementAndValue(
1248 : psOptionsNode, "Resampling",
1249 : GDALRasterIOGetResampleAlg(psOptions->eResampleAlg));
1250 :
1251 7 : if (psOptions->nThreads == -1)
1252 : {
1253 6 : CPLCreateXMLElementAndValue(psOptionsNode, "NumThreads", "ALL_CPUS");
1254 : }
1255 1 : else if (psOptions->nThreads > 1)
1256 : {
1257 0 : CPLCreateXMLElementAndValue(psOptionsNode, "NumThreads",
1258 : CPLSPrintf("%d", psOptions->nThreads));
1259 : }
1260 :
1261 7 : if (psOptions->nBitDepth)
1262 0 : CPLCreateXMLElementAndValue(psOptionsNode, "BitDepth",
1263 : CPLSPrintf("%d", psOptions->nBitDepth));
1264 :
1265 7 : const char *pszAdjust = nullptr;
1266 7 : switch (m_eGTAdjustment)
1267 : {
1268 7 : case GTAdjust_Union:
1269 7 : pszAdjust = "Union";
1270 7 : break;
1271 0 : case GTAdjust_Intersection:
1272 0 : pszAdjust = "Intersection";
1273 0 : break;
1274 0 : case GTAdjust_None:
1275 0 : pszAdjust = "None";
1276 0 : break;
1277 0 : case GTAdjust_NoneWithoutWarning:
1278 0 : pszAdjust = "NoneWithoutWarning";
1279 0 : break;
1280 0 : default:
1281 0 : break;
1282 : }
1283 :
1284 7 : if (psOptions->bHasNoData)
1285 : {
1286 1 : CPLCreateXMLElementAndValue(psOptionsNode, "NoData",
1287 : CPLSPrintf("%.16g", psOptions->dfNoData));
1288 : }
1289 6 : else if (m_bNoDataDisabled)
1290 : {
1291 0 : CPLCreateXMLElementAndValue(psOptionsNode, "NoData", "None");
1292 : }
1293 :
1294 7 : if (pszAdjust)
1295 7 : CPLCreateXMLElementAndValue(psOptionsNode, "SpatialExtentAdjustment",
1296 : pszAdjust);
1297 :
1298 7 : if (psOptions->hPanchroBand)
1299 : {
1300 : CPLXMLNode *psBand =
1301 7 : CPLCreateXMLNode(psOptionsNode, CXT_Element, "PanchroBand");
1302 : GDALRasterBand *poBand =
1303 7 : GDALRasterBand::FromHandle(psOptions->hPanchroBand);
1304 7 : if (poBand->GetDataset())
1305 : {
1306 7 : auto poPanchoDS = poBand->GetDataset();
1307 : std::map<CPLString, CPLString>::iterator oIter =
1308 7 : m_oMapToRelativeFilenames.find(poPanchoDS->GetDescription());
1309 7 : if (oIter == m_oMapToRelativeFilenames.end())
1310 : {
1311 6 : CPLCreateXMLElementAndValue(psBand, "SourceFilename",
1312 6 : poPanchoDS->GetDescription());
1313 : }
1314 : else
1315 : {
1316 1 : CPLXMLNode *psSourceFilename = CPLCreateXMLElementAndValue(
1317 1 : psBand, "SourceFilename", oIter->second);
1318 1 : CPLCreateXMLNode(CPLCreateXMLNode(psSourceFilename,
1319 : CXT_Attribute,
1320 : "relativeToVRT"),
1321 : CXT_Text, "1");
1322 : }
1323 :
1324 7 : GDALSerializeOpenOptionsToXML(psBand, poPanchoDS->GetOpenOptions());
1325 :
1326 7 : CPLCreateXMLElementAndValue(psBand, "SourceBand",
1327 : CPLSPrintf("%d", poBand->GetBand()));
1328 : }
1329 : }
1330 27 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
1331 : {
1332 : CPLXMLNode *psBand =
1333 20 : CPLCreateXMLNode(psOptionsNode, CXT_Element, "SpectralBand");
1334 :
1335 39 : for (int j = 0; j < psOptions->nOutPansharpenedBands; j++)
1336 : {
1337 39 : if (psOptions->panOutPansharpenedBands[j] == i)
1338 : {
1339 40 : for (int k = 0; k < nBands; k++)
1340 : {
1341 40 : if (static_cast<VRTRasterBand *>(GetRasterBand(k + 1))
1342 40 : ->IsPansharpenRasterBand())
1343 : {
1344 39 : if (static_cast<VRTPansharpenedRasterBand *>(
1345 39 : GetRasterBand(k + 1))
1346 39 : ->GetIndexAsPansharpenedBand() == j)
1347 : {
1348 20 : CPLCreateXMLNode(CPLCreateXMLNode(psBand,
1349 : CXT_Attribute,
1350 : "dstBand"),
1351 : CXT_Text, CPLSPrintf("%d", k + 1));
1352 20 : break;
1353 : }
1354 : }
1355 : }
1356 20 : break;
1357 : }
1358 : }
1359 :
1360 : GDALRasterBand *poBand =
1361 20 : GDALRasterBand::FromHandle(psOptions->pahInputSpectralBands[i]);
1362 20 : if (poBand->GetDataset())
1363 : {
1364 20 : auto poSpectralDS = poBand->GetDataset();
1365 : std::map<CPLString, CPLString>::iterator oIter =
1366 20 : m_oMapToRelativeFilenames.find(poSpectralDS->GetDescription());
1367 20 : if (oIter == m_oMapToRelativeFilenames.end())
1368 : {
1369 20 : CPLCreateXMLElementAndValue(psBand, "SourceFilename",
1370 20 : poSpectralDS->GetDescription());
1371 : }
1372 : else
1373 : {
1374 0 : CPLXMLNode *psSourceFilename = CPLCreateXMLElementAndValue(
1375 0 : psBand, "SourceFilename", oIter->second);
1376 0 : CPLCreateXMLNode(CPLCreateXMLNode(psSourceFilename,
1377 : CXT_Attribute,
1378 : "relativeToVRT"),
1379 : CXT_Text, "1");
1380 : }
1381 :
1382 20 : GDALSerializeOpenOptionsToXML(psBand,
1383 20 : poSpectralDS->GetOpenOptions());
1384 :
1385 20 : CPLCreateXMLElementAndValue(psBand, "SourceBand",
1386 : CPLSPrintf("%d", poBand->GetBand()));
1387 : }
1388 : }
1389 :
1390 7 : return psTree;
1391 : }
1392 :
1393 : /************************************************************************/
1394 : /* GetBlockSize() */
1395 : /************************************************************************/
1396 :
1397 192 : void VRTPansharpenedDataset::GetBlockSize(int *pnBlockXSize,
1398 : int *pnBlockYSize) const
1399 :
1400 : {
1401 192 : assert(nullptr != pnBlockXSize);
1402 192 : assert(nullptr != pnBlockYSize);
1403 :
1404 192 : *pnBlockXSize = m_nBlockXSize;
1405 192 : *pnBlockYSize = m_nBlockYSize;
1406 192 : }
1407 :
1408 : /************************************************************************/
1409 : /* AddBand() */
1410 : /************************************************************************/
1411 :
1412 1 : CPLErr VRTPansharpenedDataset::AddBand(CPL_UNUSED GDALDataType eType,
1413 : CPL_UNUSED char **papszOptions)
1414 :
1415 : {
1416 1 : CPLError(CE_Failure, CPLE_NotSupported, "AddBand() not supported");
1417 :
1418 1 : return CE_Failure;
1419 : }
1420 :
1421 : /************************************************************************/
1422 : /* IRasterIO() */
1423 : /************************************************************************/
1424 :
1425 21 : CPLErr VRTPansharpenedDataset::IRasterIO(
1426 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
1427 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
1428 : int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
1429 : GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
1430 : {
1431 21 : if (eRWFlag == GF_Write)
1432 0 : return CE_Failure;
1433 :
1434 : /* Try to pass the request to the most appropriate overview dataset */
1435 21 : if (nBufXSize < nXSize && nBufYSize < nYSize)
1436 : {
1437 : int bTried;
1438 1 : CPLErr eErr = TryOverviewRasterIO(
1439 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1440 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
1441 : nBandSpace, psExtraArg, &bTried);
1442 1 : if (bTried)
1443 0 : return eErr;
1444 : }
1445 :
1446 21 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eBufType);
1447 21 : if (nXSize == nBufXSize && nYSize == nBufYSize &&
1448 20 : nDataTypeSize == nPixelSpace && nLineSpace == nPixelSpace * nBufXSize &&
1449 20 : nBandSpace == nLineSpace * nBufYSize && nBandCount == nBands)
1450 : {
1451 82 : for (int i = 0; i < nBands; i++)
1452 : {
1453 124 : if (panBandMap[i] != i + 1 ||
1454 62 : !static_cast<VRTRasterBand *>(GetRasterBand(i + 1))
1455 62 : ->IsPansharpenRasterBand())
1456 : {
1457 0 : goto default_path;
1458 : }
1459 : }
1460 :
1461 : //{static int bDone = 0; if (!bDone) printf("(2)\n"); bDone = 1; }
1462 20 : return m_poPansharpener->ProcessRegion(nXOff, nYOff, nXSize, nYSize,
1463 20 : pData, eBufType);
1464 : }
1465 :
1466 1 : default_path:
1467 : //{static int bDone = 0; if (!bDone) printf("(3)\n"); bDone = 1; }
1468 1 : return VRTDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
1469 : nBufXSize, nBufYSize, eBufType, nBandCount,
1470 : panBandMap, nPixelSpace, nLineSpace,
1471 1 : nBandSpace, psExtraArg);
1472 : }
1473 :
1474 : /************************************************************************/
1475 : /* ==================================================================== */
1476 : /* VRTPansharpenedRasterBand */
1477 : /* ==================================================================== */
1478 : /************************************************************************/
1479 :
1480 : /************************************************************************/
1481 : /* VRTPansharpenedRasterBand() */
1482 : /************************************************************************/
1483 :
1484 192 : VRTPansharpenedRasterBand::VRTPansharpenedRasterBand(GDALDataset *poDSIn,
1485 : int nBandIn,
1486 192 : GDALDataType eDataTypeIn)
1487 192 : : m_nIndexAsPansharpenedBand(nBandIn - 1)
1488 : {
1489 192 : Initialize(poDSIn->GetRasterXSize(), poDSIn->GetRasterYSize());
1490 :
1491 192 : poDS = poDSIn;
1492 192 : nBand = nBandIn;
1493 192 : eAccess = GA_Update;
1494 192 : eDataType = eDataTypeIn;
1495 :
1496 192 : static_cast<VRTPansharpenedDataset *>(poDS)->GetBlockSize(&nBlockXSize,
1497 : &nBlockYSize);
1498 192 : }
1499 :
1500 : /************************************************************************/
1501 : /* ~VRTPansharpenedRasterBand() */
1502 : /************************************************************************/
1503 :
1504 384 : VRTPansharpenedRasterBand::~VRTPansharpenedRasterBand()
1505 :
1506 : {
1507 192 : FlushCache(true);
1508 384 : }
1509 :
1510 : /************************************************************************/
1511 : /* IReadBlock() */
1512 : /************************************************************************/
1513 :
1514 26 : CPLErr VRTPansharpenedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
1515 : void *pImage)
1516 :
1517 : {
1518 26 : const int nReqXOff = nBlockXOff * nBlockXSize;
1519 26 : const int nReqYOff = nBlockYOff * nBlockYSize;
1520 26 : int nReqXSize = nBlockXSize;
1521 26 : int nReqYSize = nBlockYSize;
1522 26 : if (nReqXOff + nReqXSize > nRasterXSize)
1523 13 : nReqXSize = nRasterXSize - nReqXOff;
1524 26 : if (nReqYOff + nReqYSize > nRasterYSize)
1525 22 : nReqYSize = nRasterYSize - nReqYOff;
1526 :
1527 : //{static int bDone = 0; if (!bDone) printf("(4)\n"); bDone = 1; }
1528 26 : const int nDataTypeSize = GDALGetDataTypeSize(eDataType) / 8;
1529 : GDALRasterIOExtraArg sExtraArg;
1530 26 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1531 52 : if (IRasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize, pImage,
1532 : nReqXSize, nReqYSize, eDataType, nDataTypeSize,
1533 26 : static_cast<GSpacing>(nDataTypeSize) * nReqXSize,
1534 26 : &sExtraArg) != CE_None)
1535 : {
1536 0 : return CE_Failure;
1537 : }
1538 :
1539 26 : if (nReqXSize < nBlockXSize)
1540 : {
1541 5613 : for (int j = nReqYSize - 1; j >= 0; j--)
1542 : {
1543 5600 : memmove(static_cast<GByte *>(pImage) +
1544 5600 : static_cast<size_t>(j) * nDataTypeSize * nBlockXSize,
1545 5600 : static_cast<GByte *>(pImage) +
1546 5600 : static_cast<size_t>(j) * nDataTypeSize * nReqXSize,
1547 5600 : static_cast<size_t>(nReqXSize) * nDataTypeSize);
1548 5600 : memset(static_cast<GByte *>(pImage) +
1549 5600 : (static_cast<size_t>(j) * nBlockXSize + nReqXSize) *
1550 5600 : nDataTypeSize,
1551 : 0,
1552 5600 : static_cast<size_t>(nBlockXSize - nReqXSize) *
1553 5600 : nDataTypeSize);
1554 : }
1555 : }
1556 26 : if (nReqYSize < nBlockYSize)
1557 : {
1558 22 : memset(static_cast<GByte *>(pImage) +
1559 22 : static_cast<size_t>(nReqYSize) * nBlockXSize * nDataTypeSize,
1560 : 0,
1561 22 : static_cast<size_t>(nBlockYSize - nReqYSize) * nBlockXSize *
1562 22 : nDataTypeSize);
1563 : }
1564 :
1565 : // Cache other bands
1566 26 : CPLErr eErr = CE_None;
1567 26 : VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);
1568 26 : if (poGDS->nBands != 1 && !poGDS->m_bLoadingOtherBands)
1569 : {
1570 10 : poGDS->m_bLoadingOtherBands = TRUE;
1571 :
1572 36 : for (int iOtherBand = 1; iOtherBand <= poGDS->nBands; iOtherBand++)
1573 : {
1574 26 : if (iOtherBand == nBand)
1575 10 : continue;
1576 :
1577 : GDALRasterBlock *poBlock =
1578 16 : poGDS->GetRasterBand(iOtherBand)
1579 16 : ->GetLockedBlockRef(nBlockXOff, nBlockYOff);
1580 16 : if (poBlock == nullptr)
1581 : {
1582 0 : eErr = CE_Failure;
1583 0 : break;
1584 : }
1585 16 : poBlock->DropLock();
1586 : }
1587 :
1588 10 : poGDS->m_bLoadingOtherBands = FALSE;
1589 : }
1590 :
1591 26 : return eErr;
1592 : }
1593 :
1594 : /************************************************************************/
1595 : /* IRasterIO() */
1596 : /************************************************************************/
1597 :
1598 127 : CPLErr VRTPansharpenedRasterBand::IRasterIO(
1599 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
1600 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
1601 : GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
1602 : {
1603 127 : if (eRWFlag == GF_Write)
1604 0 : return CE_Failure;
1605 :
1606 127 : VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);
1607 :
1608 : /* Try to pass the request to the most appropriate overview dataset */
1609 127 : if (nBufXSize < nXSize && nBufYSize < nYSize)
1610 : {
1611 : int bTried;
1612 3 : CPLErr eErr = TryOverviewRasterIO(
1613 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1614 : eBufType, nPixelSpace, nLineSpace, psExtraArg, &bTried);
1615 3 : if (bTried)
1616 0 : return eErr;
1617 : }
1618 :
1619 127 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eBufType);
1620 127 : if (nDataTypeSize > 0 && nXSize == nBufXSize && nYSize == nBufYSize &&
1621 124 : nDataTypeSize == nPixelSpace && nLineSpace == nPixelSpace * nBufXSize)
1622 : {
1623 : GDALPansharpenOptions *psOptions =
1624 124 : poGDS->m_poPansharpener->GetOptions();
1625 :
1626 : // Have we already done this request for another band ?
1627 : // If so use the cached result
1628 124 : const size_t nBufferSizePerBand =
1629 124 : static_cast<size_t>(nXSize) * nYSize * nDataTypeSize;
1630 124 : if (nXOff == poGDS->m_nLastBandRasterIOXOff &&
1631 118 : nYOff >= poGDS->m_nLastBandRasterIOYOff &&
1632 111 : nXSize == poGDS->m_nLastBandRasterIOXSize &&
1633 69 : nYOff + nYSize <= poGDS->m_nLastBandRasterIOYOff +
1634 69 : poGDS->m_nLastBandRasterIOYSize &&
1635 59 : eBufType == poGDS->m_eLastBandRasterIODataType)
1636 : {
1637 : //{static int bDone = 0; if (!bDone) printf("(6)\n"); bDone = 1; }
1638 53 : if (poGDS->m_pabyLastBufferBandRasterIO == nullptr)
1639 0 : return CE_Failure;
1640 53 : const size_t nBufferSizePerBandCached =
1641 53 : static_cast<size_t>(nXSize) * poGDS->m_nLastBandRasterIOYSize *
1642 53 : nDataTypeSize;
1643 53 : memcpy(pData,
1644 53 : poGDS->m_pabyLastBufferBandRasterIO +
1645 53 : nBufferSizePerBandCached * m_nIndexAsPansharpenedBand +
1646 53 : static_cast<size_t>(nYOff -
1647 53 : poGDS->m_nLastBandRasterIOYOff) *
1648 53 : nXSize * nDataTypeSize,
1649 : nBufferSizePerBand);
1650 53 : return CE_None;
1651 : }
1652 :
1653 71 : int nYSizeToCache = nYSize;
1654 71 : if (nYSize == 1 && nXSize == nRasterXSize)
1655 : {
1656 : //{static int bDone = 0; if (!bDone) printf("(7)\n"); bDone = 1; }
1657 : // For efficiency, try to cache at leak 256 K
1658 0 : nYSizeToCache = (256 * 1024) / nXSize / nDataTypeSize;
1659 0 : if (nYSizeToCache == 0)
1660 0 : nYSizeToCache = 1;
1661 0 : else if (nYOff + nYSizeToCache > nRasterYSize)
1662 0 : nYSizeToCache = nRasterYSize - nYOff;
1663 : }
1664 71 : const GUIntBig nBufferSize = static_cast<GUIntBig>(nXSize) *
1665 71 : nYSizeToCache * nDataTypeSize *
1666 71 : psOptions->nOutPansharpenedBands;
1667 : // Check the we don't overflow (for 32 bit platforms)
1668 71 : if (nBufferSize > std::numeric_limits<size_t>::max() / 2)
1669 : {
1670 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
1671 : "Out of memory error while allocating working buffers");
1672 0 : return CE_Failure;
1673 : }
1674 : GByte *pabyTemp = static_cast<GByte *>(
1675 71 : VSI_REALLOC_VERBOSE(poGDS->m_pabyLastBufferBandRasterIO,
1676 : static_cast<size_t>(nBufferSize)));
1677 71 : if (pabyTemp == nullptr)
1678 : {
1679 0 : return CE_Failure;
1680 : }
1681 71 : poGDS->m_nLastBandRasterIOXOff = nXOff;
1682 71 : poGDS->m_nLastBandRasterIOYOff = nYOff;
1683 71 : poGDS->m_nLastBandRasterIOXSize = nXSize;
1684 71 : poGDS->m_nLastBandRasterIOYSize = nYSizeToCache;
1685 71 : poGDS->m_eLastBandRasterIODataType = eBufType;
1686 71 : poGDS->m_pabyLastBufferBandRasterIO = pabyTemp;
1687 :
1688 142 : CPLErr eErr = poGDS->m_poPansharpener->ProcessRegion(
1689 : nXOff, nYOff, nXSize, nYSizeToCache,
1690 71 : poGDS->m_pabyLastBufferBandRasterIO, eBufType);
1691 71 : if (eErr == CE_None)
1692 : {
1693 : //{static int bDone = 0; if (!bDone) printf("(8)\n"); bDone = 1; }
1694 71 : size_t nBufferSizePerBandCached = static_cast<size_t>(nXSize) *
1695 71 : poGDS->m_nLastBandRasterIOYSize *
1696 71 : nDataTypeSize;
1697 71 : memcpy(pData,
1698 71 : poGDS->m_pabyLastBufferBandRasterIO +
1699 71 : nBufferSizePerBandCached * m_nIndexAsPansharpenedBand,
1700 : nBufferSizePerBand);
1701 : }
1702 : else
1703 : {
1704 0 : VSIFree(poGDS->m_pabyLastBufferBandRasterIO);
1705 0 : poGDS->m_pabyLastBufferBandRasterIO = nullptr;
1706 : }
1707 71 : return eErr;
1708 : }
1709 :
1710 : //{static int bDone = 0; if (!bDone) printf("(9)\n"); bDone = 1; }
1711 3 : CPLErr eErr = VRTRasterBand::IRasterIO(
1712 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1713 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
1714 :
1715 3 : return eErr;
1716 : }
1717 :
1718 : /************************************************************************/
1719 : /* SerializeToXML() */
1720 : /************************************************************************/
1721 :
1722 : CPLXMLNode *
1723 20 : VRTPansharpenedRasterBand::SerializeToXML(const char *pszVRTPathIn,
1724 : bool &bHasWarnedAboutRAMUsage,
1725 : size_t &nAccRAMUsage)
1726 :
1727 : {
1728 20 : CPLXMLNode *psTree = VRTRasterBand::SerializeToXML(
1729 : pszVRTPathIn, bHasWarnedAboutRAMUsage, nAccRAMUsage);
1730 :
1731 : /* -------------------------------------------------------------------- */
1732 : /* Set subclass. */
1733 : /* -------------------------------------------------------------------- */
1734 20 : CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
1735 : CXT_Text, "VRTPansharpenedRasterBand");
1736 :
1737 20 : return psTree;
1738 : }
1739 :
1740 : /************************************************************************/
1741 : /* GetOverviewCount() */
1742 : /************************************************************************/
1743 :
1744 23 : int VRTPansharpenedRasterBand::GetOverviewCount()
1745 : {
1746 23 : VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);
1747 :
1748 : // Build on-the-fly overviews from overviews of pan and spectral bands
1749 23 : if (poGDS->m_poPansharpener != nullptr &&
1750 46 : poGDS->m_apoOverviewDatasets.empty() && poGDS->m_poMainDataset == poGDS)
1751 : {
1752 : GDALPansharpenOptions *psOptions =
1753 19 : poGDS->m_poPansharpener->GetOptions();
1754 :
1755 : GDALRasterBand *poPanBand =
1756 19 : GDALRasterBand::FromHandle(psOptions->hPanchroBand);
1757 19 : const int nPanOvrCount = poPanBand->GetOverviewCount();
1758 19 : if (nPanOvrCount > 0)
1759 : {
1760 14 : for (int i = 0; i < poGDS->GetRasterCount(); i++)
1761 : {
1762 10 : if (!static_cast<VRTRasterBand *>(poGDS->GetRasterBand(i + 1))
1763 10 : ->IsPansharpenRasterBand())
1764 : {
1765 0 : return 0;
1766 : }
1767 : }
1768 :
1769 : int nSpectralOvrCount =
1770 4 : GDALRasterBand::FromHandle(psOptions->pahInputSpectralBands[0])
1771 4 : ->GetOverviewCount();
1772 10 : for (int i = 1; i < psOptions->nInputSpectralBands; i++)
1773 : {
1774 6 : auto poSpectralBand = GDALRasterBand::FromHandle(
1775 6 : psOptions->pahInputSpectralBands[i]);
1776 6 : if (poSpectralBand->GetOverviewCount() != nSpectralOvrCount)
1777 : {
1778 0 : return 0;
1779 : }
1780 : }
1781 4 : auto poPanBandDS = poPanBand->GetDataset();
1782 7 : for (int j = 0; j < std::min(nPanOvrCount, nSpectralOvrCount); j++)
1783 : {
1784 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>
1785 3 : poPanOvrDS(GDALCreateOverviewDataset(poPanBandDS, j, true));
1786 3 : if (!poPanOvrDS)
1787 : {
1788 0 : CPLError(CE_Warning, CPLE_AppDefined,
1789 : "GDALCreateOverviewDataset(poPanBandDS, %d, true) "
1790 : "failed",
1791 : j);
1792 0 : break;
1793 : }
1794 : GDALRasterBand *poPanOvrBand =
1795 3 : poPanOvrDS->GetRasterBand(poPanBand->GetBand());
1796 : auto poOvrDS = std::make_unique<VRTPansharpenedDataset>(
1797 3 : poPanOvrBand->GetXSize(), poPanOvrBand->GetYSize());
1798 6 : poOvrDS->m_apoDatasetsToReleaseRef.push_back(
1799 3 : std::move(poPanOvrDS));
1800 10 : for (int i = 0; i < poGDS->GetRasterCount(); i++)
1801 : {
1802 7 : GDALRasterBand *poSrcBand = poGDS->GetRasterBand(i + 1);
1803 : GDALRasterBand *poBand = new VRTPansharpenedRasterBand(
1804 7 : poOvrDS.get(), i + 1, poSrcBand->GetRasterDataType());
1805 : const char *pszNBITS =
1806 7 : poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1807 7 : if (pszNBITS)
1808 0 : poBand->SetMetadataItem("NBITS", pszNBITS,
1809 0 : "IMAGE_STRUCTURE");
1810 7 : poOvrDS->SetBand(i + 1, poBand);
1811 : }
1812 :
1813 : GDALPansharpenOptions *psPanOvrOptions =
1814 3 : GDALClonePansharpenOptions(psOptions);
1815 3 : psPanOvrOptions->hPanchroBand = poPanOvrBand;
1816 10 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
1817 : {
1818 7 : auto poSpectralBand = GDALRasterBand::FromHandle(
1819 7 : psOptions->pahInputSpectralBands[i]);
1820 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>
1821 : poSpectralOvrDS(GDALCreateOverviewDataset(
1822 7 : poSpectralBand->GetDataset(), j, true));
1823 7 : if (!poSpectralOvrDS)
1824 : {
1825 0 : CPLError(CE_Warning, CPLE_AppDefined,
1826 : "GDALCreateOverviewDataset(poSpectralBand->"
1827 : "GetDataset(), %d, true) failed",
1828 : j);
1829 0 : GDALDestroyPansharpenOptions(psPanOvrOptions);
1830 0 : return 0;
1831 : }
1832 14 : psPanOvrOptions->pahInputSpectralBands[i] =
1833 7 : poSpectralOvrDS->GetRasterBand(
1834 : poSpectralBand->GetBand());
1835 14 : poOvrDS->m_apoDatasetsToReleaseRef.push_back(
1836 7 : std::move(poSpectralOvrDS));
1837 : }
1838 3 : poOvrDS->m_poPansharpener = new GDALPansharpenOperation();
1839 3 : if (poOvrDS->m_poPansharpener->Initialize(psPanOvrOptions) !=
1840 : CE_None)
1841 : {
1842 0 : CPLError(CE_Warning, CPLE_AppDefined,
1843 : "Unable to initialize pansharpener.");
1844 0 : GDALDestroyPansharpenOptions(psPanOvrOptions);
1845 0 : return 0;
1846 : }
1847 3 : GDALDestroyPansharpenOptions(psPanOvrOptions);
1848 :
1849 3 : poOvrDS->m_poMainDataset = poGDS;
1850 3 : poOvrDS->SetMetadataItem("INTERLEAVE", "PIXEL",
1851 3 : "IMAGE_STRUCTURE");
1852 :
1853 3 : poGDS->m_apoOverviewDatasets.push_back(std::move(poOvrDS));
1854 : }
1855 : }
1856 : }
1857 23 : return static_cast<int>(poGDS->m_apoOverviewDatasets.size());
1858 : }
1859 :
1860 : /************************************************************************/
1861 : /* GetOverview() */
1862 : /************************************************************************/
1863 :
1864 6 : GDALRasterBand *VRTPansharpenedRasterBand::GetOverview(int iOvr)
1865 : {
1866 6 : if (iOvr < 0 || iOvr >= GetOverviewCount())
1867 2 : return nullptr;
1868 :
1869 4 : VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);
1870 :
1871 4 : return poGDS->m_apoOverviewDatasets[iOvr]->GetRasterBand(nBand);
1872 : }
1873 :
1874 : /*! @endcond */
|
