Index: /trunk/psLib/src/imageops/psImageInterpolate.c =================================================================== --- /trunk/psLib/src/imageops/psImageInterpolate.c (revision 20305) +++ /trunk/psLib/src/imageops/psImageInterpolate.c (revision 20306) @@ -7,6 +7,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-10-16 22:55:45 $ + * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-22 02:10:37 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -31,46 +31,215 @@ #include "psImageInterpolate.h" -static void imageInterpolateOptionsFree(psImageInterpolateOptions *options) +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Interpolation kernels +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Kernel sizes as a function of interpolation mode: probably faster than a 'switch' +static int kernelSizes[] = { 0, // NONE + 0, // FLAT + 2, // BILINEAR + 3, // BIQUADRATIC + 3, // GAUSS + 4, // LANCZOS2 + 6, // LANCZOS3 + 8 // LANCZOS4 +}; + +/// Generate a linear interpolation kernel +static inline void interpolationKernelBilinear(psF32 *kernel, // Kernel vector to populate + float frac // Fraction of pixel + ) +{ + kernel[0] = 1.0 - frac; + kernel[1] = frac; +} + +/// Generate a biquadratic interpolation kernel +/// This reduces Gene's original made-up 2D kernel to 1D +static inline void interpolationKernelBiquadratic(psF32 *kernel, // Kernel vector to populate + float frac // Fraction of pixel + ) +{ + float x = frac / 6.0; + float xx = frac * x; + kernel[0] = -x + xx; + kernel[1] = 1.0/3.0 - 2.0 * xx; + kernel[2] = x + xx; +} + +#if 0 +/// Generate a bicubic interpolation kernel +/// "cubic Hermite spline" has a=-1/2 for: +/// W(x) = 1 when x == 0 +/// = (a+2)|x|^3 - (a+3)|x|^2 + 1 when 0 < |x| < 1 +/// = a|x|^3 - 5a|x|^2 + 8a|x| - 4a when 1 < |x| < 2 +/// = 0 otherwise +static inline void interpolationKernelBicubic(psF32 *kernel, // Kernel vector to populate + float frac // Fraction of pixel + ) +{ + float frac2 = PS_SQR(frac); + float frac3 = frac2 * frac; + kernel[0] = 0.5 * frac + frac2 - 0.5 * frac3; // x = -(1 + frac) + kernel[1] = 1.5 * frac3 - 2.5 * frac2 + 1.0; // x = -frac + kernel[2] = 0.5 * frac + 2.0 * frac2 - 1.5 * frac3; // x = 1 - frac + kernel[3] = -0.5 * frac2 - frac3; // x = 2 - frac +} +#endif + +// Generate Lanczos interpolation kernel +static inline void interpolationKernelLanczos(psF32 *kernel, // Kernel vector to populate + int size, // Size of kernel + float frac // Fraction of pixel + ) +{ + if (fabs(frac) < FLT_EPSILON) { + // No real shift + for (int i = 0; i < (size - 1) / 2; i++) { + kernel[i] = 0.0; + } + kernel[(size - 1) / 2] = 1.0; + for (int i = (size - 1) / 2 + 1; i < size; i++) { + kernel[i] = 0.0; + } + return; + } + + float norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos + float norm2 = M_PI * 4.0 / (float)size; // Normalisation for sinc function 1 + float norm3 = M_PI_2 * 4.0 / (float)size; // Normalisation for sinc function 2 + float pos = - (size - 1)/2 - frac; // Position of interest + for (int i = 0; i < size; i++, pos += 1.0) { + if (fabs(pos) < FLT_EPSILON) { + kernel[i] = 1.0; + } else { + kernel[i] = norm1 * sinf(pos * norm2) * sinf(pos * norm3) / PS_SQR(pos); + } + } + return; +} + +// Generate Gauss interpolation kernel +static inline void interpolationKernelGauss(psF32 *kernel, // Kernel vector to populate + int size, // Size of kernel + float sigma, // Width of kernel + float frac // Fraction of pixel + ) +{ + for (int i = 0, pos = - (size - 1) / 2; i < size; i++, pos++) { + kernel[i] = expf(-0.5 / PS_SQR(sigma) * PS_SQR(pos - frac)); + } + // XXX Normalise? + return; +} + + +// Generate interpolation kernel +static inline void interpolationKernel(psF32 *kernel, // Kernel vector to populate + float frac, // Fraction of pixel + psImageInterpolateMode mode // Mode of interpolation + ) +{ + psAssert(frac >= 0.0 && frac < 1.0, "Pixel fraction is %f", frac); + switch (mode) { + case PS_INTERPOLATE_FLAT: + // Nothing to pre-compute + break; + case PS_INTERPOLATE_BILINEAR: + interpolationKernelBilinear(kernel, frac); + break; + case PS_INTERPOLATE_BIQUADRATIC: + interpolationKernelBiquadratic(kernel, frac); + break; + case PS_INTERPOLATE_GAUSS: + interpolationKernelGauss(kernel, kernelSizes[mode], 0.5, frac); + break; + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: + interpolationKernelLanczos(kernel, kernelSizes[mode], frac); + break; + default: + psAbort("Unsupported interpolation mode: %x", mode); + } + + return; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +static void imageInterpolationFree(psImageInterpolation *interp) { // Casting away const - psFree((psImage*)options->image); - psFree((psImage*)options->mask); - psFree((psImage*)options->variance); -} - -psImageInterpolateOptions *psImageInterpolateOptionsAlloc(psImageInterpolateMode mode, - const psImage *image, const psImage *variance, - const psImage *mask, psMaskType maskVal, - double badImage, double badVariance, - psMaskType badMask, psMaskType poorMask, - float poorFrac) -{ - psImageInterpolateOptions *options = psAlloc(sizeof(psImageInterpolateOptions)); // Options, to return - psMemSetDeallocator(options, (psFreeFunc)imageInterpolateOptionsFree); - - options->mode = mode; - // Casting away const to add to options - options->image = psMemIncrRefCounter((psImage*)image); - options->variance = psMemIncrRefCounter((psImage*)variance); - options->mask = psMemIncrRefCounter((psImage*)mask); - options->maskVal = maskVal; - options->badImage = badImage; - options->badVariance = badVariance; - options->badMask = badMask; - options->poorMask = poorMask; - options->poorFrac = poorFrac; - options->shifting = true; - - return options; -} + psFree((psImage*)interp->image); + psFree((psImage*)interp->mask); + psFree((psImage*)interp->variance); + psFree((psImage*)interp->kernel); + psFree((psImage*)interp->kernel2); + psFree((psVector*)interp->sumKernel2); +} + +psImageInterpolation *psImageInterpolationAlloc(psImageInterpolateMode mode, + const psImage *image, const psImage *variance, + const psImage *mask, psMaskType maskVal, + double badImage, double badVariance, + psMaskType badMask, psMaskType poorMask, + float poorFrac, int numKernels) +{ + psImageInterpolation *interp = psAlloc(sizeof(psImageInterpolation)); // Options, to return + psMemSetDeallocator(interp, (psFreeFunc)imageInterpolationFree); + + interp->mode = mode; + // Casting away const to add to interp + interp->image = psMemIncrRefCounter((psImage*)image); + interp->variance = psMemIncrRefCounter((psImage*)variance); + interp->mask = psMemIncrRefCounter((psImage*)mask); + interp->maskVal = maskVal; + interp->badImage = badImage; + interp->badVariance = badVariance; + interp->badMask = badMask; + interp->poorMask = poorMask; + interp->poorFrac = poorFrac; + interp->shifting = true; + + // Pre-calculate interpolation kernels + psImage *kernel = NULL; // Kernel + psImage *kernel2 = NULL; // Kernel^2 + psVector *sumKernel2 = NULL; // Sum of kernel^2 + if (numKernels > 0) { + int size = kernelSizes[mode]; // Size of kernel + + kernel = psImageAlloc(size, numKernels, PS_TYPE_F32); // Kernel + kernel2 = psImageAlloc(size, numKernels, PS_TYPE_F32); // Kernel^2 + sumKernel2 = psVectorAlloc(numKernels, PS_TYPE_F32); // Sum of kernel^2 + + for (int i = 0; i < numKernels; i++) { + float frac = i / (float)numKernels; // Fraction of shift + interpolationKernel(kernel->data.F32[i], frac, mode); + + float sum = 0.0; // Sum of kernel^2 + for (int j = 0; j < size; j++) { + sum += kernel2->data.F32[i][j] = PS_SQR(kernel->data.F32[i][j]); + } + sumKernel2->data.F32[i] = sum; + } + } + interp->kernel = kernel; + interp->kernel2 = kernel2; + interp->sumKernel2 = sumKernel2; + interp->numKernels = numKernels; + + return interp; +} + // Interpolation engine for flat mode (nearest pixel) static inline psImageInterpolateStatus interpolateFlat(double *imageValue, double *varianceValue, psMaskType *maskValue, float x, float y, - const psImageInterpolateOptions *options) + const psImageInterpolation *interp) { // Parameters have been checked by psImageInterpolate() - const psImage *image = options->image; // Image of interest + const psImage *image = interp->image; // Image of interest int xInt = round(x - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in x int yInt = round(y - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in y @@ -81,11 +250,11 @@ // At least one pixel of the interpolation kernel is off the image if (imageValue) { - *imageValue = options->badImage; + *imageValue = interp->badImage; } if (varianceValue) { - *varianceValue = options->badVariance; + *varianceValue = interp->badVariance; } if (maskValue) { - *maskValue = options->badMask; + *maskValue = interp->badMask; } @@ -97,13 +266,13 @@ case PS_TYPE_##TYPE: { \ if (imageValue) { \ - *imageValue = options->image->data.TYPE[yInt][xInt]; \ + *imageValue = interp->image->data.TYPE[yInt][xInt]; \ } \ - if (varianceValue && options->variance) { \ - *varianceValue = options->variance->data.TYPE[yInt][xInt]; \ + if (varianceValue && interp->variance) { \ + *varianceValue = interp->variance->data.TYPE[yInt][xInt]; \ } \ break; \ } - switch (options->image->type.type) { + switch (interp->image->type.type) { FLAT_CASE(U8); FLAT_CASE(U16); @@ -118,11 +287,11 @@ default: psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", - options->image->type.type); + interp->image->type.type); return PS_INTERPOLATE_STATUS_ERROR; } if (maskValue) { - if (options->mask) { - *maskValue = options->mask->data.PS_TYPE_MASK_DATA[yInt][xInt]; + if (interp->mask) { + *maskValue = interp->mask->data.PS_TYPE_MASK_DATA[yInt][xInt]; } else { *maskValue = 0; @@ -133,245 +302,207 @@ } -// Setup for interpolation by kernel -static inline void interpolateKernelSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned - int *xCentral, int *yCentral, // Central pixel of convolution - float x, float y, // Coordinates of interest - psImageInterpolateMode mode // Mode for interpolation - ) -{ - switch (mode) { - case PS_INTERPOLATE_BILINEAR: - *xNum = *yNum = 2; - // Central pixel is the pixel below the point of interest - *xCentral = floor(x - 0.5 + FLT_EPSILON); - *yCentral = floor(y - 0.5 + FLT_EPSILON); - break; - case PS_INTERPOLATE_BICUBE: - case PS_INTERPOLATE_GAUSS: - *xNum = *yNum = 3; - // Central pixel is the closest pixel to the point of interest - *xCentral = x; - *yCentral = y; - break; - case PS_INTERPOLATE_FLAT: - case PS_INTERPOLATE_LANCZOS2: - case PS_INTERPOLATE_LANCZOS3: - case PS_INTERPOLATE_LANCZOS4: - default: - psAbort("Invalid interpolation mode."); - } -} - -// Generate the interpolation kernel -// No need to normalise to unity -static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel - float kernel[xNum][yNum], // Kernel, to be set - bool *xExact, bool *yExact, // Exact shift? - int xCentral, int yCentral, // Central pixel of convolution - float x, float y, // Coordinates of interest - psImageInterpolateMode mode, // Mode for interpolation - bool allowExact // Allow exact shifts? - ) -{ - float xFrac = x - 0.5 - xCentral; // Fraction of pixel in x - float yFrac = y - 0.5 - yCentral; // Fraction of pixel in y - *xExact = (allowExact && fabs(xFrac) < DBL_EPSILON); - *yExact = (allowExact && fabs(yFrac) < DBL_EPSILON); - - switch (mode) { - case PS_INTERPOLATE_BILINEAR: { - kernel[0][0] = (1.0 - xFrac) * (1.0 - yFrac); - kernel[0][1] = xFrac * (1.0 - yFrac); - kernel[1][0] = yFrac * (1.0 - xFrac); - kernel[1][1] = xFrac * yFrac; - break; - } - case PS_INTERPOLATE_BICUBE: { - float xFrac = x - 0.5 - xCentral; // Fraction of pixel in x - float yFrac = y - 0.5 - yCentral; // Fraction of pixel in y - // Calculation variables - double xxFrac = xFrac * xFrac / 6.0; - double yyFrac = yFrac * yFrac / 6.0; - double xyFrac = 0.25 * xFrac * yFrac; - xFrac /= 6.0; - yFrac /= 6.0; - kernel[0][0] = - 1.0/9.0 - xFrac - yFrac + xxFrac + yyFrac + xyFrac; - kernel[0][1] = 2.0/9.0 - yFrac - 2.0 * xxFrac + yyFrac; - kernel[0][2] = - 1.0/9.0 + xFrac - yFrac + xxFrac + yyFrac - xyFrac; - kernel[1][0] = 2.0/9.0 - xFrac + xxFrac - 2.0 * yyFrac; - kernel[1][1] = 5.0/9.0 - 2.0 * xxFrac - 2.0 * yyFrac; - kernel[1][2] = 2.0/9.0 + xFrac + xxFrac - 2.0 * yyFrac; - kernel[2][0] = - 1.0/9.0 - xFrac + yFrac + xxFrac + yyFrac - xyFrac; - kernel[2][1] = 2.0/9.0 + yFrac - 2.0 * xxFrac + yyFrac; - kernel[2][2] = - 1.0/9.0 + xFrac + yFrac + xxFrac + yyFrac + xyFrac; - break; - } - case PS_INTERPOLATE_GAUSS: { - float xFrac = x - xCentral - 0.5; // Fraction of pixel in x - float yFrac = y - yCentral - 0.5; // Fraction of pixel in y - float sigma = 0.5; // Gaussian sigma - for (int j = 0, yPos = - (yNum - 1) / 2; j < yNum; j++, yPos++) { - for (int i = 0, xPos = - (xNum - 1) / 2; i < xNum; i++, xPos++) { - kernel[j][i] = expf(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) + - PS_SQR(yPos - yFrac))); - } - } - break; - } - case PS_INTERPOLATE_FLAT: - case PS_INTERPOLATE_LANCZOS2: - case PS_INTERPOLATE_LANCZOS3: - case PS_INTERPOLATE_LANCZOS4: - default: - psAbort("Invalid interpolation mode."); - } - return; -} // Can't interpolate: kernel is off the image #define INTERPOLATE_OFF() { \ - *imageValue = options->badImage; \ + *imageValue = interp->badImage; \ if (varianceValue) { \ - *varianceValue = options->badVariance; \ + *varianceValue = interp->badVariance; \ } \ if (maskValue) { \ - *maskValue = options->badMask; \ + *maskValue = interp->badMask; \ } \ return PS_INTERPOLATE_STATUS_OFF; \ } -// Set up and check interpolation bounds -// This macro defines many useful values -#define INTERPOLATE_BOUNDS() \ - int xLast = image->numCols - 1, yLast = image->numRows - 1; /* Last pixels of image */ \ - /* Start and stop of kernel on image */ \ - int xStart = xCentral - (xNum - 1) / 2, xStop = xCentral + xNum / 2; \ - int yStart = yCentral - (yNum - 1) / 2, yStop = yCentral + yNum / 2; \ - if (xStart > xLast || xStop < 0 || yStart > yLast || yStop < 0) { \ - /* Interpolation kernel is entirely off the image */ \ - INTERPOLATE_OFF(); \ - } \ - int xMin, xMax, yMin, yMax; /* Minimum and maximum valid pixels on image */ \ - int iMin, iMax, jMin, jMax; /* Minimum and maximum valid pixels on kernel */ \ - bool offImage = false; /* At least one pixel of the kernel is off the image */ \ - if (xStart < 0) { \ - xMin = 0; \ - iMin = -xStart; \ - offImage = true; \ - } else { \ - xMin = xStart; \ - iMin = 0; \ - } \ - if (xStop > xLast) { \ - xMax = xLast; \ - iMax = xLast - xStart + 1; \ - offImage = true; \ - } else { \ - xMax = xStop; \ - iMax = xNum; \ - } \ - if (yStart < 0) { \ - yMin = 0; \ - jMin = -yStart; \ - offImage = true; \ - } else { \ - yMin = yStart; \ - jMin = 0; \ - } \ - if (yStop > yLast) { \ - yMax = yLast; \ - jMax = yLast - yStart + 1; \ - offImage = true; \ - } else { \ - yMax = yStop; \ - jMax = yNum; \ - } - -// Refine limits if the interpolation is exact -#define INTERPOLATE_EXACT() { \ - if (xExact) { \ - if (iMin > 0 || iMax < 1) { \ - INTERPOLATE_OFF(); /* Returns */ \ - } \ - iMin = (xNum - 1) / 2; \ - iMax = iMin + 1; \ - } \ - if (yExact) { \ - if (jMin > 0 || jMax < 1) { \ - INTERPOLATE_OFF(); /* Returns */ \ - } \ - jMin = (yNum - 1) / 2; \ - jMax = jMin + 1; \ - } \ - if (xExact && yExact) { \ - return interpolateFlat(imageValue, varianceValue, maskValue, x, y, options); \ - } \ -} - -// Interpolation engine using interpolation kernel +// Set up the kernel parameters; defines some useful values +#define INTERPOLATE_KERNEL_SETUP(X, Y) \ + /* Central pixel is the pixel below the point of interest */ \ + int xCentral = floor((X) - 0.5), yCentral = floor((Y) - 0.5); /* Central pixel */ \ + float xFrac = (X) - xCentral - 0.5, yFrac = (Y) - yCentral - 0.5; /* Fraction of pixel */ + + +// Interpolation engine for (separable) interpolation kernels static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue, psMaskType *maskValue, float x, float y, - const psImageInterpolateOptions *options) + const psImageInterpolation *interp) { // Parameters have been checked by psImageInterpolate() - int xNum, yNum; // Number of interpolation kernel pixels - int xCentral, yCentral; // Central pixel of the convolution - interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode); - - const psImage *image = options->image; // Image of interest - const psImage *mask = options->mask; // Image mask - const psImage *variance = options->variance; // Image variance - - INTERPOLATE_BOUNDS(); - - float kernel[yNum][xNum]; // Interpolation kernel for straight interpolation - bool xExact, yExact; // Exact shifts? - interpolateKernelGenerate(xNum, yNum, kernel, &xExact, &yExact, xCentral, yCentral, - x, y, options->mode, options->shifting); + psImageInterpolateMode mode = interp->mode; // Mode of interpolation + const psImage *image = interp->image; // Image of interest + const psImage *mask = interp->mask; // Image mask + const psImage *variance = interp->variance; // Image variance + psMaskType maskVal = interp->maskVal; // Value to mask + bool wantVariance = variance && varianceValue; // Does the user want the variance value? + bool haveMask = mask && maskVal; // Does the user want the variance value? + + // Kernel basics + int size = kernelSizes[mode]; // Size of kernel + INTERPOLATE_KERNEL_SETUP(x, y); + + // Extent of the kernel on the image + bool xExact = fabsf(xFrac) < FLT_EPSILON, yExact = fabsf(yFrac) < FLT_EPSILON; // Are shifts exact? + if (xExact && yExact) { + // Both shifts are exact + return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp); + } + + int xLast = image->numCols - 1, yLast = image->numRows - 1; // Last pixels of image + int xStart = xCentral - (size - 1) / 2, xStop = xCentral + size / 2; // Start and stop of kernel on image + int yStart = yCentral - (size - 1) / 2, yStop = yCentral + size / 2; // Start and stop of kernel on image + if (xStart > xLast || xStop < 0 || yStart > yLast || yStop < 0) { + // Interpolation kernel is entirely off the image + INTERPOLATE_OFF(); // Returns + } + int xMin, xMax, yMin, yMax; // Minimum and maximum valid pixels on image + int iMin, iMax, jMin, jMax; // Minimum and maximum valid pixels on kernel + bool offImage = false; // At least one pixel of the kernel is off the image + + // XXX Haven't got single dimension exact shifts implemented properly yet + // XXX When it is ready, note that the check below limits + xExact = yExact = false; + if (xExact) { + if (iMin > 0 || iMax < 1) { + INTERPOLATE_OFF(); // Returns + } + iMin = (size - 1) / 2; + iMax = iMin + 1; + } else { + if (xStart < 0) { + xMin = 0; + iMin = -xStart; + offImage = true; + } else { + xMin = xStart; + iMin = 0; + } + if (xStop > xLast) { + xMax = xLast; + iMax = xLast - xStart + 1; + offImage = true; + } else { + xMax = xStop; + iMax = size; + } + } + if (yExact) { + if (jMin > 0 || jMax < 1) { + INTERPOLATE_OFF(); // Returns + } + jMin = (size - 1) / 2; + jMax = jMin + 1; + } else { + if (yStart < 0) { + yMin = 0; + jMin = -yStart; + offImage = true; + } else { + yMin = yStart; + jMin = 0; + } + if (yStop > yLast) { + yMax = yLast; + jMax = yLast - yStart + 1; + offImage = true; + } else { + yMax = yStop; + jMax = size; + } + } + + // Get the appropriate kernels + psVector *xKernelNew = NULL, *yKernelNew = NULL; // New kernels if not pre-calculated + psVector *xKernel2New = NULL, *yKernel2New = NULL;// New kernel^2 if not pre-calculated + psF32 *xKernel, *yKernel; // Kernel values + psF32 *xKernel2, *yKernel2; // Kernel^2 values + float xSumKernel2, ySumKernel2; // Sum of kernel^2 + int numKernels = interp->numKernels; // Number of pre-calculated kernels + if (numKernels > 0) { + const psImage *kernel = interp->kernel; // Pre-calculated kernels + const psImage *kernel2 = interp->kernel2; // Pre-calculated kernel^2 + const psVector *sumKernel2 = interp->sumKernel2; // Pre-calculated kernel^2 + int xIndex = xFrac * numKernels + 0.5; // Index of closest pre-calculated kernel + xKernel = kernel->data.F32[xIndex]; + xKernel2 = kernel2->data.F32[xIndex]; + xSumKernel2 = sumKernel2->data.F32[xIndex]; + int yIndex = yFrac * numKernels + 0.5; // Index of closest pre-calculated kernel + yKernel = kernel->data.F32[yIndex]; + yKernel2 = kernel2->data.F32[yIndex]; + ySumKernel2 = sumKernel2->data.F32[yIndex]; + } else { + xKernelNew = psVectorAlloc(size, PS_TYPE_F32); + xKernel = xKernelNew->data.F32; + interpolationKernel(xKernel, xFrac, mode); + yKernelNew = psVectorAlloc(size, PS_TYPE_F32); + yKernel = yKernelNew->data.F32; + interpolationKernel(yKernel, yFrac, mode); + + if (haveMask || wantVariance || offImage) { + xKernel2New = psVectorAlloc(size, PS_TYPE_F32); + xKernel2 = xKernel2New->data.F32; + interpolationKernel(xKernel2, xFrac, mode); + yKernel2New = psVectorAlloc(size, PS_TYPE_F32); + yKernel2 = yKernel2New->data.F32; + interpolationKernel(yKernel2, yFrac, mode); + xSumKernel2 = 0.0; + ySumKernel2 = 0.0; + for (int i = 0; i < size; i++) { + xSumKernel2 += xKernel2[i] = PS_SQR(xKernel2[i]); + ySumKernel2 += yKernel2[i] = PS_SQR(yKernel2[i]); + } + } + } float sumKernel = 0.0; // Sum of the kernel - double sumKernel2 = 0.0; // Sum of the kernel-squared - float sumBad = 0.0; // Sum of bad kernel-squared contributions - psMaskType maskVal = options->maskVal; // Value to mask + double sumBad = 0.0; // Sum of bad kernel-squared contributions double sumImage = 0.0; // Sum of image multiplied by kernel double sumVariance = 0.0; // Sum of variance multiplied by kernel-squared - - bool wantVariance = variance && varianceValue; // Does the user want the variance value? - bool haveMask = mask && maskVal; // Does the user want the variance value? - - INTERPOLATE_EXACT(); - + float sumKernel2 = xSumKernel2 * ySumKernel2; // Sum of kernel^2 + + // Measure kernel contribution from outside image + + if (offImage) { // Add contributions in an area outside the image -#define INTERPOLATE_KERNEL_ADD_OFFIMAGE() { \ - sumBad += PS_SQR(kernel[j][i]); \ - } - - // Measure kernel contribution from outside image - if (offImage) { +#define INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL() \ + float xSumBad = 0.0; +#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL() { \ + xSumBad += xKernel2[i]; \ + } +#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW() { \ + sumBad += yKernel2[j] * xSumBad; \ + } + + // Bottom rows if (!yExact) { - // Bottom rows for (int j = 0; j < jMin; j++) { - for (int i = 0; i < xNum; i++) { - INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); + for (int i = 0; i < size; i++) { + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); } + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); } } + // Two sides if (!xExact) { - // Two sides for (int j = jMin; j < jMax; j++) { + INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); for (int i = 0; i < iMin; i++) { - INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); } - for (int i = iMax + 1; i < xNum; i++) { - INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + for (int i = iMax; i < size; i++) { + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); } + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); } } + // Top rows if (!yExact) { - // Top rows - for (int j = jMax + 1; j < yNum; j++) { - for (int i = 0; i < xNum; i++) { - INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + for (int j = jMax; j < size; j++) { + if (!xExact) { + INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); + for (int i = 0; i < size; i++) { + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); + } + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); } } @@ -384,55 +515,112 @@ if (haveMask) { \ /* Variance and mask */ \ - for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ - for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ - float kernelValue = kernel[j][i]; /* Value of kernel */ \ - float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ - if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \ - sumBad += kernelValue2; \ + const psF32 *yKernelData = yKernel; \ + const psF32 *yKernel2Data = yKernel2; \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++, yKernel2Data++) { \ + double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ + double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \ + float xSumKernel = 0.0; /* Sum of kernel in x */ \ + float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \ + /* Dereferenced versions of inputs */ \ + const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ + const ps##TYPE *varianceData = &variance->data.TYPE[yPix][xMin]; \ + const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[yPix][xMin]; \ + const psF32 *xKernelData = xKernel; \ + const psF32 *xKernel2Data = xKernel2; \ + for (int i = iMin, xPix = xMin; i < iMax; \ + i++, xPix++, imageData++, varianceData++, maskData++, \ + xKernelData++, xKernel2Data++) { \ + float kernelValue = *xKernelData; /* Value of kernel in x */ \ + float kernelValue2 = *xKernel2Data; /* Square of kernel in x */ \ + if (*maskData & maskVal) { \ + xSumBad += kernelValue2; \ } else { \ - sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ - sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ - sumKernel += kernelValue; \ - sumKernel2 += kernelValue2; \ + xSumImage += kernelValue * *imageData; \ + xSumVariance += kernelValue2 * *varianceData; \ + xSumKernel += kernelValue; \ } \ } \ + float kernelValue = *yKernel; /* Value of kernel in y */ \ + float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \ + sumImage += kernelValue * xSumImage; \ + sumVariance += kernelValue2 * xSumVariance; \ + sumBad += kernelValue2 * xSumBad; \ + sumKernel += kernelValue * xSumKernel; \ } \ - \ } else { \ /* Variance, no mask */ \ - for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ - for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ - float kernelValue = kernel[j][i]; /* Value of kernel */ \ - float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ - sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ - sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ - sumKernel += kernelValue; \ - sumKernel2 += kernelValue2; \ + const psF32 *yKernelData = yKernel; \ + const psF32 *yKernel2Data = yKernel2; \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++, yKernel2Data++) { \ + double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ + double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \ + float xSumKernel = 0.0; /* Sum of kernel in x */ \ + /* Dereferenced versions of inputs */ \ + const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ + const ps##TYPE *varianceData = &variance->data.TYPE[yPix][xMin]; \ + const psF32 *xKernelData = xKernel; \ + const psF32 *xKernel2Data = xKernel2; \ + for (int i = iMin, xPix = xMin; i < iMax; \ + i++, xPix++, imageData++, varianceData++, xKernelData++, xKernel2Data++) { \ + float kernelValue = *xKernelData; /* Value of kernel */ \ + float kernelValue2 = *xKernel2Data; /* Square of kernel */ \ + xSumImage += kernelValue * *imageData; \ + xSumVariance += kernelValue2 * *varianceData; \ + xSumKernel += kernelValue; \ } \ + float kernelValue = *yKernelData; /* Value of kernel in y */ \ + float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \ + sumImage += kernelValue * xSumImage; \ + sumVariance += kernelValue2 * xSumVariance; \ + sumKernel += kernelValue * xSumKernel; \ } \ } \ } else if (haveMask) { \ /* Mask, no variance */ \ - for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ - for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ - float kernelValue = kernel[j][i]; /* Value of kernel */ \ - float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ - if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \ - sumBad += kernelValue2; \ + const psF32 *yKernelData = yKernel; \ + const psF32 *yKernel2Data = yKernel2; \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++, yKernel2Data++) { \ + double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ + float xSumKernel = 0.0; /* Sum of kernel in x */ \ + float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \ + /* Dereferenced versions of inputs */ \ + const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ + const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[yPix][xMin]; \ + const psF32 *xKernelData = xKernel; \ + const psF32 *xKernel2Data = xKernel2; \ + for (int i = iMin, xPix = xMin; i < iMax; \ + i++, xPix++, imageData++, maskData++, xKernelData++, xKernel2Data++) { \ + float kernelValue = *xKernelData; /* Value of kernel */ \ + float kernelValue2 = *xKernel2Data; /* Value of kernel-squared */ \ + if (*maskData & maskVal) { \ + xSumBad += kernelValue2; \ } else { \ - sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ - sumKernel += kernelValue; \ - sumKernel2 += kernelValue2; \ + xSumImage += kernelValue * *imageData; \ + xSumKernel += kernelValue; \ } \ } \ + float kernelValue = *yKernelData; /* Value of kernel in y */ \ + float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \ + sumImage += kernelValue * xSumImage; \ + sumBad += kernelValue2 * xSumBad; \ + sumKernel += kernelValue * xSumKernel; \ } \ - } else { \ + } else {\ /* Neither variance nor mask */ \ - for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ - for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ - float kernelValue = kernel[j][i]; /* Value of kernel */ \ - sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ - sumKernel += kernelValue; \ + const psF32 *yKernelData = yKernel; \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++) { \ + double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ + float xSumKernel = 0.0; /* Sum of kernel in x */ \ + /* Dereferenced versions of inputs */ \ + const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ + const psF32 *xKernelData = xKernel; \ + for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++, imageData++, xKernelData++) { \ + float kernelValue = *xKernelData; /* Value of kernel */ \ + xSumImage += kernelValue * *imageData; \ + xSumKernel += kernelValue; \ } \ + float kernelValue = *yKernelData; /* Value of kernel in y */ \ + sumImage += kernelValue * xSumImage; \ + sumKernel += kernelValue * xSumKernel; \ } \ } \ @@ -466,8 +654,8 @@ // Completely good pixel status = PS_INTERPOLATE_STATUS_GOOD; - } else if (sumBad < PS_SQR(options->poorFrac) * (sumBad + sumKernel2)) { + } else if (sumBad < PS_SQR(interp->poorFrac) * (sumBad + sumKernel2) ) { // Some pixels masked: poor pixel if (haveMask && maskValue) { - *maskValue |= options->poorMask; + *maskValue |= interp->poorMask; } status = PS_INTERPOLATE_STATUS_POOR; @@ -475,349 +663,28 @@ // Many pixels (or a few important pixels) masked: bad pixel if (haveMask && maskValue) { - *maskValue |= options->badMask; + *maskValue |= interp->badMask; } status = PS_INTERPOLATE_STATUS_BAD; } + psFree(xKernelNew); + psFree(yKernelNew); + psFree(xKernel2New); + psFree(yKernel2New); + return status; } -// Generate Lanczos interpolation kernels -static void lanczos(float values[], // Interpolation kernel to generate - bool *exact, // Exact shift? - int num, // Number of values in the kernel - float frac, // Sub-pixel position - bool allowExact // Allow exact shift? - ) -{ - if (allowExact && fabs(frac) < DBL_EPSILON) { - *exact = true; - // No real shift - for (int i = 0; i < (num - 1) / 2; i++) { - values[i] = 0.0; - } - values[(num - 1) / 2] = 1.0; - for (int i = (num - 1) / 2 + 1; i < num; i++) { - values[i] = 0.0; - } - } else { - *exact = false; - float norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos - float norm2 = M_PI * 4.0 / (float)num; // Normalisation for sinc function 1 - float norm3 = M_PI_2 * 4.0 / (float)num; // Normalisation for sinc function 2 - float pos = - (num - 1)/2 - frac; // Position of interest - for (int i = 0; i < num; i++, pos += 1.0) { - if (fabs(pos) < DBL_EPSILON) { - values[i] = 1.0; - } else { - values[i] = norm1 * sinf(pos * norm2) * sinf(pos * norm3) / PS_SQR(pos); - } - } - } - - return; -} - -// Setup for interpolation by separable kernels -static inline void interpolateSeparateSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned - int *xCentral, int *yCentral, // Central pixel of convolution - float x, float y, // Coordinates of interest - psImageInterpolateMode mode // Mode for interpolation - ) -{ - // Central pixel is the pixel below the point of interest - *xCentral = floor(x - 0.5); - *yCentral = floor(y - 0.5); - switch (mode) { - case PS_INTERPOLATE_LANCZOS2: - *xNum = *yNum = 4; - break; - case PS_INTERPOLATE_LANCZOS3: - *xNum = *yNum = 6; - break; - case PS_INTERPOLATE_LANCZOS4: - *xNum = *yNum = 8; - break; - case PS_INTERPOLATE_FLAT: - case PS_INTERPOLATE_BILINEAR: - case PS_INTERPOLATE_BICUBE: - case PS_INTERPOLATE_GAUSS: - default: - psAbort("Invalid interpolation mode."); - } -} - -// Generate the interpolation kernels for separable case; they should be normalised to unity -static inline void interpolateSeparateGenerate(int xNum, int yNum, // Size of interpolation kernel - float xKernel[xNum], float yKernel[yNum], // Kernels - bool *xExact, bool *yExact, // Exact shifts? - int xCentral, int yCentral, // Central pixel of convolution - float x, float y, // Coordinates of interest - psImageInterpolateMode mode, // Mode for interpolation - bool allowExact // Allow an exact shift? - ) -{ - // XXX Could put in an "exact shift" (i.e., xFrac = 0.0) version - switch (mode) { - case PS_INTERPOLATE_LANCZOS2: - case PS_INTERPOLATE_LANCZOS3: - case PS_INTERPOLATE_LANCZOS4: { - float xFrac = x - xCentral - 0.5; // Fraction of pixel in x - lanczos(xKernel, xExact, xNum, xFrac, allowExact); - float yFrac = y - yCentral - 0.5; // Fraction of pixel in y - lanczos(yKernel, yExact, yNum, yFrac, allowExact); - break; - } - case PS_INTERPOLATE_FLAT: - case PS_INTERPOLATE_BILINEAR: - case PS_INTERPOLATE_BICUBE: - case PS_INTERPOLATE_GAUSS: - default: - psAbort("Invalid interpolation mode."); - } -} - -// Interpolation engine for separable interpolation kernels (either for good reasons or for practical reasons) -static psImageInterpolateStatus interpolateSeparate(double *imageValue, double *varianceValue, - psMaskType *maskValue, float x, float y, - const psImageInterpolateOptions *options) -{ - // Parameters have been checked by psImageInterpolate() - - int xNum, yNum; // Number of interpolation kernel pixels - int xCentral, yCentral; // Central pixel of the convolution - interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode); - - const psImage *image = options->image; // Image of interest - const psImage *mask = options->mask; // Image mask - const psImage *variance = options->variance; // Image variance - - INTERPOLATE_BOUNDS(); - - float xKernel[xNum], yKernel[yNum]; // Interpolation kernels - bool xExact, yExact; // Exact shift? - interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, &xExact, &yExact, - xCentral, yCentral, x, y, options->mode, options->shifting); - - float sumKernel = 0.0; // Sum of the kernel - double sumKernel2 = 0.0; // Sum of the kernel-squared - double sumBad = 0.0; // Sum of bad kernel-squared contributions - psMaskType maskVal = options->maskVal; // Value to mask - double sumImage = 0.0; // Sum of image multiplied by kernel - double sumVariance = 0.0; // Sum of variance multiplied by kernel-squared - - bool wantVariance = variance && varianceValue; // Does the user want the variance value? - bool haveMask = mask && maskVal; // Does the user want the variance value? - - INTERPOLATE_EXACT(); - - // Add contributions in an area outside the image -#define INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL() \ - float xSumBad = 0.0; - -#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL() { \ - xSumBad += PS_SQR(xKernel[i]); \ - } - -#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW() { \ - sumBad += PS_SQR(yKernel[j]) * xSumBad; \ - } - - // Measure kernel contribution from outside image - if (offImage) { - // Bottom rows - if (!yExact) { - for (int j = 0; j < jMin; j++) { - INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); - for (int i = 0; i < xNum; i++) { - INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); - } - INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); - } - } - // Two sides - if (!xExact) { - for (int j = jMin; j < jMax; j++) { - INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); - for (int i = 0; i < iMin; i++) { - INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); - } - for (int i = iMax + 1; i < xNum; i++) { - INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); - } - INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); - } - } - // Top rows - if (!yExact) { - for (int j = jMax + 1; j < yNum; j++) { - if (!xExact) { - INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); - for (int i = 0; i < xNum; i++) { - INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); - } - INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); - } - } - } - } - -#define INTERPOLATE_SEPARATE_CASE(TYPE) \ - case PS_TYPE_##TYPE: { \ - if (wantVariance) { \ - if (haveMask) { \ - /* Variance and mask */ \ - for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ - double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ - double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \ - float xSumKernel = 0.0; /* Sum of kernel in x */ \ - double xSumKernel2 = 0.0; /* Sum of kernel-squared in x */ \ - float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \ - /* Dereferenced versions of inputs */ \ - const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ - const ps##TYPE *varianceData = &variance->data.TYPE[yPix][xMin]; \ - const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[yPix][xMin]; \ - for (int i = iMin, xPix = xMin; i < iMax; \ - i++, xPix++, imageData++, varianceData++, maskData++) { \ - float kernelValue = xKernel[i]; /* Value of kernel in x */ \ - double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel in x */ \ - if (*maskData & maskVal) { \ - xSumBad += kernelValue2; \ - } else { \ - xSumImage += kernelValue * *imageData; \ - xSumVariance += kernelValue2 * *varianceData; \ - xSumKernel += kernelValue; \ - xSumKernel2 += kernelValue2; \ - } \ - } \ - float kernelValue = yKernel[j]; /* Value of kernel in y */ \ - double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared in y */ \ - sumImage += kernelValue * xSumImage; \ - sumVariance += kernelValue2 * xSumVariance; \ - sumBad += kernelValue2 * xSumBad; \ - sumKernel += kernelValue * xSumKernel; \ - sumKernel2 += kernelValue2 * xSumKernel2; \ - } \ - } else { \ - /* Variance, no mask */ \ - for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ - double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ - double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \ - float xSumKernel = 0.0; /* Sum of kernel in x */ \ - double xSumKernel2 = 0.0; /* Sum of kernel-squared in x */ \ - for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ - float kernelValue = xKernel[i]; /* Value of kernel */ \ - double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ - xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ - xSumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ - xSumKernel += kernelValue; \ - xSumKernel2 += kernelValue2; \ - } \ - float kernelValue = yKernel[j]; /* Value of kernel in y */ \ - double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared in y */ \ - sumImage += kernelValue * xSumImage; \ - sumVariance += kernelValue2 * xSumVariance; \ - sumKernel += kernelValue * xSumKernel; \ - sumKernel2 += kernelValue2 * xSumKernel2; \ - } \ - } \ - } else if (haveMask) { \ - /* Mask, no variance */ \ - for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ - double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ - float xSumKernel = 0.0; /* Sum of kernel in x */ \ - double xSumKernel2 = 0.0; /* Sum of kernel-squared in x */ \ - float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \ - for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ - float kernelValue = xKernel[i]; /* Value of kernel */ \ - double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared */ \ - if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \ - xSumBad += kernelValue2; \ - } else { \ - xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ - xSumKernel += kernelValue; \ - xSumKernel2 += kernelValue2; \ - } \ - } \ - float kernelValue = yKernel[j]; /* Value of kernel in y */ \ - double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared in y */ \ - sumImage += kernelValue * xSumImage; \ - sumBad += kernelValue2 * xSumBad; \ - sumKernel += kernelValue * xSumKernel; \ - sumKernel2 += kernelValue2 * xSumKernel2; \ - } \ - } else {\ - /* Neither variance nor mask */ \ - for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ - double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ - float xSumKernel = 0.0; /* Sum of kernel in x */ \ - for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ - float kernelValue = xKernel[i]; /* Value of kernel */ \ - xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ - xSumKernel += kernelValue; \ - } \ - float kernelValue = yKernel[j]; /* Value of kernel in y */ \ - sumImage += kernelValue * xSumImage; \ - sumKernel += kernelValue * xSumKernel; \ - } \ - } \ - } \ - break; - - - switch (image->type.type) { - INTERPOLATE_SEPARATE_CASE(U8); - INTERPOLATE_SEPARATE_CASE(U16); - INTERPOLATE_SEPARATE_CASE(U32); - INTERPOLATE_SEPARATE_CASE(U64); - INTERPOLATE_SEPARATE_CASE(S8); - INTERPOLATE_SEPARATE_CASE(S16); - INTERPOLATE_SEPARATE_CASE(S32); - INTERPOLATE_SEPARATE_CASE(S64); - INTERPOLATE_SEPARATE_CASE(F32); - INTERPOLATE_SEPARATE_CASE(F64); - default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", - image->type.type); - return PS_INTERPOLATE_STATUS_ERROR; - } - - psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; // Status of interpolation - *imageValue = sumImage / sumKernel; - if (wantVariance) { - *varianceValue = sumVariance / sumKernel2; - } - if (sumBad == 0) { - // Completely good pixel - status = PS_INTERPOLATE_STATUS_GOOD; - } else if (sumBad < PS_SQR(options->poorFrac) * (sumBad + sumKernel2) ) { - // Some pixels masked: poor pixel - if (haveMask && maskValue) { - *maskValue |= options->poorMask; - } - status = PS_INTERPOLATE_STATUS_POOR; - } else { - // Many pixels (or a few important pixels) masked: bad pixel - if (haveMask && maskValue) { - *maskValue |= options->badMask; - } - status = PS_INTERPOLATE_STATUS_BAD; - } - - return status; -} - - psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, psMaskType *maskValue, - float x, float y, const psImageInterpolateOptions *options) + float x, float y, const psImageInterpolation *interp) { PS_ASSERT_PTR_NON_NULL(imageValue, PS_INTERPOLATE_STATUS_ERROR); - PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR); - - const psImage *image = options->image; // Image to interpolate - const psImage *mask = options->mask; // Mask to interpolate - const psImage *variance = options->variance; // Variance to interpolate + PS_ASSERT_PTR_NON_NULL(interp, PS_INTERPOLATE_STATUS_ERROR); + + const psImage *image = interp->image; // Image to interpolate + const psImage *mask = interp->mask; // Mask to interpolate + const psImage *variance = interp->variance; // Variance to interpolate PS_ASSERT_IMAGE_NON_NULL(image, PS_INTERPOLATE_STATUS_ERROR); @@ -834,21 +701,20 @@ } - PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(options->poorFrac, 0.0, PS_INTERPOLATE_STATUS_ERROR); - - switch (options->mode) { + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(interp->poorFrac, 0.0, PS_INTERPOLATE_STATUS_ERROR); + + switch (interp->mode) { case PS_INTERPOLATE_FLAT: - return interpolateFlat(imageValue, varianceValue, maskValue, x, y, options); + return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp); case PS_INTERPOLATE_BILINEAR: - case PS_INTERPOLATE_BICUBE: + case PS_INTERPOLATE_BIQUADRATIC: case PS_INTERPOLATE_GAUSS: - return interpolateKernel(imageValue, varianceValue, maskValue, x, y, options); case PS_INTERPOLATE_LANCZOS2: case PS_INTERPOLATE_LANCZOS3: case PS_INTERPOLATE_LANCZOS4: - return interpolateSeparate(imageValue, varianceValue, maskValue, x, y, options); + return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Specified interpolation method (%d) is not supported."), - options->mode); + interp->mode); return PS_INTERPOLATE_STATUS_ERROR; } @@ -859,108 +725,44 @@ -static float varianceFactorFlat(float x, float y, const psImageInterpolateOptions *options) -{ - // There's no smearing - return 1.0; -} - -static float varianceFactorKernel(float x, float y, const psImageInterpolateOptions *options) -{ - int xNum, yNum; // Number of interpolation kernel pixels - int xCentral, yCentral; // Central pixel of the convolution - interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode); - - const psImage *image = options->image; // Image of interest - int xLast = image->numCols - 1; // Last pixel in x - int yLast = image->numRows - 1; // Last pixel in y - - if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast || - yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) { - // At least one pixel of the interpolation kernel is off the image - return NAN; - } - - float kernel[yNum][xNum]; // Interpolation kernel for straight interpolation - bool xExact, yExact; // Exact shift? - // Pretend we're not shifting pixels (i.e., we don't care about any exact shifting) because on the off - // chance that a shift is integral, we don't want to get a trivial kernel back, but something that sort of - // reflects what's going on. - interpolateKernelGenerate(xNum, yNum, kernel, &xExact, &yExact, xCentral, yCentral, - x, y, options->mode, false); - - float sumKernel = 0.0; // Sum of kernel - double sumKernel2 = 0.0; // Sum of kernel squares - for (int j = 0; j < yNum; j++) { - for (int i = 0; i < xNum; i++) { - float value = kernel[j][i]; // Value of kernel - sumKernel += value; - sumKernel2 += PS_SQR(value); - } - } - - return sumKernel2 / PS_SQR(sumKernel); -} - -static float varianceFactorSeparate(float x, float y, const psImageInterpolateOptions *options) -{ - int xNum, yNum; // Number of interpolation kernel pixels - int xCentral, yCentral; // Central pixel of the convolution - interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode); - - const psImage *image = options->image; // Image of interest - int xLast = image->numCols - 1; // Last pixel in x - int yLast = image->numRows - 1; // Last pixel in y - - if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast || - yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) { - // At least one pixel of the interpolation kernel is off the image - return NAN; - } - - float xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation - bool xExact, yExact; // Exact shift? - // Pretend we're not shifting pixels (i.e., we don't care about any exact shifting) because on the off - // chance that a shift is integral, we don't want to get a trivial kernel back, but something that sort of - // reflects what's going on. - interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, &xExact, &yExact, - xCentral, yCentral, x, y, options->mode, false); - - float ySumKernel = 0.0; // Sum of kernel in y - double ySumKernel2 = 0.0; // Sum of kernel squared in y - for (int j = 0; j < yNum; j++) { - float xSumKernel = 0.0; // Sum of kernel in x - double xSumKernel2 = 0.0; // Sum of kernel squared in x - for (int i = 0; i < xNum; i++) { - float value = xKernel[i]; // Value of kernel - xSumKernel += value; - xSumKernel2 += PS_SQR(value); - } - float value = yKernel[j]; // Value of kernel - ySumKernel += xSumKernel * value; - ySumKernel2 += xSumKernel2 * PS_SQR(value); - } - - return ySumKernel2 / PS_SQR(ySumKernel); -} - -float psImageInterpolateVarianceFactor(float x, float y, const psImageInterpolateOptions *options) -{ - PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR); - - switch (options->mode) { +static float varianceFactorKernel(float x, float y, psImageInterpolateMode mode) +{ + int size = kernelSizes[mode]; // Size of kernel + + // Kernel basics + INTERPOLATE_KERNEL_SETUP(x, y); + + psF32 xKernel[size], yKernel[size]; // Interpolation kernels + interpolationKernel(xKernel, xFrac, mode); + interpolationKernel(yKernel, yFrac, mode); + + float xSumKernel = 0.0, ySumKernel = 0.0; // Sum of kernel + float xSumKernel2 = 0.0, ySumKernel2 = 0.0; // Sum of kernel^2 + for (int i = 0; i < size; i++) { + xSumKernel += xKernel[i]; + xSumKernel2 += PS_SQR(xKernel[i]); + ySumKernel += yKernel[i]; + ySumKernel2 += PS_SQR(yKernel[i]); + } + + return (xSumKernel2 * ySumKernel2) / PS_SQR(xSumKernel * ySumKernel); +} + +float psImageInterpolateVarianceFactor(float x, float y, psImageInterpolateMode mode) +{ + switch (mode) { case PS_INTERPOLATE_FLAT: - return varianceFactorFlat(x, y, options); + // No smearing by design + return 1.0; case PS_INTERPOLATE_BILINEAR: - case PS_INTERPOLATE_BICUBE: + case PS_INTERPOLATE_BIQUADRATIC: case PS_INTERPOLATE_GAUSS: - return varianceFactorKernel(x, y, options); case PS_INTERPOLATE_LANCZOS2: case PS_INTERPOLATE_LANCZOS3: case PS_INTERPOLATE_LANCZOS4: - return varianceFactorSeparate(x, y, options); + return varianceFactorKernel(x, y, mode); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Specified interpolation method (%d) is not supported."), - options->mode); + mode); return NAN; } @@ -977,5 +779,5 @@ if (!strcasecmp(name, "FLAT")) return PS_INTERPOLATE_FLAT; if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR; - if (!strcasecmp(name, "BICUBE")) return PS_INTERPOLATE_BICUBE; + if (!strcasecmp(name, "BIQUADRATIC")) return PS_INTERPOLATE_BIQUADRATIC; if (!strcasecmp(name, "GAUSS")) return PS_INTERPOLATE_GAUSS; if (!strcasecmp(name, "LANCZOS2")) return PS_INTERPOLATE_LANCZOS2; Index: /trunk/psLib/src/imageops/psImageInterpolate.h =================================================================== --- /trunk/psLib/src/imageops/psImageInterpolate.h (revision 20305) +++ /trunk/psLib/src/imageops/psImageInterpolate.h (revision 20306) @@ -7,6 +7,6 @@ * @author Paul Price, Institute for Astronomy * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-21 01:12:44 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-22 02:10:37 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii */ @@ -18,12 +18,12 @@ /// Enumeration of options in interpolation typedef enum { - PS_INTERPOLATE_NONE, ///< no interpolate defined (error state) - PS_INTERPOLATE_FLAT, ///< Flat interpolation (nearest pixel) - PS_INTERPOLATE_BILINEAR, ///< Bilinear interpolation - PS_INTERPOLATE_BICUBE, ///< Bicubic interpolation with 3x3 region - PS_INTERPOLATE_GAUSS, ///< Gaussian inteprolation with 3x3 region - PS_INTERPOLATE_LANCZOS2, ///< Sinc interpolation with 4x4 pixel kernel - PS_INTERPOLATE_LANCZOS3, ///< Sinc interpolation with 6x6 pixel kernel - PS_INTERPOLATE_LANCZOS4, ///< Sinc interpolation with 8x8 pixel kernel + PS_INTERPOLATE_NONE = 0, ///< No interpolate defined (error state) + PS_INTERPOLATE_FLAT, ///< Flat interpolation (nearest pixel) + PS_INTERPOLATE_BILINEAR, ///< Bilinear interpolation + PS_INTERPOLATE_BIQUADRATIC, ///< Biquadratic interpolation with 3x3 region + PS_INTERPOLATE_GAUSS, ///< Gaussian inteprolation with 3x3 region + PS_INTERPOLATE_LANCZOS2, ///< Sinc interpolation with 4x4 pixel kernel + PS_INTERPOLATE_LANCZOS3, ///< Sinc interpolation with 6x6 pixel kernel + PS_INTERPOLATE_LANCZOS4, ///< Sinc interpolation with 8x8 pixel kernel } psImageInterpolateMode; @@ -54,30 +54,34 @@ float poorFrac; ///< Fraction of flux in bad pixels before output is marked bad bool shifting; ///< Shifting images? Don't interpolate if the shift is exact. -} psImageInterpolateOptions; + int numKernels; ///< Number of pre-calculated kernels + const psImage *kernel, *kernel2; ///< 1D interpolation kernel and kernel^2 (row) for each spacing + const psVector *sumKernel2; ///< Sum of kernel^2 for each spacing +} psImageInterpolation; /// Allocator -psImageInterpolateOptions *psImageInterpolateOptionsAlloc(psImageInterpolateMode mode, // Interpolation mode - const psImage *image, // Input image - const psImage *variance, // Variance image - const psImage *mask, // Mask image - psMaskType maskVal, // Value to mask - double badImage, // Value for image if bad - double badVariance, // Value for variance if bad - psMaskType badMask, // Mask value for bad pixels - psMaskType poorMask, // Mask value for poor pixels - float poorFrac // Fraction of flux for question +psImageInterpolation *psImageInterpolationAlloc( + psImageInterpolateMode mode, // Interpolation mode + const psImage *image, // Input image + const psImage *variance, // Variance image + const psImage *mask, // Mask image + psMaskType maskVal, // Value to mask + double badImage, // Value for image if bad + double badVariance, // Value for variance if bad + psMaskType badMask, // Mask value for bad pixels + psMaskType poorMask, // Mask value for poor pixels + float poorFrac, // Fraction of flux for question + int numKernels // Number of interpolation kernels to pre-calculate ) PS_ATTR_MALLOC; - /// Interpolate image pixel value given floating point coordinates. -psImageInterpolateStatus psImageInterpolate(double *imageValue, ///< Return value for image - double *varianceValue, ///< Return value for variance - psMaskType *maskValue, ///< Return value for mask - float x, float y, ///< Location to which to interpolate - const psImageInterpolateOptions *options ///< Options +psImageInterpolateStatus psImageInterpolate( + double *imageValue, ///< Return value for image + double *varianceValue, ///< Return value for variance + psMaskType *maskValue, ///< Return value for mask + float x, float y, ///< Location to which to interpolate + const psImageInterpolation *options ///< Options ); - // Return the appropriate interpolation mode given a char string name for that mode @@ -91,5 +95,5 @@ /// factor. float psImageInterpolateVarianceFactor(float x, float y, ///< Position of interest - const psImageInterpolateOptions *options ///< Interpolation options + psImageInterpolateMode mode ///< Interpolation mode ); Index: /trunk/psLib/src/imageops/psImagePixelExtract.c =================================================================== --- /trunk/psLib/src/imageops/psImagePixelExtract.c (revision 20305) +++ /trunk/psLib/src/imageops/psImagePixelExtract.c (revision 20306) @@ -8,6 +8,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-04-04 22:42:02 $ + * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-22 02:10:37 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -679,6 +679,6 @@ float dY = (endRow - startRow) / (float)(nSamples-1); - psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, input, NULL, mask, maskVal, - 0, 0, 0, 0, 0); + psImageInterpolation *interp = psImageInterpolationAlloc(mode, input, NULL, mask, maskVal, + 0, 0, 0, 0, 0, 0); #define LINEAR_CUT_CASE(TYPE) \ Index: /trunk/psLib/test/imageops/tap_psImageInterpolate2.c =================================================================== --- /trunk/psLib/test/imageops/tap_psImageInterpolate2.c (revision 20305) +++ /trunk/psLib/test/imageops/tap_psImageInterpolate2.c (revision 20306) @@ -79,7 +79,6 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 12345); - psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, image, - variance, mask, 0, NAN, NAN, - 0, 0, 0.0); + psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, mask, 0, NAN, NAN, + 0, 0, 0.0, 0); ok(interp, "Interpolation options set"); skip_start(!interp, 2 * num, "Interpolation options failed"); @@ -93,28 +92,19 @@ psMaskType maskVal; - psImageInterpolateStatus status = psImageInterpolate(&imageVal, &varianceVal, &maskVal, x, y, interp); - ok(status != PS_INTERPOLATE_STATUS_ERROR, "Interpolation at %.2f,%.2f (%x)", x, y, status); - skip_start(status == PS_INTERPOLATE_STATUS_ERROR, 1, "Interpolation failed"); + psImageInterpolateStatus interpOK = psImageInterpolate(&imageVal, &varianceVal, &maskVal, + x, y, interp); + ok(interpOK != PS_INTERPOLATE_STATUS_ERROR, "Interpolation at %.2f,%.2f (%x)", x, y, interpOK); + skip_start(interpOK == PS_INTERPOLATE_STATUS_ERROR, 1, "Interpolation failed"); - int xCentral, yCentral; // Central pixel of interpolation - switch (mode) { - case PS_INTERPOLATE_BICUBE: - case PS_INTERPOLATE_GAUSS: - xCentral = x; - yCentral = y; - break; - default: - xCentral = floor(x - 0.5); - yCentral = floor(y - 0.5); - break; - } + int xCentral = x - 0.5, yCentral = y - 0.5; // Central pixel of interpolation if (xCentral - (xKernel - 1) / 2 < 0 || xCentral + xKernel / 2 > xSize - 1 || yCentral - (yKernel - 1) / 2 < 0 || yCentral + yKernel / 2 > ySize - 1) { - ok(status == PS_INTERPOLATE_STATUS_BAD || status == PS_INTERPOLATE_STATUS_POOR, - "Interpolation at border"); + ok(interpOK == PS_INTERPOLATE_STATUS_OFF || interpOK == PS_INTERPOLATE_STATUS_BAD || + interpOK == PS_INTERPOLATE_STATUS_POOR, + "Interpolation at %.2f,%.2f (border): %x", x, y, interpOK); } else { - is_double_tol(imageVal, imageFunc(x, y), tol, "Interpolation = %f vs %f", - imageVal, imageFunc(x, y)); + is_double_tol(imageVal, imageFunc(x, y), tol, "Interpolation = %f vs %f (%d)", + imageVal, imageFunc(x, y), interpOK); } @@ -146,6 +136,6 @@ check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_BILINEAR, 2, 2, 1.0e-6); // 66 tests - check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_BICUBE, 3, 3, 1.0e-6); // 66 tests - check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_BICUBE, 3, 3, 1.0e-6); // 66 tests + check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_BIQUADRATIC, 3, 3, 1.0e-6); // 66 tests + check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_BIQUADRATIC, 3, 3, 1.0e-6); // 66 tests check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_GAUSS, 3, 3, 1.0e-3); // 66 tests