Re: [Qgis-developer] Raster resampling suggestion

2014-02-18 Thread Nick Dolezal
Hey Marco,

Here's the code I have for a bilinear-based implementation (attached),
though I've omitted the actual NN and bilinear resampling code itself for
brevity.  This does have a number of implementation-specific assumptions
which would need to be reworked for the general case.  (the general case in
GIS is hard and makes my brain bleed.)  So it's not really close to the
point of being ready for a pull request.

It could be adapted to work on RGBA output, though with some modifications
mostly involving using the alpha channel with a NN resample in place of the
NODATA values.

It should be all ANSI C with no external dependencies. (except for one
irrelevant Accelerate framework call)

Hope this is useful -- let me know if anything else is needed.


On Tue, Feb 18, 2014 at 1:14 AM, Marco Hugentobler <
marco.hugentob...@sourcepole.ch> wrote:

> Hi Nick
>
> Resampling on data level does not work in case of colormap rasters,
> therefore the resampling filter currently comes after assigning rgba. One
> possibility would be to provide a user choice if resampling should be done
> on color values or on data.
>
> I think it is best if you post your code (as a pull request on github or
> via mailing list), then others can give it a try and provide feedback. I'm
> not really a raster power user myself (mostly working with rasters as
> background images for maps).
>
> Regards,
> Marco
// 
// GetZoomedTile_BilinearPunch:
// 
//
// Wrapper around bilinear resize for a raster that preserves NODATA values
// without artifacts.
//
// Evenly resizes (zooming in only) a planar 32-bit float tile of varying resolution,
//normally 256x256 - 384x384.
//
// It will work on other widths and heights, but some of the logic expecting evenly
//divisible srcWH and resampleBufferWH may need to change, or need a ceilf.
//
// Preserves NODATA values vaguely (but not really) derived from PVR textures; 
//a "punch through" technique is used based upon NN resampling.
//
// To avoid interpolating NODATAs, a neighborhood raster stats mean function is used.
//
// Currently assumes NODATA = 0.0F and that interpolated values will not be transformed 
//into NODATA.
//
void gbImage_GetZoomedTile_BilinearPunch_v2(float*   src,
const size_t srcWidth,
const size_t srcHeight,
const size_t srcWidthUnpadded,
const size_t srcHeightUnpadded,
const size_t srcDataCount,
float*   dest,
const size_t destWidth,
const size_t destHeight,
float*   resampleBuffer,
const size_t resampleBufferWidth,
const size_t resampleBufferHeight,
float*   tempBuffer,
const size_t tempBufferBytes,
const intinterpolationTypeId,
const bool   doNotTile,
const float  NODATA)
{
const bool useBitmask  = srcDataCount < srcWidth * srcHeight;
const bool cleanOutput = destWidth == resampleBufferWidth && destHeight == resampleBufferHeight;

// first, get NN resampled mask to filter resampled output with later
if (useBitmask && interpolationTypeId == kGB_TileEngine_Interp_Bilinear)
{
	// returns NN resampled raster tile of 0 or 1, which is somewhat more efficient than using conditionals later, but will *not* work in the general case.
// for the general case, do a normal NN resize here, and use a conditional in the final filtering step instead of vmul
gbImage_GetZoomedTile_NN_FromCrop_Bitmask(src, srcWidthUnpadded, srcHeightUnpadded, srcWidth, dest, destWidth, destHeight, NODATA);

_gbImage_FillNODATA(src, srcWidth, srcHeight, srcWidthUnpadded, srcHeightUnpadded, srcDataCount, tempBuffer, tempBufferBytes, interpolationTypeId, NODATA);
}//if

// always do NoNODATA if possible, other version is branchy and slower due to inline NODATA conditionals
if (interpolationTypeId == kGB_TileEngine_Interp_Bilinear || !useBitmask)
{
gbImage_BilinearVectorized_NoNODATA(src,
srcWidth,
srcHeight,
cleanOutput ? dest : resampleBuffer,
resampleBufferWidth,
resampleBufferHeight,
true);
   

[Qgis-developer] Raster resampling suggestion

2014-02-17 Thread 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!
___
Qgis-developer mailing list
Qgis-developer@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/qgis-developer