Index: trunk/psLib/src/imageops/psImageInterpolate.c =================================================================== --- trunk/psLib/src/imageops/psImageInterpolate.c (revision 19128) +++ trunk/psLib/src/imageops/psImageInterpolate.c (revision 19129) @@ -7,6 +7,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-06-24 02:04:55 $ + * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-20 04:06:31 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -162,5 +162,6 @@ } -// Generate the interpolation kernel; it should be normalised to unity +// Generate the interpolation kernel +// No need to normalise to unity static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel double kernel[xNum][yNum], // Kernel, to be set @@ -204,15 +205,8 @@ double yFrac = y - yCentral - 0.5; // Fraction of pixel in y double sigma = 0.5; // Gaussian sigma - double norm = 0.0; // Normalisation for (int j = 0, yPos = - (yNum - 1) / 2; j < yNum; j++, yPos++) { for (int i = 0, xPos = - (xNum - 1) / 2; i < xNum; i++, xPos++) { - norm += kernel[j][i] = exp(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) + - PS_SQR(yPos - yFrac))); - } - } - norm = 1.0 / norm; - for (int j = 0; j < yNum; j++) { - for (int i = 0; i < xNum; i++) { - kernel[j][i] *= norm; + kernel[j][i] = exp(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) + + PS_SQR(yPos - yFrac))); } } @@ -228,4 +222,58 @@ } +// 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 */ \ + *imageValue = options->badImage; \ + if (varianceValue) { \ + *varianceValue = options->badVariance; \ + } \ + if (maskValue) { \ + *maskValue = options->badMask; \ + } \ + return PS_INTERPOLATE_STATUS_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 = xStop - xLast; \ + 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 = yStop - yLast; \ + offImage = true; \ + } else { \ + yMax = yStop; \ + jMax = yNum; \ + } + // Interpolation engine using interpolation kernel static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue, @@ -240,126 +288,153 @@ 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 - if (imageValue) { - *imageValue = options->badImage; - } - if (varianceValue) { - *varianceValue = options->badVariance; - } - if (maskValue) { - *maskValue = options->badMask; - } - return PS_INTERPOLATE_STATUS_OFF; - } + const psImage *mask = options->mask; // Image mask + const psImage *variance = options->variance; // Image variance + + INTERPOLATE_BOUNDS(); double kernel[yNum][xNum]; // Interpolation kernel for straight interpolation interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode); - // Image interpolation, according to image type - #define KERNEL_IMAGE_CASE(TYPE) \ - case PS_TYPE_##TYPE: { \ - for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \ - for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \ - *imageValue += kernel[j][i] * image->data.TYPE[yPix][xPix]; \ - } \ - } \ - break; \ - } - - // Calculate the value for the image - if (imageValue) { - *imageValue = 0.0; - switch (image->type.type) { - KERNEL_IMAGE_CASE(U8); - KERNEL_IMAGE_CASE(U16); - KERNEL_IMAGE_CASE(U32); - KERNEL_IMAGE_CASE(U64); - KERNEL_IMAGE_CASE(S8); - KERNEL_IMAGE_CASE(S16); - KERNEL_IMAGE_CASE(S32); - KERNEL_IMAGE_CASE(S64); - KERNEL_IMAGE_CASE(F32); - KERNEL_IMAGE_CASE(F64); - default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", - image->type.type); - return PS_INTERPOLATE_STATUS_ERROR; - } - } - - // Check the mask value - const psImage *mask = options->mask; // Image mask - psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_GOOD; // Status of interpolation - if (maskValue) { - *maskValue = 0; - if (mask) { - psMaskType maskVal = options->maskVal; // Mask value - double badContrib = 0.0; // Amount of kernel on bad pixels - double totContrib = 0.0; // Total amount of kernel - bool badPix = false; // Is there a bad pixel? - for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { - for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { - if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { - badContrib += fabs(kernel[j][i]); - badPix = true; - } - *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix]; - totContrib += fabs(kernel[j][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 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? + + // Add contributions in an area outside the image +#define INTERPOLATE_KERNEL_ADD_OFFIMAGE() { \ + float kernelValue = kernel[j][i]; /* Value of kernel */ \ + double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ + sumBad += kernelValue2; \ + sumKernel2 += kernelValue2; \ + } + + // Measure kernel contribution from outside image + if (offImage) { + // Bottom rows + for (int j = 0; j < jMin; j++) { + for (int i = 0; i < xNum; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); } - - if (badPix) { - if (badContrib / totContrib >= options->poorFrac) { - *maskValue |= options->badMask; - status = PS_INTERPOLATE_STATUS_BAD; - } else { - *maskValue |= options->poorMask; - status = PS_INTERPOLATE_STATUS_POOR; - } + } + // Two sides + for (int j = jMin; j < jMax; j++) { + for (int i = 0; i < iMin; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); } - } - } - - // Finally, the variance - const psImage *variance = options->variance; // Image variance - if (varianceValue && variance) { - *varianceValue = 0.0; - - // Variance interpolation, according to image type - #define KERNEL_VARIANCE_CASE(TYPE) \ - case PS_TYPE_##TYPE: { \ - double sumKernel2 = 0.0; /* Sum of kernel squares */ \ - for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \ - for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \ - double kernel2 = PS_SQR(kernel[j][i]); /* Kernel squared */ \ - sumKernel2 += kernel2; \ - *varianceValue += kernel2 * variance->data.TYPE[yPix][xPix]; \ - } \ - } \ - *varianceValue /= sumKernel2; /* Normalise so that sum of kernel squares is unity */ \ - break; \ + for (int i = iMax + 1; i < xNum; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); } - - switch (variance->type.type) { - KERNEL_VARIANCE_CASE(U8); - KERNEL_VARIANCE_CASE(U16); - KERNEL_VARIANCE_CASE(U32); - KERNEL_VARIANCE_CASE(U64); - KERNEL_VARIANCE_CASE(S8); - KERNEL_VARIANCE_CASE(S16); - KERNEL_VARIANCE_CASE(S32); - KERNEL_VARIANCE_CASE(S64); - KERNEL_VARIANCE_CASE(F32); - KERNEL_VARIANCE_CASE(F64); - default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", - image->type.type); - return PS_INTERPOLATE_STATUS_ERROR; - } + } + // Top rows + for (int j = jMax + 1; j < yNum; j++) { + for (int i = 0; i < xNum; 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 */ \ + double 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 += PS_SQR(kernelValue); \ + } \ + } \ + \ + } 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 */ \ + double 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 */ \ + if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \ + sumBad += PS_SQR(kernelValue); \ + } else { \ + sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ + sumKernel += kernelValue; \ + } \ + sumKernel2 += PS_SQR(kernelValue); \ + } \ + } \ + } 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); + INTERPOLATE_KERNEL_CASE(U32); + INTERPOLATE_KERNEL_CASE(U64); + INTERPOLATE_KERNEL_CASE(S8); + INTERPOLATE_KERNEL_CASE(S16); + INTERPOLATE_KERNEL_CASE(S32); + INTERPOLATE_KERNEL_CASE(S64); + INTERPOLATE_KERNEL_CASE(F32); + INTERPOLATE_KERNEL_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 * PS_SQR(sumKernel) / sumKernel2; + } + if (sumBad == 0) { + // Completely good pixel + status = PS_INTERPOLATE_STATUS_GOOD; + } else if (sumBad / sumKernel2 < PS_SQR(options->poorFrac)) { + // 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; } @@ -472,139 +547,205 @@ 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 - if (imageValue) { - *imageValue = options->badImage; - } - if (varianceValue) { - *varianceValue = options->badVariance; - } - if (maskValue) { - *maskValue = options->badMask; - } - return PS_INTERPOLATE_STATUS_OFF; - } + const psImage *mask = options->mask; // Image mask + const psImage *variance = options->variance; // Image variance + + INTERPOLATE_BOUNDS(); double xKernel[xNum], yKernel[yNum]; // Interpolation kernels interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode); - // Image interpolation, according to image type - #define SEPARATE_IMAGE_CASE(TYPE) \ - case PS_TYPE_##TYPE: { \ - for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \ - double xInterpValue = 0.0; /* Interpolation in x */ \ - for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \ - xInterpValue += xKernel[i] * image->data.TYPE[yPix][xPix]; \ - } \ - *imageValue += yKernel[j] * xInterpValue; /* Interpolating in y */ \ - } \ - break; \ - } - - // Calculate the value for the image - if (imageValue) { - *imageValue = 0.0; - switch (image->type.type) { - SEPARATE_IMAGE_CASE(U8); - SEPARATE_IMAGE_CASE(U16); - SEPARATE_IMAGE_CASE(U32); - SEPARATE_IMAGE_CASE(U64); - SEPARATE_IMAGE_CASE(S8); - SEPARATE_IMAGE_CASE(S16); - SEPARATE_IMAGE_CASE(S32); - SEPARATE_IMAGE_CASE(S64); - SEPARATE_IMAGE_CASE(F32); - SEPARATE_IMAGE_CASE(F64); - default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", - image->type.type); - return PS_INTERPOLATE_STATUS_ERROR; - } - } - - // Check the mask value - const psImage *mask = options->mask; // Image mask - psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_GOOD; // Status of interpolation - if (maskValue) { - *maskValue = 0; - if (mask) { - psMaskType maskVal = options->maskVal; // Mask value - double badContrib = 0.0; // Amount of kernel on bad pixels - double totContrib = 0.0; // Total amount of kernel - bool badPix = false; // Number of bad pixels - for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { - // Interpolation in x - double xBadContrib = 0.0; // Amount of x kernel on bad pixels - double xTotContrib = 0.0; // Total amount of x kernel - for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { - if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { - xBadContrib += fabs(xKernel[i]); - badPix = true; - } - *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix]; - xTotContrib += fabs(xKernel[i]); - } - // Interpolating in y - badContrib += fabs(yKernel[j]) * xBadContrib; - totContrib += fabs(yKernel[j]) * xTotContrib; + 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 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? + + // Add contributions in an area outside the image +#define INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL() \ + float xSumBad = 0.0; \ + double xSumKernel2 = 0.0; +#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL() { \ + float kernelValue = xKernel[i]; /* Value of kernel */ \ + double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ + xSumBad += kernelValue2; \ + xSumKernel2 += kernelValue2; \ + } +#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW() { \ + float kernelValue = yKernel[j]; /* Value of kernel */ \ + double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ + sumBad += kernelValue2 * xSumBad; \ + sumKernel2 += kernelValue2 * xSumKernel2; \ + } + + // Measure kernel contribution from outside image + if (offImage) { + // Bottom rows + for (int j = 0; j < jMin; j++) { + INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); + for (int i = 0; i < xNum; i++) { + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); } - - if (badPix) { - if (badContrib / totContrib >= options->poorFrac) { - *maskValue |= options->badMask; - status = PS_INTERPOLATE_STATUS_BAD; - } else { - *maskValue |= options->poorMask; - status = PS_INTERPOLATE_STATUS_POOR; - } + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); + } + // Two sides + for (int j = jMin; j < jMax; j++) { + INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); + for (int i = 0; i < iMin; i++) { + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); } - } - } - - // Finally, the variance - const psImage *variance = options->variance; // Image variance - if (varianceValue && variance) { - *varianceValue = 0.0; - - // Variance interpolation, according to image type - #define SEPARATE_VARIANCE_CASE(TYPE) \ - case PS_TYPE_##TYPE: { \ - double ySumKernel2 = 0.0; /* Sum of kernel squared in y */ \ - for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \ - double xSumKernel2 = 0.0; /* Sum of kernel squared in x */ \ - double xInterpValue = 0.0; /* Interpolation in x */ \ - for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \ - double kernel2 = PS_SQR(xKernel[i]); /* Kernel squared */ \ - xSumKernel2 += kernel2; \ - xInterpValue += kernel2 * variance->data.TYPE[yPix][xPix]; \ - } \ - double kernel2 = PS_SQR(yKernel[j]); /* Kernel squared */ \ - ySumKernel2 += xSumKernel2 * kernel2; \ - *varianceValue += xInterpValue * kernel2; /* Interpolating in y */ \ - } \ - *varianceValue /= ySumKernel2; \ - break; \ - } - - switch (variance->type.type) { - SEPARATE_VARIANCE_CASE(U8); - SEPARATE_VARIANCE_CASE(U16); - SEPARATE_VARIANCE_CASE(U32); - SEPARATE_VARIANCE_CASE(U64); - SEPARATE_VARIANCE_CASE(S8); - SEPARATE_VARIANCE_CASE(S16); - SEPARATE_VARIANCE_CASE(S32); - SEPARATE_VARIANCE_CASE(S64); - SEPARATE_VARIANCE_CASE(F32); - SEPARATE_VARIANCE_CASE(F64); - default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", - variance->type.type); - return PS_INTERPOLATE_STATUS_ERROR; - } + for (int i = iMax + 1; i < xNum; i++) { + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); + } + INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); + } + // Top rows + for (int j = jMax + 1; j < yNum; j++) { + 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 */ \ + for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ + float kernelValue = xKernel[i]; /* Value of kernel in x */ \ + double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel in x */ \ + if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \ + xSumBad += kernelValue2; \ + } else { \ + 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; \ + 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 * PS_SQR(sumKernel) / sumKernel2; + } + if (sumBad == 0) { + // Completely good pixel + status = PS_INTERPOLATE_STATUS_GOOD; + } else if (sumBad / sumKernel2 < PS_SQR(options->poorFrac)) { + // 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; } @@ -617,4 +758,5 @@ float x, float y, const psImageInterpolateOptions *options) { + PS_ASSERT_PTR_NON_NULL(imageValue, PS_INTERPOLATE_STATUS_ERROR); PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR); Index: trunk/psLib/test/imageops/tap_psImageInterpolate2.c =================================================================== --- trunk/psLib/test/imageops/tap_psImageInterpolate2.c (revision 19128) +++ trunk/psLib/test/imageops/tap_psImageInterpolate2.c (revision 19129) @@ -86,6 +86,6 @@ for (int i = 0; i < num; i++) { - float x = psRandomUniform(rng) * xSize; - float y = psRandomUniform(rng) * ySize; + float x = psRandomUniform(rng) * (xSize - 1); + float y = psRandomUniform(rng) * (ySize - 1); // Values from interpolation @@ -93,7 +93,7 @@ psMaskType maskVal; - bool success = psImageInterpolate(&imageVal, &varianceVal, &maskVal, x, y, interp); - ok(success, "Interpolation at %.2f,%.2f", x, y); - skip_start(!success, 1, "Interpolation failed"); + 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"); int xCentral, yCentral; // Central pixel of interpolation @@ -112,5 +112,6 @@ if (xCentral - (xKernel - 1) / 2 < 0 || xCentral + xKernel / 2 > xSize - 1 || yCentral - (yKernel - 1) / 2 < 0 || yCentral + yKernel / 2 > ySize - 1) { - ok(isnan(imageVal), "Interpolation = %f vs NAN (border)", imageVal); + ok(status == PS_INTERPOLATE_STATUS_BAD || status == PS_INTERPOLATE_STATUS_POOR, + "Interpolation at border"); } else { is_double_tol(imageVal, imageFunc(x, y), tol, "Interpolation = %f vs %f",