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