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 int interpolationTypeId,
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);
}//if
else
{
// strict bilinear implementation with 3rd pass and inline NODATA detection with NN / horizontal linear fallback
// generally, looks worse and performs about the same
gbImage_BilinearVectorized(src,
srcWidth,
srcHeight,
resampleBuffer,
resampleBufferWidth,
resampleBufferHeight,
width*height-1,
(srcWidthUnpadded * srcHeightUnpadded - srcDataCount),
tempBuffer,
tempBufferBytes,
doNotTile);
}//else
if (useBitmask || !cleanOutput)
{
_gbImage_FilterWithBitmaskByRow_Any(cleanOutput ? dest : resampleBuffer,
resampleBufferWidth,
useBitmask ? tempBuffer : NULL,
useBitmask ? destHeight : 0,
dest,
destWidth,
destHeight);
}//if
}//gbImage_GetZoomedTile_BilinearPunch_v2
static inline void _gbImage_FillNODATA(float* src,
const size_t srcWidth,
const size_t srcHeight,
const size_t srcWidthUnpadded,
const size_t srcHeightUnpadded,
const size_t srcDataCount,
float* tempBuffer,
const size_t tempBufferBytes,
const int interpolationTypeId,
const float NODATA)
{
if ( srcDataCount >= srcWidth * srcHeight
|| interpolationTypeId == kGB_TileEngine_Interp_NN
|| interpolationTypeId == kGB_TileEngine_Interp_BilinearCareful)
{
return;
}//if
size_t cellSize;
if (interpolationTypeId == kGB_TileEngine_Interp_Lanczos3x3 || interpolationTypeId == kGB_TileEngine_Interp_Lanczos5x5)
{
const size_t largestDiff = MAX(srcWidth - srcWidthUnpadded, srcHeight-srcHeightUnpadded);
cellSize = MAX(MIN(largestDiff * 2 + 2, srcWidth), 1);
gbImage_ExtendEdgesRightAndDownOnly(src, srcWidth, srcHeight, srcWidthUnpadded, srcHeightUnpadded);
}//if
else
{
cellSize = 3;
}//else
gbImage_RasterNeighborhood_Mean(src, (int)cellSize, srcWidth*srcHeight, 0, 0, srcWidth-1, srcHeight-1, srcWidth, srcHeight, true, tempBuffer, tempBufferBytes);
}//_gbImage_FillNODATA
// this did not vectorize well, at all.
void gbImage_RasterNeighborhood_Mean(float* src,
const int cellSize,
const size_t dataCount,
const size_t rwx0,
const size_t rwy0,
const size_t rwx1,
const size_t rwy1,
const size_t width,
const size_t height,
const bool smoothing,
float* tempBuffer,
const size_t tempBufferBytes)
{
const size_t n = width * height;
size_t dataFound = 0;
bool didMalloc = false;
if (tempBuffer == NULL)
{
didMalloc = true;
tempBuffer = malloc(sizeof(float) * n);
}//if
size_t _cellSize = (size_t)cellSize;
_cellSize = _cellSize > width && _cellSize > height ? MAX(width, height) : _cellSize;
const size_t cellOffsetBack = _cellSize % 2 == 0 ? (_cellSize >> 1) - 1 : (_cellSize - 1) >> 1;
const size_t cellOffsetFwd = _cellSize % 2 == 0 ? _cellSize >> 1 : (_cellSize - 1) >> 1;
memcpy(tempBuffer, src, sizeof(float) * n);
size_t startAX;
size_t startAY;
size_t endAX;
size_t endAY;
size_t y_width;
size_t aY_width;
size_t y;
size_t x;
size_t aY;
size_t aX;
float f_mean;
float f_n;
// don't use additional smoothing if windowing
const bool useSmoothing = smoothing && rwx0 == 0 && rwy0 == 0 && rwx1 == (width-1) && rwy1 == (height-1);
for (y = rwy0; y <= rwy1; y++)
{
startAY = y >= cellOffsetBack ? y - cellOffsetBack : rwy0;
endAY = y <= rwy1 - cellOffsetFwd ? y + cellOffsetFwd : rwy1;
y_width = y * width;
for (x = rwx0; x <= rwx1; x++)
{
if (tempBuffer[y_width+x] != 0.0F)
{
dataFound++;
startAX = x >= cellOffsetBack ? x - cellOffsetBack : rwx0;
endAX = x <= rwx1 - cellOffsetFwd ? x + cellOffsetFwd : rwx1;
f_mean = 0.0F;
f_n = 0.0F;
// 1. get mean value
for (aY = startAY; aY <= endAY; aY++)
{
aY_width = aY * width;
for (aX = startAX; aX <= endAX; aX++)
{
if (tempBuffer[aY_width + aX] != 0.0F)
{
f_mean += tempBuffer[aY_width + aX];
f_n += 1.0F;
}//if
else if (useSmoothing && src[aY_width + aX] != 0.0F) // optional smoothing
{
f_mean += src[aY_width + aX];
f_n += 1.0F;
}//else if
}//for AX
}//for AY
f_mean = f_n != 0.0F ? f_mean / f_n : 0.0F;
// 2. set NODATAs to mean value
for (aY = startAY; aY <= endAY; aY++)
{
aY_width = aY * width;
for (aX = startAX; aX <= endAX; aX++)
{
if (tempBuffer[aY_width + aX] == 0.0F)
{
src[aY_width + aX] = f_mean;
}//if
}//for AX
}//for AY
}//if
}//for X
if (dataFound >= dataCount)
{
break;
}//if
}//for Y
if (didMalloc)
{
free(tempBuffer);
tempBuffer = NULL;
}//if
}//gbImage_RasterNeighborhood_Mean
// useful for extending edges for Lanczos kernel resampling of continuous
// raster data without needing to load another tile. mitigates, but does
// not eliminate, tile boundary artifacts in such a case.
void gbImage_ExtendEdgesRightAndDownOnly(float* src,
const size_t width,
const size_t height,
const size_t dataWidth,
const size_t dataHeight)
{
size_t x;
size_t y;
size_t y_width;
size_t y_width_x;
float dataVal;
if (dataWidth > 0 && width > dataWidth)
{
for (y=0; y<dataHeight; y++)
{
y_width = y * width;
dataVal = src[y_width + dataWidth - 1];
for (x=dataWidth; x<width; x++)
{
y_width_x = y_width + x;
if (src[y_width_x] == 0.0F)
{
src[y_width_x] = dataVal;
}//if
}//for
}//for
}//if
if (dataHeight > 0 && height > dataHeight)
{
size_t ysb1_width;
for (y=dataHeight; y<height; y++)
{
y_width = y * width;
ysb1_width = (y-1) * width;
for (x=0; x<width; x++)
{
y_width_x = y_width + x;
if (src[y_width_x] == 0.0F)
{
src[y_width_x] = src[ysb1_width + x];
}//if
}//for
}//for
}//if
}//ExtendEdgesRightAndDownOnly
// this both filters resampled output with a NODATA mask of 0.0F,
// and copies rows into dest from src which are typically padded
// for resampling due to the kernel extent.
//
// this is for bilinear interpolation -- Lanczos (or cubic) needs
// a clip or some analog to clean up data values turned into NODATAs
// by the negative kernel lobe, or to clamp to the original data range
// if desired.
static inline void _gbImage_FilterWithBitmaskByRow_Any(const float* src,
const size_t srcWidth,
const float* bitmask,
const size_t bitmaskWidth,
float* dest,
const size_t destWidth,
const size_t destHeight)
{
const size_t procWidth = bitmask != NULL ? MIN(destWidth, bitmaskWidth) : destWidth;
size_t y_bitmaskWidth;
size_t y_srcWidth;
size_t y_destWidth;
if (bitmask != NULL)
{
for (size_t y=0; y<destHeight; y++)
{
y_bitmaskWidth = y * bitmaskWidth;
y_srcWidth = y * srcWidth;
y_destWidth = y * destWidth;
// oh noes, the Accelerate framework
vDSP_vmul(src + y_srcWidth, 1,
bitmask + y_bitmaskWidth, 1,
dest + y_destWidth, 1,
procWidth);
}//for
}//if
else
{
const size_t rowBytes = sizeof(float) * destWidth;
for (size_t y=0; y<destHeight; y++)
{
y_srcWidth = y * srcWidth;
y_destWidth = y * destWidth;
memcpy(dest + y_destWidth, src + y_srcWidth, rowBytes);
}//for
}//else
}//_gbImage_FilterTargetWithBitmaskByRow_Any_______________________________________________
Qgis-developer mailing list
Qgis-developer@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/qgis-developer
