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