/****************************************************************************** * * 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. * * Note that no reference will be taken on the passed bands. Consequently, * they or their dataset to which they belong to must 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_poPansharpener(nullptr), 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() { if (m_poMainDataset == nullptr) return FALSE; VRTPansharpenedDataset *poMainDatasetLocal = m_poMainDataset; m_poMainDataset = nullptr; 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_apoDatasetsToClose. for (size_t i = 0; i < m_apoOverviewDatasets.size(); i++) { bHasDroppedRef = TRUE; delete m_apoOverviewDatasets[i]; } m_apoOverviewDatasets.resize(0); if (m_poPansharpener != nullptr) { // Delete the pansharper object before closing the dataset // because it may have warped the bands into an intermediate VRT delete m_poPansharpener; m_poPansharpener = nullptr; // Close in reverse order (VRT firsts and real datasets after) for (int i = static_cast(m_apoDatasetsToClose.size()) - 1; i >= 0; i--) { bHasDroppedRef = TRUE; GDALClose(m_apoDatasetsToClose[i]); } m_apoDatasetsToClose.resize(0); } if (poMainDatasetLocal != this) { // To avoid killing us for (size_t i = 0; i < poMainDatasetLocal->m_apoOverviewDatasets.size(); i++) { if (poMainDatasetLocal->m_apoOverviewDatasets[i] == this) { poMainDatasetLocal->m_apoOverviewDatasets[i] = nullptr; break; } } bHasDroppedRef |= poMainDatasetLocal->CloseDependentDatasets(); } return bHasDroppedRef; } /************************************************************************/ /* GetFileList() */ /************************************************************************/ char **VRTPansharpenedDataset::GetFileList() { char **papszFileList = GDALDataset::GetFileList(); if (m_poPansharpener != nullptr) { GDALPansharpenOptions *psOptions = m_poPansharpener->GetOptions(); if (psOptions != nullptr) { std::set oSetNames; if (psOptions->hPanchroBand != nullptr) { GDALDatasetH hDS = GDALGetBandDataset(psOptions->hPanchroBand); if (hDS != nullptr) { papszFileList = CSLAddString(papszFileList, GDALGetDescription(hDS)); oSetNames.insert(GDALGetDescription(hDS)); } } for (int i = 0; i < psOptions->nInputSpectralBands; i++) { if (psOptions->pahInputSpectralBands[i] != nullptr) { GDALDatasetH hDS = GDALGetBandDataset(psOptions->pahInputSpectralBands[i]); if (hDS != nullptr && oSetNames.find(GDALGetDescription( hDS)) == oSetNames.end()) { papszFileList = CSLAddString(papszFileList, GDALGetDescription(hDS)); oSetNames.insert(GDALGetDescription(hDS)); } } } } } 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; GDALPansharpenOptions *psPanOptions; /* -------------------------------------------------------------------- */ /* 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; GDALDataset *poPanDatasetToClose = 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; } poPanDatasetToClose = poPanDataset; 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(poPanDatasetToClose); return CE_Failure; } oMapNamesToDataset[osSourceFilename] = poPanDataset; m_apoDatasetsToClose.push_back(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; } // 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); GDALClose(poPanDatasetToClose); m_apoDatasetsToClose.resize(0); 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; double adfPanGT[6] = {0, 0, 0, 0, 0, 0}; const bool bPanGeoTransformValid = (poPanDataset->GetGeoTransform(adfPanGT) == CE_None); if (!bPanGeoTransformValid) { CPLError(CE_Failure, CPLE_AppDefined, "Panchromatic band has no associated geotransform"); GDALClose(poPanDatasetToClose); m_apoDatasetsToClose.resize(0); return CE_Failure; } int nPanXSize = poPanBand->GetXSize(); int nPanYSize = poPanBand->GetYSize(); int bHasNoData = FALSE; double dfNoData = poPanBand->GetNoDataValue(&bHasNoData); double dfLRPanX = adfPanGT[0] + nPanXSize * adfPanGT[1] + nPanYSize * adfPanGT[2]; double dfLRPanY = adfPanGT[3] + nPanXSize * adfPanGT[4] + nPanYSize * adfPanGT[5]; bool bFoundRotatingTerms = (adfPanGT[2] != 0.0 || adfPanGT[4] != 0.0); double dfMinX = adfPanGT[0]; double dfMaxX = dfLRPanX; double dfMaxY = adfPanGT[3]; double dfMinY = dfLRPanY; if (dfMinY > dfMaxY) std::swap(dfMinY, dfMaxY); const double dfPanMinY = dfMinY; const double dfPanMaxY = dfMaxY; CPLString osPanProjection, osPanProjectionProj4; if (poPanDataset->GetProjectionRef()) { osPanProjection = poPanDataset->GetProjectionRef(); char *pszProj4 = nullptr; OGRSpatialReference oSRS(osPanProjection); if (oSRS.exportToProj4(&pszProj4) == OGRERR_NONE) osPanProjectionProj4 = pszProj4; CPLFree(pszProj4); } /* -------------------------------------------------------------------- */ /* 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; } 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; m_apoDatasetsToClose.push_back(poDataset); } } // 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. CPLString osProjection; if (poDataset->GetProjectionRef()) osProjection = poDataset->GetProjectionRef(); if (!osPanProjection.empty()) { if (!osProjection.empty()) { if (osPanProjection != osProjection) { CPLString osProjectionProj4; char *pszProj4 = nullptr; OGRSpatialReference oSRS(osProjection); if (oSRS.exportToProj4(&pszProj4) == OGRERR_NONE) osProjectionProj4 = pszProj4; CPLFree(pszProj4); if (osPanProjectionProj4 != osProjectionProj4) { 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 (!osProjection.empty()) { CPLError(CE_Warning, CPLE_AppDefined, "Pan dataset has no projection, whereas %s has " "one. Results might be incorrect", osSourceFilename.c_str()); } double adfSpectralGeoTransform[6]; if (poDataset->GetGeoTransform(adfSpectralGeoTransform) != CE_None) { CPLError(CE_Failure, CPLE_AppDefined, "Spectral band has no associated geotransform"); goto error; } if (adfSpectralGeoTransform[3] * adfPanGT[3] < 0) { CPLError(CE_Failure, CPLE_NotSupported, "Spectral band vertical orientation is " "different from pan one"); goto error; } int bIsThisOneNonMatching = FALSE; double dfPixelSize = std::max(adfSpectralGeoTransform[1], std::abs(adfSpectralGeoTransform[5])); if (std::abs(adfPanGT[0] - adfSpectralGeoTransform[0]) > dfPixelSize || std::abs(adfPanGT[3] - adfSpectralGeoTransform[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 || (adfSpectralGeoTransform[2] != 0.0 || adfSpectralGeoTransform[4] != 0.0); double dfLRSpectralX = adfSpectralGeoTransform[0] + poDataset->GetRasterXSize() * adfSpectralGeoTransform[1] + poDataset->GetRasterYSize() * adfSpectralGeoTransform[2]; double dfLRSpectralY = adfSpectralGeoTransform[3] + poDataset->GetRasterXSize() * adfSpectralGeoTransform[4] + poDataset->GetRasterYSize() * adfSpectralGeoTransform[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 = adfSpectralGeoTransform[3]; if (dfThisMinY > dfThisMaxY) std::swap(dfThisMinY, dfThisMaxY); if (bIsThisOneNonMatching && m_eGTAdjustment == GTAdjust_Union) { dfMinX = std::min(dfMinX, adfSpectralGeoTransform[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, adfSpectralGeoTransform[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(adfPanGT[1], std::abs(adfPanGT[5])); if (std::abs(adfPanGT[0] - dfMinX) > dfPixelSize || std::abs(dfPanMaxY - dfMaxY) > dfPixelSize || std::abs(dfLRPanX - dfMaxX) > dfPixelSize || std::abs(dfPanMinY - dfMinY) > dfPixelSize) { dfMinX = adfPanGT[0] + std::floor((dfMinX - adfPanGT[0]) / adfPanGT[1] + 0.5) * adfPanGT[1]; dfMaxX = dfLRPanX + std::floor((dfMaxX - dfLRPanX) / adfPanGT[1] + 0.5) * adfPanGT[1]; if (adfPanGT[5] > 0) { dfMinY = dfPanMinY + std::floor((dfMinY - dfPanMinY) / adfPanGT[5] + 0.5) * adfPanGT[5]; dfMaxY = dfPanMinY + std::floor((dfMaxY - dfPanMinY) / adfPanGT[5] + 0.5) * adfPanGT[5]; } else { dfMinY = dfPanMaxY + std::floor((dfMinY - dfPanMaxY) / adfPanGT[5] + 0.5) * adfPanGT[5]; dfMaxY = dfPanMaxY + std::floor((dfMaxY - dfPanMaxY) / adfPanGT[5] + 0.5) * adfPanGT[5]; } } std::map::iterator oIter = oMapNamesToDataset.begin(); for (; oIter != oMapNamesToDataset.end(); ++oIter) { GDALDataset *poSrcDS = oIter->second; double adfGT[6]; if (poSrcDS->GetGeoTransform(adfGT) != CE_None) continue; // Check if this dataset needs adjustments dfPixelSize = std::max(adfGT[1], std::abs(adfGT[5])); dfPixelSize = std::max(adfPanGT[1], dfPixelSize); dfPixelSize = std::max(std::abs(adfPanGT[5]), dfPixelSize); double dfThisMinX = adfGT[0]; double dfThisMaxX = adfGT[0] + poSrcDS->GetRasterXSize() * adfGT[1]; double dfThisMaxY = adfGT[3]; double dfThisMinY = adfGT[3] + poSrcDS->GetRasterYSize() * adfGT[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; } double adfAdjustedGT[6]; adfAdjustedGT[0] = dfMinX; adfAdjustedGT[1] = adfGT[1]; adfAdjustedGT[2] = 0; adfAdjustedGT[3] = adfGT[5] > 0 ? dfMinY : dfMaxY; adfAdjustedGT[4] = 0; adfAdjustedGT[5] = adfGT[5]; int nAdjustRasterXSize = static_cast(0.5 + (dfMaxX - dfMinX) / adfAdjustedGT[1]); int nAdjustRasterYSize = static_cast( 0.5 + (dfMaxY - dfMinY) / std::abs(adfAdjustedGT[5])); VRTDataset *poVDS = new VRTDataset(nAdjustRasterXSize, nAdjustRasterYSize); poVDS->SetWritable(FALSE); poVDS->SetDescription(std::string("Shifted ") .append(poSrcDS->GetDescription()) .c_str()); poVDS->SetGeoTransform(adfAdjustedGT); 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 - adfGT[0]) / adfGT[1] + 0.001)), adfGT[5] < 0 ? static_cast(std::floor( (dfMaxY - dfThisMaxY) / adfGT[5] + 0.001)) : static_cast(std::floor( (dfMinY - dfThisMinY) / adfGT[5] + 0.001)), static_cast(0.5 + (dfMaxX - dfMinX) / adfGT[1]), static_cast(0.5 + (dfMaxY - dfMinY) / std::abs(adfGT[5])), 0, 0, nAdjustRasterXSize, nAdjustRasterYSize); poVRTBand->AddSource(poSimpleSource); } oIter->second = poVDS; if (poSrcDS == poPanDataset) { memcpy(adfPanGT, adfAdjustedGT, 6 * sizeof(double)); poPanDataset = poVDS; poPanBand = poPanDataset->GetRasterBand(nPanBand); nPanXSize = poPanDataset->GetRasterXSize(); nPanYSize = poPanDataset->GetRasterYSize(); } m_apoDatasetsToClose.push_back(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. */ /* -------------------------------------------------------------------- */ { double adfOutGT[6]; if (GetGeoTransform(adfOutGT) != CE_None && GetGCPCount() == 0 && GetProjectionRef()[0] == '\0') { if (bPanGeoTransformValid) { SetGeoTransform(adfPanGT); } if (poPanDataset->GetProjectionRef() != nullptr && poPanDataset->GetProjectionRef()[0] != '\0') { SetProjection(poPanDataset->GetProjectionRef()); } } } /* -------------------------------------------------------------------- */ /* Parse rest of PansharpeningOptions */ /* -------------------------------------------------------------------- */ iSpectralBand = 0; for (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 = 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 = new GDALPansharpenOperation(); eErr = m_poPansharpener->Initialize(psPanOptions); if (eErr != CE_None) { // Delete the pansharper object before closing the dataset // because it may have warped the bands into an intermediate VRT delete m_poPansharpener; m_poPansharpener = nullptr; // Close in reverse order (VRT firsts and real datasets after) for (int i = static_cast(m_apoDatasetsToClose.size()) - 1; i >= 0; i--) { GDALClose(m_apoDatasetsToClose[i]); } m_apoDatasetsToClose.resize(0); } GDALDestroyPansharpenOptions(psPanOptions); return eErr; error: // Close in reverse order (VRT firsts and real datasets after) for (int i = static_cast(m_apoDatasetsToClose.size()) - 1; i >= 0; i--) { GDALClose(m_apoDatasetsToClose[i]); } m_apoDatasetsToClose.resize(0); 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; 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 poPanchoDS = poBand->GetDataset(); std::map::iterator oIter = m_oMapToRelativeFilenames.find(poPanchoDS->GetDescription()); if (oIter == m_oMapToRelativeFilenames.end()) { CPLCreateXMLElementAndValue(psBand, "SourceFilename", poPanchoDS->GetDescription()); } else { CPLXMLNode *psSourceFilename = CPLCreateXMLElementAndValue( psBand, "SourceFilename", oIter->second); CPLCreateXMLNode(CPLCreateXMLNode(psSourceFilename, CXT_Attribute, "relativeToVRT"), CXT_Text, "1"); } GDALSerializeOpenOptionsToXML(psBand, poPanchoDS->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()) { CPLCreateXMLElementAndValue(psBand, "SourceFilename", poSpectralDS->GetDescription()); } 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; int nReqXSize = nBlockXSize; int nReqYSize = nBlockYSize; if (nReqXOff + nReqXSize > nRasterXSize) nReqXSize = nRasterXSize - nReqXOff; if (nReqYOff + nReqYSize > nRasterYSize) nReqYSize = nRasterYSize - nReqYOff; //{static int bDone = 0; if (!bDone) printf("(4)\n"); bDone = 1; } const int nDataTypeSize = GDALGetDataTypeSize(eDataType) / 8; 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) { 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) { 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++) { auto 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()); VRTPansharpenedDataset *poOvrDS = new VRTPansharpenedDataset( poPanOvrBand->GetXSize(), poPanOvrBand->GetYSize()); poOvrDS->m_apoDatasetsToClose.push_back(poPanOvrDS); for (int i = 0; i < poGDS->GetRasterCount(); i++) { GDALRasterBand *poSrcBand = poGDS->GetRasterBand(i + 1); GDALRasterBand *poBand = new VRTPansharpenedRasterBand( poOvrDS, i + 1, poSrcBand->GetRasterDataType()); const char *pszNBITS = poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE"); if (pszNBITS) poBand->SetMetadataItem("NBITS", pszNBITS, "IMAGE_STRUCTURE"); poOvrDS->SetBand(i + 1, poBand); } GDALPansharpenOptions *psPanOvrOptions = GDALClonePansharpenOptions(psOptions); psPanOvrOptions->hPanchroBand = poPanOvrBand; for (int i = 0; i < psOptions->nInputSpectralBands; i++) { auto poSpectralBand = GDALRasterBand::FromHandle( psOptions->pahInputSpectralBands[i]); auto poSpectralOvrDS = GDALCreateOverviewDataset( poSpectralBand->GetDataset(), j, true); if (!poSpectralOvrDS) { CPLError(CE_Warning, CPLE_AppDefined, "GDALCreateOverviewDataset(poSpectralBand->" "GetDataset(), %d, true) failed", j); delete poOvrDS; GDALDestroyPansharpenOptions(psPanOvrOptions); return 0; } psPanOvrOptions->pahInputSpectralBands[i] = poSpectralOvrDS->GetRasterBand( poSpectralBand->GetBand()); poOvrDS->m_apoDatasetsToClose.push_back(poSpectralOvrDS); } poOvrDS->m_poPansharpener = new GDALPansharpenOperation(); if (poOvrDS->m_poPansharpener->Initialize(psPanOvrOptions) != CE_None) { CPLError(CE_Warning, CPLE_AppDefined, "Unable to initialize pansharpener."); delete poOvrDS; GDALDestroyPansharpenOptions(psPanOvrOptions); return 0; } GDALDestroyPansharpenOptions(psPanOvrOptions); poOvrDS->m_poMainDataset = poGDS; poOvrDS->SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE"); poGDS->m_apoOverviewDatasets.push_back(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 */