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