Index: trunk/psLib/src/imageops/psImageInterpolate.c =================================================================== --- trunk/psLib/src/imageops/psImageInterpolate.c (revision 20104) +++ trunk/psLib/src/imageops/psImageInterpolate.c (revision 20207) @@ -7,6 +7,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-10-13 21:47:41 $ + * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-16 22:55:45 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -166,5 +166,5 @@ // 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 + float kernel[xNum][yNum], // Kernel, to be set bool *xExact, bool *yExact, // Exact shift? int xCentral, int yCentral, // Central pixel of convolution @@ -174,6 +174,6 @@ ) { - double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x - double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y + 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); @@ -188,6 +188,6 @@ } case PS_INTERPOLATE_BICUBE: { - double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x - double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y + 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; @@ -208,11 +208,11 @@ } case PS_INTERPOLATE_GAUSS: { - double xFrac = x - xCentral - 0.5; // Fraction of pixel in x - double yFrac = y - yCentral - 0.5; // Fraction of pixel in y - double sigma = 0.5; // Gaussian sigma + 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] = exp(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) + - PS_SQR(yPos - yFrac))); + kernel[j][i] = expf(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) + + PS_SQR(yPos - yFrac))); } } @@ -326,5 +326,5 @@ INTERPOLATE_BOUNDS(); - double kernel[yNum][xNum]; // Interpolation kernel for straight interpolation + float kernel[yNum][xNum]; // Interpolation kernel for straight interpolation bool xExact, yExact; // Exact shifts? interpolateKernelGenerate(xNum, yNum, kernel, &xExact, &yExact, xCentral, yCentral, @@ -345,8 +345,5 @@ // 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; \ + sumBad += PS_SQR(kernel[j][i]); \ } @@ -390,5 +387,5 @@ 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 */ \ + float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \ sumBad += kernelValue2; \ @@ -397,6 +394,6 @@ sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ sumKernel += kernelValue; \ + sumKernel2 += kernelValue2; \ } \ - sumKernel2 += PS_SQR(kernelValue); \ } \ } \ @@ -407,5 +404,5 @@ 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 */ \ + float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ @@ -420,11 +417,12 @@ 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 += PS_SQR(kernelValue); \ + sumBad += kernelValue2; \ } else { \ sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ sumKernel += kernelValue; \ + sumKernel2 += kernelValue2; \ } \ - sumKernel2 += PS_SQR(kernelValue); \ } \ } \ @@ -463,10 +461,10 @@ *imageValue = sumImage / sumKernel; if (wantVariance) { - *varianceValue = sumVariance * PS_SQR(sumKernel) / sumKernel2; + *varianceValue = sumVariance / sumKernel2; } if (sumBad == 0) { // Completely good pixel status = PS_INTERPOLATE_STATUS_GOOD; - } else if (sumBad / sumKernel2 < PS_SQR(options->poorFrac)) { + } else if (sumBad < PS_SQR(options->poorFrac) * (sumBad + sumKernel2)) { // Some pixels masked: poor pixel if (haveMask && maskValue) { @@ -487,5 +485,5 @@ // Generate Lanczos interpolation kernels -static void lanczos(double values[], // Interpolation kernel to generate +static void lanczos(float values[], // Interpolation kernel to generate bool *exact, // Exact shift? int num, // Number of values in the kernel @@ -506,13 +504,13 @@ } else { *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 - double norm3 =M_PI_2 * 4.0 / (float)num; // Normalisation for sinc function 2 - double pos = - (num - 1)/2 - frac; // Position of interest + 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 * sin(pos * norm2) * sin(pos * norm3) / PS_SQR(pos); + values[i] = norm1 * sinf(pos * norm2) * sinf(pos * norm3) / PS_SQR(pos); } } @@ -553,5 +551,5 @@ // 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 - double xKernel[xNum], double yKernel[yNum], // Kernels + float xKernel[xNum], float yKernel[yNum], // Kernels bool *xExact, bool *yExact, // Exact shifts? int xCentral, int yCentral, // Central pixel of convolution @@ -566,7 +564,7 @@ case PS_INTERPOLATE_LANCZOS3: case PS_INTERPOLATE_LANCZOS4: { - double xFrac = x - xCentral - 0.5; // Fraction of pixel in x + float xFrac = x - xCentral - 0.5; // Fraction of pixel in x lanczos(xKernel, xExact, xNum, xFrac, allowExact); - double yFrac = y - yCentral - 0.5; // Fraction of pixel in y + float yFrac = y - yCentral - 0.5; // Fraction of pixel in y lanczos(yKernel, yExact, yNum, yFrac, allowExact); break; @@ -598,5 +596,5 @@ INTERPOLATE_BOUNDS(); - double xKernel[xNum], yKernel[yNum]; // Interpolation kernels + float xKernel[xNum], yKernel[yNum]; // Interpolation kernels bool xExact, yExact; // Exact shift? interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, &xExact, &yExact, @@ -605,5 +603,5 @@ 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 + 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 @@ -617,17 +615,12 @@ // Add contributions in an area outside the image #define INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL() \ - float xSumBad = 0.0; \ - double xSumKernel2 = 0.0; + float xSumBad = 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; \ - } + xSumBad += PS_SQR(xKernel[i]); \ + } + #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; \ + sumBad += PS_SQR(yKernel[j]) * xSumBad; \ } @@ -696,6 +689,6 @@ xSumVariance += kernelValue2 * *varianceData; \ xSumKernel += kernelValue; \ + xSumKernel2 += kernelValue2; \ } \ - xSumKernel2 += kernelValue2; \ } \ float kernelValue = yKernel[j]; /* Value of kernel in y */ \ @@ -745,6 +738,6 @@ xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ xSumKernel += kernelValue; \ + xSumKernel2 += kernelValue2; \ } \ - xSumKernel2 += kernelValue2; \ } \ float kernelValue = yKernel[j]; /* Value of kernel in y */ \ @@ -794,10 +787,10 @@ *imageValue = sumImage / sumKernel; if (wantVariance) { - *varianceValue = sumVariance * PS_SQR(sumKernel) / sumKernel2; + *varianceValue = sumVariance / sumKernel2; } if (sumBad == 0) { // Completely good pixel status = PS_INTERPOLATE_STATUS_GOOD; - } else if (sumBad / sumKernel2 < PS_SQR(options->poorFrac)) { + } else if (sumBad < PS_SQR(options->poorFrac) * (sumBad + sumKernel2) ) { // Some pixels masked: poor pixel if (haveMask && maskValue) { @@ -888,5 +881,5 @@ } - double kernel[yNum][xNum]; // Interpolation kernel for straight interpolation + 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 @@ -896,5 +889,5 @@ x, y, options->mode, false); - double sumKernel = 0.0; // Sum of kernel + float sumKernel = 0.0; // Sum of kernel double sumKernel2 = 0.0; // Sum of kernel squares for (int j = 0; j < yNum; j++) { @@ -925,5 +918,5 @@ } - double xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation + 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 @@ -933,8 +926,8 @@ xCentral, yCentral, x, y, options->mode, false); - double ySumKernel = 0.0; // Sum of kernel in y + 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++) { - double xSumKernel = 0.0; // Sum of kernel in x + 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++) {