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 11 : std::abs(dfThisMaxY - dfMaxY) <= dfPixelSize &&
739 37 : 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 : const char *pszNBITS =
787 25 : poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
788 25 : if (pszNBITS)
789 0 : poVRTBand->SetMetadataItem("NBITS", pszNBITS,
790 0 : "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 59 : pszBitDepth = GDALRasterBand::FromHandle(ahSpectralBands[0])
1038 59 : ->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1039 72 : if (pszBitDepth)
1040 15 : nBitDepth = atoi(pszBitDepth);
1041 72 : if (nBitDepth)
1042 : {
1043 44 : for (int i = 0; i < nBands; i++)
1044 : {
1045 29 : if (!static_cast<VRTRasterBand *>(GetRasterBand(i + 1))
1046 29 : ->IsPansharpenRasterBand())
1047 2 : continue;
1048 27 : if (GetRasterBand(i + 1)->GetMetadataItem(
1049 27 : "NBITS", "IMAGE_STRUCTURE") == nullptr)
1050 : {
1051 27 : if (nBitDepth != 8 && nBitDepth != 16 && nBitDepth != 32)
1052 : {
1053 12 : GetRasterBand(i + 1)->SetMetadataItem(
1054 : "NBITS", CPLSPrintf("%d", nBitDepth),
1055 6 : "IMAGE_STRUCTURE");
1056 : }
1057 : }
1058 0 : else if (nBitDepth == 8 || nBitDepth == 16 || nBitDepth == 32)
1059 : {
1060 0 : GetRasterBand(i + 1)->SetMetadataItem("NBITS", nullptr,
1061 0 : "IMAGE_STRUCTURE");
1062 : }
1063 : }
1064 : }
1065 : }
1066 :
1067 72 : if (GDALGetRasterBandXSize(ahSpectralBands[0]) >
1068 144 : GDALGetRasterBandXSize(poPanBand) ||
1069 72 : GDALGetRasterBandYSize(ahSpectralBands[0]) >
1070 72 : GDALGetRasterBandYSize(poPanBand))
1071 : {
1072 0 : CPLError(CE_Warning, CPLE_AppDefined,
1073 : "Dimensions of spectral band larger than panchro band");
1074 : }
1075 :
1076 72 : aMapDstBandToSpectralBandIter = aMapDstBandToSpectralBand.begin();
1077 242 : for (; aMapDstBandToSpectralBandIter != aMapDstBandToSpectralBand.end();
1078 170 : ++aMapDstBandToSpectralBandIter)
1079 : {
1080 171 : const int nDstBand = 1 + aMapDstBandToSpectralBandIter->first;
1081 341 : if (nDstBand > nBands ||
1082 170 : !static_cast<VRTRasterBand *>(GetRasterBand(nDstBand))
1083 170 : ->IsPansharpenRasterBand())
1084 : {
1085 1 : CPLError(CE_Failure, CPLE_AppDefined,
1086 : "SpectralBand.dstBand = '%d' invalid", nDstBand);
1087 1 : goto error;
1088 : }
1089 : }
1090 :
1091 71 : if (adfWeights.empty())
1092 : {
1093 160 : for (int i = 0; i < static_cast<int>(ahSpectralBands.size()); i++)
1094 : {
1095 115 : adfWeights.push_back(1.0 / ahSpectralBands.size());
1096 : }
1097 : }
1098 26 : else if (adfWeights.size() != ahSpectralBands.size())
1099 : {
1100 1 : CPLError(CE_Failure, CPLE_AppDefined,
1101 : "%d weights defined, but %d input spectral bands",
1102 1 : static_cast<int>(adfWeights.size()),
1103 1 : static_cast<int>(ahSpectralBands.size()));
1104 1 : goto error;
1105 : }
1106 :
1107 70 : if (aMapDstBandToSpectralBand.empty())
1108 : {
1109 1 : CPLError(CE_Warning, CPLE_AppDefined,
1110 : "No spectral band is mapped to an output band");
1111 : }
1112 :
1113 : /* -------------------------------------------------------------------- */
1114 : /* Instantiate poPansharpener */
1115 : /* -------------------------------------------------------------------- */
1116 70 : psPanOptions.reset(GDALCreatePansharpenOptions());
1117 70 : psPanOptions->ePansharpenAlg = GDAL_PSH_WEIGHTED_BROVEY;
1118 70 : psPanOptions->eResampleAlg = eResampleAlg;
1119 70 : psPanOptions->nBitDepth = nBitDepth;
1120 70 : psPanOptions->nWeightCount = static_cast<int>(adfWeights.size());
1121 140 : psPanOptions->padfWeights =
1122 70 : static_cast<double *>(CPLMalloc(sizeof(double) * adfWeights.size()));
1123 70 : memcpy(psPanOptions->padfWeights, &adfWeights[0],
1124 70 : sizeof(double) * adfWeights.size());
1125 70 : psPanOptions->hPanchroBand = poPanBand;
1126 70 : psPanOptions->nInputSpectralBands =
1127 70 : static_cast<int>(ahSpectralBands.size());
1128 140 : psPanOptions->pahInputSpectralBands = static_cast<GDALRasterBandH *>(
1129 70 : CPLMalloc(sizeof(GDALRasterBandH) * ahSpectralBands.size()));
1130 70 : memcpy(psPanOptions->pahInputSpectralBands, &ahSpectralBands[0],
1131 70 : sizeof(GDALRasterBandH) * ahSpectralBands.size());
1132 70 : psPanOptions->nOutPansharpenedBands =
1133 70 : static_cast<int>(aMapDstBandToSpectralBand.size());
1134 140 : psPanOptions->panOutPansharpenedBands = static_cast<int *>(
1135 70 : CPLMalloc(sizeof(int) * aMapDstBandToSpectralBand.size()));
1136 70 : aMapDstBandToSpectralBandIter = aMapDstBandToSpectralBand.begin();
1137 235 : for (int i = 0;
1138 235 : aMapDstBandToSpectralBandIter != aMapDstBandToSpectralBand.end();
1139 165 : ++aMapDstBandToSpectralBandIter, ++i)
1140 : {
1141 165 : psPanOptions->panOutPansharpenedBands[i] =
1142 165 : aMapDstBandToSpectralBandIter->second;
1143 : }
1144 70 : psPanOptions->bHasNoData = bHasNoData;
1145 70 : psPanOptions->dfNoData = dfNoData;
1146 70 : psPanOptions->nThreads = nThreads;
1147 :
1148 70 : if (nBands == psPanOptions->nOutPansharpenedBands)
1149 66 : SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
1150 :
1151 70 : m_poPansharpener = std::make_unique<GDALPansharpenOperation>();
1152 70 : eErr = m_poPansharpener->Initialize(psPanOptions.get());
1153 70 : if (eErr != CE_None)
1154 : {
1155 : // Delete the pansharper object before closing the dataset
1156 : // because it may have warped the bands into an intermediate VRT
1157 1 : m_poPansharpener.reset();
1158 :
1159 : // Close in reverse order (VRT firsts and real datasets after)
1160 4 : for (int i = static_cast<int>(m_apoDatasetsToReleaseRef.size()) - 1;
1161 4 : i >= 0; i--)
1162 : {
1163 3 : m_apoDatasetsToReleaseRef[i].reset();
1164 : }
1165 1 : m_apoDatasetsToReleaseRef.clear();
1166 : }
1167 :
1168 70 : return eErr;
1169 :
1170 16 : error:
1171 : // Close in reverse order (VRT firsts and real datasets after)
1172 66 : for (int i = static_cast<int>(m_apoDatasetsToReleaseRef.size()) - 1; i >= 0;
1173 : i--)
1174 : {
1175 50 : m_apoDatasetsToReleaseRef[i].reset();
1176 : }
1177 16 : m_apoDatasetsToReleaseRef.clear();
1178 16 : return CE_Failure;
1179 : }
1180 :
1181 : /************************************************************************/
1182 : /* SerializeToXML() */
1183 : /************************************************************************/
1184 :
1185 8 : CPLXMLNode *VRTPansharpenedDataset::SerializeToXML(const char *pszVRTPathIn)
1186 :
1187 : {
1188 8 : CPLXMLNode *psTree = VRTDataset::SerializeToXML(pszVRTPathIn);
1189 :
1190 8 : if (psTree == nullptr)
1191 0 : return psTree;
1192 :
1193 : /* -------------------------------------------------------------------- */
1194 : /* Set subclass. */
1195 : /* -------------------------------------------------------------------- */
1196 8 : CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
1197 : CXT_Text, "VRTPansharpenedDataset");
1198 :
1199 : /* -------------------------------------------------------------------- */
1200 : /* Serialize the block size. */
1201 : /* -------------------------------------------------------------------- */
1202 8 : CPLCreateXMLElementAndValue(psTree, "BlockXSize",
1203 : CPLSPrintf("%d", m_nBlockXSize));
1204 8 : CPLCreateXMLElementAndValue(psTree, "BlockYSize",
1205 : CPLSPrintf("%d", m_nBlockYSize));
1206 :
1207 : /* -------------------------------------------------------------------- */
1208 : /* Serialize the options. */
1209 : /* -------------------------------------------------------------------- */
1210 8 : if (m_poPansharpener == nullptr)
1211 0 : return psTree;
1212 8 : const GDALPansharpenOptions *psOptions = m_poPansharpener->GetOptions();
1213 8 : if (psOptions == nullptr)
1214 0 : return psTree;
1215 :
1216 : CPLXMLNode *psOptionsNode =
1217 8 : CPLCreateXMLNode(psTree, CXT_Element, "PansharpeningOptions");
1218 :
1219 8 : if (psOptions->ePansharpenAlg == GDAL_PSH_WEIGHTED_BROVEY)
1220 : {
1221 8 : CPLCreateXMLElementAndValue(psOptionsNode, "Algorithm",
1222 : "WeightedBrovey");
1223 : }
1224 : else
1225 : {
1226 0 : CPLAssert(false);
1227 : }
1228 8 : if (psOptions->nWeightCount)
1229 : {
1230 16 : CPLString osWeights;
1231 29 : for (int i = 0; i < psOptions->nWeightCount; i++)
1232 : {
1233 21 : if (i)
1234 13 : osWeights += ",";
1235 21 : osWeights += CPLSPrintf("%.16g", psOptions->padfWeights[i]);
1236 : }
1237 8 : CPLCreateXMLElementAndValue(
1238 : CPLCreateXMLNode(psOptionsNode, CXT_Element, "AlgorithmOptions"),
1239 : "Weights", osWeights.c_str());
1240 : }
1241 8 : CPLCreateXMLElementAndValue(
1242 : psOptionsNode, "Resampling",
1243 8 : GDALRasterIOGetResampleAlg(psOptions->eResampleAlg));
1244 :
1245 8 : if (psOptions->nThreads == -1)
1246 : {
1247 7 : CPLCreateXMLElementAndValue(psOptionsNode, "NumThreads", "ALL_CPUS");
1248 : }
1249 1 : else if (psOptions->nThreads > 1)
1250 : {
1251 0 : CPLCreateXMLElementAndValue(psOptionsNode, "NumThreads",
1252 0 : CPLSPrintf("%d", psOptions->nThreads));
1253 : }
1254 :
1255 8 : if (psOptions->nBitDepth)
1256 0 : CPLCreateXMLElementAndValue(psOptionsNode, "BitDepth",
1257 0 : CPLSPrintf("%d", psOptions->nBitDepth));
1258 :
1259 8 : const char *pszAdjust = nullptr;
1260 8 : switch (m_eGTAdjustment)
1261 : {
1262 8 : case GTAdjust_Union:
1263 8 : pszAdjust = "Union";
1264 8 : break;
1265 0 : case GTAdjust_Intersection:
1266 0 : pszAdjust = "Intersection";
1267 0 : break;
1268 0 : case GTAdjust_None:
1269 0 : pszAdjust = "None";
1270 0 : break;
1271 0 : case GTAdjust_NoneWithoutWarning:
1272 0 : pszAdjust = "NoneWithoutWarning";
1273 0 : break;
1274 0 : default:
1275 0 : break;
1276 : }
1277 :
1278 8 : if (psOptions->bHasNoData)
1279 : {
1280 1 : CPLCreateXMLElementAndValue(psOptionsNode, "NoData",
1281 1 : CPLSPrintf("%.16g", psOptions->dfNoData));
1282 : }
1283 7 : else if (m_bNoDataDisabled)
1284 : {
1285 0 : CPLCreateXMLElementAndValue(psOptionsNode, "NoData", "None");
1286 : }
1287 :
1288 8 : if (pszAdjust)
1289 8 : CPLCreateXMLElementAndValue(psOptionsNode, "SpatialExtentAdjustment",
1290 : pszAdjust);
1291 :
1292 8 : if (psOptions->hPanchroBand)
1293 : {
1294 : CPLXMLNode *psBand =
1295 8 : CPLCreateXMLNode(psOptionsNode, CXT_Element, "PanchroBand");
1296 : GDALRasterBand *poBand =
1297 8 : GDALRasterBand::FromHandle(psOptions->hPanchroBand);
1298 8 : if (poBand->GetDataset())
1299 : {
1300 8 : auto poPanchroDS = poBand->GetDataset();
1301 : std::map<CPLString, CPLString>::iterator oIter =
1302 8 : m_oMapToRelativeFilenames.find(poPanchroDS->GetDescription());
1303 8 : if (oIter == m_oMapToRelativeFilenames.end())
1304 : {
1305 : const char *pszUnderlyingDS =
1306 7 : poPanchroDS->GetMetadataItem("UNDERLYING_DATASET");
1307 7 : if (!pszUnderlyingDS)
1308 6 : pszUnderlyingDS = poPanchroDS->GetDescription();
1309 7 : CPLCreateXMLElementAndValue(psBand, "SourceFilename",
1310 : pszUnderlyingDS);
1311 : }
1312 : else
1313 : {
1314 1 : CPLXMLNode *psSourceFilename = CPLCreateXMLElementAndValue(
1315 1 : psBand, "SourceFilename", oIter->second);
1316 1 : CPLCreateXMLNode(CPLCreateXMLNode(psSourceFilename,
1317 : CXT_Attribute,
1318 : "relativeToVRT"),
1319 : CXT_Text, "1");
1320 : }
1321 :
1322 8 : GDALSerializeOpenOptionsToXML(psBand,
1323 8 : poPanchroDS->GetOpenOptions());
1324 :
1325 8 : CPLCreateXMLElementAndValue(psBand, "SourceBand",
1326 : CPLSPrintf("%d", poBand->GetBand()));
1327 : }
1328 : }
1329 29 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
1330 : {
1331 : CPLXMLNode *psBand =
1332 21 : CPLCreateXMLNode(psOptionsNode, CXT_Element, "SpectralBand");
1333 :
1334 40 : for (int j = 0; j < psOptions->nOutPansharpenedBands; j++)
1335 : {
1336 40 : if (psOptions->panOutPansharpenedBands[j] == i)
1337 : {
1338 41 : for (int k = 0; k < nBands; k++)
1339 : {
1340 41 : if (static_cast<VRTRasterBand *>(GetRasterBand(k + 1))
1341 41 : ->IsPansharpenRasterBand())
1342 : {
1343 40 : if (static_cast<VRTPansharpenedRasterBand *>(
1344 40 : GetRasterBand(k + 1))
1345 40 : ->GetIndexAsPansharpenedBand() == j)
1346 : {
1347 21 : CPLCreateXMLNode(CPLCreateXMLNode(psBand,
1348 : CXT_Attribute,
1349 : "dstBand"),
1350 : CXT_Text, CPLSPrintf("%d", k + 1));
1351 21 : break;
1352 : }
1353 : }
1354 : }
1355 21 : break;
1356 : }
1357 : }
1358 :
1359 : GDALRasterBand *poBand =
1360 21 : GDALRasterBand::FromHandle(psOptions->pahInputSpectralBands[i]);
1361 21 : if (poBand->GetDataset())
1362 : {
1363 21 : auto poSpectralDS = poBand->GetDataset();
1364 : std::map<CPLString, CPLString>::iterator oIter =
1365 21 : m_oMapToRelativeFilenames.find(poSpectralDS->GetDescription());
1366 21 : if (oIter == m_oMapToRelativeFilenames.end())
1367 : {
1368 : const char *pszUnderlyingDS =
1369 21 : poSpectralDS->GetMetadataItem("UNDERLYING_DATASET");
1370 21 : if (!pszUnderlyingDS)
1371 20 : pszUnderlyingDS = poSpectralDS->GetDescription();
1372 21 : CPLCreateXMLElementAndValue(psBand, "SourceFilename",
1373 : pszUnderlyingDS);
1374 : }
1375 : else
1376 : {
1377 0 : CPLXMLNode *psSourceFilename = CPLCreateXMLElementAndValue(
1378 0 : psBand, "SourceFilename", oIter->second);
1379 0 : CPLCreateXMLNode(CPLCreateXMLNode(psSourceFilename,
1380 : CXT_Attribute,
1381 : "relativeToVRT"),
1382 : CXT_Text, "1");
1383 : }
1384 :
1385 21 : GDALSerializeOpenOptionsToXML(psBand,
1386 21 : poSpectralDS->GetOpenOptions());
1387 :
1388 21 : CPLCreateXMLElementAndValue(psBand, "SourceBand",
1389 : CPLSPrintf("%d", poBand->GetBand()));
1390 : }
1391 : }
1392 :
1393 8 : return psTree;
1394 : }
1395 :
1396 : /************************************************************************/
1397 : /* GetBlockSize() */
1398 : /************************************************************************/
1399 :
1400 195 : void VRTPansharpenedDataset::GetBlockSize(int *pnBlockXSize,
1401 : int *pnBlockYSize) const
1402 :
1403 : {
1404 195 : assert(nullptr != pnBlockXSize);
1405 195 : assert(nullptr != pnBlockYSize);
1406 :
1407 195 : *pnBlockXSize = m_nBlockXSize;
1408 195 : *pnBlockYSize = m_nBlockYSize;
1409 195 : }
1410 :
1411 : /************************************************************************/
1412 : /* AddBand() */
1413 : /************************************************************************/
1414 :
1415 1 : CPLErr VRTPansharpenedDataset::AddBand(CPL_UNUSED GDALDataType eType,
1416 : CPL_UNUSED CSLConstList papszOptions)
1417 :
1418 : {
1419 1 : CPLError(CE_Failure, CPLE_NotSupported, "AddBand() not supported");
1420 :
1421 1 : return CE_Failure;
1422 : }
1423 :
1424 : /************************************************************************/
1425 : /* IRasterIO() */
1426 : /************************************************************************/
1427 :
1428 21 : CPLErr VRTPansharpenedDataset::IRasterIO(
1429 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
1430 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
1431 : int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
1432 : GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
1433 : {
1434 21 : if (eRWFlag == GF_Write)
1435 0 : return CE_Failure;
1436 :
1437 : /* Try to pass the request to the most appropriate overview dataset */
1438 21 : if (nBufXSize < nXSize && nBufYSize < nYSize)
1439 : {
1440 : int bTried;
1441 1 : CPLErr eErr = TryOverviewRasterIO(
1442 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1443 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
1444 : nBandSpace, psExtraArg, &bTried);
1445 1 : if (bTried)
1446 0 : return eErr;
1447 : }
1448 :
1449 21 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eBufType);
1450 21 : if (nXSize == nBufXSize && nYSize == nBufYSize &&
1451 20 : nDataTypeSize == nPixelSpace && nLineSpace == nPixelSpace * nBufXSize &&
1452 20 : nBandSpace == nLineSpace * nBufYSize && nBandCount == nBands)
1453 : {
1454 82 : for (int i = 0; i < nBands; i++)
1455 : {
1456 124 : if (panBandMap[i] != i + 1 ||
1457 62 : !static_cast<VRTRasterBand *>(GetRasterBand(i + 1))
1458 62 : ->IsPansharpenRasterBand())
1459 : {
1460 0 : goto default_path;
1461 : }
1462 : }
1463 :
1464 : //{static int bDone = 0; if (!bDone) printf("(2)\n"); bDone = 1; }
1465 20 : return m_poPansharpener->ProcessRegion(nXOff, nYOff, nXSize, nYSize,
1466 20 : pData, eBufType);
1467 : }
1468 :
1469 1 : default_path:
1470 : //{static int bDone = 0; if (!bDone) printf("(3)\n"); bDone = 1; }
1471 1 : return VRTDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
1472 : nBufXSize, nBufYSize, eBufType, nBandCount,
1473 : panBandMap, nPixelSpace, nLineSpace,
1474 1 : nBandSpace, psExtraArg);
1475 : }
1476 :
1477 : /************************************************************************/
1478 : /* ==================================================================== */
1479 : /* VRTPansharpenedRasterBand */
1480 : /* ==================================================================== */
1481 : /************************************************************************/
1482 :
1483 : /************************************************************************/
1484 : /* VRTPansharpenedRasterBand() */
1485 : /************************************************************************/
1486 :
1487 195 : VRTPansharpenedRasterBand::VRTPansharpenedRasterBand(GDALDataset *poDSIn,
1488 : int nBandIn,
1489 195 : GDALDataType eDataTypeIn)
1490 195 : : m_nIndexAsPansharpenedBand(nBandIn - 1)
1491 : {
1492 195 : Initialize(poDSIn->GetRasterXSize(), poDSIn->GetRasterYSize());
1493 :
1494 195 : poDS = poDSIn;
1495 195 : nBand = nBandIn;
1496 195 : eAccess = GA_Update;
1497 195 : eDataType = eDataTypeIn;
1498 :
1499 195 : static_cast<VRTPansharpenedDataset *>(poDS)->GetBlockSize(&nBlockXSize,
1500 : &nBlockYSize);
1501 195 : }
1502 :
1503 : /************************************************************************/
1504 : /* ~VRTPansharpenedRasterBand() */
1505 : /************************************************************************/
1506 :
1507 390 : VRTPansharpenedRasterBand::~VRTPansharpenedRasterBand()
1508 :
1509 : {
1510 195 : FlushCache(true);
1511 390 : }
1512 :
1513 : /************************************************************************/
1514 : /* IReadBlock() */
1515 : /************************************************************************/
1516 :
1517 26 : CPLErr VRTPansharpenedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
1518 : void *pImage)
1519 :
1520 : {
1521 26 : const int nReqXOff = nBlockXOff * nBlockXSize;
1522 26 : const int nReqYOff = nBlockYOff * nBlockYSize;
1523 26 : const int nReqXSize = std::min(nBlockXSize, nRasterXSize - nReqXOff);
1524 26 : const int nReqYSize = std::min(nBlockYSize, nRasterYSize - nReqYOff);
1525 :
1526 : //{static int bDone = 0; if (!bDone) printf("(4)\n"); bDone = 1; }
1527 26 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eDataType);
1528 : GDALRasterIOExtraArg sExtraArg;
1529 26 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1530 52 : if (IRasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize, pImage,
1531 : nReqXSize, nReqYSize, eDataType, nDataTypeSize,
1532 26 : static_cast<GSpacing>(nDataTypeSize) * nReqXSize,
1533 26 : &sExtraArg) != CE_None)
1534 : {
1535 0 : return CE_Failure;
1536 : }
1537 :
1538 26 : if (nReqXSize < nBlockXSize)
1539 : {
1540 5613 : for (int j = nReqYSize - 1; j >= 0; j--)
1541 : {
1542 5600 : memmove(static_cast<GByte *>(pImage) +
1543 5600 : static_cast<size_t>(j) * nDataTypeSize * nBlockXSize,
1544 5600 : static_cast<GByte *>(pImage) +
1545 5600 : static_cast<size_t>(j) * nDataTypeSize * nReqXSize,
1546 5600 : static_cast<size_t>(nReqXSize) * nDataTypeSize);
1547 5600 : memset(static_cast<GByte *>(pImage) +
1548 5600 : (static_cast<size_t>(j) * nBlockXSize + nReqXSize) *
1549 5600 : nDataTypeSize,
1550 : 0,
1551 5600 : static_cast<size_t>(nBlockXSize - nReqXSize) *
1552 5600 : nDataTypeSize);
1553 : }
1554 : }
1555 26 : if (nReqYSize < nBlockYSize)
1556 : {
1557 22 : memset(static_cast<GByte *>(pImage) +
1558 22 : static_cast<size_t>(nReqYSize) * nBlockXSize * nDataTypeSize,
1559 : 0,
1560 22 : static_cast<size_t>(nBlockYSize - nReqYSize) * nBlockXSize *
1561 22 : nDataTypeSize);
1562 : }
1563 :
1564 : // Cache other bands
1565 26 : CPLErr eErr = CE_None;
1566 26 : VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);
1567 26 : if (poGDS->nBands != 1 && !poGDS->m_bLoadingOtherBands)
1568 : {
1569 10 : poGDS->m_bLoadingOtherBands = TRUE;
1570 :
1571 36 : for (int iOtherBand = 1; iOtherBand <= poGDS->nBands; iOtherBand++)
1572 : {
1573 26 : if (iOtherBand == nBand)
1574 10 : continue;
1575 :
1576 : GDALRasterBlock *poBlock =
1577 16 : poGDS->GetRasterBand(iOtherBand)
1578 16 : ->GetLockedBlockRef(nBlockXOff, nBlockYOff);
1579 16 : if (poBlock == nullptr)
1580 : {
1581 0 : eErr = CE_Failure;
1582 0 : break;
1583 : }
1584 16 : poBlock->DropLock();
1585 : }
1586 :
1587 10 : poGDS->m_bLoadingOtherBands = FALSE;
1588 : }
1589 :
1590 26 : return eErr;
1591 : }
1592 :
1593 : /************************************************************************/
1594 : /* IRasterIO() */
1595 : /************************************************************************/
1596 :
1597 127 : CPLErr VRTPansharpenedRasterBand::IRasterIO(
1598 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
1599 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
1600 : GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
1601 : {
1602 127 : if (eRWFlag == GF_Write)
1603 0 : return CE_Failure;
1604 :
1605 127 : VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);
1606 :
1607 : /* Try to pass the request to the most appropriate overview dataset */
1608 127 : if (nBufXSize < nXSize && nBufYSize < nYSize)
1609 : {
1610 : int bTried;
1611 3 : CPLErr eErr = TryOverviewRasterIO(
1612 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1613 : eBufType, nPixelSpace, nLineSpace, psExtraArg, &bTried);
1614 3 : if (bTried)
1615 0 : return eErr;
1616 : }
1617 :
1618 127 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eBufType);
1619 127 : if (nDataTypeSize > 0 && nXSize == nBufXSize && nYSize == nBufYSize &&
1620 124 : nDataTypeSize == nPixelSpace && nLineSpace == nPixelSpace * nBufXSize)
1621 : {
1622 : const GDALPansharpenOptions *psOptions =
1623 124 : poGDS->m_poPansharpener->GetOptions();
1624 :
1625 : // Have we already done this request for another band ?
1626 : // If so use the cached result
1627 124 : const size_t nBufferSizePerBand =
1628 124 : static_cast<size_t>(nXSize) * nYSize * nDataTypeSize;
1629 124 : if (nXOff == poGDS->m_nLastBandRasterIOXOff &&
1630 118 : nYOff >= poGDS->m_nLastBandRasterIOYOff &&
1631 111 : nXSize == poGDS->m_nLastBandRasterIOXSize &&
1632 69 : nYOff + nYSize <= poGDS->m_nLastBandRasterIOYOff +
1633 69 : poGDS->m_nLastBandRasterIOYSize &&
1634 59 : eBufType == poGDS->m_eLastBandRasterIODataType)
1635 : {
1636 : //{static int bDone = 0; if (!bDone) printf("(6)\n"); bDone = 1; }
1637 53 : if (poGDS->m_pabyLastBufferBandRasterIO == nullptr)
1638 0 : return CE_Failure;
1639 53 : const size_t nBufferSizePerBandCached =
1640 53 : static_cast<size_t>(nXSize) * poGDS->m_nLastBandRasterIOYSize *
1641 53 : nDataTypeSize;
1642 53 : memcpy(pData,
1643 53 : poGDS->m_pabyLastBufferBandRasterIO +
1644 53 : nBufferSizePerBandCached * m_nIndexAsPansharpenedBand +
1645 53 : static_cast<size_t>(nYOff -
1646 53 : poGDS->m_nLastBandRasterIOYOff) *
1647 53 : nXSize * nDataTypeSize,
1648 : nBufferSizePerBand);
1649 53 : return CE_None;
1650 : }
1651 :
1652 71 : int nYSizeToCache = nYSize;
1653 71 : if (nYSize == 1 && nXSize == nRasterXSize)
1654 : {
1655 : //{static int bDone = 0; if (!bDone) printf("(7)\n"); bDone = 1; }
1656 : // For efficiency, try to cache at leak 256 K
1657 0 : nYSizeToCache = (256 * 1024) / nXSize / nDataTypeSize;
1658 0 : if (nYSizeToCache == 0)
1659 0 : nYSizeToCache = 1;
1660 0 : else if (nYOff + nYSizeToCache > nRasterYSize)
1661 0 : nYSizeToCache = nRasterYSize - nYOff;
1662 : }
1663 71 : const GUIntBig nBufferSize = static_cast<GUIntBig>(nXSize) *
1664 71 : nYSizeToCache * nDataTypeSize *
1665 71 : psOptions->nOutPansharpenedBands;
1666 : // Check the we don't overflow (for 32 bit platforms)
1667 71 : if (nBufferSize > std::numeric_limits<size_t>::max() / 2)
1668 : {
1669 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
1670 : "Out of memory error while allocating working buffers");
1671 0 : return CE_Failure;
1672 : }
1673 : GByte *pabyTemp = static_cast<GByte *>(
1674 71 : VSI_REALLOC_VERBOSE(poGDS->m_pabyLastBufferBandRasterIO,
1675 : static_cast<size_t>(nBufferSize)));
1676 71 : if (pabyTemp == nullptr)
1677 : {
1678 0 : return CE_Failure;
1679 : }
1680 71 : poGDS->m_nLastBandRasterIOXOff = nXOff;
1681 71 : poGDS->m_nLastBandRasterIOYOff = nYOff;
1682 71 : poGDS->m_nLastBandRasterIOXSize = nXSize;
1683 71 : poGDS->m_nLastBandRasterIOYSize = nYSizeToCache;
1684 71 : poGDS->m_eLastBandRasterIODataType = eBufType;
1685 71 : poGDS->m_pabyLastBufferBandRasterIO = pabyTemp;
1686 :
1687 71 : CPLErr eErr = poGDS->m_poPansharpener->ProcessRegion(
1688 : nXOff, nYOff, nXSize, nYSizeToCache,
1689 71 : poGDS->m_pabyLastBufferBandRasterIO, eBufType);
1690 71 : if (eErr == CE_None)
1691 : {
1692 : //{static int bDone = 0; if (!bDone) printf("(8)\n"); bDone = 1; }
1693 71 : size_t nBufferSizePerBandCached = static_cast<size_t>(nXSize) *
1694 71 : poGDS->m_nLastBandRasterIOYSize *
1695 71 : nDataTypeSize;
1696 71 : memcpy(pData,
1697 71 : poGDS->m_pabyLastBufferBandRasterIO +
1698 71 : nBufferSizePerBandCached * m_nIndexAsPansharpenedBand,
1699 : nBufferSizePerBand);
1700 : }
1701 : else
1702 : {
1703 0 : VSIFree(poGDS->m_pabyLastBufferBandRasterIO);
1704 0 : poGDS->m_pabyLastBufferBandRasterIO = nullptr;
1705 : }
1706 71 : return eErr;
1707 : }
1708 :
1709 : //{static int bDone = 0; if (!bDone) printf("(9)\n"); bDone = 1; }
1710 3 : CPLErr eErr = VRTRasterBand::IRasterIO(
1711 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1712 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
1713 :
1714 3 : return eErr;
1715 : }
1716 :
1717 : /************************************************************************/
1718 : /* SerializeToXML() */
1719 : /************************************************************************/
1720 :
1721 : CPLXMLNode *
1722 21 : VRTPansharpenedRasterBand::SerializeToXML(const char *pszVRTPathIn,
1723 : bool &bHasWarnedAboutRAMUsage,
1724 : size_t &nAccRAMUsage)
1725 :
1726 : {
1727 21 : CPLXMLNode *psTree = VRTRasterBand::SerializeToXML(
1728 : pszVRTPathIn, bHasWarnedAboutRAMUsage, nAccRAMUsage);
1729 :
1730 : /* -------------------------------------------------------------------- */
1731 : /* Set subclass. */
1732 : /* -------------------------------------------------------------------- */
1733 21 : CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
1734 : CXT_Text, "VRTPansharpenedRasterBand");
1735 :
1736 21 : return psTree;
1737 : }
1738 :
1739 : /************************************************************************/
1740 : /* GetOverviewCount() */
1741 : /************************************************************************/
1742 :
1743 25 : int VRTPansharpenedRasterBand::GetOverviewCount()
1744 : {
1745 25 : VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);
1746 :
1747 : // Build on-the-fly overviews from overviews of pan and spectral bands
1748 50 : if (poGDS->m_poPansharpener != nullptr &&
1749 50 : poGDS->m_apoOverviewDatasets.empty() && poGDS->m_poMainDataset == poGDS)
1750 : {
1751 : const GDALPansharpenOptions *psOptions =
1752 20 : poGDS->m_poPansharpener->GetOptions();
1753 :
1754 : GDALRasterBand *poPanBand =
1755 20 : GDALRasterBand::FromHandle(psOptions->hPanchroBand);
1756 20 : int nOvrCount = poPanBand->GetOverviewCount();
1757 20 : if (nOvrCount > 0)
1758 : {
1759 14 : for (int i = 0; i < poGDS->GetRasterCount(); i++)
1760 : {
1761 10 : if (!static_cast<VRTRasterBand *>(poGDS->GetRasterBand(i + 1))
1762 10 : ->IsPansharpenRasterBand())
1763 : {
1764 0 : return 0;
1765 : }
1766 : }
1767 :
1768 : // Limit number of overviews of the VRTPansharpenedRasterBand to
1769 : // the minimum number of overviews of the pan and spectral source
1770 : // bands, and also make sure all spectral bands have overviews
1771 : // of the same dimension for a given level.
1772 4 : std::vector<std::pair<int, int>> sizeSpectralOverviews;
1773 14 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
1774 : {
1775 10 : auto poSpectralBand = GDALRasterBand::FromHandle(
1776 10 : psOptions->pahInputSpectralBands[i]);
1777 10 : nOvrCount =
1778 10 : std::min(nOvrCount, poSpectralBand->GetOverviewCount());
1779 10 : if (i == 0)
1780 : {
1781 7 : for (int iOvr = 0; iOvr < nOvrCount; ++iOvr)
1782 : {
1783 3 : auto poOvrBand = poSpectralBand->GetOverview(iOvr);
1784 : sizeSpectralOverviews.emplace_back(
1785 3 : poOvrBand->GetXSize(), poOvrBand->GetYSize());
1786 : }
1787 : }
1788 : else
1789 : {
1790 10 : for (int iOvr = 0; iOvr < nOvrCount; ++iOvr)
1791 : {
1792 4 : auto poOvrBand = poSpectralBand->GetOverview(iOvr);
1793 4 : if (sizeSpectralOverviews[iOvr].first !=
1794 8 : poOvrBand->GetXSize() ||
1795 4 : sizeSpectralOverviews[iOvr].second !=
1796 4 : poOvrBand->GetYSize())
1797 : {
1798 0 : nOvrCount = iOvr;
1799 0 : break;
1800 : }
1801 : }
1802 : }
1803 : }
1804 :
1805 4 : auto poPanBandDS = poPanBand->GetDataset();
1806 7 : for (int j = 0; j < nOvrCount; j++)
1807 : {
1808 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>
1809 3 : poPanOvrDS(GDALCreateOverviewDataset(poPanBandDS, j, true));
1810 3 : if (!poPanOvrDS)
1811 : {
1812 0 : CPLError(CE_Warning, CPLE_AppDefined,
1813 : "GDALCreateOverviewDataset(poPanBandDS, %d, true) "
1814 : "failed",
1815 : j);
1816 0 : break;
1817 : }
1818 : GDALRasterBand *poPanOvrBand =
1819 3 : poPanOvrDS->GetRasterBand(poPanBand->GetBand());
1820 : auto poOvrDS = std::make_unique<VRTPansharpenedDataset>(
1821 3 : poPanOvrBand->GetXSize(), poPanOvrBand->GetYSize());
1822 6 : poOvrDS->m_apoDatasetsToReleaseRef.push_back(
1823 3 : std::move(poPanOvrDS));
1824 10 : for (int i = 0; i < poGDS->GetRasterCount(); i++)
1825 : {
1826 7 : GDALRasterBand *poSrcBand = poGDS->GetRasterBand(i + 1);
1827 : auto poBand = std::make_unique<VRTPansharpenedRasterBand>(
1828 7 : poOvrDS.get(), i + 1, poSrcBand->GetRasterDataType());
1829 : const char *pszNBITS =
1830 7 : poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
1831 7 : if (pszNBITS)
1832 0 : poBand->SetMetadataItem("NBITS", pszNBITS,
1833 0 : "IMAGE_STRUCTURE");
1834 7 : int bHasNoData = FALSE;
1835 : const double dfNoData =
1836 7 : poSrcBand->GetNoDataValue(&bHasNoData);
1837 7 : if (bHasNoData)
1838 1 : poBand->SetNoDataValue(dfNoData);
1839 14 : poBand->SetColorInterpretation(
1840 7 : poSrcBand->GetColorInterpretation());
1841 7 : poOvrDS->SetBand(i + 1, std::move(poBand));
1842 : }
1843 :
1844 : std::unique_ptr<GDALPansharpenOptions,
1845 : decltype(&GDALDestroyPansharpenOptions)>
1846 : psPanOvrOptions(GDALClonePansharpenOptions(psOptions),
1847 3 : GDALDestroyPansharpenOptions);
1848 3 : psPanOvrOptions->hPanchroBand = poPanOvrBand;
1849 10 : for (int i = 0; i < psOptions->nInputSpectralBands; i++)
1850 : {
1851 7 : auto poSpectralBand = GDALRasterBand::FromHandle(
1852 7 : psOptions->pahInputSpectralBands[i]);
1853 : std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>
1854 : poSpectralOvrDS(GDALCreateOverviewDataset(
1855 7 : poSpectralBand->GetDataset(), j, true));
1856 7 : if (!poSpectralOvrDS)
1857 : {
1858 0 : CPLError(CE_Warning, CPLE_AppDefined,
1859 : "GDALCreateOverviewDataset(poSpectralBand->"
1860 : "GetDataset(), %d, true) failed",
1861 : j);
1862 0 : return 0;
1863 : }
1864 14 : psPanOvrOptions->pahInputSpectralBands[i] =
1865 7 : poSpectralOvrDS->GetRasterBand(
1866 : poSpectralBand->GetBand());
1867 14 : poOvrDS->m_apoDatasetsToReleaseRef.push_back(
1868 7 : std::move(poSpectralOvrDS));
1869 : }
1870 3 : poOvrDS->m_poPansharpener =
1871 6 : std::make_unique<GDALPansharpenOperation>();
1872 6 : if (poOvrDS->m_poPansharpener->Initialize(
1873 6 : psPanOvrOptions.get()) != CE_None)
1874 : {
1875 0 : CPLError(CE_Warning, CPLE_AppDefined,
1876 : "Unable to initialize pansharpener.");
1877 0 : return 0;
1878 : }
1879 :
1880 3 : poOvrDS->m_poMainDataset = poGDS;
1881 3 : poOvrDS->SetMetadataItem("INTERLEAVE", "PIXEL",
1882 3 : "IMAGE_STRUCTURE");
1883 :
1884 3 : poGDS->m_apoOverviewDatasets.push_back(std::move(poOvrDS));
1885 : }
1886 : }
1887 : }
1888 25 : return static_cast<int>(poGDS->m_apoOverviewDatasets.size());
1889 : }
1890 :
1891 : /************************************************************************/
1892 : /* GetOverview() */
1893 : /************************************************************************/
1894 :
1895 7 : GDALRasterBand *VRTPansharpenedRasterBand::GetOverview(int iOvr)
1896 : {
1897 7 : if (iOvr < 0 || iOvr >= GetOverviewCount())
1898 2 : return nullptr;
1899 :
1900 5 : VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);
1901 :
1902 5 : return poGDS->m_apoOverviewDatasets[iOvr]->GetRasterBand(nBand);
1903 : }
1904 :
1905 : /*! @endcond */
|
