Index: /trunk/psLib/src/imageops/psImageInterpolate.c =================================================================== --- /trunk/psLib/src/imageops/psImageInterpolate.c (revision 19139) +++ /trunk/psLib/src/imageops/psImageInterpolate.c (revision 19140) @@ -7,6 +7,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-20 04:06:31 $ + * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-21 01:12:44 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -60,4 +60,5 @@ options->poorMask = poorMask; options->poorFrac = poorFrac; + options->shifting = true; return options; @@ -166,13 +167,18 @@ static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel double 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 + psImageInterpolateMode mode, // Mode for interpolation + bool allowExact // Allow exact shifts? ) { + double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x + double 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: { - double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x - double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y kernel[0][0] = (1.0 - xFrac) * (1.0 - yFrac); kernel[0][1] = xFrac * (1.0 - yFrac); @@ -220,5 +226,18 @@ psAbort("Invalid interpolation mode."); } -} + return; +} + +// Can't interpolate: kernel is off the image +#define INTERPOLATE_OFF() { \ + *imageValue = options->badImage; \ + if (varianceValue) { \ + *varianceValue = options->badVariance; \ + } \ + if (maskValue) { \ + *maskValue = options->badMask; \ + } \ + return PS_INTERPOLATE_STATUS_OFF; \ + } // Set up and check interpolation bounds @@ -231,12 +250,5 @@ 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; \ + INTERPOLATE_OFF(); \ } \ int xMin, xMax, yMin, yMax; /* Minimum and maximum valid pixels on image */ \ @@ -276,4 +288,25 @@ } +// 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 static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue, @@ -294,5 +327,7 @@ double kernel[yNum][xNum]; // Interpolation kernel for straight interpolation - interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode); + bool xExact, yExact; // Exact shifts? + interpolateKernelGenerate(xNum, yNum, kernel, &xExact, &yExact, xCentral, yCentral, + x, y, options->mode, options->shifting); float sumKernel = 0.0; // Sum of the kernel @@ -306,4 +341,6 @@ bool haveMask = mask && maskVal; // Does the user want the variance value? + INTERPOLATE_EXACT(); + // Add contributions in an area outside the image #define INTERPOLATE_KERNEL_ADD_OFFIMAGE() { \ @@ -316,23 +353,29 @@ // 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 (!yExact) { + // Bottom rows + for (int j = 0; j < jMin; j++) { + for (int i = 0; i < xNum; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } } } - // Two sides - for (int j = jMin; j < jMax; j++) { - for (int i = 0; i < iMin; 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 + 1; i < xNum; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } } - for (int i = iMax + 1; i < xNum; i++) { - INTERPOLATE_KERNEL_ADD_OFFIMAGE(); - } - } - // Top rows - for (int j = jMax + 1; j < yNum; j++) { - for (int i = 0; i < xNum; i++) { - INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } + if (!yExact) { + // Top rows + for (int j = jMax + 1; j < yNum; j++) { + for (int i = 0; i < xNum; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } } } @@ -445,11 +488,12 @@ // Generate Lanczos interpolation kernels static void lanczos(double values[], // Interpolation kernel to generate + bool *exact, // Exact shift? int num, // Number of values in the kernel - float frac // Sub-pixel position + float frac, // Sub-pixel position + bool allowExact // Allow exact shift? ) { - // XXX: Instead of generating a convolution kernel that does no shifting, need to set a boolean that says - // we can do an exact shift. - if (fabs(frac) < DBL_EPSILON) { + if (allowExact && fabs(frac) < DBL_EPSILON) { + *exact = true; // No real shift for (int i = 0; i < (num - 1) / 2; i++) { @@ -461,5 +505,5 @@ } } else { - double norm4 = 0.0; // Normalisation + *exact = false; double norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos double norm2 = M_PI * 4.0 / (float)num; // Normalisation for sinc function 1 @@ -467,9 +511,9 @@ double pos = - (num - 1)/2 - frac; // Position of interest for (int i = 0; i < num; i++, pos += 1.0) { - norm4 += values[i] = norm1 * sin(pos * norm2) * sin(pos * norm3) / PS_SQR(pos); - } - norm4 = 1.0 / norm4; - for (int i = 0; i < num; i++) { - values[i] *= norm4; + if (fabs(pos) < DBL_EPSILON) { + values[i] = 1.0; + } else { + values[i] = norm1 * sin(pos * norm2) * sin(pos * norm3) / PS_SQR(pos); + } } } @@ -510,7 +554,9 @@ static inline void interpolateSeparateGenerate(int xNum, int yNum, // Size of interpolation kernel double xKernel[xNum], double 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 + psImageInterpolateMode mode, // Mode for interpolation + bool allowExact // Allow an exact shift? ) { @@ -521,7 +567,7 @@ case PS_INTERPOLATE_LANCZOS4: { double xFrac = x - xCentral - 0.5; // Fraction of pixel in x - lanczos(xKernel, xNum, xFrac); + lanczos(xKernel, xExact, xNum, xFrac, allowExact); double yFrac = y - yCentral - 0.5; // Fraction of pixel in y - lanczos(yKernel, yNum, yFrac); + lanczos(yKernel, yExact, yNum, yFrac, allowExact); break; } @@ -553,5 +599,7 @@ double xKernel[xNum], yKernel[yNum]; // Interpolation kernels - interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode); + 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 @@ -564,4 +612,6 @@ 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 @@ -585,29 +635,37 @@ 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 (!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(); } - 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(); + 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(); } - for (int i = iMax + 1; i < xNum; i++) { - INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); + } + // 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(); + } } - 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(); } } @@ -826,5 +884,10 @@ double kernel[yNum][xNum]; // Interpolation kernel for straight interpolation - interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode); + 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); double sumKernel2 = 0.0; // Sum of kernel squares @@ -855,5 +918,10 @@ double xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation - interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode); + 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); double ySumKernel2 = 0.0; // Sum of kernel squared in y Index: /trunk/psLib/src/imageops/psImageInterpolate.h =================================================================== --- /trunk/psLib/src/imageops/psImageInterpolate.h (revision 19139) +++ /trunk/psLib/src/imageops/psImageInterpolate.h (revision 19140) @@ -7,6 +7,6 @@ * @author Paul Price, Institute for Astronomy * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-06-17 21:24:58 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-21 01:12:44 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii */ @@ -53,4 +53,5 @@ psMaskType poorMask; ///< Mask value to give poor pixels 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;