Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

VRT: avoid artifacts on source boundaries when the VRT has an alpha or mask band, and non-nearest neighbour resampling is involved #11585

Merged
merged 1 commit into from
Jan 9, 2025
Merged

VRT: avoid artifacts on source boundaries when the VRT has an alpha or mask band, and non-nearest neighbour resampling is involved #11585

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
110 changes: 110 additions & 0 deletions autotest/gcore/vrt_read.py
View file
Open in desktop
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
Open in desktop
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
Loading