/****************************************************************************** * * Project: Virtual GDAL Datasets * Purpose: Implementation of VRTPansharpenedRasterBand and *VRTPansharpenedDataset. Author: Even Rouault * ****************************************************************************** * Copyright (c) 2015, Even Rouault * * SPDX-License-Identifier: MIT ****************************************************************************/ #include "cpl_port.h" #include "gdal_vrt.h" #include "vrtdataset.h" #include #include #include #include #include #include #include #include #include #include #include #include #include "cpl_conv.h" #include "cpl_error.h" #include "cpl_minixml.h" #include "cpl_string.h" #include "cpl_vsi.h" #include "gdal.h" #include "gdal_priv.h" #include "gdalpansharpen.h" #include "ogr_core.h" #include "ogr_spatialref.h" /************************************************************************/ /* GDALCreatePansharpenedVRT() */ /************************************************************************/ /** * Create a virtual pansharpened dataset. * * This function will create a virtual pansharpened dataset. * * Starting with GDAL 3.12, a reference will be taken on the dataset owing * each input bands. * Before 3.12, no reference were taken on the passed bands. Consequently, * they or their dataset to which they belong had to be kept open until * this virtual pansharpened dataset is closed. * * The returned dataset will have no associated filename for itself. If you * want to write the virtual dataset description to a file, use the * GDALSetDescription() function (or SetDescription() method) on the dataset * to assign a filename before it is closed. * * @param pszXML Pansharpened VRT XML where <SpectralBand> elements have * no explicit SourceFilename and SourceBand. The spectral bands in the XML will * be assigned the successive values of the pahInputSpectralBands array. Must * not be NULL. * * @param hPanchroBand Panchromatic band. Must not be NULL. * * @param nInputSpectralBands Number of input spectral bands. Must be greater * than zero. * * @param pahInputSpectralBands Array of nInputSpectralBands spectral bands. * * @return NULL on failure, or a new virtual dataset handle on success to be * closed with GDALClose(). * * @since GDAL 2.1 */ GDALDatasetH GDALCreatePansharpenedVRT(const char *pszXML, GDALRasterBandH hPanchroBand, int nInputSpectralBands, GDALRasterBandH *pahInputSpectralBands) { VALIDATE_POINTER1(pszXML, "GDALCreatePansharpenedVRT", nullptr); VALIDATE_POINTER1(hPanchroBand, "GDALCreatePansharpenedVRT", nullptr); VALIDATE_POINTER1(pahInputSpectralBands, "GDALCreatePansharpenedVRT", nullptr); CPLXMLNode *psTree = CPLParseXMLString(pszXML); if (psTree == nullptr) return nullptr; VRTPansharpenedDataset *poDS = new VRTPansharpenedDataset(0, 0); CPLErr eErr = poDS->XMLInit(psTree, nullptr, hPanchroBand, nInputSpectralBands, pahInputSpectralBands); CPLDestroyXMLNode(psTree); if (eErr != CE_None) { delete poDS; return nullptr; } return GDALDataset::ToHandle(poDS); } /*! @cond Doxygen_Suppress */ /************************************************************************/ /* ==================================================================== */ /* VRTPansharpenedDataset */ /* ==================================================================== */ /************************************************************************/ /************************************************************************/ /* VRTPansharpenedDataset() */ /************************************************************************/ VRTPansharpenedDataset::VRTPansharpenedDataset(int nXSize, int nYSize, int nBlockXSize, int nBlockYSize) : VRTDataset(nXSize, nYSize, nBlockXSize > 0 ? nBlockXSize : std::min(nXSize, 512), nBlockYSize > 0 ? nBlockYSize : std::min(nYSize, 512)), m_poMainDataset(nullptr), m_bLoadingOtherBands(FALSE), m_pabyLastBufferBandRasterIO(nullptr), m_nLastBandRasterIOXOff(0), m_nLastBandRasterIOYOff(0), m_nLastBandRasterIOXSize(0), m_nLastBandRasterIOYSize(0), m_eLastBandRasterIODataType(GDT_Unknown), m_eGTAdjustment(GTAdjust_Union), m_bNoDataDisabled(FALSE) { eAccess = GA_Update; m_poMainDataset = this; } /************************************************************************/ /* ~VRTPansharpenedDataset() */ /************************************************************************/ VRTPansharpenedDataset::~VRTPansharpenedDataset() { VRTPansharpenedDataset::FlushCache(true); VRTPansharpenedDataset::CloseDependentDatasets(); CPLFree(m_pabyLastBufferBandRasterIO); } /************************************************************************/ /* CloseDependentDatasets() */ /************************************************************************/ int VRTPansharpenedDataset::CloseDependentDatasets() { int bHasDroppedRef = VRTDataset::CloseDependentDatasets(); /* -------------------------------------------------------------------- */ /* Destroy the raster bands if they exist. */ /* -------------------------------------------------------------------- */ for (int iBand = 0; iBand < nBands; iBand++) { delete papoBands[iBand]; } nBands = 0; // Destroy the overviews before m_poPansharpener as they might reference // files that are in m_apoDatasetsToReleaseRef. bHasDroppedRef |= !m_apoOverviewDatasets.empty(); m_apoOverviewDatasets.clear(); if (m_poPansharpener != nullptr) { // Delete the pansharper object before closing the dataset // because it may have warped the bands into an intermediate VRT m_poPansharpener.reset(); // Close in reverse order (VRT firsts and real datasets after) for (int i = static_cast(m_apoDatasetsToReleaseRef.size()) - 1; i >= 0; i--) { bHasDroppedRef = TRUE; m_apoDatasetsToReleaseRef[i].reset(); } m_apoDatasetsToReleaseRef.clear(); } return bHasDroppedRef; } /************************************************************************/ /* GetFileList() */ /************************************************************************/ char **VRTPansharpenedDataset::GetFileList() { char **papszFileList = GDALDataset::GetFileList(); if (m_poPansharpener != nullptr) { const GDALPansharpenOptions *psOptions = m_poPansharpener->GetOptions(); if (psOptions != nullptr) { std::set oSetNames; if (psOptions->hPanchroBand != nullptr) { GDALDataset *poDS = GDALDataset::FromHandle( GDALGetBandDataset(psOptions->hPanchroBand)); if (poDS != nullptr) { const char *pszName = poDS->GetMetadataItem("UNDERLYING_DATASET"); if (!pszName) pszName = poDS->GetDescription(); papszFileList = CSLAddString(papszFileList, pszName); oSetNames.insert(pszName); } } for (int i = 0; i < psOptions->nInputSpectralBands; i++) { if (psOptions->pahInputSpectralBands[i] != nullptr) { GDALDataset *poDS = GDALDataset::FromHandle(GDALGetBandDataset( psOptions->pahInputSpectralBands[i])); if (poDS != nullptr) { const char *pszName = poDS->GetMetadataItem("UNDERLYING_DATASET"); if (!pszName) pszName = poDS->GetDescription(); if (!cpl::contains(oSetNames, pszName)) { papszFileList = CSLAddString(papszFileList, pszName); oSetNames.insert(pszName); } } } } } } return papszFileList; } /************************************************************************/ /* XMLInit() */ /************************************************************************/ CPLErr VRTPansharpenedDataset::XMLInit(const CPLXMLNode *psTree, const char *pszVRTPathIn) { return XMLInit(psTree, pszVRTPathIn, nullptr, 0, nullptr); } CPLErr VRTPansharpenedDataset::XMLInit(const CPLXMLNode *psTree, const char *pszVRTPathIn, GDALRasterBandH hPanchroBandIn, int nInputSpectralBandsIn, GDALRasterBandH *pahInputSpectralBandsIn) { CPLErr eErr; std::unique_ptr psPanOptions(nullptr, GDALDestroyPansharpenOptions); /* -------------------------------------------------------------------- */ /* Initialize blocksize before calling sub-init so that the */ /* band initializers can get it from the dataset object when */ /* they are created. */ /* -------------------------------------------------------------------- */ m_nBlockXSize = atoi(CPLGetXMLValue(psTree, "BlockXSize", "512")); m_nBlockYSize = atoi(CPLGetXMLValue(psTree, "BlockYSize", "512")); /* -------------------------------------------------------------------- */ /* Parse PansharpeningOptions */ /* -------------------------------------------------------------------- */ const CPLXMLNode *psOptions = CPLGetXMLNode(psTree, "PansharpeningOptions"); if (psOptions == nullptr) { CPLError(CE_Failure, CPLE_AppDefined, "Missing PansharpeningOptions"); return CE_Failure; } if (CPLGetXMLValue(psOptions, "MSShiftX", nullptr) != nullptr || CPLGetXMLValue(psOptions, "MSShiftY", nullptr) != nullptr) { CPLError(CE_Failure, CPLE_AppDefined, "Options MSShiftX and MSShiftY are no longer supported. " "Spatial adjustment between panchromatic and multispectral " "bands is done through their geotransform."); return CE_Failure; } CPLString osSourceFilename; GDALDataset *poPanDataset = nullptr; GDALRasterBand *poPanBand = nullptr; std::map oMapNamesToDataset; int nPanBand; if (hPanchroBandIn == nullptr) { const CPLXMLNode *psPanchroBand = CPLGetXMLNode(psOptions, "PanchroBand"); if (psPanchroBand == nullptr) { CPLError(CE_Failure, CPLE_AppDefined, "PanchroBand missing"); return CE_Failure; } osSourceFilename = CPLGetXMLValue(psPanchroBand, "SourceFilename", ""); if (osSourceFilename.empty()) { CPLError(CE_Failure, CPLE_AppDefined, "PanchroBand.SourceFilename missing"); return CE_Failure; } const bool bRelativeToVRT = CPL_TO_BOOL(atoi(CPLGetXMLValue( psPanchroBand, "SourceFilename.relativetoVRT", "0"))); if (bRelativeToVRT) { const std::string osAbs = CPLProjectRelativeFilenameSafe( pszVRTPathIn, osSourceFilename.c_str()); m_oMapToRelativeFilenames[osAbs] = osSourceFilename; osSourceFilename = osAbs; } const CPLStringList aosOpenOptions( GDALDeserializeOpenOptionsFromXML(psPanchroBand)); poPanDataset = GDALDataset::Open(osSourceFilename, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr, aosOpenOptions.List()); if (poPanDataset == nullptr) { return CE_Failure; } const char *pszSourceBand = CPLGetXMLValue(psPanchroBand, "SourceBand", "1"); nPanBand = atoi(pszSourceBand); if (poPanBand == nullptr) poPanBand = poPanDataset->GetRasterBand(nPanBand); if (poPanBand == nullptr) { CPLError(CE_Failure, CPLE_AppDefined, "%s invalid band of %s", pszSourceBand, osSourceFilename.c_str()); GDALClose(poPanDataset); return CE_Failure; } oMapNamesToDataset[osSourceFilename] = poPanDataset; } else { poPanBand = GDALRasterBand::FromHandle(hPanchroBandIn); nPanBand = poPanBand->GetBand(); poPanDataset = poPanBand->GetDataset(); if (poPanDataset == nullptr) { CPLError( CE_Failure, CPLE_AppDefined, "Cannot retrieve dataset associated with panchromatic band"); return CE_Failure; } oMapNamesToDataset[CPLSPrintf("%p", poPanDataset)] = poPanDataset; poPanDataset->Reference(); } m_apoDatasetsToReleaseRef.push_back( std::unique_ptr( poPanDataset)); poPanDataset = m_apoDatasetsToReleaseRef.back().get(); // Figure out which kind of adjustment we should do if the pan and spectral // bands do not share the same geotransform const char *pszGTAdjustment = CPLGetXMLValue(psOptions, "SpatialExtentAdjustment", "Union"); if (EQUAL(pszGTAdjustment, "Union")) m_eGTAdjustment = GTAdjust_Union; else if (EQUAL(pszGTAdjustment, "Intersection")) m_eGTAdjustment = GTAdjust_Intersection; else if (EQUAL(pszGTAdjustment, "None")) m_eGTAdjustment = GTAdjust_None; else if (EQUAL(pszGTAdjustment, "NoneWithoutWarning")) m_eGTAdjustment = GTAdjust_NoneWithoutWarning; else { m_eGTAdjustment = GTAdjust_Union; CPLError(CE_Failure, CPLE_AppDefined, "Unsupported value for GeoTransformAdjustment. Defaulting to " "Union"); } const char *pszNumThreads = CPLGetXMLValue(psOptions, "NumThreads", nullptr); int nThreads = 0; if (pszNumThreads != nullptr) { if (EQUAL(pszNumThreads, "ALL_CPUS")) nThreads = -1; else nThreads = atoi(pszNumThreads); } const char *pszAlgorithm = CPLGetXMLValue(psOptions, "Algorithm", "WeightedBrovey"); if (!EQUAL(pszAlgorithm, "WeightedBrovey")) { CPLError(CE_Failure, CPLE_AppDefined, "Algorithm %s unsupported", pszAlgorithm); return CE_Failure; } std::vector adfWeights; if (const CPLXMLNode *psAlgOptions = CPLGetXMLNode(psOptions, "AlgorithmOptions")) { const char *pszWeights = CPLGetXMLValue(psAlgOptions, "Weights", nullptr); if (pszWeights != nullptr) { char **papszTokens = CSLTokenizeString2(pszWeights, " ,", 0); for (int i = 0; papszTokens && papszTokens[i]; i++) adfWeights.push_back(CPLAtof(papszTokens[i])); CSLDestroy(papszTokens); } } GDALRIOResampleAlg eResampleAlg = GDALRasterIOGetResampleAlg( CPLGetXMLValue(psOptions, "Resampling", "Cubic")); std::vector ahSpectralBands; std::map aMapDstBandToSpectralBand; std::map::iterator aMapDstBandToSpectralBandIter; int nBitDepth = 0; bool bFoundNonMatchingGT = false; GDALGeoTransform panGT; const bool bPanGeoTransformValid = (poPanDataset->GetGeoTransform(panGT) == CE_None); if (!bPanGeoTransformValid) { CPLError(CE_Failure, CPLE_AppDefined, "Panchromatic band has no associated geotransform"); return CE_Failure; } int nPanXSize = poPanBand->GetXSize(); int nPanYSize = poPanBand->GetYSize(); int bHasNoData = FALSE; double dfNoData = poPanBand->GetNoDataValue(&bHasNoData); double dfLRPanX = panGT[0] + nPanXSize * panGT[1] + nPanYSize * panGT[2]; double dfLRPanY = panGT[3] + nPanXSize * panGT[4] + nPanYSize * panGT[5]; bool bFoundRotatingTerms = (panGT[2] != 0.0 || panGT[4] != 0.0); double dfMinX = panGT[0]; double dfMaxX = dfLRPanX; double dfMaxY = panGT[3]; double dfMinY = dfLRPanY; if (dfMinY > dfMaxY) std::swap(dfMinY, dfMaxY); const double dfPanMinY = dfMinY; const double dfPanMaxY = dfMaxY; const auto poPanSRS = poPanDataset->GetSpatialRef(); /* -------------------------------------------------------------------- */ /* First pass on spectral datasets to check their georeferencing. */ /* -------------------------------------------------------------------- */ int iSpectralBand = 0; for (const CPLXMLNode *psIter = psOptions->psChild; psIter; psIter = psIter->psNext) { GDALDataset *poDataset; if (psIter->eType != CXT_Element || !EQUAL(psIter->pszValue, "SpectralBand")) continue; if (nInputSpectralBandsIn && pahInputSpectralBandsIn != nullptr) { if (iSpectralBand == nInputSpectralBandsIn) { CPLError(CE_Failure, CPLE_AppDefined, "More SpectralBand elements than in source array"); goto error; } poDataset = GDALRasterBand::FromHandle( pahInputSpectralBandsIn[iSpectralBand]) ->GetDataset(); if (poDataset == nullptr) { CPLError(CE_Failure, CPLE_AppDefined, "SpectralBand has no associated dataset"); goto error; } osSourceFilename = poDataset->GetDescription(); oMapNamesToDataset[CPLSPrintf("%p", poDataset)] = poDataset; poDataset->Reference(); } else { osSourceFilename = CPLGetXMLValue(psIter, "SourceFilename", ""); if (osSourceFilename.empty()) { CPLError(CE_Failure, CPLE_AppDefined, "SpectralBand.SourceFilename missing"); goto error; } int bRelativeToVRT = atoi( CPLGetXMLValue(psIter, "SourceFilename.relativetoVRT", "0")); if (bRelativeToVRT) { const std::string osAbs = CPLProjectRelativeFilenameSafe( pszVRTPathIn, osSourceFilename.c_str()); m_oMapToRelativeFilenames[osAbs] = osSourceFilename; osSourceFilename = osAbs; } poDataset = oMapNamesToDataset[osSourceFilename]; if (poDataset == nullptr) { const CPLStringList aosOpenOptions( GDALDeserializeOpenOptionsFromXML(psIter)); poDataset = GDALDataset::Open( osSourceFilename, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr, aosOpenOptions.List()); if (poDataset == nullptr) { goto error; } oMapNamesToDataset[osSourceFilename] = poDataset; } else { poDataset->Reference(); } } m_apoDatasetsToReleaseRef.push_back( std::unique_ptr( poDataset)); poDataset = m_apoDatasetsToReleaseRef.back().get(); // Check that the spectral band has a georeferencing consistent // of the pan band. Allow an error of at most the size of one pixel // of the spectral band. const auto poSpectralSRS = poDataset->GetSpatialRef(); if (poPanSRS) { if (poSpectralSRS) { if (!poPanSRS->IsSame(poSpectralSRS)) { CPLError(CE_Warning, CPLE_AppDefined, "Pan dataset and %s do not seem to " "have same projection. Results might " "be incorrect", osSourceFilename.c_str()); } } else { CPLError(CE_Warning, CPLE_AppDefined, "Pan dataset has a projection, whereas %s " "not. Results might be incorrect", osSourceFilename.c_str()); } } else if (poSpectralSRS) { CPLError(CE_Warning, CPLE_AppDefined, "Pan dataset has no projection, whereas %s has " "one. Results might be incorrect", osSourceFilename.c_str()); } GDALGeoTransform spectralGT; if (poDataset->GetGeoTransform(spectralGT) != CE_None) { CPLError(CE_Failure, CPLE_AppDefined, "Spectral band has no associated geotransform"); goto error; } if (spectralGT[5] * panGT[5] < 0) { CPLError(CE_Failure, CPLE_NotSupported, "Spectral band vertical orientation is " "different from pan one"); goto error; } int bIsThisOneNonMatching = FALSE; double dfPixelSize = std::max(spectralGT[1], std::abs(spectralGT[5])); if (std::abs(panGT[0] - spectralGT[0]) > dfPixelSize || std::abs(panGT[3] - spectralGT[3]) > dfPixelSize) { bIsThisOneNonMatching = TRUE; if (m_eGTAdjustment == GTAdjust_None) { CPLError(CE_Warning, CPLE_AppDefined, "Georeferencing of top-left corner of pan " "dataset and %s do not match", osSourceFilename.c_str()); } } bFoundRotatingTerms = bFoundRotatingTerms || (spectralGT[2] != 0.0 || spectralGT[4] != 0.0); double dfLRSpectralX = spectralGT[0] + poDataset->GetRasterXSize() * spectralGT[1] + poDataset->GetRasterYSize() * spectralGT[2]; double dfLRSpectralY = spectralGT[3] + poDataset->GetRasterXSize() * spectralGT[4] + poDataset->GetRasterYSize() * spectralGT[5]; if (std::abs(dfLRPanX - dfLRSpectralX) > dfPixelSize || std::abs(dfLRPanY - dfLRSpectralY) > dfPixelSize) { bIsThisOneNonMatching = TRUE; if (m_eGTAdjustment == GTAdjust_None) { CPLError(CE_Warning, CPLE_AppDefined, "Georeferencing of bottom-right corner of " "pan dataset and %s do not match", osSourceFilename.c_str()); } } double dfThisMinY = dfLRSpectralY; double dfThisMaxY = spectralGT[3]; if (dfThisMinY > dfThisMaxY) std::swap(dfThisMinY, dfThisMaxY); if (bIsThisOneNonMatching && m_eGTAdjustment == GTAdjust_Union) { dfMinX = std::min(dfMinX, spectralGT[0]); dfMinY = std::min(dfMinY, dfThisMinY); dfMaxX = std::max(dfMaxX, dfLRSpectralX); dfMaxY = std::max(dfMaxY, dfThisMaxY); } else if (bIsThisOneNonMatching && m_eGTAdjustment == GTAdjust_Intersection) { dfMinX = std::max(dfMinX, spectralGT[0]); dfMinY = std::max(dfMinY, dfThisMinY); dfMaxX = std::min(dfMaxX, dfLRSpectralX); dfMaxY = std::min(dfMaxY, dfThisMaxY); } bFoundNonMatchingGT = bFoundNonMatchingGT || bIsThisOneNonMatching; iSpectralBand++; } if (nInputSpectralBandsIn && iSpectralBand != nInputSpectralBandsIn) { CPLError(CE_Failure, CPLE_AppDefined, "Less SpectralBand elements than in source array"); goto error; } /* -------------------------------------------------------------------- */ /* On-the-fly spatial extent adjustment if needed and asked. */ /* -------------------------------------------------------------------- */ if (bFoundNonMatchingGT && (m_eGTAdjustment == GTAdjust_Union || m_eGTAdjustment == GTAdjust_Intersection)) { if (bFoundRotatingTerms) { CPLError( CE_Failure, CPLE_NotSupported, "One of the panchromatic or spectral datasets has rotating " "terms in their geotransform matrix. Adjustment not possible"); goto error; } if (m_eGTAdjustment == GTAdjust_Intersection && (dfMinX >= dfMaxX || dfMinY >= dfMaxY)) { CPLError( CE_Failure, CPLE_NotSupported, "One of the panchromatic or spectral datasets has rotating " "terms in their geotransform matrix. Adjustment not possible"); goto error; } if (m_eGTAdjustment == GTAdjust_Union) CPLDebug("VRT", "Do union of bounding box of panchromatic and " "spectral datasets"); else CPLDebug("VRT", "Do intersection of bounding box of panchromatic " "and spectral datasets"); // If the pandataset needs adjustments, make sure the coordinates of the // union/intersection properly align with the grid of the pandataset // to avoid annoying sub-pixel shifts on the panchro band. double dfPixelSize = std::max(panGT[1], std::abs(panGT[5])); if (std::abs(panGT[0] - dfMinX) > dfPixelSize || std::abs(dfPanMaxY - dfMaxY) > dfPixelSize || std::abs(dfLRPanX - dfMaxX) > dfPixelSize || std::abs(dfPanMinY - dfMinY) > dfPixelSize) { dfMinX = panGT[0] + std::floor((dfMinX - panGT[0]) / panGT[1] + 0.5) * panGT[1]; dfMaxX = dfLRPanX + std::floor((dfMaxX - dfLRPanX) / panGT[1] + 0.5) * panGT[1]; if (panGT[5] > 0) { dfMinY = dfPanMinY + std::floor((dfMinY - dfPanMinY) / panGT[5] + 0.5) * panGT[5]; dfMaxY = dfPanMinY + std::floor((dfMaxY - dfPanMinY) / panGT[5] + 0.5) * panGT[5]; } else { dfMinY = dfPanMaxY + std::floor((dfMinY - dfPanMaxY) / panGT[5] + 0.5) * panGT[5]; dfMaxY = dfPanMaxY + std::floor((dfMaxY - dfPanMaxY) / panGT[5] + 0.5) * panGT[5]; } } std::map::iterator oIter = oMapNamesToDataset.begin(); for (; oIter != oMapNamesToDataset.end(); ++oIter) { GDALDataset *poSrcDS = oIter->second; GDALGeoTransform gt; if (poSrcDS->GetGeoTransform(gt) != CE_None) continue; // Check if this dataset needs adjustments dfPixelSize = std::max(gt[1], std::abs(gt[5])); dfPixelSize = std::max(panGT[1], dfPixelSize); dfPixelSize = std::max(std::abs(panGT[5]), dfPixelSize); double dfThisMinX = gt[0]; double dfThisMaxX = gt[0] + poSrcDS->GetRasterXSize() * gt[1]; double dfThisMaxY = gt[3]; double dfThisMinY = gt[3] + poSrcDS->GetRasterYSize() * gt[5]; if (dfThisMinY > dfThisMaxY) std::swap(dfThisMinY, dfThisMaxY); if (std::abs(dfThisMinX - dfMinX) <= dfPixelSize && std::abs(dfThisMaxY - dfMaxY) <= dfPixelSize && std::abs(dfThisMaxX - dfMaxX) <= dfPixelSize && std::abs(dfThisMinY - dfMinY) <= dfPixelSize) { continue; } GDALGeoTransform adjustedGT; adjustedGT[0] = dfMinX; adjustedGT[1] = gt[1]; adjustedGT[2] = 0; adjustedGT[3] = gt[5] > 0 ? dfMinY : dfMaxY; adjustedGT[4] = 0; adjustedGT[5] = gt[5]; const double dfAdjustedRasterXSize = 0.5 + (dfMaxX - dfMinX) / adjustedGT[1]; const double dfAdjustedRasterYSize = 0.5 + (dfMaxY - dfMinY) / std::abs(adjustedGT[5]); if (!(dfAdjustedRasterXSize >= 1 && dfAdjustedRasterXSize < INT_MAX && dfAdjustedRasterYSize >= 1 && dfAdjustedRasterYSize < INT_MAX)) { CPLError(CE_Failure, CPLE_AppDefined, "Datasets are too disjoint"); goto error; } const int nAdjustRasterXSize = static_cast(dfAdjustedRasterXSize); const int nAdjustRasterYSize = static_cast(dfAdjustedRasterYSize); VRTDataset *poVDS = new VRTDataset(nAdjustRasterXSize, nAdjustRasterYSize); poVDS->SetWritable(FALSE); poVDS->SetMetadataItem("UNDERLYING_DATASET", poSrcDS->GetDescription()); poVDS->SetGeoTransform(adjustedGT); poVDS->SetProjection(poPanDataset->GetProjectionRef()); for (int i = 0; i < poSrcDS->GetRasterCount(); i++) { GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(i + 1); poVDS->AddBand(poSrcBand->GetRasterDataType(), nullptr); VRTSourcedRasterBand *poVRTBand = static_cast( poVDS->GetRasterBand(i + 1)); const char *pszNBITS = poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE"); if (pszNBITS) poVRTBand->SetMetadataItem("NBITS", pszNBITS, "IMAGE_STRUCTURE"); VRTSimpleSource *poSimpleSource = new VRTSimpleSource(); poVRTBand->ConfigureSource( poSimpleSource, poSrcBand, FALSE, static_cast( std::floor((dfMinX - gt[0]) / gt[1] + 0.001)), gt[5] < 0 ? static_cast(std::floor( (dfMaxY - dfThisMaxY) / gt[5] + 0.001)) : static_cast(std::floor( (dfMinY - dfThisMinY) / gt[5] + 0.001)), static_cast(0.5 + (dfMaxX - dfMinX) / gt[1]), static_cast(0.5 + (dfMaxY - dfMinY) / std::abs(gt[5])), 0, 0, nAdjustRasterXSize, nAdjustRasterYSize); poVRTBand->AddSource(poSimpleSource); } oIter->second = poVDS; if (poSrcDS == poPanDataset) { panGT = adjustedGT; poPanDataset = poVDS; poPanBand = poPanDataset->GetRasterBand(nPanBand); nPanXSize = poPanDataset->GetRasterXSize(); nPanYSize = poPanDataset->GetRasterYSize(); } m_apoDatasetsToReleaseRef.push_back( std::unique_ptr( poVDS)); } } if (nRasterXSize == 0 && nRasterYSize == 0) { nRasterXSize = nPanXSize; nRasterYSize = nPanYSize; } else if (nRasterXSize != nPanXSize || nRasterYSize != nPanYSize) { CPLError(CE_Failure, CPLE_AppDefined, "Inconsistent declared VRT dimensions with panchro dataset"); goto error; } /* -------------------------------------------------------------------- */ /* Initialize all the general VRT stuff. This will even */ /* create the VRTPansharpenedRasterBands and initialize them. */ /* -------------------------------------------------------------------- */ eErr = VRTDataset::XMLInit(psTree, pszVRTPathIn); if (eErr != CE_None) { goto error; } /* -------------------------------------------------------------------- */ /* Inherit georeferencing info from panchro band if not defined */ /* in VRT. */ /* -------------------------------------------------------------------- */ { GDALGeoTransform outGT; if (GetGeoTransform(outGT) != CE_None && GetGCPCount() == 0 && GetSpatialRef() == nullptr) { if (bPanGeoTransformValid) { SetGeoTransform(panGT); } if (poPanSRS) { SetSpatialRef(poPanSRS); } } } /* -------------------------------------------------------------------- */ /* Parse rest of PansharpeningOptions */ /* -------------------------------------------------------------------- */ iSpectralBand = 0; for (const CPLXMLNode *psIter = psOptions->psChild; psIter; psIter = psIter->psNext) { if (psIter->eType != CXT_Element || !EQUAL(psIter->pszValue, "SpectralBand")) continue; GDALDataset *poDataset; GDALRasterBand *poBand; const char *pszDstBand = CPLGetXMLValue(psIter, "dstBand", nullptr); int nDstBand = -1; if (pszDstBand != nullptr) { nDstBand = atoi(pszDstBand); if (nDstBand <= 0) { CPLError(CE_Failure, CPLE_AppDefined, "SpectralBand.dstBand = '%s' invalid", pszDstBand); goto error; } } if (nInputSpectralBandsIn && pahInputSpectralBandsIn != nullptr) { poBand = GDALRasterBand::FromHandle( pahInputSpectralBandsIn[iSpectralBand]); poDataset = poBand->GetDataset(); if (poDataset) { poDataset = oMapNamesToDataset[CPLSPrintf("%p", poDataset)]; CPLAssert(poDataset); if (poDataset) poBand = poDataset->GetRasterBand(poBand->GetBand()); } } else { osSourceFilename = CPLGetXMLValue(psIter, "SourceFilename", ""); const bool bRelativeToVRT = CPL_TO_BOOL(atoi( CPLGetXMLValue(psIter, "SourceFilename.relativetoVRT", "0"))); if (bRelativeToVRT) osSourceFilename = CPLProjectRelativeFilenameSafe( pszVRTPathIn, osSourceFilename.c_str()); poDataset = oMapNamesToDataset[osSourceFilename]; CPLAssert(poDataset); const char *pszSourceBand = CPLGetXMLValue(psIter, "SourceBand", "1"); const int nBand = atoi(pszSourceBand); poBand = poDataset->GetRasterBand(nBand); if (poBand == nullptr) { CPLError(CE_Failure, CPLE_AppDefined, "%s invalid band of %s", pszSourceBand, osSourceFilename.c_str()); goto error; } } if (bHasNoData) { double dfSpectralNoData = poPanBand->GetNoDataValue(&bHasNoData); if (bHasNoData && dfSpectralNoData != dfNoData) bHasNoData = FALSE; } ahSpectralBands.push_back(poBand); if (nDstBand >= 1) { if (aMapDstBandToSpectralBand.find(nDstBand - 1) != aMapDstBandToSpectralBand.end()) { CPLError( CE_Failure, CPLE_AppDefined, "Another spectral band is already mapped to output band %d", nDstBand); goto error; } aMapDstBandToSpectralBand[nDstBand - 1] = static_cast(ahSpectralBands.size() - 1); } iSpectralBand++; } if (ahSpectralBands.empty()) { CPLError(CE_Failure, CPLE_AppDefined, "No spectral band defined"); goto error; } { const char *pszNoData = CPLGetXMLValue(psOptions, "NoData", nullptr); if (pszNoData) { if (EQUAL(pszNoData, "NONE")) { m_bNoDataDisabled = TRUE; bHasNoData = FALSE; } else { bHasNoData = TRUE; dfNoData = CPLAtof(pszNoData); } } } if (GetRasterCount() == 0) { for (int i = 0; i < static_cast(aMapDstBandToSpectralBand.size()); i++) { if (aMapDstBandToSpectralBand.find(i) == aMapDstBandToSpectralBand.end()) { CPLError(CE_Failure, CPLE_AppDefined, "Hole in SpectralBand.dstBand numbering"); goto error; } GDALRasterBand *poInputBand = GDALRasterBand::FromHandle( ahSpectralBands[aMapDstBandToSpectralBand[i]]); GDALRasterBand *poBand = new VRTPansharpenedRasterBand( this, i + 1, poInputBand->GetRasterDataType()); poBand->SetColorInterpretation( poInputBand->GetColorInterpretation()); if (bHasNoData) poBand->SetNoDataValue(dfNoData); SetBand(i + 1, poBand); } } else { int nIdxAsPansharpenedBand = 0; for (int i = 0; i < nBands; i++) { if (static_cast(GetRasterBand(i + 1)) ->IsPansharpenRasterBand()) { if (aMapDstBandToSpectralBand.find(i) == aMapDstBandToSpectralBand.end()) { CPLError(CE_Failure, CPLE_AppDefined, "Band %d of type VRTPansharpenedRasterBand, but " "no corresponding SpectralBand", i + 1); goto error; } else { static_cast( GetRasterBand(i + 1)) ->SetIndexAsPansharpenedBand(nIdxAsPansharpenedBand); nIdxAsPansharpenedBand++; } } } } // Figure out bit depth { const char *pszBitDepth = CPLGetXMLValue(psOptions, "BitDepth", nullptr); if (pszBitDepth == nullptr) pszBitDepth = GDALRasterBand::FromHandle(ahSpectralBands[0]) ->GetMetadataItem("NBITS", "IMAGE_STRUCTURE"); if (pszBitDepth) nBitDepth = atoi(pszBitDepth); if (nBitDepth) { for (int i = 0; i < nBands; i++) { if (!static_cast(GetRasterBand(i + 1)) ->IsPansharpenRasterBand()) continue; if (GetRasterBand(i + 1)->GetMetadataItem( "NBITS", "IMAGE_STRUCTURE") == nullptr) { if (nBitDepth != 8 && nBitDepth != 16 && nBitDepth != 32) { GetRasterBand(i + 1)->SetMetadataItem( "NBITS", CPLSPrintf("%d", nBitDepth), "IMAGE_STRUCTURE"); } } else if (nBitDepth == 8 || nBitDepth == 16 || nBitDepth == 32) { GetRasterBand(i + 1)->SetMetadataItem("NBITS", nullptr, "IMAGE_STRUCTURE"); } } } } if (GDALGetRasterBandXSize(ahSpectralBands[0]) > GDALGetRasterBandXSize(poPanBand) || GDALGetRasterBandYSize(ahSpectralBands[0]) > GDALGetRasterBandYSize(poPanBand)) { CPLError(CE_Warning, CPLE_AppDefined, "Dimensions of spectral band larger than panchro band"); } aMapDstBandToSpectralBandIter = aMapDstBandToSpectralBand.begin(); for (; aMapDstBandToSpectralBandIter != aMapDstBandToSpectralBand.end(); ++aMapDstBandToSpectralBandIter) { const int nDstBand = 1 + aMapDstBandToSpectralBandIter->first; if (nDstBand > nBands || !static_cast(GetRasterBand(nDstBand)) ->IsPansharpenRasterBand()) { CPLError(CE_Failure, CPLE_AppDefined, "SpectralBand.dstBand = '%d' invalid", nDstBand); goto error; } } if (adfWeights.empty()) { for (int i = 0; i < static_cast(ahSpectralBands.size()); i++) { adfWeights.push_back(1.0 / ahSpectralBands.size()); } } else if (adfWeights.size() != ahSpectralBands.size()) { CPLError(CE_Failure, CPLE_AppDefined, "%d weights defined, but %d input spectral bands", static_cast(adfWeights.size()), static_cast(ahSpectralBands.size())); goto error; } if (aMapDstBandToSpectralBand.empty()) { CPLError(CE_Warning, CPLE_AppDefined, "No spectral band is mapped to an output band"); } /* -------------------------------------------------------------------- */ /* Instantiate poPansharpener */ /* -------------------------------------------------------------------- */ psPanOptions.reset(GDALCreatePansharpenOptions()); psPanOptions->ePansharpenAlg = GDAL_PSH_WEIGHTED_BROVEY; psPanOptions->eResampleAlg = eResampleAlg; psPanOptions->nBitDepth = nBitDepth; psPanOptions->nWeightCount = static_cast(adfWeights.size()); psPanOptions->padfWeights = static_cast(CPLMalloc(sizeof(double) * adfWeights.size())); memcpy(psPanOptions->padfWeights, &adfWeights[0], sizeof(double) * adfWeights.size()); psPanOptions->hPanchroBand = poPanBand; psPanOptions->nInputSpectralBands = static_cast(ahSpectralBands.size()); psPanOptions->pahInputSpectralBands = static_cast( CPLMalloc(sizeof(GDALRasterBandH) * ahSpectralBands.size())); memcpy(psPanOptions->pahInputSpectralBands, &ahSpectralBands[0], sizeof(GDALRasterBandH) * ahSpectralBands.size()); psPanOptions->nOutPansharpenedBands = static_cast(aMapDstBandToSpectralBand.size()); psPanOptions->panOutPansharpenedBands = static_cast( CPLMalloc(sizeof(int) * aMapDstBandToSpectralBand.size())); aMapDstBandToSpectralBandIter = aMapDstBandToSpectralBand.begin(); for (int i = 0; aMapDstBandToSpectralBandIter != aMapDstBandToSpectralBand.end(); ++aMapDstBandToSpectralBandIter, ++i) { psPanOptions->panOutPansharpenedBands[i] = aMapDstBandToSpectralBandIter->second; } psPanOptions->bHasNoData = bHasNoData; psPanOptions->dfNoData = dfNoData; psPanOptions->nThreads = nThreads; if (nBands == psPanOptions->nOutPansharpenedBands) SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE"); m_poPansharpener = std::make_unique(); eErr = m_poPansharpener->Initialize(psPanOptions.get()); if (eErr != CE_None) { // Delete the pansharper object before closing the dataset // because it may have warped the bands into an intermediate VRT m_poPansharpener.reset(); // Close in reverse order (VRT firsts and real datasets after) for (int i = static_cast(m_apoDatasetsToReleaseRef.size()) - 1; i >= 0; i--) { m_apoDatasetsToReleaseRef[i].reset(); } m_apoDatasetsToReleaseRef.clear(); } return eErr; error: // Close in reverse order (VRT firsts and real datasets after) for (int i = static_cast(m_apoDatasetsToReleaseRef.size()) - 1; i >= 0; i--) { m_apoDatasetsToReleaseRef[i].reset(); } m_apoDatasetsToReleaseRef.clear(); return CE_Failure; } /************************************************************************/ /* SerializeToXML() */ /************************************************************************/ CPLXMLNode *VRTPansharpenedDataset::SerializeToXML(const char *pszVRTPathIn) { CPLXMLNode *psTree = VRTDataset::SerializeToXML(pszVRTPathIn); if (psTree == nullptr) return psTree; /* -------------------------------------------------------------------- */ /* Set subclass. */ /* -------------------------------------------------------------------- */ CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"), CXT_Text, "VRTPansharpenedDataset"); /* -------------------------------------------------------------------- */ /* Serialize the block size. */ /* -------------------------------------------------------------------- */ CPLCreateXMLElementAndValue(psTree, "BlockXSize", CPLSPrintf("%d", m_nBlockXSize)); CPLCreateXMLElementAndValue(psTree, "BlockYSize", CPLSPrintf("%d", m_nBlockYSize)); /* -------------------------------------------------------------------- */ /* Serialize the options. */ /* -------------------------------------------------------------------- */ if (m_poPansharpener == nullptr) return psTree; const GDALPansharpenOptions *psOptions = m_poPansharpener->GetOptions(); if (psOptions == nullptr) return psTree; CPLXMLNode *psOptionsNode = CPLCreateXMLNode(psTree, CXT_Element, "PansharpeningOptions"); if (psOptions->ePansharpenAlg == GDAL_PSH_WEIGHTED_BROVEY) { CPLCreateXMLElementAndValue(psOptionsNode, "Algorithm", "WeightedBrovey"); } else { CPLAssert(false); } if (psOptions->nWeightCount) { CPLString osWeights; for (int i = 0; i < psOptions->nWeightCount; i++) { if (i) osWeights += ","; osWeights += CPLSPrintf("%.16g", psOptions->padfWeights[i]); } CPLCreateXMLElementAndValue( CPLCreateXMLNode(psOptionsNode, CXT_Element, "AlgorithmOptions"), "Weights", osWeights.c_str()); } CPLCreateXMLElementAndValue( psOptionsNode, "Resampling", GDALRasterIOGetResampleAlg(psOptions->eResampleAlg)); if (psOptions->nThreads == -1) { CPLCreateXMLElementAndValue(psOptionsNode, "NumThreads", "ALL_CPUS"); } else if (psOptions->nThreads > 1) { CPLCreateXMLElementAndValue(psOptionsNode, "NumThreads", CPLSPrintf("%d", psOptions->nThreads)); } if (psOptions->nBitDepth) CPLCreateXMLElementAndValue(psOptionsNode, "BitDepth", CPLSPrintf("%d", psOptions->nBitDepth)); const char *pszAdjust = nullptr; switch (m_eGTAdjustment) { case GTAdjust_Union: pszAdjust = "Union"; break; case GTAdjust_Intersection: pszAdjust = "Intersection"; break; case GTAdjust_None: pszAdjust = "None"; break; case GTAdjust_NoneWithoutWarning: pszAdjust = "NoneWithoutWarning"; break; default: break; } if (psOptions->bHasNoData) { CPLCreateXMLElementAndValue(psOptionsNode, "NoData", CPLSPrintf("%.16g", psOptions->dfNoData)); } else if (m_bNoDataDisabled) { CPLCreateXMLElementAndValue(psOptionsNode, "NoData", "None"); } if (pszAdjust) CPLCreateXMLElementAndValue(psOptionsNode, "SpatialExtentAdjustment", pszAdjust); if (psOptions->hPanchroBand) { CPLXMLNode *psBand = CPLCreateXMLNode(psOptionsNode, CXT_Element, "PanchroBand"); GDALRasterBand *poBand = GDALRasterBand::FromHandle(psOptions->hPanchroBand); if (poBand->GetDataset()) { auto poPanchroDS = poBand->GetDataset(); std::map::iterator oIter = m_oMapToRelativeFilenames.find(poPanchroDS->GetDescription()); if (oIter == m_oMapToRelativeFilenames.end()) { const char *pszUnderlyingDS = poPanchroDS->GetMetadataItem("UNDERLYING_DATASET"); if (!pszUnderlyingDS) pszUnderlyingDS = poPanchroDS->GetDescription(); CPLCreateXMLElementAndValue(psBand, "SourceFilename", pszUnderlyingDS); } else { CPLXMLNode *psSourceFilename = CPLCreateXMLElementAndValue( psBand, "SourceFilename", oIter->second); CPLCreateXMLNode(CPLCreateXMLNode(psSourceFilename, CXT_Attribute, "relativeToVRT"), CXT_Text, "1"); } GDALSerializeOpenOptionsToXML(psBand, poPanchroDS->GetOpenOptions()); CPLCreateXMLElementAndValue(psBand, "SourceBand", CPLSPrintf("%d", poBand->GetBand())); } } for (int i = 0; i < psOptions->nInputSpectralBands; i++) { CPLXMLNode *psBand = CPLCreateXMLNode(psOptionsNode, CXT_Element, "SpectralBand"); for (int j = 0; j < psOptions->nOutPansharpenedBands; j++) { if (psOptions->panOutPansharpenedBands[j] == i) { for (int k = 0; k < nBands; k++) { if (static_cast(GetRasterBand(k + 1)) ->IsPansharpenRasterBand()) { if (static_cast( GetRasterBand(k + 1)) ->GetIndexAsPansharpenedBand() == j) { CPLCreateXMLNode(CPLCreateXMLNode(psBand, CXT_Attribute, "dstBand"), CXT_Text, CPLSPrintf("%d", k + 1)); break; } } } break; } } GDALRasterBand *poBand = GDALRasterBand::FromHandle(psOptions->pahInputSpectralBands[i]); if (poBand->GetDataset()) { auto poSpectralDS = poBand->GetDataset(); std::map::iterator oIter = m_oMapToRelativeFilenames.find(poSpectralDS->GetDescription()); if (oIter == m_oMapToRelativeFilenames.end()) { const char *pszUnderlyingDS = poSpectralDS->GetMetadataItem("UNDERLYING_DATASET"); if (!pszUnderlyingDS) pszUnderlyingDS = poSpectralDS->GetDescription(); CPLCreateXMLElementAndValue(psBand, "SourceFilename", pszUnderlyingDS); } else { CPLXMLNode *psSourceFilename = CPLCreateXMLElementAndValue( psBand, "SourceFilename", oIter->second); CPLCreateXMLNode(CPLCreateXMLNode(psSourceFilename, CXT_Attribute, "relativeToVRT"), CXT_Text, "1"); } GDALSerializeOpenOptionsToXML(psBand, poSpectralDS->GetOpenOptions()); CPLCreateXMLElementAndValue(psBand, "SourceBand", CPLSPrintf("%d", poBand->GetBand())); } } return psTree; } /************************************************************************/ /* GetBlockSize() */ /************************************************************************/ void VRTPansharpenedDataset::GetBlockSize(int *pnBlockXSize, int *pnBlockYSize) const { assert(nullptr != pnBlockXSize); assert(nullptr != pnBlockYSize); *pnBlockXSize = m_nBlockXSize; *pnBlockYSize = m_nBlockYSize; } /************************************************************************/ /* AddBand() */ /************************************************************************/ CPLErr VRTPansharpenedDataset::AddBand(CPL_UNUSED GDALDataType eType, CPL_UNUSED char **papszOptions) { CPLError(CE_Failure, CPLE_NotSupported, "AddBand() not supported"); return CE_Failure; } /************************************************************************/ /* IRasterIO() */ /************************************************************************/ CPLErr VRTPansharpenedDataset::IRasterIO( GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize, void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType, int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace, GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg) { if (eRWFlag == GF_Write) return CE_Failure; /* Try to pass the request to the most appropriate overview dataset */ if (nBufXSize < nXSize && nBufYSize < nYSize) { int bTried; CPLErr eErr = TryOverviewRasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize, eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace, nBandSpace, psExtraArg, &bTried); if (bTried) return eErr; } const int nDataTypeSize = GDALGetDataTypeSizeBytes(eBufType); if (nXSize == nBufXSize && nYSize == nBufYSize && nDataTypeSize == nPixelSpace && nLineSpace == nPixelSpace * nBufXSize && nBandSpace == nLineSpace * nBufYSize && nBandCount == nBands) { for (int i = 0; i < nBands; i++) { if (panBandMap[i] != i + 1 || !static_cast(GetRasterBand(i + 1)) ->IsPansharpenRasterBand()) { goto default_path; } } //{static int bDone = 0; if (!bDone) printf("(2)\n"); bDone = 1; } return m_poPansharpener->ProcessRegion(nXOff, nYOff, nXSize, nYSize, pData, eBufType); } default_path: //{static int bDone = 0; if (!bDone) printf("(3)\n"); bDone = 1; } return VRTDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize, eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace, nBandSpace, psExtraArg); } /************************************************************************/ /* ==================================================================== */ /* VRTPansharpenedRasterBand */ /* ==================================================================== */ /************************************************************************/ /************************************************************************/ /* VRTPansharpenedRasterBand() */ /************************************************************************/ VRTPansharpenedRasterBand::VRTPansharpenedRasterBand(GDALDataset *poDSIn, int nBandIn, GDALDataType eDataTypeIn) : m_nIndexAsPansharpenedBand(nBandIn - 1) { Initialize(poDSIn->GetRasterXSize(), poDSIn->GetRasterYSize()); poDS = poDSIn; nBand = nBandIn; eAccess = GA_Update; eDataType = eDataTypeIn; static_cast(poDS)->GetBlockSize(&nBlockXSize, &nBlockYSize); } /************************************************************************/ /* ~VRTPansharpenedRasterBand() */ /************************************************************************/ VRTPansharpenedRasterBand::~VRTPansharpenedRasterBand() { FlushCache(true); } /************************************************************************/ /* IReadBlock() */ /************************************************************************/ CPLErr VRTPansharpenedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage) { const int nReqXOff = nBlockXOff * nBlockXSize; const int nReqYOff = nBlockYOff * nBlockYSize; const int nReqXSize = std::min(nBlockXSize, nRasterXSize - nReqXOff); const int nReqYSize = std::min(nBlockYSize, nRasterYSize - nReqYOff); //{static int bDone = 0; if (!bDone) printf("(4)\n"); bDone = 1; } const int nDataTypeSize = GDALGetDataTypeSizeBytes(eDataType); GDALRasterIOExtraArg sExtraArg; INIT_RASTERIO_EXTRA_ARG(sExtraArg); if (IRasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize, pImage, nReqXSize, nReqYSize, eDataType, nDataTypeSize, static_cast(nDataTypeSize) * nReqXSize, &sExtraArg) != CE_None) { return CE_Failure; } if (nReqXSize < nBlockXSize) { for (int j = nReqYSize - 1; j >= 0; j--) { memmove(static_cast(pImage) + static_cast(j) * nDataTypeSize * nBlockXSize, static_cast(pImage) + static_cast(j) * nDataTypeSize * nReqXSize, static_cast(nReqXSize) * nDataTypeSize); memset(static_cast(pImage) + (static_cast(j) * nBlockXSize + nReqXSize) * nDataTypeSize, 0, static_cast(nBlockXSize - nReqXSize) * nDataTypeSize); } } if (nReqYSize < nBlockYSize) { memset(static_cast(pImage) + static_cast(nReqYSize) * nBlockXSize * nDataTypeSize, 0, static_cast(nBlockYSize - nReqYSize) * nBlockXSize * nDataTypeSize); } // Cache other bands CPLErr eErr = CE_None; VRTPansharpenedDataset *poGDS = static_cast(poDS); if (poGDS->nBands != 1 && !poGDS->m_bLoadingOtherBands) { poGDS->m_bLoadingOtherBands = TRUE; for (int iOtherBand = 1; iOtherBand <= poGDS->nBands; iOtherBand++) { if (iOtherBand == nBand) continue; GDALRasterBlock *poBlock = poGDS->GetRasterBand(iOtherBand) ->GetLockedBlockRef(nBlockXOff, nBlockYOff); if (poBlock == nullptr) { eErr = CE_Failure; break; } poBlock->DropLock(); } poGDS->m_bLoadingOtherBands = FALSE; } return eErr; } /************************************************************************/ /* IRasterIO() */ /************************************************************************/ CPLErr VRTPansharpenedRasterBand::IRasterIO( GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize, void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType, GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg) { if (eRWFlag == GF_Write) return CE_Failure; VRTPansharpenedDataset *poGDS = static_cast(poDS); /* Try to pass the request to the most appropriate overview dataset */ if (nBufXSize < nXSize && nBufYSize < nYSize) { int bTried; CPLErr eErr = TryOverviewRasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg, &bTried); if (bTried) return eErr; } const int nDataTypeSize = GDALGetDataTypeSizeBytes(eBufType); if (nDataTypeSize > 0 && nXSize == nBufXSize && nYSize == nBufYSize && nDataTypeSize == nPixelSpace && nLineSpace == nPixelSpace * nBufXSize) { const GDALPansharpenOptions *psOptions = poGDS->m_poPansharpener->GetOptions(); // Have we already done this request for another band ? // If so use the cached result const size_t nBufferSizePerBand = static_cast(nXSize) * nYSize * nDataTypeSize; if (nXOff == poGDS->m_nLastBandRasterIOXOff && nYOff >= poGDS->m_nLastBandRasterIOYOff && nXSize == poGDS->m_nLastBandRasterIOXSize && nYOff + nYSize <= poGDS->m_nLastBandRasterIOYOff + poGDS->m_nLastBandRasterIOYSize && eBufType == poGDS->m_eLastBandRasterIODataType) { //{static int bDone = 0; if (!bDone) printf("(6)\n"); bDone = 1; } if (poGDS->m_pabyLastBufferBandRasterIO == nullptr) return CE_Failure; const size_t nBufferSizePerBandCached = static_cast(nXSize) * poGDS->m_nLastBandRasterIOYSize * nDataTypeSize; memcpy(pData, poGDS->m_pabyLastBufferBandRasterIO + nBufferSizePerBandCached * m_nIndexAsPansharpenedBand + static_cast(nYOff - poGDS->m_nLastBandRasterIOYOff) * nXSize * nDataTypeSize, nBufferSizePerBand); return CE_None; } int nYSizeToCache = nYSize; if (nYSize == 1 && nXSize == nRasterXSize) { //{static int bDone = 0; if (!bDone) printf("(7)\n"); bDone = 1; } // For efficiency, try to cache at leak 256 K nYSizeToCache = (256 * 1024) / nXSize / nDataTypeSize; if (nYSizeToCache == 0) nYSizeToCache = 1; else if (nYOff + nYSizeToCache > nRasterYSize) nYSizeToCache = nRasterYSize - nYOff; } const GUIntBig nBufferSize = static_cast(nXSize) * nYSizeToCache * nDataTypeSize * psOptions->nOutPansharpenedBands; // Check the we don't overflow (for 32 bit platforms) if (nBufferSize > std::numeric_limits::max() / 2) { CPLError(CE_Failure, CPLE_OutOfMemory, "Out of memory error while allocating working buffers"); return CE_Failure; } GByte *pabyTemp = static_cast( VSI_REALLOC_VERBOSE(poGDS->m_pabyLastBufferBandRasterIO, static_cast(nBufferSize))); if (pabyTemp == nullptr) { return CE_Failure; } poGDS->m_nLastBandRasterIOXOff = nXOff; poGDS->m_nLastBandRasterIOYOff = nYOff; poGDS->m_nLastBandRasterIOXSize = nXSize; poGDS->m_nLastBandRasterIOYSize = nYSizeToCache; poGDS->m_eLastBandRasterIODataType = eBufType; poGDS->m_pabyLastBufferBandRasterIO = pabyTemp; CPLErr eErr = poGDS->m_poPansharpener->ProcessRegion( nXOff, nYOff, nXSize, nYSizeToCache, poGDS->m_pabyLastBufferBandRasterIO, eBufType); if (eErr == CE_None) { //{static int bDone = 0; if (!bDone) printf("(8)\n"); bDone = 1; } size_t nBufferSizePerBandCached = static_cast(nXSize) * poGDS->m_nLastBandRasterIOYSize * nDataTypeSize; memcpy(pData, poGDS->m_pabyLastBufferBandRasterIO + nBufferSizePerBandCached * m_nIndexAsPansharpenedBand, nBufferSizePerBand); } else { VSIFree(poGDS->m_pabyLastBufferBandRasterIO); poGDS->m_pabyLastBufferBandRasterIO = nullptr; } return eErr; } //{static int bDone = 0; if (!bDone) printf("(9)\n"); bDone = 1; } CPLErr eErr = VRTRasterBand::IRasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg); return eErr; } /************************************************************************/ /* SerializeToXML() */ /************************************************************************/ CPLXMLNode * VRTPansharpenedRasterBand::SerializeToXML(const char *pszVRTPathIn, bool &bHasWarnedAboutRAMUsage, size_t &nAccRAMUsage) { CPLXMLNode *psTree = VRTRasterBand::SerializeToXML( pszVRTPathIn, bHasWarnedAboutRAMUsage, nAccRAMUsage); /* -------------------------------------------------------------------- */ /* Set subclass. */ /* -------------------------------------------------------------------- */ CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"), CXT_Text, "VRTPansharpenedRasterBand"); return psTree; } /************************************************************************/ /* GetOverviewCount() */ /************************************************************************/ int VRTPansharpenedRasterBand::GetOverviewCount() { VRTPansharpenedDataset *poGDS = static_cast(poDS); // Build on-the-fly overviews from overviews of pan and spectral bands if (poGDS->m_poPansharpener != nullptr && poGDS->m_apoOverviewDatasets.empty() && poGDS->m_poMainDataset == poGDS) { const GDALPansharpenOptions *psOptions = poGDS->m_poPansharpener->GetOptions(); GDALRasterBand *poPanBand = GDALRasterBand::FromHandle(psOptions->hPanchroBand); int nOvrCount = poPanBand->GetOverviewCount(); if (nOvrCount > 0) { for (int i = 0; i < poGDS->GetRasterCount(); i++) { if (!static_cast(poGDS->GetRasterBand(i + 1)) ->IsPansharpenRasterBand()) { return 0; } } // Limit number of overviews of the VRTPansharpenedRasterBand to // the minimum number of overviews of the pan and spectral source // bands, and also make sure all spectral bands have overviews // of the same dimension for a given level. std::vector> sizeSpectralOverviews; for (int i = 0; i < psOptions->nInputSpectralBands; i++) { auto poSpectralBand = GDALRasterBand::FromHandle( psOptions->pahInputSpectralBands[i]); nOvrCount = std::min(nOvrCount, poSpectralBand->GetOverviewCount()); if (i == 0) { for (int iOvr = 0; iOvr < nOvrCount; ++iOvr) { auto poOvrBand = poSpectralBand->GetOverview(iOvr); sizeSpectralOverviews.emplace_back( poOvrBand->GetXSize(), poOvrBand->GetYSize()); } } else { for (int iOvr = 0; iOvr < nOvrCount; ++iOvr) { auto poOvrBand = poSpectralBand->GetOverview(iOvr); if (sizeSpectralOverviews[iOvr].first != poOvrBand->GetXSize() || sizeSpectralOverviews[iOvr].second != poOvrBand->GetYSize()) { nOvrCount = iOvr; break; } } } } auto poPanBandDS = poPanBand->GetDataset(); for (int j = 0; j < nOvrCount; j++) { std::unique_ptr poPanOvrDS(GDALCreateOverviewDataset(poPanBandDS, j, true)); if (!poPanOvrDS) { CPLError(CE_Warning, CPLE_AppDefined, "GDALCreateOverviewDataset(poPanBandDS, %d, true) " "failed", j); break; } GDALRasterBand *poPanOvrBand = poPanOvrDS->GetRasterBand(poPanBand->GetBand()); auto poOvrDS = std::make_unique( poPanOvrBand->GetXSize(), poPanOvrBand->GetYSize()); poOvrDS->m_apoDatasetsToReleaseRef.push_back( std::move(poPanOvrDS)); for (int i = 0; i < poGDS->GetRasterCount(); i++) { GDALRasterBand *poSrcBand = poGDS->GetRasterBand(i + 1); auto poBand = std::make_unique( poOvrDS.get(), i + 1, poSrcBand->GetRasterDataType()); const char *pszNBITS = poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE"); if (pszNBITS) poBand->SetMetadataItem("NBITS", pszNBITS, "IMAGE_STRUCTURE"); int bHasNoData = FALSE; const double dfNoData = poSrcBand->GetNoDataValue(&bHasNoData); if (bHasNoData) poBand->SetNoDataValue(dfNoData); poBand->SetColorInterpretation( poSrcBand->GetColorInterpretation()); poOvrDS->SetBand(i + 1, std::move(poBand)); } std::unique_ptr psPanOvrOptions(GDALClonePansharpenOptions(psOptions), GDALDestroyPansharpenOptions); psPanOvrOptions->hPanchroBand = poPanOvrBand; for (int i = 0; i < psOptions->nInputSpectralBands; i++) { auto poSpectralBand = GDALRasterBand::FromHandle( psOptions->pahInputSpectralBands[i]); std::unique_ptr poSpectralOvrDS(GDALCreateOverviewDataset( poSpectralBand->GetDataset(), j, true)); if (!poSpectralOvrDS) { CPLError(CE_Warning, CPLE_AppDefined, "GDALCreateOverviewDataset(poSpectralBand->" "GetDataset(), %d, true) failed", j); return 0; } psPanOvrOptions->pahInputSpectralBands[i] = poSpectralOvrDS->GetRasterBand( poSpectralBand->GetBand()); poOvrDS->m_apoDatasetsToReleaseRef.push_back( std::move(poSpectralOvrDS)); } poOvrDS->m_poPansharpener = std::make_unique(); if (poOvrDS->m_poPansharpener->Initialize( psPanOvrOptions.get()) != CE_None) { CPLError(CE_Warning, CPLE_AppDefined, "Unable to initialize pansharpener."); return 0; } poOvrDS->m_poMainDataset = poGDS; poOvrDS->SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE"); poGDS->m_apoOverviewDatasets.push_back(std::move(poOvrDS)); } } } return static_cast(poGDS->m_apoOverviewDatasets.size()); } /************************************************************************/ /* GetOverview() */ /************************************************************************/ GDALRasterBand *VRTPansharpenedRasterBand::GetOverview(int iOvr) { if (iOvr < 0 || iOvr >= GetOverviewCount()) return nullptr; VRTPansharpenedDataset *poGDS = static_cast(poDS); return poGDS->m_apoOverviewDatasets[iOvr]->GetRasterBand(nBand); } /*! @endcond */