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