diff --git a/autotest/gcore/vrt_read.py b/autotest/gcore/vrt_read.py index 4f24787e176e..5fb40487e4e3 100755 --- a/autotest/gcore/vrt_read.py +++ b/autotest/gcore/vrt_read.py @@ -2862,3 +2862,113 @@ def test_vrt_read_float32_complex_source_from_cfloat32(): ds = gdal.Open("data/vrt/complex_non_zero_real_zero_imag_as_float32.vrt") assert struct.unpack("f" * 4, ds.ReadRaster()) == (1, 1, 1, 1) + + +############################################################################### + + +def test_vrt_resampling_with_mask_and_overviews(tmp_vsimem): + + filename1 = str(tmp_vsimem / "in1.tif") + ds = gdal.GetDriverByName("GTiff").Create(filename1, 100, 10) + ds.SetGeoTransform([0, 1, 0, 0, 0, -1]) + ds.GetRasterBand(1).WriteRaster(0, 0, 48, 10, b"\xFF" * (48 * 10)) + ds.GetRasterBand(1).WriteRaster(48, 0, 2, 10, b"\xF0" * (2 * 10)) + ds.CreateMaskBand(gdal.GMF_PER_DATASET) + ds.GetRasterBand(1).GetMaskBand().WriteRaster(0, 0, 52, 10, b"\xFF" * (52 * 10)) + ds.BuildOverviews("NEAR", [2]) + ds.GetRasterBand(1).GetOverview(0).Fill( + 127 + ) # to demonstrate we ignore overviews unfortunately + ds = None + + filename2 = str(tmp_vsimem / "in2.tif") + ds = gdal.GetDriverByName("GTiff").Create(filename2, 100, 10) + ds.SetGeoTransform([0, 1, 0, 0, 0, -1]) + ds.GetRasterBand(1).WriteRaster(48, 0, 52, 10, b"\xF0" * (52 * 10)) + ds.CreateMaskBand(gdal.GMF_PER_DATASET) + ds.GetRasterBand(1).GetMaskBand().WriteRaster(48, 0, 52, 10, b"\xFF" * (52 * 10)) + ds.BuildOverviews("NEAR", [2]) + ds.GetRasterBand(1).GetOverview(0).Fill( + 127 + ) # to demonstrate we ignore overviews unfortunately + ds = None + + vrt_filename = str(tmp_vsimem / "test.vrt") + gdal.BuildVRT(vrt_filename, [filename1, filename2], resampleAlg=gdal.GRIORA_Cubic) + + ds = gdal.Open(vrt_filename) + assert ( + ds.ReadRaster(buf_xsize=10, buf_ysize=1) + == b"\xFF\xFF\xFF\xFF\xFC\xF0\xF0\xF0\xF0\xF0" + ) + assert ( + ds.GetRasterBand(1).GetMaskBand().ReadRaster(buf_xsize=10, buf_ysize=1) + == b"\xFF" * 10 + ) + + vrt_filename = str(tmp_vsimem / "test.vrt") + gdal.BuildVRT( + vrt_filename, [filename1, filename2], resampleAlg=gdal.GRIORA_Bilinear + ) + + ds = gdal.Open(vrt_filename) + assert ( + ds.ReadRaster(buf_xsize=10, buf_ysize=1) + == b"\xFF\xFF\xFF\xFF\xFB\xF1\xF0\xF0\xF0\xF0" + ) + assert ( + ds.GetRasterBand(1).GetMaskBand().ReadRaster(buf_xsize=10, buf_ysize=1) + == b"\xFF" * 10 + ) + + +############################################################################### + + +def test_vrt_resampling_with_alpha_and_overviews(tmp_vsimem): + + filename1 = str(tmp_vsimem / "in1.tif") + ds = gdal.GetDriverByName("GTiff").Create( + filename1, 100, 10, 2, options=["ALPHA=YES"] + ) + ds.SetGeoTransform([0, 1, 0, 0, 0, -1]) + ds.GetRasterBand(1).WriteRaster(0, 0, 48, 10, b"\xFF" * (48 * 10)) + ds.GetRasterBand(1).WriteRaster(48, 0, 2, 10, b"\xF0" * (2 * 10)) + ds.GetRasterBand(2).WriteRaster(0, 0, 52, 10, b"\xFF" * (52 * 10)) + ds.BuildOverviews("NEAR", [2]) + ds.GetRasterBand(1).GetOverview(0).Fill( + 127 + ) # to demonstrate we ignore overviews unfortunately + ds = None + + filename2 = str(tmp_vsimem / "in2.tif") + ds = gdal.GetDriverByName("GTiff").Create( + filename2, 100, 10, 2, options=["ALPHA=YES"] + ) + ds.SetGeoTransform([0, 1, 0, 0, 0, -1]) + ds.GetRasterBand(1).WriteRaster(48, 0, 52, 10, b"\xF0" * (52 * 10)) + ds.GetRasterBand(2).WriteRaster(48, 0, 52, 10, b"\xFF" * (52 * 10)) + ds.BuildOverviews("NEAR", [2]) + ds.GetRasterBand(1).GetOverview(0).Fill( + 127 + ) # to demonstrate we ignore overviews unfortunately + ds = None + + vrt_filename = str(tmp_vsimem / "test.vrt") + gdal.BuildVRT(vrt_filename, [filename1, filename2], resampleAlg=gdal.GRIORA_Cubic) + + ds = gdal.Open(vrt_filename) + assert ds.ReadRaster( + buf_xsize=10, buf_ysize=1 + ) == b"\xFF\xFF\xFF\xFF\xFC\xF0\xF0\xF0\xF0\xF0" + (b"\xFF" * 10) + + vrt_filename = str(tmp_vsimem / "test.vrt") + gdal.BuildVRT( + vrt_filename, [filename1, filename2], resampleAlg=gdal.GRIORA_Bilinear + ) + + ds = gdal.Open(vrt_filename) + assert ds.ReadRaster( + buf_xsize=10, buf_ysize=1 + ) == b"\xFF\xFF\xFF\xFF\xFB\xF1\xF0\xF0\xF0\xF0" + (b"\xFF" * 10) diff --git a/frmts/vrt/vrtsourcedrasterband.cpp b/frmts/vrt/vrtsourcedrasterband.cpp index cc3c0e6f6698..05bd3111979c 100644 --- a/frmts/vrt/vrtsourcedrasterband.cpp +++ b/frmts/vrt/vrtsourcedrasterband.cpp @@ -124,6 +124,25 @@ bool VRTSourcedRasterBand::CanIRasterIOBeForwardedToEachSource( GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize, int nBufXSize, int nBufYSize, GDALRasterIOExtraArg *psExtraArg) const { + const auto IsNonNearestInvolved = [this, psExtraArg] + { + if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour) + { + return true; + } + for (int i = 0; i < nSources; i++) + { + if (papoSources[i]->GetType() == VRTComplexSource::GetTypeStatic()) + { + auto *const poSource = + static_cast(papoSources[i]); + if (!poSource->GetResampling().empty()) + return true; + } + } + return false; + }; + // If resampling with non-nearest neighbour, we need to be careful // if the VRT band exposes a nodata value, but the sources do not have it. // To also avoid edge effects on sources when downsampling, use the @@ -131,7 +150,7 @@ bool VRTSourcedRasterBand::CanIRasterIOBeForwardedToEachSource( // nominal resolution, and then downsampling), but only if none of the // contributing sources have overviews. if (eRWFlag == GF_Read && (nXSize != nBufXSize || nYSize != nBufYSize) && - psExtraArg->eResampleAlg != GRIORA_NearestNeighbour && nSources != 0) + nSources != 0 && IsNonNearestInvolved()) { bool bSourceHasOverviews = false; const bool bIsDownsampling = (nBufXSize < nXSize && nBufYSize < nYSize); @@ -148,6 +167,27 @@ bool VRTSourcedRasterBand::CanIRasterIOBeForwardedToEachSource( VRTSimpleSource *const poSource = static_cast(papoSources[i]); + if (poSource->GetType() == VRTComplexSource::GetTypeStatic()) + { + auto *const poComplexSource = + static_cast(poSource); + if (!poComplexSource->GetResampling().empty()) + { + const int lMaskFlags = + const_cast(this) + ->GetMaskFlags(); + if ((lMaskFlags != GMF_ALL_VALID && + lMaskFlags != GMF_NODATA) || + IsMaskBand()) + { + // Unfortunately this will prevent using overviews + // of the sources, but it is unpractical to use + // them without serious implementation complications + return false; + } + } + } + double dfXOff = nXOff; double dfYOff = nYOff; double dfXSize = nXSize; @@ -388,9 +428,41 @@ CPLErr VRTSourcedRasterBand::IRasterIO( // recursion l_poDS->SetEnableOverviews(false); } + + const auto eResampleAlgBackup = psExtraArg->eResampleAlg; + if (psExtraArg->eResampleAlg == GRIORA_NearestNeighbour) + { + std::string osResampling; + for (int i = 0; i < nSources; i++) + { + if (papoSources[i]->GetType() == + VRTComplexSource::GetTypeStatic()) + { + auto *const poComplexSource = + static_cast(papoSources[i]); + if (!poComplexSource->GetResampling().empty()) + { + if (i == 0) + osResampling = poComplexSource->GetResampling(); + else if (osResampling != + poComplexSource->GetResampling()) + { + osResampling.clear(); + break; + } + } + } + } + if (!osResampling.empty()) + psExtraArg->eResampleAlg = + GDALRasterIOGetResampleAlg(osResampling.c_str()); + } + const auto eErr = GDALRasterBand::IRasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg); + + psExtraArg->eResampleAlg = eResampleAlgBackup; l_poDS->SetEnableOverviews(bBackupEnabledOverviews); return eErr; }