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 : if (poPanBand == nullptr)
340 72 : poPanBand = poPanDataset->GetRasterBand(nPanBand);
341 72 : if (poPanBand == nullptr)
342 : {
343 1 : CPLError(CE_Failure, CPLE_AppDefined, "%s invalid band of %s",
344 : pszSourceBand, osSourceFilename.c_str());
345 1 : GDALClose(poPanDataset);
346 1 : return CE_Failure;
347 : }
348 71 : oMapNamesToDataset[osSourceFilename] = poPanDataset;
349 : }
350 : else
351 : {
352 17 : poPanBand = GDALRasterBand::FromHandle(hPanchroBandIn);
353 17 : nPanBand = poPanBand->GetBand();
354 17 : poPanDataset = poPanBand->GetDataset();
355 17 : if (poPanDataset == nullptr)
356 : {
357 0 : CPLError(
358 : CE_Failure, CPLE_AppDefined,
359 : "Cannot retrieve dataset associated with panchromatic band");
360 0 : return CE_Failure;
361 : }
362 17 : oMapNamesToDataset[CPLSPrintf("%p", poPanDataset)] = poPanDataset;
363 17 : poPanDataset->Reference();
364 : }
365 88 : m_apoDatasetsToReleaseRef.push_back(
366 176 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>(
367 : poPanDataset));
368 88 : poPanDataset = m_apoDatasetsToReleaseRef.back().get();
369 : // Figure out which kind of adjustment we should do if the pan and spectral
370 : // bands do not share the same geotransform
371 : const char *pszGTAdjustment =
372 88 : CPLGetXMLValue(psOptions, "SpatialExtentAdjustment", "Union");
373 88 : if (EQUAL(pszGTAdjustment, "Union"))
374 82 : m_eGTAdjustment = GTAdjust_Union;
375 6 : else if (EQUAL(pszGTAdjustment, "Intersection"))
376 2 : m_eGTAdjustment = GTAdjust_Intersection;
377 4 : else if (EQUAL(pszGTAdjustment, "None"))
378 2 : m_eGTAdjustment = GTAdjust_None;
379 2 : else if (EQUAL(pszGTAdjustment, "NoneWithoutWarning"))
380 1 : m_eGTAdjustment = GTAdjust_NoneWithoutWarning;
381 : else
382 : {
383 1 : m_eGTAdjustment = GTAdjust_Union;
384 1 : CPLError(CE_Failure, CPLE_AppDefined,
385 : "Unsupported value for GeoTransformAdjustment. Defaulting to "
386 : "Union");
387 : }
388 :
389 : const char *pszNumThreads =
390 88 : CPLGetXMLValue(psOptions, "NumThreads", nullptr);
391 88 : int nThreads = 0;
392 88 : if (pszNumThreads != nullptr)
393 : {
394 41 : if (EQUAL(pszNumThreads, "ALL_CPUS"))
395 41 : nThreads = -1;
396 : else
397 0 : nThreads = atoi(pszNumThreads);
398 : }
399 :
400 : const char *pszAlgorithm =
401 88 : CPLGetXMLValue(psOptions, "Algorithm", "WeightedBrovey");
402 88 : if (!EQUAL(pszAlgorithm, "WeightedBrovey"))
403 : {
404 1 : CPLError(CE_Failure, CPLE_AppDefined, "Algorithm %s unsupported",
405 : pszAlgorithm);
406 1 : return CE_Failure;
407 : }
408 :
409 174 : std::vector<double> adfWeights;
410 87 : if (const CPLXMLNode *psAlgOptions =
411 87 : CPLGetXMLNode(psOptions, "AlgorithmOptions"))
412 : {
413 : const char *pszWeights =
414 37 : CPLGetXMLValue(psAlgOptions, "Weights", nullptr);
415 37 : if (pszWeights != nullptr)
416 : {
417 37 : char **papszTokens = CSLTokenizeString2(pszWeights, " ,", 0);
418 128 : for (int i = 0; papszTokens && papszTokens[i]; i++)
419 91 : adfWeights.push_back(CPLAtof(papszTokens[i]));
420 37 : CSLDestroy(papszTokens);
421 : }
422 : }
423 :
424 87 : GDALRIOResampleAlg eResampleAlg = GDALRasterIOGetResampleAlg(
425 : CPLGetXMLValue(psOptions, "Resampling", "Cubic"));
426 :
427 174 : std::vector<GDALRasterBand *> ahSpectralBands;
428 174 : std::map<int, int> aMapDstBandToSpectralBand;
429 87 : std::map<int, int>::iterator aMapDstBandToSpectralBandIter;
430 87 : int nBitDepth = 0;
431 87 : bool bFoundNonMatchingGT = false;
432 87 : GDALGeoTransform panGT;
433 : const bool bPanGeoTransformValid =
434 87 : (poPanDataset->GetGeoTransform(panGT) == CE_None);
435 87 : if (!bPanGeoTransformValid)
436 : {
437 1 : CPLError(CE_Failure, CPLE_AppDefined,
438 : "Panchromatic band has no associated geotransform");
439 1 : return CE_Failure;
440 : }
441 86 : int nPanXSize = poPanBand->GetXSize();
442 86 : int nPanYSize = poPanBand->GetYSize();
443 86 : int bHasNoData = FALSE;
444 86 : double dfNoData = poPanBand->GetNoDataValue(&bHasNoData);
445 86 : double dfLRPanX = panGT[0] + nPanXSize * panGT[1] + nPanYSize * panGT[2];
446 86 : double dfLRPanY = panGT[3] + nPanXSize * panGT[4] + nPanYSize * panGT[5];
447 86 : bool bFoundRotatingTerms = (panGT[2] != 0.0 || panGT[4] != 0.0);
448 86 : double dfMinX = panGT[0];
449 86 : double dfMaxX = dfLRPanX;
450 86 : double dfMaxY = panGT[3];
451 86 : double dfMinY = dfLRPanY;
452 86 : if (dfMinY > dfMaxY)
453 11 : std::swap(dfMinY, dfMaxY);
454 86 : const double dfPanMinY = dfMinY;
455 86 : const double dfPanMaxY = dfMaxY;
456 :
457 86 : const auto poPanSRS = poPanDataset->GetSpatialRef();
458 :
459 : /* -------------------------------------------------------------------- */
460 : /* First pass on spectral datasets to check their georeferencing. */
461 : /* -------------------------------------------------------------------- */
462 86 : int iSpectralBand = 0;
463 555 : for (const CPLXMLNode *psIter = psOptions->psChild; psIter;
464 469 : psIter = psIter->psNext)
465 : {
466 : GDALDataset *poDataset;
467 472 : if (psIter->eType != CXT_Element ||
468 472 : !EQUAL(psIter->pszValue, "SpectralBand"))
469 264 : continue;
470 :
471 208 : if (nInputSpectralBandsIn && pahInputSpectralBandsIn != nullptr)
472 : {
473 42 : if (iSpectralBand == nInputSpectralBandsIn)
474 : {
475 1 : CPLError(CE_Failure, CPLE_AppDefined,
476 : "More SpectralBand elements than in source array");
477 3 : goto error;
478 : }
479 : poDataset = GDALRasterBand::FromHandle(
480 41 : pahInputSpectralBandsIn[iSpectralBand])
481 41 : ->GetDataset();
482 41 : if (poDataset == nullptr)
483 : {
484 0 : CPLError(CE_Failure, CPLE_AppDefined,
485 : "SpectralBand has no associated dataset");
486 0 : goto error;
487 : }
488 41 : osSourceFilename = poDataset->GetDescription();
489 :
490 41 : oMapNamesToDataset[CPLSPrintf("%p", poDataset)] = poDataset;
491 41 : poDataset->Reference();
492 : }
493 : else
494 : {
495 166 : osSourceFilename = CPLGetXMLValue(psIter, "SourceFilename", "");
496 166 : if (osSourceFilename.empty())
497 : {
498 1 : CPLError(CE_Failure, CPLE_AppDefined,
499 : "SpectralBand.SourceFilename missing");
500 1 : goto error;
501 : }
502 165 : int bRelativeToVRT = atoi(
503 : CPLGetXMLValue(psIter, "SourceFilename.relativetoVRT", "0"));
504 165 : if (bRelativeToVRT)
505 : {
506 : const std::string osAbs = CPLProjectRelativeFilenameSafe(
507 107 : pszVRTPathIn, osSourceFilename.c_str());
508 107 : m_oMapToRelativeFilenames[osAbs] = osSourceFilename;
509 107 : osSourceFilename = osAbs;
510 : }
511 165 : poDataset = oMapNamesToDataset[osSourceFilename];
512 165 : if (poDataset == nullptr)
513 : {
514 : const CPLStringList aosOpenOptions(
515 70 : GDALDeserializeOpenOptionsFromXML(psIter));
516 :
517 70 : poDataset = GDALDataset::Open(
518 : osSourceFilename, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
519 : nullptr, aosOpenOptions.List());
520 70 : if (poDataset == nullptr)
521 : {
522 1 : goto error;
523 : }
524 69 : oMapNamesToDataset[osSourceFilename] = poDataset;
525 : }
526 : else
527 : {
528 95 : poDataset->Reference();
529 : }
530 : }
531 205 : m_apoDatasetsToReleaseRef.push_back(
532 410 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>(
533 : poDataset));
534 205 : poDataset = m_apoDatasetsToReleaseRef.back().get();
535 :
536 : // Check that the spectral band has a georeferencing consistent
537 : // of the pan band. Allow an error of at most the size of one pixel
538 : // of the spectral band.
539 205 : const auto poSpectralSRS = poDataset->GetSpatialRef();
540 205 : if (poPanSRS)
541 : {
542 164 : if (poSpectralSRS)
543 : {
544 164 : if (!poPanSRS->IsSame(poSpectralSRS))
545 : {
546 1 : CPLError(CE_Warning, CPLE_AppDefined,
547 : "Pan dataset and %s do not seem to "
548 : "have same projection. Results might "
549 : "be incorrect",
550 : osSourceFilename.c_str());
551 : }
552 : }
553 : else
554 : {
555 0 : CPLError(CE_Warning, CPLE_AppDefined,
556 : "Pan dataset has a projection, whereas %s "
557 : "not. Results might be incorrect",
558 : osSourceFilename.c_str());
559 : }
560 : }
561 41 : else if (poSpectralSRS)
562 : {
563 0 : CPLError(CE_Warning, CPLE_AppDefined,
564 : "Pan dataset has no projection, whereas %s has "
565 : "one. Results might be incorrect",
566 : osSourceFilename.c_str());
567 : }
568 :
569 205 : GDALGeoTransform spectralGT;
570 205 : if (poDataset->GetGeoTransform(spectralGT) != CE_None)
571 : {
572 0 : CPLError(CE_Failure, CPLE_AppDefined,
573 : "Spectral band has no associated geotransform");
574 0 : goto error;
575 : }
576 205 : if (spectralGT[5] * panGT[5] < 0)
577 : {
578 0 : CPLError(CE_Failure, CPLE_NotSupported,
579 : "Spectral band vertical orientation is "
580 : "different from pan one");
581 0 : goto error;
582 : }
583 :
584 205 : int bIsThisOneNonMatching = FALSE;
585 205 : double dfPixelSize = std::max(spectralGT[1], std::abs(spectralGT[5]));
586 381 : if (std::abs(panGT[0] - spectralGT[0]) > dfPixelSize ||
587 176 : std::abs(panGT[3] - spectralGT[3]) > dfPixelSize)
588 : {
589 29 : bIsThisOneNonMatching = TRUE;
590 29 : if (m_eGTAdjustment == GTAdjust_None)
591 : {
592 2 : CPLError(CE_Warning, CPLE_AppDefined,
593 : "Georeferencing of top-left corner of pan "
594 : "dataset and %s do not match",
595 : osSourceFilename.c_str());
596 : }
597 : }
598 410 : bFoundRotatingTerms = bFoundRotatingTerms ||
599 205 : (spectralGT[2] != 0.0 || spectralGT[4] != 0.0);
600 205 : double dfLRSpectralX = spectralGT[0] +
601 205 : poDataset->GetRasterXSize() * spectralGT[1] +
602 205 : poDataset->GetRasterYSize() * spectralGT[2];
603 205 : double dfLRSpectralY = spectralGT[3] +
604 205 : poDataset->GetRasterXSize() * spectralGT[4] +
605 205 : poDataset->GetRasterYSize() * spectralGT[5];
606 381 : if (std::abs(dfLRPanX - dfLRSpectralX) > dfPixelSize ||
607 176 : std::abs(dfLRPanY - dfLRSpectralY) > dfPixelSize)
608 : {
609 29 : bIsThisOneNonMatching = TRUE;
610 29 : if (m_eGTAdjustment == GTAdjust_None)
611 : {
612 2 : CPLError(CE_Warning, CPLE_AppDefined,
613 : "Georeferencing of bottom-right corner of "
614 : "pan dataset and %s do not match",
615 : osSourceFilename.c_str());
616 : }
617 : }
618 :
619 205 : double dfThisMinY = dfLRSpectralY;
620 205 : double dfThisMaxY = spectralGT[3];
621 205 : if (dfThisMinY > dfThisMaxY)
622 23 : std::swap(dfThisMinY, dfThisMaxY);
623 205 : if (bIsThisOneNonMatching && m_eGTAdjustment == GTAdjust_Union)
624 : {
625 24 : dfMinX = std::min(dfMinX, spectralGT[0]);
626 24 : dfMinY = std::min(dfMinY, dfThisMinY);
627 24 : dfMaxX = std::max(dfMaxX, dfLRSpectralX);
628 24 : dfMaxY = std::max(dfMaxY, dfThisMaxY);
629 : }
630 181 : else if (bIsThisOneNonMatching &&
631 5 : m_eGTAdjustment == GTAdjust_Intersection)
632 : {
633 2 : dfMinX = std::max(dfMinX, spectralGT[0]);
634 2 : dfMinY = std::max(dfMinY, dfThisMinY);
635 2 : dfMaxX = std::min(dfMaxX, dfLRSpectralX);
636 2 : dfMaxY = std::min(dfMaxY, dfThisMaxY);
637 : }
638 :
639 205 : bFoundNonMatchingGT = bFoundNonMatchingGT || bIsThisOneNonMatching;
640 :
641 205 : iSpectralBand++;
642 : }
643 :
644 83 : if (nInputSpectralBandsIn && iSpectralBand != nInputSpectralBandsIn)
645 : {
646 1 : CPLError(CE_Failure, CPLE_AppDefined,
647 : "Less SpectralBand elements than in source array");
648 1 : goto error;
649 : }
650 :
651 : /* -------------------------------------------------------------------- */
652 : /* On-the-fly spatial extent adjustment if needed and asked. */
653 : /* -------------------------------------------------------------------- */
654 :
655 82 : if (bFoundNonMatchingGT && (m_eGTAdjustment == GTAdjust_Union ||
656 5 : m_eGTAdjustment == GTAdjust_Intersection))
657 : {
658 13 : if (bFoundRotatingTerms)
659 : {
660 0 : CPLError(
661 : CE_Failure, CPLE_NotSupported,
662 : "One of the panchromatic or spectral datasets has rotating "
663 : "terms in their geotransform matrix. Adjustment not possible");
664 2 : goto error;
665 : }
666 13 : if (m_eGTAdjustment == GTAdjust_Intersection &&
667 2 : (dfMinX >= dfMaxX || dfMinY >= dfMaxY))
668 : {
669 1 : CPLError(
670 : CE_Failure, CPLE_NotSupported,
671 : "One of the panchromatic or spectral datasets has rotating "
672 : "terms in their geotransform matrix. Adjustment not possible");
673 1 : goto error;
674 : }
675 12 : if (m_eGTAdjustment == GTAdjust_Union)
676 11 : CPLDebug("VRT", "Do union of bounding box of panchromatic and "
677 : "spectral datasets");
678 : else
679 1 : CPLDebug("VRT", "Do intersection of bounding box of panchromatic "
680 : "and spectral datasets");
681 :
682 : // If the pandataset needs adjustments, make sure the coordinates of the
683 : // union/intersection properly align with the grid of the pandataset
684 : // to avoid annoying sub-pixel shifts on the panchro band.
685 12 : double dfPixelSize = std::max(panGT[1], std::abs(panGT[5]));
686 12 : if (std::abs(panGT[0] - dfMinX) > dfPixelSize ||
687 7 : std::abs(dfPanMaxY - dfMaxY) > dfPixelSize ||
688 21 : std::abs(dfLRPanX - dfMaxX) > dfPixelSize ||
689 2 : std::abs(dfPanMinY - dfMinY) > dfPixelSize)
690 : {
691 10 : dfMinX =
692 10 : panGT[0] +
693 10 : std::floor((dfMinX - panGT[0]) / panGT[1] + 0.5) * panGT[1];
694 10 : dfMaxX =
695 10 : dfLRPanX +
696 10 : std::floor((dfMaxX - dfLRPanX) / panGT[1] + 0.5) * panGT[1];
697 10 : if (panGT[5] > 0)
698 : {
699 0 : dfMinY = dfPanMinY +
700 0 : std::floor((dfMinY - dfPanMinY) / panGT[5] + 0.5) *
701 0 : panGT[5];
702 0 : dfMaxY = dfPanMinY +
703 0 : std::floor((dfMaxY - dfPanMinY) / panGT[5] + 0.5) *
704 0 : panGT[5];
705 : }
706 : else
707 : {
708 10 : dfMinY = dfPanMaxY +
709 10 : std::floor((dfMinY - dfPanMaxY) / panGT[5] + 0.5) *
710 10 : panGT[5];
711 10 : dfMaxY = dfPanMaxY +
712 20 : std::floor((dfMaxY - dfPanMaxY) / panGT[5] + 0.5) *
713 10 : panGT[5];
714 : }
715 : }
716 :
717 : std::map<CPLString, GDALDataset *>::iterator oIter =
718 12 : oMapNamesToDataset.begin();
719 34 : for (; oIter != oMapNamesToDataset.end(); ++oIter)
720 : {
721 23 : GDALDataset *poSrcDS = oIter->second;
722 23 : GDALGeoTransform gt;
723 23 : if (poSrcDS->GetGeoTransform(gt) != CE_None)
724 3 : continue;
725 :
726 : // Check if this dataset needs adjustments
727 23 : dfPixelSize = std::max(gt.xscale, std::abs(gt.yscale));
728 23 : dfPixelSize = std::max(panGT[1], dfPixelSize);
729 23 : dfPixelSize = std::max(std::abs(panGT[5]), dfPixelSize);
730 23 : double dfThisMinX = gt.xorig;
731 : double dfThisMaxX =
732 23 : gt.xorig + poSrcDS->GetRasterXSize() * gt.xscale;
733 23 : double dfThisMaxY = gt.yorig;
734 : double dfThisMinY =
735 23 : gt.yorig + poSrcDS->GetRasterYSize() * gt.yscale;
736 23 : if (dfThisMinY > dfThisMaxY)
737 2 : std::swap(dfThisMinY, dfThisMaxY);
738 23 : if (std::abs(dfThisMinX - dfMinX) <= dfPixelSize &&
739 11 : std::abs(dfThisMaxY - dfMaxY) <= dfPixelSize &&
740 37 : std::abs(dfThisMaxX - dfMaxX) <= dfPixelSize &&
741 3 : std::abs(dfThisMinY - dfMinY) <= dfPixelSize)
742 : {
743 3 : continue;
744 : }
745 :
746 20 : GDALGeoTransform adjustedGT;
747 20 : adjustedGT[0] = dfMinX;
748 20 : adjustedGT[1] = gt.xscale;
749 20 : adjustedGT[2] = 0;
750 20 : adjustedGT[3] = gt.yscale > 0 ? dfMinY : dfMaxY;
751 20 : adjustedGT[4] = 0;
752 20 : adjustedGT[5] = gt.yscale;
753 : const double dfAdjustedRasterXSize =
754 20 : 0.5 + (dfMaxX - dfMinX) / adjustedGT[1];
755 : const double dfAdjustedRasterYSize =
756 20 : 0.5 + (dfMaxY - dfMinY) / std::abs(adjustedGT[5]);
757 20 : if (!(dfAdjustedRasterXSize >= 1 &&
758 20 : dfAdjustedRasterXSize < INT_MAX &&
759 19 : dfAdjustedRasterYSize >= 1 &&
760 19 : dfAdjustedRasterYSize < INT_MAX))
761 : {
762 1 : CPLError(CE_Failure, CPLE_AppDefined,
763 : "Datasets are too disjoint");
764 1 : goto error;
765 : }
766 19 : const int nAdjustRasterXSize =
767 : static_cast<int>(dfAdjustedRasterXSize);
768 19 : const int nAdjustRasterYSize =
769 : static_cast<int>(dfAdjustedRasterYSize);
770 :
771 : VRTDataset *poVDS =
772 19 : new VRTDataset(nAdjustRasterXSize, nAdjustRasterYSize);
773 19 : poVDS->SetWritable(FALSE);
774 19 : poVDS->SetMetadataItem("UNDERLYING_DATASET",
775 19 : poSrcDS->GetDescription());
776 19 : poVDS->SetGeoTransform(adjustedGT);
777 19 : poVDS->SetProjection(poPanDataset->GetProjectionRef());
778 :
779 44 : for (int i = 0; i < poSrcDS->GetRasterCount(); i++)
780 : {
781 25 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(i + 1);
782 25 : poVDS->AddBand(poSrcBand->GetRasterDataType(), nullptr);
783 : VRTSourcedRasterBand *poVRTBand =
784 : static_cast<VRTSourcedRasterBand *>(
785 25 : poVDS->GetRasterBand(i + 1));
786 :
787 : const char *pszNBITS =
788 25 : poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
789 25 : if (pszNBITS)
790 0 : poVRTBand->SetMetadataItem("NBITS", pszNBITS,
791 0 : "IMAGE_STRUCTURE");
792 :
793 25 : VRTSimpleSource *poSimpleSource = new VRTSimpleSource();
794 100 : poVRTBand->ConfigureSource(
795 : poSimpleSource, poSrcBand, FALSE,
796 25 : static_cast<int>(
797 25 : std::floor((dfMinX - gt.xorig) / gt.xscale + 0.001)),
798 25 : gt.yscale < 0
799 22 : ? static_cast<int>(std::floor(
800 22 : (dfMaxY - dfThisMaxY) / gt.yscale + 0.001))
801 3 : : static_cast<int>(std::floor(
802 3 : (dfMinY - dfThisMinY) / gt.yscale + 0.001)),
803 25 : static_cast<int>(0.5 + (dfMaxX - dfMinX) / gt.xscale),
804 25 : static_cast<int>(0.5 +
805 25 : (dfMaxY - dfMinY) / std::abs(gt.yscale)),
806 : 0, 0, nAdjustRasterXSize, nAdjustRasterYSize);
807 :
808 25 : poVRTBand->AddSource(poSimpleSource);
809 : }
810 :
811 19 : oIter->second = poVDS;
812 19 : if (poSrcDS == poPanDataset)
813 : {
814 9 : panGT = adjustedGT;
815 9 : poPanDataset = poVDS;
816 9 : poPanBand = poPanDataset->GetRasterBand(nPanBand);
817 9 : nPanXSize = poPanDataset->GetRasterXSize();
818 9 : nPanYSize = poPanDataset->GetRasterYSize();
819 : }
820 :
821 19 : m_apoDatasetsToReleaseRef.push_back(
822 38 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>(
823 : poVDS));
824 : }
825 : }
826 :
827 80 : if (nRasterXSize == 0 && nRasterYSize == 0)
828 : {
829 58 : nRasterXSize = nPanXSize;
830 58 : nRasterYSize = nPanYSize;
831 : }
832 22 : else if (nRasterXSize != nPanXSize || nRasterYSize != nPanYSize)
833 : {
834 1 : CPLError(CE_Failure, CPLE_AppDefined,
835 : "Inconsistent declared VRT dimensions with panchro dataset");
836 1 : goto error;
837 : }
838 :
839 : /* -------------------------------------------------------------------- */
840 : /* Initialize all the general VRT stuff. This will even */
841 : /* create the VRTPansharpenedRasterBands and initialize them. */
842 : /* -------------------------------------------------------------------- */
843 79 : eErr = VRTDataset::XMLInit(psTree, pszVRTPathIn);
844 :
845 79 : if (eErr != CE_None)
846 : {
847 1 : goto error;
848 : }
849 :
850 : /* -------------------------------------------------------------------- */
851 : /* Inherit georeferencing info from panchro band if not defined */
852 : /* in VRT. */
853 : /* -------------------------------------------------------------------- */
854 :
855 : {
856 78 : GDALGeoTransform outGT;
857 144 : if (GetGeoTransform(outGT) != CE_None && GetGCPCount() == 0 &&
858 66 : GetSpatialRef() == nullptr)
859 : {
860 66 : if (bPanGeoTransformValid)
861 : {
862 66 : SetGeoTransform(panGT);
863 : }
864 66 : if (poPanSRS)
865 : {
866 47 : SetSpatialRef(poPanSRS);
867 : }
868 : }
869 : }
870 :
871 : /* -------------------------------------------------------------------- */
872 : /* Parse rest of PansharpeningOptions */
873 : /* -------------------------------------------------------------------- */
874 78 : iSpectralBand = 0;
875 497 : for (const CPLXMLNode *psIter = psOptions->psChild; psIter;
876 419 : psIter = psIter->psNext)
877 : {
878 422 : if (psIter->eType != CXT_Element ||
879 422 : !EQUAL(psIter->pszValue, "SpectralBand"))
880 235 : continue;
881 :
882 : GDALDataset *poDataset;
883 : GDALRasterBand *poBand;
884 :
885 187 : const char *pszDstBand = CPLGetXMLValue(psIter, "dstBand", nullptr);
886 187 : int nDstBand = -1;
887 187 : if (pszDstBand != nullptr)
888 : {
889 181 : nDstBand = atoi(pszDstBand);
890 181 : if (nDstBand <= 0)
891 : {
892 1 : CPLError(CE_Failure, CPLE_AppDefined,
893 : "SpectralBand.dstBand = '%s' invalid", pszDstBand);
894 3 : goto error;
895 : }
896 : }
897 :
898 186 : if (nInputSpectralBandsIn && pahInputSpectralBandsIn != nullptr)
899 : {
900 70 : poBand = GDALRasterBand::FromHandle(
901 35 : pahInputSpectralBandsIn[iSpectralBand]);
902 35 : poDataset = poBand->GetDataset();
903 35 : if (poDataset)
904 : {
905 35 : poDataset = oMapNamesToDataset[CPLSPrintf("%p", poDataset)];
906 35 : CPLAssert(poDataset);
907 35 : if (poDataset)
908 35 : poBand = poDataset->GetRasterBand(poBand->GetBand());
909 : }
910 : }
911 : else
912 : {
913 151 : osSourceFilename = CPLGetXMLValue(psIter, "SourceFilename", "");
914 151 : const bool bRelativeToVRT = CPL_TO_BOOL(atoi(
915 : CPLGetXMLValue(psIter, "SourceFilename.relativetoVRT", "0")));
916 151 : if (bRelativeToVRT)
917 188 : osSourceFilename = CPLProjectRelativeFilenameSafe(
918 94 : pszVRTPathIn, osSourceFilename.c_str());
919 151 : poDataset = oMapNamesToDataset[osSourceFilename];
920 151 : CPLAssert(poDataset);
921 : const char *pszSourceBand =
922 151 : CPLGetXMLValue(psIter, "SourceBand", "1");
923 151 : const int nBand = atoi(pszSourceBand);
924 151 : poBand = poDataset->GetRasterBand(nBand);
925 151 : if (poBand == nullptr)
926 : {
927 1 : CPLError(CE_Failure, CPLE_AppDefined, "%s invalid band of %s",
928 : pszSourceBand, osSourceFilename.c_str());
929 1 : goto error;
930 : }
931 : }
932 :
933 185 : if (bHasNoData)
934 : {
935 3 : double dfSpectralNoData = poPanBand->GetNoDataValue(&bHasNoData);
936 3 : if (bHasNoData && dfSpectralNoData != dfNoData)
937 0 : bHasNoData = FALSE;
938 : }
939 :
940 185 : ahSpectralBands.push_back(poBand);
941 185 : if (nDstBand >= 1)
942 : {
943 179 : if (aMapDstBandToSpectralBand.find(nDstBand - 1) !=
944 358 : aMapDstBandToSpectralBand.end())
945 : {
946 1 : CPLError(
947 : CE_Failure, CPLE_AppDefined,
948 : "Another spectral band is already mapped to output band %d",
949 : nDstBand);
950 1 : goto error;
951 : }
952 178 : aMapDstBandToSpectralBand[nDstBand - 1] =
953 178 : static_cast<int>(ahSpectralBands.size() - 1);
954 : }
955 :
956 184 : iSpectralBand++;
957 : }
958 :
959 75 : if (ahSpectralBands.empty())
960 : {
961 1 : CPLError(CE_Failure, CPLE_AppDefined, "No spectral band defined");
962 1 : goto error;
963 : }
964 :
965 : {
966 74 : const char *pszNoData = CPLGetXMLValue(psOptions, "NoData", nullptr);
967 74 : if (pszNoData)
968 : {
969 7 : if (EQUAL(pszNoData, "NONE"))
970 : {
971 0 : m_bNoDataDisabled = TRUE;
972 0 : bHasNoData = FALSE;
973 : }
974 : else
975 : {
976 7 : bHasNoData = TRUE;
977 7 : dfNoData = CPLAtof(pszNoData);
978 : }
979 : }
980 : }
981 :
982 74 : if (GetRasterCount() == 0)
983 : {
984 156 : for (int i = 0; i < static_cast<int>(aMapDstBandToSpectralBand.size());
985 : i++)
986 : {
987 109 : if (aMapDstBandToSpectralBand.find(i) ==
988 218 : aMapDstBandToSpectralBand.end())
989 : {
990 1 : CPLError(CE_Failure, CPLE_AppDefined,
991 : "Hole in SpectralBand.dstBand numbering");
992 1 : goto error;
993 : }
994 108 : GDALRasterBand *poInputBand = GDALRasterBand::FromHandle(
995 108 : ahSpectralBands[aMapDstBandToSpectralBand[i]]);
996 : GDALRasterBand *poBand = new VRTPansharpenedRasterBand(
997 108 : this, i + 1, poInputBand->GetRasterDataType());
998 108 : poBand->SetColorInterpretation(
999 108 : poInputBand->GetColorInterpretation());
1000 108 : if (bHasNoData)
1001 13 : poBand->SetNoDataValue(dfNoData);
1002 108 : SetBand(i + 1, poBand);
1003 : }
1004 : }
1005 : else
1006 : {
1007 26 : int nIdxAsPansharpenedBand = 0;
1008 97 : for (int i = 0; i < nBands; i++)
1009 : {
1010 72 : if (static_cast<VRTRasterBand *>(GetRasterBand(i + 1))
1011 72 : ->IsPansharpenRasterBand())
1012 : {
1013 68 : if (aMapDstBandToSpectralBand.find(i) ==
1014 136 : aMapDstBandToSpectralBand.end())
1015 : {
1016 1 : CPLError(CE_Failure, CPLE_AppDefined,
1017 : "Band %d of type VRTPansharpenedRasterBand, but "
1018 : "no corresponding SpectralBand",
1019 : i + 1);
1020 1 : goto error;
1021 : }
1022 : else
1023 : {
1024 : static_cast<VRTPansharpenedRasterBand *>(
1025 67 : GetRasterBand(i + 1))
1026 67 : ->SetIndexAsPansharpenedBand(nIdxAsPansharpenedBand);
1027 67 : nIdxAsPansharpenedBand++;
1028 : }
1029 : }
1030 : }
1031 : }
1032 :
1033 : // Figure out bit depth
1034 : {
1035 : const char *pszBitDepth =
1036 72 : CPLGetXMLValue(psOptions, "BitDepth", nullptr);
1037 72 : if (pszBitDepth == nullptr)
1038 59 : pszBitDepth = GDALRasterBand::FromHandle(ahSpectralBands[0])
1039 59 : ->GetMetadataItem("NBITS", "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 : "NBITS", "IMAGE_STRUCTURE") == nullptr)
1051 : {
1052 27 : if (nBitDepth != 8 && nBitDepth != 16 && nBitDepth != 32)
1053 : {
1054 12 : GetRasterBand(i + 1)->SetMetadataItem(
1055 : "NBITS", CPLSPrintf("%d", nBitDepth),
1056 6 : "IMAGE_STRUCTURE");
1057 : }
1058 : }
1059 0 : else if (nBitDepth == 8 || nBitDepth == 16 || nBitDepth == 32)
1060 : {
1061 0 : GetRasterBand(i + 1)->SetMetadataItem("NBITS", nullptr,
1062 0 : "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("INTERLEAVE", "PIXEL", "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 : const char *pszNBITS =
1831 7 : poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1832 7 : if (pszNBITS)
1833 0 : poBand->SetMetadataItem("NBITS", pszNBITS,
1834 0 : "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("INTERLEAVE", "PIXEL",
1883 3 : "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 */
|
