gdalwarp and cubic spline resampling: broken nodata handling; cubic resamp is Ok

[maphew]

There is a bug in gdalwarp and the Cubic Spline resampler. See rcs-screen.png for the broken output and rc-screen.png for non-broken using the same input files. They are screen captures from Openev showing the resultant mosaics from gdalwarp, with "Windowed Raster Re-enhancement" applied after zooming in to bottom left corner.

The files in the bug package were processed with fwtools 1.2.2, but I've replicated the problem with "GDAL 1.5dev, FWTools 1.3.0, released 2007/04/06"

---

Commands Used

 gdalwarp --version
GDAL 1.4.0.0, FWTools 1.2.2, released 2007/02/22

Mosaic using Cubic Spline resampling, downsize output 90% for bug report

 gdalwarp -rcs -tr 30 -multi dem_115b_30m ak_115b_dem dem_gd-merge-rcs.tif 
 gdal_translate -outsize 10% 10% -co compress=lzw dem_gd-merge-rcs.tif \
 gdal-rcs-bug\rcs-out.tif 

Mosaic using Cubic resampling, downsize 90%

 gdalwarp -rc -tr 30 -multi dem_115b_30m ak_115b_dem dem_gd-merge-rc.tif 
 gdal_translate -outsize -co compress=lzw 10% 10% dem_gd-merge-rcs.tif \
 gdal-rcs-bug\rc-out.tif 

the "version" field in Ticket Properties doesn't have 1.5 listed

[maphew]

Screen dump of cubic resamp, it's okay

[maphew]

Screen dump of cubic spline resamp, it's broken

[maphew]

package of data and commands to replicate the bug

[maphew]

.tif files in the attached package were reduced a further 25% in order to get under the 1mb attach limit.

[warmerdam]

Andrey,

I've added you as a cc: in case this might relate to the work you have done recently, or perhaps you would like to look into it.

[dron]

The problem arises when the density component calculated. It is due to the floating point computation errors. GWKCubicSplineResample() returns dfDensity value != 1.0 even if the all input densities are equal to 1.0. Sometimes it is slightly less than 1.0, sometimes greater than 1.0. When dfDensity is less than 1.0 GWKSetPixelValue() applies density mask and we can see results above.

Possible solutions:

1. Test all values in the filter window for equality. If values are the same then no filtering should be applied. This solution increases computation time.
2. Rework resampling logic to avoid this problem at all. I have no idea how it can be done.

Frank, do you have any thoughts?

Best regards, Andrey

[warmerdam]

Andrey,

I'm sorry for dropping the ball on this. I didn't notice your request for thoughts.

I haven't examined this very closely, but I imagine we will need some sort of epsilon test to see if we are very close to a density of 1.0 and if so to treat it as exactly 1.0.

[warmerdam]

Taking this one back to work on the epsilon handling.

Matt, I downloaded the gdalwarp-rcs-bug.zip file but it isn't clear to me whether it is enough to reproduce the problem. Don't I need the input data you use in the gdalwarp commands to reproduce the problem?

I'd prefer to reproduce the result myself before trying to fix.

[maphew]

The .zip has everything necessary to reproduce the bug, just add .tif to the input files. Sorry for being unclear. The commands posted earlier were from before putting the package together.

gdalwarp -rc -tr 30 -multi dem_115b_30m.tif ak_115b_dem.tif dem_gd-merge-rc.tif
gdalwarp -rcs -tr 30 -multi dem_115b_30m.tif ak_115b_dem.tif dem_gd-merge-rcs.tif

[Even Rouault]

Duplicated in #1945

[Even Rouault]

The commandline is incorrect indeed. -tr requires 2 values, not 1. I've added a check for that in trunk in r12553. However this has nothing to do with the core of the problem.

[Even Rouault]

To reproduce the problem, I used the following command line:

gdalwarp -srcnodata -3.4028234663852886e+38  -dstnodata -3.4028234663852886e+38 -rcs gdalwarp-rcs-bug/dem_115b_30m.tif gdalwarp-rcs-bug/ak_115b_dem.tif dem_gd-merge-rc.tif

The following patch suggested by Frank improves the situation a lot :

Index: alg/gdalwarpkernel.cpp
===================================================================
--- alg/gdalwarpkernel.cpp      (révision 12552)
+++ alg/gdalwarpkernel.cpp      (copie de travail)
@@ -703,7 +703,7 @@
 /*      existing destination value, and mix it with the source to       */
 /*      get the new "to apply" value.  Also compute composite density.  */
 /* -------------------------------------------------------------------- */
-    if( dfDensity < 1.0 )
+    if( dfDensity < 1.0 - 1e-5 )
     {
         double dfDstReal, dfDstImag, dfDstDensity = 1.0;

However, this is not sufficient. There are still problems at the edges of the source images where the computation at line 783 of alg/gdalwarpkernel.cpp tries to average the nodata value (dfDstReal, dfDstDensity=1.0) with a valid value (dfReal, dfDensity=.95), leading to a big negative number. There's probably something missing on the destination to test for nodata values.

We do not observe this problem when using -rc (cubic sampling) since when one of the source pixel is nodata, GWKCubicResample returns FALSE. The side effect is that on the above test case, we have an edge of nodata value between the two source images...

[warmerdam]

OK, I'm seeing a few issues here. The -rc has some funky overshoot effects on edges which can produce a white fringe in the output in a case like this:

gdalwarp -ot Int16 -rc dem_115b_30m.tif out.tif

This seems to be handled properly with -rc if the source nodata value is known. eg.

gdalwarp -ot Int16 -rc dem_115b_30m.tif out.tif -srcnodata -32768

It appears that gdalwarp ignores nodata values set in the source file, so they have to be specified on the commandline.

On the other hand, as in Even's case above, the -rcs interpolator does not seem to properly honour nodata settings, and ends up mixing in values as demonstrated with:

gdalwarp -ot Int16 -rc dem_115b_30m.tif out.tif -srcnodata -32768

I believe that the earlier discussions of epsilon and so forth are missing the main point that some interpolators are ignoring source nodata masking when collecting values for interpolation. I am going to treat that as the main issue of this bug report.

[warmerdam]

I have created #1949 to address issue of source nodata being ignored unless specified on the commandline for gdalwarp.

[warmerdam]

Actually, the command to demonstrate the problem with the spline interpolator is:

gdalwarp -ot Int16 -rcs dem_115b_30m.tif out.tif -srcnodata -32768

[warmerdam]

Problem corrected in trunk (r12558) with this comment:

Fixed Cubic Spline and Lanczos. interpolators (functions GWKCubicSplineResample() and GWKLanczosResample()) to renormalize the results if the accumulated weights are less than 1.0 due to some no-data pixels in the kernel. (#1610)

Also back ported into 1.4 branch as r12559.

Patch also fixes a different problem with cubic resampler on edges near nodata.

[warmerdam]

OK, some additional problems. Mixing with background nodata values was occuring and caused dramatic results when operating in full floating point. r12560 in trunk avoid this for the most part by not doing mixing when dfDensity is very near 1.0 in GWKSetPixel().

But, when mosaicing the two files there were still some pixels getting mixed due to densities further from 1 (0.996) and so digging deeper I see that panDstValid[] was not getting created when -dstnodata was set, nor was it being checked in GWKSetPixel(). Both issues corrected as r12561 in trunk.

Reopening pending merging these back into 1.4 branch.

[warmerdam]

Changes merged into 1.4 branch as r12562. Closing again.

[maphew]

thank you frank and roualt!

[warmerdam]

Note, this change was botched by me, and introduced #1964 in 1.4.3.