Skip to content

Commit

Permalink
Merge pull request #11585 from rouault/fix_11583
Browse files Browse the repository at this point in the history 
VRT: avoid artifacts on source boundaries when the VRT has an alpha or mask band, and non-nearest neighbour resampling is involved
  • Loading branch information
@rouault
rouault authored Jan 9, 2025
2 parents eb84213 + 4dc6363 commit d4c58d2
Show file tree
Hide file tree
Showing 2 changed files with 183 additions and 1 deletion.
110 changes: 110 additions & 0 deletions autotest/gcore/vrt_read.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
74 changes: 73 additions & 1 deletion frmts/vrt/vrtsourcedrasterband.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,14 +124,33 @@ 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<VRTComplexSource *>(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
// base implementation of IRasterIO() (that is acquiring sources at their
// 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);
Expand All @@ -148,6 +167,27 @@ bool VRTSourcedRasterBand::CanIRasterIOBeForwardedToEachSource(
VRTSimpleSource *const poSource =
static_cast<VRTSimpleSource *>(papoSources[i]);

if (poSource->GetType() == VRTComplexSource::GetTypeStatic())
{
auto *const poComplexSource =
static_cast<VRTComplexSource *>(poSource);
if (!poComplexSource->GetResampling().empty())
{
const int lMaskFlags =
const_cast<VRTSourcedRasterBand *>(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;
Expand Down Expand Up @@ -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<VRTComplexSource *>(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;
}
Expand Down

0 comments on commit d4c58d2

Please sign in to comment.