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