gdalwarp with lanczos creates aliased result

zdila commented 1 month ago

What is the bug?

Gdalwarp with -r lanczos can cause aliased edges.

lanczos:

cubic:

I have three input vrt files but I can reproduce it only on single one.

Steps to reproduce the issue

gdalwarp -s_srs 'EPSG:8353' -t_srs 'EPSG:3857' -tr 0.29858214173896972949 0.29858214173896972949  -r lanczos -of GTiff -co COMPRESS=WEBP -co WEBP_LEVEL=85 -co NUM_THREADS=ALL_CPUS -wo NUM_THREADS=ALL_CPUS -multi -te 2118692.6 6294091.5 2120157.9 6295276.7 Ortofoto_2021_stred_jtsk_rgb/all.vrt stred-extract-warped.tif

Versions and provenance

Tried with 2 versions:

Debian linux testing: GDAL 3.9.2, released 2024/08/13
Official osged/gdal Docker image: GDAL 3.9.3, released 2024/10/07

Additional context

Tried also with no success:

JPEG instead fo WEBP
removing -co TILED=YES -co BIGTIFF=YES
removing alpha channel from vrt

rouault commented 1 month ago

Gdalwarp with -r lanczos can cause aliased edges.

you'll have to explain more the issue with some highlighting to show where you see something wrong. I can't even notice a difference with your 2 screeshots... And provide the input dataset

zdila commented 1 month ago

As animated gif:

aliasing

And provide the input dataset

It has 160GB so I will provide a description where to download it and create it.

rouault commented 1 month ago

It has 160GB so I will provide a description where to download it and create it.

you can also try to use gdal_translate -projwin or -srcwin to create a minimum reproducer. I don't have 160 GB available on my hard disk!

zdila commented 1 month ago

If I extract small area and transform it, I can't reproduce the issue.

gdal_translate -projwin -422444.0 -1185375.5 -422211.0 -1185563.9 Ortofoto_2021_stred_jtsk_rgb/all.vrt extracted.tif
gdalwarp -s_srs 'EPSG:8353' -t_srs 'EPSG:3857' -tr 0.29858214173896972949 0.29858214173896972949 -r lanczos -of GTiff -co COMPRESS=WEBP -co WEBP_LEVEL=85 -co NUM_THREADS=ALL_CPUS -wo NUM_THREADS=ALL_CPUS -multi extracted.tif stred-extract-warped.tif

zdila commented 1 month ago

Anyway, to reproduce it:

download 2 zipfiles Central part 2021 (RGB) from https://www.geoportal.sk/en/zbgis/orthophotomosaic/2nd-cycle/
extract them to single directory
run gdalbuildvrt -a_srs EPSG:8353 -r lanczos all.vrt *.tif
run gdalwarp -s_srs 'EPSG:8353' -t_srs 'EPSG:3857' -tr 0.29858214173896972949 0.29858214173896972949 -r lanczos -of GTiff -co COMPRESS=WEBP -co WEBP_LEVEL=85 -co NUM_THREADS=ALL_CPUS -wo NUM_THREADS=ALL_CPUS -multi -te 2118692.6 6294091.5 2120157.9 6295276.7 all.vrt stred-extract-warped.tif

zdila commented 1 month ago

Created small reproducible testcase:

download https://drive.google.com/file/d/13mmQS9nqyLY0D0oa5c0rdJNESO8wdZ_T/view?usp=sharing
run gdalwarp -s_srs 'EPSG:8353' -t_srs 'EPSG:3857' -tr 0.29858214173896972949 0.29858214173896972949 -r lanczos -of GTiff -co COMPRESS=WEBP -co WEBP_LEVEL=85 -co NUM_THREADS=ALL_CPUS -wo NUM_THREADS=ALL_CPUS -multi extracted.tif warped.tif

jratike80 commented 1 month ago

There is lot of extra in the command. This simplified one makes almost similar result: gdalwarp -s_srs EPSG:8353 -t_srs EPSG:3857 -r lanczos -of GTiff extracted.tif warped2.tif

I think that maybe the special input SRS EPSG:8353 plays some role. For testing the hypothesis I tried the following command and the result looks better to my eyes.

gdalwarp -s_srs EPSG:3067 -t_srs EPSG:3857 -r lanczos -of GTiff extracted.tif warped3.tif

zdila commented 1 month ago

Now please try with https://drive.google.com/file/d/1KOAjpfTy98o0eiVEnLJgnw3BPAOOW0bv/view?usp=sharing and same command and the issue will not reproduce. I wonder what difference between these two files could cause that aliasing.

rouault commented 1 month ago

and the issue will not reproduce.

actually it does, for example when comparing the cubicspline vs lanczos around coordinate 2365033.68, 6240237.04 in EPSG:3857

rouault commented 1 month ago

I've made a significant breakthrough in debugging this. The following patch disabling the Lanczos optimized code path (used since 12 years...) seems to fix the aliasing issue.

diff --git a/alg/gdalwarpkernel.cpp b/alg/gdalwarpkernel.cpp
index 565a6657a1..8f687abbd4 100644
--- a/alg/gdalwarpkernel.cpp
+++ b/alg/gdalwarpkernel.cpp
@@ -3566,7 +3566,7 @@ static GWKResampleWrkStruct *GWKResampleCreateWrkStruct(GDALWarpKernel *poWK)
     psWrkStruct->padfRowImag =
         static_cast<double *>(CPLCalloc(nXDist, sizeof(double)));

-    if (poWK->eResample == GRA_Lanczos)
+    if (0 && poWK->eResample == GRA_Lanczos)
     {
         psWrkStruct->pfnGWKResample = GWKResampleOptimizedLanczos;

zdila commented 1 month ago

actually it does, for example when comparing the cubicspline vs lanczos around coordinate 2365033.68, 6240237.04 in EPSG:3857

Hmm... In my case I see no problem there. I compare warped2.tif with extracted2.tif and the only relevant difference i see there is a global 1px offset.

Please post a screenshot so I could see your result. Thank you.

Mine is like this:

rouault commented 1 month ago

Please post a screenshot so I could see your result

Lanczos:

vs cubic spline:

zdila commented 1 month ago

Thank you. Strange that I can't reproduce it with that extracted2.tif file (and you can). I used Debian linux testing: GDAL 3.9.2, released 2024/08/13. It would be more logical for me if I could reproduce it on both files as you can. Will try more tomorrow.

rouault commented 1 month ago

Fix in #11046

rouault commented 1 month ago

Workaround: if just reprojecting without trying to downsample/oversample, adding -wo XSCALE=1 -wo YSCALE=1 will avoid the issue

zdila commented 1 month ago

Wow! That's a lot of code!

I've compiled your fork/branch (GDAL 3.11.0dev-7727dbaa84, released 2024/10/18) and for me it works. Also the workaround (-wo XSCALE=1 -wo YSCALE=1) does with the distro gdal.

Big thank you!

OSGeo / gdal