[Qgis-developer] Raster resampling suggestion

[Nick Dolezal]

Hey,

I've been using qgis and have been pretty impressed with it so far.  I would however like to offer a suggestion for improvement of the raster resampling.

Right now, it appears the resampling is done on the final RGBA output, rather than on the data values.  Further, there are issues with both negative kernel lobe artifacts from the cubic convolution kernel and NODATA interpolation artifacts from the bilinear resampling.

My suggestion on how to deal with this is:
1.  Resample the data values, not the RGBA output
2.  To remove the NODATA artifacts, create two rasters while resampling:
2a.  First, create a nearest-neighbor only resampled image from the original.
2b.  Second, before interpolating, use a neighborhood mean that only fills NODATA cells, edge extend, or combination of the two to prevent NODATA interpolation artifacts.  This must be done to a minimum of (kernelExtent - 1.0) cells.  Then, interpolate normally.
3.  Use the NN raster as a mask to gate data values from the resampled image.  ie:  dest[i] = nnMask[i] != NODATA ? cubicResampledSrc[i] : NODATA;  Should be 2-4 cycles per pixel SIMD depending on dual issue.  If the NODATA value is 0, you can just clip everything to either 0 or 1 in the NN raster and do lots of nice memory-aligned vmuls.
3B. If using cubic (or Lanczos), you'll need to integrate or do another clipping pass to make sure the evil negative kernel lobe fairy didn't turn any data into NODATA.
4.  If you're maintaining the alpha channel at this point, I would suggest always using NN to resize it.

I have implemented this myself, and it works quite well, for both bilinear and Lanczos interpolation.  Not terribly creative, it's derived from how the alpha channel in PowerVR textures work.  I did not find the neighborhood mean to significantly distort anything, compared to a very strict bilinear implementation that never interpolated NODATA but rather fell back on NN.

Net results are:
1.  Higher quality interpolation, if the source raster uses a data type with more width than uint8_t.  You can actually zoom in crazy amounts on even bilinear interpolated rasters as long as you're interpolating something like 32-bit floats.
2.  No interpolation into NODATA values.
3.  With bicubic, elimination of negative kernel lobe black ringing artifacts in NODATA areas.  (ringing artifacts will be present to some degree still.  This can mitigated by transforming to the logarithmic domain and back but that's relatively expensive)
4.  Possibly faster performance if you were interpolating 4 channels of RGBA interleaved data which could be done with 2x planar interpolations instead.

If it would help, I can provide source code of a working implementation and release it without copyright, although as this stuff is simple but implementation-dependent, I'm not sure it'd be more useful than the above.  Just pad any NODATA around the data values, interpolate and gate / clip with the 2nd NN raster.

Thanks!