Index: trunk/psLib/src/imageops/psImageInterpolate.c =================================================================== --- trunk/psLib/src/imageops/psImageInterpolate.c (revision 20306) +++ trunk/psLib/src/imageops/psImageInterpolate.c (revision 20311) @@ -7,6 +7,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-10-22 02:10:37 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-22 02:48:50 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -55,4 +55,5 @@ } +#if 0 /// Generate a biquadratic interpolation kernel /// This reduces Gene's original made-up 2D kernel to 1D @@ -67,4 +68,28 @@ kernel[2] = x + xx; } +#else +/// Generate a biquadratic interpolation kernel +static inline void interpolationKernelBiquadratic( + psF32 kernel[kernelSizes[PS_INTERPOLATE_BIQUADRATIC]][kernelSizes[PS_INTERPOLATE_BIQUADRATIC]], // Kernel + float xFrac, float yFrac // Fraction of pixel + ) +{ + 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; +} +#endif + #if 0 @@ -135,5 +160,5 @@ -// Generate interpolation kernel +// Generate 1D interpolation kernel static inline void interpolationKernel(psF32 *kernel, // Kernel vector to populate float frac, // Fraction of pixel @@ -148,7 +173,4 @@ case PS_INTERPOLATE_BILINEAR: interpolationKernelBilinear(kernel, frac); - break; - case PS_INTERPOLATE_BIQUADRATIC: - interpolationKernelBiquadratic(kernel, frac); break; case PS_INTERPOLATE_GAUSS: @@ -160,4 +182,5 @@ interpolationKernelLanczos(kernel, kernelSizes[mode], frac); break; + case PS_INTERPOLATE_BIQUADRATIC: // 2D kernel default: psAbort("Unsupported interpolation mode: %x", mode); @@ -207,22 +230,37 @@ 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]); + + switch (mode) { + case PS_INTERPOLATE_BILINEAR: + case PS_INTERPOLATE_GAUSS: + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: + 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; } - sumKernel2->data.F32[i] = sum; - } - } + } + break; + case PS_INTERPOLATE_BIQUADRATIC: + // 2D kernel, would cost too much memory to pre-calculate + break; + default: + psAbort("Unsupported interpolation mode: %x", mode); + } + interp->kernel = kernel; interp->kernel2 = kernel2; @@ -316,14 +354,104 @@ // Set up the kernel parameters; defines some useful values -#define INTERPOLATE_KERNEL_SETUP(X, Y) \ +#define INTERPOLATE_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 psImageInterpolation *interp) + float xFrac = (X) - xCentral - 0.5, yFrac = (Y) - yCentral - 0.5; /* Fraction of pixel */ \ + bool xExact = fabsf(xFrac) < FLT_EPSILON, yExact = fabsf(yFrac) < FLT_EPSILON; /* Are shifts exact? */ \ + +// Determine the range of the kernel; defines some useful values +#define INTERPOLATE_RANGE() \ + int xLast = image->numCols - 1, yLast = image->numRows - 1; /* Last pixels of image */ \ + /* Proper range of kernel on image; not all may be defined */ \ + int xStart = xCentral - (size - 1) / 2, xStop = xCentral + size / 2; \ + int yStart = yCentral - (size - 1) / 2, yStop = yCentral + size / 2; \ + 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 limit checks ([ij]{Min,Max}) below are wrong: */ \ + /* should probably refer to [xy]{Start,Stop}. */ \ + xExact = yExact = false; \ + if (xExact) { \ + if (iMin > 0 || iMax < 1) { /* This is wrong */ \ + 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) { /* This is wrong */ \ + 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; \ + } \ + } + +// Determine the result of the interpolation after all the math has been done +#define INTERPOLATE_RESULT() \ + 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(interp->poorFrac) * (sumBad + sumKernel2) ) { \ + /* Some pixels masked: poor pixel */ \ + if (haveMask && maskValue) { \ + *maskValue |= interp->poorMask; \ + } \ + status = PS_INTERPOLATE_STATUS_POOR; \ + } else { \ + /* Many pixels (or a few important pixels) masked: bad pixel */ \ + if (haveMask && maskValue) { \ + *maskValue |= interp->badMask; \ + } \ + status = PS_INTERPOLATE_STATUS_BAD; \ + } + +// Interpolation engine for separable interpolation kernels +static psImageInterpolateStatus interpolateSeparable(double *imageValue, double *varianceValue, + psMaskType *maskValue, float x, float y, + const psImageInterpolation *interp) { // Parameters have been checked by psImageInterpolate() @@ -339,75 +467,9 @@ // 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? + INTERPOLATE_SETUP(x, y); 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; - } - } + INTERPOLATE_RANGE(); // Get the appropriate kernels @@ -510,5 +572,5 @@ } -#define INTERPOLATE_KERNEL_CASE(TYPE) \ +#define INTERPOLATE_SEPARATE_CASE(TYPE) \ case PS_TYPE_##TYPE: { \ if (wantVariance) { \ @@ -630,4 +692,162 @@ 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; + } + + INTERPOLATE_RESULT(); + + psFree(xKernelNew); + psFree(yKernelNew); + psFree(xKernel2New); + psFree(yKernel2New); + + return status; +} + +// Interpolation engine for (separable) interpolation kernels +static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue, + psMaskType *maskValue, float x, float y, + const psImageInterpolation *interp) +{ + // Parameters have been checked by psImageInterpolate() + + 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_SETUP(x, y); + if (xExact && yExact) { + return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp); + } + INTERPOLATE_RANGE(); + + // Get the appropriate kernels + psF32 kernel[size][size]; + psAssert(mode == PS_INTERPOLATE_BIQUADRATIC, "Mode is %x", mode); + interpolationKernelBiquadratic(kernel, xFrac, yFrac); + + float sumKernel = 0.0; // Sum of the kernel + 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 + float sumKernel2 = 0.0; // Sum of kernel^2 + + // 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) { + if (!yExact) { + // Bottom rows + for (int j = 0; j < jMin; j++) { + for (int i = 0; i < size; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } + } + } + if (!xExact) { + // Two sides + for (int j = jMin; j < jMax; j++) { + for (int i = 0; i < iMin; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } + for (int i = iMax; i < size; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } + } + } + if (!yExact) { + // Top rows + for (int j = jMax; j < size; j++) { + for (int i = 0; i < size; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } + } + } + } + +#define INTERPOLATE_KERNEL_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + if (wantVariance) { \ + 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; \ + } else { \ + sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ + sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ + sumKernel += kernelValue; \ + sumKernel2 += kernelValue2; \ + } \ + } \ + } \ + \ + } 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; \ + } \ + } \ + } \ + } 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; \ + } else { \ + sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ + sumKernel += kernelValue; \ + sumKernel2 += kernelValue2; \ + } \ + } \ + } \ + } 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; \ + } \ + } \ + } \ + } \ + break; + + switch (image->type.type) { INTERPOLATE_KERNEL_CASE(U8); INTERPOLATE_KERNEL_CASE(U16); @@ -646,30 +866,5 @@ } - 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(interp->poorFrac) * (sumBad + sumKernel2) ) { - // Some pixels masked: poor pixel - if (haveMask && maskValue) { - *maskValue |= interp->poorMask; - } - status = PS_INTERPOLATE_STATUS_POOR; - } else { - // Many pixels (or a few important pixels) masked: bad pixel - if (haveMask && maskValue) { - *maskValue |= interp->badMask; - } - status = PS_INTERPOLATE_STATUS_BAD; - } - - psFree(xKernelNew); - psFree(yKernelNew); - psFree(xKernel2New); - psFree(yKernel2New); + INTERPOLATE_RESULT(); return status; @@ -706,11 +901,12 @@ case PS_INTERPOLATE_FLAT: return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp); + case PS_INTERPOLATE_BIQUADRATIC: + return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp); case PS_INTERPOLATE_BILINEAR: - case PS_INTERPOLATE_BIQUADRATIC: case PS_INTERPOLATE_GAUSS: case PS_INTERPOLATE_LANCZOS2: case PS_INTERPOLATE_LANCZOS3: case PS_INTERPOLATE_LANCZOS4: - return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp); + return interpolateSeparable(imageValue, varianceValue, maskValue, x, y, interp); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, @@ -730,5 +926,6 @@ // Kernel basics - INTERPOLATE_KERNEL_SETUP(x, y); + INTERPOLATE_SETUP(x, y); + xExact = yExact = false; psF32 xKernel[size], yKernel[size]; // Interpolation kernels